热–力耦合作用下的冷冻消融手术模拟研究

A Simulation Study of Cryoablation Surgery under the Effect of Heat-Force Coupling

摘要: 本文通过COMSOL软件建立了三维肝肿瘤模型，研究了热–力耦合作用下的冷冻消融手术过程。结果表明，在线弹性模型为前提的情况下，只有探针周围极小区域的应力超过了屈服应力，与不考虑力场的数值模拟相比，该模型中组织相变导致的变形对冷冻范围产生了约10%的影响。

Abstract: In this paper, a three-dimensional liver tumor model was established by COMSOL software to study the cryoablation procedure under thermal-force coupling. The results show that only a very small area around the probe exceeds the yield stress under the premise of the in-line elastic model, and the deformation due to the tissue phase change in this model exerts about 10% influence on the cryoablation extent compared with the numerical simulation without considering the force field.

1. 引言

肿瘤已成为一个重要的健康问题，影响着人们的健康。虽然良性肿瘤不会对人类健康构成直接威胁，但仍有可能令人担忧。相反，恶性肿瘤，即通常所说的癌症，对人们的生活和整体健康有着深远的影响。癌症是导致死亡的主要原因，其发病率在各国都很高，阻碍了各国提高预期寿命的努力[1]。根据世界卫生组织(WHO)对2019年的估计，在全球183个国家中，有112个国家的70岁以下人群的第一或第二大死因是癌症。此外，还有23个国家将癌症列为第三或第四大死因[2]。在各种癌症中，肝癌已成为许多国家最常见的恶性肿瘤之一。就死亡率而言，它排在第四位(男性第五位，女性第九位)，就发病率而言，它排在第六位(男性第二位，女性第六位) [3]。因此，治疗肝癌对于保障全球公众健康和实施有效的癌症管理措施至关重要。

在目前的临床实践中，针对肿瘤的治疗方法多种多样。这些方法包括放疗、化疗、手术切除、射频消融、冷冻消融和药物治疗。冷冻消融术又称冷冻手术，是一种新型的微创治疗方法。这种技术是在组织中切开一个小切口，然后将冷冻探针从切口插入，到达受影响的区域。然后，冷冻探针与目标组织直接接触，利用极低的温度引起冷冻，进而破坏病变细胞。近年来，冷冻消融术在治疗乳腺癌、肝癌、肺癌、肾癌、前列腺癌和软组织局部恶性肿瘤的过程中，作为手术切除的替代疗法受到了广泛欢迎[4]。

低温消融手术中常用的制冷探头包括液氮、氩气、氧化亚氮和二氧化碳。理论上，这些物质的最低温度可能分别达到−196℃、−187℃、−89.5℃和−78.5℃ [5]。生物组织和细胞对低温的反应可以用阿伦尼斯方程来量化，该方程表明，较低的温度会延长这些生物成分的存活时间。过快或过慢的冷却速度都会导致组织损伤，进而导致组织死亡[6]。此外，复温速度也会对组织损伤程度产生重大影响[7]。学术界广泛研究的重点是确定冷冻组织的最佳时间和温度。人们普遍认为，与组织损伤相关的温度范围通常在−20℃至−40℃之间。一般来说，细胞在−20℃至−30℃之间会明显受到破坏，如果温度降至−40℃至−50℃以下，细胞极有可能受到破坏[8] [9]。组织的低温冷冻会造成血管损伤，破坏循环系统，并最终因冰晶的形成而导致细胞损伤[10]。

