劉 瑩， 章德發(fā)， 殷艷飛， 張智亮

腦卒中已成為危害人類健康的全球重大公共衛(wèi)生問題，所致死亡人數(shù)2/3發(fā)生于發(fā)展中國家。我國心血管疾病患者達2.3億，其中腦卒中患者超過700萬［1］，每年用于治療腦卒中費用達數(shù)百億元，對我國社會經(jīng)濟發(fā)展構成嚴重負擔。腦卒中是一種受環(huán)境與遺傳因素共同作用而緩慢形成的復雜疾病，20%由頸動脈粥樣斑塊性狹窄引起［2］。高血壓、糖尿病和吸煙等因素已證實對粥樣斑塊形成具有重要影響［3-4］，但這些因素難以說明粥樣斑塊分布部位的特殊性。粥樣斑塊好發(fā)于血管幾何結構復雜部位，表明局部血流動力學參數(shù)在粥樣斑塊形成與發(fā)展中起到十分重要作用［5-6］。劉桂英等［7］通過計算流體力學方法研究右冠狀動脈內血流動力學與粥樣硬化斑塊形成的關系，結果表明局部低壁面切應力和低壁面壓力區(qū)域與粥樣硬化斑塊高發(fā)區(qū)域重合。Assemat等［8］對血流動力學參數(shù)與粥樣斑塊形成、發(fā)展和破裂過程進行分析，發(fā)現(xiàn)復雜血流和低壁面切應力有利于粥樣斑塊形成。目前多數(shù)血流動力學研究未考慮血管壁作用，分析非穩(wěn)態(tài)血流時管壁具有吸收和釋放血流能量的功能，有利于血液循環(huán)流動。因此，在流固耦合作用下充分了解頸動脈粥樣斑塊與頸動脈內血流動力學分布關系，對腦卒中早期預防與治療具有關鍵意義。

本文采用計算流體力學方法對流固耦合作用下頸動脈內一個心動周期典型時刻點血流動力學進行數(shù)值分析，從血流動力學特征參數(shù)方面分析心動收縮期和心動舒張期對血流分布特性的影響，探索頸動脈粥樣斑塊形成與發(fā)展機制。

1 材料與方法

1.1 實驗材料

采用UG NX 8.0軟件構建頸動脈幾何模型（圖1），包括頸總動脈（CCA）、頸內動脈（ICA）和頸外動脈（ECA）；CCA內徑為5.3 mm，ICA內徑為3.8mm，ECA內徑為3.2mm。頸動脈具有人體生理參數(shù)的幾何彎曲與錐度特征，頸動脈內血液密度為1.055×103kg/m3，流場內生理壓力為 1.333 9× 104Pa［9-10］。非穩(wěn)態(tài)血液沿CCA軸向流入，入口血液速度隨時間變化，ICA與ECA軸向為血流出口，其相對壓力設為零［11-12］，選擇Carreau-Yasuda黏度模型作為非牛頓血液模型［13］：

上式中，μ0＝22×10-3Pa·s，μ∞＝2.2×10-3Pa·s，λ＝0.110 s，n＝0.392，z＝0.644，γ為切應變率，黏度單位為Pa.s。頸動脈管壁設為均勻、各項同性線彈性材料，管壁密度為1.15×103kg/m3，壁厚均為0.5 mm，彈性模量為0.5mPa，泊松比為0.45，管壁兩端面軸向位移假設為零［14］。采用ICEM軟件對頸動脈血液和管壁進行體網(wǎng)格劃分。為提高流固耦合面計算精度，對流固耦合面血液和管壁分別采用逐層加密的3棱柱網(wǎng)格劃分，其余部分采用4面體網(wǎng)格劃分。

圖1 頸動脈模型

1.2 實驗方法

1.2.1 血流控制方程 血液為黏性、不可壓縮的非牛頓流體，血流控制方程為三維非穩(wěn)態(tài)血液流動Navier-Stokes方程［15］：

上式中，u為血流速度，ρ為血液密度，P為頸動脈流場內壓力，τ為應力張量。

1.2.2 頸動脈管壁控制方程 頸動脈管壁運動控制方程［16］如下：

上式中，ρs，ɑs，σs分別為頸動脈管壁密度、管壁質點加速度和管壁應力張量。血液與頸動脈管壁流固耦合面應滿足以下條件：

上式中，d、n、σ分別為位移、邊界法向和應力張量，f和s分別表示血液與頸動脈管壁。

選取頸動脈上6個特征點進行分析（圖1），P1、P2、P3分別為分叉部位附近CCA、ICA、ECA管徑軸線上一點，P4、P5、P6分別為分叉部位附近CCA、ICA、ECA外側流固耦合面上的一點。一個心動周期為0.8 s，非穩(wěn)態(tài)迭代計算時間步長為0.02 s，得到一個心動周期典型時刻血流動力學分布。

