基于多孔介質(zhì)模型的膜式氧合器內(nèi)部流場分析

葉非華，廖虎，易國斌

（1廣東工業(yè)大學(xué)輕工化工學(xué)院，廣東廣州510006；2廣東順德工業(yè)設(shè)計研究院（廣東順德創(chuàng)新設(shè)計研究院），廣東佛山528300；3武漢科技大學(xué)機械自動化學(xué)院，湖北武漢430081）

膜式氧合器是心肺手術(shù)中心肺機器的重要組成部分，可用于體外生命支持[1-2]。利用計算流體力學(xué)（CFD）對膜式氧合器中的血流進(jìn)行數(shù)值模擬，可以預(yù)測其各項性能（例如壓力分布、灌注動力學(xué)和氣體傳輸性質(zhì)等）[3-8]。由于氧合器纖維膜的不透明性，難以通過實驗直接觀察血液灌注動力學(xué)，相比于傳統(tǒng)實驗，CFD對氧合器性能的預(yù)測更為經(jīng)濟(jì)有效，是研究和優(yōu)化氧合器設(shè)計的重要方法之一[9-10]。本文以一款市售國產(chǎn)膜式氧合器為研究對象，采用多孔介質(zhì)模型，通過數(shù)值模擬氧合器內(nèi)部流場特性，詳細(xì)分析了流體的速度分布、壓力分布、湍流強度分布等，并結(jié)合快速溶血數(shù)值預(yù)估模型計算得到了氧合器的標(biāo)準(zhǔn)溶血指數(shù)NIH。研究結(jié)果有利于產(chǎn)品開發(fā)人員了解膜式氧合器內(nèi)部流體運動特性及其對產(chǎn)品性能的影響，為后續(xù)優(yōu)化設(shè)計氧合器，進(jìn)一步改善其性能提供一定的理論指導(dǎo)。

1 膜式氧合器工作原理

圖1 所示為國產(chǎn)某型號膜式氧合器的工作原理。該膜式氧合器主要包含變溫室、氧合室、中間管道，以及血液、溫水和氧氣的進(jìn)出口等。兩個腔室內(nèi)部均由中空纖維膜填充，其中變溫室用于對血液的溫度進(jìn)行調(diào)節(jié)，氧合室主要起排出血液二氧化碳，增加氧氣含量的作用。膜式氧合器工作時，血液從變溫室下底的進(jìn)血口進(jìn)入氧合器，由下而上從中空纖維管外通過變溫室，與中空纖維膜管內(nèi)從上往下流的溫水在纖維表面進(jìn)行熱交換。變溫后的血液通過中間管道進(jìn)入氧合室，再由上往下從中空纖維膜管外通過氧合室，與中空纖維膜管內(nèi)從上往下流的氣體通過纖維壁上的膜孔進(jìn)行氣體交換。最后變溫、氧合后的血液通過出血口離開氧合器。膜式氧合器避免了血液和水、氧氣直接接觸，生物相容性更好。

圖1 膜式氧合器工作原理

2 數(shù)值模型

2.1 多孔介質(zhì)模型

多孔介質(zhì)模型是用于計算膜式氧合器內(nèi)部流場的主要CFD 模型之一，即假設(shè)氧合器中的纖維束為多孔介質(zhì)，通過滲透率和孔隙率估算纖維束對流體產(chǎn)生的動量損失，這種方法已被證明能夠以有效的方式對氧合器性能進(jìn)行模擬[11-15]。本研究采用Navier-Stokes方程和連續(xù)性方程描述笛卡爾坐標(biāo)系中不可壓縮黏性流體的穩(wěn)態(tài)運動，如式(1)。

式中，ρ為密度；ui為流體速度；P為壓力；τij為應(yīng)力張量；gi為重力；S為源項。

多孔介質(zhì)模型是在Navier-Stokes方程上疊加了一個動量源項S，來描述多孔介質(zhì)區(qū)域中流體的動量損失，i方向上的源項為式(3)。

