基于多相流的軸流血泵流場分析及溶血指數預測

王帶領，譚建平，喻哲欽

?

基于多相流的軸流血泵流場分析及溶血指數預測

王帶領，譚建平，喻哲欽

(中南大學 機電工程學院，湖南 長沙，410083)

為研究血泵內部血細胞分布規律及溶血預測方法，以自制軸流血泵為例，應用多相流分析方法，采用多重參考坐標系法(MRF)等技術建立數值分析模型，并通過體外循環實驗驗證模型的有效性。進一步分析血泵內部血細胞濃度、速度、壓力等的分布規律，得到血細胞破壞區域和一般規律。根據優化的溶血模型對血泵的溶血性能進行評估，在此基礎上提出溶血實驗指標標準溶血指標(NIH)與溶血預測值之間的轉換關系。研究結果表明：血細胞在血泵內部不是均勻分布，葉輪處有明顯分離現象；血泵內部剪切應力在 200 Pa 以下的區域的體積分數約為98.6%，高剪切應力主要位于葉輪頂部和外殼的間隙處；在剪切應力為常數時，優化的溶血預測模型和Giersiepen實驗得到的溶血模型相符合；采用優化模型計算得到血泵平均溶血預測值0.005 7，具有較好的溶血性能。

血泵；多相流；標準溶血指標；溶血模型

血泵已經成為幫助心臟功能障礙患者的一種有效方法。血泵經歷了從容積式到懸浮技術的變革，其在驅動控制、生物相容性、等方面都有了極大的發展[1]。在血泵的設計開發中，計算流體力學(CFD)得到了普遍應用。基于血泵虛擬樣機的血液動力學性能數值分析能夠極大地縮短設計周期，節約成本[2?3]，但同時伴隨著一些問題，其中溶血和血栓尤為突出[4?5]：因此，研究數值模擬血泵內部流場特性和溶血預測算法對優化血泵設計具有重要意義，國內外學者為尋找一種有效的溶血預測方法進行了廣泛研究[6]。YANO等[7]采用拉格朗日粒子追蹤法對血泵的溶血性能進行了數值分析，并與實驗結果進行對比，驗證了模型的合理性。GARON等[8]通過歐拉法建立了快速溶血模型，通過體外溶血實驗對模型進行了驗證。TASKIN等[9]通過設計特定的溶血實驗設備，比較了拉格朗日幾種冪律模型和歐拉模型的區別，實驗結果表明歐拉模型溶血預測誤差較大。LI等[5]使用多相流模型對軸流泵流場進行了仿真分析，并提出了基于剪切應力的溶血估計算法，通過實驗進行驗證。周冰晶等[10?11]通過拉格朗日粒子追蹤法對自主研發的血泵進行了溶血預測，從而對血泵結構優化提供了理論指導。壽宸等[12]通過比較4種不同葉輪形式的離心血泵仿真結果，并采用快速溶血模型進行溶血預測，發現螺旋線葉輪性能較優。以往的CFD分析中大都把血液近似為單一牛頓流體進行數值計算。采用單一流體一方面不能得到血細胞在血泵內部的分布規律；另一方面，簡化的血液模型與實際結果有差異，不能反映血細胞對流場的影響，難以揭示真實血液的流動規律。鑒于血液組成的復雜性，本文作者通過采用多相流模型，兼顧血細胞與血漿的相互作用，使模擬流動與實際流動更接近，分析血細胞體積濃度、速度，壓力等的分布規律。并在原有的拉格朗日粒子追蹤法的基礎上，考慮加載歷史對溶血的影響，提出優化溶血模型，并對兩者進行比較分析。人們普遍采用標準溶血指標(NIH)對血泵的溶血特性進行評估[5]，考慮到NIH與溶血預測值這2個指標的意義和計算方法的區別，提出標準溶血指標(NIH)與溶血預測值之間的轉換關系。

1 模型及算法

1.1 血泵模型及網格劃分

