基于CTA的冠狀動脈分叉病變血流動力學(xué)仿真分析

郭金興 陳廣新 包婉秋 張春霞 張美樂 劉陽

摘? 要： 針對不同分叉病變區(qū)域長度的個體化冠狀動脈的血流動力學(xué)指標(biāo)分布特征，進行個體化建模，探討不同分叉病變區(qū)域長度對冠狀動脈的影響。應(yīng)用計算流體力學(xué)仿真計算，獲得不同分叉病變區(qū)域長度的冠狀動脈血流動力學(xué)指標(biāo)分布特征。利用CFX CCL語言實現(xiàn)不同的血流動力學(xué)指標(biāo)參數(shù)。通過計算機建模與仿真，可實現(xiàn)CTA-STL模型-CFD網(wǎng)格-CFD仿真結(jié)果，結(jié)果為冠狀動脈分叉病變的血流動力學(xué)變化受分叉病變的長度影響。

關(guān)鍵詞： CTA;冠狀動脈分叉病變;血流動力學(xué)仿真

中圖分類號： TP391.4 ???文獻標(biāo)識碼： A??? DOI：10.3969/j.issn.1003-6970.2020.07.024

本文著錄格式：郭金興，陳廣新，包婉秋，等. 基于CTA的冠狀動脈分叉病變血流動力學(xué)仿真分析[J]. 軟件，2020，41（07）：120-125

CTA Based Simulation Analysis of Coronary Bifurcation Lesions

GUO Jin-xing¹， CHEN Guang-xin²， BAO Wan-qiu²， ZHANG Chun-xia²， ZHANG Mei-le²， LIU Yang^3*

（1. Hongqi Hospital of Mudanjiang Medical University， Mudanjiang， 157011， China; 2. Medical Imaging College of Mudanjiang MedicalUniversity， Mudanjiang， 157011， China; 3. Registrars office of Mudanjiang Medical University， Mudanjiang， 157011， China）

【Abstract】： According to the distribution characteristics of individual coronary artery Hemodynamics index with different length of bifurcation lesion area， the individual model was established to explore the influence of different length of bifurcation lesion area on coronary artery. The distribution characteristics of coronary artery Hemodynamics index were obtained by using computational fluid dynamics simulation calculation. Using CFX CCL language to implement different Hemodynamics parameters. Cta-stl model-CFD grid-CFD simulation results are realized by computer modeling and simulation. The results show that the Hemodynamics of coronary bifurcation lesions is affected by the length of bifurcation lesions.

【Key words】： CTA; Coronary bifurcation lesions; Hemodynamics simulation

0? 引言

冠狀動脈狹窄與心肌缺血是否存在相關(guān)性一直沒有定論^[1-3]，中等程度額的冠狀動脈狹窄病變可能會引起心肌缺血，冠狀動脈分叉病變的影響因素包括狹窄率、狹窄區(qū)域長度，所以冠狀動脈的形態(tài)學(xué)研究對弄清冠狀動脈的血流動力學(xué)特征變化具有重要的意義^[4-5]。本研究即基于不同狹窄區(qū)域長度的三維模型，應(yīng)用血流動力學(xué)的指標(biāo)包括壁面剪切應(yīng)力、

震蕩剪切因子、避免剪切應(yīng)力梯度、時間平均壁面剪切應(yīng)力等對冠狀動脈不同病變長度的LAD進行分析，探討病變長度對冠狀動脈的的血流動力學(xué)因素的影響關(guān)系。

1? 個體化冠狀動脈三維模型構(gòu)建