其中等式右邊第一項代表黏性損失，第二項代表慣性損失。在流量Q＜6.00L/min時，慣性損失可以忽略不計[11,15-16]。因此多孔介質(zhì)區(qū)域的動量損失僅取決于利用達(dá)西定律計算得到的黏性阻力，如式(4)。

式中，1/α為黏性阻力系數(shù)；A為橫截面面積；ΔP為流體穿過纖維束的壓降；μ為流體黏度；Q為流量；L為纖維束長度。

2.2 三維快速溶血數(shù)值預(yù)估模型

傳統(tǒng)的溶血數(shù)值預(yù)估方法是基于拉格朗日粒子軌跡追蹤的方法來進(jìn)行溶血計算[17]，需要在設(shè)備入口處注入粒子并監(jiān)測其在流場中的瞬時切應(yīng)力和停留時間，累積計算沿著運動軌跡的所有時間步中的溶血值，最后求得總?cè)苎禂?shù)[18]。由于多孔介質(zhì)模型無法反應(yīng)單個纖維周圍的局部速度分布和切應(yīng)力分布，使用拉格朗日法計算的溶血值會遠(yuǎn)大于實驗值[19]。新型的溶血數(shù)值預(yù)估方法為三維快速溶血數(shù)值預(yù)估模型，這種模型考慮了流場整體平均效應(yīng)，是對整個穩(wěn)定流場的速度域與切應(yīng)力域進(jìn)行線性場積分計算[20-21]。

2.2.1 流體剪切應(yīng)力計算模型

在湍流狀態(tài)下，血液中的剪切應(yīng)力包括由流體黏性引起的分子切應(yīng)力和湍流引起的雷諾切應(yīng)力。其張量表達(dá)式為式(5)，其中雷諾切應(yīng)力張量為式(6)。

式 中，μτ為 湍 流 黏 度；k為 湍 流 動 能；δij為Kronecker 數(shù)。利用米澤斯屈服準(zhǔn)則（mises yield criterion）將剪切應(yīng)力的張量形式簡化為剪切應(yīng)力標(biāo)量形式，如式(7)。

2.2.2 溶血數(shù)值預(yù)估模型

Giersiepen 等[17]通過大量溶血實驗，歸納了溶血值D與剪切應(yīng)力τ以及暴露時間t的冪函數(shù)關(guān)系式，如式(8)。

式中，D為溶血值，是衡量溶血量的參數(shù)；ΔHb代表被損傷的血紅蛋白濃度，g/L；Hb代表總血紅蛋白濃度（Hb=140g/L）；τ為剪切應(yīng)力；t為暴露時間；C、a、b是通過試驗數(shù)據(jù)回歸得到的經(jīng)驗常數(shù)。Giersiepen 等通過試驗得到的一組常數(shù)值：C=3.62×10-7、a=2.416、b=0.785，被眾多文獻(xiàn)所引用，過往研究也證明了這組常數(shù)應(yīng)用于溶血指數(shù)計算時的準(zhǔn)確性。因此，在本研究中也采用這組常數(shù)數(shù)值進(jìn)行溶血指數(shù)計算。

Garon 等[20]基于雙曲輸運方程開發(fā)了一種三維快速溶血數(shù)值模擬的方法來預(yù)測溶血。雙曲輸運方程如式(9)。

式中，V為速度矢量；DI為線性溶血指數(shù)；源項σ為單位時間溶血破壞率，具體計算如式(10)、式(11)。

在計算獲得一個穩(wěn)定流場后，膜式氧合器內(nèi)部線性平均溶血指數(shù)為式(12)。

式中，Q為體積流量。通過指數(shù)換算得到溶血值，如式(13)。

將溶血值轉(zhuǎn)換為最終的標(biāo)準(zhǔn)溶血參數(shù)值（NIH），如式(14)。

3 實驗方案