因此，在使用低温技术时，既要最大限度地冻结目标组织，又要尽量减少对周围正常组织的伤害，这一点至关重要。在这种情况下，既能彻底清除异常肿瘤细胞，又不会对邻近的健康组织造成重大损害，具有重要意义。事实证明，数值模拟是低温消融过程中温度场分布建模的重要工具，因为现场实验往往不切实际。研究人员利用数值模拟来研究和优化冷冻消融程序。Handler [11]等利用有限单元法对心室组织的冷冻消融程序进行了模拟。Burkov [12]等结合热源边界、人体组织随温度变化的热物理性质以及可变的冷冻手术参数，研究了多探针冷冻消融过程中的温度分布。Jaberzadeh [13]等比较了利用多个探针进行肝脏冷冻消融手术的不同术前计划，他们提出了一种新型模拟算法，旨在优化冷冻范围，同时考虑精确计算冷冻区域。Nabaei和Karimi [14]进行了一项模拟低温消融程序的研究，以调查不同血管大小和位置对肿瘤的影响，结果表明，动脉会对冷冻手术的结果产生重大影响。Nazemian和Nabaei [15]等从CT图像中重建患者特定的肝脏和门静脉网络几何图形，并模拟肝脏肿瘤的冷冻消融。Keangin [16]等建立了一个将电磁波方程、生物热方程和机械变形方程相结合的数学模型，他们使用有限单元法对该模型进行了求解，发现与不包含变形的模型相比，包含变形的模型精度更高。

在冷冻消融手术中，相变现象是关键因素之一。人体组织的含水量通常超过70%，冷冻过程中的相变有可能对组织造成机械损伤。因此，必须考虑相变导致的体积变化以及温度引起的应力。然而，关于冷冻消融过程中的热应力和应变的研究很少。本研究旨在系统研究有变形和无变形数值模型对温度分布和应力分布的影响。所提出的数值模型将生物热方程与机械变形方程相结合，并使用有限元方法进行求解。通过将形变纳入模型，可以全面了解低温消融过程中的温度分布和应力分布。

2. 数值模型的建立及参数的选定

2.1. 几何模型

本研究使用AutoCAD 2023创建了肝组织模型，X、Y和Z方向的尺寸分别约为25 cm、13 cm和10 cm。该模型不考虑肝组织内的血管结构。然后将肝组织模型导入COMSOL Multiphysics 6.0进行模拟仿真，网格划分如图1所示。研究中使用的冷冻探针是一个直径为3 mm的圆柱体，向肝组织内的肿瘤提供冷冻能量。冷冻探针采用高压氮气进行热交换，工作温度为−196℃。模型中还包括导管部分，但导管部分是隔热的。

Figure 1. Geometric model

图1. 几何模型

表1详细列出了低温探针的位置和尺寸以及网格细化区域。

Table 1. Location and size of probe and unit refinement areas

表1. 探针和网格细化区域的位置和尺寸

类型	坐标(中心点)	尺寸
探针	x = −175 mm, y = 70 mm, z = 86 mm	半径r = 1.5 mm, 长度h = 12 mm
网格细化区	x = −175 mm, y = 70 mm, z = 86 mm	半径r = 25 mm

2.2. 参数选定

模拟中相关参数选取于多项相关研究[14] [15] [17]-[20]，具体值见表2。假定在相变的温度范围内，物理特性与温度呈线性关系。

Table 2. Material parameters [14] [15] [17]-[20]

表2. 材料参数[14] [15] [17]-[20]

参数	值	单位	描述
Tu	−1	℃	相变开始温度
Tf	−8	℃	相变结束温度
Kt	0.42 (未冻结)/1.56 (冻结)	Wm^−^1∙℃^−^1	组织导热系数
Ct	3600 (未冻结)/1800 (冻结)	Jkg^−^1∙℃^−^1	组织比热容
Q1	250	MJ∙m^−^3	相变潜热
Qmet	4200 (组织)/42,000 (肿瘤)	W∙m^−^3	代谢热
Wb	0.0005	s^−^1	血液灌注速率
ρb	1.06	kg∙m^−^3	血液密度
Cp,b	3600	Jkg^−^1∙℃^−^1	血液比热容
ρt	1000 (未冻结)/942 (冻结)	kg∙m^−^3	组织密度
Fw	0.711	1	组织水含量
e	0.087	1	组织水相变体应变
μt	0.4 (未冻结)/0.33 (冻结)	1	泊松比
Et	12 (未冻结)/14.8e6 (冻结)	kPa	弹性模量
σs	203 (未冻结)/53,500 (冻结)	kPa	极限强度

3. 理论模型

3.1. 温度模型

