不同边界条件对主动脉血流仿真结果的影响
曹欣荣1, 王婕1, 王荣品2, 张先文1, 唐劲天1
1. 清华大学 工程物理系, 粒子技术与辐射成像教育部重点实验室, 北京 100084
2. 贵州省人民医院 放射科, 贵阳 550002
唐劲天, 研究员, E-mail:tangjt@mail.tsinghua.edu.cn

作者简介: 曹欣荣(1988—), 男(汉), 山东, 博士研究生。

摘要

主动脉血流仿真可应用于血流运动和部分心血管疾病形成原因分析和风险评价。该文基于核磁共振成像扫描重建的主动脉模型,通过磁共振相位对比成像测量主动脉及主要支脉的时间-流量曲线及血液动力学参数,讨论了不同边界条件下时间-流量曲线、时间-壁面切应力曲线、平均壁面切应力和震荡剪切因子的影响。结果表明: 出口边界相对压力为0、 采用拟合血流方程的入口条件等边界下计算结果与测量的实际血流曲线作为出口和入口边界条件计算结果差异较大,其中入口参数对壁面切应力、平均壁面切应力和震荡剪切因子计算结果影响较大,出口参数对出口时间-流量曲线影响较大。

关键词: 计算流体力学; 主动脉; 磁共振相位对比成像; 平均壁面切应力; 震荡剪切因子
中图分类号:Q66 文献标志码:A 文章编号:1000-0054(2014)06-0700-06
Influence of the boundary conditions on aorta blood flow simulations
Xinrong CAO1, Jie WANG1, Rongpin WANG2, Xianwen ZHANG1, Jintian TANG1
1. Key Laboratory of Particle and Radiation Imaging of Ministry of Education, Department of Engineering Physics, Tsinghua University, Beijing 100084, China
2. Department of Radiology, Guizhou Provincial People’s Hospital, Guiyang 550002, China
Abstract

Simulations of blood flows in the aorta area provide a better understanding of the hemodynamic blood flow mechanism for evaluating the risk of some cardiovascular diseases. An aorta model was reconstructed from magnetic resonance imaging (MRI) images while the blood flow and hemodynamic parameters in the ascending aorta, the descending aorta and the main peripheral arteries were measured using phase contrast MRI (PC-MRI). The boundary conditions affected the flow rate, wall shear stress, averaged wall shear stress and the oscillating shear index. The results show that simulations with zero pressure outlets, simulations with fitted input flows and zero pressure outlets, and simulations with fitted input flow and outflow discharges give quite different results. The input flow variations are mainly related to the local and average wall shear stresses and the oscillating shear index. The outlet flows are strongly related to the outlet boundary conditions.

Keyword: computational fluid dynamics; aorta; phase contrast magnetic resonance imaging (PC-MRI); averaged wall shear stress; oscillating shear index

根据世界卫生组织统计,心血管疾病已经成为仅次于癌症的威胁人类健康的第二杀手。通过计算流体力学(computational fluid dynamics, CFD)进行血流流动仿真应用,了解血液流动过程中与管壁的相互作用是探究心血管疾病的有效方式。胸主动脉处血管结构复杂,了解其血液流动对部分心血管疾病(动脉粥样硬化病变、动脉瘤和主动脉狭窄等)的形成原因和潜在危险评价等有重要作用[1,2]

目前,主动脉模型的计算主要集中于血液流动情况和壁面切应力(wall shear stress, WSS)的分析。主要分析参数有平均壁面切应力(averaged wall shear stress, AWSS)和震荡剪切因子(oscillatory shear index, OSI)等[3,4]。研究表明, AWSS和OSI在血管瘤风险预测[1]和动脉粥样硬化诊疗[5]等方面有着重要意义。

国内一部分学者利用几何构造的模型对主动脉弓处的血液流动进行了仿真[6], 一部分学者利用影像重建的主动脉模型进行了仿真计算[7]。但是,大部分采用的血流模型入口出口参数与实际中血液流动存在一定的差异: 部分学者采用拟合的血流方程模拟入口血流,在出口处部分采用出口相对压力为0的边界条件进行仿真计算[8]; 部分学者采用流量分数比例进行计算[9]。目前,在利用实际测量血流流量或者时间-流速/流量曲线进行仿真计算这方面研究较少。

