頸動(dòng)脈分叉的非穩(wěn)態(tài)數(shù)值模擬分析

王汝良，胡霖霖，郭金興，張凱旋，張臘梅，戴志穎，陳廣新*

?

頸動(dòng)脈分叉的非穩(wěn)態(tài)數(shù)值模擬分析

王汝良1，胡霖霖2，郭金興3，張凱旋2，張臘梅2，戴志穎2，陳廣新2*

（1. 牡丹江醫(yī)學(xué)院附屬紅旗醫(yī)院影像科，黑龍江 牡丹江 157011；2. 牡丹江醫(yī)學(xué)院附屬紅旗醫(yī)院介入科， 黑龍江 牡丹江 157011；3. 牡丹江醫(yī)學(xué)院醫(yī)學(xué)影像學(xué)院，黑龍江 牡丹江 157011）

本課題旨在研究頸動(dòng)脈分叉心動(dòng)周期內(nèi)的血液動(dòng)力學(xué)參數(shù)，如頸動(dòng)脈的壁切應(yīng)力、壓力、血流速度及渦流變化及各力學(xué)參數(shù)對(duì)頸動(dòng)脈斑塊形成的影響。選取患者的頸動(dòng)脈CT血管成像數(shù)據(jù)構(gòu)建有限元模型，通過(guò)計(jì)算流體力學(xué)軟件仿真計(jì)算。在設(shè)定了初始條件的情況下，獲得了頸動(dòng)脈的血液動(dòng)力學(xué)參數(shù)，得出了頸動(dòng)脈壁面剪切應(yīng)力及壁面壓力的的血流規(guī)律，并對(duì)頸動(dòng)脈的壁切應(yīng)力、壓力、血流速度及渦流變化進(jìn)行了詳細(xì)的分析。頸動(dòng)脈竇部位的壁面剪切應(yīng)力最小，血管壁的壁面壓力存在“負(fù)壓”效應(yīng)，負(fù)壓效應(yīng)導(dǎo)致血流速度趨緩、血管直徑減小、血流量減小、血流阻力增加。在脈動(dòng)周期一定時(shí)刻血流在頸動(dòng)脈的分叉附近會(huì)形成渦流，渦流區(qū)內(nèi)血液低速流動(dòng)且存在回流及二次流等復(fù)雜血流，實(shí)驗(yàn)結(jié)果符合血液的流體力學(xué)運(yùn)動(dòng)方程的理論分析。

頸動(dòng)脈；有限元模型；渦流；血流動(dòng)力學(xué)分析

0 引言

心腦血管疾病嚴(yán)重威脅人類健康，給社會(huì)和家庭造成了巨大的生活與財(cái)務(wù)負(fù)擔(dān)[1]，是老年人致死的主要病因之一，危害性僅次于惡性腫瘤。我國(guó)腦卒發(fā)病率的增長(zhǎng)速度高達(dá)8.7%，據(jù)文獻(xiàn)記載，10~20%的腦卒中是由于血管內(nèi)粥樣硬化的斑塊脫落所引起的[2,3]，研究發(fā)現(xiàn)，動(dòng)脈硬化及血管內(nèi)膜的增厚，血液運(yùn)輸通道變狹窄，血液流速增加，對(duì)管壁的沖擊增大，管壁破損后形成斑塊，最終斑塊 脫落會(huì)形成血栓，而頸動(dòng)脈處斑塊脫落可導(dǎo)致腦卒中[4,5]。病理解剖發(fā)現(xiàn)斑塊具有病灶選擇性，好發(fā)于頸動(dòng)脈竇及頸動(dòng)脈分叉等部位。管壁壓力、血流速度、回流、渦流、二次流、震蕩剪切因子（Oscillatory shear factor，OSI）及低平均壁面剪切應(yīng)力（Low average wall shear stress，TAWSS）是斑塊的血流動(dòng)力學(xué)成因。粒子滯留時(shí)間（Particle retention time，RRT）的長(zhǎng)短將會(huì)直接影響動(dòng)脈粥樣硬化的發(fā)生和發(fā)展[6,10]。Malek等人對(duì)血管的壁面切應(yīng)力進(jìn)行了詳細(xì)的研究，研究結(jié)果顯示壁面切應(yīng)力的臨界值低于0.4 Pa時(shí)會(huì)發(fā)生動(dòng)脈粥樣硬化[8,9]。因此，研究頸動(dòng)脈的血液流體力學(xué)因素具有重要的臨床意義。