選取牡丹江醫(yī)學(xué)院附屬紅旗醫(yī)院冠狀動脈患者病例一例，CTA斷層影像數(shù)據(jù)總計402張，層厚0.5?mm，使用比利時醫(yī)學(xué)交互式影像控制系統(tǒng)（Mater?ialists Interative Medical Image Control System，MIMICS）高級分割工具（ADVANCED SEGMENT）模塊中的Coronary專用分割工具進行冠狀動脈分割（圖1），提取出左冠狀動脈蒙版（Mask），并經(jīng)三維模型計算（Calculate Part生成初步的三維模型并以stl格式導(dǎo)入正向工程軟件3-matic中進行光順、三維模型修復(fù)、三角面片劃分，最終獲取冠狀動脈三維模型，此模型病變長度為0（圖2），圖2中標(biāo)示了模型的入口、出口位置和病變區(qū)域。在此動脈模型基礎(chǔ)上構(gòu)建10 mm，15 mm病變（狹窄）模型。MIMICS是醫(yī)學(xué)逆向工程軟件，該軟件可實現(xiàn)斷層數(shù)據(jù)的提取、三位建模、測量等，3-matic medical版本是基于正向工程技術(shù)的一款三維模型的構(gòu)建、修復(fù)處理軟件。

2 ?計算方法

2.1 ?材料屬性與邊界條件

本研究的仿真計算是基于ANSYS CFX軟件，Ansys CFX是一款高性能計算流體動力學(xué)（CFD）軟件工具，能快速穩(wěn)健地提供準(zhǔn)確可靠的解決方案，適用于眾多CFD和多物理場應(yīng)用，作為世界上唯一采用全隱式耦合算法的大型商業(yè)軟件^[6-8]。算法上的先進性，豐富的物理模型和前后處理的完善性使ANSYS CFX在結(jié)果精確性，計算穩(wěn)定性，計算速度和靈活性上都有優(yōu)異的表現(xiàn)。假設(shè)血流密度為1066 kg/m³，血液為牛頓流體，動力粘度為0.0035 。冠狀動脈入口采用速度入口，入口速度曲線如圖3所示，出口壓力曲線如圖4所示，不考慮重力的影響。計算兩個周期，每個周期的時長為0.8 s，取最后一個周期的結(jié)果進行研究。由于血管內(nèi)壓力很小，設(shè)置流體域的壓力為0。雷諾數(shù)為Re=1430<2300，因此采用層流^[9-10]。

2.3 ?網(wǎng)格劃分

各模型的網(wǎng)格劃分采用非結(jié)構(gòu)化的四面體網(wǎng)格（ANSYS FLUENT MESHING劃分），為保證計算精度，邊界層采用5層加密（圖5所示），網(wǎng)格劃分后進行網(wǎng)格獨立性驗證，滿足精度要求。ANSYS FLUENT MESHING是全新的基于Ribbon風(fēng)格的界面，提高了操作的便捷性，改善了用戶的體驗，同時提供了基于包面方法的全自動腳本生成網(wǎng)格、基于ANSA集成FLUENT MESHING的網(wǎng)格生成、基于SCDM結(jié)合FLUENT MESHING等多種網(wǎng)格生成流程。應(yīng)用該軟件劃分的網(wǎng)格質(zhì)量較高，并提供獨有的多面體網(wǎng)格技術(shù)^[11-13]。

3 ?仿真計算結(jié)果分析

本計算結(jié)果圖像后處理分析采用ENSIGHT進行分析，ENSIGHT是美國CEI公司開發(fā)的一款具有尖端的科學(xué)工程可視化的圖像后處理軟件。

3.1 ?血流動力學(xué)參數(shù)選擇

各參數(shù)的實現(xiàn)采用CCL編碼開發(fā)，ccl語言通俗易懂，可開發(fā)性較高，本實驗程序編碼部分如下：

inrad1 = sqrt（area（）@INLET1/ pi）

inrad2 = sqrt（area（）@INLET2/ pi）

invel1 = max（0 [m s^-1]， invmax1 * （1.0 - （zxradius1/inrad1）^2））

invel2 = max（0 [m s^-1]， invmax2 * （1.0 - （zxradius2/inrad2）^2））

invmax1 = 2 * INLET1f / （ pi *（inrad1^2）*areaAve（density）@INLET1）

invmax2 = 2 * INLET2f / （ pi *（inrad2^2）*areaAve（density）@INLET2）

