基于兩相流歐拉方法的中空纖維管血流動力學數值模擬

朱 橋 輝, 許 少 鋒, 王 子 恒, 陸 俊 杰, 張 學 昌

(1.浙江大學 機械工程學院, 浙江 杭州 310027;2.浙大寧波理工學院 機電與能源工程學院, 浙江 寧波 315100;3.浙江理工大學 機械與自動控制學院, 浙江 杭州 310018 )

0 引 言

人工肝支持系統(artificial liver support system,ALSS)簡稱人工肝,可將人體血液引出至體外機械裝置幫助凈化人體血液內各種有害物質、改善血液內環境,暫時替代和輔助肝臟代謝、解毒等功能[1]．人工肝在臨床上可用于肝病重癥、肝衰竭患者的相關治療[2]．在對危重癥肺炎患者的治療中[3],人工肝可用于清除炎癥介質等有害物質,補充白蛋白等有益物質,從而阻斷細胞因子風暴,幫助患者改善呼吸功能,減輕肺炎癥狀,提高危重癥肺炎患者的救治率．中空纖維組件是人工肝的核心裝置,其中有大量的中空纖維管,當血液流經中空纖維管時,血液內的有害物質通過中空纖維膜進入外腔,從而達到凈化血液的目的,研究血液在中空纖維管內的流動行為對研究人工肝具有科學意義．

Baskurt等[4]的研究表明血液的表觀黏度取決于現有的剪切力,并由紅細胞比體積、血漿黏度、紅細胞聚集和紅細胞的力學特性決定．Shariatkhah等[5]使用非線性黏彈性模型研究了毛細血管中周期性發展的血流,并根據實驗數據與Giesekus 模型推導出了流動因子和零剪切速率黏度的最佳值,結果表明只有非線性黏彈性模型才能準確描述毛細血管中血流的實驗數據．Apostolidis等[6]研究了動脈血流動中的非牛頓血液流變學效應,結果顯示采用牛頓模型與非牛頓模型描述血液黏度時模擬結果存在顯著差異,這些差異主要歸因于流動中低剪切速率區域和高剪切速率區域之間存在的耦合．Hund等[7]在Krieger懸浮液模型的基礎上,提出了新的血液黏度模型,解決了先前模型在紅細胞比體積和剪切速率方面表現出的不連續性問題．Wu等[8]使用混合理論模擬了微通道的血流,紅細胞被建模為具有剪切依賴性黏度的非線性流體,并考慮了紅細胞比體積的影響,結果顯示在速度場和紅細胞的體積分數分布方面與實驗觀察結果非常吻合．Marhefka等[9]對血液在微尺度流動進行了實驗觀察,結果顯示血液在微管或小管(直徑小于0.3 mm)內流動時在靠近管壁的地方會出現一層薄薄的血漿層,稱為紅細胞的耗竭層,在耗竭層下游的分支血管中會出現血漿撇渣效應,導致在較小的血管分支中紅細胞的減少,可能會損害血液在遠端毛細血管中氧氣的運輸．其他學者還發現了F?hr?us-Lindqvist效應、塞子流等現象[10]．

血液的微觀流動已經得到了大量研究,但對血液中紅細胞受力的作用機理和遷移規律的研究還遠遠不夠．為了真實反映出人工肝中空纖維管內的血流動力學和紅細胞遷移規律,本文對中空纖維管內血液兩相流動進行三維數值模擬,構建中空纖維管的三維模型,采用歐拉雙流體模型,并使用用戶自定義函數(UDF)描述紅細胞升力方程,模擬人工肝中空纖維管內血液的流動情況,討論壁面升力作用對紅細胞徑向移動的影響,分析中空纖維管內血液兩相流體的速度分布、紅細胞徑向體積分數分布、黏度分布以及壁面條件對紅細胞體積分數分布的影響,并根據相關數值結果分析紅細胞所受升力的作用機理．研究結果旨在探討人工肝內血液微觀流動行為規律,以期為中空纖維組件的選材與優化設計提供一定的理論指導,提高研發效率,節約研發成本．

1 模型建立