本文主要通过磁共振相位对比成像(phase contrast magnetic resonance imaging, PC-MRI)技术测量了主动脉和其主要支脉的时间-流量曲线,利用3维影像重建获取了主动脉模型,利用流体仿真软件Fluent分析以下影响因素: 1) 测量血流流量作为出口入口条件; 2) 测量血流流量为入口条件,出口相对压力为0; 3) 拟合血流流量作为入口边界条件,相对压力为0作为出口边界条件; 4) 拟合血流流量作为入口边界条件,出口采用相对流量边界条件。考察各种边界条件对模型出口血流速度和壁面切应力以及AWSS和OSI的影响。

1 材料与方法

选择一位健康成年男性志愿者,通过MRI进行断层扫描重建其主动脉,利用PC-MRI测量其在升主动脉、降主动脉、左右锁骨下动脉和左右颈总动脉的时间-流量曲线和血液动力学参数,通过计算流体力学软件对重建的主动脉模型进行仿真。

1.1 模型重建

通过MRI对该志愿者进行冠状位扫描胸部影像(Siemens, SimphonyTim, 30层,间距2 mm)。获取影像后,通过区域增长方法选择主动脉以及相应支脉的血液区域进行重建,对重建模型平滑处理,最终获取的主动脉模型见图1。模型的血流入口为升主动脉(asending aorta, AA), 出口包括降主动脉(desending aorta, DA)末端、左右锁骨下动脉(left and right subclavian, LS和RS)末端、左右颈总动脉(left and right carotid, LC和RC)末端,与志愿者实际测量部位一致。

图1 重建的主动脉模型 (图中点1、 2、 3、 4为选择计算AWSS和OSI的兴趣点)

1.2 血流测量与模型

PC-MRI是一种无创的血流测量方法,有标准的定位方法和扫描参数,其测量结果不受骨骼及其他周边组织的影响[10]。本研究中主要测量了升主动脉、降主动脉、左右锁骨下动脉和左右颈总动脉的时间-流量曲线和血液动力学参数。

具体测量如下: 升主动脉选择距主动脉瓣膜2 cm处,降主动脉选择隔肌平面,左右锁骨下动脉选择距开口2~3 cm处,左右颈总动脉选择距开口2~3 cm处,采用电影法快速相位对比序列屏气下扫描。扫描参数: 脉冲序列重复时间/回波时间(TR/TE)=自动选择最小重复时间/Min full, 翻转角20°, 带宽31.25 kHz, 视野40 cm×40 cm, 矩阵256×128, 扫描层厚5 mm, 激励次数1, 速度编码值150 cm/s, 编码方向为Slice。采用回顾性心电门控触发分段K空间成像,依次对升主动脉、降主动脉、左右锁骨下动脉及左右颈总动脉垂直于血管长轴方向进行单层20个时相屏气扫描。使用跨平面技术测量,测量部位垂直于血流方向。通过Siemens argus血流测量软件进行后处理,得到一个心动周期的时间-流量曲线(图2)[11]。为了满足仿真计算步长要求,对获取的时间 -流量曲线通过3次样条差值获取测量时间外的血流流量值,差值后相邻数据点时间间隔为1 ms。

图2 Siemens argus血流测量软件测量获取的升主动脉时间流量图

1.3 仿真计算

将重建的主动脉模型导入ICEM(ANSYS Inc.)进行表面分割,分为血流的入口和5个出口以及血管壁。网格划分时,血管壁采用角柱网格,内部采用四面体网格。网格数为244 053, 节点数为66 881。

仿真模型主要根据Navier-Stokes方程求解,

ρut+ρ(u·)u=·[-pI+τ]+F.(1)

其中: ρ表示血液密度, u为血流速度矢量, p为压力, I为单位张量, F为体积力,

τ=μ((u+u)T)).(2)

计算中血液密度选择为1 060 kg/m3[12], 黏度系数μ为0.003 5 kg/(m·s)[13]

本文中采用Ansys Workbench中的Fluent模块通过有限体积法进行计算。利用Fluent的用户定义函数(user defined function, UDF)将时变的入口或者出口血流速度作为相应的边界条件。计算时间间隔为0.001 s, 计算中每10步保存一次计算结果,心动周期为0.8 s, 计算长度为4.0 s, 选取计算结果稳定的第5个周期计算结果作为最终分析结果。

