基于VLBM的髂外動脈血流動力學分析

戎晨彬,陳 柔,嚴微微,俞慧丹,許友生

(1.浙江科技學院 機械與能源工程學院,杭州 310023;2.中國計量大學 計量測量工程學院,杭州 310018;3.印第安納大學-普渡大學印第安納波利斯分校 機械工程系,印第安納波利斯 46202)

下肢動脈粥樣硬化狹窄是一種常見的動脈狹窄疾病,其主要特征是由下肢動脈內(nèi)壁上形成沉積[1](斑塊)而導致的動脈彈性下降,管徑狹窄。臨床表現(xiàn)早期有間歇性跛行、下肢發(fā)涼麻木和靜息性下肢疼痛等,晚期則伴隨著動脈狹窄可能會出現(xiàn)發(fā)紺、潰瘍和壞疽等[2-3]。該病的發(fā)病率逐年增高,且隨年齡的增長而增長,據(jù)調(diào)查[4-5],45歲以上人群的發(fā)病率約為5%,而70歲以上人群的發(fā)病率則大于15%。

髂外動脈是髂總動脈的終支之一,向下肢供應血液。準確的動脈血壓測量對判斷動脈粥樣硬化狹窄程度十分重要,如脈搏波傳導速度,該指標是測量動脈硬化的可靠方法,但測量數(shù)值與血壓(舒張壓和收縮壓)的相關性較強[6-7]。有創(chuàng)動脈血壓(invasive blood pressure,IBP)是監(jiān)測探頭經(jīng)動脈穿刺置管后直接測量動脈血壓,是測量血壓的金標準,測量值準確可靠;但并發(fā)癥較多,容易造成血腫、血栓、感染、血管損傷等,并且對穿刺技術要求極高[8]。計算機斷層掃描血管造影(computed tomography angiography,CTA)是一種臨床常用的無創(chuàng)動脈診斷方法,但它無法提供相關的流體力學信息。近年來,隨計算機硬件水平的快速發(fā)展,CTA與計算流體力學(computational fluid dynamics,CFD)相結合來診斷動脈疾病的技術也愈發(fā)成熟,尤其在冠狀動脈領域[9-10]已取得一定的成效,而在髂動脈領域鮮有相關研究。作為基于介觀模擬尺度的CFD方法,格子玻爾茲曼方法(lattice Boltzmann method,LBM)具有一些獨特的優(yōu)勢,例如易處理復雜的幾何形狀和邊界[11-12],非常適合計算多相流等復雜流體[13-14],以及易實現(xiàn)大規(guī)模的圖形處理器(graphics processing unit,GPU)并行計算[15]。體積格子玻爾茲曼方法(volumetric lattice Boltzmann method,VLBM)是Yu等[16]基于LBM提出的,相比流體粒子位于格點的傳統(tǒng)LBM,VLBM假設粒子均勻分布在格子單元內(nèi),隨碰撞將遷移至相鄰格子單元,并將格子單元按流體、固體和邊界進行分類。該方法可以更好地處理任意彎曲的復雜邊界,并提高了GPU并行計算能力。

在上述研究的基礎上,本研究提出一種基于GPU并行架構開發(fā)的VLBM,可在短時間內(nèi)獲得血流動力學分析結果,從而輔助醫(yī)生臨床診斷。為驗證此方法的可行性,我們對狹窄髂外動脈的血流動力學進行了研究,首先介紹VLBM及其在模擬中的應用;然后對比模擬結果和IBP結果,驗證VLBM在狹窄髂外動脈血流動力學分析中的可行性,并通過將支架植入狹窄處進行修復,獲得恢復動脈的血流動力學信息;最后總結分析狹窄對血流的影響。本研究中狹窄髂外動脈的CTA影像、IBP測量結果等數(shù)據(jù)均由杭州市第一人民醫(yī)院于2017年提供。

1 原理與方法

1.1 體積格子玻爾茲曼方法

體積格子玻爾茲曼方法(VLBM)假設粒子均勻分布在格子中,并將格子單元分成流體、固體和邊界3種類型,如圖1所示。引入?yún)?shù)P(x,t)=Vs(x,t)/V定義每個格子單元中固體體積的百分比,其中,P(x,t)=0的格子單元為流體單元,P(x,t)=1的為固態(tài)單元,以及0

圖1 VLBM中格子單元的分類

