針對肝移植術前清洗器官的灌注環節所引發的血管內流體動力學行為進行研究, 為術前清洗的相關操作提供理論指導。構建了帶有異物的一級直血管、彎曲血管實體模型和血管內液體湍流的控制方程。實驗測定醫用灌注液的物理參數, 提出了估算灌注技術參數的方法, 實現了不同條件下的灌注仿真。基于構建的流體控制方程和血管模型進行計算, 其初步結果符合醫療現場實際操作值, 結果較為清晰地顯示了灌注過程中(尤其是管內異物周圍)的流體動力學行為, 如回流滯留現象等。同時, 結果也顯示了血管內各處流速場大小和管壁面壓強值分布。隨著入口初始灌注速度的增加, 兩種類型的血管模型內的壓強值和流速場均增大, 在異物處產生壁透壓和回流滯留區, 同時由于彎曲血管形狀的影響, 其內部流體動力學行為比直血管更加復雜。
引用本文： 劉軍, 范勇, 劉懿禾. 簡易肝血管模型灌注過程中的流體動力學行為初探. 生物醫學工程學雜志, 2016, 33(2): 260-267. doi: 10.7507/1001-5515.20160045 復制
引言
目前國內外不論是發達地區還是欠發達地區肝病誘發的死亡率都很高，肝移植是而目前治療終末期肝病的唯一有效手段。肝供體質量高低是影響肝移植手術成功率的關鍵環節之一，術前對肝供體的清洗(灌注)是必須操作的醫療過程。對肝臟內部血管的流體動力學特性，尤其是灌注過程的動力學行為的研究，有助于提高肝臟供體的灌注質量和對肝供體的保護。對血管建模并進行相應的力學分析，將有助于降低肝供體灌注過程中的主觀因素對清洗質量的影響，同時為規范操作提供理論依據。
劉向明等[1]對生物醫學中的數學建模問題進行了綜述，闡明了其分類和各自特點以及研究目標。王振全等[2]概述了目前國內肝臟血流動力學臨床應用的研究進展。肝臟血管最多可向下分至8級，而目前，國內外對于肝臟血管機制的研究最多能達到前4級。為了構建肝臟血管的結構模型，在血管組織及其結構的可視化研究基礎上，國內外大多利用超聲成像[3-5]、核磁共振[6-7]和計算機斷層成像(computed tomography, CT)[8-10]等技術來獲得血管系統的二維或三維實際圖像數據，然后將數據導入三維造型軟件中進行肝臟血管模型的建立。也有研究者直接利用三維造型技術進行血管系統的模型構建，將重點放在血管和血液的流-固耦合及其相互影響和作用的研究上。近年來，研究者也提出了一些新的原理和方法，如Hoganson等[11]應用仿生原理進行肝臟血管模型的優化及設計，開發了肝臟血管的實體模型。Matsumoto等[12]利用二維灰階血流成像技術(two-dimensional blood flow, B-Flow)對肝臟血管進行圖像分析，并將該技術與彩色多普勒超聲技術進行比較，探討了二者的異同點和優缺點。
對肝移植手術中的醫學問題進行定量分析涉及有關工學領域的知識，如肝臟血管數學模型、傳感器檢測和血流動力學計算等方面。國內外對肝移植的研究主要是建立在核磁共振或CT技術提供的肝臟血流分布信息上，并依此來進行相關研究。Debbaut等[13]利用CT技術對肝臟血管進行圖像建模，使用文獻數據模擬肝臟內低溫機械灌注情況下的液體流動狀況，得出灌注壓力及流體黏度對于血管壁面剪切應力(wall shear stress, WSS)的影響。但并未給出血管內流體的詳細數學模型，且其模擬實驗使用文獻數據，可靠性沒有保障。近年來，一些研究雖然對于探究肝臟血管流體動力學有一定的指導意義，但尚無采用醫學-工學結合方法對灌注液的物質參數、運輸特性、灌注技術參數及其與血管相互間的影響在肝臟血管數學模型上進行理論計算與分析。
本文針對上述問題，首先嘗試使用三維實體建模軟件建立肝臟血管帶有異物(如纖維素沉積、膠原纖維變性、內皮細胞及肌細胞壞死等原因造成的異物堆積)的一級簡易模型(包括直管模型和彎曲血管模型)，其次對肝臟血管內的流體以及血管壁進行數學建模，給出微分方程。最后，通過實驗對專用灌注液的物理特性參數進行測量，結合手術中使用的相關參數，通過流-固耦合數據仿真，模擬不同初始條件下灌注過程中的運輸特性，以及在灌注過程中流體對血管壁和異物表面的沖擊作用等現象，并使用方程對沖擊力進行量化。
1 肝血管的實體建模
首先建立簡易小口徑直、彎曲血管模型。由于肝內靜脈、動脈的血管直徑不統一，為了便于計算和分析，本文血管模型內直徑為3 mm，血管外直徑為5 mm，血管壁厚1 mm，血管壁設為各向同性，現設定其物理參數[14-16]如下：彈性模量為1.2 MPa，泊松比為0.49，密度為1 120 kg/m3。
1.1 簡易直血管模型的建立
建立一級簡易直血管的三維實體模型，并模擬血管內異物(如纖維素沉積、膠原纖維變性、內皮細胞及肌細胞壞死等原因造成的異物堆積)。其模型如圖 1(a)所示。為計算方便，血管長設為20 mm，管內異物為半球形，為后續計算方便設定其半徑為1.5 mm。
圖1
帶有異物的直血管和彎曲血管模型示意圖
(a)正常直血管模型; (b)正常彎曲血管模型
Figure1. Models of straight and curved vessels with foreign matter(a) normal straight vessel model; (b) normal curved vessel model
1.2 簡易彎曲血管模型的建立
建立一級帶異物簡易彎曲血管模型如圖 1(b)所示。為了便于后續分析計算，假設血管內異物位于彎曲段的頂部。
2 建立管內流體動力學和管壁的數學模型
2.1 笛卡爾坐標系下血管內不可壓流體湍流瞬時控制方程
根據醫療現場使用數據計算，得出灌注液體在血管內是以湍流的形式存在的。一般情況下，無論湍流模型多么復雜，非穩態的連續方程和Navier-Stokes方程對于湍流的瞬時運動仍然是適用的[17]。不可壓縮流體湍流瞬時控制方程如下：
| $ div\boldsymbol{u}=0 $ |
| $ \frac{\partial u}{\partial t}+div\left(u\cdot \boldsymbol{u} \right)=-\frac{1}{\rho }\cdot \frac{\partial p}{\partial x}+\boldsymbol{v}\cdot div\left(grad\ u \right) $ |
| $ \frac{\partial v}{\partial t}+div\left(v\cdot \boldsymbol{u} \right)=-\frac{1}{\rho }\cdot \frac{\partial p}{\partial x}+\boldsymbol{v}\cdot div\left(grad\ v \right) $ |
| $ \frac{\partial w}{\partial t}+div\left(w\cdot \boldsymbol{u} \right)=-\frac{1}{\rho }\cdot \frac{\partial p}{\partial x}+\boldsymbol{v}\cdot div\left(grad\ w \right) $ |
其中p表示管內壓強, ρ表示管內液體密度, u、v和w分別為速度矢量u在x、y和z方向的分量, v為運動黏度。
為了考察管內脈動的影響，采用時間平均法，并引入Reynolds平均法[16], 用平均值與脈動值之和代替流動變量則有：
| $ \begin{matrix} \boldsymbol{u}=\overline{\boldsymbol{u}}+\boldsymbol{u}', u=\overline{u}+u', v=\overline{v}+v', w=\overline{w}+w' \ p=\overline{p}+p' \\ \end{matrix} $ |
將公式(5)代入公式(1)和公式(2)中，并對時間取平均，忽略密度脈動的影響，寫出可壓湍流平均流動控制方程[17]如下(為了方便起見，除脈動值外，下式中去掉了表示時均值的上劃線符號“-”)。
連續方程：
| $ \frac{\partial p}{\partial t}+div\left(\rho \ \ \boldsymbol{u} \right)=0 $ |
動量方程及運動方程：
| $ \frac{\partial \left(\rho u \right)}{\partial t}=div\left(\rho u\ \ \boldsymbol{u} \right)=-\frac{\partial p}{\partial x}+div\left(\mu \cdot grad\ u \right)+\left[-\frac{\partial \left(\rho \overline{u{{'}^{2}}} \right)}{\partial x}-\frac{\partial \left(\rho \overline{u'v'} \right)}{\partial y}-\frac{\partial \left(\rho \overline{u'w'} \right)}{\partial z} \right]+{{S}_{u}} $ |
| $ \frac{\partial \left(\rho v \right)}{\partial t}=div\left(\rho v\ \ \boldsymbol{u} \right)=-\frac{\partial p}{\partial y}+div\left(\mu \cdot grad\ v \right)+\left[-\frac{\partial \left(\rho \overline{u'v'} \right)}{\partial x}-\frac{\partial \left(\rho \overline{v{{'}^{2}}} \right)}{\partial y}-\frac{\partial \left(\rho \overline{v'w'} \right)}{\partial z} \right]+{{S}_{v}} $ |
| $ \frac{\partial \left(\rho w \right)}{\partial t}=div\left(\rho w\ \ \boldsymbol{u} \right)=-\frac{\partial p}{\partial z}+div\left(\mu \cdot grad\ w \right)+\left[-\frac{\partial \left(\rho \overline{u'w'} \right)}{\partial x}-\frac{\partial \left(\rho \overline{u'v'} \right)}{\partial y}-\frac{\partial \left(\rho \overline{w{{'}^{2}}} \right)}{\partial z} \right]+{{S}_{w}} $ |
其中μ為流體的動力黏度，Su，Sv，Sw是動量守恒方程的廣義源項。
Su=Fx+sx，Sv=Fy+sy，Sw=Fz+sz，對于黏性為常數的不可壓流體，sx=sy=sz=0，且Fx，Fy，Fz為微元上的體力，若體力只有重力，則Fx=Fy=0，Fz=-ρ·g。
2.2 邊界條件的設定
壁面條件：因為壁面上流體的速度為零，管壁面滿足無滑移條件，即：u=v=w=0，且在壁面處流體的速度與壁面平行，即:u⊥n。其中：n代表壁面的單位法向量。由于血管壁的徑向速度非常小，故將壁面的徑向速度簡化設定為零。
出口條件：根據流動充分發展條件，在血管末端出口處令其相對壓力為0 Pa。
設管壁異物處血管半徑為剩余通暢部分的y方向距離的函數。則可建立管壁的運動控制方程[18]：
| $ \nabla \cdot {{\sigma }_{\text{s}}}={{\rho }_{\text{s}}}{{a}_{\text{s}}} $ |
式中：σs是血管壁應力張量；ρs是血管壁密度；as是血管壁加速度。
2.3 流-固耦合交界面
交界面的處理對于流-固耦合是非常重要的，它是用來傳遞速度和位移。一般情況，交界面應滿足下列條件：
| $ {{d}_{\text{s}}}={{d}_{\text{f}}} $ |
| $ {{\sigma }_{\text{s}}}\cdot {{\boldsymbol{n}}_{\text{s}}}={{\sigma }_{\text{f}}}\cdot {{\boldsymbol{n}}_{\text{f}}} $ |
| $ {{u}_{\text{s}}}={{u}_{\text{f}}} $ |
其中：d為位移，n為邊界法向；下標s、f分別表示固體和流體。
在分析流-固耦合影響的管內液體流動問題時，對液體流動的計算和血管變形的計算應分別進行，并將各自的計算作為對方下次迭代計算的邊界條件或者初始條件。
3 灌注液的物理參數的測定
本研究采用羥乙基淀粉130/0.4氯化鈉注射液，其摩爾質量為150 000 g/mol[19]，摩爾取代度為0.38~0.45，本研究關注的灌注液參數為其密度ρ和動力黏度η。使用密度計來測試其密度如圖 2中左側所示。本灌注液為不可壓縮牛頓流體，多次采樣求取平均值得出其密度為1.029 g/cm3。對于其黏度使用回轉式液體黏度計進行測量，如圖 2中右側所示。在測定的過程中不斷提高轉子速度，即分別在6、12、30、60 r/min的轉速下測量3次，最后求取平均數，得出該種灌注液的動力黏度為5.762 5 mPa·s。后續的數據仿真將使用上述參數進行仿真計算。
圖2
測量灌注液密度及黏度所用實驗儀器
Figure2.
Experimental apparatus for density measurement and viscosity measurement of perfusion fluid
4 流-固耦合動力學仿真計算
4.1 網格劃分，初始條件及邊界條件設定
使用計算流體力學(computational fluid dynamics, CFD)軟件進行實體模型的的網格劃分；在Ansys12.0-CFX軟件中進行數據仿真，分別設定入口流速為0.3、0.5、0.8 m/s。離散格式采用混合格式，混合因子設定為0.75，最大殘差值設定為1×10-5，出口相對壓力設定為0 Pa，壁面無滑移，數值模擬時使用標準k-ε湍流模型進行描述。
4.2 仿真原理及結果
血管內灌注液沖擊異物時速度變化所引起的繞流現象如圖 3所示。當流速較小時，異物表面流場行為如圖 3左側所示，流線緊貼異物表面，由層流脫離表面變為湍流的分離點較靠后；當流速較大時，異物表面流場行為如圖 3右側所示，分離點前移，表面流線更早地由層流轉換為湍流。不可壓縮黏性流體繞流物體流動時，在物體表面由于應力而引起摩擦阻力Ff。并且，由于流線(體)在分離點與物體表面分離，變為湍流而產生能量損失，使得流線(體)繞曲壁面流動時，物體前后總壓強不同而產生壓差阻力Fp。它們的和為黏性流體繞流物體的阻力之和，稱為繞流阻力Fd[20]，為便于分析計算，參考文獻和有關實驗設其大小等于異物和血管壁的接觸摩擦力F阻。因而只要沖擊力(壓強)的大小超過繞流阻力或大于沖擊表面壓強p，即可認為異物被沖走(灌注沖洗成功)。
摩擦阻力Ff是黏性直接作用的結果，其計算過程在分離點之前的異物表面層流部分進行。壓差阻力Fp是黏性間接作用的結果，受分離層的影響很大，它與被繞流物體的形狀密切相關。以上兩種力表示為如下形式[20]：
| $ {{F}_{f}}={{C}_{f}}\frac{\rho {{V}^{2}}}{2}{{A}_{f}} $ |
| $ {{F}_{p}}={{C}_{p}}\frac{\rho {{V}^{2}}}{2}{{A}_{p}} $ |
圖3
灌注液遇血管內異物時動力學行為
Figure3.
Dynamic activities of the perfusion fluid with foreign matter in the vessel flowing around
則有
| $ {{F}_{d}}=\left({{C}_{f}}{{A}_{f}}+{{C}_{p}}{{A}_{p}} \right)\frac{\rho {{V}^{2}}}{2}={{C}_{d}}\frac{\rho {{V}^{2}}}{2}A $ |
式中：Cf和Cp分別為摩擦阻力系數和壓差阻力系數；Af為切應力作用的面積，Ap=A為物體向圖 3中X=0面即沖擊速度的法面的投影面積，若血管內異物為不規則物體，則需借助其他醫學手段測明其向沖擊速度的法面的投影面積。本模型中A=2.763 8×10-6 m2；Cd稱為繞流系數，主要取決于雷諾數，也與異物表面的粗糙度有關，本模型取值0.4。為了驗證分析，設異物重量m為0.05 g，血管內壁摩擦系數μ=0.75，被沖擊面的面積近似求得S=5.24×10-6 m2, 則有F阻=Fd=mgμ=0.000 05 kg×9.8 m/s2×0.75=0.000 367 5 N，得出沖擊表面的壓強近似為p=70 Pa(設標準大氣壓為0 Pa)。由公式(16)可得理論沖擊速度為V=0.803 806 6 m/s。為便于計算和分析，忽略由于流體的黏性造成的水頭損失，同時不考慮灌注液和血管壁的熱傳導的影響。現在使用該速度V作為入口初始速度進行直血管及彎曲血管模型的灌注仿真。
4.2.1 直血管內的仿真結果
本文分別采用V=0.3 m/s、V=0.5 m/s、V=0.8 m/s進行對比灌注仿真，異物表面的壓強分布結果如圖 4上排所示，數值結果如表 1所示：圖中灌注液從左側流向右側，為便于分析，此圖只顯示異物表面的壓強分布。三種速度下異物迎流面和背流面壓差達均大于p，說明用該速度沖擊異物成功。隨著灌注速度的增大，異物表面壓強也相應增大。
圖4
不同的灌注速度下直血管內異物表面壓強圖(上排)和管內異物周圍的流場速度云圖(下排)
Figure4.
Pressure distribution on foreign matter surface (upper figures) and flow velocity contours around foreign matter (lower figures) in straight vessel under the different perfusion velocities
表1
直血管不同灌注速度下異物表面壓力場和速度場
Table1.
Pressure field and velocity field on the foreign matter surface under different perfusion velocities in the straight vessel
圖 4下排顯示了在上述三種灌注速度下異物周圍的流場速度云圖，管內最大速度均出現在異物頂部，最小速度均出現在異物背流面底部。
4.2.2 彎曲血管內的仿真結果
圖 5上排依次為彎管血管內在三種灌注速度(0.3、0.5、0.8 m/s)下物體表面壓強分布圖，灌注液從左側流向右側。數值結果如表 2所示。異物的迎流面和背流面的壓差估值均大于前文所述p，即沖洗成功。
圖5
不同的灌注速度下彎曲血管內異物表面壓強圖(上排)和管內異物周圍的流場速度云圖(下排)
Figure5.
Pressure distribution on foreign matter surface (upper figures) and flow velocity contours around foreign matter (lower figures) in curved vessel under the different perfusion velocities
表2
彎曲血管不同灌注速度下異物表面壓力場和速度場
Table2.
Pressure field and velocity field on the foreign matter surface under different perfusion velocities in the blend vessel
圖 5下排為不同灌注速度下彎管內異物周圍的流場速度云圖。速度最大區域出現在異物斜上方，同時血管彎曲處速度場變化較大。
圖 6所示為在灌注速度0.8 m/s且為流固雙向耦合的情況下，直血管內壁的XY面正應力分布和血管內壁變形量。上圖為正應力分布絕對值：0~254.2 Pa。下圖為變形量分布：0~3.779×10-6 m，變形集中于異物處周圍。
圖6
流固雙向耦合且灌注速度0.8 m/s下的直血管壁的XY面正應力分布圖(上)以及血管壁變形圖(下)
Figure6.
Normal stress distribution (upper) on the XY plane of the straight vessel wall and the vessel wall deformation (lower) under the perfusion velocity 0.8 m/s with the condition of fluid-solid coupling
5 討論
5.1 血管內異物的表面壓強分布
從圖 4上排和圖 5上排對比分析得出：不論在直血管或彎曲血管中，隨著灌注速度的增大，管內異物表面的壓強變化范圍隨之增大。且在三種速度下表面壓強均呈層狀變化分布，壓強高的位置正是灌注液速度較低處，與流體力學行為相符。異物迎流面和背流面的壓差也與灌注速度呈正相關，這在血管承受范圍內對于沖擊異物是有利的。
由于血管內異物頂部速度最大，兩種類型的血管在該處均出現了負的壁透壓。這是由于黏性流體在繞流物體時，頂部速度較高，壓強較小，使異物內部壓強大于外部壓強，因此產生一個向外的升力，且該值大小與入口初始灌注速度正相關，其計算形式類似于繞流阻力。同樣，壁透壓在血管承受范圍內有利于異物沖洗，但直管內該壁透壓值大于彎曲血管，這是由于彎曲血管內水頭損失大于直管而造成的能量損失。
5.2 血管內流場分布
兩種類型的血管隨著灌注速度的增加，血管內的流速場均增大。由于管內液體有黏性，需要克服內部的黏性阻力來維持流動，從而造成能量損失。在直管內主要發生沿程損失，在彎管內還會發生由于血管形狀變化造成的局部水頭損失，因此會導致彎曲血管內流速場整體小于直血管。同時在彎曲血管的異物周圍和兩個血管拐點處存在較大的流速場變化。三種灌注速度下均在異物頂部出現了最大流速場，直血管內該場位于異物正上方，彎曲血管內該場位置稍稍靠后，導致這里異物的表面壁透壓(向外升力)也是斜向的，比正向外的升力更加有利于沖洗，會減小對血管的損傷。
兩種類型的血管在不同灌注速度下，異物的背流面底部均出現了回流滯留區。該區域的出現使得原本的異物后部的沖擊速度降低，在這一部分流速場大小介于異物頂部和迎流面之間。該回流滯留區的出現不利于異物清洗，而且有可能導致壁面剪切力較低，造成流動淤滯區[21]，形成新的阻塞。由圖中分析可以得出，沖擊速度越大，這一部分的回流現象越明顯，因此可以在不影響沖洗成功的前提下適當降低灌注速度來避免這種情況的出現。同樣，在血管彎曲處也會產生回流滯留區。綜上所述，在血管內的凸處的背流面、凹處的中部均會出現低壁面剪切力的回流區，這些地方有可能形成新的阻塞等。
同時，在三種灌注速度下，如圖 4下排和圖 5下排均存在流體緊貼異物表面的部分流速較小的現象。這是由于黏性流體在繞流時貼近物體表面的部分會形成很窄的一層層流，而這層流狀態會隨著流動在分離點消失，在物體后半部成為湍流，而當加大灌注速度，分離現象會更加明顯，分離點會“前移”，如圖 3(b)所示。
5.3 血管壁面壓強分布以及對于術前灌洗過程的理論指導
根據天津市第一中心醫院術前清洗醫療現場的數據，灌注液設置在距灌注入口高1 m處通過管道注入肝供體進行清洗，即以V=4.5 m/s的速度進行灌注。該參數值是依托經驗和動物實驗數據，還沒有理論分析數據的支撐。
本研究針對以上問題提出了一種根據血管內異物的大小和形狀來初步計算灌注速度的方法，同時考慮壁面壓強。由造影分析血管中異物的大小，并由公式(16)算出灌注速度值，用此速度能夠成功沖洗異物。本研究參照手術參數選取V=0.5 m/s和V=0.3 m/s進行對比試驗。從圖 6可以得出，在估算得出的灌注速度下，壓強值均小于人體血管正常壓強值(收縮壓90~120 mm Hg，舒張壓60~90 mm Hg)，故不會對血管造成破壞。根據有關文獻的報道[22]，人體正常的壁面剪切力極值在40 Pa左右，過大的剪切力會導致血管壁受損。最合理的方法是在保證沖洗能夠成功的前提下，使血管內壁面承受的壓強和壁面剪切力盡可能不大于其極限值。
6 結論
本研究針對肝臟血管內流體動力學行為進行了定量的數據仿真和機制分析。結合流-固耦合的方法對灌注液在管內清洗中的動力學行為進行了較為詳細的闡述。由于灌注液在血管內沖洗異物時的流動為湍流，其流動動力學行為比較復雜，在異物背流面易形成回流滯留區而導致壁面剪切力較低，進而有可能產生新的阻塞。而且流體經過異物時都會產生一定的繞流現象。受血管形狀的影響，彎曲血管內的流體動力學行為比直血管中更加復雜，壁透壓的絕對值均小于直血管，其管內流速場變化更大。本研究提供了一種可行的估算灌注速度場大小的方法，可以為改進供肝灌注流程與質量提供有價值的理論指導。
引言
目前國內外不論是發達地區還是欠發達地區肝病誘發的死亡率都很高，肝移植是而目前治療終末期肝病的唯一有效手段。肝供體質量高低是影響肝移植手術成功率的關鍵環節之一，術前對肝供體的清洗(灌注)是必須操作的醫療過程。對肝臟內部血管的流體動力學特性，尤其是灌注過程的動力學行為的研究，有助于提高肝臟供體的灌注質量和對肝供體的保護。對血管建模并進行相應的力學分析，將有助于降低肝供體灌注過程中的主觀因素對清洗質量的影響，同時為規范操作提供理論依據。
劉向明等[1]對生物醫學中的數學建模問題進行了綜述，闡明了其分類和各自特點以及研究目標。王振全等[2]概述了目前國內肝臟血流動力學臨床應用的研究進展。肝臟血管最多可向下分至8級，而目前，國內外對于肝臟血管機制的研究最多能達到前4級。為了構建肝臟血管的結構模型，在血管組織及其結構的可視化研究基礎上，國內外大多利用超聲成像[3-5]、核磁共振[6-7]和計算機斷層成像(computed tomography, CT)[8-10]等技術來獲得血管系統的二維或三維實際圖像數據，然后將數據導入三維造型軟件中進行肝臟血管模型的建立。也有研究者直接利用三維造型技術進行血管系統的模型構建，將重點放在血管和血液的流-固耦合及其相互影響和作用的研究上。近年來，研究者也提出了一些新的原理和方法，如Hoganson等[11]應用仿生原理進行肝臟血管模型的優化及設計，開發了肝臟血管的實體模型。Matsumoto等[12]利用二維灰階血流成像技術(two-dimensional blood flow, B-Flow)對肝臟血管進行圖像分析，并將該技術與彩色多普勒超聲技術進行比較，探討了二者的異同點和優缺點。
對肝移植手術中的醫學問題進行定量分析涉及有關工學領域的知識，如肝臟血管數學模型、傳感器檢測和血流動力學計算等方面。國內外對肝移植的研究主要是建立在核磁共振或CT技術提供的肝臟血流分布信息上，并依此來進行相關研究。Debbaut等[13]利用CT技術對肝臟血管進行圖像建模，使用文獻數據模擬肝臟內低溫機械灌注情況下的液體流動狀況，得出灌注壓力及流體黏度對于血管壁面剪切應力(wall shear stress, WSS)的影響。但并未給出血管內流體的詳細數學模型，且其模擬實驗使用文獻數據，可靠性沒有保障。近年來，一些研究雖然對于探究肝臟血管流體動力學有一定的指導意義，但尚無采用醫學-工學結合方法對灌注液的物質參數、運輸特性、灌注技術參數及其與血管相互間的影響在肝臟血管數學模型上進行理論計算與分析。
本文針對上述問題，首先嘗試使用三維實體建模軟件建立肝臟血管帶有異物(如纖維素沉積、膠原纖維變性、內皮細胞及肌細胞壞死等原因造成的異物堆積)的一級簡易模型(包括直管模型和彎曲血管模型)，其次對肝臟血管內的流體以及血管壁進行數學建模，給出微分方程。最后，通過實驗對專用灌注液的物理特性參數進行測量，結合手術中使用的相關參數，通過流-固耦合數據仿真，模擬不同初始條件下灌注過程中的運輸特性，以及在灌注過程中流體對血管壁和異物表面的沖擊作用等現象，并使用方程對沖擊力進行量化。
1 肝血管的實體建模
首先建立簡易小口徑直、彎曲血管模型。由于肝內靜脈、動脈的血管直徑不統一，為了便于計算和分析，本文血管模型內直徑為3 mm，血管外直徑為5 mm，血管壁厚1 mm，血管壁設為各向同性，現設定其物理參數[14-16]如下：彈性模量為1.2 MPa，泊松比為0.49，密度為1 120 kg/m3。
1.1 簡易直血管模型的建立
建立一級簡易直血管的三維實體模型，并模擬血管內異物(如纖維素沉積、膠原纖維變性、內皮細胞及肌細胞壞死等原因造成的異物堆積)。其模型如圖 1(a)所示。為計算方便，血管長設為20 mm，管內異物為半球形，為后續計算方便設定其半徑為1.5 mm。
圖1
帶有異物的直血管和彎曲血管模型示意圖
(a)正常直血管模型; (b)正常彎曲血管模型
Figure1. Models of straight and curved vessels with foreign matter(a) normal straight vessel model; (b) normal curved vessel model
1.2 簡易彎曲血管模型的建立
建立一級帶異物簡易彎曲血管模型如圖 1(b)所示。為了便于后續分析計算，假設血管內異物位于彎曲段的頂部。
2 建立管內流體動力學和管壁的數學模型
2.1 笛卡爾坐標系下血管內不可壓流體湍流瞬時控制方程
根據醫療現場使用數據計算，得出灌注液體在血管內是以湍流的形式存在的。一般情況下，無論湍流模型多么復雜，非穩態的連續方程和Navier-Stokes方程對于湍流的瞬時運動仍然是適用的[17]。不可壓縮流體湍流瞬時控制方程如下：
| $ div\boldsymbol{u}=0 $ |
| $ \frac{\partial u}{\partial t}+div\left(u\cdot \boldsymbol{u} \right)=-\frac{1}{\rho }\cdot \frac{\partial p}{\partial x}+\boldsymbol{v}\cdot div\left(grad\ u \right) $ |
| $ \frac{\partial v}{\partial t}+div\left(v\cdot \boldsymbol{u} \right)=-\frac{1}{\rho }\cdot \frac{\partial p}{\partial x}+\boldsymbol{v}\cdot div\left(grad\ v \right) $ |
| $ \frac{\partial w}{\partial t}+div\left(w\cdot \boldsymbol{u} \right)=-\frac{1}{\rho }\cdot \frac{\partial p}{\partial x}+\boldsymbol{v}\cdot div\left(grad\ w \right) $ |
其中p表示管內壓強, ρ表示管內液體密度, u、v和w分別為速度矢量u在x、y和z方向的分量, v為運動黏度。
為了考察管內脈動的影響，采用時間平均法，并引入Reynolds平均法[16], 用平均值與脈動值之和代替流動變量則有：
| $ \begin{matrix} \boldsymbol{u}=\overline{\boldsymbol{u}}+\boldsymbol{u}', u=\overline{u}+u', v=\overline{v}+v', w=\overline{w}+w' \ p=\overline{p}+p' \\ \end{matrix} $ |
將公式(5)代入公式(1)和公式(2)中，并對時間取平均，忽略密度脈動的影響，寫出可壓湍流平均流動控制方程[17]如下(為了方便起見，除脈動值外，下式中去掉了表示時均值的上劃線符號“-”)。
連續方程：
| $ \frac{\partial p}{\partial t}+div\left(\rho \ \ \boldsymbol{u} \right)=0 $ |
動量方程及運動方程：
| $ \frac{\partial \left(\rho u \right)}{\partial t}=div\left(\rho u\ \ \boldsymbol{u} \right)=-\frac{\partial p}{\partial x}+div\left(\mu \cdot grad\ u \right)+\left[-\frac{\partial \left(\rho \overline{u{{'}^{2}}} \right)}{\partial x}-\frac{\partial \left(\rho \overline{u'v'} \right)}{\partial y}-\frac{\partial \left(\rho \overline{u'w'} \right)}{\partial z} \right]+{{S}_{u}} $ |
| $ \frac{\partial \left(\rho v \right)}{\partial t}=div\left(\rho v\ \ \boldsymbol{u} \right)=-\frac{\partial p}{\partial y}+div\left(\mu \cdot grad\ v \right)+\left[-\frac{\partial \left(\rho \overline{u'v'} \right)}{\partial x}-\frac{\partial \left(\rho \overline{v{{'}^{2}}} \right)}{\partial y}-\frac{\partial \left(\rho \overline{v'w'} \right)}{\partial z} \right]+{{S}_{v}} $ |
| $ \frac{\partial \left(\rho w \right)}{\partial t}=div\left(\rho w\ \ \boldsymbol{u} \right)=-\frac{\partial p}{\partial z}+div\left(\mu \cdot grad\ w \right)+\left[-\frac{\partial \left(\rho \overline{u'w'} \right)}{\partial x}-\frac{\partial \left(\rho \overline{u'v'} \right)}{\partial y}-\frac{\partial \left(\rho \overline{w{{'}^{2}}} \right)}{\partial z} \right]+{{S}_{w}} $ |
其中μ為流體的動力黏度，Su，Sv，Sw是動量守恒方程的廣義源項。
Su=Fx+sx，Sv=Fy+sy，Sw=Fz+sz，對于黏性為常數的不可壓流體，sx=sy=sz=0，且Fx，Fy，Fz為微元上的體力，若體力只有重力，則Fx=Fy=0，Fz=-ρ·g。
2.2 邊界條件的設定
壁面條件：因為壁面上流體的速度為零，管壁面滿足無滑移條件，即：u=v=w=0，且在壁面處流體的速度與壁面平行，即:u⊥n。其中：n代表壁面的單位法向量。由于血管壁的徑向速度非常小，故將壁面的徑向速度簡化設定為零。
出口條件：根據流動充分發展條件，在血管末端出口處令其相對壓力為0 Pa。
設管壁異物處血管半徑為剩余通暢部分的y方向距離的函數。則可建立管壁的運動控制方程[18]：
| $ \nabla \cdot {{\sigma }_{\text{s}}}={{\rho }_{\text{s}}}{{a}_{\text{s}}} $ |
式中：σs是血管壁應力張量；ρs是血管壁密度；as是血管壁加速度。
2.3 流-固耦合交界面
交界面的處理對于流-固耦合是非常重要的，它是用來傳遞速度和位移。一般情況，交界面應滿足下列條件：
| $ {{d}_{\text{s}}}={{d}_{\text{f}}} $ |
| $ {{\sigma }_{\text{s}}}\cdot {{\boldsymbol{n}}_{\text{s}}}={{\sigma }_{\text{f}}}\cdot {{\boldsymbol{n}}_{\text{f}}} $ |
| $ {{u}_{\text{s}}}={{u}_{\text{f}}} $ |
其中：d為位移，n為邊界法向；下標s、f分別表示固體和流體。
在分析流-固耦合影響的管內液體流動問題時，對液體流動的計算和血管變形的計算應分別進行，并將各自的計算作為對方下次迭代計算的邊界條件或者初始條件。
3 灌注液的物理參數的測定
本研究采用羥乙基淀粉130/0.4氯化鈉注射液，其摩爾質量為150 000 g/mol[19]，摩爾取代度為0.38~0.45，本研究關注的灌注液參數為其密度ρ和動力黏度η。使用密度計來測試其密度如圖 2中左側所示。本灌注液為不可壓縮牛頓流體，多次采樣求取平均值得出其密度為1.029 g/cm3。對于其黏度使用回轉式液體黏度計進行測量，如圖 2中右側所示。在測定的過程中不斷提高轉子速度，即分別在6、12、30、60 r/min的轉速下測量3次，最后求取平均數，得出該種灌注液的動力黏度為5.762 5 mPa·s。后續的數據仿真將使用上述參數進行仿真計算。
圖2
測量灌注液密度及黏度所用實驗儀器
Figure2.
Experimental apparatus for density measurement and viscosity measurement of perfusion fluid
4 流-固耦合動力學仿真計算
4.1 網格劃分，初始條件及邊界條件設定
使用計算流體力學(computational fluid dynamics, CFD)軟件進行實體模型的的網格劃分；在Ansys12.0-CFX軟件中進行數據仿真，分別設定入口流速為0.3、0.5、0.8 m/s。離散格式采用混合格式，混合因子設定為0.75，最大殘差值設定為1×10-5，出口相對壓力設定為0 Pa，壁面無滑移，數值模擬時使用標準k-ε湍流模型進行描述。
4.2 仿真原理及結果
血管內灌注液沖擊異物時速度變化所引起的繞流現象如圖 3所示。當流速較小時，異物表面流場行為如圖 3左側所示，流線緊貼異物表面，由層流脫離表面變為湍流的分離點較靠后；當流速較大時，異物表面流場行為如圖 3右側所示，分離點前移，表面流線更早地由層流轉換為湍流。不可壓縮黏性流體繞流物體流動時，在物體表面由于應力而引起摩擦阻力Ff。并且，由于流線(體)在分離點與物體表面分離，變為湍流而產生能量損失，使得流線(體)繞曲壁面流動時，物體前后總壓強不同而產生壓差阻力Fp。它們的和為黏性流體繞流物體的阻力之和，稱為繞流阻力Fd[20]，為便于分析計算，參考文獻和有關實驗設其大小等于異物和血管壁的接觸摩擦力F阻。因而只要沖擊力(壓強)的大小超過繞流阻力或大于沖擊表面壓強p，即可認為異物被沖走(灌注沖洗成功)。
摩擦阻力Ff是黏性直接作用的結果，其計算過程在分離點之前的異物表面層流部分進行。壓差阻力Fp是黏性間接作用的結果，受分離層的影響很大，它與被繞流物體的形狀密切相關。以上兩種力表示為如下形式[20]：
| $ {{F}_{f}}={{C}_{f}}\frac{\rho {{V}^{2}}}{2}{{A}_{f}} $ |
| $ {{F}_{p}}={{C}_{p}}\frac{\rho {{V}^{2}}}{2}{{A}_{p}} $ |
圖3
灌注液遇血管內異物時動力學行為
Figure3.
Dynamic activities of the perfusion fluid with foreign matter in the vessel flowing around
則有
| $ {{F}_{d}}=\left({{C}_{f}}{{A}_{f}}+{{C}_{p}}{{A}_{p}} \right)\frac{\rho {{V}^{2}}}{2}={{C}_{d}}\frac{\rho {{V}^{2}}}{2}A $ |
式中：Cf和Cp分別為摩擦阻力系數和壓差阻力系數；Af為切應力作用的面積，Ap=A為物體向圖 3中X=0面即沖擊速度的法面的投影面積，若血管內異物為不規則物體，則需借助其他醫學手段測明其向沖擊速度的法面的投影面積。本模型中A=2.763 8×10-6 m2；Cd稱為繞流系數，主要取決于雷諾數，也與異物表面的粗糙度有關，本模型取值0.4。為了驗證分析，設異物重量m為0.05 g，血管內壁摩擦系數μ=0.75，被沖擊面的面積近似求得S=5.24×10-6 m2, 則有F阻=Fd=mgμ=0.000 05 kg×9.8 m/s2×0.75=0.000 367 5 N，得出沖擊表面的壓強近似為p=70 Pa(設標準大氣壓為0 Pa)。由公式(16)可得理論沖擊速度為V=0.803 806 6 m/s。為便于計算和分析，忽略由于流體的黏性造成的水頭損失，同時不考慮灌注液和血管壁的熱傳導的影響。現在使用該速度V作為入口初始速度進行直血管及彎曲血管模型的灌注仿真。
4.2.1 直血管內的仿真結果
本文分別采用V=0.3 m/s、V=0.5 m/s、V=0.8 m/s進行對比灌注仿真，異物表面的壓強分布結果如圖 4上排所示，數值結果如表 1所示：圖中灌注液從左側流向右側，為便于分析，此圖只顯示異物表面的壓強分布。三種速度下異物迎流面和背流面壓差達均大于p，說明用該速度沖擊異物成功。隨著灌注速度的增大，異物表面壓強也相應增大。
圖4
不同的灌注速度下直血管內異物表面壓強圖(上排)和管內異物周圍的流場速度云圖(下排)
Figure4.
Pressure distribution on foreign matter surface (upper figures) and flow velocity contours around foreign matter (lower figures) in straight vessel under the different perfusion velocities
表1
直血管不同灌注速度下異物表面壓力場和速度場
Table1.
Pressure field and velocity field on the foreign matter surface under different perfusion velocities in the straight vessel
圖 4下排顯示了在上述三種灌注速度下異物周圍的流場速度云圖，管內最大速度均出現在異物頂部，最小速度均出現在異物背流面底部。
4.2.2 彎曲血管內的仿真結果
圖 5上排依次為彎管血管內在三種灌注速度(0.3、0.5、0.8 m/s)下物體表面壓強分布圖，灌注液從左側流向右側。數值結果如表 2所示。異物的迎流面和背流面的壓差估值均大于前文所述p，即沖洗成功。
圖5
不同的灌注速度下彎曲血管內異物表面壓強圖(上排)和管內異物周圍的流場速度云圖(下排)
Figure5.
Pressure distribution on foreign matter surface (upper figures) and flow velocity contours around foreign matter (lower figures) in curved vessel under the different perfusion velocities
表2
彎曲血管不同灌注速度下異物表面壓力場和速度場
Table2.
Pressure field and velocity field on the foreign matter surface under different perfusion velocities in the blend vessel
圖 5下排為不同灌注速度下彎管內異物周圍的流場速度云圖。速度最大區域出現在異物斜上方，同時血管彎曲處速度場變化較大。
圖 6所示為在灌注速度0.8 m/s且為流固雙向耦合的情況下，直血管內壁的XY面正應力分布和血管內壁變形量。上圖為正應力分布絕對值：0~254.2 Pa。下圖為變形量分布：0~3.779×10-6 m，變形集中于異物處周圍。
圖6
流固雙向耦合且灌注速度0.8 m/s下的直血管壁的XY面正應力分布圖(上)以及血管壁變形圖(下)
Figure6.
Normal stress distribution (upper) on the XY plane of the straight vessel wall and the vessel wall deformation (lower) under the perfusion velocity 0.8 m/s with the condition of fluid-solid coupling
5 討論
5.1 血管內異物的表面壓強分布
從圖 4上排和圖 5上排對比分析得出：不論在直血管或彎曲血管中，隨著灌注速度的增大，管內異物表面的壓強變化范圍隨之增大。且在三種速度下表面壓強均呈層狀變化分布，壓強高的位置正是灌注液速度較低處，與流體力學行為相符。異物迎流面和背流面的壓差也與灌注速度呈正相關，這在血管承受范圍內對于沖擊異物是有利的。
由于血管內異物頂部速度最大，兩種類型的血管在該處均出現了負的壁透壓。這是由于黏性流體在繞流物體時，頂部速度較高，壓強較小，使異物內部壓強大于外部壓強，因此產生一個向外的升力，且該值大小與入口初始灌注速度正相關，其計算形式類似于繞流阻力。同樣，壁透壓在血管承受范圍內有利于異物沖洗，但直管內該壁透壓值大于彎曲血管，這是由于彎曲血管內水頭損失大于直管而造成的能量損失。
5.2 血管內流場分布
兩種類型的血管隨著灌注速度的增加，血管內的流速場均增大。由于管內液體有黏性，需要克服內部的黏性阻力來維持流動，從而造成能量損失。在直管內主要發生沿程損失，在彎管內還會發生由于血管形狀變化造成的局部水頭損失，因此會導致彎曲血管內流速場整體小于直血管。同時在彎曲血管的異物周圍和兩個血管拐點處存在較大的流速場變化。三種灌注速度下均在異物頂部出現了最大流速場，直血管內該場位于異物正上方，彎曲血管內該場位置稍稍靠后，導致這里異物的表面壁透壓(向外升力)也是斜向的，比正向外的升力更加有利于沖洗，會減小對血管的損傷。
兩種類型的血管在不同灌注速度下，異物的背流面底部均出現了回流滯留區。該區域的出現使得原本的異物后部的沖擊速度降低，在這一部分流速場大小介于異物頂部和迎流面之間。該回流滯留區的出現不利于異物清洗，而且有可能導致壁面剪切力較低，造成流動淤滯區[21]，形成新的阻塞。由圖中分析可以得出，沖擊速度越大，這一部分的回流現象越明顯，因此可以在不影響沖洗成功的前提下適當降低灌注速度來避免這種情況的出現。同樣，在血管彎曲處也會產生回流滯留區。綜上所述，在血管內的凸處的背流面、凹處的中部均會出現低壁面剪切力的回流區，這些地方有可能形成新的阻塞等。
同時，在三種灌注速度下，如圖 4下排和圖 5下排均存在流體緊貼異物表面的部分流速較小的現象。這是由于黏性流體在繞流時貼近物體表面的部分會形成很窄的一層層流，而這層流狀態會隨著流動在分離點消失，在物體后半部成為湍流，而當加大灌注速度，分離現象會更加明顯，分離點會“前移”，如圖 3(b)所示。
5.3 血管壁面壓強分布以及對于術前灌洗過程的理論指導
根據天津市第一中心醫院術前清洗醫療現場的數據，灌注液設置在距灌注入口高1 m處通過管道注入肝供體進行清洗，即以V=4.5 m/s的速度進行灌注。該參數值是依托經驗和動物實驗數據，還沒有理論分析數據的支撐。
本研究針對以上問題提出了一種根據血管內異物的大小和形狀來初步計算灌注速度的方法，同時考慮壁面壓強。由造影分析血管中異物的大小，并由公式(16)算出灌注速度值，用此速度能夠成功沖洗異物。本研究參照手術參數選取V=0.5 m/s和V=0.3 m/s進行對比試驗。從圖 6可以得出，在估算得出的灌注速度下，壓強值均小于人體血管正常壓強值(收縮壓90~120 mm Hg，舒張壓60~90 mm Hg)，故不會對血管造成破壞。根據有關文獻的報道[22]，人體正常的壁面剪切力極值在40 Pa左右，過大的剪切力會導致血管壁受損。最合理的方法是在保證沖洗能夠成功的前提下，使血管內壁面承受的壓強和壁面剪切力盡可能不大于其極限值。
6 結論
本研究針對肝臟血管內流體動力學行為進行了定量的數據仿真和機制分析。結合流-固耦合的方法對灌注液在管內清洗中的動力學行為進行了較為詳細的闡述。由于灌注液在血管內沖洗異物時的流動為湍流，其流動動力學行為比較復雜，在異物背流面易形成回流滯留區而導致壁面剪切力較低，進而有可能產生新的阻塞。而且流體經過異物時都會產生一定的繞流現象。受血管形狀的影響，彎曲血管內的流體動力學行為比直血管中更加復雜，壁透壓的絕對值均小于直血管，其管內流速場變化更大。本研究提供了一種可行的估算灌注速度場大小的方法，可以為改進供肝灌注流程與質量提供有價值的理論指導。