課題組自主設計了一個由磁力驅動的微型軸流血泵，包含前導、葉輪、后導、永磁體，軸承等結構。前導葉和后導葉有導流作用，葉輪采用螺旋結構，其高速旋轉促使血液運動并獲得壓能，從而起到輔助人體心臟泵血的功能；永磁體放置在葉輪部分，分為N和S兩極，通過與由單片機控制的三齒槽定子驅動線圈的相互作用使血泵的葉輪持續轉動[13]。血泵進出口直徑為16.4 mm，葉輪輪緣直徑為16.0 mm；模型總長81.0 mm。圖1所示為血泵結構示意圖。

將血泵幾何文件導入到ICEM中劃分網格。由于流體域各部分運動狀況不同，把模型劃分為3部分即前導葉區，葉輪區，后導葉區，使用多重參考坐標系法(MRF)，分別建立動靜坐標系，并在動靜區域創建交界面。因為葉輪、導葉結構較復雜，使用非結構化網格技術進行處理，共得到2 410 710個網格單元。圖2所示為網格劃分示意圖。

圖1 血泵結構示意圖

圖2 網格劃分示意圖

1.2 歐拉多相流模型

紅細胞是血液的重要構成成分，約占血液體積的40%。由于血液成分復雜，因此，選用多相流進行模擬。在多相流模型中，歐拉模型能較好地模擬固相體積分數超過10%的情況，且與其他模型相比具有較高的求解精度，因此，選用歐拉模型進行仿真計算。

將血漿定義為液相，其密度為1.03 t/m3，黏度為1.60 mPa?s；紅細胞定義為固相，其密度為1.09 t/m3，黏度為6 mPa?s[14]。人體中成熟的紅細胞為雙凹面圓盤狀，直徑為7~9 μm，為了簡化模型，固相(紅細胞)假設為球形，直徑均勻，計算中取值為8 μm。

1.3 邊界條件及算法

為了滿足人體基本的生理需求，血泵設計要求進出口壓差約為13.3 kPa，流量為5 L/min。入口使用速度進口邊界，根據流量和截面積，計算出進口速度為0.4 m/s。出口使用壓力出口條件，壓力為13.3 kPa。血泵葉輪區域的壁面設置為旋轉邊界，其轉速與葉輪區的轉速相同，該血泵設計工況為8 000 r/min。其余壁面定義為無滑移固壁邊界。進口紅細胞體積分數設置為40%。

模型采用標準的?湍流模型，使用simple算法，離散格式設置為二階迎風。

1.4 仿真模型驗證

為驗證仿真模型的可靠性，搭建血泵體外循環實驗系統，設計工況下進行體外循環性能實驗。該實驗系統主要由儲液槽、血泵、驅動線圈、壓力傳感器、渦輪流量計、電動閥及各支架和管道等組成。壓力傳感器分別測量血泵進出口壓力，渦輪流量計測量血泵輸出流量。實驗中循環介質采用體積分數為37%的甘油和水的混合溶液，其在常溫條件下，動力黏度和血液的匹配。設定葉輪轉速為8 000 r/min，并調節血泵進出口壓差為13.3 kPa。

在上述工況下進行血泵性能實驗。表1所示為體外循環實驗中實驗流量與仿真結果之間對比。由表1可知：在8 000 r/min下血泵的平均流量約為4.8 L/min，基本滿足人體生理需求；其中測量值與仿真值之間的最大誤差為4.2%，均小于5.0%，誤差在允許范圍內。實驗結果在一定程度上驗證了血泵數值仿真模型的有效性，為后續血泵流場分布規律的研究提供了保證。

表1 實驗流量與仿真流量比較

2 仿真結果分析

2.1 固相體積分數分布

