動靜區(qū)域耦合方法對血泵模擬結(jié)果的影響研究

林子祺， 于 勇， 張 宇

(1.北京理工大學(xué) 宇航學(xué)院流體力學(xué)實驗室， 北京 100081；2.清華大學(xué) 醫(yī)學(xué)院航天泰心科技有限公司人工心臟聯(lián)合研究中心， 北京 100091)

引言

心室輔助裝置(Ventricular Assist Device,VAD)以醫(yī)工結(jié)合的方式幫助病人維持血液循環(huán)。對于VAD此類復(fù)雜泵體，計算流體力學(xué)(Computational Fluid Dynamics,CFD)能減少實驗開銷，實現(xiàn)對設(shè)備的評估和優(yōu)化；但同時需要考慮CFD準(zhǔn)確度的問題[1-3]。為方便進行CFD準(zhǔn)確度的研究，美國食品與藥品監(jiān)管局(Food and Drug Administration, FDA)提出了“關(guān)鍵路徑計劃”(Critical Path Initiative, CPI)，發(fā)布了一個基準(zhǔn)模型及圍繞該模型獲得的多個工況下的實驗結(jié)果，以供與CFD仿真結(jié)果進行比對。2012年SANDY等[4]給出了該模型的詳細幾何參數(shù)。相較現(xiàn)有產(chǎn)品設(shè)計(如Impella、新一代HeartMate等)，標(biāo)準(zhǔn)血泵對紅細胞的破壞性更強。但標(biāo)準(zhǔn)模型有詳細的流動和溶血指數(shù)實驗數(shù)據(jù)，是評價CFD仿真VAD準(zhǔn)確性的合適選擇[5]。

2017年，MALINAUSKAS等[6]總結(jié)了來自全球24個CFD仿真小組針對標(biāo)準(zhǔn)模型使用不同網(wǎng)格、湍流模型計算得出的結(jié)果，并與實驗數(shù)據(jù)進行對比。所有小組對葉輪附近的流速仿真較為準(zhǔn)確，但對泵出口擴散段的流速仿真結(jié)果較差。對于溶血估計的計算結(jié)果所有小組之間的差別很大，同時與實驗差別也較大。

對常規(guī)工業(yè)泵的CFD模擬，所選湍流模型和定常/非定常方法的不同會影響CFD的結(jié)果[7-8]。在對標(biāo)準(zhǔn)血泵的模擬中，MALINAUSKAS[6]總結(jié)的24個小組中有14組明確使用定常仿真方法。近幾年針對VAD的CFD仿真也大多使用定常方法，但是在出口擴散段流速與溶血率絕對值的估計準(zhǔn)確度有待提升。使用非定常仿真一般能獲得比定常仿真更準(zhǔn)確的結(jié)果，但是計算量大，是否值得在工程上推廣應(yīng)用是一個值得討論的問題。

本研究針對FDA標(biāo)準(zhǔn)血泵模型，對比了定常SRF和非定常滑移網(wǎng)格的仿真方法，評估了兩種方法對速度、溶血率仿真準(zhǔn)確性的影響，來討論非定常計算的必要性。

1 材料與方法

1.1 FDA標(biāo)模設(shè)計與工況

研究所用三維模型引自網(wǎng)站https://nciphub.org/wiki/FDA_CFD。轉(zhuǎn)子直徑52 mm，泵體內(nèi)徑60 mm，轉(zhuǎn)子背面與泵體、轉(zhuǎn)子葉刀頂部與泵體的間隔皆為1 mm。標(biāo)準(zhǔn)泵模型及對應(yīng)的實驗裝置如圖1所示。具體幾何尺寸可下載網(wǎng)站中的Blood_Pump.zip文件。

圖1 FDA發(fā)布的CPI實驗裝置

FDA為此標(biāo)準(zhǔn)血泵做了6個工況的實驗[7]。本研究選用了其中實驗數(shù)據(jù)完整的4個工況，工況參數(shù)見表1。研究所提及的“實驗數(shù)據(jù)”皆指由FDA公布的流場或溶血數(shù)據(jù)。

表1 FDA血泵模擬工況參數(shù)