本實(shí)驗(yàn)利用實(shí)際病人的頸動(dòng)脈CT血管成像（CT angiography，CTA）數(shù)據(jù)構(gòu)建有限元模型，通過(guò)計(jì)算流體力學(xué)（Computational fluid mechanics，CFD）軟件進(jìn)行仿真計(jì)算，獲得頸動(dòng)脈在心動(dòng)周期內(nèi)的血液動(dòng)力學(xué)參數(shù)，分析頸動(dòng)脈的壁切應(yīng)力、壓力、血流速度及渦流變化情況。

1 材料與方法

1.1 資料

圖像數(shù)據(jù)：數(shù)據(jù)來(lái)源于一男性64歲患者，選取80張斷層厚度為0.5 mm，圖像矩陣為512×512的頸部CTA斷層掃描Dicom原始數(shù)據(jù)。儀器應(yīng)用東芝64排螺旋CT，CT掃描參數(shù)：電壓120 V，電流250 mA。

圖像后處理工作站：戴爾Precision T7810，Xeon E5-2609 v3 CPU、nVIDIA Quadro k2200圖像顯卡；圖像后處理及計(jì)算流體力學(xué)分析軟件：Mimics 20.0、3-matic 12.0、Ansys 16.0。

1.2 方法

1.2.1 建立有限元模型

將80張Dicom格式CTA數(shù)據(jù)文件導(dǎo)入Mimics中應(yīng)用閾值分割、動(dòng)態(tài)區(qū)域增長(zhǎng)等方法，初步重建建頸動(dòng)脈三維模型。將重建后的三維模型導(dǎo)入3-matic軟件中進(jìn)行光滑處理，去除細(xì)小分支，切平出口、入口面，最后完成頸動(dòng)脈三維幾何模型（圖1），該模型包括頸總動(dòng)脈（Common carotid artery，CCA）、頸內(nèi)動(dòng)脈（Internal carotid artery，ICA）、頸外動(dòng)脈(External carotid artery，ECA)。將該模型以stl格式導(dǎo)入Ansys icem cfd劃分網(wǎng)格，邊界層實(shí)施五層棱柱加密以保證計(jì)算的精度。網(wǎng)格類型為四面體，網(wǎng)格單元節(jié)點(diǎn)共計(jì)218524個(gè)，網(wǎng)格單元為121872個(gè)。將劃分之后的網(wǎng)格導(dǎo)入Ansys fluent 16.0中轉(zhuǎn)換為多面體網(wǎng)格，可以進(jìn)一步提升計(jì)算精度。圖1中A、B、C、D、E、F為選定的位置觀測(cè)點(diǎn)。

圖1 頸動(dòng)脈分叉幾何模型

1.2.2 血流控制方程

本研究假定血液為不可壓縮、粘性的牛頓流體，血管壁為剛性壁，控制方程為：

式中為速度矢量，p為流場(chǎng)壓力，ρ為血流密度，μ為動(dòng)力粘度。本文中，μ = 0.00356 Pa×s，ρ = 1060 kg/m3。

1.2.3 邊界條件

本文研究采用瞬態(tài)計(jì)算，入口速度為隨時(shí)間變化的脈動(dòng)曲線（圖2）。出口壓力設(shè)定為0，血管壁無(wú)滑動(dòng)，各速度分量為0，不考慮重力影響。計(jì)算指定時(shí)間步長(zhǎng)為0.01 s，為獲得穩(wěn)定周期結(jié)果，共進(jìn)行三個(gè)心動(dòng)周期計(jì)算，對(duì)第三個(gè)周期的計(jì)算結(jié)果進(jìn)行分析。

圖2 入口的血流波形

2 結(jié)果分析

2.1 壁面剪切應(yīng)力（WSS）分析

研究結(jié)果表明[1,3]，動(dòng)脈粥樣硬化一般發(fā)生在血管壁面的剪切應(yīng)力低于1 Pa的區(qū)域。圖3是位置觀測(cè)點(diǎn)A、B、C、D的壁面切應(yīng)力曲線圖。圖3結(jié)果顯示，頸動(dòng)脈竇處（A點(diǎn)）的壁面剪切應(yīng)力震蕩范圍為0~0.12 Pa，始終保持在較低的波動(dòng)范圍。頸動(dòng)脈B點(diǎn)的壁面剪切應(yīng)力在整個(gè)心動(dòng)周期內(nèi)，始終遠(yuǎn)高于A點(diǎn)。頸動(dòng)脈D點(diǎn)的壁面剪切應(yīng)力在整個(gè)心動(dòng)周期內(nèi)，始終遠(yuǎn)高于C點(diǎn)。B、C、D、E各點(diǎn)的壁切應(yīng)力變化趨勢(shì)與A點(diǎn)保持一致，并隨時(shí)間震蕩。圖4是一個(gè)心動(dòng)周期內(nèi)典型時(shí)刻壁面剪切力的變化情況，圖4表明，頸動(dòng)脈竇部位的壁面剪切應(yīng)力分布始終是整個(gè)頸動(dòng)脈區(qū)域最小的，而臨床研究已經(jīng)證實(shí)，動(dòng)脈粥樣硬化易發(fā)生在低壁面切應(yīng)力處[7]。

