椎动脉开窗血流动力学仿真分析

2020-07-10张雪莲

软件 2020年5期

张雪莲，黄伟

（1. 牡丹江医学院医学影像学院，黑龙江牡丹江 157011；2. 牡丹江医学院教务处，黑龙江牡丹江 157011）

0 引言

椎动脉开窗属于罕见的孔型脑动脉，即开窗畸形变异[1-2]。在影像学上脑动脉开窗畸形的发现率为0.3%-0.9%，常见于临近椎基底动脉交界处。椎动脉开窗常伴有脑动脉瘤发生[3-4]。临床上常根据脑动脉瘤的形成推断动脉存在开窗畸形。应用计算流体力学方法（computational fluid dynamics，CFD）数值模拟动脉的血流动力学指标，可对椎动脉开窗形态的血流状态进行定量分析[6-8]，且可对动脉瘤的形成的危险因素进行预测、分析，对临床治疗进行指导。本文旨在通过对一例椎动脉开窗畸形影像数据进行三维重建，应用计算机数值模拟仿真，分析开窗畸形脑动脉的血流动力学指标及对动脉瘤的形成影响[9-12]，对临床诊断与治疗预后提供借鉴和指导。

1 材料与方法

1.1 实验设备与软件

本实验分别采用医学交互式影像控制系统、CFD网格划分软件FLUENT MESHING、CFD仿真计算软件CFX、图像后处理软件ENSIGHT进行建模、网格划分、计算、结果分析（表 1）。应用TOSHIBA/AQUILION 64排 CT扫描，选取电流250 mA，电压120 kv，层厚0.5 mm；前臂肘静脉团注射造影剂，注射速度3.5 ml/s，总量为90 ml；CT数据以 DICOM3.0 格式存储。

表1 实验使用软件Tab.1 Experimental software

1.2 三维模型重建与网格划分

由牡丹江医学院附属红旗医院提供的椎动脉开窗畸形患者的 CTA影像数据，患者女性，68岁，突发头痛，CTA显示椎动脉开窗畸形。将采集的CTA影像数据以 DICOM格式导入MIMICS21.0软件中进行图像处理，采用阈值分割、计算三维等算法，生成初步的三维模型，然后导入3-matic软件进行表面网格光顺、除去细小分支，并经多次光滑、重画网格，切平出入口（2进4出），最后检查三角面片质量，最后将编辑完毕的椎动脉模型以STL格式导出保存。使用ANSYS FLUENT MESHING软件对构建好的STL格式三维模型进行网格划分，由于模型结构比较复杂，按照非结构化的四面体进行网格划分，为保证计算精度，边界层采用5层加密。

1.3 计算设置

1.3.1 边界设置

入口采用速度入口，为保证计算快速收敛，两个入口速度值存在较小的差别。入口速度曲线如图1-2所示。出口给予脉动压力，如图3所示。出口、入口位置如图4所示。

图1 入口1速度曲线Fig.1 Velocity of inlet1

图2 入口2速度曲线Fig.2 Velocity of inlet2

图3 出口压力曲线Fig.3 Pressure of outlet

图4 椎动脉入口、出口位置图Fig.4 Location of vertebral artery inlet and outlet

1.3.2 求解设置

根据雷诺数计算设置采用层流，牛顿、不可压缩流体，不考虑重力影响。本计算采用瞬态计算，计算三个周期。步长设为0.0008 s，总步数为1000步。为保证计算收敛稳定，取计算结果第三个周期进行分析。

1.4 观察指标

血流动力学参数选取：时间平均壁面切应力（time average wall shear stress，TAWSS）、平均壁面切应力梯度（time average wall shear stress grade，TAWSSG）、剪切震荡系数（oscillatory shear index，OSI）、动脉瘤形成指标（aneurysm formation index，AFI）），梯度震荡数值（GON）等参数作为观察指标。由于 CFX无法直接输出 TAWSS、TAWSSG、OSI、AFI、GON等指标，需要ccl语言编程。

2 结果

2.1 OSI云图

由OSI云图（图5）可见在椎动脉开窗处，存在高OSI区域，其他区域表现低OSI区域。表明由于椎动脉存在开窗畸形导致不同的椎动脉区域血流震荡，高低OSI震荡易导致动脉瘤的发生。

图5 椎动脉开窗OSI云图Fig.5 Vertebral artery fenestration OSI cloud

2.2 TAWSS云图

图6 为椎动脉TAWSS云图，TAWSS为脉动流在一个心脏周期内用每个节点上积分 WSS量的值来计算 TAWSS，TAWSS云图中椎动脉各区域中呈现出分布不均，TAWSS呈现红色，表明TAWSS较大，TAWSS呈现蓝色表明 TAWSS较小，图中高TAWSS区域分布于椎动脉的出口附近区域：主干、分叉。由于主干区域的血流速度较大，该区域TAWSS较其他区域集中，数值较大。