實驗使用了粒子圖像測速法(Particle Image Velocimetry, PIV)技術(shù)測量血泵內(nèi)特定剖面的流場。實驗用血為豬血液，密度1056 kg/m3，動力黏度0.0034 kg/(m·s)。實驗時以雷諾數(shù)Re和流動系數(shù)φ為基準(zhǔn)對實驗時使用的流體流量與轉(zhuǎn)速進行適當(dāng)調(diào)整[9]。所用雷諾數(shù)的定義是

(1)

式中，ρ—— 血液密度

ω—— 轉(zhuǎn)子轉(zhuǎn)速

D—— 轉(zhuǎn)子直徑

μ—— 血液動力黏度

流動系數(shù)的定義是：

(2)

式中，Q—— 入口處血液質(zhì)量流量

實驗中定義坐標(biāo)原點為轉(zhuǎn)子無葉片一面的中心、z軸正方向與流體流入方向相反、x軸正方向與出口流體流動方向相同。同時定義了若干條采用線段(在使用PIV拍攝的剖面內(nèi))，采集線段上的流速，采樣線段的命名沿用文獻[7]，CS2，CS3，CS5分別表示不同PIV的拍攝平面，在本研究中B表示沿象限平分線的采樣線段，L1～L4表示擴散段不同的采樣線段。坐標(biāo)軸、拍攝剖面示意及采樣線段所在位置如圖2所示。采樣線段端點坐標(biāo)見表2。

表2 仿真流速采集線段定義

圖2 采樣線段所在位置

1.2 流動基本方程

血液視為單相牛頓流體。數(shù)值求解RANS和URANS方程組。其中動量方程為：

(3)

式中，U—— 速度張量

τ—— 瞬時剪切應(yīng)力

f—— 動力黏度附加力

湍流模型采用RNG k-epsilon 模型：

(4)

其中,C1=1.42,C2=1.68。

湍流黏性系數(shù)

(5)

其中，Cμ=0.0845。

在計算中，Y+為10左右，壁面采用增強壁面函數(shù)(Enhanced Wall Treatment)

(6)

仿真中模型入口直徑為12 mm，入口邊界條件見表3。出口使用自由流邊界(outflow boundary)。

表3 模擬速度入口邊界參數(shù)

1.3 流場仿真可信度分析

本研究使用了3個參數(shù)衡量仿真所得到流場與實驗結(jié)果的偏差：確定系數(shù)D、極值誤差α與極值點誤差α0，其中確定系數(shù)越接近1表示仿真與實驗的整體偏差越小，極值與極值點誤差用以定量描述仿真與實驗曲線在極點處的偏離程度。

(7)

(8)

(9)

式中，vi—— PIV方法所測一個空間點在xOy平面內(nèi)的速度

vmax—— 使用PIV方法所測的特定線段上在xOy平面內(nèi)的速度最大值

x0—— 取得vmax時特定線段上的位置參數(shù)

xmax,xmin—— 分別為采樣線段坐標(biāo)的最大值與最小值

1.4 溶血預(yù)測模型與方法

CFD用于評估溶血的絕對水平非常具有挑戰(zhàn)性[10]。GIERSIPEN等[11]提出了一種基于應(yīng)力的冪律模型，該方法將血漿游離血紅蛋白的產(chǎn)生與剪切應(yīng)力大小和暴露時間聯(lián)系起來：

(10)

H是溶血率，即血漿中的自由血紅蛋白和總的血紅蛋白的百分比。Hct是血細胞比容，fHb是自由血漿中的血紅蛋白(mg/dl plasma)，Hb是總的血紅蛋白濃度(mg/dl blood)，C，a和b是經(jīng)驗系數(shù)。C的單位是Pa-1s-1，a和b是無量綱的，系數(shù)通常由在均勻剪切條件下在Couette型剪切裝置中的溶血測量確定[12]。有許多實驗團隊以此方法針對不同的動物血液整定配套的實驗參數(shù)C，a，b[13]。對于豬血，使用SONG等[14]整定的參數(shù)C=1.8×10-8，a=0.765，b=1.991。

借助離散相模型(Discrete Phase Model，DPM)對溶血率進行估計。對于計算已收斂的流場，在入口出撒入無質(zhì)量的顆粒并進行追蹤，統(tǒng)計其所經(jīng)歷的流場單元的剪切應(yīng)力，給出溶血的預(yù)測值，即統(tǒng)計第p個粒子在第i個時間段內(nèi)的溶血率：

(11)

隨后再計算單個粒子在第i個時間段之前的累計溶血率：

