关键词: 主动脉瓣；有限元；应力场；约束反力

本文用有限元法分析了人主动脉瓣在正常生理状态下的应力场、附着缘的约束反力。非线性本构关系通过实验测定，但在实用上，它可以近似地用分段线弹性来模拟。文中考虑了瓣叶高度h对邻瓣间接触面积的影响、对应力分布和附着缘反力的影响。

分类号：R318.01

STRESS FIELD ANALYSIS OF AORTIC VALVE UNDERPHYSIOLOGICAL CONDITION

Li Zhangzheng, Yu Jianhua, Li Jintang, Chen Junkai

(Department of Civil Engineering & Applied Mechanics,Sichuan University, Chengdu, 610065)

ABSTRACT：We analysed the stress field and constrained force at attachment margin in human aortic valve by using finite element method. The nonlinear constitutive relation of the tissue was determined by experiment, however, could be approximately simplified as several linear elastic segments corresponding to various Young′s moduli E. We also considered the effects of valve height on the contact area between the two neigbouring leaflets, stress distribution in the valve, and reaction at attachment margin of the tissue.

Key words:Aortic valve; Finite element; Stress field; Constrained force

1瓣叶几何形状描述

零应力时瓣叶的参考构型极为复杂。Angell等人^［1］研究证实，动脉瓣的解剖模型变化较大。然而，Mer-cer等人^［2］则指出，在心脏舒张期，瓣叶形状大致为半个抛物面。Mohamed等人^［3］详细研究了猪瓣叶，测量结果的表明，在瓣膜关闭、处于未受力状态时，单个瓣叶的形状接近半个椭圆抛物面。本文采用后者的结论，即以椭圆抛物面(图1)来描述瓣叶：

式中a、b为常数。瓣叶顶部与主动脉血管的连结点(附着点)用A、A′来表示，瓣底附着在血管的B点，该三点位于x₁、x₃平面(与竖直面xz夹θ角)内。椭圆的长轴为AA′方向，俯视图如图2所示。三个瓣叶均分圆面积，由于对称关系，应有∠ACA′=2∠BCA=120°。参数a、b由瓣的高度h(竖向)和血管半径R确定。

设z=h,即x₃=h/cosθ时，椭圆的长轴长度等于A、A′两点之间的距离，则由(1)式有如下关系存在：

已知h、R的值，就可依据(2)式来求得参数a。另一参数b,从几何上讲应满足如下不等式：

b值大小决定着瓣叶表面面积，换言之，若已知瓣叶表面积，便可用迭代的办法求出b。此处取(3)式的中间值，b=0.78a。

图沿是瓣在铅直面xz平面内的投影，它实际存在一条弧形的关闭线(图中虚线AOA′)。关闭线以上，与邻瓣相遇而形成接触面(ACA′OA),ACA′为游离缘；关闭线以下，朝向心室为非接触面。接触面COAC和COA′C分别按竖直平面计算，非接触面满足椭圆抛物面方程式(1)。

1椭圆抛物面图2瓣在水平面上的投影

2有限元单元划分

用有限元进行分析时，因为对称关系，所以只取一个瓣叶的二分之一，如图2中的BACB或图3中的BOCAB部分。瓣为空间曲面，用三角形单元进行模拟比较合适^［4］。如图4所示，将欲分析的半个瓣叶BCAB共划分成77个网格结点、121个单元。其中关闭线OPA以上的接触面分为50个单元，关闭线以下的非接触面分为71个单元。

3瓣在竖直面xz内投影图4单元划分

整体坐标取为xyz(如图1所示)，单元网格采取自动划分。接触面(竖向平面)内的结点坐标直接由平面方程在整体坐标下计算，非接触面内结点坐标由(1)式计算并按下式转换成整体坐标。

x=x₁,y=x₂cosθ+x₃sinθ,z=x₃cosθ-x₂sinθ(4)

从解剖学、生理学^［5，8］上可知，对中国人而言，R的平均取值为：R=12.5mm；h的值作为一个可变的量，计算时取不同数值。这样就可以比较h变化对接触面面积、附着缘长度及瓣内应力场等的影响。

三角形薄膜单元，位移u_i取结点值的线性插值函数，即：

式中N_k为形函数，u^k_i为结点k在i方向的位移。

3本构关系及平衡方程

作者等人对人主动脉瓣的本构关系进行了实验研究^［6］，由各向同性、不可压缩超弹性模型得到三维本构方程，其中一维拉伸应力-应变关系可表述为：

这里G=3.868 MPa为材料的剪切模量，λ=1+ε为试件的长度比(现实长度L与初始长度L₀的比值)，σ为柯西应力。由上式可定义切线模量E(λ)如下：

由(7)式可知，E是λ的非线性函数。从实用的角度出发，可以将它进行分段线性化处理^［3］。在10%的应变范围内，可均分成5个线性段，其相应的折线模量E分别为：11.6、13.5、18.0、26.6和42.0MPa。

非线性问题，可采取增量加载、分段线性的办法^［4］进行处理。对于增量应力和增量应变，仍然满足广义胡克定律。泊松比μ可取值为0.4995(真正的不可压缩材料μ=0.5)。

设时间增量Δt，对应的柯西应力、应变增量分别为σ^Δt_ij和e^Δt_ij，则由虚位移原理得平衡方程：

这里，δ表示变分，V为现时构型的体积，而W^Δt代表外力的虚功。该方程不能直接求解，因为现时构型V未知。由更新的拉格朗日方法，对(8)式进行推导，可得增量平衡方程^［3］：

［K］_T^。｛u｝=｛ΔP｝(9)

式中［K］_T为切线刚度矩阵，它是单元刚度的集成。｛u｝为增量位移矢，｛ΔP｝为增量力矢量。解(9)式，即可得增量位移，进而可求得增量应变、增量应力。

4计算结果

设瓣所受外载荷仅为血压，取正常值16.0 kPa(120mmHg)。R=12.5mm,h分别为18.0、20.0和22.0mm。附着缘上的边界结点固定，对称轴线BOC上单元结点约束x方向位移，而接触面上各结点只可沿z方向移动。

4.1接触参数

主动脉半径确定之后，关闭线的位置将只随h而变。h不同取值时，瓣间接触面积(一个完整瓣)、中心位置接触线OC的长度，如表1所示，它们随h的增大而增大。三个瓣叶游离缘的总长度(6R)及接触面总面积决定着瓣的启闭功能^［7］。接触面积大，关闭性能好；若接触面积过小(h较小)，则将会出现关闭不全的病症。

一个瓣叶附着缘ABA′的总长度S，可沿着弧长积分而求得：

部分结果见表1。(12)式过于复杂，不便于实际临床应用。在一定范围内，我们由(12)式通过电脑算出50余对数据，在三维曲面图S=S(R,h)上发现，它十分接近于一个平面。为此，进行了线性拟合，结果如下：

S=1.2572R+1.6546h+0.4052(mm)(13)

(13)式与(12)式相比，当R≤h≤2R时，其相对误差＜0.6