Pennes模型是第一个考虑到血流对组织冷冻过程影响的传热模型[21]。该模型在工程热物理学和医学领域都具有开创性意义，因为它能够对生物传热过程、血流的影响以及冷冻手术中的精确治疗进行定量分析。Pennes模型已被广泛应用于临床医学和冷冻手术研究中[12] [14] [15]。在本研究中，肝脏组织被认为是各向同性的，并且在建模中忽略了血管。

瞬态生物热方程有效地描述了热量在肝组织中的传递过程，考虑到新陈代谢热和血液灌注热的Pennes方程如下所示：

$\begin{array}{c}\frac{\partial \left(\rho cT\right)}{\partial t}=\nabla \cdot \left(k\nabla T\right)+{\rho }_b{c}_b{\omega }_b\left(T_b-T\right)+Q_{met}\end{array}$ (1)

其中，$\nabla$ 代表温度空间梯度的拉普拉斯算子。

已知纯水的相变过程发生在0℃。然而，对于非纯净物质，相变的发生是不确定的，可以认为是在一定的温度范围内发生的。尽管如此，一旦发生相变，就会产生潜热，表示在相变过程中吸收或释放能量。与相变相关的热物理问题具有一些特征，如存在两相区以及两相界面存在潜热。因此，相变区域的传热过程是非线性的，给理论分析带来了挑战。从实际角度来看，尽管相变问题具有非线性，但仍必须使用数值方法来解决。在本研究中，生物组织的初始温度设定为37℃。根据组织的冻结状态，将其分为三种不同的状态：未冻结、正在冻结和完全冻结。当温度低于−1℃时，组织开始冻结；当温度达到−8℃时，组织完全冻结。为了计算冷冻过程中的相变潜热，采用了显热容法。这种方法可以将与固相和液相之间转变相关的潜热考虑在内。通过考虑显热容，计算程序可以准确分析热传递现象。显热容法的具体计算公式如下：

$\begin{array}{c}C\left(T\right)=\{\begin{array}{cc}{C}_f;& T<T_f\\ \frac{Q_1}{\rho \left(T_u-T_f\right)}+\frac{C_f+C_u}{2};& T_f\le T\le T_u\\ C_u;& T>T_u\end{array}\end{array}$ (2)

其中，下标f代表冻结状态，下标u代表未冻结状态。

热导率和密度被认为在相变区间内线性变化，冷冻和非冷冻状态下的平均值被作为相变区间的值。具体计算方法如下：

$\begin{array}{c}K\left(T\right)=\{\begin{array}{cc}{K}_f;& T<T_f\\ \frac{K_f+K_u}{2};& T_f\le T\le T_u\\ K_u;& T>T_u\end{array}\end{array}$ (3)

$\begin{array}{c}\rho \left(T\right)=\{\begin{array}{cc}{\rho }_f;& T<T_f\\ \frac{{\rho }_f+{\rho }_u}{2};& T_f\le T\le T_u\\ {\rho }_u;& T>T_u\end{array}\end{array}$ (4)

其中，下标f代表冻结状态，下标u代表未冻结状态。

在这项研究中，假定肝脏组织在温度场中表现出各向同性。此外，在冷冻状态下，肝组织和肿瘤组织内的血液灌注率和代谢热均可忽略不计。本研究中使用的模型参数见表2。

关于冷冻探针，其初始温度设定为37℃，并在一分钟内线性降至−196℃。值得注意的是，只有冷冻探针底部的预定长度被认为在传递冷冻能量时处于活跃状态，而探针的其余部分被视为热绝缘体。

3.2. 力学模型

冷冻消融过程中，组织应力和应变的主要原因是温度的变化。然而，与材料中遇到的传统热应力不同，冷冻消融涉及相变现象。肝脏组织中大量水分的存在约占70% [22]。当发生相变时，例如在冷冻过程中从液态转变为固态，组织内部会发生巨大的体积变化。这些体积变化会导致受周围非冻结区域限制的冻结区域膨胀，从而产生压应力。反之，由于冷冻区域的膨胀，非冷冻区域也会产生压缩应力。因此，相变过程会产生不可忽略的温度变形和温度应力[19]。

Rabin [18]等对冷冻软生物组织在受到外加压应力时的机械反应进行了研究。实验结果表明，冷冻生物组织在压缩应力下表现出一种相当奇特的行为。经过近似线性弹性的初始阶段后，在恒定应变下出现了一系列突然的应力下降。Shi [19]等研究发现，在含水量高的生物材料中的水到冰的转变过程中，体积的大幅膨胀是产生巨大应力的主要原因。此外，在快速冻结的情况下，松弛效应被认为可以忽略不计，因此采用弹性模型足以描述由此产生的机械变形。

本研究使用的低温探针直径为3 mm，长度为12 mm，工作温度为−196℃。冷却过程被认为是快速的，因此需要使用弹性模型来分析组织的机械变形。为了测量人体组织在未冷冻状态下的机械特性，研究者们已经进行了大量研究，并获得了宝贵的数据，这些数据有助于描述人体组织的线性和非线性行为[20]。在本研究中，表2列出了考虑和纳入的模型参数。假设肝脏组织表现出各向同性，三维笛卡尔坐标(x、y和z)的平衡微分方程可表述如下：

$\begin{array}{c}\begin{array}{c}\frac{\partial {\sigma }_x}{\partial x}+\frac{\partial {\tau }_{yx}}{\partial y}+\frac{\partial {\tau }_{zx}}{\partial z}+f_x=0\\ \frac{\partial {\sigma }_y}{\partial y}+\frac{\partial {\tau }_{zy}}{\partial z}+\frac{\partial {\tau }_{xy}}{\partial x}+f_y=0\\ \frac{\partial {\sigma }_x}{\partial z}+\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\tau }_{yz}}{\partial y}+f_z=0\end{array}\}\end{array}$ (5)