inxcen1 = areaAve（Global X Coordinate）@INLET1

inxcen2 = areaAve（Global X Coordinate） @INLET2

inycen1 = areaAve（Global Y Coordinate） @INLET1

inycen2 = areaAve（Global Y Coordinate） @INLET2

inzcen1 = areaAve（Global Z Coordinate） @INLET1

inzcen2 = areaAve（Global Z Coordinate） @INLET2

numaneurysms = 1

numsystoliccycles = 3

peaksystole = 0.13

systoliccyclelength = 0.8

visc = viscval * 1[Pa s]

viscval = 0.0035

vortmag2 = Vorticity X^2 + Vorticity Y^2 + Vorticity Z^2

zxradius1 = sqrt（（x-inxcen1）^2 + （y-inycen1）^2 + （z-inzcen1）^2）

zxradius2 = sqrt（（x-inxcen2）^2 + （y-inycen2）^2 + （z-inzcen2）^2）

END

ADDITIONAL VARIABLE： Qvar

Option = Definition

Tensor Type = SCALAR

Units = [s^-2]

Variable Type = Unspecified

END

ADDITIONAL VARIABLE： ViscDisp

Option = Definition

Tensor Type = SCALAR

Units = [s^-2]

Variable Type = Unspecified

END

ADDITIONAL VARIABLE： WSSField

Option = Definition

Tensor Type = SCALAR

Units = [Pa]

Variable Type = Specific

END

ADDITIONAL VARIABLE： WSSG

Option = Definition

Tensor Type = VECTOR

Units = [Pa m^-1]

Variable Type = Unspecified

END

ADDITIONAL VARIABLE： WSSxF

Option = Definition

Tensor Type = SCALAR

Units = [Pa]

Variable Type = Specific

END

ADDITIONAL VARIABLE： WSSyF

Option = Definition

Tensor Type = SCALAR

Units = [Pa]

Variable Type = Specific

END

ADDITIONAL VARIABLE： WSSzF

Option = Definition

Tensor Type = SCALAR

Units = [Pa]

Variable Type = Specific

END

MATERIAL： Blood

Material Group = User

Option = Pure Substance

PROPERTIES：

Option = General Material

EQUATION OF STATE：

Density = dens

Molar Mass = 1.0 [kg kmol^-1]

Option = Value

END

DYNAMIC VISCOSITY：

Dynamic Viscosity = visc

Option = Value

END

END

END

END

FLOW： Flow Analysis 1

SOLUTION UNITS：

Angle Units = [rad]

Length Units = [m]

Mass Units = [kg]

Solid Angle Units = [sr]

Temperature Units = [K]

Time Units = [s]

END

ANALYSIS TYPE：

Option = Transient

EXTERNAL SOLVER COUPLING：

Option = None

END

INITIAL TIME：

Option = Automatic with Value

Time = 0 [s]

END

TIME DURATION：

Option = Total Time

Total Time = 2.4 [s]

END

TIME STEPS：

Option = Timesteps

Timesteps = 0.008 [s]

END

END

DOMAIN： FLUIDdom

Coord Frame = Coord 0

Domain Type = Fluid

Location = Assembly

BOUNDARY： INLET1

Boundary Type = INLET

Location = INLET1

BOUNDARY CONDITIONS：

ADDITIONAL VARIABLE： WSSField

Option = Zero Flux

END

ADDITIONAL VARIABLE： WSSxF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE： WSSyF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

ADDITIONAL VARIABLE： WSSzF

Additional Variable Value = 0 [kg m^-1 s^-2]

Option = Value

END

FLOW REGIME：

Option = Subsonic

END

MASS AND MOMENTUM：

Normal Speed = invel1

Option = Normal Speed

END

END

END

ADDITIONAL VARIABLE： OSIfield

Additional Variable Value = （1 - ＼

（WSTaveMag/（WSSField.Trnavg+ 1e-15[Pa]）））/2

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： PressGauge

Additional Variable Value = pref + Pressure

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： WSSField