通過實驗測量膜式氧合器不同流量下變溫室和氧合室的壓降值，計算纖維束的黏性阻力。開展氧合器內(nèi)部流場數(shù)值計算，驗證模擬結(jié)果的準(zhǔn)確性。實驗流程見圖2，本實驗使用市售國產(chǎn)某型號成人膜式氧合器，以水作為流體介質(zhì)，通過離心泵將水從進(jìn)血口導(dǎo)入氧合器，測量離心泵在不同轉(zhuǎn)速時氧合器進(jìn)血口，出血口的流量值及壓強值，其中流量值與壓強值通過相應(yīng)的傳感器（速度傳感器為SONOFLOW CO.56，壓力傳感器為NORa 有創(chuàng)壓力傳感器）測量得到。之后將氧合器的變溫室和氧合室分離，并分別測量這兩部分在不同流量下的壓降值。

圖2 實驗流程

利用matlab軟件將實驗數(shù)據(jù)Q-ΔP線性插值擬合，如圖3 所示，得到變溫室和氧合室區(qū)域的ΔP/Q值，將其分別代入式(4)即可求得兩區(qū)域各自的黏性阻力：氧合室黏性阻力1/α=9.210×108m-2；變溫室黏性阻力1/α=1.069×109m-2。

圖3 氧合室和變溫室壓降-流量擬合曲線

4 仿真計算設(shè)置

使用商業(yè)CAD 軟件SolidWorks 2018 繪制膜式氧合器的流道三維模型。計算模型由兩個計算域組成：環(huán)形纖維束區(qū)域（多孔介質(zhì)域）和外部流體域（纖維束與殼體之間的間隙、入口管道、出口管道以及中間管道）。為讓流動充分發(fā)展和避免回流，加長了出口管道，如圖4所示，圖5給出了氧合器各部分尺寸。

圖4 膜式氧合器幾何模型示意圖

圖5 氧合器結(jié)構(gòu)尺寸（單位：mm）

由于流體域結(jié)構(gòu)復(fù)雜，多孔介質(zhì)域與流體域皆采用非結(jié)構(gòu)網(wǎng)格劃分，兩計算域交界面設(shè)置為Interface。采用商業(yè)軟件ICEM17.0 劃分流場網(wǎng)格。為了驗證網(wǎng)格無關(guān)性以及滲透率的正確性，分別以三種不同網(wǎng)格數(shù)（190萬、296萬、420萬）計算入口流量為1.65～3.90L/min，流體介質(zhì)為水時的進(jìn)出口壓降值，在低雷諾數(shù)（Re＜6000）下采用RNGk-ε湍流模型模擬，高雷諾數(shù)（Re＞6000）下使用標(biāo)準(zhǔn)k-ε湍流模型模擬，并與實驗結(jié)果進(jìn)行比較。

由圖6可見，網(wǎng)格數(shù)逐漸增大對計算結(jié)果影響較小，296 萬網(wǎng)格與420 萬網(wǎng)格的計算結(jié)果基本一致，且與實驗測試結(jié)果之間的誤差較小，這證明了計算的準(zhǔn)確性。綜合考慮計算耗時和求解結(jié)果精度，在本論文中確定以296萬網(wǎng)格數(shù)進(jìn)行研究。

圖6 網(wǎng)格無關(guān)性驗證

圖7 管道軸面速度云圖

利用商業(yè)軟件FLUENT17.0 進(jìn)行流場模擬計算。入口采用速度入口邊界條件（對應(yīng)入口流量為1.65～6.00L/min），出口為靜壓0Pa 的壓力出口，操作壓力為默認(rèn)的大氣壓101325Pa。湍流模型選擇RNGk-ε模型，采用SIMPLE 算法進(jìn)行壓力-速度耦合，離散格式采用二階迎風(fēng)模式。流體介質(zhì)血液被視為不可壓縮的牛頓流體，具有恒定密度1055kg/m3和黏度0.0035kg/(m·s)，所有壁面均為無滑移固壁邊界，將中空纖維膜作為各向同性多孔介質(zhì)處理。當(dāng)所有速度和湍流殘差小于10-5時，計算視為收斂。

5 計算結(jié)果及分析

