解海衛(wèi) 張 艷 諸 凱

（天津商業(yè)大學(xué)機(jī)械工程學(xué)院，天津市制冷技術(shù)重點(diǎn)實(shí)驗(yàn)室，天津 300134）

引言

模擬生物溫度場(chǎng)用于輔助分析已經(jīng)成為生物醫(yī)學(xué)工程研究中的常用方法。血液與生物組織有復(fù)雜的傳熱耦合關(guān)系，準(zhǔn)確地分析血液對(duì)流換熱是模擬生物組織傳熱的必要前提。

許多研究者為分析血液對(duì)流換熱做了大量的實(shí)驗(yàn)研究和理論計(jì)算。Charm進(jìn)行了實(shí)驗(yàn)測(cè)量，并分析了血液在不銹鋼毛細(xì)管（直徑0.6mm）內(nèi)的熱交換［1］。Shah通過(guò)測(cè)量一定加熱量P下對(duì)應(yīng)的溫升T，由已知經(jīng)典的努賽爾數(shù)（Nu數(shù)）和P/T的擬合公式，確定了單管內(nèi)血液的Nu數(shù)［2］。Victor對(duì)單個(gè)血管的血液流動(dòng)和換熱進(jìn)行了模擬，計(jì)算得到了血液Nu數(shù)沿管長(zhǎng)的變化，其中壁面條件是均勻壁溫和均勻熱流兩種，入口為完全發(fā)展流動(dòng)［3］，入口條件是均勻流速、溫度兩種［4］。Luisa假設(shè)血液為牛頓流體，根據(jù)動(dòng)量和能量守恒，得出了單根血管局部加熱時(shí)的對(duì)流換熱系數(shù)公式［5］。Tungjitkusolmun假設(shè)血管壁溫為常數(shù)，分析計(jì)算了熱療時(shí)大血管對(duì)生物組織溫度場(chǎng)的影響［6］。

但是，這些研究均是以單個(gè)血管為對(duì)象，著重用不同方法、確定不同換熱條件下的血液Nu數(shù)（或h），而考慮在循環(huán)系統(tǒng)實(shí)際分支結(jié)構(gòu)的基礎(chǔ)上，分析血液換熱的文獻(xiàn)鮮見(jiàn)報(bào)道。Anil研究表明，血管的分支結(jié)構(gòu)對(duì)血液速度場(chǎng)分布有明顯影響，因此分支結(jié)構(gòu)對(duì)血液熱交換也一定存在影響［7］。根據(jù)傳熱學(xué)理論，流速和管徑也是影響管內(nèi)對(duì)流換熱的主要因素。筆者通過(guò)建立血管樹(shù)簡(jiǎn)化模型，模擬計(jì)算了分支結(jié)構(gòu)內(nèi)血液的三維速度場(chǎng)、溫度場(chǎng)和血液Nu數(shù)；從血管樹(shù)的實(shí)際結(jié)構(gòu)出發(fā)，分析了血管分支結(jié)構(gòu)參數(shù)（半徑比、分支角）、血液流速、血管半徑對(duì)血液與組織之間熱交換的影響。

1 血管分支結(jié)構(gòu)模型

研究表明，血液和組織的熱交換主要發(fā)生在動(dòng)脈和小動(dòng)脈處（R＞0.1mm），而這些動(dòng)脈一般是以上級(jí)血管的分支形式存在，且末端為二元分支結(jié)構(gòu)［8］，因此計(jì)算血液換熱應(yīng)在分支結(jié)構(gòu)中分析。建立循環(huán)系統(tǒng)分支結(jié)構(gòu)模型，如圖1所示，模型由一個(gè)主干血管和兩個(gè)分支血管組成。根據(jù)血管樹(shù)生成理論，由優(yōu)化原則建立血管樹(shù)，其主干血管和分支血管的軸線在同一平面內(nèi)［9］，即該血管樹(shù)關(guān)于此平面對(duì)稱(chēng)，因此截取分支結(jié)構(gòu)的一半作為最終建立的血管樹(shù)物理模型。