2 結果

2.1 血流分析

頸動脈內血流分布分析結果顯示，血流加速期（0.05 s）頸動脈分叉部位血流速度較大，ICA和ECA內血流速度與CCA內血流速度相近，方向與管壁平行；血流正峰值期（0.19 s）入口血流速度達到最大，頸動脈內血流速度整體達到最大，且ECA入口處左側存在血流停滯區(qū)，該區(qū)域內血流速度極低，血流流線少，血細胞不斷與管壁碰撞與摩擦，管壁細胞結構和功能發(fā)生改變，且血液中脂質等大分子易沉積，為粥樣斑塊形成提供有利生理環(huán)境；收縮期（0.3 s）血流處于心動收縮期末，血流速度快速減低，整個頸動脈內血流速度減低，且血流停滯區(qū)面積逐漸增大；血流負峰值期（0.54 s）血流處于心動舒張期，入口血流速度很低，頸動脈管壁回彈作用為血流流動提供動力，整個頸動脈內血流流線分布均勻、速度較低，血流停滯區(qū)面積很??；與心動舒張期相比，心動收縮期快速變化的血流對頸動脈管壁作用效果較大，ECA入口處左側形成的血流停滯區(qū)面積很大，對粥樣斑塊形成影響較明顯。（圖2）

圖2 典型時刻頸動脈內血流分布

下頸動脈內血流速度分析結果顯示，各特征點血流速度均隨入口血流速度增大而增大；心動收縮期P3血流速度變化不明顯，速度值均很小，而P1、P2血流速度變化趨勢一致，且血流速度變化幅度與變化率均較大；心動舒張期P1、P2、P3血流速度變化趨勢相近，速度值均較??；整個心動周期內P3血流速度均很小，與ECA入口處左側血流停滯區(qū)相符，進一步表明ECA入口處左側易形成粥樣斑塊。（圖3）

2.2 壁面壓力與切應力分析

壁面壓力分析結果顯示，P4、P5、P6相對壁面壓力在心動收縮期隨入口血流速度增大而逐漸增大，ICA和CCA（P4、P6）相對壁面壓力明顯高于ECA左側（P5）；ICA內相對壁面壓力最大，CCA內相對壁面壓力其次，ECA內相對壁面壓力最?。恍膭邮鎻埰赑4、P5、P6相對壁面壓力為負值，形成“負壓”效應，各特征點處形成的“負壓”效應相近；低壁面壓力作用下 “負壓”效應使血管內血流動力受到很大影響，血流速度變慢，引起腦部供血不足，導致缺血性腦卒中發(fā)生。（圖4）

圖3 特征點血流速度隨時間變化曲線

圖4 特征點相對壁面壓力隨時間變化曲線

壁面切應力分析結果顯示，P4、P5、P6壁面切應力隨入口血流速度增大而增大，心動收縮期ICA、CCA（P4、P6）壁面切應力明顯大于ECA左側（P5），而心動舒張期P4、P5、P6壁面切應力變化很小，且趨勢一致；與P4、P6壁面切應力相比，P5壁面切應力在整個心動周期均處于低壁面切應力范圍，且與血流停滯區(qū)位置吻合；低壁面切應力作用下近壁面血流正常流態(tài)被破壞，血細胞與管壁不斷摩擦與碰撞，并激活血小板活性，促進了管壁內膜過度增生。（圖5）

圖5 特征點壁面切應力隨時間變化曲線

2.3 頸動脈管壁形變量分析

頸動脈管壁形變量分析結果顯示，P4、P5、P6管壁形變量均隨入口血流速度增大而增大，且P4、P5、P6形變量依次減??；心動收縮期P4、P5、P6形變量變化較明顯，變化趨勢一致，且在0.19 s時達到最大，但心動舒張期P4、P5、P6形變量變化很小，且軌跡幾乎重合；與心動舒張期相比，心動收縮期ICA、ECA內血流速度大，血流對管壁作用增強，引起P4、P5管壁形變量較大。（圖6）

圖6 特征點頸動脈管壁形變量隨時間變化曲線

2.4 頸動脈管壁von Mises等效應力分析

頸動脈管壁von Mises等效應力分析結果顯示，一個心動周期內不同時刻頸動脈分叉部和ECA入口處左側von Mises等效應力始終較其它部位大，易發(fā)生應力集中，血管壁發(fā)生破裂的危險系數(shù)明顯增加；血流加速期（0.05 s）血流速度較小，頸動脈內von Mises等效應力值整體較小；血流正峰值期（0.19 s）血流速度達到最大，頸動脈內von Mises等效應力明顯增大，尤其在整個CCA、ECA、ICA入口部位特別大；心動收縮期（0.3 s）和血流負峰值期（0.54 s）血流處于心動舒張期，心臟停止對外供血，頸動脈內血流主要依靠彈性管壁回彈作用為血流提供動力，此時管壁形變量逐漸減小，管壁von Mises等效應力快速減小。（圖7）