1.1 數值模型

人體血液是由紅細胞、白細胞、血小板、血漿等組成的混合流體,血液中紅細胞與血漿的比例大于0.99,而白細胞與血小板所占比例不到0.01,一般忽略不計．因此血液被看作紅細胞與血漿組成的兩相系統,其中血漿在血液內的體積分數εp=0.55,紅細胞在血液內的體積分數εRBC=0.45[11]．對血液流動模擬計算時,血漿被看作液相,紅細胞被看作顆粒相(固相),由于紅細胞的體積分數較高,可將顆粒相紅細胞假設為擬流體,模擬采用的兩相流模型為歐拉雙流體液-固兩相流模型．

1.1.1 血液兩相流基本方程 采用歐拉雙流體模型,各相均可用統一的質量和動量守恒方程描述．質量守恒方程[12]為

(1)

式中：k為p或RBC,p、RBC分別代表血漿與紅細胞;ρ為密度;ε為體積分數;v為速度矢量．血漿和紅細胞體積分數之和恒為1,可表示為

(2)

動量守恒方程[12]為

(3)

式中：p為流場壓力,τk為應力張量,Fk為虛擬質量力、浮升力等力源項,β為相間動量交換系數,g為重力．

1.1.2 本構方程 為使上述基本方程組封閉,需給出密度、黏度等參數對應的本構方程．混合密度ρmix是由血漿和紅細胞的密度按其所占的體積分數加權計算得到,表達式為

(4)

1.1.3 血液黏度模型 血液黏度模型是模擬血液動力學流動的關鍵因素．非牛頓剪切稀化模型能夠描述中性稠密浮力懸浮液中高體積分數顆粒的流動．采用Carreau-Yasuda模型來模擬血液流動的非牛頓特性,血液黏度會隨著紅細胞體積分數和剪切速率的變化而變化,根據人體血液的實驗流變學數據,計算得到混合物的具體量綱一黏度表達式[13-14]為

(5)

(6)

(7)

正常生理條件下,血漿是牛頓流體,血漿黏度[16]為0.001 2 kg/(m·s),而紅細胞黏度會隨其體積分數的變化而變化,呈現出剪切稀化特性[17]．

液相(血漿)應力張量以牛頓形式表示,本構方程為

(8)

式中：κp表示血漿的體積黏度,I表示單位應力張量．一般條件下,血漿的體積黏度κp設為0．

固相(紅細胞)應力張量的本構方程[18]為

τRBC=-pRBCδ+εRBCμRBC(?v+(?v)T)+

(9)

式中：δ表示Kronecker函數;pRBC表示紅細胞顆粒間因發生相互碰撞、排斥作用而產生的壓力,由于紅細胞的比體積范圍內顆粒間產生的壓力微乎其微,可忽略不計;κRBC表示紅細胞的體積黏度,其數值為0;紅細胞黏度μRBC采用Carreau-Yasuda模型計算得到．

1.1.4 曳力系數 顆粒相與液相(血漿)之間的曳力采用Schiller-Naumann模型來表示相間動量交換系數β,其表達式[12]為

(10)

其中曳力系數Cd表達式為

(11)

(12)

式中：dRBC表示紅細胞顆粒的直徑,φ表示形狀因子．紅細胞顆粒的直徑被假定為8 μm,由于不考慮紅細胞變形影響,其形狀因子φ為1[15]．

1.1.5 虛擬質量力 當加速的顆粒相與液相(血漿)發生碰撞時,由于具有一定慣性作用,顆粒相會受到虛擬質量力,其表達式[19]為

(13)

1.1.6 升力方程 受限剪切流場中,在壁面條件和徑向速度梯度的作用下,考慮到紅細胞變形作用[20],紅細胞會受到經過其軸心的非慣性升力作用[21-22],該作用力會導致紅細胞發生流場內的徑向運動,對血液的流場分布、紅細胞的徑向分布、壁面剪切力產生影響．

兩相流模型中紅細胞所受的非慣性升力作用,當紅細胞的軸心距離壁面小于5 μm時,升力作用由壁面條件決定;當紅細胞的軸心距離壁面大于5 μm時,升力作用由剪切流場中的速度梯度引起,其表達式[22-23]為

