基于流固耦合的介入機器人診療時血流動力學(xué)分析

朱宗銘 季蘇強 王 浩 唐蒲華 梁 亮

1.長沙學(xué)院機電工程學(xué)院,長沙,4100222.湘潭大學(xué)機械工程與力學(xué)學(xué)院,湘潭,4111053.湖南師范大學(xué)工程與設(shè)計學(xué)院,長沙,410012

0 引言

區(qū)別于傳統(tǒng)的血管介入手術(shù),介入機器人可在血管內(nèi)自由移動,并完成診斷血管、疏通血管和定點投藥等工作[1]。血流動力學(xué)參數(shù)包括血流速度、血壓、血流阻力(血管壁面剪切應(yīng)力)等指標(biāo)。血液在血管內(nèi)流動時,血液流動和血管變形存在流固耦合效應(yīng)。人們常引入計算流體力學(xué)、流固耦合等方法來預(yù)測血流動力學(xué)參數(shù)的變化。

血管內(nèi)無介入裝置時,國內(nèi)外學(xué)者考慮流固耦合效應(yīng),對不同部位血管(胸動脈[2]、頸動脈[3]、椎動脈[4])、不同形狀血管(直管[5]、彎管[6]和分叉管[7])和不同病變血管(含斑塊[8]和囊腫[9])的血流動力學(xué)參數(shù)進行了數(shù)值計算和實驗研究。

血管內(nèi)存在介入裝置時,血管的結(jié)構(gòu)將發(fā)生變化,血液流動也隨之改變。PEWOWARUK等[10]對導(dǎo)管介入前后的豬肺動脈分支狹窄模型中的血液動力學(xué)進行了數(shù)值模擬預(yù)測。JOHARI等[11]基于流固耦合效應(yīng),研究了專用支架植入頸動脈分叉血管后的管內(nèi)血流動力學(xué)。DAI等[12]探討了多個重疊非覆膜支架的孔隙率對主動脈夾層血流動力學(xué)的影響。WILLIAMSON等[13]采用粒子圖像測速(particle image velocimetry,PIV)技術(shù)實驗研究了主動脈髂分叉管中植入支架前后的分叉管近端和遠端的血流動力學(xué)參數(shù)。江帆等[14]采用流固耦合方法,數(shù)值計算了彈性血管和剛性血管中介入機器人在固定位置旋轉(zhuǎn)時的血流速度和血壓等參數(shù)。LIANG等[15]研究了平移速度和轉(zhuǎn)速等運行參數(shù)對剛性血管內(nèi)運行的介入機器人綜合性能的影響。

血流動力學(xué)研究分為血管內(nèi)無介入裝置和有介入裝置兩種。無介入裝置時,考慮血管彈性(流固耦合作用)的研究主要涉及血管的部位、幾何形狀和病變結(jié)構(gòu)。有介入裝置時,研究主要涉及介入裝置的類型(導(dǎo)管、支架、機器人等)、介入裝置的運動狀態(tài)(靜止、旋轉(zhuǎn)、移動等)、血管的力學(xué)性質(zhì)(彈性和剛性)等。但針對彈性血管的研究未考慮介入裝置的診療移動,針對介入裝置移動的研究又假設(shè)在剛性血管內(nèi),即缺乏綜合分析介入裝置運動、血液流動和血管變形的耦合關(guān)系。本文基于血液和血管的雙向流固耦合作用,對脈動血流環(huán)境中,介入機器人在血管內(nèi)不同高度旋進診療時的血流特性開展數(shù)值計算和實驗研究。

1 血管介入機器人的結(jié)構(gòu)

介入機器人系統(tǒng)采用永磁體驅(qū)動,如圖1所示,外部驅(qū)動永磁體為徑向磁化的空心圓環(huán)。介入機器人外表面帶有螺旋肋,內(nèi)部有一個徑向磁化的實心圓柱磁鐵(與外部永磁體的異極相吸)?？刂仆獠坑来朋w的平移和旋轉(zhuǎn)(旋進)運動來驅(qū)動介入機器人做同樣的旋進運動,清除血管內(nèi)的血栓。

圖1 血管介入機器人系統(tǒng)結(jié)構(gòu)示意圖

2 數(shù)值計算模型與方法

2.1 動力學(xué)控制方程

假設(shè)血液流體不受溫度影響且不可壓縮,則流體域控制方程包括連續(xù)性方程和動量守恒方程[16]：