其中f代表单位体积表面力，τ代表切向应力。

几何方程为：

$\begin{array}{c}\begin{array}{c}{\epsilon }_x=\frac{\partial u}{\partial x},{\epsilon }_y=\frac{\partial v}{\partial y},{\epsilon }_x=\frac{\partial w}{\partial z},\\ {\gamma }_{yz}=\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z},{\gamma }_{xx}=\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x},{\gamma }_{xy}=\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y},\end{array}\}\end{array}$ (6)

其中ɛ代表应变，ɤ代表剪应变，u，v，w分别代表x，y，z方向的位移。

物理方程为：

$\begin{array}{c}\begin{array}{c}{\epsilon }_x=\frac{1}{E}\left[{\sigma }_x-\mu \left({\sigma }_y+{\sigma }_x\right)\right]+{\epsilon }_x^{th}\\ {\epsilon }_y=\frac{1}{E}\left[{\sigma }_y-\mu \left({\sigma }_x+{\sigma }_x\right)\right]+{\epsilon }_y^{th}\\ {ϵ}_x=\frac{1}{E}\left[{\sigma }_x-\mu \left({\sigma }_x+{\sigma }_y\right)\right]+{\epsilon }_z^{th}\\ {\gamma }_{yx}=\frac{2\left(1+\mu \right)}{E}{\tau }_{yx}\\ {\gamma }_{xx}=\frac{2\left(1+\mu \right)}{E}{\tau }_{xx}\\ {\gamma }_{xy}=\frac{2\left(1+\mu \right)}{E}{\tau }_{xy}\end{array}\}\end{array}$ (7)

其中ɛ^th代表温度应变。其表达式为：

(8)

其中T_R代表初始温度，T代表当前温度，α是热膨胀系数，F_w是组织水含量，e组织水体积应变。F(T)表达式如下：