我們以Qian等[17]于1992年提出的DdQm(d維空間,m個離散速度)系列LBM基礎模型中的D2Q9模型為例,來討論傳統(tǒng)LBM和VLBM的差異。圖2(a)為傳統(tǒng)LBM模型,粒子在格點處碰撞后,沿其速度方向往相鄰格點遷移,fi(x,t)為粒子密度分布函數(shù)。而在VLBM中,假設粒子均勻分布在格子單元內(nèi),隨碰撞遷移至相鄰格子單元中,如圖2(b)所示。與傳統(tǒng)LBM的fi(x,t)類似,VLBM引入ni(x,t)表示離散速度為ei的粒子在時間t占據(jù)格子單元的分布函數(shù),i為離散速度數(shù)量,fi(x,t)與ni(x,t)的關系為

圖2 傳統(tǒng)LBM和VLBM粒子分布示意圖

ni(x,t)=fi(x,t)[1-P(x,t)]V。

(1)

式(1)中:V取單位體積(V=1)。與LBM相似,VLBM的碰撞方程如下:

(2)

式(2)中:n′i(x,t)為碰撞后的偏分布函數(shù),其中i=0,1,2,…,b;τ為無量綱弛豫時間;ni,eq(x,t)為平衡態(tài)分布函數(shù),

(3)

式(3)中:ωi為權重系數(shù);cs為格子聲速;U為宏觀速度;N(x,t)=∑ni(x,t)為當前格子單元密度。本研究選用i=0,1,2,…,18的D3Q19格子模型進行三維模擬,離散速度ei在i=0時為(0,0,0),i=1～6時為(±1,0,0)、(0,±1,0)、(0,0,±1),以及i=7～18時為(±1,±1,0)、(0,±1,±1)、(±1,0,±1)。相對應的權重系數(shù)在i=0時ωi=1/3,i=1～6時ωi=1/18,以及i=7～18時ωi=1/36。為將非滑移邊界條件整合到格子模型中,流動部分分為兩部分:1)從迎風方向的鄰近格子單元中遷移出的粒子數(shù)[1-P(x,t)]n′i(x-eiδt,t);2)從逆風方向的格子單元反彈的粒子數(shù)P(x-eiδt,t)n′i*(x,t)。此處δt為時間步長,n′i*(x,t)為對應n′i(x,t)從非滑移壁面反彈后的粒子分布函數(shù)。因此,遷移方程為

n″i(x,t+δt)=[1-P(x,t)]n′i(x-eiδt,t)+P(x-eiδt,t)n′i*(x,t)。

(4)

宏觀參數(shù)密度和速度分別從粒子分布函數(shù)的0、1階中獲得:

ρ(x,t)=∑ni(x,t)/[1-P(x,t)];

(5)

u(x,t)=∑eini(x,t)/∑ni(x,t)。

(6)

壓強可從密度中獲得:

(7)

式(7)中:p0為參考壓強;ρ0為參考密度。

1.2 計算模型設置

通過對CTA影像進行處理,提取狹窄髂外動脈形態(tài)結構模型。CTA是一種醫(yī)學檢查方法,它將CT掃描與注射造影劑相結合,產(chǎn)生人體特定部位的多層等距二維數(shù)字圖像,圖像由大量不同灰度值的像素按矩陣排列方式構成。3D Slicer軟件可以設定灰度閾值從而分割出二維圖像中髂動脈部分,提取該部分的輪廓線,并通過幾何單元拼接來擬合動脈壁,最終達到重構3D動脈模型的效果。重構支架植入前的3D髂外動脈狹窄模型如圖3(a)所示,該髂外動脈為重度狹窄,狹窄率約為70%。狹窄率表示動脈的狹窄程度,常用的計算公式為:狹窄率=(狹窄遠端正常直徑-狹窄段的最窄直徑)/狹窄遠端正常直徑。支架植入后如圖3(b)所示,髂動脈模型通過3D Slicer軟件拉伸狹窄處壁面后導出,將狹窄處恢復成健康動脈。

圖3 3D髂外動脈結構模型

為簡便起見,假設血液為不可壓縮牛頓流體,密度為1 060 kg/m3,運動黏度為3.3×10-6m2/s,動脈壁為剛性壁,且壁面采用無滑移條件。將脈動流作為髂總動脈剖面的速度入口邊界條件,在髂內(nèi)外動脈末端設置壓力出口邊界條件,具體邊界條件方案如下:

