现在,心脏瓣膜疾病已经成为流行性疾病之一。全世界每年约有20万例患者需要进行主动脉瓣的手术治疗^[1],主要采用人工瓣膜进行替换。人工心脏瓣膜分为生物瓣膜和机械瓣膜,其中机械瓣膜具有耐久性好、寿命长等优点,缺点是需要终身服用抗凝药物,可能出现血栓等风险^[2]。生物瓣膜的优点是具有较好的生物相容性和良好的血流动力学性能,且随着微创介入疗法的成熟生物瓣膜将广泛应用。但是,生物瓣膜在体内服役时易出现钙化或纤维结构破坏,因此其服役寿命有限^[3]。

前期研究结果表明,生物瓣膜的钙化、破损、撕裂等与其材料结构和体内服役时经历的复杂应力应变有密切的关系^[4]。为了揭示瓣膜在体内运行环境下的微观受力情况,国内外学者针对生物瓣膜进行三维建模及有限元分析。袁泉等^[5]对不同型面的瓣膜进行静力学分析,认为球型瓣叶应力分布较为合理。马雪洁等^[6]建立了瓣膜与血液间的耦合计算模型,证明了积分法分析瓣膜应力的准确性。Weinberg等^[7]利用流固耦合有限元模型分析了不同年龄段人群的心脏瓣膜钙化、组织失效等情况。

人工瓣膜植入体内后很难检测其流体动力学特性,而脉动流检测能有效模拟体内血流环境,成为检测心脏瓣膜流体动力学性能的重要手段^[8,9]。因此,利用体外脉动流检测结合有限元分析方法可有效分析人工瓣膜的结构—力学性能间关系。本文利用人工生物瓣膜体外脉动检测数据,结合三维建模和有限元分析方法,研究生物瓣膜的结构—力学性能间的关系。

1 体外脉动流的实验检测