Kinematic Diffusivity = 1e-15 [m^2 s^-1]

Option = Poisson Equation

END

ADDITIONAL VARIABLE： WSSG

Option = Vector Algebraic Equation

Vector xValue = -（1.0-Normal X*Normal X）*WSSField.Gradient X ＼

-（0.0-Normal X*Normal Y）*WSS?Field.Gradient Y -（0.0-Normal X*Normal ＼

Z）*WSSField.Gradient Z

Vector yValue = -（0.0-Normal Y* Normal X）*WSSField.Gradient X ＼

-（1.0-Normal Y*Normal Y）*WSS?Field.Gradient Y -（0.0-Normal Y*Normal ＼

Z）*WSSField.Gradient Z

Vector zValue = -（0.0-Normal Z*Nor?mal X）*WSSField.Gradient X ＼

-（0.0-Normal Z*Normal Y）*WSS?Field.Gradient Y -（1.0-Normal Z*Normal ＼

Z）*WSSField.Gradient Z

END

ADDITIONAL VARIABLE： WSSxF

Kinematic Diffusivity = 1e-15 [m^2 s^-1]

Option = Poisson Equation

END

ADDITIONAL VARIABLE： WSSyF

Kinematic Diffusivity = 1e-15 [m^2 s^-1]

Option = Poisson Equation

END

ADDITIONAL VARIABLE： WSSzF

Kinematic Diffusivity = 1e-15 [m^2 s^-1]

Option = Poisson Equation

END

COMBUSTION MODEL：

Option = None

END

HEAT TRANSFER MODEL：

Option = None

END

THERMAL RADIATION MODEL：

Option = None

END

TURBULENCE MODEL：

Option = Laminar

END

END

END

本文選擇血流速度、時間平均壁面切應(yīng)力（time average wall shear stress，TAWSS）、平均壁面切應(yīng)力梯度（time average wall shear stress grade，TAWSSG）、剪切震蕩系數(shù)（oscillatory shear index，OSI）、壁面切應(yīng)力（wall shear stress，WSS）是指血液流動時在血管壁表面上引起的切向的動態(tài)摩擦力。對于脈動流，在一個心臟周期內(nèi)用每個節(jié)點上積分WSS量的值來計算TAWSS：

其中wss_i是瞬時剪切應(yīng)力矢量，T是周期的持續(xù)時間。同時提出了壁面切應(yīng)力梯度（wall shear stress，WSSG），可以更明顯觀測WSS數(shù)值變化，WSSG的時間平均值即平均壁面切應(yīng)力梯度（TAWSSG）：

在一個心動周期內(nèi)，震蕩剪切系數(shù)（OSI）可以描述WSS方向的變化程度，OSI值介于（0，0.5）之間，OSI數(shù)值越大，表示W(wǎng)SS方向的變化也就越大，但其值的大小與WSS沒有必然聯(lián)系，其表達式為：

3.2 ?病變長度對TAWSS的影響

本研究提取了最后一個周期的TAWSS云圖 （圖6），由TAWSS分布云圖可見，在分叉病變區(qū)域，無病變、10 mm病變長度，15 mm病變長度的病變部位高TAWSS區(qū)域依次遞增。在病變區(qū)域，高TAWSS區(qū)域的值都較小，普遍都在0.6左右，而非病變區(qū)域的分叉部位病變長度10 mm、15 mm的模型高TAWSS區(qū)域要大于無病變者。

Fig.6? TAWSS distribution nephogram （from left to right： no lesion （a）， 10 mm lesion length model （b），15 mm lesion length model （c））

3.3 ?病變長度對OSI的影響

Zhang等人的那研究證實低的OSI能夠降低血管內(nèi)膜增生的可能性。由圖7可見，無病變模型、10 mm病變長度模型、15 mm病變長度模型OSI分布差別較小。

3.4 ?病變長度對TAWSSG的影響

由圖8可見，各模型的分叉病變高TWSSG區(qū)域按無病變、10 mm、15 mm順序依次遞增。在管壁分叉處，TAWSSG值都較高，有血管損傷的風(fēng)險。

