液氮冷刀共軛傳熱模型消融腫瘤的研究

公維芳，曹秒，王彬

（長春理工大學 生命科學技術學院，長春 130022）

惡性腫瘤是全球第二大死因，在未來的二十年中全世界癌癥病例數可能會增加60%[1]。盡管外科手術是腫瘤治療的首選方法，但是大部分病人在明確診斷時腫瘤進展已處于中晚期，能獲得手術機會的病人約為20%～30%。近年來，臨床上逐漸開始應用新興的腫瘤療法[2]，低溫冷凍手術作為一種微創物理療法，具有創傷小、恢復快等優點[3]，使得一些不耐受手術的多發性腫瘤轉移的患者獲得延長生存期甚至獲得根治腫瘤的機會。

冷凍消融手術兼顧創傷小和治療效果是非常困難的，因為監測凍結過程這一技術難題限制了冷凍消融手術在臨床應用上的推廣，因此模擬研究就顯得格外重要。Kumar等人[4]通過雙曲線和拋物線型兩種模型模擬了肺癌手術數值仿真過程中的相變傳熱過程。舒春等人[5]采用浸入邊界法來研究冷凍過程中血管的熱效應，并將仿真結果與已發表的結果對比驗證，結果吻合。CHAPAL S M等人[6]提出了瞬態三維兩相制冷模型來評估前列腺手術最佳冷凍效果，通過數值計算研究了組織內特定點的溫度分布。GOLKAR E[7]提出基于GPU的三維冰球快速建模方法預測冰球直徑，仿真實驗驗證結果是誤差僅在冰球直徑的5.8。SHARMA等人[8]仿真過程中加入導熱系數低的全氟乙烷以達到快速凍結的目的。本文采用液氮為一相，液氮蒸發氣化為二相的雙相流共軛傳熱模型對冷刀的物理行為進行建模，并在氣泡填充法基礎上進行了計算優化，此模型可用于分析冷刀插入深度、結構參數及放置位置對消融結果的影響。

1 求解生物傳熱方程

冷刀共軛傳熱模型為應用于微小圓管通道的流體模型，即流體流經一個空間位置固定的微團，因此其連續性方程可以寫為如下形式：

其中，ρm為液氮混合相密度；v為液氮平均流速；ρ1為液氮一相密度；φ1為一相體積分數；ρ2為二相密度；φ2為二相體積分數。

生物組織間的耦合傳熱可由Pennes生物傳熱方程描述[9]：

其中，T為組織溫度；ρ為密度；cb為血液的比熱；為血液灌注率；體積代謝熱（腫瘤：，正常組織：）；下標b表示血液。

在組織凍結過程中，相變往往發生在很小的時間范圍（Tl，Ts）內，其中上轉變溫度Tl=1 ℃，下轉變溫度Ts=-8℃，引入比熱和導熱系數方程，作為腫瘤內部熱源：

其中，Cs=1977J/(kg·℃)是凍結組織的比熱；Cl=3688J/(kg·℃)是未凍結組織的比熱；Ks=1W/（m·℃）是凍結組織的導熱系數；Kl=0.42 W/（m·℃）是未凍結組織的導熱系數，Q=250 MJ/Kg是組織潛熱。

2 共軛傳熱模型建立及數值計算過程

液氮冷凍刀內部的流體流動和傳熱是不穩定的瞬態相變現象，從而使冷凍過程無法控制，無法實時精確監測冷凍效果，而數值計算和仿真研究可以在計算機上通過設定條件來模擬實驗所需環境和實驗過程，精確的物理模型可以幫助研究者通過冷刀和生物組織的相互作用來分析和監測冷凍消融手術過程中的溫度分布和消融體積變化情況。

單個冷凍刀僅可完全摧毀徑向長度小于2 cm的腫瘤，當腫瘤大于等于3 cm時就需要多個冷刀相互協助消融[10]。而三個直徑為3 mm的冷刀可在相對較短的冷凍時間內殺死腫瘤細胞，并盡可能減少對健康細胞的損傷[11]。建立的模型如圖1所示，包括圓柱形健康組織，球形腫瘤組織，以及由內外同心管組成的冷刀，冷刀采用不銹鋼材質內外管底端距離5 mm，內管直徑為可調節的，外管直徑為3 mm，絕熱段42 mm傳熱段8 mm，其中組織頂部暴露在空氣中進行對流換熱，模型周圍設定為恒定體溫37℃。為了研究溫度隨時間變化的趨勢，取P2、P3、P4點分別距離冷刀5 mm、15 mm、10 mm，在仿真過程中記錄測溫點的溫度變化，用來評估冷刀的制冷效果。

圖1 腫瘤消融三維模型和冷刀結構模型

采用有限元離散軟件Fluent對液氮冷刀模型進行數值仿真計算，采用GAMBIT對腫瘤冷刀模型網格劃分[12]，由于網格疏密對數值計算結果有一定的影響，計算結果的精確度隨網格數量的增加而提高，與此同時會造成計算量增大、計算周期增長的問題。因此做以下兩點處理：首先進行網格無關性驗證，所謂網格無關性即在計算精度和計算開銷間尋求合適的點，這個點所處的位置就達到了網格無關性的閾值；其次在需要著重研究的冷刀傳熱段采用較密的劃分，腫瘤區域次之，健康組織采用較稀的網格數量。如表1所示，網格數量60 W和80 W消融體積差小于1%，結合計算效率，選取60 W網格對應數值來進行仿真計算。

表1 網格無關性驗證

求解過程選擇基于壓力的瞬態求解模式，并應用Eulerian模型來計算雙相流問題，選擇能量方程作為控制方程，冷刀入口處雷諾數為20 000，K-epsilon Realizable黏度模型實現對液氮冷刀的數值求解。該模型的湍動能及耗散率運輸方程為：