5.1 流場分析

為研究膜式氧合器內(nèi)部流體運動特性，利用速度云圖、矢量圖、流線圖進(jìn)行同步分析。圖7顯示了入口流量為4.50L/min 時的速度分布云圖，其中A、B、C三個截面分別為膜式氧合器入口管道、中間管道和出口管道軸中間面。圖8為全局速度流線圖。

圖8 全局速度流線圖

從圖7 和圖8 可見，由于入口和流動腔室設(shè)計，流體在進(jìn)入變溫室的纖維束之前，軸向射流受到纖維束區(qū)域的流動阻力，先分流到流體域中，之后流體在膜前后壓差的作用下突破流動阻力徑向穿過纖維束區(qū)域，向上匯集到中間管道，由中間管道進(jìn)入氧合室，最后在出口處重新匯集流出氧合器。提取不同流量下中間管道和出口管道的軸切面平均速度變化，如圖9和圖10所示，在中間管道起點處由于流體的匯集，流速激增。隨著流體的繼續(xù)流動，平均流速逐漸降低，直至到達(dá)管道中部流速才保持穩(wěn)定，但在管道終點處由于流體的分流，使得切面平均流速有所增加。出口管道流速分布與中間管道類似。

圖9 中間管道軸切面平均速度變化

圖10 出口管道軸切面平均速度變化

為對比不同流量下流場分布情況，進(jìn)一步研究了膜式氧合器徑向中間截面的速度云圖和出口管道區(qū)域的速度矢量圖。如圖11 所示，流量在1.65～6.00L/min 時，隨著入口流量的增加，膜式氧合器內(nèi)部速度梯度分布形式基本保持不變。從圖12 可見，在出口管道上壁面出現(xiàn)較大渦旋，其渦流區(qū)約占流道的1/6～1/3，過大的渦旋產(chǎn)生了血液再循環(huán)流動現(xiàn)象，會導(dǎo)致高剪切應(yīng)力的產(chǎn)生，同時增加了血紅細(xì)胞在氧合器中的暴露時間，加劇了血紅細(xì)胞破壞的可能性。隨著流量的增加，渦旋區(qū)域面積逐漸增大，且渦旋中心逐漸向出口方向移動，回流區(qū)域面積增大，使得倒流量逐漸增加。

圖11 氧合器徑向中間截面的速度云圖

圖12 出口管道速度矢量圖

圖13 氧合器壓力分布

圖13顯示了入口流速為0.955m/s（流量4.50L/min）時全局和徑向中間截面的總壓云圖。由于入口管道、中間管道和出口管道的非軸對稱設(shè)計，膜式氧合器內(nèi)部壓力分布不同于Gage 等[11]和Pelosi 等[16]預(yù)測的同心均勻模式，而是呈傾斜狀態(tài)逐漸減小。出入口總壓降ΔP 約為36677Pa，由于纖維束的高黏性阻力，大部分壓力損失位于纖維束內(nèi)，其中53.3%位于氧合室，42.6%位于變溫室。入口管道和中間管道對壓降的影響可忽略不計，出口管道提供的壓力損失約占4.1%，且由于血液進(jìn)入出口管道時流通面積有較大變化，在出口管道內(nèi)產(chǎn)生了負(fù)壓區(qū)，負(fù)壓區(qū)與渦旋區(qū)一致。

圖14 壁面剪切應(yīng)力分布云圖

5.2 血液損傷分析

本文將三維快速溶血數(shù)值預(yù)估模型通過用戶自定義函數(shù)（UDF）編譯到FLUENT 17.0中，在獲得穩(wěn)定的流場后計算膜式氧合器的標(biāo)準(zhǔn)溶血指標(biāo)，同時將壁面剪切應(yīng)力、湍流強度和紅細(xì)胞停留時間作為輔助參數(shù)，用來分析氧合器中血液的損傷程度。