按照上述參數和算法設置，對血泵內部兩相分布情況進行數值分析。圖3所示為血細胞體積分布云圖。從圖3可見：血細胞在血泵內部并非均勻分布，體積分數為35%~42%。其中，前導葉區域由于沒有運動部件，流動較平穩，血細胞體積分數分布均勻，約為40%。隨著進一步流動，在葉輪區域，由于葉輪的高速旋轉和因離心力與慣性的影響，血細胞的體積分數逐漸開始發生變化，其中，在葉輪前端部分差異不明顯，隨著輪轂增大，血細胞逐漸偏向葉輪輪緣處。特別是在葉輪分流葉片處，輪轂處附近血細胞體積分數減小，葉輪工作面及輪緣處血細胞濃度有明顯增大，達到42%。當血細胞運動到后導葉區域，由于導葉的導流作用，切向速度減小，軸向速度增加，流動再次趨于平穩，血細胞濃度分布重新變得均勻。由于血細胞和血漿密度接近，并且血細胞直徑較小，從而對流場的影響有限，各部分體積分數差異不是很大，僅在葉輪區域存在小程度的分離現象。

(a) 血泵內部體積分數；(b) 軸截面體積分數

2.2 速度與壓力分布

圖4所示為軸截面單相與兩相速度分布。由圖4可見：兩者速度分布基本相同；在血泵進口位置，速度較小；在葉輪部分，由于軸流血泵的葉輪具有極高的轉速，從輪轂沿著徑向速度逐漸增大，特別是在葉輪邊緣具有較大的速度梯度，對細胞破壞作用大；在出口部位，由于導葉的作用，速度較均勻。多相流下液體的最大速度略大。這是因為在葉輪邊緣血細胞體積濃度較高，血細胞與血漿的相互作用，血漿速度增大。云忠等[15]發現紅細胞其垂直撞擊速度超過6 m/s，有發生破碎的可能。從圖4可見：血液最大流速達 6.95 m/s，特別是在輪緣與外殼之間的區域，紅細胞有破碎的危險。

(a) 兩相流速度云圖；(b) 單相流速度云圖

圖5所示為血泵軸向截面總壓分布。由圖5可見：入口處壓力較小；流道內沿著流動方向，壓力逐漸增大；分流葉片處具有較大的壓力梯度，對紅細胞的損傷較嚴重；在出口區域，由于后導葉的作用，速度轉換為壓力能，壓力分布均勻，進出口壓差約為 13.3 kPa，能滿足人體基本需求。

圖5 軸截面總壓分布云圖

3 溶血預測模型

3.1 通用溶血預測模型

GIERSIEPEN等[16]通過實驗得到了溶血公式。該模型指出溶血程度與流場剪切應力、紅細胞在剪切應力中的暴露時間呈冪函數關系：

血泵內的實際流動為不規則的湍流運動，剪切應力為黏性剪切應力與湍流剪切應力的合成。黏性剪切應力與湍流剪切應力計算公式見文獻[10]。剪切應力標量可以表示如下：

溶血算法主要有歐拉法和拉格朗日法，均基于式(1)所示的溶血實驗公式進行推導。其中拉格朗日算法核心思想是：當單個紅細胞流過血泵時，由于運動軌跡上各點的剪切應力是不斷變化的，將粒子軌跡線分成個微小段，在任意微小時間段(=1~)內，計算溶血值，進行疊加。

YANO等[7]采用溶血算法計算任意小段時間內的溶血值：

跡線上某一時刻的累計溶血值D,i可以通過D,i?1遞推計算得到，遞推公式如下：

經計算，可以得到整條跡線的溶血值，對條跡線進行平均處理，可以得到血泵溶血預測值。

3.2 優化溶血模型

式(1)是基于切應力不變的情況下得出的，在剪切應力不變情況下，據式(1)，任意時間步的溶血貢獻值d,i如下：

比較式(3)與(5)，YANO等[7]將每個時間步內用Δ0.785代替Δ(0.785)，因為0.785不符合線性疊加且一階導數單調遞減，即

導致溶血預測值偏大。考慮到方程中時間非線性以及加載歷史的影響對算法進行改進，引入等效時間。在第1個時間段d1內，若它的溶血預測值d,1為

