張雪蓮，黃 偉

（1. 牡丹江醫(yī)學(xué)院醫(yī)學(xué)影像學(xué)院，黑龍江 牡丹江 157011；2. 牡丹江醫(yī)學(xué)院教務(wù)處，黑龍江 牡丹江 157011）

0 引言

椎動(dòng)脈開(kāi)窗屬于罕見(jiàn)的孔型腦動(dòng)脈，即開(kāi)窗畸形變異[1-2]。在影像學(xué)上腦動(dòng)脈開(kāi)窗畸形的發(fā)現(xiàn)率為0.3%-0.9%，常見(jiàn)于臨近椎基底動(dòng)脈交界處。椎動(dòng)脈開(kāi)窗常伴有腦動(dòng)脈瘤發(fā)生[3-4]。臨床上常根據(jù)腦動(dòng)脈瘤的形成推斷動(dòng)脈存在開(kāi)窗畸形。應(yīng)用計(jì)算流體力學(xué)方法（computational fluid dynamics，CFD）數(shù)值模擬動(dòng)脈的血流動(dòng)力學(xué)指標(biāo)，可對(duì)椎動(dòng)脈開(kāi)窗形態(tài)的血流狀態(tài)進(jìn)行定量分析[6-8]，且可對(duì)動(dòng)脈瘤的形成的危險(xiǎn)因素進(jìn)行預(yù)測(cè)、分析，對(duì)臨床治療進(jìn)行指導(dǎo)。本文旨在通過(guò)對(duì)一例椎動(dòng)脈開(kāi)窗畸形影像數(shù)據(jù)進(jìn)行三維重建，應(yīng)用計(jì)算機(jī)數(shù)值模擬仿真，分析開(kāi)窗畸形腦動(dòng)脈的血流動(dòng)力學(xué)指標(biāo)及對(duì)動(dòng)脈瘤的形成影響[9-12]，對(duì)臨床診斷與治療預(yù)后提供借鑒和指導(dǎo)。

1 材料與方法

1.1 實(shí)驗(yàn)設(shè)備與軟件

本實(shí)驗(yàn)分別采用醫(yī)學(xué)交互式影像控制系統(tǒng)、CFD網(wǎng)格劃分軟件FLUENT MESHING、CFD仿真計(jì)算軟件CFX、圖像后處理軟件ENSIGHT進(jìn)行建模、網(wǎng)格劃分、計(jì)算、結(jié)果分析（表 1）。應(yīng)用TOSHIBA/AQUILION 64排 CT掃描，選取電流250 mA，電壓120 kv，層厚0.5 mm；前臂肘靜脈團(tuán)注射造影劑，注射速度3.5 ml/s，總量為90 ml；CT數(shù)據(jù)以 DICOM3.0 格式存儲(chǔ)。

表1 實(shí)驗(yàn)使用軟件Tab.1 Experimental software

1.2 三維模型重建與網(wǎng)格劃分

由牡丹江醫(yī)學(xué)院附屬紅旗醫(yī)院提供的椎動(dòng)脈開(kāi)窗畸形患者的 CTA影像數(shù)據(jù)，患者女性，68歲，突發(fā)頭痛，CTA顯示椎動(dòng)脈開(kāi)窗畸形。將采集的CTA影像數(shù)據(jù)以 DICOM格式導(dǎo)入MIMICS21.0軟件中進(jìn)行圖像處理，采用閾值分割、計(jì)算三維等算法，生成初步的三維模型，然后導(dǎo)入3-matic軟件進(jìn)行表面網(wǎng)格光順、除去細(xì)小分支，并經(jīng)多次光滑、重畫網(wǎng)格，切平出入口（2進(jìn)4出），最后檢查三角面片質(zhì)量，最后將編輯完畢的椎動(dòng)脈模型以STL格式導(dǎo)出保存。使用ANSYS FLUENT MESHING軟件對(duì)構(gòu)建好的STL格式三維模型進(jìn)行網(wǎng)格劃分，由于模型結(jié)構(gòu)比較復(fù)雜，按照非結(jié)構(gòu)化的四面體進(jìn)行網(wǎng)格劃分，為保證計(jì)算精度，邊界層采用5層加密。