圖3 A、B、C、D各點(diǎn)壁面剪切應(yīng)力

圖4 不同時(shí)刻剪切應(yīng)力變化情況

2.2 壁面壓力分析

圖5為頸動(dòng)脈B、E、F壁面壓強(qiáng)隨時(shí)間的變化曲線，圖6為整個(gè)頸動(dòng)脈各個(gè)典型時(shí)刻的壓強(qiáng)對(duì)時(shí)間的變化曲線。實(shí)驗(yàn)結(jié)果顯示，在時(shí)間段=0.08 s到=0.3 s（心動(dòng)收縮減速期），血管壁包括B、E、F點(diǎn)的壁面壓力存在負(fù)值現(xiàn)象即“負(fù)壓”效應(yīng)，與入口血流速度降低正相關(guān)。在低壁面壓力作用下，這種“負(fù)壓”效應(yīng)導(dǎo)致血流速度趨緩、血管直徑減小、血流量減小、血流阻力增加及腦部供血不足，繼而引發(fā)缺血性腦卒中，或誘發(fā)局部血管痙攣，致使腦缺血進(jìn)一步加重。此外，血流速度慢將導(dǎo)致血液中的脂質(zhì)附著血管壁時(shí)間增長(zhǎng)，可促使粥樣斑塊的進(jìn)一步發(fā)展[12]，當(dāng)初始條件完全設(shè)定的情況下該實(shí)驗(yàn)結(jié)論與血液的流體力學(xué)運(yùn)動(dòng)方程的理論分析和求解完全吻合。

圖5 B、E、F各點(diǎn)壁面壓力變化情況

圖6 各典型時(shí)刻壓力變化情況

2.3 血流分析

頸動(dòng)脈血流線圖結(jié)果顯示（圖5），在心動(dòng)收縮加速期，頸總動(dòng)脈（CCA）、頸內(nèi)動(dòng)脈（ICA）、頸外動(dòng)脈（ECA）的血流速度相差不多，方向與軸向平行。在=0.08 s時(shí)，入口速度達(dá)到峰值，CCA、ICA、ECA的血流速度都達(dá)到最高值。在ICA的A、B（分叉處）、C點(diǎn)附近區(qū)域速度較低，血流停滯面積較小，血細(xì)胞與管壁之間存在碰撞和摩擦。在心動(dòng)收縮期內(nèi)血流速度劇烈變化，對(duì)頸動(dòng)脈結(jié)構(gòu)和功能會(huì)產(chǎn)生影響，有利于血液中脂質(zhì)等成分沉積，為斑塊形成提供了條件[12]。在脈動(dòng)收縮期末（=0.3 s），頸動(dòng)脈整體血流速度最低；在舒張期，血流速度保持管壁產(chǎn)生作用，尤其是頸動(dòng)脈竇附近區(qū)域、分叉部位，更有利于斑塊的形成、發(fā)展[14,15]。

圖7是不同典型時(shí)刻的血流流線圖。=0.06 s，=0.08 s為心動(dòng)周期中的收縮上升期，=0.09s，= 0.1 s，=0.12s，=0.16 s，=0.21 s，=0.3 s為心動(dòng)周期的收縮下降期。其余時(shí)刻為心動(dòng)周期的舒張期。從收縮下降期開(kāi)始，渦流逐漸形成，在0.3s時(shí)刻渦流最大，渦流位置發(fā)生在頸動(dòng)脈竇部附近。從舒張期開(kāi)始，渦流強(qiáng)度逐漸降低。研究結(jié)果顯示，渦流區(qū)內(nèi)血流速度小且存在回流及二次流等復(fù)雜血流，血液中脂質(zhì)及纖維蛋白易沉積于此，有利于斑塊的形成。此外，血流速度大小、方向的頻繁改變對(duì)斑塊的形成與發(fā)展也有促進(jìn)作用[13,16]。

3 結(jié)論