選取主干血管半徑為0.7mm，屬于熱重要血管的范圍［10］，長(zhǎng)度為20mm。分支Ⅰ與分支Ⅱ的長(zhǎng)度為40mm；分支Ⅰ與主干血管之間的夾角（即分支角α（定義見(jiàn)圖1），分別為90°、120°、150°、180°，分支血管Ⅰ的半徑與主干血管的半徑之比（即半徑比）分別為0.75、0.8、0.85、0.9；分支血管Ⅱ的半徑由Murray定律確定。

圖1 血管分支結(jié)構(gòu)Fig.1 Schematic of the furcated configuration

2 數(shù)學(xué)模型的建立

2.1 模型假設(shè)［3，11］

忽略血管的膨脹性，將流動(dòng)視為穩(wěn)態(tài)、層流；血液為常物性，密度、比熱不隨溫度變化；忽略血液的代謝熱；忽略近壁面處血漿層的影響；血液為不可壓非牛頓流體，其非牛頓性用冪律模型描述，有

式中，γ為切變率，1/s；τ為切應(yīng)力，Pa；k為常數(shù)，取0.014 67（Pa·sn）；n為非牛頓指數(shù)，取0.775 5［12］。

2.2 控制方程

計(jì)算血液與組織的對(duì)流換熱，需耦合求解血液的流場(chǎng)和溫度場(chǎng)，即求解以下控制方程為

質(zhì)量方程（連續(xù)性方程）

(1)看皮下脂肪層的厚度。“瘦肉精”豬皮下脂肪層明顯變薄，通常不足1cm;而正常豬肉在皮層和瘦肉之間會(huì)有一層脂肪，肥膘約為1～2cm厚。因此，遇到脂肪層太薄、太松軟的豬肉要格外小心。

動(dòng)量方程

能量方程

式中，U為流體速度，m/s；S為廣義源項(xiàng)，Pa/m；Ψ為黏性耗散函數(shù)，1/s2；ρ為密度，kg/m3；P為壓力，Pa；μ為動(dòng)力黏度，Pa·s；T為溫度，K；c為定壓比熱，J/（kg·K）；λ為導(dǎo)熱系數(shù)，W/（m·K）。

血液各參數(shù)取值見(jiàn)表1。

表1 血液熱物性參數(shù)取值［13］Tab.1 Properties of blood

2.3 邊界條件

1）主干血管入口條件：假設(shè)主干血管入口邊界條件為均勻速度和均勻溫度，Tin=310K。

2）分支血管出口條件：Pout=0（相對(duì)壓力），qout=0。

3）血管軸線對(duì)稱(chēng)面條件：對(duì)稱(chēng)面絕熱，即q=0；流動(dòng)邊界條件為在垂直于對(duì)稱(chēng)面的方向上流速為0。

4）血管壁面條件：壁面無(wú)滑移，即Uw=0；熱邊界條件為均勻壁溫條件，Tin=309K。

2.4 網(wǎng)格劃分和數(shù)值模擬

圖2 有限元網(wǎng)格。（a）直管部分網(wǎng)格；（b）分支處網(wǎng)格；（c）血管截面網(wǎng)格Fig.2 Finite element mesh of the vascular tree.（a）mesh on straight part of vessel；（b）mesh on crosssection；（c）mesh at the bifurcation point

利用ANSYS有限元軟件中的FLOTRAN模塊，分析血管分支結(jié)構(gòu)內(nèi)的速度場(chǎng)和溫度場(chǎng)，得到血管壁面處血液與壁面的換熱量，進(jìn)而根據(jù)努塞爾數(shù)定義，計(jì)算得出分支血管Ⅰ管長(zhǎng)x處（如圖1所示）血管壁面任意位置ф處的血液努塞爾數(shù)Nux，Ф數(shù)，以及x處的截面平均努賽爾數(shù)Nux。

3 分支結(jié)構(gòu)內(nèi)的速度場(chǎng)和溫度場(chǎng)