1.3 計(jì)算設(shè)置

1.3.1 邊界設(shè)置

入口采用速度入口，為保證計(jì)算快速收斂，兩個(gè)入口速度值存在較小的差別。入口速度曲線如圖1-2所示。出口給予脈動(dòng)壓力，如圖3所示。出口、入口位置如圖4所示。

圖1 入口1速度曲線Fig.1 Velocity of inlet1

圖2 入口2速度曲線Fig.2 Velocity of inlet2

圖3 出口壓力曲線Fig.3 Pressure of outlet

圖4 椎動(dòng)脈入口、出口位置圖Fig.4 Location of vertebral artery inlet and outlet

1.3.2 求解設(shè)置

根據(jù)雷諾數(shù)計(jì)算設(shè)置采用層流，牛頓、不可壓縮流體，不考慮重力影響。本計(jì)算采用瞬態(tài)計(jì)算，計(jì)算三個(gè)周期。步長(zhǎng)設(shè)為0.0008 s，總步數(shù)為1000步。為保證計(jì)算收斂穩(wěn)定，取計(jì)算結(jié)果第三個(gè)周期進(jìn)行分析。

1.4 觀察指標(biāo)

血流動(dòng)力學(xué)參數(shù)選取：時(shí)間平均壁面切應(yīng)力（time average wall shear stress，TAWSS）、平均壁面切應(yīng)力梯度（time average wall shear stress grade，TAWSSG）、剪切震蕩系數(shù)（oscillatory shear index，OSI）、動(dòng)脈瘤形成指標(biāo)（aneurysm formation index，AFI）），梯度震蕩數(shù)值（GON）等參數(shù)作為觀察指標(biāo)。由于 CFX無(wú)法直接輸出 TAWSS、TAWSSG、OSI、AFI、GON等指標(biāo)，需要ccl語(yǔ)言編程。

2 結(jié)果

2.1 OSI云圖

由OSI云圖（圖5）可見(jiàn)在椎動(dòng)脈開(kāi)窗處，存在高OSI區(qū)域，其他區(qū)域表現(xiàn)低OSI區(qū)域。表明由于椎動(dòng)脈存在開(kāi)窗畸形導(dǎo)致不同的椎動(dòng)脈區(qū)域血流震蕩，高低OSI震蕩易導(dǎo)致動(dòng)脈瘤的發(fā)生。

圖5 椎動(dòng)脈開(kāi)窗OSI云圖Fig.5 Vertebral artery fenestration OSI cloud

2.2 TAWSS云圖

圖6 為椎動(dòng)脈TAWSS云圖，TAWSS為脈動(dòng)流在一個(gè)心臟周期內(nèi)用每個(gè)節(jié)點(diǎn)上積分 WSS量的值來(lái)計(jì)算 TAWSS，TAWSS云圖中椎動(dòng)脈各區(qū)域中呈現(xiàn)出分布不均，TAWSS呈現(xiàn)紅色，表明TAWSS較大，TAWSS呈現(xiàn)藍(lán)色表明 TAWSS較小，圖中高TAWSS區(qū)域分布于椎動(dòng)脈的出口附近區(qū)域：主干、分叉。由于主干區(qū)域的血流速度較大，該區(qū)域TAWSS較其他區(qū)域集中，數(shù)值較大。

圖6 TAWSS云圖Fig.6 TAWSS Distribution Nephogram

2.3 椎動(dòng)脈GON云圖