3 結論

一個心動周期內，流固耦合作用下ECA入口處左側（P5）存在血流停滯區(qū)，且心動收縮期內血流停滯區(qū)面積明顯大于心動舒張期，血流停滯區(qū)內血流速度極低，血流流線少，血細胞不斷與管壁碰撞與摩擦，管壁細胞結構和功能發(fā)生改變，且血液中脂質等大分子易沉積，為粥樣斑塊形成提供有利生理環(huán)境。

圖7 典型時刻頸動脈von Mises等效應力分布

整個心動周期內，ECA入口處左側（P5）存在低壁面壓力和低壁面切應力區(qū)，與血流停滯區(qū)吻合。低壁面壓力作用下形成了“負壓”效應，血流動力受到影響，血流速度變慢，引起腦部供血不足，導致缺血性腦卒中發(fā)生。低壁面切應力作用下，近壁面血流的正常流態(tài)被破壞，血小板活性增強，促進內膜過度增生。

流固耦合作用下，與心動舒張期相比，心動收縮期內頸動脈管壁形變量和von Mises等效應力值均較大，且局部存在突變。ECA（P5）管壁形變量和von Mises等效應力值較P4和P6大，易引起應力集中，增大血管破裂風險。

［1］ Wan LH,Zhao J,Zhang XP,et al.Stroke prevention knowledge and prestroke health behaviors among hypertensive stroke patients in mainland China［J］.JCardiovasc Nurs,2014,29:E1-E9.

［2］ Linfante I,Andreone V,Akkawi N,et al.Internal carotid artery stenting in patients over 80 years of age:single-center experience and review of the literature［J］.JNeuroimaging,2009,19:158-163.

［3］ Wentzel JJ,Chatzizisis YS,Gijsen FJ,et al.Endothelial shear stress in the evolution of coronary atherosclerotic plaque and vascular remodelling:current understanding and remaining questions［J］.Cardiovasc Res,2012,96:234-243.

［4］ 彭應龍，宋 莉，佟小強，等.頸動脈支架成形術血流動力學改變的影響因素分析［J］.介入放射學雜志，2013，22：535-539.

［5］ 王曉曦，劉宏斌，胡小忠，等.冠狀動脈局部血流動力學參數(shù)分析及斑塊預警研究［J］.解放軍醫(yī)學院學報，2013，34：617-620.

［6］ 江國權，方興根，徐善水，等.顱內動脈瘤破裂的血流動力學和形態(tài)學因素［J］.介入放射學雜志，2014，23：1109-1113.

［7］ 劉桂英.右冠狀動脈的血流動力學數(shù)值模擬［D］.廣州：南方醫(yī)科大學，2014.

［8］ Assemat P,Armitage JA,Siu KK,et al.Three-dimensional numerical simulation of blood flow in mouse aortic arch around atherosclerotic plaques［J］.Appl Math Model,2014,38:4175-4185.

［9］ Chaichana T,Sun Z,Jewkes J.Computation of hemodynamics in the left coronary artery with variable angulations［J］.JBiomech, 2011,44:1869-1878.

［10］Blagojevic M,Nikolic A,Zivkovic M,et al.Role of oscillatory shear index in predicting the occurrence and development of plaque［J］.JSerbian Soc ComputMech,2013,7:29-37.

［11］楊金有，徐躍平，俞 航，等.人體主動脈弓內非牛頓血液流動數(shù)值模擬分析［J］.中國醫(yī)學物理學雜志，2011，28：2422-2425.

［12］Karmonik C,Bismuth JX,Davies MG,et al.Computational hemodynamics in the human aorta:a computational fluid dynamics study of three cases with patient-specific geometries and inflow rates［J］.Technol Health Care,2008,16:343-354.

［13］Chen J,Lu XY.Numerical investigation of the non-Newtonian pulsatile blood flow in a bifurcation model with a non-planar branch［J］.JBiomech,2006,39:818-832.

［14］Yang JY,Hong Y.Numerical simulations of the non-newtonian blood blow in human abdominal artery based on reverse engineering［J］.Appl Mech Mater,2012,140:195-199.

［15］Lee CJ,Zhang Y,Takao H,et al.The influence of elastic upstream artery length on fluid-structure interaction modeling:A comparative study using patient-specific cerebral aneurysm［J］. Med Eng Phys,2013,35:1377-1384.

［16］Valencia A,Baeza F.Numerical simulation of fluid-structure interaction in stenotic arteries considering two layer nonlinear anisotropic structuralmodel［J］.Int Communi Heat Mass Transf, 2009,36:137-142.