則在時間步d2，需要重新設定時間起點。為此，引入等效時間deq：

重復以上步驟，將下一個時間步重新設定的時間起點。經計算，可以得到整條跡線的溶血預測值。對條跡線進行平均處理，可以得到標準溶血預測值。假設細胞在運動過程中受到的剪切應力為常數，加載時間為0.5 s，切應力為100 Pa。圖6所示為常剪切應力下不同模型的溶血預測。從圖6可見：優化算法在剪切應力恒定時與GIERSIEPEN公式所得計算結果相吻合，溶血值為1.4×10?2，YANO等[7]的算法溶血值為2.0×10?2，與公式值相比偏大，分析結果一致。

3.3 溶血預測結果

標量切應力分布影響血泵的溶血性能，特別是切應力大于200 Pa的區域對紅細胞有較大破壞作用[17]。圖7所示為血泵內部剪切應力云圖，圖8所示為不同剪切應力區域在總體區域中所占的比例。從圖7可見：內部最大剪切應力約為877 Pa，主要分布在葉輪和后導葉入口區域。從圖8可見：剪切應力在 200 Pa 以下的區域所占比例約為98.6%，超過200 Pa的區域約為1.4%，可見使紅細胞產生較大損傷的高剪切應力區域較小，這些部分主要位于外殼和葉輪頂部的間隙以及后導葉入口處，此處血細胞的濃度也較高，血細胞可能因為剪切應力較大，與葉輪劇烈碰撞而破碎。

在FLUENT中通過UDF，根據切應力標量公式定義函數。導出每條跡線上切應力與時間的對應數據文件；利用MATLAB讀取數據文件并通過程序進行數據處理從而得到每條跡線的溶血預測值，進一步得到平均溶血預測值。圖9所示為其中1條流線的剪切應力和溶血隨時間變化規律。從圖9可見：紅細胞在開始階段即前導葉區域，剪切應力較小，溶血預測值累計值基本為0；當紅細胞進入葉輪區域時，速度梯度較大，剪切應力出現高峰區，溶血預測值急劇增加；當細胞經過后導葉時，流動逐漸平穩，剪切應力減小，溶血預測值基本保持不變。這2種算法所得溶血預測值變化趨勢相同，YANO等[7]的算法溶血預測值較高(1.0×10?3)，優化算法溶血預測值為2.8×10?4。

(a) 常剪切應力加載；(b) 不同算法溶血值

圖7 剪切應力分布云圖

圖8 不同剪切應力區域占比

(a) 剪切應力變化；(b) 溶血變化

選取117條流線，利用優化算法對117條流線溶血值進行平均，得到血泵的平均溶血值為0.005 7，溶血性能較好，可為血泵溶血性能的數值預測提供參考。

3.4 標準溶血指標與預測值的轉換

國際上，標準溶血指標(NIH)被廣泛用來評估體外輔助裝置對紅細胞的損傷程度[5]。計算式可表示為

式中：ct為紅細胞壓積百分數，%；ΔHB為測試間隔時間內血漿游離血紅蛋白的增量，g/L；為總循環體積，L；為輔助裝置輸出流量，L/min；為測試時間間隔，min。

溶血預測值為游離血紅蛋白與血液中全部血紅蛋白的比值即溶血率，為量綱一量。兩者不能直接進行比較，需要通過公式建立聯系來驗證模型。

根據上述溶血模型單次循環溶血預測值，進行次循環的溶血估計值為D，其中0=0。D可以通過下式計算得到：

進行次循環的血紅蛋白增量值ΔHB可以計算如下：

式中：GB為血紅蛋白質量濃度，g/L。

通過ΔHB，利用式(10)可以計算得到標準溶血指數預測值IH，與實驗結果NIH進行對比，從而評估溶血模型。

4 結論

