基于ABAQUS的升主动脉置换术术后血流数值模拟

孙凯顺，张洪明，张桂敏，夏健明，孙毅

1.昆明理工大学工程力学系，云南昆明650500；2.云南省阜外心血管病医院，云南昆明650032

前言

ABAQUS是一款具有强大非线性分析能力和模拟复杂系统高可靠性的有限元软件。在计算流体软件中，ABAQUS/CFD模块有着强大的模拟计算能力，能模拟层流与湍流，稳态与瞬态的内部和外部流动。雷诺数范围很宽，对复杂几何有着很强的适应性。有限元法和混合有限体积法的求解方法，使其在不可压缩流的层流与湍流问题上有着较高的求解精度［1］。主动脉夹层（Aortic Dissection,AD）是一种严重的心血管疾病，主动脉内壁的内膜撕裂口发生后，血液流入主动脉中膜外层或中外膜交界处，可引起急性主动脉破裂而致患者死亡。或因夹层真腔受到假腔的挤压而致腔体狭窄和闭塞，肠管、肾脏、下肢等由真腔供血的重要脏器发生缺血性改变，并伴随严重并发症［2］。Stanford分类法将主动脉夹层分为A型和B型。Stanford A型：夹层破口累及升主动脉，无论第一破口源自何处；Stanford B型：夹层第一破口不累及升主动脉。当今Stanford A型发病率显著升高，严重威胁人民生命健康。主动脉外科手术的术式种类繁多，目前推行按照内膜破口位置、数量和夹层涉及范围制定个性化手术方案［3］，人工血管置换术就是其中一种有效治疗主动脉夹层的术式［4］。深入研究主动脉血流动力学规律的意义重大，通过体内实测获得一系列主动脉血流动力学指标的分布较为困难。现如今医学影像学、计算流体动力学的发展以及高性能计算机并行计算的应用，使数值模拟仿真的研究方法实现体外模拟主动脉血流成为高效手段，进而通过后处理呈现相关流域的血流动力学指标［5］。主动脉血流数值模拟的研究已有较多进展。边界条件的差异直接影响模拟计算的收敛性［6］。关于主动脉入口流速，较多的研究采用理想化流速［7］。研究表明理想化的流速分布会使模型输出与体内真实血液动力学参数产生一定偏差［8］。流型假设上，层流被认为适用于人类动脉中的血流［9-10］，且被广泛运用到主动脉血流模拟［11-12］，实际上主动脉血流存在湍流［13］。本文研究对象为一Stanford A型主动脉夹层病例，破口位于升主动脉，治疗方式采用人工升主动脉置换术。获取患者术后胸部CT图像中的AD数据以构建流域几何模型［14］，主动脉入口流速使用术后实测流速，湍流模型采用Spalart-Allmaras，该湍流模型收敛性较好，可呈现二次流等基本流动特征［15］。经模型网格划分后进行数值模拟。

1 计算模型和数值方法

1.1 控制方程

假定血液流动为不可压缩牛顿流体的粘性流动，其连续性和动量方程的无量纲形式可表示如下［16］：

其中，u为流速矢量，p为流体压力，Re为流体雷诺数。

1.2 数值模拟

采集本研究病例的术后胸部CT图像AD数据，在高性能计算机上应用MIMICS软件进行三维重构，经过自动化阈值分割和手动分离，删除所需模型以外的组织，得到包括主动脉根部、升主动脉、主动脉弓、降主动脉以及主动脉弓上部的3大分支（头臂干、左颈总动脉、左锁骨下动脉），再进行模型表面光顺化，获得术后流场CAD模型。

划分网格是流体数值模拟的重要一步，精细的网格使计算准确性较高，但会增加计算机硬件成本和时间成本。将CAD模型文件导入网格划分软件，如图1所示，划分高度复杂几何体网格有如下3步主要步骤：（1）几何预处理，将几何体抽壳后留下壁面块状曲面，将壁面块状曲面之间的共享边转化为抑制边，可实现网格划分时壁面块状曲面被视为整体连续曲面。（2）在连续曲面的基础上进行2D网格划分和网格质量控制。（3）在2D网格的基础上进行CFD网格划分，并通过网格质量检查。边界层网格为三菱柱单元，非边界层网格为四面体单元［17］。为了在保证良好收敛性的同时降低计算量，选用如图2所示的网格模型，共计单元数为292 012。

图1 主动脉流域几何模型Fig.1 Geometric model of aortic drainage area

图2 主动脉流域网格模型Fig.2 Mesh model of aortic drainage area