3.5 ?病變長度對血流速度分布的影響

為考察病變長度對血流速度的影響，選取血管病變位置的截面速度進行比較研究，圖9為截面的位置。三個不同病變長度的截面血流速度分布如圖10所示，由圖可見，在脈動流周期的入口血流速度峰值時刻，無病變、10 mm病變長度、15 mm病變長度模型的截面血流速度依次遞增。

4 ?討論

本研究構(gòu)建了冠狀動脈的LAD的不同長度狹窄區(qū)域的血管3D模型，并在此基礎(chǔ)上使用ANSYS CFX流體仿真軟件進行仿真計算，獲得冠狀動脈不同狹窄區(qū)域長度的血流動力學(xué)指標(biāo)分布特征。冠狀動脈WSS是引起血管發(fā)生病變及病變惡化的重要影響因素，WSS高、低區(qū)域的震蕩是引起血管損傷的重要因素，其低WSS區(qū)域或震蕩可能擴大病變狹窄程度，而高WSS可能會造成斑塊的軟化和不穩(wěn)定;高OSI可能會增加內(nèi)皮細胞的功能紊亂和血管內(nèi)膜增生^[12-13]。

參考文獻

1. Mozaffarian D， Benjamin EJ， Go AS， et al. Heart disease and stroke statistics-2015 update： a report from the american heart association. Circulation， 2015， 131： e29-e322.

2. Zarins CK， Taylor CA， Min JK. Computed fractional flow reserve（FFTCT） derived from coronary CT angiography. Journal of Cardiovascular Translational Research， 2013， 6（5）： 708-714.

3. Zhang JM， Zhong L， Luo T， et al. Numerical simulation and clinical implications of stenosis in coronary blood flow. In： BioMed Research International， 2014： 1-10.

4. Taylor CA， Fonte TA， Min JK. Computational fluid dyna?mics applied to cardiac computed tomography for noninva?sive quantification of fractional flow reserve. Journal of the American College of Cardiology， 2013， 61（22）： 2233-2241.

5. Koo BK， Erglis A， Doh JH， et al. Diagnosis of ischemia- causing coronary stenoses by noninvasive fractional flow reserve computed from coronary computed tomographic angiograms. Journal of the American College of Cardiology， 2011， 58（19）： 1889-1996.

6. Yong AS， Ng AC， Brieger D， et al. Three-dimensional and two dimensional quantitative coronary angiography， and their prediction of reduced fractional flow reserve. European Heart Journal， 2011， 32（3）： 345-353.

7. Iguchi T， Hasegawa T， Nishimura S， et al. Impact of lesion length on functional significance in intermediate coronary lesions. Clinical Investigations， 2013， 36（3）： 172-177.

8. Kristensen TS， Engstr?mb T， Kelb?k H， et al. Correlation between coronary computed tomographic angiography and fractional flow reserve. International Journal of Cardiology， 2010， 144（2）： 200-205.

9. Alghamdi A， Balgaith M， Alkhaldi A. Influence of the length of coronary artery lesions on fractional flow reserve across intermediate coronary obstruction. European Heart Journal Supplements， 2014， 16（Supplement B）： 76-79.

10. 隋國慶， 張培新， 楊國柱， 等. 冠狀動脈的數(shù)值模擬分析及在支架介入的應(yīng)用研究[J]. 軟件. 2020（1）.

11. 陳廣新， 趙東良， 郭金興， 等. 基于CTA的個體化腦動脈瘤的流固耦合分析及其臨床應(yīng)用[J]. 2020（1）.

12. 王汝良， 胡霖霖， 郭金興， 等. 頸動脈分叉的非穩(wěn)態(tài)數(shù)值模擬分析[J]. 2018（10）.

13. 張凱旋， 陳廣新， 邱收， 等. 椎動脈阻斷術(shù)前后基底動脈瘤的血流動力學(xué)數(shù)值模擬分析[J]. 2019（6）.