Dp,i=Dp,i-1+(1-Dp,i-1)dp,i

(12)

最后，根據(jù)所有粒子的溶血率進行統(tǒng)計平均

(13)

其中，N為布置的DPM粒子個數(shù)，P為粒子編號。參考Richard的處理方式[7]，相關(guān)溶血系數(shù)將獲得的絕對溶血率與工況C5的絕對溶血率相除，以呈現(xiàn)同一方法對不同工況溶血率變化趨勢的估計能力。

1.5 網(wǎng)格及無關(guān)性驗證

定常與非定常仿真使用相同的網(wǎng)格進行。網(wǎng)格采用組合的方式，由轉(zhuǎn)子周圍流體域及殼體壁面附近流體域組成，不同流體域之間使用interface進行流場信息交換，如圖3所示。圖3中外部流場網(wǎng)格以半透明灰色表示，近轉(zhuǎn)子流場網(wǎng)格以顯示面網(wǎng)格來表示，紫色部分為近轉(zhuǎn)子網(wǎng)格的內(nèi)部固壁邊界(轉(zhuǎn)子)。

圖3 網(wǎng)格組合方式

對于定常仿真，對轉(zhuǎn)子附近的流體域定義旋轉(zhuǎn)坐標(biāo)系；對非定常仿真，轉(zhuǎn)子附近流體域有旋轉(zhuǎn)運動，其他流體域靜止，動、靜區(qū)域通過interface進行耦合，在保證庫朗數(shù)小于1的前提下選擇合適的時間步長：為保證進行C6工況(流速最大)仿真時庫朗數(shù)小于1，選取的時間步長為5e-5 s，在進行C1/C2/C5的仿真時沿用此時間步長。

進行網(wǎng)格無關(guān)性驗證的不同網(wǎng)格屬性及進行部分工況的定常計算結(jié)果見表4。

使用表4中的網(wǎng)格對C5工況也進行了仿真，獲取表3中各個線段的流況數(shù)據(jù)并與實驗數(shù)據(jù)進行對比，得到確定系數(shù)相差小于0.001。三種網(wǎng)格對工況C1/C2進出口壓頭的計算結(jié)果偏差在2%以內(nèi)。特別地，根據(jù)FDA實驗所獲得的C1工況壓頭為172±12 mmHg，C2工況壓頭為353±18 mmHg，網(wǎng)格M2與M3獲得的壓頭均在此范圍內(nèi)且十分接近。綜合考慮計算資源的消耗與計算資源的提升，使用M2網(wǎng)格進行后續(xù)的計算。

表4 多種參與測試的網(wǎng)格屬性

2 結(jié)果

2.1 流場

先后使用SRF與滑移網(wǎng)格方法對表2中4個工況進行仿真，在表3所示部分采樣位置下分析流場速度，按照1.3可信度分析方法獲得的結(jié)果見表5。

表5 仿真的定量比較記錄

圖4表示了針對采樣線段CS2B，四種工況下分別使用SRF與滑移網(wǎng)格方法獲得的仿真數(shù)據(jù)與實驗數(shù)據(jù)比較結(jié)果。可以看出，滑移網(wǎng)格方法對實驗結(jié)果的擬合比SRF方法更加準(zhǔn)確。

圖4 采樣線段CS2B四個工況下定常SRF與非定常滑移網(wǎng)格仿真與CPI實驗流速對比

圖5表示了工況C2下，對出口擴散段的4條采樣線段，使用SRF與滑移網(wǎng)格方法進行的仿真數(shù)據(jù)與實驗數(shù)據(jù)比較結(jié)果。相比于泵內(nèi)的流動，出口擴散段的仿真距離實驗結(jié)果仍有一定差距。但是對工況C2的CS5L3/CS5L4兩條采樣線段而言，滑移網(wǎng)格的擬合效果有明顯改善，但在別的工況、別的出口擴散段采樣線段中，滑移網(wǎng)格方法并沒有獲得更接近實驗數(shù)據(jù)的結(jié)果。

圖5 工況C2出口擴散段采樣線段CPI實驗和仿真結(jié)果對比

2.2 溶血系數(shù)

圖6所示為仿真獲得的相關(guān)溶血系數(shù)σ與實驗數(shù)據(jù)C的對比。前人仿真數(shù)據(jù)引自MALINAUSKAS的整理工作[6]。可以看出在使用基于應(yīng)力的冪律模型預(yù)測溶血時，無論使用定常SRF還是非定常滑移網(wǎng)格方法計算都不能獲得較準(zhǔn)確的效果。