本課題對(duì)頸動(dòng)脈在心動(dòng)周期內(nèi)的血液動(dòng)力學(xué)參數(shù)，如頸動(dòng)脈的壁切應(yīng)力、壓力、血流速度及渦流變化進(jìn)行了研究，同時(shí)探討了各力學(xué)參數(shù)對(duì)頸動(dòng)脈斑塊形成的影響，實(shí)驗(yàn)結(jié)果與血液的流體力學(xué)運(yùn)動(dòng)方程的理論分析和求解完全吻合。研究結(jié)果證實(shí)，頸動(dòng)脈竇部位的壁面剪切應(yīng)力分布是頸動(dòng)脈區(qū)域最小的，動(dòng)脈粥樣硬化易發(fā)生在低壁面切應(yīng)力處。血管壁的壁面壓力存在負(fù)值現(xiàn)象即“負(fù)壓”效應(yīng)。這種“負(fù)壓”效應(yīng)導(dǎo)致血流速度趨緩、血管直徑減小、血流量減小、血流阻力增加及腦部供血不足，繼而引發(fā)缺血性腦卒中，或誘發(fā)局部血管痙攣。在心動(dòng)收縮期內(nèi)血流速度劇烈變化，對(duì)頸動(dòng)脈結(jié)構(gòu)和功能會(huì)產(chǎn)生影響，為斑塊形成提供了條件。研究結(jié)果顯示，渦流區(qū)內(nèi)血流速度小且存在回流及二次流等復(fù)雜血流，血流速度大小、方向的頻繁改變對(duì)斑塊的形成與發(fā)展有一定的促進(jìn)作用。

圖7 不同時(shí)刻流線圖

[1] 袁煒, 陳忠利. 頸動(dòng)脈血液動(dòng)力學(xué)的數(shù)值模擬[J]. 中國(guó)組織工程研究. 2014, 18(42): 6785-6786.

[2] Razavi A, Shirani E, Sadeghi MR. Numerical simulation of blood pulsatile flow in a stenosed carotid artery using different rheological models[J]. J Biomech. 2011; 44(11): 2021-2030.

[3] 張耀楠, 李艷, 康雁. 基于CT影像的頸動(dòng)脈分叉血流動(dòng)力學(xué)特性分析[J]. 東北大學(xué)學(xué)報(bào)(自然科學(xué)版), 2014, 35(3): 356-357.

[4] 周志尊, 王佳美, 胡明成, 陳廣新, 劉陽(yáng), 董默, 孫強(qiáng), 孫鵬, 劉佳維. 頸動(dòng)脈三維重建及臨床應(yīng)用的研究[J]. 軟件, 2017, 38(8): 32-35.

[5] Bai-Nan X, Fu-Yu W, Lei L, Xiao-Jun Z, Hai-Yue J. Hemodynamics model of fluid-solid interaction in internal carotid artery aneurysms[J]. Neurosurg Rev, 2011, 34(1): 39-47.

[6] Karmonik C, Yen C, Diaz O, Klucznik R, Grossman RG, Benndorf G. Temporal variations of wall shear stress parameters in intracranial aneurysms-importance of patient-specific inflow waveforms for CFD calculations[J]. Acta neurochirurgica, 2010, 152(8): 1391-8.

[7] Watanabe T, Isoda H, Takehara Y, Terada M, Naito T, Kosugi T, Onishi Y, Tanoi C, Izumi T. Hemodynamic vascular biomarkers for initiation of paraclinoid internalcarotidarteryaneurysms using patient-specific computational fluid dynamic simulation based on magnetic resonance imaging[J]. Neuroradiology, 2018, 60(5): 545-555.

[8] Xu B, Zhong H, Duan S. Modeling of internalcarotidarteryaneurysm and blood flow simulation[J]. Technology and health care: official journal of the European Society for Engineering and Medicine, 2015, 23(1): 43-8.

[9] Ren G, Cao X, Wang D, Jiang P, Li Y, Pei B. hemodynamic numerical simulationofcarotidartery aneurysm before and after surgery based on CT date[J]. Journal of biomedical engineering, 2014, 31(2): 341-6.

[10] Wang F, Xu B, Sun Z, Wu C, Zhang X. Wall shear stress in intracranial aneurysms and adjacent arteries. Neural regeneration research[J]. 2013, 8(11): 1007-15.

[11] Lal BK, Beach KW, Sumner DS. Intracranial collateralization determineshemodynamicforces forcarotid plaque disruption. Journal of vascular surgery, 2011, 54(5): 1461-71.

[12] De Wilde D, Trachet B, Debusschere N, Iannaccone F, Swillens A, Degroote J, Vierendeels J, De Meyer GRY, Segers P. Assessment of shear stress related parameters in thecarotidbifurcation using mouse-specific FSI simulations[J]. Journal of biomechanics, 2016, 49(11): 2135-2142.