(1)

(2)

固體域控制方程可由牛頓第二定律導(dǎo)出:

(3)

式中,σs為血管壁應(yīng)力張量;Fs為血管壁體積力;ρs為血管壁密度;ds為血管壁區(qū)域的局部加速度。

式(1)、式(2)為黏性流體動力學(xué)控制方程,式(3)為血管壁的固體控制方程,它們組成數(shù)值計算的數(shù)學(xué)模型。

在不考慮溫度和熱量影響的情況下,血液和血管壁耦合交界面處流體和固體滿足

(4)

式中,σ為應(yīng)力張量;n為邊界法向向量;u為位移矢量;下標(biāo)f表示流體,s表示固體。

2.2 系統(tǒng)建模

采用Pro/E軟件建立包括機器人、血管和血液流體的介入機器人系統(tǒng)三維模型,相關(guān)參數(shù)見表1。為模擬機器人的運動,在機器人表面設(shè)計了一層包裹的流體(流體區(qū)域Ⅰ)。流體區(qū)域Ⅰ的厚度為0.5 mm。血管內(nèi)剩余流體為流體區(qū)域Ⅱ。機器人系統(tǒng)的計算幾何模型如圖2所示。

表1 介入機器人系統(tǒng)相關(guān)參數(shù)

圖2 血管機器人系統(tǒng)計算幾何模型

2.3 網(wǎng)格劃分

流固耦合過程中使用分區(qū)耦合算法計算,即流體域和固體域先獨立計算,再進行數(shù)據(jù)交換,因此采用ANSYS-MESH模塊對血管和血液進行網(wǎng)格劃分。血管區(qū)域采用六面體結(jié)構(gòu)網(wǎng)格Hex,默認單元尺寸為4.2426 mm。血液區(qū)域采用以四面體網(wǎng)格Tet為主、以楔形網(wǎng)格為輔的網(wǎng)格,默認單元尺寸為4.1976 mm。

流固耦合計算時,網(wǎng)格尺寸對計算結(jié)果和計算速度有較大影響。網(wǎng)格無關(guān)性分析的參考指標(biāo)設(shè)為介入機器人運行時的血流出口平均速度。根據(jù)上述網(wǎng)格的劃分方法,逐步對模型進行網(wǎng)格加密,并在其他參數(shù)相同的條件下進行數(shù)值模擬計算。采用不同網(wǎng)格數(shù)量的模型進行數(shù)值模擬,表2所示為網(wǎng)格數(shù)量與血流出口平均速度的關(guān)系,網(wǎng)格數(shù)量從387 904增加到632 711時,血流出口速度變化率(偏差)小于1%,滿足網(wǎng)格無關(guān)性要求。綜合考慮,選用第3組網(wǎng)格模型進行后續(xù)的數(shù)值模擬研究。

表2 網(wǎng)格無關(guān)性研究結(jié)果

2.4 雙向流固耦合計算方法

介入機器人在血管中運行時,血管在流動血液的壓力作用下產(chǎn)生變形,血管變形又會對血液流動造成影響,因此需構(gòu)建雙向流固耦合模型。

雙向流固耦合研究中,流體區(qū)域(血液)和固體結(jié)構(gòu)區(qū)域(血管)是并行求解的。每個子步驟中,流體域解和固體域解都需要收斂,才能繼續(xù)下一個時間步長的計算,否則需要重新計算直至結(jié)果收斂。求解過程中,數(shù)據(jù)傳遞順序如下:流體域?qū)毫π畔鬟f給固體域,固體域接收到壓力信息后產(chǎn)生位移,位移引起計算區(qū)域發(fā)生變化,進而發(fā)生網(wǎng)格重構(gòu),之后,固體域再反饋給流體域。上述步驟完成一次即完成一個時間步長的計算,計算時間遞增,直至達到計算周期結(jié)束,否則將固體域的位置與速度信息傳遞給流體域。雙向流固耦合的計算過程如圖3所示。

圖3 雙向流固耦合計算流程圖

2.5 脈動血流入口條件和求解器參數(shù)設(shè)置

心房和心室規(guī)律性的收縮和舒張實現(xiàn)了間歇射血,使得動脈血具有顯著的脈動特性。假設(shè)心臟搏動頻率為每分鐘75次,即心臟搏動周期為0.8 s。根據(jù)文獻[17],采用圖4所示的主動脈入口速度曲線擬合血流入口速度與時間的關(guān)系:

圖4 主動脈血流入口速度曲線

vinlet=

(5)

采用用戶自定義函數(shù)(userdefined function,UDF)編譯脈動血流口速度條件,出口壓力設(shè)為0。將血管的進出口端面定義為固定支撐。因為機器人的旋轉(zhuǎn)運動,血管內(nèi)流體的流動狀態(tài)為湍流,故計算時采用標(biāo)準(zhǔn)k-ε湍流模型。設(shè)置湍流強度I為5%,水力直徑D為18 mm。

根據(jù)實際情況,數(shù)值計算中設(shè)定機器人的平移速度vr為4 mm/s,運動方向與血流方向相反,機器人旋轉(zhuǎn)速度nr為120 r/min。采用標(biāo)準(zhǔn)的SIMPLE算法求解流體壓力和速度耦合方程,壓力、動量、湍流動能和耗散率的差分格式均為二階迎風(fēng)格式。為模擬機器人表面鄰近區(qū)域流體的運動,對機器人外層流體進行滑移網(wǎng)格設(shè)置,給定機器人外層流體轉(zhuǎn)速等于機器人轉(zhuǎn)速。

數(shù)值計算為非穩(wěn)態(tài)計算,假定機器人沿著X軸正方向(圖2)在管道內(nèi)做旋進運動,其運動規(guī)律通過UDF控制。機器人移動和流固耦合作用導(dǎo)致的網(wǎng)格更新通過動網(wǎng)格重構(gòu)實現(xiàn)。動網(wǎng)格區(qū)域創(chuàng)建中,流體與血管的交界面設(shè)定為系統(tǒng)耦合面,動網(wǎng)格重構(gòu)時對鄰近的邊界層進行變形處理。系統(tǒng)耦合模塊中,選擇流體與血管的交界面為數(shù)據(jù)傳輸界面,其中,力傳輸?shù)膩喫神Y因子設(shè)定為0.5。計算時間步長設(shè)置為0.01 s,計算結(jié)束時間為0.8 s。

3 數(shù)值計算結(jié)果與分析

機器人進入血管后,血液流動和血管變形都將受到影響,介入機器人在血管內(nèi)不同高度運行時的血流動力學(xué)參數(shù)也將變化。假設(shè)血管中心軸的Y坐標(biāo)為0,介入機器人在血管內(nèi)5個不同高度(y1～y5)即機器人中心點的Y坐標(biāo)分別為2.50 mm、1.25 mm、0 mm、-1.25 mm、-2.50 mm。

3.1 血流速度

圖5所示為血流速度最大正值(t=0.25 s,收縮期內(nèi))時,介入機器人在血管內(nèi)不同位置和無機器人介入時的血流速度流線,其中,血液從左側(cè)入口流向右側(cè)出口,機器人從右向左運動。由圖5可知,t=0.25 s時,介入機器人在靠近血管壁時,血流速度略微增大,在血管中心的血流速度略微偏小,但總體差別較小,可見介入機器人高度對血流速度影響較小。血管內(nèi)無介入機器人時,血流狀態(tài)為層流,血管中心軸處的血流速度最大,為0.28 m/s,明顯小于有介入機器人時的最大血流速度0.48 m/s。機器人的介入使得血管內(nèi)徑變小,血液流過機器人時,機器人周圍流體域受到壓縮,該區(qū)域血液流動速度增加,形成了高速流域。總體來說,不論介入機器人在血管哪個高度運行,血流速度均在人體的安全范圍內(nèi),不會對血管造成損傷。

(a)機器人高度為2.50 mm

3.2 血壓

圖6為血流速度最大正值(t=0.25 s,收縮期內(nèi))和最大負值(t=0.45 s,舒張期內(nèi))時,介入機器人在血管不同高度和無機器人介入時的血液對血管壁壓力(血壓)分布圖,可以看出,t=0.25 s時,機器人在不同高度的血壓變化趨勢大致相同:血管入口端的血壓最大,并沿血流方向逐漸減小;血壓最大值為120 Pa。無機器人介入時,血管內(nèi)血壓均勻變化,從入口向出口逐漸變小,血壓最大值為72 Pa。因此,機器人的介入使血管內(nèi)的血壓變化不均勻,機器人周圍的血壓分布更復(fù)雜。收縮期內(nèi),機器人附近和血管末端會形成局部負壓區(qū)。無機器人介入時,血流負壓區(qū)僅在血管末端產(chǎn)生。在心臟收縮期,介入機器人在不同高度位置運行時,血壓變化很小。