为了对比不同的血流参数对计算结果的影响,将测量的各动脉处的血流流量结果作为出口入口边界条件,计算结果记为Sim1; 入口血流采用测量数据,出口采用相对压力为0的边界条件,计算结果记为Sim2; 入口血流采用文[6]中的血流速度作为边界条件,各个出口相对压力为0, 其余计算参数不变,记为Sim3; 将各个模型的流量比值[9]作为出口参数,入口血流采用Sim3相同的入口条件进行计算,记为Sim4。

Vinlet(t)=2.40e-7.557tsin13.09t,0t0.24s;-0.344e-7.557(t-0.24)sin13.09(t-0.24),0.24s<t0.274s;0,0.274s<t0.8s.3

重建主动脉模型入口截面面积为3.83 cm2, Sim1和Sim2的入口血流平均速度为0.23 m/s, Sim3和Sim4的入口血流平均速度为0.20 m/s。Sim1和Sim2入口处平均Reynolds数为1 519, Womersley数为17; Sim3和Sim4入口处平均Reynolds数为1 316, Womersley数为17[14,15]

1.4 计算结果处理

计算Sim1—Sim4的出口、入口的流量和、沿主动脉血液流动方向依次选取了4个兴趣点(图1兴趣点1—4)的平均壁面切应力(AWSS)和震荡剪切因子值(OSI)。

流量和计算公式为

V=0TQ(t)dt.(4)

其中, Q( t)为该处的血流流量, T为心动周期时间。

兴趣点的AWSS和OSI值[3,4]主要计算公式为

AWSS=1T0Tτwdt,(5)OSI=121-0Tτwdt/0Tτwdt.(6)

2 结果分析
2.1 计算出口的血流流量对比

图3显示了Sim1—Sim4计算获取的入口和出口处的血流流量 Q随时间变化的曲线。可以看出,与Sim1相比, Sim2只有在降主动脉时间-流量曲线与测量结果接近,其他出口差异较大。Sim2和Sim3的时间-流量曲线类似。Sim4计算降主动脉和左锁骨下动脉的流量峰值最大, Sim3左右颈总动脉和右锁骨下动脉的流量峰值最大。另外, Sim3左右锁骨下动脉和左右颈总动脉射血后(0.2~0.4 s)的回流( Q<0)较大。

图3 一个心动周期内的各个出入口的血流量结果

计算获取Sim1—Sim4一个心动周期的各个入口和出口流量和,结果见表1

表1 出入口流量和结果

表1中可以看出,与Sim1相比, Sim2的降主动脉流量和较大,右锁骨下动脉和右颈总动脉流量和较小; Sim3降主动脉流量和较大,而左右锁骨下动脉和右颈总动脉结果较小; Sim4的降主动脉、左右颈总动脉流量和较小,左右所锁骨下动脉流量和较大。Sim4与Sim1的流量和结果差异表明实际测量的个体结果和统计获取的结果存在一定的差异,因此在仿真计算中采用个体测量的结果更加符合实际。

2.2 0.1 s时壁面切应力分布

图4 t=0.1 s时Sim1—Sim4的壁面切应力分布轮廓图。可以看出, Sim2和Sim3在左右锁骨下动脉和颈左右总动脉处的壁面切应力较大, Sim4降主动脉段和左锁骨下动脉壁面切应力值较大。这有可能是因为Sim2和Sim3在0.1 s时刻在左右锁骨下动脉的血流流量较大,而Sim4在降主动脉和左锁骨下动脉处的血流流量较大。

图4 t=0.1 s时刻Sim1—Sim4的壁面切应力分布

2.3 兴趣点壁面切应力时间曲线

兴趣点1—4的壁面切应力随时间变化的曲线如图5所示。可以看出, Sim1和Sim2计算的壁面切应力曲线类似,在0~0.2 s射血初期基本一致; Sim3和Sim4的壁面切应力曲线类似,兴趣点1、 2、 3在0~0.2 s射血初期结果基本一致,而兴趣点4结果差异较大。Sim1、 Sim2与Sim3、 Sim4的计算结果存在较大的差异。

图5 选取的兴趣点在一个心动周期内的时间壁面切应力曲线

根据壁面切应力定义,

τw=μuy|y=0,(7)