按照ISO-5840国际测试标准要求,采用ViVitro Labs Inc.脉动实验台对人工生物瓣膜(Edwards #2625)的血流动力学特性进行检测。测试温度为37℃,使用的流体为生理盐水。设置心循环率70 cycles/min,心输出量5.6 L/min,平均主动脉压100 mmHg,保证主动脉收缩压120 mmHg,舒张压80 mmHg,脉动前向流百分比35%。检测10个模拟周期,由测试软件ViVitest计算出生物瓣膜在上述参数条件下的脉动流特征参数,包括有效开口面积、返流百分比、平均跨瓣压差等。

2 有限元分析模型

2.1 几何模型

根据Smut^[10]研究的参数指标设计正常瓣膜的三维结构,表1列出了瓣膜模型设计的五个主要参考参数。利用UG NX 8.5(SIEMENS, Plano, USA)进行复杂曲面造型,图1给出了生成的瓣膜三维模型,并以IGES格式保存,假设此状态为零压力状态。

显示原图| 下载原图ZIP| 生成PPT

图1   生物瓣的几何模型

Fig.1   Geometric model of the bioprosthetic heart valve: (a) front view, (b) top view

2.2 网格模型

心脏瓣膜具有类似于薄膜的薄壁结构,瓣膜性能损伤的主因是弯曲变形,将瓣膜视为仅受张力的薄膜结构会限制弯曲效果。使用三维实体单元分析复杂模型的应力分布,因为三维实体单元在解决薄膜结构互相渗透时是比较复杂的,而壳元素的运动约束能力较容易,因此使用壳单元来代替三维实体单元分析瓣膜的受力作用^[11]。将心脏瓣膜的三维几何模型导入ANSYS有限元软件(ANSYS, Pittsburgh, USA)进行有限元网格划分,瓣膜由4103个壳单元(shell 163)组成,厚度为0.2 mm,模型的网格划分采用四边形网格。

2.3 材料属性、载荷设置及边界条件

真实的心脏瓣膜为各向异性超弹性材料,本文将生物瓣膜简化为各向同性线弹性材料,表2列出了主要的材料属性数值。

图2给出了边界设置示意图。生物瓣膜是固定在支架上的,其刚度相对支架较低,因此假定三个瓣叶是被约束在零位移的刚性支架上,瓣膜模型整个下边界范围(缝合边)自由度限制为零,上边界范围(自由边)是自由活动的^[12]。

显示原图| 下载原图ZIP| 生成PPT

图2   边界设置示意图

Fig.2   Boundary settings for the bioprosthetic heart valve

瓣膜承受来自左心室和主动脉两边的压力,将近似生理条件下脉动流检测得到的随时间变化的压差曲线作为有限元分析的输入压力载荷。

周围流体对心脏瓣膜的开合功能有重要影响,血液流速不仅对心动周期的打开阶段有明显的影响,且在关闭状态瓣膜因压差发生急速变化,周围液体的阻尼对瓣叶变形起主要作用。使用瑞利阻尼决定模型在液体中运动产生的粘性阻尼显示动力学分析,即

[C]=α∙M+β∙K(1)

σd=β∙Del∙ε(2)

其中, [C]为结构阻尼矩阵, M为结构质量矩阵, [K]为结构刚度矩阵, α为质量比例系数, β为刚度比例系数(粘性阻尼系数), σd为附加阻尼压力, Del为弹性刚度, ε为应变率。阻尼系数随响应频率而不断变化,刚度比例系数 β与响应频率呈线性增长关系,质量比例系数 α与响应频率成反比关系。为了避免高频响应的影响,基于Marom^[13]等对阻尼参数的分析,设置瓣叶刚度比例阻尼 β=0.01 s产生更优化的结果。

2.4 接触设置

为了保证瓣叶表面之间接触时不相互穿透,在表面间建立一种接触关系以防止在有限元分析中相互穿过。本文定义瓣叶之间接触行为为无摩擦接触,利用罚函数法来设置接触刚度在接触界面强制协调性以满足不穿透条件。通过接触面上探测点计算接触面的相互作用力^[14]。

2.5 图像处理分析

选取收缩期间三个时间点,将截取的模拟计算结果和实验图片分别导入Image J软件进行图像处理分析。先对图像进行标定,计算图片像素对应实物面积比值c。然后读取图像,将彩色图像转换为灰度图像,选定瓣膜开口阈值范围,生成二值化区域。进行开运算,通过腐蚀、膨胀算法处理使瓣膜开口与周围干扰分离,且腐蚀膨胀的程度相同不改变区域面积。根据形态学原理选取瓣膜开闭状态下的面积范围、圆度特征(瓣膜区域面积/最小外接圆面积)及瓣膜所在位置,得到瓣膜开口区域,对此区域进行面积计算,乘以比值c得到瓣膜的真实面积。

3 结果和讨论

3.1 脉动流检测结果

图3给出了一个心动周期近似生理条件的脉动压穿过瓣膜时其上下表面的压力曲线。一个心动周期的持续时间为0.85 s,平均主动脉压力保持在100 mmHg。在收缩期间,当t=0 s到t=0.07 s时左心室压力随着脉动流增大。当左心室压力增大到与主动脉相等时瓣膜开始打开,t=0.07 s到t=0.25 s为打开阶段。在舒张期间左心室压力随着脉动流减少而低于主动脉压力,关闭阶段开始,此时左心室与主动脉压差与打开阶段压差相比变化很大,特别在t=0.40 s到t=0.45 s期间压差呈现极具陡峭的梯度。随后压差改变较为平缓直到t=0.85 s,然后快速上升。在心动周期末期,瓣膜回到零压力加载状态。表3列出了生物瓣膜在体外脉动流检测下的参数指标与ISO-5840标准结果比较。有效开口面积2.0 cm²,跨瓣压差10.8 mmHg,返流百分比8.4%,均符合ISO-5840国际标准。

显示原图| 下载原图ZIP| 生成PPT

图3   心动周期过程的生物瓣膜检测曲线

Fig.3   Time-dependent blood pressure of the left ventricular and aortic pressure curve(a), and time-dependent transvalvular blood flow (b) for a bioprothetic heart valve in a cardiac cycle

3.2 有限元分析结果

本文结合体外脉动环境下测得的实验数据,为有限元分析提供近似生理条件下的参数边界条件,深入分析心动周期过程生物瓣膜应力应变分布情况。图4给出了生物瓣膜有限元分析时加载的压差曲线。生物瓣膜模型的零压力状态作为有限元分析的起始时间,周期为0.78 s,左心室和主动脉随时间变化的压力差作为显示动力学分析的压力载荷,曲线上压差的正值和负值指压力分别加载在瓣膜的左心室和主动脉一侧。图5、图6给出了心动周期中瓣膜在不同时间点下最大主应力和位移的分布情况。在打开阶段压力差从10.8 mmHg降到0,一个相对较小的压力加载在瓣膜的心室一侧,应力主要集中在瓣叶的缝合边及腹部弧形区域,最大主应力达到452 kPa。在关闭阶段,压力加载在瓣膜的主动脉一侧,压差值相较于打开阶段变得非常大,应力主要集中在缝合边两侧,即自由边与缝合边交界部分,最大主应力达到1.46 MPa。在心动周期末期,瓣膜回到零压力加载状态。

显示原图| 下载原图ZIP| 生成PPT

图4   随时间变化的压差加载曲线

Fig.4   Loading time-dependent curve of transvalvular pressure

显示原图| 下载原图ZIP| 生成PPT

图5   心动周期过程最大主应力分布情况(MPa)

Fig.5   Maximum principal stress distribution among the heart valve leaflets during the cardiac cycle (MPa)

显示原图| 下载原图ZIP| 生成PPT

图6   心动周期过程位移分布情况(mm)

Fig.6   Displacement distribution among the heart valve leaflets during the cardiac cycle (mm)

Halevi等^[15]对人体主动脉瓣膜样品钙化情况分析结果显示,瓣膜腹部是钙化的主要区域,有限元分析揭示其也是主要应力集中区域之一。植入到人体中的生物瓣膜要受到跨瓣压力、循环弯曲应力和剪切应力等多重作用。近期研究表明,人工生物瓣膜的钙化与瓣膜在血液中的受力情况密切相关^[16]。机械应力对生物瓣膜钙化的影响是一个非常复杂的过程,甚至不同方向的应力都对钙化效果产生不同的影响;不同应力的联合作用,可能是协同的或拮抗的^[17]。对人工生物瓣膜的微观受力的有限元分析,将为其体内钙化机制提供一些参考。

Mohammadi等^[18]研究在瓣叶表面的节点利用罚函数法计算推动力,解决了瓣叶之间接触问题。本文使用显示动力学中罚函数算法,适于瓣叶大变形分析,在心动过程中垂直于瓣叶间接触表面施加推动作用力,在瓣膜关闭阶段瓣叶间最大接触反力约为1.58 N,避免瓣叶之间渗透。

因为周围血液带来的阻尼效应影响动态模拟结果的准确性,Kim等^[19]设定阻尼常数为18.9 g/s,Kamensky等^[20]设定每单位面积上阻尼常数为80 g/(cm² s),以确保生物瓣膜在生理条件下正常打开。本文充分考虑惯性力和阻尼力对瓣膜动力学性能的影响,综合以上研究结果设置粘性阻尼系数为0.01 s。结果表明,生物瓣膜可在近生理条件下正常开合。

3.3 结果的比较

使用分辨率1080 p、帧频率60 fps的相机拍摄生物瓣膜在体外脉动实验运动状态视频,图7给出了截取t=0.05 s,t=0.20 s,t=0.50 s三个时间点生物瓣膜脉动流实验和有限元仿真情况的开合状态,图像处理后得到三个时间点实验条件下开口面积分别为0.27 cm²、1.32 cm²、0 cm²,仿真计算的开口面积分别为0.30 cm²、1.46 cm²、0 cm²。结果表明,实验数据与仿真计算结果相接近。

显示原图| 下载原图ZIP| 生成PPT

图7   收缩期和舒张期不同时间点瓣膜开口面积的脉动流实验和模拟仿真结果

Fig.7   Geometric orifice areas from experimental tests (a-c) and FE simulation (d-f) at different time points: (a, d) 0.05 s and (b, e) 0.20 s during the systole period, and (c, f) 0.50 s during the diastole period

在近似生理条件下人工瓣膜的有效开口面积是其重要功能指标之一,Pibarot等^[21]指出当置换瓣膜的有效开口面积低于正常人瓣膜的有效面积时,会产生瓣膜不匹配。有研究人员根据流体的定常流动情况推导瓣膜开口面积^[22],但是瓣膜打开和关闭的速度及开口面积受到血流量、跨瓣压差等因素影响,根据公式计算的结果与实验检测数据有一定的差别。Gabbay等^[23]指出,在高跨瓣压差和大的心搏出量情况下瓣膜有效开口面积计算值会偏高。Gillinov等^[24]指出,在体外测试条件下瓣膜的有效开口面积受到脉动流波形、流速等因素影响,因此难以通过理论计算给出准确的数值。本文采用脉动流实验过程中实时捕获瓣膜开闭过程的图像,经图像处理获得开口面积,结果比理论计算方法更为准确,且与实验结果相匹配,进一步验证了仿真计算方法的可靠性。

4 结论

(1) 近似生理条件下的脉动流实验能准确检测人工心脏瓣膜的流体动力学性能,本文使用ViVitro脉动实验台在脉动流条件下对Edwards人工生物瓣膜进行检测评估,主要参数均符合ISO-5840国际检测标准。

(2) 在生理压力载荷下生物瓣膜有限元模型揭示其在收缩期最大主应力达到425 kPa,应力集中在弯曲变形严重的腹部;舒张期最大主应力1.46 MPa,应力集中在瓣叶缝合边的两侧。局部应力集中和大变形对瓣膜钙化、结构失效起到重要影响。

(3) 根据体外脉动流检测病结合有限元分析方法可有效分析人工生物瓣膜的结构—力学性能间关系,能为人工生物瓣膜的结构设计和性能优化提供高效可靠的评价,为生物瓣膜在体内服役过程中出现的钙化和受力破损提供理论依据。

The authors have declared that no competing interests exist.