1) 多相流模型能夠揭示血泵內部血細胞分布的一般規律。血細胞在血泵內不是均勻分布的，進出口區域分布較均勻，葉輪處存在小程度分離現象。由于血細胞體積較小，且其密度與血漿的密度接近，離析現象不嚴重。

2) 血泵內部切應力超過200 Pa的區域約為1.4%，這些區域主要位于外殼和葉輪頂部的間隙與后導葉入口處，高剪切應力和劇烈的碰撞可能導致血細胞破碎。

3) 優化的溶血模型在剪切應力為常數時，與GIERSIEPEN的溶血實驗公式相吻合。YANO等的算法溶血預測值偏高，改進模型的計算值較優。

4) 建立了實驗指標IH與溶血預測值之間的數學聯系，為評估溶血模型提供了理論方法。

[1] 韓偉, 韓冰雪, 王漢義, 等. 基于CFD螺旋離心式血泵與離心式血泵內部流場的數值分析[J]. 蘭州理工大學學報, 2014, 40(5): 56?59. HAN Wei, HAN Bingxue, WANG Hanyi, et al. Number analysis of internal flow field in screw centrifugal blood pump and centrifugal blood pump based on CFD[J]. Journal of Lanzhou University of Technology, 2014, 40(5): 56?59.

[2] FRASER K H, TASKIN M E, GRIFFITH B P, et al. The use of computational fluid dynamics in the development of ventricular assist devices[J]. Medical Engineering & Physics, 2011, 33(3): 263?280.

[3] COWGER J A, ROMANO M A, SHAH P, et al. Hemolytic:a harbinger of adverse outcome after left ventricular assist device implant[J]. Heart Lung Transplant, 2014, 33(1): 35?43.

[4] TASKIN M E, FRASER K H, ZHANG T, et al. Computational characterization of flow and hemolytic performance of the UltraMag blood pump for circulatory support[J]. Artificial Organs, 2010, 34(12): 1099?1113.

[5] LI Tieyan, YE Liang, FW Hong, et al. The simulation of multiphase flow field in implantable blood pump and analysis of hemolytic capability[J]. Journal of Hydrodynamics: Ser B, 2013, 25(4): 606?615.

[6] GOUBERGRITS L. Numerical modeling of blood damage: current status, challenges and future prospects[J]. Expert Review of Medical Devices, 2006, 3(5): 527?531.

[7] YANO T, SEKINE K, MITOH A, et al. An estimation method of hemolysis within an axial flow blood pump by computational fluid dynamics analysis[J]. Artificial Organs, 2003, 27(10): 920?925.

[8] GARON A, FARINAS M I. Fast three-dimensional numerical hemolysis approximation[J]. Artificial Organs, 2004, 28(11): 1016?1025.

[9] TASKIN M E, FRASER K H, ZHANG T, et al. Evaluation of eulerian and lagrangian models for hemolysis estimation[J]. ASAIO Journal, 2012, 58(4): 363?372.

[10] 周冰晶, 荊騰, 王芳群, 等. 兩級軸流血泵基于血液損傷的數值分析[J]. 生物醫學工程學雜志, 2016, 33(4): 686?690. ZHOU Bingjing, JING Teng, WANG Fangqun, et al. Numerical analysis of two-stage axial blood pump based on blood damage[J]. Journal of Biomedical Engineering, 2013, 25(4): 606?615.

[11] 吳華春, 龔高, 王子彥, 等. 軸流式磁懸浮血泵流場數值模擬及溶血預測[J]. 中國機械工程, 2013, 24(3): 399?403. WU Huachun, GONG Gao, WANG Ziyan, et al. Flow field numerical simulation and hemolytic for axial flow maglev blood pump[J]. China Mechanical Engineering, 2013, 24(3): 339?403.

[12] 壽宸, 郭勇君, 蘇磊, 等. 基于快速溶血預估模型的離心血泵葉輪特性數值分析[J]. 生物醫學工程學雜志, 2014, 31(6): 1260?1264. SHOU Chen, GUO Yongjun, SU Lei, et al. Numerical assessment of impeller features of centrifugal blood pump based on fast hemolysis approximation model[J]. Journal of Biomedical Engineering, 2014, 31(6): 1260?1264.