将上述网格文件导入高性能计算机的ABAQUS。主要前处理设置如下：（1）应用幅值曲线命令（Amplitude）的Smooth类型，将患者术后主动脉进口实测流速进行瞬态加载，如图3所示。（2）根据文献［18］，血液密度设定为1 060 kg/m3，血液运动粘度设定为3.71 mPa·s。（3）所有壁面假定为无滑移边界条件［19］。（4）本文的流体域设定为不可压缩牛顿流体［20］。（5）根据流动充分发展条件将主动脉弓3大分支出口和主动脉出口设为自由出流。（6）定义雷诺数Re=ρvd/μ，其中ρ为血液密度，v为主动脉进口流速，d为进口截面的特征长度，μ为血液运动粘度，得平均雷诺数为9 900，可知该流动存在湍流。使用上述湍流模型，进行10个心动周期的计算，共计6.6 s，采用自动分析步长，得到稳定收敛解，取第10个周期进行分析。

图3 主动脉进口流速曲线Fig.3 Aortic inlet velocity

2 结果与讨论

经过10个周期的模拟计算，获得主动脉壁面压力云图和主动脉纵向剖面流速矢量云图。取第10个周期6个关键时刻的计算结果进行分析，即t1=0 s、t2=0.105 6 s、t3=0.217 8 s、t4=0.231 0 s、t5=0.257 4 s、t6=0.369 6 s，这些时刻分别是开始射血时刻、射血早期时刻、射血加速最快时刻、射血峰值时刻、射血速度平稳时刻、射血降速最快时刻。

2.1 血管壁面压力分布

6个关键时刻的壁面压力云图如图4所示。主动脉壁面压力与边界几何形状和血流速度密切相关。在心脏收缩期，研究对象的壁面压力随着入口血流速度增大而迅速升高。在心跳周期不同时间点，流域壁面压力整体呈现从主动脉近心端到远心端的下降趋势，局部压力有差异分布。并且在流速较高时段，由于壁面几何改变，血流流动使升主动脉流域和主动脉弓流域的外侧压力分布明显高于内侧。当入口血流速度到达0.231 s的最高流速2 m/s时，流域整体壁面压力到达峰值。随后进入舒张期，整体壁面压力伴随入口血流速度的减少而降低。

图4 主动脉壁面压力分布（Pa）Fig.4 Wall pressure distribution at the aorta(Pa)

2.2 流域流速分布

6个关键时刻的主动脉纵向剖面流速矢量云图如图5所示。本例的入流流速存在病理性，从最低流速0.031 7 m/s到峰值流速2 m/s的循环中呈现非线性上升下降，入流平均流速为1.008 m/s。在心动收缩期初期，流域整体流速不高，升主动脉与主动脉弓流域流速较小，并伴随着如图所示的二次流，分别为升主动脉中部和降主动脉起始处内侧的二次回流和二次环流。心脏收缩期属于快速射血期，随着入口流速的提高，流域整体流速逐渐提高，升主动脉内侧流速与主动脉弓流速涨幅较小，且存在由于壁面复杂几何和二次流状态而表现出的局部性和波动性，升主动脉的二次流出现了向主动脉内侧偏移的趋势。到心室最大收缩期0.231 s，入口流速到达峰值时，流域整体流速保持上升，不同位置点的流速高峰时刻不同，这说明流场流速与入口流速没有同步性。当进入心室减慢射血期和舒张期，流域流速持续降低，升主动脉涡流和回流重新回靠。在几何关系上，涡流较多出现在局部几何复杂或腔体狭窄的动脉中，如心脏瓣膜、静脉瓣膜、血管分支处以及血管急转弯的内侧边，可知本例存在两处典型急转弯，为主动脉根部与升主动脉交界处内侧和主动脉弓与降主动脉交界处内侧。流速而论，流速对于涡流形成的影响呈现双向作用。当血流流速较低时，血液粘滞性可致流动阻力较大，轴流不明显而边流作用增加，涡流形成可能性提高。当流速较高时，血流剪切力增加且血流冲击内壁，也可致涡流出现的几率上升。

图5 主动脉纵向剖面流速矢量分布（m/s）Fig.5 Blood velocity vector distribution at aortic vertical section(m/s)

3 结论

血流力学行为与模型几何形状和进出口边界条件息息相关。本文介绍了ABAQUS/CFD在主动脉血流数值模拟上的应用，得到主动脉壁面压力云图和剖面流速云图。相比健康成年人体，本例的壁面压力分布和流速分布整体上偏大，需要术后进一步治疗调理。在不同时刻，从近心端到远心端，壁面压力整体呈现下降趋势，局部压力有差异分布，较大压力位置为主动脉根部内侧和升主动脉起始处外侧。二次流的主要分布位置为升主动脉内侧和降主动脉起始处内侧。心血管疾病的发病机理和病理过程与血流力学环境的异常密不可分。局部力学属性的突变，包括血液的湍流、二次环流、二次回流、血管内壁局部高压等，会使血细胞、脂粒等附壁堆积（或导致动脉粥样硬化类疾病）、血管内膜受损、血管局部扩张等现象发生。本例的血流数值模拟可为人工血管置换术治疗Stanford A型主动脉夹层后的病情发展研究提供参考。