[13] Saho T, Onishi H. Evaluation of the impact ofcarotidarterybifurcationangle on hemodynamics by use of computational fluid dynamics: asimulationand volunteer study[J]. Radio-logical physics and technology, 2016, 9(2): 277-85.

[14] Gharahi H, Zambrano BA, Zhu DC, DeMarco JK, Baek S. Computational fluid dynamicsimulationof humancarotidarterybifurcationbased on anatomy and volumetric blood flow rate measured with magnetic resonance imaging[J]. International journal of advances in engineering sciences and applied math-ematics, 2016, 8(1): 40-60.

[15] Younis HF, Kaazempur-Mofrad MR, Chan RC, Isasi AG, Hinton DP, Chau AH, Kim LA, Kamm RD. Hemodynamics and wall mechanics in humancarotidbifurcationand its conseq-uences for atherogenesis: investigation of inter-individual vari-ation[J]. Bio-m-echanics and modeling in mechanobiology 2004, 3(1): 17-32.

[16] Younis HF, Kaazempur-Mofrad MR, Chung C, Chan RC, Kamm RD. Computational analysis of the effects of exercise on hemodynamics in thecarotidbifurcation[J]. Annals of biomed-ical engineering, 2003, 31(8): 995-1006.

Carotid artery Bifurcation Unsteady State Numerical Simulation Analysis

WANG Ru-liang1, HU Lin-lin2, GUO Jin-xing3, ZHANG Kai-xuan2, ZHANG La-mei2, DAI Zhi-ying2, CHEN Guang-xin2*

(1. Department of radiology, Hongqi Hospital Affiliated to Mudanjiang medical University, Mudanjiang Heilongjiang 157011; 2. Medical image college of Mudanjiang medical university，Mudanjiang Heilongjiang 157011; 3. Department of intervention, Hongqi Hospital Affiliated to Mudanjiang medical University, Mudanjiang Heilongjiang 157011)

The purpose of this project is to study carotid artery bifurcation hemodynamic parameters in the carotid cardiac cycle, such as wall shear stress, pressure, blood flow velocity and eddy current velocity and the influence of various mechanical parameters on the formation of carotid plaque. A finite element model was constructed based on patient’s CT angiographic data of the carotid artery. Through the simulation calculation of fluid mechanics software, the flow line was designed in the experimental to observe the blood turbulence.Results: Under the initial conditions, the hemodynamic parameters of the carotid artery were obtained, the shear stress of the carotid wall and the blood flow pattern were obtained, and the wall shear stress, pressure, blood flow velocity and eddy flow changes of the carotid artery were analyzed in detail. The wall shear stress of the carotid sinus is the smallest, and the wall pressure of the vascular wall has the effect of "negative pressure", which slow down the blood flow rate, decrease blood vessel diameter, and increase blood flow resistance. When the velocity reaches a specific value，the blood flow forms vortex near the carotid bifurcation. Eddy blood flow velocity in this area is small and there is reflux and secondary flow. The experimental results conform to the fluid mechanics equations of the blood of theoretical analysis.

Carotid artery; Finite element model; Eddy current; Hemodynamic analysis

TP319

A

10.3969/j.issn.1003-6970.2018.10.008

黑龍江省自然基金項(xiàng)目(No. H2015082)；黑龍江省教育廳項(xiàng)目(No.12541846)；黑龍江省大學(xué)生創(chuàng)業(yè)項(xiàng)目(No.201646)

王汝良(1969-)，教授、主任醫(yī)師，研究方向：影像與核醫(yī)學(xué)；胡霖霖(1993-)，碩士研究生，研究方向：醫(yī)學(xué)圖像處理；郭金興(1984-)，主管護(hù)師，研究方向：介入治療；張凱旋(1993-)，碩士研究生，研究方向：醫(yī)學(xué)圖像處理；張臘梅(1993-)，碩士研究生，研究方向：超聲診斷；戴志穎(1996-)，學(xué)生，研究方向：醫(yī)學(xué)圖像處理。

陳廣新(1978-)，講師，研究方向：人工智能與醫(yī)學(xué)圖像處理。牡丹江醫(yī)學(xué)院醫(yī)學(xué)影像學(xué)院醫(yī)學(xué)圖像處理教研室。

王汝良，胡霖霖，郭金興，等. 頸動(dòng)脈分叉的非穩(wěn)態(tài)數(shù)值模擬分析[J]. 軟件，2018，39（10）：36-41