基于多孔介質模型的膜式氧合器內部流場分析

葉非華，廖虎，易國斌

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

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

1 膜式氧合器工作原理

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

圖1 膜式氧合器工作原理

2 數值模型

2.1 多孔介質模型

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

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

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

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

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

2.2 三維快速溶血數值預估模型

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

2.2.1 流體剪切應力計算模型

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

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

2.2.2 溶血數值預估模型

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

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

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

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

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

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

將溶血值轉換為最終的標準溶血參數值（NIH），如式(14)。

3 實驗方案

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

圖2 實驗流程

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

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

4 仿真計算設置

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

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

圖5 氧合器結構尺寸（單位：mm）

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

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

圖6 網格無關性驗證

圖7 管道軸面速度云圖

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

5 計算結果及分析

5.1 流場分析

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

圖8 全局速度流線圖

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

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

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

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

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

圖12 出口管道速度矢量圖

圖13 氧合器壓力分布

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

圖14 壁面剪切應力分布云圖

5.2 血液損傷分析

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

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

圖15 湍流強度分布云圖

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

圖17 較大流量時紅細胞停留時間分布

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

圖18 標準溶血指標NIH隨流量的變化

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

6 結論

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

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

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

（3）氧合器血液的入口及出口位置是血液損傷的高發區域，同時該氧合器因分離式單通道結構設計等原因，導致紅細胞停留時間偏長，不利于氧合器的長期使用。

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