(a)機器人高度為2.50 mm(t=0.25 s)

t=0.45 s時,血液流動方向發(fā)生變化,部分血液回流,血壓為負值(負壓),血壓絕對值從入口端向出口端遞減。無機器人介入時,血壓變化相對均勻,并且小于機器人介入的情況。介入機器人在不同高度的血壓的最大負值和最小負值均相差較小。

因此,在心臟收縮和舒張期內(nèi),介入機器人在血管內(nèi)運行高度的變化對血壓影響均較小,且壓力(-92 ～120 Pa)在人體安全范圍內(nèi),不會對血管壁產(chǎn)生損傷,這說明介入機器人在血管內(nèi)診療運行時是安全的。

3.3 血管壁面剪切應(yīng)力

正常情況下,血液流動的剪切力對血管壁有益,能維持血管內(nèi)皮細胞的功能。剪切應(yīng)力過大或過小時,血管內(nèi)皮細胞功能發(fā)生異常,影響血管壁的健康。因此,維持適當(dāng)?shù)谋诿婕羟袘?yīng)力對維護心血管健康至關(guān)重要。

壁面剪切應(yīng)力與速度梯度有關(guān),但較難通過臨床直接測量獲得。圖7所示為一個脈動血流周期內(nèi)的平均血管壁剪切應(yīng)力-時間曲線。圖7中,縱坐標(biāo)為正值表示壁面剪切應(yīng)力方向向右,負值表示壁面剪切應(yīng)力方向向左。無機器人介入時,血管壁所受的最大剪切應(yīng)力為0.076 Pa。血管壁面剪切應(yīng)力長時間小于0.4 Pa會導(dǎo)致養(yǎng)分供給不足,使該區(qū)域易產(chǎn)生粥樣硬化斑塊。機器人介入后,血管壁面剪切應(yīng)力明顯增大,無論機器人在哪個高度運行,血管壁面剪切應(yīng)力均在正常范圍內(nèi)波動,更有利于維持血管健康。

圖7 平均血管壁剪切應(yīng)力-時間曲線

由圖7可以看出,介入機器人在不同高度的血管壁面剪切應(yīng)力的變化趨勢大致相同,并且與主動脈血流入口速度的變化規(guī)律基本一致。機器人越靠近血管壁,血管壁面剪切應(yīng)力越大。因此,在機器人介入血管時,可調(diào)整機器人運行位置來改善壁面剪切應(yīng)力以維持血管健康。

3.4 血管性質(zhì)

上述研究考慮真實血管的彈性,采用流固耦合的分析方法,但血管性質(zhì)對計算結(jié)果影響的大小如何,有待進一步分析。

表3所示為脈動血流收縮期內(nèi)最大血流入口速度(t=0.25 s)時的血流動力學(xué)參數(shù),可以看出,彈性血管模型下的最大血流速度、最大血壓和血管壁面剪切應(yīng)力均略小于剛性血管模型的結(jié)果。這是因為血管壁的彈性對血液流動起到一定的阻礙作用,使血流速度略微降低,導(dǎo)致血壓和壁面剪切應(yīng)力都減小。綜合來看,考慮血管彈性的流固耦合計算結(jié)果更接近真實血管。

表3 彈性血管模型和剛性血管模型的血流動力學(xué)參數(shù)(t=0.25 s)

4 介入機器人管內(nèi)脈動流場實驗臺

基于磁控機器人管內(nèi)恒定流場測量系統(tǒng)[18],搭建了介入機器人管內(nèi)脈動流場測量系統(tǒng),如圖8所示。介入機器人驅(qū)動部分主要由六自由度機器手、旋轉(zhuǎn)電機、工作臺、外部驅(qū)動永磁體、試驗管道、螺旋介入機器人等組成。介入機器人系統(tǒng)相關(guān)參數(shù)與表1相同。實驗中,外部永磁體平移速度和轉(zhuǎn)速與數(shù)值計算中的數(shù)值相同。

(a)全景圖