Sim1、 Sim2与Sim3、 Sim4模型计算的该处速度切向变化不同。Sim2和Sim3均采用压力为0的出口边界条件,说明不同的入口边界条件会引起壁面处的速度不同,进而影响壁面的壁面切应力分布。对比Sim1和Sim2的结果或者Sim3和Sim4的结果可以说明,选定的4种不同出口边界条件对壁面切应力影响较小,在计算壁面切应力时入口边界条件更为重要。兴趣点4在Sim3、 Sim4的壁面切应力结果差异较大有可能是因为降主动脉处血流流速结果差异较大。因此,在AWSS和OSI分析中,最好采用实际测量的时间-流量曲线作为边界条件。

2.4 兴趣点AWSS和OSI值

选择计算了兴趣点1和4的AWSS值和OSI值,并且计算Sim2、 Sim3和Sim4的计算值与Sim1值的相对误差(见表2和3)。相对误差为

Error=AWSSSimi-AWSSSim1AWSSSim1×100%.(8)

表2和3表明, Sim1、 Sim2和Sim4兴趣点1的AWSS和OSI计算结果接近,兴趣点4的AWSS计算结果Sim1、 Sim2和Sim3接近。这说明 AWSS 和OSI计算结果变化趋势可能不一致,同时壁面切应力曲线差异与AWSS、 OSI计算结果并不一致。例如,点1处Sim1和Sim4的AWSS 和OSI值接近,但是时间-壁面切应力曲线差异较大(图5), 说明AWSS、 OSI与壁面切应力曲线差异不直接相关。

表2 兴趣点1的AWSS和OSI值计算结果
表3 兴趣点4的AWSS和OSI值计算结果
3 讨 论

本研究结果表明在仿真计算主动脉处血流时,采用不同的入口、出口边界条件下出口时间-流量曲线、时间-壁面切应力曲线、 AWSS和OSI等存在较大差异。目前,主要计算采用的边界条件与Sim3、 Sim4边界相同,没有考虑个体差异性,因此有可能计算结果与实际人体结果存在较大的差异。

在实际血液流动中,不同支脉上游下游动脉阻力和顺应性随血流运动改变不同,血压随时间变化曲线也不相同[16], 因此采用相对压力为0的出口边界条件会导致较多的误差。

根据实际测量的时间-流量曲线进行计算,结果会更加真实反映动脉血液流动情况。

4 结论与展望

本文基于MRI影像重建的胸主动脉模型,在不同边界条件下计算讨论了出口入口边界条件对时间-流量曲线、壁面切应力、 AWSS和OSI值的影响。研究结果表明: 在不同的边界条件下主动脉内血液流动情况差异较大,这说明在血流仿真计算中应该尽量采用符合实际测量的血流流动参数,这样才能保证仿真结果与真实结果的吻合度。

本文模型在计算中没有考虑血管管壁弹性和运动,今后将进一步采用流固耦合(fluid-structure interaction, FSI)计算方法进行仿真计算,以获取更加符合实际血液流动的结果。

The authors have declared that no competing interests exist.