圖3表示血管連接處速度和溫度分布，其中（a）和（b）分別為模擬得到的主干血管血液流速為80mm/s、主干血管半徑為0.7mm、分支角為120°、分支Ⅰ半徑比為0.8時(shí)，血管連接區(qū)域的速度和溫度分布等值線圖。由圖3（a）可以看出，血液從主干血管流入分支血管時(shí)，速度減小，流向改變，截面高流速區(qū)向一側(cè)偏移，速度場(chǎng)由主干血管的軸對(duì)稱(chēng)分布變?yōu)榉种а苋肟谔幍姆禽S對(duì)稱(chēng)分布，然后逐漸恢復(fù)為軸對(duì)稱(chēng)分布狀態(tài)；由圖3（b）可以看出，分支血管入口處的溫度分布也呈現(xiàn)與速度場(chǎng)類(lèi)似的變化，原主干血管軸心處的低溫區(qū)，隨流動(dòng)偏移到分支血管位置1一側(cè)（見(jiàn)圖2）。此外，兩個(gè)分支血管的速度場(chǎng)和溫度場(chǎng)的偏心程度不同，這是由于兩個(gè)分支血管的分支角度和半徑比不同而引起的。

從圖3還可以看出，除主干血管出口區(qū)域速度場(chǎng)有一個(gè)較小的擾動(dòng)外，兩個(gè)分支血管的存在對(duì)主干血管的速度場(chǎng)和溫度場(chǎng)基本沒(méi)有影響，主干血管的速度和溫度等值線基本呈平行狀態(tài)進(jìn)入血管分叉區(qū)域，分支對(duì)主干血管內(nèi)流動(dòng)換熱的影響可以忽略，即在循環(huán)系統(tǒng)中，血管主要受上級(jí)來(lái)流主干血管的影響，而其下級(jí)分支血管的影響可以忽略。由此可見(jiàn)，在模型中，分支血管Ⅰ、Ⅱ內(nèi)血液換熱主要受主干血管的影響，而分支Ⅰ、Ⅱ的下級(jí)分支對(duì)其影響可忽略，所以血管樹(shù)模型無(wú)需建立分支Ⅰ、Ⅱ的下級(jí)分支。

圖3 血管連接處速度和溫度分布。（a）速度場(chǎng)分布等值線（單位：m/s）；（b）溫度場(chǎng)分布等值線（單位：℃）Fig.3 Contour plot of the numerical solutions at the bifurcation position of the vessel system.（a）contour plot of velocity distribution（Unit：m/s）；（b）contour plot of temperature distribution（Unit：℃）

4 分支血管內(nèi)Nu數(shù)變化規(guī)律與分析

由上述可知，分支結(jié)構(gòu)的存在會(huì)影響分支血管內(nèi)的速度場(chǎng)和溫度場(chǎng)，使其在入口處呈現(xiàn)非軸對(duì)稱(chēng)分布，相應(yīng)地也會(huì)影響分支血管內(nèi)的Nu數(shù)，使之呈現(xiàn)與單個(gè)直管內(nèi)Nu數(shù)不同的分布規(guī)律。因此，筆者以分支血管Ⅰ為研究對(duì)象，計(jì)算了血管分支結(jié)構(gòu)中分支血管Ⅰ內(nèi)任意截面x（即管長(zhǎng)x處）的平均努賽爾數(shù)Nux和徑向位置1、2、3處（如圖2所示）的局部努賽爾數(shù)（Nux，1、Nux，2、Nux，3），并分析了分支結(jié)構(gòu)參數(shù)（半徑比、分支角）、分支血液流速、分支血管半徑的變化對(duì)截面平均努賽爾數(shù)Nux的影響。

4.1 分支血管內(nèi)努賽爾數(shù)分布規(guī)律

圖4為分支血管Ⅰ半徑為0.56mm、平均流速為30mm/s，分支角為120°，半徑比為0.8時(shí)，其內(nèi)截面平均Nux與相同管徑、流速下單個(gè)直管內(nèi)截面平均Nux的對(duì)比曲線。從該圖可以看出，在血管入口段，分支血管的Nux明顯小于單個(gè)直管的Nux，且收斂速度大于直管。這是因?yàn)檠毫鹘?jīng)主干血管后，由入口處的速度、溫度均勻的流體變?yōu)槌隹谔幩俣?、溫度完全發(fā)展的流體，雖然流入分支后，由于方向的改變?cè)斐伤俣葓?chǎng)和溫度場(chǎng)的變化，但相對(duì)于流入單個(gè)直管的血液，更接近于熱完全發(fā)展階段，因此Nux較小，且收斂速度快。在血管出口處，分支血管和單個(gè)直管的Nux都達(dá)到收斂值，約3.75。