其中，Gk表示由于平均速度梯度引起的湍動能產生；Gb是浮力引起的湍動能產生；YM是可壓縮湍流脈動膨脹對總的耗散率的影響；C2和C1ε是常數，分別為1.9和1.44；σk和σε分別是湍動能及其耗散率的湍流普朗特數，作為默認值常數分別為1.0和1.2。

通過將數值計算的預測結果與張鑫等人[13]公布的實驗結果進行比較，驗證了兩相流共軛傳熱冷刀模型預測生物組織內部溫度分布的有效性。

圖2 測溫點溫度對比圖

3 氣泡填充法的優化

氣泡填充法是一種基于計算機的規劃方法，用于在指定區域內使一定數量的氣泡受力不斷運動直至均勻分布的物理方法[14]。該方法首先在給定區域內形成球形氣泡，假設氣泡間受到類似范德華力的作用，通過施加物理意義上的張弛找到緊密的力平衡的氣泡排列。此外沿著目標區域的輪廓還放置了一組靜止的氣泡作為邊界氣泡。根據文獻[15]，簡化的范德華力模型如下：

邊界條件如下：

其中，l是兩個相互作用的氣泡間的實際距離；l0是相互作用的氣泡間的穩定距離；k0是距離為l0時的線性常數；α、β、γ、ε是簡化范德華力模型的系數。

在目前研究中，所有氣泡的大小和形狀在趨于平衡的過程中是保持不變的[16]，氣泡的尺寸是根據目標區域的大小、正在使用的低溫探針數量來選擇的。因此，氣泡體積如下：

其中，Vp是氣泡的體積大小；Vt是目標區域的體積大小；n為此次消融的氣泡數量，即對應的冷刀數量；RP是氣泡軸向半徑；1.596為此模型通過反復實驗得出，可用于所提出的肝部腫瘤模型。由于肝部腫瘤模型內無特殊熱源，所以冷刀橫向布局應為均勻分布，即θ=60°。若想達到最佳消融效果，既要得到最大的冰球覆蓋面積又要保全健康組織，即令冷凍治療范圍超過腫瘤2 mm左右。已知腫瘤半徑Rt，所以根據三角函數公式計算得到冷刀間距d：

4 實驗結果與分析

4.1 冷刀插入深度及長徑比影響分析

冷凍手術的第一步就是要確定冷刀插入深度，由實驗可得冷刀在冷凍過程中形成的消融體積是橢球形，其中單個冷刀消融的軸向和徑向直徑可達到9 mm和11.6 mm。隨著消融時間的增加，冷凍范圍逐漸逼近健康組織，所以確定插入深度可以有效減少健康組織的損傷。由圖2可以直觀看出，傳熱段8 mm的冷刀插入腫瘤10 mm在消融至300 s時，冷凍范圍距健康組織過近會導致組織不可逆的損傷，與插入腫瘤12 mm和14 mm模型對比可知插入深度14 mm冷凍范圍居腫瘤中心，為適宜的插入深度。

圖3 冷刀插入深度溫度分布圖

冷刀內管是液氮流入的通道，管徑的大小影響冷刀頭處的液氮流速，而液氮流速的改變進一步影響組織間的換熱過程和組織中產生的溫度分布，因此做了以下三組模型進行對比分析。圖3為直徑不同的三組冷刀模型流速矢量圖，由圖可以看出內管直徑越小形成的湍流窩越小，液氮流下撞擊管壁造成冷刀震蕩越強烈，比較后得出內管直徑1.8 mm的冷刀模型速度矢量圖形成的湍流較為規則。消融情況由圖4驗證可得，內管直徑1.8 mm的冷刀模型消融體積最大，300 s時單冷刀消融體積達到0.681 5 cm3，而內管直徑繼續增加會對回流產生巨大阻礙，因此選用1.8 mm的內管直徑做進一步的研究。

圖4 不同內管直徑的冷刀頭液氮流速矢量圖

4.2 冷刀布局對肝腫瘤的消融結果

由式（10）計算得出最佳冷刀間距為14 mm，圖5是其在-22℃下消融體積隨時間變化圖。由圖5知150 s時模型熱耦合區域（中心區域）未見明顯空隙，這是冷量疊加的效果，即在一定距離范圍內，冷刀間存在協同作用。對間距為13 mm和15 mm的冷刀布局模型仿真計算做進一步驗證，其消融體積及各測溫點溫度如圖6所示，間距為14 mm的冷刀模型消融效果最佳。冷刀間距越小，中心熱耦合區域溫度越低，隨冷刀間距增大，近冷刀端p3溫度越低，而遠冷刀端p4溫度差異不大。實驗表明隨著治療時間的增加，冷刀間距增大會加重對健康組織的損傷，而冷刀間距為14 mm模型的遠冷刀端溫度最低，在得到最大消融體積的同時滿足對健康組織損傷較小的要求。

圖5 不同內管直徑模型300 s消融體積圖

圖6 冷刀間距為14 mm的模型-22℃等溫線體積圖

圖7 二維模型消融體積及溫度變化曲線

5 結論

基于有限元數值仿真方法建立了雙相流共軛傳熱模型，并應用了改進的氣泡填充法，分析了三個液氮冷刀應用于肝臟腫瘤的三組模型的消融結果并對其進行比較。隨著消融時間的增加，可形成橢球形消融區域，為減少健康組織損傷并保證消融效果，計算并比較了不同的冷刀插入深度、內管直徑和冷刀間距在目標區域及健康組織的溫度分布情況。最終得出插入深度為14 mm，冷刀內管直徑為1.8 mm，間距為14 mm時的多冷刀布局具有很好的消融效果。