圖6 相對溶血系數(shù)對比圖

3 討論

研究使用的網(wǎng)格及計算方法進行過網(wǎng)格無關(guān)性驗證，在泵內(nèi)流場(采樣線段CS2B/CS3B)的仿真上對比實驗數(shù)據(jù)的確定系數(shù)在0.9以上，C6工況有略低的確定系數(shù)，而滑移網(wǎng)格方法能獲得比SRF方法更高確定系數(shù)的結(jié)果。

在出口擴散段流場(采樣線段CS5L1～4)信息的仿真上，離出口越遠，仿真數(shù)據(jù)與實驗數(shù)據(jù)的差距越大，在C6工況差距最大，且在此處非定常滑移網(wǎng)格方法也不能保證比定常SRF方法更接近實驗數(shù)據(jù)。對于溶血的預(yù)測亦然。針對溶血預(yù)測的準(zhǔn)確性，WU等人認(rèn)為基于應(yīng)力的冪律模型對溶血的預(yù)測本身并不準(zhǔn)確[17]，如果采用基于WU等人提出的能量耗散冪律模型會準(zhǔn)確很多[18]，這值得進一步關(guān)注。從目前的結(jié)果看來，非定常滑移網(wǎng)格方法在溶血方面沒有顯著的優(yōu)勢可能與此有關(guān)。

在數(shù)值仿真中，也在同樣的部位出現(xiàn)了與實驗值差異較大與溶血預(yù)測不準(zhǔn)確的問題[10]。且對于偏離設(shè)計的C6工況，各項參數(shù)也會出現(xiàn)數(shù)據(jù)離散[5]。如何把出口擴散段的流場信息仿真準(zhǔn)確及獲得準(zhǔn)確的溶血預(yù)測結(jié)果，是接下來的工作難點。對于CFD仿真的準(zhǔn)確度問題，使用定常方法還是非定常方法只是一個討論的側(cè)面。HUO等人也指出湍流模型、網(wǎng)格類型與密度、時間步長等因素對FDA血泵流場的預(yù)測結(jié)果都有很大影響[16]。在改變這些因素之后，在難以準(zhǔn)確預(yù)測的區(qū)域非定常方法是否會有比定常方法更顯著的優(yōu)勢，有待進一步研究。

SRF方法不改變網(wǎng)格結(jié)構(gòu)，通過施加旋轉(zhuǎn)坐標(biāo)系來進行定常的迭代計算。此方法在模擬軸流泵(如HearMate II)時[15]，由于出口法向量與轉(zhuǎn)子轉(zhuǎn)軸同向、流體總是向著轉(zhuǎn)軸運動，可以認(rèn)為是與實際情況吻合的。但對于本研究中的FDA標(biāo)準(zhǔn)血泵，出口法向量垂直于轉(zhuǎn)子轉(zhuǎn)軸，SRF方法模擬的準(zhǔn)確性尚待討論。而滑移網(wǎng)格方法通過改變轉(zhuǎn)子角位置進行非定常仿真，可以保證與實驗實際情況吻合，但會消耗更多的計算資源。本研究亦表明，滑移網(wǎng)格方法確實在一定程度上獲得了比SRF方法更準(zhǔn)確的仿真結(jié)果。

4 結(jié)論

非定常滑移網(wǎng)格方法同定常SRF方法相比，結(jié)果更貼合實驗結(jié)果。對于泵內(nèi)流速場的求解，非定常方法能有更高的確定系數(shù)，尤其對于峰值速度、對取得峰值所在極值點的仿真比定常方法明顯更加準(zhǔn)確；對于出口擴散段的求解，兩種方法對于越靠外側(cè)的特征線段求解與實驗結(jié)果差異越大，其中非定常方法在外側(cè)的求解能優(yōu)于定常結(jié)果，但仍然不能從本質(zhì)上改善計算準(zhǔn)確性。

綜上所述，非定常求解結(jié)果在一定程度上會比定常求解更準(zhǔn)確，但考慮到工程應(yīng)用時的計算效率問題，尤其是在溶血估計上的不良表現(xiàn)，可以考慮選擇定常方法在損失少量精度的情況下較快地獲得計算結(jié)果。