(14)

式中：U表示量綱一滑移速度,與紅細胞內部、外部流體黏度比有關,在本式中U恒為常數,即U=0.03[19,24];h表示紅細胞軸心與壁面的距離;f(s)是以s為自變量的量綱一函數[25-26],s是指與紅細胞具有相同表面積的球體相比紅細胞的減少量,假設紅細胞在h=5 μm處不再和壁面發生接觸,由壁面條件所引起的升力[25-26]為0,并假設f(s)隨著h呈線性變化,當h=0 μm時,f(s)=2;當h=5 μm時,f(s)=0．

1.2 幾何模型與求解設置

由于1.1節中用于血液流動模擬的兩相流數學模型比較復雜,在模擬中空纖維管內血液流動時,為了降低數值求解過程對計算機性能的要求以及相關計算參數的求解,中空纖維管模型采用三維軸對稱結構,模型示意圖如圖1所示,中空纖維管截面直徑Ф為100 μm,長度L為5 mm．

圖1 中空纖維管模型示意圖

采用ANSYS ICEM 2020R2對三維結構離散化,模型網格劃分采用六面體網格法,為準確模擬血液在壁面附近的流動情況,通過對壁面邊界層進行加密,共設置了11個邊界層．其中,模擬時設置了3種不同網格數量的中空纖維管模型進行數值求解,并選用距入口0.5L截面處的紅細胞徑向體積分數分布結果進行比較,表1給出了不同網格數量模型下的紅細胞在壁面鄰近處的最大體積分數比較．由表1可知,不同網格數量模型下紅細胞在壁面鄰近處的最大體積分數存在差異,網格數量為80 735的模型與網格數量為200 901的模型計算結果相差3.404%,而網格數量為200 901的模型與網格數量為429 600的模型計算結果相差0.488%．因此在保證計算效率和計算精度的前提下,本文選用網格數量為200 901的模型進行數值求解．

表1 不同網格數量模型下紅細胞在距入口0.5L截面處的最大體積分數比較

模擬時血液混合密度為1 059.7 kg/m3．血漿密度為1 030 kg/m3,黏度為0.001 2 kg/(m·s)．紅細胞顆粒的直徑為8 μm,密度為1 096 kg/m3．紅細胞的體積分數在入口處恒有εRBC=0.45,液相(血漿)和固相(紅細胞)在入口處保持相同速度,模擬過程中流速始終保持不變,均為0.2 m/s．出口邊界條件設置為開放出口,相對壓強為0．根據中空纖維管直徑與血液流速,計算得到平均雷諾數Re為17.66,可看作層流狀態．使用計算流體模擬軟件ANSYS FLUENT 2020R2進行模擬,迭代計算采用經典相間耦合Simple算法．

2 結果與討論

2.1 血液兩相流體的速度分布

對中空纖維管內血液兩相流體的速度分布特性進行分析,在流場作用下,入口邊界條件對血液各相沿徑向速度分布的影響較小,模擬發現距入口不同位置處,血液兩相流體沿徑向速度分布曲線基本相同．圖2給出了紅細胞在距入口0.5L截面處的速度分布．由圖2可知,中空纖維管軸心處速度最大,管壁處速度為0,符合層流流動特性．根據圖2的相關結果,圖3給出了紅細胞和血漿在距入口0.5L截面處的速度分布曲線,同時給出了將血液作為單相牛頓流體的速度分布曲線,橫坐標r/r0表示量綱一半徑,r表示紅細胞或血漿與管內中心的距離,r0表示中空纖維管半

圖2 紅細胞在距入口0.5L截面處的速度分布

圖3 中空纖維管內血液兩相流體沿徑向速度

徑．由圖3可知,血液各相的速度呈近似拋物線分布,各相在中空纖維管軸心上的流速大于壁面流速,但血液兩相流模型中血液各相的速度小于血液流動(單相牛頓流體)所形成的拋物線速度．血液表現出一定的非牛頓特性,由于顆粒相紅細胞的存在,顆粒相會拖曳液相(血漿)的流動,減緩其流動速度．其中,徑向各位置處紅細胞流動速度略小于血漿流動速度,該結果與Chandran等[27]采用擴散通量模型得到的研究結論一致．