流場測量部分主要由粒子圖像測速系統(tǒng)、蠕動泵、壓力計、流量計、玻璃水槽、水槽支架和計算機等組成。脈動流由蠕動泵產(chǎn)生,通過流量計和壓力計對流體數(shù)據(jù)進行監(jiān)控。

結(jié)合流量計,對蠕動泵輸出流體速度的幅值和脈動頻率進行調(diào)節(jié),使一個周期的流體速度與圖4中的血流入口速度曲線紅色部分相近,并將輸出流體速度繪制成曲線,如圖9所示。

圖9 蠕動泵輸出流體速度曲線

5 實驗結(jié)果與分析

圖10所示為蠕動泵輸出流體速度達到最大值時,機器人在高度2.50 mm處向左運動時周圍的流體速度,可以看出,機器人向左運動時,機器人下方存在一個流體高速區(qū),機器人前后方為流體中速區(qū),機器人尾部出現(xiàn)流體低速區(qū)。機器人介入后,機器人下方管道變窄,流體壓力增大,流速隨之增加,形成流體高速區(qū)?？傮w來看,機器人周圍流體速度的實驗值與計算值大體相似。實驗時,外磁體的運動驅(qū)動機器人運動,機器人運動的滯后性導(dǎo)致機器人所受磁吸力不沿水平方向,機器人前進時與管壁反復(fù)碰撞,因此機器人運動軌跡呈波浪形。數(shù)值計算時為了簡化計算,假設(shè)機器人沿水平方向運動,因此實驗測量的管內(nèi)流體速度沒有數(shù)值計算結(jié)果均勻。

(a)實驗測量

為定量比較實驗測量和數(shù)值計算的結(jié)果,建立圖11所示的參考坐標(biāo)系,其中,管道中心軸為X軸,通過機器人中心的垂直線為Y軸,vf為流體速度,vr為機器人速度。

圖11 流體速度定量比較的參考坐標(biāo)系

圖12所示為蠕動泵輸出流體速度達到最大值時,介入機器人旋進和無介入機器人時的Y軸上管內(nèi)流體速度的實驗值(多次實驗測量結(jié)果的平均值)和計算值,可以看出,無介入機器人時,對于脈動流入口,從管道底部往上,管內(nèi)流體速度先增大,再相對平穩(wěn),最后減小,數(shù)值計算結(jié)果變化趨勢與實驗測量結(jié)果基本一致。管內(nèi)中間段流體區(qū)域的數(shù)值計算結(jié)果與實驗測量結(jié)果的最大誤差約為13.1%。有機器人介入時,從管道底部往上,管內(nèi)流體速度先快速增大,然后快速減小,機器人位置處的流體速度接近0,數(shù)值計算結(jié)果與實驗測量結(jié)果的變化趨勢基本一致。機器人和管道之間流體區(qū)域的數(shù)值計算結(jié)果與實驗測量的結(jié)果最大誤差約為19.7%。造成數(shù)值計算結(jié)果與實驗測量結(jié)果誤差的主要原因是,實驗中機器人沿X軸運行時存在一定的上下波動。對比來看,機器人的介入使管內(nèi)機器人下方流體速度增大。

(a)無介入機器人

6 結(jié)論

使用雙向流固耦合方法獲得了介入機器人在彈性血管內(nèi)運行時的管內(nèi)流場特性。對比分析了無介入機器人和介入機器人在血管中不同高度時的血流流線和血壓的變化差異,討論了血管壁面剪切應(yīng)力與血流入口速度的關(guān)系。研究發(fā)現(xiàn),機器人的介入對血流速度、血壓、血管壁面剪切應(yīng)力等血流動力學(xué)參數(shù)有一定的增大作用,但數(shù)值均在安全范圍以內(nèi),不會對血管造成損傷;機器人運行高度的變化對血管壁面剪切應(yīng)力略有影響;彈性血管模型下的血流速度、血壓、血管壁面剪切應(yīng)力等血流動力學(xué)參數(shù)均低于剛性血管模型下結(jié)果。

搭建了介入機器人管內(nèi)脈動流場測量系統(tǒng),在彈性管道中實驗測量了介入機器人旋進和無介入機器人的管內(nèi)機流體速度。結(jié)果表明:實驗測得的機器人周圍流體速度與數(shù)值計算結(jié)果基本一致,兩者誤差在20%以內(nèi)。這對介入機器人運動時管內(nèi)復(fù)雜流場的計算是可以接受的。