参考文献
[1] Buchanan J, Kleinstreuer C, Hyun S, et al. Hemodynamics simulation and identification of susceptible sites of atherosclerotic lesion formation in a model abdominal aorta[J]. Journal of Biomechanics, 2003, 36(8): 1185-1196. [本文引用:2] [JCR: 2.496]
[2] 付文宇, 乔爱科. 基于 CT 图像的胸主动脉瘤模型数值模拟[J]. 科技导报, 2009, 27(20): 27-31.
FU Wenyu, QIAO Aike. Numerical simulation of human thoracic aorta based on CT images[J]. Sciense & Technology Revew, 2009, 27(20): 27-31. (in Chinese) [本文引用:1]
[3] Ku D N, Giddens D P, Zarins C K, et al. Pulsatile flow and atherosclerosis in the human carotid bifurcation: Positive correlation between plaque location and low oscillating shear stress[J]. Arteriosclerosis, Thrombosis and Vascular Biology, 1985, 5(3): 293-302. [本文引用:2] [JCR: 5.533]
[4] Soulis J V, Fytanidis D K, Papaioannou V Cet al. Oscillating shear index, wall shear stress and low density lipoprotein accumulation in human RCAs [C]// 3rd Micro and Nano Flows Conference. Thessaloniki, Greece, 2011: 22-24. [本文引用:2]
[5] Caro C, Fitz-Gerald J, Schroter R. Atheroma and arterial wall shear observation, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis[J]Proceedings of the Royal Society of London. Series B: Biological Sciences, 1971, 177(1046): 109-133. [本文引用:1]
[6] 向亚菲, 李功文, 汤敬东, . 人胸主动脉血液流动的三维数值分析 [J]. 同济大学学报: 医学版, 2010, 31(4): 19-23.
XIANG Yafei, LI Gongwen, TANG Jingdong, et al. Three-dimensional simulation of blood flow in human thoracic aorta[J]. Journal of Tongji University: Medical Science, 2010, 31(4): 19-23. (in Chinese) [本文引用:1] [CJCR: 0.675]
[7] 杨金有, 俞航, 单晶心, . 基于计算流体力学方法进行不同个体的主动脉弓内血流模拟对比分析[J]. 中国医学物理学杂志, 2010, 27(4): 2059-2062.
YANG Jinyou, YU Hang, SHAN Jingxin, et al. Comparative analysis of different individuals of the aortic arch blood flow simulation based on computational fluid dynamics methods[J]. Chinese Journal of Medical Physics, 2010, 27(4): 2059-2062. (in Chinese) [本文引用:1] [CJCR: 0.604]
[8] Shahcheraghi N, Dwyer H A, Cheer A Y, et al. Unsteady and three-dimensional simulation of blood flow in the human aortic arch[J]. Journal of Biomechanical Engineering, 2002, 124(4): 378-387. [本文引用:1] [JCR: 1.748]
[9] Fytanidis D K, Soulis J V, Papaioannou V C, et al. Oscillatory flow in human arteries [C]// 2010 10th IEEE International Conference on Information Technology and Applications in Biomedicine (ITAB). Hong Kong, China, 2010: 1-4. [本文引用:2]
[10] 王荣品, 梁长虹, 黄美萍, . 3. 0T MR 相位对比法测量正常人肺、 体循环血流[J]. 中国医学影像技术, 2011, 27(1): 61-65.
WANG Rongpin, LIANG Changhong, HUANG Meiping, et al. Quantitative measurement of hemodynamics of pulmonary and systemic circulation in healthy volunteers with phase-contrast MR imaging on 3. 0T system[J]. Chinese Journal of Medical Imaging Technology, 2011, 27(1): 61-65. (in Chinese) [本文引用:1] [CJCR: 0.943]
[11] Wang R P, Liang C H, Huang M P, et al. Assessment of aortopulmonary collateral flow and pulmonary vascular growth using a 3. 0T magnetic resonance imaging system in patients who underwent bidirectional Glenn shunting[J]. European Journal of Cardio-Thoracic Surgery, 2012, 41(6): 146-153. [本文引用:1] [JCR: 2.814]
[12] Laganà K, Balossino R, Migliavacca F K, et al. Multiscale modeling of the cardiovascular system: Application to the study of pulmonary and coronary perfusions in the univentricular circulation[J]. Journal of Biomechanics, 2005, 38(5): 1129-1141. [本文引用:1] [JCR: 2.496]
[13] Gallo D, De Santis G, Negri F D, et al. On the use of in vivo measured flow rates as boundary conditions for image-based hemodynamic models of the human aorta: Implications for indicators of abnormal flow[J]. Annals of Biomedical Engineering, 2012, 40(3): 729-741. [本文引用:1] [JCR: 3.231]
[14] Mori D, Yamaguchi T. Computational fluid dynamics modeling and analysis of the effect of 3-D distortion of the human aortic arch[J]. Computer Methods in Biomechanics & Biomedical Engineering, 2002, 5(3): 249-260. [本文引用:1] [JCR: 1.793]
[15] 杨金有. 胸主动脉内血液流动的计算流体力学方法研究 [D]. 沈阳: 中国医科大学, 2010.
YANG Jinyou. Computational Fluid Dynamics Method Analysis for Thoracic Aorta Blood Flow [D]. Shenyang: Chinese Medical University, 2010. (in Chinese) [本文引用:1]
[16] McDonald D A. Blood Flow in Arteries[M]. Baltimore, MD: Williams & Wilkins Company, 1974. [本文引用:1]