2.2 紅細胞徑向體積分數分布

對中空纖維管內紅細胞在距入口不同位置處沿徑向的體積分數分布進行統計,共選取了3個截面,其位置與入口的距離分別為0.2L、0.5L、0.8L,圖4給出了紅細胞在距入口0.5L截面處沿徑向的體積分數分布．結果顯示,在徑向上,紅細胞在壁面處受壁面產生的升力作用較大,紅細胞遠離壁面．與入口處紅細胞體積分數εRBC=0.45相比,壁面處紅細胞體積分數減小并小于入口處

圖4 在距入口0.5L截面處紅細胞沿徑向的體積分數分布

紅細胞體積分數,說明紅細胞向遠離壁面的方向發生遷移,與Geislinger等[28]的研究結論相符．在鄰近壁面處(距離壁面一定距離處),紅細胞所受的壁面升力作用與流場剪切梯度對紅細胞產生的升力作用兩者相互平衡[21],此時紅細胞體積分數出現最大值并高于入口處紅細胞體積分數．在管內軸心附近處,紅細胞體積分數略小于入口處紅細胞體積分數．

根據不同截面上紅細胞徑向體積分數分布結果,圖5給出了紅細胞在不同截面上沿徑向的體積分數分布曲線．結果顯示,在管內軸心附近,紅細胞體積分數略小于入口處紅細胞體積分數,而在鄰近壁面處出現最大值,紅細胞沿徑向體積分數分布呈雙峰狀．從鄰近壁面處至管內軸心處,紅細胞體積分數由最大值逐漸減小至局部最小值;而從鄰近壁面處至壁面處,紅細胞體積分數急劇減?。酥?對比各截面上紅細胞在鄰近壁面處的體積分數分布情況,距入口0.2L截面上,紅細胞體積分數的最大值比其他兩個截面上的最大值大,隨著各截面與入口距離的增加,各截面上紅細胞體積分數最大值逐漸減小;而在壁面處,隨著截面與入口距離增加,紅細胞體積分數逐漸增大,距入口0.8L截面上紅細胞體積分數最大．

圖5 紅細胞在不同截面上沿徑向的體積分數分布曲線

2.3 黏度分布

采用血液兩相流模型模擬時,顆粒相紅細胞被看作擬流體,模擬過程采用Carreau-Yasuda黏度模型描述紅細胞黏度．圖6給出了紅細胞在距入口0.5L截面處的黏度分布．紅細胞在壁面處流動速度小,速度梯度大,紅細胞黏度低;在軸心附近處流動速度大,速度梯度小,紅細胞黏度高．根據圖6的黏度分布結果,圖7給出了血液兩相流體沿徑向各位置處的黏度分布曲線．由圖7可知,紅細胞呈現出非牛頓特性,具有一定的剪切稀化特性．紅細胞沿徑向黏度變化最大,由于血漿是一種牛頓流體,其黏度不隨體積分數與流速的變化而變化,模擬過程中黏度始終保持不變．血液主要是由紅細胞與血漿組成的混合溶液,即血液由牛頓流體與非牛頓流體混合而成,整體上表現出一定的非牛頓特性,血液黏度會隨紅細胞體積分數與剪切速率的變化而變化．

圖6 紅細胞在距入口0.5L截面處的黏度分布

圖7 血液各相沿徑向的黏度分布曲線

2.4 壁面條件對紅細胞體積分數分布的影響