图6 TAWSS云图Fig.6 TAWSS Distribution Nephogram

2.3 椎动脉GON云图

图7 为GON云图，由图可见，椎动脉GON分布极其不均匀，在椎动脉开窗附近的主血管区域呈现较大的GON值，表明椎动脉表面存在强烈的震荡力和压缩力，原因是血液在进入瘤腔后形成涡流，导致动脉壁震荡，这种冲击对血管壁造成膨胀或者扩张。

图7 GON云图Fig.7 GON distribution nephogram

2.4 椎动脉AFI云图

图8 中椎动脉AFI值普遍较大，在椎动脉的分叉处，存在相对较小的AFI值，说明血流不断冲击管壁，形成涡流，血流运动极其不稳定，容易导致脑动脉瘤的形成。

图8 AFI云图Fig.8 AFI distribution nephogram

2.5 TAWSSG云图

图9 为椎动脉开窗畸形的TAWSSG云图，由图中可见，高 TAWSSG均分布于分叉、开窗区域，TAWSSG较大易造成血管壁损伤。椎动脉开窗畸形导致血流在狭窄、分叉区域速度较快，冲击血管壁，导致血管壁局部区域TAWSSG较高。

图9 TAWSSG云图Fig.9 TAWSSG distribution nephogram

3 结论

椎动脉开窗畸形引起血液异常流动，在此基础上，AFI表征动脉瘤形成指数[13-15]，椎动脉开窗畸形血管形成动脉瘤的风险较大。同时，OSI震荡表明，血液在流动过程中，局部区域速度忽高忽低，来回震荡，形成动脉瘤易发因素。此外，在其他血流动力学指标上，椎动脉开窗畸形也导致了异常情况的发生，发展。通过计算机三维建模与计算流体力学仿真，可获取椎动脉开窗畸形的内部血液流变指标，进而可为临床应用提供指导与借鉴。

4 讨论

椎动脉瘤与椎动脉开窗畸形有一定的因果关系，一些学者[16-20]已通过宏观血流动力学特征证实，椎动脉开窗畸形动脉瘤与两支血管之间形成囊动脉瘤[21-28]是由于血流动力学参数剧烈变化的导致的结构改变，血液流动时，椎动脉左血管管径较粗，形成优势支，其血流量、流速均大于右支[29-33]在血管分叉处，形成最大压力。

5 ccl部分源代码

ADDITIONAL VARIABLE：

ADDITIONAL VARIABLE： Gmag

Additional Variable Value = sqrt(Gx *Gx + Gy * Gy + Gz * Gz)

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： Gx

Additional Variable Value = -(1.0-Normal X*Normal X)*WSSxF.Gradient

X -(0.0-Normal X*Normal Y)*WSSxF.Gradient Y -(0.0-Normal X*Normal

Z)*WSSxF.Gradient Z

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： Gy

Additional Variable Value = -(0.0-Normal Y*Normal X)*WSSyF.Gradient

X -(1.0-Normal Y*Normal Y)*WSSyF.Gradient Y -(0.0-Normal Y*Normal

Z)*WSSyF.Gradient Z

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： Gz

Additional Variable Value = -(0.0-Normal Z*Normal X)*WSSzF.Gradient

X -(0.0-Normal Z*Normal Y)*WSSzF.Gradient Y -(1.0-Normal Z*Normal

Z)*WSSzF.Gradient Z

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： PressGauge

Additional Variable Value = pref +Pressure

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： Qvar

Additional Variable Value = -0.5*(0.5*Shear Strain Rate^2 -

0.5 *vortmag2)

Option = Algebraic Equation

END

ADDITIONAL VARIABLE： ViscDisp

Additional Variable Value = (Shear Strain Rate)^2

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)*WSSField.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)*WSSField.Gradient Y -(0.0-Normal Y*Normal

Z)*WSSField.Gradient Z

Vector zValue = -(0.0-Normal Z*Normal X)*WSSField.Gradient X

-(0.0-Normal Z*Normal Y)*WSSField.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

INITIALISATION：

Option = Automatic

INITIAL CONDITIONS：

Velocity Type = Cartesian

ADDITIONAL VARIABLE： WSSField

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

Option = Automatic with Value

END

ADDITIONAL VARIABLE： WSSxF

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

Option = Automatic with Value

END

ADDITIONAL VARIABLE： WSSyF

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

Option = Automatic with Value

END

ADDITIONAL VARIABLE： WSSzF

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

Option = Automatic with Value

END

CARTESIAN VELOCITY COMPONENTS：

Option = Automatic with Value

U = 0 [m s^-1]

V = 0 [m s^-1]

W = 0 [m s^-1]

END

STATIC PRESSURE：

Option = Automatic with Value

Relative Pressure = 0 [Pa]

END

OUTPUT CONTROL：

RESULTS：

File Compression Level = Default

Option = Standard

END

TRANSIENT RESULTS： Transient Results 1

File Compression Level = Default

Include Mesh = No