如圖14 所示，流量為4.50L/min 時的最大壁面剪切應(yīng)力約為118Pa，位于氧合室入口流體域外壁面。高壁面剪切應(yīng)力區(qū)域存在于氧合室和變溫室的入口及出口部分，與高流速區(qū)對應(yīng)，而其他流域的剪切應(yīng)力較低，不太可能產(chǎn)生過量流動引起血液損傷。圖15 顯示流體域中湍流強度較大，主要出現(xiàn)在分流或回流處，而多孔介質(zhì)域中也存在小區(qū)域的湍流強度較大區(qū)，位于與流體域的相連之處。氧合器纖維束內(nèi)的血流速度較小，施加在血液上的剪切應(yīng)力低，對血細(xì)胞破壞小。而氧合器流體域，尤其是氧合室和變溫室的入口及出口部分，皆出現(xiàn)較高的速度、壁面剪切應(yīng)力和湍流強度分布，這些區(qū)域是血細(xì)胞容易遭到破壞的高發(fā)區(qū)域。

圖15 湍流強度分布云圖

圖16 流量1.65～6.00L/min時紅細(xì)胞停留時間

圖17 較大流量時紅細(xì)胞停留時間分布

計算收斂后，將直徑為7μm的球形顆粒代表紅細(xì)胞，利用離散相模型（DPM）計算粒子從氧合器入口到出口所經(jīng)歷的時間。圖16 顯示隨著流量的增加，紅細(xì)胞的停留時間及其標(biāo)準(zhǔn)差逐漸減小。圖17 顯示在較大流量時，紅細(xì)胞在氧合器中的停留時間分布趨于平穩(wěn)。當(dāng)流量為3.50L/min 時，紅細(xì)胞平均停留時間依然達(dá)到9.13s，停留時間過長。這是由于該膜式氧合器采用的是分離式單通道結(jié)構(gòu)設(shè)計，血液預(yù)存量大，通過氧合器流程長。此外流場中形成的渦旋、回流等不規(guī)則流動現(xiàn)象也是造成紅細(xì)胞停留時間變長的原因。紅細(xì)胞在氧合器停留時間長，血液損傷的風(fēng)險高，不利于氧合器的長期使用。

圖18 標(biāo)準(zhǔn)溶血指標(biāo)NIH隨流量的變化

基于三維快速溶血數(shù)值預(yù)估模型，對流量在1.65～6.00L/min下膜式氧合器內(nèi)部流場進(jìn)行溶血估算，得到該氧合器的標(biāo)準(zhǔn)溶血指標(biāo)NIH為0.0084～0.0835g/100L（見圖18），結(jié)合人體生理允許的最大溶血指標(biāo)0.1g/100L，說明該氧合器的溶血性能滿足使用要求。

6 結(jié)論

本文基于各向同性多孔介質(zhì)模型，通過計算流體力學(xué)對膜式氧合器中的血流進(jìn)行數(shù)值模擬，采用RNGk-ε湍流模型和溶血數(shù)值預(yù)估模型，對不同流量下氧合器內(nèi)部流體運動特性及其對性能的影響進(jìn)行了研究，得出以下結(jié)論。

（1）低流量1.65～3.00L/min 下，各向同性多孔介質(zhì)模型能準(zhǔn)確模擬氧合器內(nèi)部血液流動，計算值與實驗值相差較小，但隨著流量的增加，兩者之間的偏差逐漸增大。

（2）膜式氧合器內(nèi)部壓力分布呈傾斜狀態(tài)且逐漸減小，大部分壓力損失位于纖維束內(nèi)，其中53.3%位于氧合室，42.6%位于變溫室。

（3）氧合器血液的入口及出口位置是血液損傷的高發(fā)區(qū)域，同時該氧合器因分離式單通道結(jié)構(gòu)設(shè)計等原因，導(dǎo)致紅細(xì)胞停留時間偏長，不利于氧合器的長期使用。

（4）流量為1.65～6.00L/min 時，該氧合器標(biāo)準(zhǔn)溶血指標(biāo)NIH 為0.0084～0.0835g/100L，滿足人體生理允許的使用范圍。