4.3.2 主动脉瓣流固耦合模型构建
如前文所述，浸没边界法在处理有大变形、大位移，以及流体域可能出现被固体域隔断的流固耦合问题时，有不需要构造贴体网格、计算稳健可靠的特点，因此本研究也选择了浸没边界法作为流固耦合模型构建的方法。在浸没边界法中，流体域通常被处理成多材料单元类型的静态实体单元网格，允许同一网格中有多种材料并存；而固体域采用拉格朗日方法，使用的实体单元或壳单元、膜单元网格可以直接嵌入流体网格中自由移动而无须在流固耦合接面上构造贴体网格；流固耦合界面之间的计算过程采用罚函数法，仿真过程中求解器会搜索固体网格在流体网格中发生穿透的位移变化，然后给固体网格的节点施加弹性回复力，以试图将其“推回原位”，计算过程中固体域节点的质量和动能也将按照距离远近的权重重新分配到周围液体的网格节点上，通过这样的办法对加速度和速度进行耦合。
1）计算网格
在本研究中，结构域（主动脉瓣及主动脉根部血管）和流体域（血液）可分别构建计算网格，而无须用血管和瓣膜几何对流体域进行切割和划分边界层。由于在本研究中，LSDYNA（Version R11.1.0，LSTC，美国；Ansys，美国）被选用为求解器，因此网格类型和一些前处理的方法均按照其格式要求完成。
结构域的构建主要依靠上一节中逆向重建完成的组织解剖结构模型。完成重建的患者个体化主动脉根和主动脉瓣几何被导入Hypermesh（Altair Engineering，美国）进行网格化处理。为了降低计算成本，在本研究中，瓣叶和血管都被简化成用三维全积分壳单元处理；而为了避免绑定接触（Tie）可能导致的局部应力过大问题，瓣叶和主动脉根部的网格进行共节点处理。
不同于主动脉根部和瓣叶的网格模型，血液部分的网格模型直接在Hypermesh中进行构建。首先将主动脉根部血管模型两端封闭并在入口端面划分四边形面网格，然后运用Hypermesh三维网格模块中的实体网格匹配（Solid map）工具，选取主动脉根部血管模型出入口分别作为源面和目标面，入口端面上的四边形网格作为拉伸网格，血管壁外表面曲面作为扫略控制面，网格尺度设置为约1 mm，拖曳出血管内血液六面体网格。由于心动周期过程中，主动脉窦和主动脉可能会反复舒张或收缩甚至有较大位移，为了保证所有结构在心动周期内都完整“浸没”在血液内而不会“逃出”流体域而导致血管内形成真空或塌陷，对上一步生成的血管网格沿径向继续偏置两层六面体网格作为流体域网格。
LS-DYNA中无法直接将压力载荷施加在流体域的实体单元上，而是使用与流体域连通的受压流体（pressurized fluid）层给流体域内的流体间接传递静水压力或者动态压力。因此，对血液模型位于主动脉根部血管两端的网格，分别额外地偏置出一层单元作为入口和出口处受压流体层，以添加载荷。
结构域网格和流体域网格在同一坐标系下被分别构建，再各自存为单独的DYNA关键字文件——K文件，以便根据几何或者材料属性的变化做局部修改或调整。当结构域和流体域K文件被计算主程序调用时，即自动装配成相互嵌套的流固耦合网格模型（图4-12）。
图4-12 采用浸没边界法构建的手术组患者个体化主动脉瓣流固耦合模型示意图
2）材料属性
本研究所涉及的生物材料包括原生瓣膜、主动脉根部血管和心包膜。这三种生物材料都由多层软组织构成，并由特定方向排列的胶原纤维增强，显示出高度的各向异性。原生瓣叶由包括心室层、海绵层和纤维层的三层结构组成，在纤维层中有宏观上正交于主动脉轴向的胶原纤维成束排列。自体心包膜由浆膜层和纤维层组成，在纤维层中亦有宏观上正交于心底-心尖连线的胶原纤维平行排列。主动脉血管壁也是三层结构，内层是由内皮细胞和纵向排列的结缔组织构成的内膜，中层是环形排列着弹性纤维和胶原纤维的中膜，外层是由纵向排列的结缔组织构成的外膜。三种不同组织在微观解剖下的结构特征如图4-13所示。
在本研究中，心包膜、原生瓣叶和主动脉血管都被认为是正交各向异性材料，其周向（circumferential）、径向（radial）和轴线（longitudinal）方向上的弹性模量E、剪切模量G、泊松比μ、密度和厚度信息来自Grand-Allen和Luo等的体外力学测试，具体数值见表4-1。
图4-13 原生瓣叶、原生自体心包膜及主动脉内的胶原纤维宏观排列
表4-1 正交各向异性的生物组织材料属性
注：Ecirc—周向弹性模量；Erad—径向弹性模量；Elong—轴向弹性模量；Gxy—xy平面内的剪切模量；Gyz—yz平面内的剪切模量；Gzx—zx平面内的剪切模量；μxy—xy平面内的泊松比；μyz—yz平面内的泊松比；μzx—zx平面内的泊松比。
材料的正交各向异性属性需要以正确的方向映射到有限元模型中各组织的网格单元上。在LS-DYNA中壳单元的正交各向异性主要通过定义材料局部坐标系来实现（图4-14），其中材料局部坐标系的a、b和c方向对应于材料属性表中的x、y、z方向，分别通过以下方法来定义：
图4-14 LS-DYNA中壳单元定义材料局部坐标系的示意图［13］
（1）沿主动脉轴线方向定义向量V，壳单元的面法线方向为n，则主动脉原生瓣叶、心包膜瓣叶替代物和主动脉根血管的a方向都可根据右手定则定义为向量V和其各自结构上壳单元的面法线方向的叉乘a=V×n。
（2）各壳单元材料的b方向定义为a方向和面法线方向的叉乘b=n×a。
（3）各壳单元材料的c方向则为c和b方向的叉乘方向。
（4）额外的，对于利用自体心包进行主动脉瓣修复的模型，可定义β角度来使整个坐标系绕着c方向旋转，这使控制自体心包主动脉瓣叶中纤维朝向相对于原生瓣膜纤维的周向排列发生偏转成为可能。
如前所述，本研究中的血液被设为不可压缩牛顿流体，其密度为1 050 kg/m3，动力学黏度为0.003 5 Pa·s。需要特别注意的是，为了提高LS-DYNA显式计算的速度，本研究中血液中的声速被定义为150 m/s。
3）载荷与边界条件
如前所述，在本研究中，变化的心室压力和主动脉压力需要通过流体域两端的受压流体层间接地施加到流体域上。在LS-DYNA中可使用状态方程（equation of state，EOS）来将受压流体层中的材料与它们所受到的压力、密度、体积和内能联系起来。
在本研究中，受压流体层的压力主要是由Gruneisen状态方程来控制的。Gruneisen状态方程是一种用来表示固体材料在特定温度条件下体积和压力之间关系的方程，但工程上研究凝聚态液体材料的冲击或高压问题时，往往也采用固体形式的高压方程。LSDYNA中常用的Gruneisen状态方程的形式如下，定义为可压缩材料上受到的压力：
式中 E——初始单位体积的内能，量纲与压强相同；
C——常温常压无扰动状态声速；(https://www.daowen.com)
α——常数，γ0的一阶体积修正系数；
S1，S2，S3——待定的常系数；
ρ0——初始密度；
μ——压缩状态系数，由下式确定：
本研究引用了正常的左心室压力和主动脉压力波形作为模板，并根据中山医院心外科术中测量的患者个体化的高低血压为依据进行缩放，从而得到个体化的左心室压力和主动脉压力波形（图4-15）作为载荷输入，分别施加在前述的受压流体层上。
图4-15 分别作为入口和出口脉动载荷的左心室压力（Pv）和主动脉压力（Pao）波形图
本研究中的流体域和固体域的边界条件做如下设置：对于血管部分，约束主动脉根部血管两端的节点；对于血液部分，约束流体域所有实体单元形成的整体外表面所有节点，并消除压力波在流体网格边界上发生折射、反射而影响计算结果。
另外，为了更稳健地计算瓣膜开闭过程中三片瓣叶反复相互接触和释放的复杂过程，本研究定义了带软约束（soft constraint）的点面接触（node-to-surface contact），并通过阻尼控制（damping control）来降低瓣膜在高速开闭过程中可能的震荡带来的计算不稳定性。
4）虚拟手术模型
本研究需要构建基于患者个体化主动脉根部几何的疾病瓣膜模型、正常模型和手术模型。疾病模型主要是对4.3.1节基于患者CT影像逆向重建的个体化主动脉瓣几何直接网格化获得。而正常模型和手术模型则通过对患者个体化主动脉根模型做出结构和/或材料上的调整获得。
通过对异常的主动脉根和主动脉瓣叶形态进行修复从而构建出正常模型。主动脉根的解剖结构异常被认为会影响瓣膜正常功能。对于上述患者个体化主动脉瓣疾病模型，从主动脉侧往下可观察到主动脉根部有异常膨出（图4-16），主动脉瓣瓣环也因此有一定的扭曲。所以，正常模型首先即在几何上继续利用Geomagic Studio 12.0对畸形的主动脉根部进行平滑处理，消除异常膨出，并光顺主动脉瓣瓣环的B样条曲线，以符合正常主动脉窦圆润饱满、瓣环曲线光滑简洁的解剖特征；然后再利用新的瓣环的空间曲线与前述定义的瓣叶自由边及脊线B样条曲线，构建正常瓣叶的几何形态。
在本研究中，手术模型的几何网格主要采用了虚拟手术的方法进行建模，以在仿真模型上模拟临床实际采用的自体心包主动脉瓣修复术的过程（图4-17a）。该虚拟手术主要分成四步（图4-17b），由于是单纯的结构力学仿真，所以采用Abaqus Standard（Version 2017，Dassault System，法国）完成：
图4-16 疾病模型和正常模型的构建
图4-17 自体心包主动脉瓣修复术的虚拟手术模拟
（1）首先按照临床常用的瓣叶裁剪模板构建合适该患者主动脉根部尺寸的扇贝状瓣叶片体。
（2）然后将该片体模型弯折、牵拉，使瓣叶片体根部与主动脉根曲线重合，瓣叶两侧直边与主动脉根部相邻两窦连合部位的交界棱线重合。
（3）接着对瓣叶根部和两侧网格边缘的节点进行位移约束，但是释放其转动约束，并从主动脉侧对三个瓣叶施加600 Pa均布压力载荷（大约4%的瓣膜关闭状态时的最大跨瓣压力），让三个瓣叶轻柔地接触交叠，形成略微关闭的形态，模拟主动脉瓣修复手术过程最后零应力状态下的瓣叶形态。
（4）最后将变形后的瓣叶网格导出，转化单元类型并保存为LS-DYNA的K文件格式，作为手术模型的瓣叶初始网格模型。
本研究也设计制作了3D打印的瓣叶快速装夹装置以验证上述虚拟手术方法各步骤的合理性（图4-18），以及最后形成的零应力瓣叶形态的准确性。该瓣叶快速装夹装置由内外镶嵌的套筒组成，内外套筒都呈三尖皇冠形状；内侧套筒壁为圆周完整的衬垫；外侧套筒的三尖特征凸起成柱，柱上有等分的键槽，瓣叶假体（图中用纸替代）被卡在内外套筒之间，瓣叶连合部位可穿过外侧套筒键槽至工装外部固定。使用该工装可实践上述虚拟手术关于瓣叶裁剪、弯折、根部固定、向内轻压和定型的过程。由于其操作简便、结果直观的优点，非常便于仿真模型与之对照。通过仿真模型与瓣叶快速装夹实践对比显示，本研究中所提出的手工缝制瓣膜的虚拟手术方法是可行可靠的。另外值得指出的是，对该瓣叶快速装夹工装可以进行进一步优化和改造，使其在未来的工作中也可用于快速装夹不同瓣叶设计进行体外流体力学测试，以简化制作完整瓣膜假体进行测试的复杂程度，并加速瓣膜设计开发的进程。
图4-18 用3D打印的瓣叶快速装夹装置验证虚拟手术方法
关于材料属性，正常模型、疾病模型和手术模型采用了不同的材料局部坐标系定义方式。对于疾病模型和正常模型，直接采用前文所述的结合沿血管轴线方向的向量及瓣叶各网格单元法向的方式，定义原生瓣叶上正交于血管轴线并呈环形排列的主材料方向。对于手术模型，则通过定义偏转角β分别为0°、45°和90°，让瓣叶假体材料的主方向与原生瓣叶环形排列的材料主方向发生偏转，以模拟实际进行自体心包主动脉瓣修复术时由于对心包膜任意方向的裁剪而带来的三种自体心包膜胶纤维固有方向与原生瓣叶胶原纤维优势方向之间可能发生的失拟偏斜现象（图4-19）。0°模型意味着心包膜内的胶原纤维排列方向和原生瓣叶纤维束优势方向完全相同的理想匹配状态；45°模型意味着心包膜内胶原纤维斜行偏转于原生瓣纤维束优势方向；90°模型则说明心包膜内胶原纤维方向垂直于理想方向而呈沿血管轴线方向排列。
图4-19 自体心包主动脉瓣修复术中瓣叶假体与胶原纤维排列方向的夹角
5）求解设置和求解环境
对心血管的流体力学相关仿真一般都始于零应力状态或者静息状态，在显式动态分析时，数值模型对周期性载荷在刚开始加载时的瞬时冲击的响应往往与多周期后消除了惯性效应的响应不太一致。特别对于心脏瓣膜的分析，瓣叶的运动形态及流场中血液动力学特征等的分布，也往往需要若干周期才能达到稳定或与正常生理条件下接近。在本研究中，其判断标准为主要分析变量前后两个周期内不发生显著变化。因此，为了避免这种显式动态分析下由于初始加载瞬时启动可能带来的惯性偏差，基于前期试算的结果对流固耦合模型加载两个心动周期的周期性脉动压力载荷，可使瓣叶运动学特征和流场形态达到周期性稳定，而在后处理时仅对第二个周期的结果进行分析。
本研究的流固耦合计算使用LS-DYNA的显式求解器进行。计算的硬件环境为如下配置的高性能服务器：CPU为10核心20线程的Intel Xeon，E5-2680 Xeon CPU主频为2.8 GHz，内存为32 GB。