由2.2節可知,紅細胞在各截面沿徑向的體積分數分布受到壁面升力作用以及血液流動對其產生的流體力學作用的影響．根據Wu等[8]的相關研究,壁面對囊泡狀物體的升力會受到壁面剪切速率與囊泡變形作用兩者共同的影響．為了分析紅細胞在不同壁面條件下所受升力作用對其體積分數分布的影響,本節共設置了3組紅細胞軸心與壁面距離h=0 μm時的f(s),以模擬血液在中空纖維管不同壁面彈性條件下的流動情況．假設f(s)隨著紅細胞距壁面的距離h呈線性變化,f(s)越小表征壁面所產生的壁面剪切速率對紅細胞的變形作用越小時,即當紅細胞距壁面的距離越小,受到來自壁面對紅細胞產生的非慣性升力越大．圖8給出了不同壁面條件下紅細胞在距入口0.5L截面處沿徑向體積分數分布曲線．由圖8可知,隨著壁面處的剪切速率與壁面對紅細胞的變形作用越來越小(即壁面處f(s)越來越小),紅細胞在鄰近壁面處的體積分數逐漸降低,而在壁面處的體積分數越來越高,表明紅細胞遠離壁面的趨勢減弱．

圖8 不同壁面條件下紅細胞在距入口0.5L截面處沿徑向體積分數分布曲線

為了說明非慣性升力對紅細胞徑向遷移的影響,圖9給出了紅細胞在徑向上的遷移機理示意圖,圖中Fl表示壁面對紅細胞產生的非慣性升力,Fs表示流場剪切梯度對紅細胞產生的升力．如圖所示,當紅細胞在壁面處時,所受的升力Fl>Fs,此時紅細胞遠離壁面;當紅細胞在遠離壁面處于管內中心線附近時,主要受流場剪切梯度產生的升力Fs,此時紅細胞有沿徑向向壁面遷移的趨勢;當紅細胞在緊鄰壁面處(離壁面距離約為5 μm)時,紅細胞受到的非慣性升力與剪切梯度產生的升力兩者達到平衡,即Fl=Fs,此時紅細胞會在鄰近壁面處積聚,在該處紅細胞體積分數出現最大值．當壁面處對應的f(s)減小時,紅細胞在壁面處受到的非慣性升力Fl變小,導致徑向上紅細胞向管內中心的遷移趨勢減弱,紅細胞在鄰近壁面處積聚程度下降,因此紅細胞在壁面處的體積分數增大．該結果與Abkarian等[25]模擬的數值結果吻合．

圖9 紅細胞所受非慣性升力與剪切梯度產生的升力作用機理示意圖

3 結 論

(1)在層流流場作用下,血液中紅細胞與血漿的流速呈近似拋物線分布,各相在軸心處流速均大于在壁面處流速,兩相流體速度小于血液單相流動形成的拋物線速度,沿徑向各位置處的紅細胞流速略大于血漿流速．

(2)紅細胞黏度會隨體積分數與剪切速率變化而變化,壁面處速度梯度大,紅細胞黏度低;軸心附近處速度梯度小,紅細胞黏度高,血液整體上表現出非牛頓特性,呈現剪切稀化特性．

(3)紅細胞在中空纖維管內流動時,壁面對紅細胞產生升力作用,使得紅細胞遠離壁面,在距離壁面一定位置處,紅細胞所受壁面升力作用與流場產生的流體力學作用達到平衡,紅細胞體積分數在該位置處出現最大值．隨著截面與入口距離的增加,該截面處紅細胞體積分數最大值逐漸減小,在壁面處紅細胞體積分數則逐漸增大．整體上看,紅細胞沿徑向體積分數分布呈雙峰狀．

(4)壁面對紅細胞的升力作用會對紅細胞沿徑向體積分數分布產生影響,隨壁面處剪切速率與壁面對紅細胞變形作用的減小,紅細胞在距離壁面一定位置處的體積分數減小,而在壁面處的體積分數增大,紅細胞由壁面向管內中心運動的趨勢減弱．

本文采用兩相流歐拉方法描述了真實條件下血液在微觀尺度下的流動規律,準確設置了血液兩相流體的黏度模型、曳力模型、紅細胞升力方程,相關參數根據表達式使用UDF描述．但由于本文主要采用升力方程來模擬紅細胞在流動時的受力以反映管壁對紅細胞遷移的影響,后續相關研究還可通過設置中空纖維管彈性條件,采用流固耦合等方法來模擬管壁對血液流動的影響．