1)速度入口邊界條件。多普勒超聲(Doppler ultrasound,DUS)可測量一個心動周期內(nèi)的血流速度曲線。在模擬中,需將速度曲線插值成1 000個時間點,如圖4(a)中的藍點所示。本研究選取了6個具有代表性的紅點作為采樣時間點,以表示血流在各個階段(舒張期和收縮期)流動時的不同現(xiàn)象。在指定時刻對髂總動脈剖面設定速度邊界時,將DUS同一時刻測得的速度作為髂總動脈中心處的最大速度,然后用拋物線分布到壁面上,以此貼近實際情況,如圖4(b)所示。值得注意的是,速度值是脈動的,每個時間點的入口邊界應按上述方法設置。為將速度邊界條件引入VLBM,我們采用Guo等[18]于2002年提出的非平衡態(tài)外推格式,該方法具有運算簡單、數(shù)值穩(wěn)定等優(yōu)點。

圖4 髂總動脈入口血流速度

2)壓力出口邊界條件。壓力出口使用與電路相似的三元windkessel模型[19],圖5為該模型應用于動脈中時的示意圖。血流速度(Q)由近端血管阻力(r)、遠端血管順應性(C)和每個出口下游遠端血管阻力(R)的綜合效應控制,因此,三元windkessel模型也被稱為RCR模型。基于血流和電路之間的相似性,流速Q(mào)和平均壓力p之間的關系滿足以下常微分方程[20]:

圖5 動脈中三元windkessel模型示意圖

(8)

求解式(8)可得到壓力的解析解

(9)

式(9)中:p0為出口的初始壓力。R、C和r這3個參數(shù)的關系和取值范圍可通過參考文獻[21]來確定。同樣,采用非平衡態(tài)外推格式將壓力邊界條件引入VLBM。

1.3 算法流程設計

使用C++編程語言開發(fā)基于GPU并行架構的VLBM算法,算法的設計流程具體可分為以下幾個步驟:

1)讀取髂外動脈模型,根據(jù)模型尺寸確定計算區(qū)域;

2)設置格子單元長度l,對計算區(qū)域進行格子劃分,計算每個格子單元中固體體積的百分比P(x,t),并設置邊界條件及給定所有格點初始宏觀參量;

3)使用引入VLBM的D3Q19模型的控制方程,通過給定的初始宏觀參量計算出所有格點各個速度方向的平衡態(tài)分布函數(shù),作為模擬的初場;

4)通過粒子碰撞遷移規(guī)則求解控制方程;

5)使用非平衡態(tài)外推格式對相應邊界處格點進行邊界處理;

6)計算所有格點的宏觀參量,如速度、密度、壓力等;

7)判斷宏觀參量是否收斂;

8)如果宏觀參量收斂,輸出計算結果;否則重回第4步繼續(xù)計算,直至計算結果收斂。

使用該算法模擬狹窄髂外動脈,將狹窄處上下端的血壓設置為計算輸出結果進行格子無關性檢驗,發(fā)現(xiàn)在達到格子單元長度l=2.99×10-4m后,計算輸出結果的變化很小,因此可用該格子單元長度模擬狹窄髂外動脈,此時對應的格子數(shù)量為142×198×258。對于一個典型的髂外動脈模擬任務,在計算機CPU為Intel Xeon E5-2683 v3,GPU為NVIDIA Tesla P100的配置下,需要約20 min的時間可完成模擬任務。

2 結果與討論

2.1 驗證VLBM的可行性

驗證VLBM用于狹窄髂外動脈模擬的可行性是第一步,也是至關重要的一步。我們將模擬獲得的狹窄處上下端血壓與IBP結果做對比以進行驗證,IBP結果來自臨床測量,如圖6的紅色虛線所示。在使用計算結果之前,需要確認時間收斂性,當宏觀參數(shù)為每個相鄰的心動周期收斂時,設置收斂條件。圖6中的藍色實線為使用自編的VLBM算法計算的髂外動脈狹窄處上下端血壓曲線。從對比結果可以看出,2條代表模擬血壓和IBP的曲線具有較高的吻合度。定量地檢測血壓的具體參數(shù),即收縮壓Ps、舒張壓Pd和平均血壓Pmean。收縮壓是血壓曲線中的最大值,舒張壓是血壓曲線中的最小值,平均血壓定義為Pmean=(Ps+2Pd)/3[22]。表1為髂外動脈狹窄處上下端IBP與VLBM血壓參數(shù)對比及相對誤差,其中相對誤差為0.455%～1.585%。因此,無論是定性還是定量的數(shù)據(jù)結果均表明將VLBM用于狹窄髂外動脈的血流動力學分析是可行的。