[13] 譚卓, 譚建平, 劉云龍, 等. 大間隙磁力驅動軸流式血泵的電磁特性[J]. 中南大學學報(自然科學版), 2015, 46(1): 99?106. TAN Zhuo, TAN Jianping, LIU Yunlong, et al. Electromagnetic characteristics of large gap magnetic driving axial flow blood pump[J]. Journal of Central South University (Science and Technology), 2015, 46(1): 99?106.

[14] 云忠, 向闖, 蔡超, 等. 邊界振動流場對單個紅細胞損傷影響分析[J]. 生物醫學工程學雜志, 2016, 33(1): 78?82. YUN Zhong, XIANG Chuang, CAI Chao, et al. Analysis of influence on single erythrocyte injury caused by oscillating boundary flow field[J]. Journal of Biomedical Engineering, 2016, 33(1): 78?82.

[15] 云忠, 譚建平. 基于血液撞擊損傷機理的高速螺旋血泵仿真分析[J]. 中南大學學報(自然科學版), 2008, 39(1): 135?142. YUN Zhong, TAN Jianping. Simulation analysis of high-speed spiral blood pump based on impacting injure principle of blood[J]. Journal of Central South University (Science and Technology), 2015, 46(1): 99?106.

[16] GIERSIEPEN M, WURZINGER L, OPITZ R, et al. Estimation of shear stress-related blood damage in heart valve prostheses: in vitro comparison of 25 aortic valves[J]. The International Journal of Artificial Organs, 1990, 13(5): 300?306.

[17] 吳廣輝, 藺嫦燕, 張錫文, 等. 磁驅動血泵溶血分析[J]. 中國生物醫學工程學報, 2010, 29(3): 473?476. WU Guanghui, LIN Changyan, ZHANG Xiwen, et al. Hemolysis test of an axial flow ventricular assistant device[J]. Chinese Journal of Biomedical Engineering, 2010, 29(3): 473?476.

(編輯 陳燦華)

Flow field analysis and hemolytic prediction of axial blood pump based on multiphase flow

WANG Dailing, TAN Jianping, YU Zheqin

(School of Mechanical and Electrical Engineering, Central South University, Changsha 410083, China)

In order to study the distribution of blood cells in blood pump and the method of hemolytic prediction, the self-made axial flow blood pump was taken as an example, analytical method of multiphase flow and multiple reference frame method (MRF) technology were used, numerical analytical model was established. The validity of the model was verified by cardiopulmonary bypass. The internal distribution of blood cell concentration, velocity, pressure and so on in blood pump was analyzed. In order to get the blood cell damage area and the general law, the improved particle tracking method was used to evaluate the hemolytic performance of blood pump, and the conversion relationship between the standard hemolytic index (NIH) and hemolytic prediction value was obtained. The results show that the blood cells in the flow channel are not evenly distributed, and there is separation phenomenon at the impeller blood pump. When the internal stress in the area is below 200 Pa, volume fraction is about 98.6%, high shear stress is mainly located in the gap between the impeller’s top and the outer shell. Under the condition that shear stress is a constant, the improved hemolytic prediction model can match the hemolytic formula obtained by Giersiepen test well. The modified hemolytic model is used to estimate the hemolysis value of the blood pump and the average hemolytic value is 0.005 7, which shows that the modified hemolytic model has good hemolytic capability.

blood pump; multiphase flow; standard hemolytic index; hemolytic model

TH312

A

1672?7207(2018)08?1929?07

10.11817/j.issn.1672?7207.2018.08.013

2017?08?10；

2017?10?21

國家自然科學基金資助項目(51475477，3167099)(Projects(51475477, 3167099) supported by the National Natural Science Foundation of China)

譚建平，博士，教授，從事電液控制及血泵研究；E-mail：jptan@163.com