在上述條件下，分支Ⅰ內(nèi)3個(gè)典型圓周位置處的壁面局部Nux，Ф（Ф=1，2，3）隨x的變化曲線?？梢钥闯觯涸诜种а苋肟诙危煌瑘A周位置之間Nux，Ф差異較大，這是由于血液流入分支血管時(shí)，流向改變引起溫度出現(xiàn)非軸對(duì)稱(chēng)分布，不同圓周位置處的溫度梯度不同，使Nux，Ф出現(xiàn)差異；而在出口截面段，3個(gè)圓周位置處的Nux，Ф基本相同，約為3.75，這與分支血管內(nèi)溫度場(chǎng)由非軸對(duì)稱(chēng)分布恢復(fù)為軸對(duì)稱(chēng)分布的變化過(guò)程一致。

圖4 分支血管與單個(gè)直管截面平均Nux數(shù)的對(duì)比曲線Fig.4 The comparison between the mean Nuxat daughter vessel I and at a straight vessel

圖5 不同圓周位置處的局部Nux，Ф數(shù)隨管長(zhǎng)變化的曲線Fig.5 The local Nux，Ф at three different positions 1，2 and 3 in daughter vessel I versusx

4.2 分支結(jié)構(gòu)參數(shù)對(duì)截面平均努賽爾數(shù)的影響

當(dāng)分支角分別是90°、120°、150°、180°時(shí)，其他參數(shù)與4.1節(jié)條件相同，分支Ⅰ內(nèi)Nux隨x變化的曲線如圖6所示；當(dāng)血管半徑比分別是0.75、0.8、0.85、0.9時(shí)，分支Ⅰ內(nèi)Nux隨x變化的曲線如圖7所示。由這兩圖可以看出：分支角、半徑比越小，入口段截面平均Nux越大；分支角、半徑比變化對(duì)分支血管截面平均Nux的收斂速度（熱入口段長(zhǎng)度）沒(méi)有影響，收斂值都約為3.75。這主要是由于分支角、半徑比越小，分支血管入口處速度場(chǎng)和溫度場(chǎng)的偏心程度越嚴(yán)重，相應(yīng)Nux也就會(huì)越大。

4.3 分支血管管徑、流速對(duì)Nux的影響

當(dāng)分支血管半徑為0.48、0.56、0.64、0.8mm時(shí)，或分支血管血液流速為20、30、40、50mm/s時(shí)，其他參數(shù)與4.1節(jié)條件相同，分支Ⅰ內(nèi)Nux隨x變化的曲線如圖8和圖9所示?？梢钥闯觯悍种а馨霃?、血液流速越大，入口段截面平均Nux越大；血管半徑和血液流速的變化對(duì)分支血管截面平均Nux的收斂速度也有影響，血管半徑或血液流速越小，Nux的收斂速度越快，收斂值都約為3.75。這主要是由于血管半徑或流速越小，血液流動(dòng)中管壁對(duì)流體的黏滯力與血流的慣性力之比越大，血流從入口處的偏心流動(dòng)恢復(fù)為軸對(duì)稱(chēng)流動(dòng)的速度越快，相應(yīng)熱入口段的長(zhǎng)度也越短，Nux的收斂速度越快。

圖6 不同分支角度對(duì)應(yīng)的Nux隨x變化的曲線Fig.6 The mean Nuxwith different furcating angels

圖7 不同半徑比對(duì)應(yīng)的Nux隨x變化的曲線Fig.7 The mean Nuxwith different bifurcation ratios

圖8 不同分支血管半徑對(duì)應(yīng)的Nux隨x變化的曲線Fig.8 The mean Nuxwith different vessel radii

圖9 不同分支血管流速對(duì)應(yīng)的Nux隨x變化的曲線Fig.9 The mean Nuxwith different blood velocities

5 結(jié)論