Option = Selected Variables

Output Variables List = Absolute Pressure,Total Pressure,Pressure,PressGauge,Velocity,Velocity u,Velocity v,Velocity

w,Vorticity,Vorticity X,Vorticity Y,Vorticity Z,Wall Shear,Wall

Shear X,Wall Shear Y,Wall Shear Z,Shear Strain

Rate,Qvar,ViscDisp,WSSField,WSSG,AFI

OUTPUT FREQUENCY：

Option = Time List

Time List = 1.6 [s], 1.608 [s], 1.616 [s],1.624 [s], 1.632 [s],

1.64 [s], 1.648 [s], 1.656 [s], 1.664[s], 1.672 [s], 1.68 [s],

1.688 [s], 1.696 [s], 1.704 [s], 1.712[s], 1.72 [s], 1.728 [s],

1.736 [s], 1.744 [s], 1.752 [s], 1.76[s], 1.768 [s], 1.776 [s],

1.784 [s], 1.792 [s], 1.8 [s], 1.808 [s],1.816 [s], 1.824 [s],

1.832 [s], 1.84 [s], 1.848 [s], 1.856[s], 1.864 [s], 1.872 [s],

1.88 [s], 1.888 [s], 1.896 [s], 1.904[s], 1.912 [s], 1.92 [s],

1.928 [s], 1.936 [s], 1.944 [s], 1.952[s], 1.96 [s], 1.968 [s],

1.976 [s], 1.984 [s], 1.992 [s], 2 [s],2.008 [s], 2.016 [s],

2.024 [s], 2.032 [s], 2.04 [s], 2.048[s], 2.056 [s], 2.064 [s],

2.072 [s], 2.08 [s], 2.088 [s], 2.096[s], 2.104 [s], 2.112 [s],

2.12 [s], 2.128 [s], 2.136 [s], 2.144[s], 2.152 [s], 2.16 [s],

2.168 [s], 2.176 [s], 2.184 [s], 2.192[s], 2.2 [s], 2.208 [s],

2.216 [s], 2.224 [s], 2.232 [s], 2.24[s], 2.248 [s], 2.256 [s],

2.264 [s], 2.272 [s], 2.28 [s], 2.288[s], 2.296 [s], 2.304 [s],

2.312 [s], 2.32 [s], 2.328 [s], 2.336[s], 2.344 [s], 2.352 [s],

2.36 [s], 2.368 [s], 2.376 [s], 2.384[s], 2.392 [s], 2.4 [s]

END

TRANSIENT STATISTICS： Trans2

Option = Maximum

Output Variables List = OSIfield,WSSField,WSSG

Start Iteration List = 201

Stop Iteration List = 300

END

TRANSIENT STATISTICS： Transient-StatsMin

Option = Minimum

Output Variables List = AFI

Start Iteration List = 201

Stop Iteration List = 300

END

TRANSIENT STATISTICS： trnstat1

Option = Arithmetic Average

Output Variables List = Velocity,Wall Shear,WSSField,AFI,Gx,Gy,Gz,Gmag

Start Iteration List = 201

Stop Iteration List = 300

END

SOLVER CONTROL：

ADVECTION SCHEME：

Blend Factor = 1.0

Option = Specified Blend Factor

END

CONVERGENCE CONTROL：

Maximum Number of Coefficient Loops= 5

Minimum Number of Coefficient Loops = 1

Timescale Control = Coefficient Loops

END

CONVERGENCE CRITERIA：

Residual Target = 0.005

Residual Type = MAX

END

EQUATION CLASS： av

ADVECTION SCHEME：

Option = High Resolution

END

TRANSIENT SCHEME：

Option = Second Order Backward Euler

TIMESTEP INITIALISATION：

Option = Extrapolation

END

EXPERT PARAMETERS：

linearly exact numerics = t

END

INTERPOLATOR STEP CONTROL：

Runtime Priority = Standard

MEMORY CONTROL：

Memory Allocation Factor = 1.0

END

PARALLEL HOST LIBRARY：

HOST DEFINITION： desktop0jjvpii

Remote Host Name = DESKTOP-0JJVPII

Host Architecture String = winntamd64

Installation Root = D：Program FilesANSYS Incv%vCFX

END

PARTITIONER STEP CONTROL：

Multidomain Option = Automatic

Runtime Priority = Standard

MEMORY CONTROL：

Memory Allocation Factor = 1.0

END

PARTITION SMOOTHING：

Maximum Partition Smoothing Sweeps= 100

Option = Smooth

END

PARTITIONING TYPE：

MeTiS Type = k-way

Option = MeTiS

Partition Size Rule = Automatic

END

RUN DEFINITION：

Run Mode = Full

SOLVER STEP CONTROL：

Runtime Priority = Standard

MEMORY CONTROL：

Memory Allocation Factor = 1.0

END

PARALLEL ENVIRONMENT：

Number of Processes = 1

Start Method = Serial

END

猜你喜欢