圖6 髂外動脈狹窄處上下端IBP與VLBM血壓對比

表1 髂外動脈狹窄處上下端IBP與VLBM血壓參數(shù)對比及相對誤差

2.2 支架植入前后血流動力學分析

分別選取支架植入前狹窄率為70%和支架植入后還原的髂外動脈流場,觀察狹窄對流場的影響。圖4(a)的6個抽樣點中,t=0.20 s時的入口血流流速最大,狹窄處上下端的流場變化最明顯,可以更好地顯示流場的變化。因此以t=0.20 s為例,圖7顯示了支架植入前后髂外動脈的流線,將狹窄處放大后做支架植入前后的速度截面圖如圖8所示,從中可以看出,與支架植入后還原的髂外動脈相比,血流流經(jīng)狹窄處后,流速明顯增大并向動脈中心處匯聚。與髂外動脈狹窄處血流速度有明顯的增加相比,髂內(nèi)動脈血流速度的變化可以忽略不計。

圖7 髂外動脈支架植入前后的流線圖

圖8 髂外動脈支架植入前后的血流速度截面圖

壁面切應力(wall shear stress,WSS)是血流作用于動脈內(nèi)壁產(chǎn)生的機械力之一[23],對動脈粥樣硬化狹窄的形成起著重要作用。有研究表明[24-27],在早期低WSS血流分布區(qū)域更易于動脈粥樣硬化斑塊的形成,這可能導致動脈管腔逐漸變窄,斑塊侵入腔內(nèi),WSS會先恢復正常,之后開始升高,WSS越高,對動脈壁施加的拉伸力越大,從而增大動脈破裂的風險。將圖7中髂外動脈狹窄(還原)部分放大后如圖9和圖10所示。圖9(a)顯示了狹窄髂外動脈的WSS分布,其中時間點取自圖4(a)中的6個采樣時間點。狹窄部位WSS較大,高WSS的分布范圍隨血流速度的變化而變化。當t=0.20 s血流接近收縮期峰值時,高WSS的分布范圍更廣,隨后逐漸減小,與一個心動周期內(nèi)血流的變化趨勢一致。支架植入后,每個時間點的WSS均減小。綜上所述,狹窄將導致高WSS,增加動脈破裂的風險。

圖9 髂外動脈支架植入前后的WSS分布

圖10 髂外動脈支架植入前后的壓力云圖

圖10為髂外動脈支架植入前后的壓力云圖。支架植入后血壓隨時間變化,還原處上下端壓降約為0.047 6 mmHg,可以忽略不計。相反,支架植入前的狹窄處有一個明顯的壓降,此外,血流量越大,壓降越大。經(jīng)計算,支架植入前狹窄處的平均壓降約為8.838 5 mmHg。從定量角度來分析,如圖11所示,支架植入后,收縮壓在狹窄(還原)處上端降低,同時在下端升高,以此減小壓降,而支架植入前后舒張壓變化不大。因此,狹窄將導致壓降增大,尤其是在流速較高的收縮期。

圖11 髂外動脈支架植入前后血壓對比曲線

3 結 論

本研究使用自編的GPU并行加速VLBM算法,對狹窄髂外動脈進行了數(shù)值模擬和血流動力學分析,利用CTA影像重構了髂外動脈結構模型。計算得到的髂外動脈血壓與IBP結果具有較高的吻合度,驗證了VLBM模擬狹窄髂外動脈的可行性。此外,將支架植入狹窄處,還原動脈,經(jīng)對比分析得出以下結論:

1)與支架植入后恢復的動脈相比,血流流經(jīng)狹窄處后,流速明顯增大并向動脈中心處匯聚。髂外動脈狹窄對髂內(nèi)動脈血流速度影響可忽略不計。

2)狹窄處WSS較高,WSS高分布區(qū)的大小隨血流速度的變化而變化,增加了動脈破裂的風險。

3)狹窄處有明顯的壓降,尤其是在流速較高的收縮期,因此壓降會隨著血流速度的增大而增大。

上述研究結果對探討血壓、WSS、狹窄率之間的關系及相關的動脈模擬研究具有參考意義,本研究中的髂外動脈血流動力學分析可為相關髂外動脈疾病的臨床無創(chuàng)計算輔助決策提供數(shù)據(jù)支持和理論依據(jù)。