以反映循環(huán)系統(tǒng)實(shí)際結(jié)構(gòu)的血管分支物理模型為平臺(tái)，用有限元方法數(shù)值模擬了血液的三維速度場(chǎng)、溫度場(chǎng)，得到了分支血管中典型位置處Nux，Ф和截面平均Nux，分析了血管半徑、血液流速和血管樹(shù)結(jié)構(gòu)參數(shù)（半徑比、分支角）對(duì)截面平均Nux的影響，結(jié)果表明：

1）當(dāng)血液由主干血管流入分支血管時(shí)，血管截面速度場(chǎng)和溫度場(chǎng)由主干血管出口的軸對(duì)稱(chēng)分布變?yōu)榉种а苋肟谔幍姆禽S對(duì)稱(chēng)分布，然后隨著在分支血管內(nèi)的流動(dòng)又逐漸恢復(fù)為軸對(duì)稱(chēng)狀態(tài)，而分支血管對(duì)主干血管內(nèi)流動(dòng)和換熱的影響較小。

2）分支血管入口段的Nux明顯小于單個(gè)直管的Nux，且收斂速度大于直管，收斂值均約為3.75。分支血管入口處血液的Nux，Ф為明顯非軸對(duì)稱(chēng)分布，且隨血液沿血管繼續(xù)流動(dòng)，其不對(duì)稱(chēng)程度逐漸變小。

3）分支角越小、半徑比越小、分支血管半徑越大、分支血管流速越大時(shí)，入口段截面平均Nux越大。各參數(shù)對(duì)入口段截面平均Nux的影響程度為：血液流速＞半徑比，血管半徑＞分支角。分支血管半徑、流速越小，截面平均Nux的收斂速度越快。

［1］Charm S，Paltiel B，Kurland GS.Heat transfer coefficients in blood flow［J］.Biorheology，1968，5（1）：133-145.

［2］Shah J，Dos SI，Haemmerich D，et al.Instrument to measure the heat convection coefficient on the endothelial surface of arteries and veins［J］.Medical＆ Biological Engineering and Computing，2005，43（4）：522-527.

［3］Victor SA，Shah VL.Steady state heat transfer to blood flowing in the entrance region of a tube［J］.International Journal of Heat and Mass Transfer，1976，19（6）：777-783.

［4］Victor SA，Shah VL.Heat transfer to blood flowing in a tube［J］.Biorheology，1975，123（3）：361-368.

［5］Luisa C，Dos SI，Haemmerich D.Theoretical analysis of the heat convection coefficient in large vessels and the significance for thermal ablative therapies［J］.Physics in Medicine and Biology，2003，48（12）：4125-4134.

［6］Tungjitkusolmun S，Staelin T，Haemmerich D，et al.Threedimensional finite-element analysis for radio-frequency hepatic tumor ablation［J］.IEEE Transactions on Biomedical Engineering，2002，49（1）：3-9.

［7］Anil K，Varshney CL，Sharma GC.Computational technique for flow in blood vessels with porous effects［J］.Applied Mathematics and Mechanics，2005，26（1）：63-72.

［8］Marek K，Yan R，Johanne B，et al.Physiologically based modeling of 3-D vascular networks and CT scan angiography［J］.IEEE Transactions on Medical Imaging，2003，22（2）：248-257.

［9］張艷，張于峰，諸凱，等.三維動(dòng)脈血管樹(shù)的模擬與分析［J］.天津大學(xué)學(xué)報(bào)，2007，40（9）：1071-1076.

［10］Chato JC.Heat transfer to blood vessels［J］.Journal of Biomechanical Engineering，1980，102（5）：110-118.

［11］Barozzi GS，Dumas A.Convective heat transfer coefficients in the circulation［J］.Journal Biomechanical Engineering，1991，113（3）：308-313.

［12］Walburn FJ，Schneck DJ.A constitutive equation for whole human blood［J］.Biorheology，1976，13（1）：201-214.

［13］陳琦，白景峰，陳亞珠.伴行血管在高強(qiáng)度聚焦超聲下對(duì)溫度場(chǎng)的影響［J］.上海交通大學(xué)學(xué)報(bào)，2004，38（1）：130-134.