$\begin{array}{c}F\left(T\right)=\{\begin{array}{l}0,T_u\ge T\hfill \\ \left(T-T_u\right)/\left(T_f-T_u\right),T_f<T<T_u\hfill \\ 1,T\le T_f\hfill \end{array}\end{array}$ (9)

与取弹性模量的方法相同，由于猪肝和人肝的成分相似，所以用猪肝的数据代替了人肝从−8℃到−196℃的温度膨胀系数。Rabin [23]测量了猪肝在−10℃至−180℃区间的热膨胀系数变化。

$\begin{array}{c}\alpha =C_0+C_{1}T+C_2{T}^2+C_3{T}^3\end{array}$ (10)

其中C₀ = 7.323e−5，C₁ = 4.344e−7，C₂ = 6.105e−10，C₃ = 0。

在这项特殊研究中，假定冷冻组织紧贴探针表面。这意味着组织和探针表面接触的位置的所有位移分量(u、v、w)都等于零。

4. 计算结果分析

本研究采用有限元法(FEM)分析与低温消融手术相关的瞬态问题。计算方案包括构建有限元模型，并利用组织特性进行传热计算，以确定局部发热项。为确保精确近似，在相关敏感区域指定了网格细化。使用COMSOL Multiphysics 6.0中的有限元模型解决了生物热和机械变形的耦合方程。这样就能对低温消融过程中肝脏组织内发生的现象进行全面分析。组织内各点的热特性首先被计算出来，并用于求解随时间变化的温度分布。随后，利用这些温度求解机械变形方程，该方程描述了肝组织的机械行为。然后利用获得的肝组织变形重新计算温度场。记录每个时间步的温度，直到达到稳定状态。为了将有限元模型离散化，本研究使用了标准的四面体元素，并使用拉格朗日二次方程来逼近每个元素内温度和应变分布的变化。

本研究利用COMSOL Multiphysics 6.0中的固体力学和固体传热模块进行了耦合计算，肝组织的材料数据如，密度、弹性模量、比热容、导热系数、线膨胀系数等根据−1℃和−8℃两个冻结开始和冻结完成的关键节点设置为分段函数进行赋值。其中，−1℃至−8℃的区间内假设材料特性参数是连续线性变化的，对分段赋值函数进行了线性连续处理。

根据临床数据，模拟了持续900秒的整个冷冻手术过程。本研究进行了网格敏感性测试，以尽量减少数值计算中的尺寸效应。图2显示了这些测试得出的收敛曲线，它展示了距离探针10 mm处的敏感点的温度与时间之间的关系。最终方案使用了121,034个单元。

关于本研究计算结果的准确性方面，由于目前实验条件的限制，生物组织冷冻的实验仅限于离体生物组织的热–力耦合场验证[24] [25]，Rabin和Vispute [24] [25]在冷冻生物保存领域提出了多物理场耦合模型和验证过程。但不同于其研究，本研究是针对于冷冻消融手术过程的仿真，肝组织在体内，温度场方面受到血液灌注热和生物代谢热的影响较大，这是与离体生物组织不同之处，因此存在着伦理道德的限制无法进行实验验证。此外，纯温度场的冷冻手术模拟与验证已有很多研究[14] [15] [26]，这些研究证明了本研究所采用的温度场相关参数在Pennes传热方程下进行冷冻手术温度场模拟的准确性。因此本研究的计算结果虽然缺乏具体实验验证，但仍具有一定的准确性。

Figure 2. Temperature versus time for simulations with different number of cells

图2. 不同单元数量模拟的温度与时间关系

4.1. 温度分布

图3显示了有变形模型中肝脏组织的温度分布，YZ平面的坐标为x = −175 mm。图4显示了无形变模型中的温度分布。在这两个模型中，垂直于探针表面方向的温度显然较低。探针周围的温度分布呈椭圆形，最低温度值出现在探针表面附近，随着与探针距离的增加而逐渐升高。在冻结的初始阶段，温度迅速下降，然后达到峰值，最终接近稳定状态。值得注意的是，在300秒到900秒这段时间内，温度场的变化并不明显。这可归因于几个因素：首先，低温探针与组织的接触面积保持固定，限制了冷冻能量的传递。其次，随着时间的推移，冷却表面的尺寸变大。图3和图4比较了两种模型的温度分布。在没有变形的模型中，见图4，肝组织内的温度分布呈现明显的椭圆形。这是因为探针的圆柱形作为温度源，导致温度分布呈椭圆形。

(a) t = 300 s (b) t = 900 s

Figure 3. YZ-plane temperature distribution (heat-force coupling)

图3. YZ平面温度分布(热–力耦合)

(a) t = 300 s (b) t = 900 s

Figure 4. YZ-plane temperature distribution (pure temperature field)

图4. YZ平面温度分布(纯温度场)

另一方面，在有变形的模型中，见图3，与无变形的模型相比，肝组织内的温度分布呈现出较弱的椭圆形。此外，针头和针尾附近的分布也不对称。这种不对称可能受到温度传导和组织变形的共同影响。首先，由于探针的尖端部分可以进行热交换，而尾部和导管部分是隔热的，因此冷冻探针两侧和尖端方向的温度传递过程最快。其次，肝组织的变形受到探针存在的限制，探针受到周围肝组织的挤压并与之粘合。如前所述，位移边界条件的这一假设已在前面的章节中说明。由于这种限制，肝组织的变形主要发生在探针的侧向和尖端。温度传导和组织变形的共同作用，使变形模型中肝脏组织内部的温度分布更像圆形，且不对称。在这种分布中还可以观察到针尾变细的现象。此外，我们还注意到，在有变形的模型中，探针横向和针向的温度传导速度明显快于无变形的模型。这可归因于低温探针在肝组织内产生的极低温度与肝组织的机械特性之间的相关性。由此产生的组织变形减缓了温度衰减过程，导致温度传导速度慢于无变形模型。

4.2. 应力分布

图5描述了冻结时间为300时YZ平面上的von Mises应力分布。在本研究中，冷冻肝组织的机械参数是以兔肝为基础的，因为兔肝的成分与人类肝脏相似。Rabin [18]研究了冷冻兔肝组织在压缩条件下的力学特性。他们报告的极限应力为53.5 MPa，超过该值，冷冻组织就会出现微裂缝。在当前的研究中，探针表面的von Mises应力似乎超过了53.5 MPa，主要集中在距离探针表面约3 mm的区域。这种局部高应力集中可能是由于在探针表面施加了固定位移边界条件而导致的应力集中。另外，冻结区域的直径约为13 mm。此外，最大von Mises应力出现在针尖处，达到约1000 MPa。这一峰值应力可能是针尖尖锐的几何轮廓造成的应力集中所致。值得注意的是，这里假定屈服应力为53.5 MPa，任何高于53.5 MPa的应力都是不现实的，因为组织中形成的微裂缝会释放如此高的应力。因此，在未来的研究中考虑到塑性变形或微裂纹形成的应力模型对于冷冻手术中的精确应力建模是必要的。

Figure 5. YZ-plane stress distribution

图5. YZ平面应力分布

5. 结论

本研究使用COMSOL Multiphysics 6.0建立了多物理场有限元热应力模型，以模拟低温消融过程。 对人体肝脏组织进行了实体建模，同时考虑了组织的变形。对有变形和无变形的模型进行了比较，以分析温度场的分布。得出了以下结论：

1) 在是否考虑热应力和热应变的两个模型中，最低温度值都出现在探针表面附近，温度值随着与探针距离的增加而逐渐升高。在冻结的初始阶段，温度迅速下降，然后达到峰值，最终接近稳定状态。

2) 在没有变形的模型中，肝组织内的温度分布呈现明显的椭圆形。而在有变形的模型中肝组织的温度分布由于热应力和热应变的作用而不呈椭圆形分布。

3) 在有变形的模型中，探针横向和针向的温度传导速度由于变形的影响而明显快于无变形的模型。

在有变形探针表面的von Mises应力超过了假定屈服应力53.5 MPa，主要集中在距离探针表面约3 mm的区域。

4) 在有变形的模型中，肝脏组织表面由于探针和组织冻结后的约束作用呈现出不均匀的变形模式。

这表明了在模拟中考虑组织变形的重要性。这些发现为低温消融手术的实际治疗提供了宝贵的见解。此外，有必要进行进一步的实验来研究组织在冷冻状态下的机械特性，尤其是其塑性特性。了解冷冻组织的机械行为对于开发更精确的应力模型和改进冷冻手术的整体模拟至关重要。

NOTES

^*第一作者。