圖7 為GON云圖，由圖可見(jiàn)，椎動(dòng)脈GON分布極其不均勻，在椎動(dòng)脈開(kāi)窗附近的主血管區(qū)域呈現(xiàn)較大的GON值，表明椎動(dòng)脈表面存在強(qiáng)烈的震蕩力和壓縮力，原因是血液在進(jìn)入瘤腔后形成渦流，導(dǎo)致動(dòng)脈壁震蕩，這種沖擊對(duì)血管壁造成膨脹或者擴(kuò)張。

圖7 GON云圖Fig.7 GON distribution nephogram

2.4 椎動(dòng)脈AFI云圖

圖8 中椎動(dòng)脈AFI值普遍較大，在椎動(dòng)脈的分叉處，存在相對(duì)較小的AFI值，說(shuō)明血流不斷沖擊管壁，形成渦流，血流運(yùn)動(dòng)極其不穩(wěn)定，容易導(dǎo)致腦動(dòng)脈瘤的形成。

圖8 AFI云圖Fig.8 AFI distribution nephogram

2.5 TAWSSG云圖

圖9 為椎動(dòng)脈開(kāi)窗畸形的TAWSSG云圖，由圖中可見(jiàn)，高 TAWSSG均分布于分叉、開(kāi)窗區(qū)域，TAWSSG較大易造成血管壁損傷。椎動(dòng)脈開(kāi)窗畸形導(dǎo)致血流在狹窄、分叉區(qū)域速度較快，沖擊血管壁，導(dǎo)致血管壁局部區(qū)域TAWSSG較高。

圖9 TAWSSG云圖Fig.9 TAWSSG distribution nephogram

3 結(jié)論

椎動(dòng)脈開(kāi)窗畸形引起血液異常流動(dòng)，在此基礎(chǔ)上，AFI表征動(dòng)脈瘤形成指數(shù)[13-15]，椎動(dòng)脈開(kāi)窗畸形血管形成動(dòng)脈瘤的風(fēng)險(xiǎn)較大。同時(shí)，OSI震蕩表明，血液在流動(dòng)過(guò)程中，局部區(qū)域速度忽高忽低，來(lái)回震蕩，形成動(dòng)脈瘤易發(fā)因素。此外，在其他血流動(dòng)力學(xué)指標(biāo)上，椎動(dòng)脈開(kāi)窗畸形也導(dǎo)致了異常情況的發(fā)生，發(fā)展。通過(guò)計(jì)算機(jī)三維建模與計(jì)算流體力學(xué)仿真，可獲取椎動(dòng)脈開(kāi)窗畸形的內(nèi)部血液流變指標(biāo)，進(jìn)而可為臨床應(yīng)用提供指導(dǎo)與借鑒。

4 討論

椎動(dòng)脈瘤與椎動(dòng)脈開(kāi)窗畸形有一定的因果關(guān)系，一些學(xué)者[16-20]已通過(guò)宏觀血流動(dòng)力學(xué)特征證實(shí)，椎動(dòng)脈開(kāi)窗畸形動(dòng)脈瘤與兩支血管之間形成囊動(dòng)脈瘤[21-28]是由于血流動(dòng)力學(xué)參數(shù)劇烈變化的導(dǎo)致的結(jié)構(gòu)改變，血液流動(dòng)時(shí)，椎動(dòng)脈左血管管徑較粗，形成優(yōu)勢(shì)支，其血流量、流速均大于右支[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

END

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

END

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

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

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

END

TRANSIENT SCHEME：

Option = Second Order Backward Euler

TIMESTEP INITIALISATION：

Option = Extrapolation

END

END

END

EXPERT PARAMETERS：

linearly exact numerics = t

END

END

INTERPOLATOR STEP CONTROL：

Runtime Priority = Standard

MEMORY CONTROL：

Memory Allocation Factor = 1.0

END

END

PARALLEL HOST LIBRARY：

HOST DEFINITION： desktop0jjvpii

Remote Host Name = DESKTOP-0JJVPII

Host Architecture String = winntamd64

Installation Root = D：Program FilesANSYS Incv%vCFX

END

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

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

END

END

END