基于混合有限元的心肌纤维动力学仿真方法

徐宇　袁永峰　王宽全

摘 要：心脏的节律性收缩是其泵血功能的动力学基础，研究其收缩行为可以为心脏循环系统运转提供科学指导。心脏的收缩是心脏被动力和主动力共同作用的结果。被动力行为与心肌组织的材料性质有关，主动力行为心脏受到电生理调控。本文基于混合有限元方法建立心肌纤维弹性力学模型，仿真实验表明该模型可以有效仿真心脏组织在体力、牵引力和主动力作用下的形态学变化。综上本文提出的方法可以为今后建立电力耦合模型和心脏动力学研究提供理论指导。

关键词：心肌纤维力学建模；有限弹性力学；混合有限元

中图分类号：TP391.41 文献标识码：A 文章编号：2095-2163（2015）05-

Abstract： The power of heart pumping origins from rhythmic contraction of the myocardial tissue. Research of heart contraction is benefit to understand the circulation system. Passive property of heart， which is determined by materials features of myocardial issue， in combination with the active force instigated by electrical activity that leads to its contraction. In this paper， it provides that based on mixed finite element methods myocardial fiber mechanic model is developed to simulate elastic features of heart contraction. The simulated results suggest that this model can effectively simulate mechanic actions of heart under body force， traction and active force. The method provide a theoretical guide to set up electro-mechanic coupling model of heart and cardiac mechanic behaviors.

Keywords： Myocardial Fiber Mechanic Model； Finite Elasticity； Mixed Finite Element Methods

0引 言

心臟的计算模型可以分为电生理模型和力学模型两大类。心脏电生理模型的研究始于Noble[1]，其中是将Hodgkin[2]关于电兴奋在枪乌贼神经上传导的模型运用到心脏的浦氏纤维中。此后，心脏电生理建模取得了长足的发展。然而，心脏的力学模型研究滞后于电生理模型[3-4]，在心脏力学建模领域存在着诸多尚未解决的难题。目前，国外的文献虽然提出了心脏力学的建模方法，但很多文献存在不一致和不准确的地方，比如在处理不可压缩约束时，某些文献在应力表达式加入了项 ，有些文献加入了项 ；有些文献将弱形式和基于最小能量泛函的变分形式混用。另外，在模型中考虑心脏的不可压缩性[5]会给模型的数值求解带来难度，研究人员主要通过罚函数法[6-8]和多域有限元法[9，10]来解决该问题。罚函数法主要用于基于位移的有限元框架中，其形式简单，只需要计算位移变量，然而在材料趋向于不可压缩时，该方法会出现锁定现象[11]。多域有限元方法能够很好地避免锁定现象，对心脏不可压缩性质进行良好建模。基于以上考虑，本文首先严格按照有限弹性力学的方法建立了心脏的力学模型，然后给出了基于混合有限元方法的模型求解流程，最后通过实验对方法进行了测试。

1基于有限弹性力学的心脏力学模型

心脏由左右心房、左右心室四个腔室构成，左心室是四个腔室中体积最大的，左心室收缩时产生的高压将血液泵送至左心房[12]。为了能承受高压，左心室的心室壁比较厚，分别由心内层，心肌层和心外层构成。本文的力学模型主要针对左心室，并在左心室的一个六面体形状的切块上进行仿真分析。

图1为心脏切块形变示意图，用 表示未发生形变的心脏切块，其上的点记为 ，称该坐标为材料坐标[13]；用 表示发生形变后的心脏切块，其上的点记为 ，称该坐标为空间坐标[13]。将初始构形中的质点 映射到当前构形中的质点 的函数称为形变函数：

接着需要考虑如何度量心脏形变程度，本文使用有限弹性力学中的Cauthy应变张量 来衡量心脏形变的大小，Cauthy应变张量 的表达式为：

公式（6）中的 为Cauthy应变张量， 为二阶单位张量。

心肌组织发生变形后，组织内各部分之间产生了相互作用的内力，本文将引进有限弹性力学中的应力来描述心肌组织间的内力。心肌组织在内力和外力的共同作用下达到动态平衡，在本文的模型中，则只是考虑具有两个状态的拟静态平衡。平衡方程是联系心脏组织的运动学行为与心脏组织的材料性质的桥梁，根据动量守恒定律可以建立心脏组织的应力平衡偏微分方程：

公式（7）中， 表示散度算子， 为形变梯度张量， 为第二Piola-Kirchhoff应力张量[14]， 为心脏单位体积受到的外力。另外，根据角动量守恒定律可以知道，第二Piola-Kirchhoff应力张量 为对称张量，利用该性质可以在模型的求解中避免不必要的张量运算，节省时间开销。

本文在心脏的材料坐标系下进行模型求解，所以问题的求解区域为心脏的初始构形 。研究将 的边界划分为互不相交的两部分 和 ，其中 表示狄利克雷边界，该边界上的质点位移值为一给定的函数； 表示纽依曼边界，该边界上的质点受到沿着边界法向的牵引力，牵引力与心脏组织的形变无关。边界条件的数学描述如下：

公式（10）中， 為形变梯度张量的行列式。

2 混合有限元求解方法

混合有限元方法通过在应变能函数中引进拉格朗日乘数项来满足不可压缩的限制，并且把静压力作为与位移域独立的未知变量，这样可以有效地避免“锁定现象”[8]的发生。该方法对应的应变能函数为：

公式（19）中， 表示关于初始位置 的梯度， 表示张量的对称化运算， 为弹性张量[17]，其它符号含义同上文。对于线性系统（17），本文使用迭代法[17]求解。

3 数值实验

本文的数值实验使用C++语言和有限元库deal.II[15]实现。在数值实验中，选用 （cm）的心脏切块，切块被划分成 的网格，如图2（a）所示，总的自由度为2 443，切块的 边界为狄利克雷边界，规定该边界上的位移为0，除此之外没有其它边界条件。选用文献[11]中的指数应变能函数进行仿真，其数学形式可表述为：

本文设计了如下三组实验：

（1）体力 心脏切块受到的体力 N/cm3

（2）牵引力 牵引力作用于 平面，大小为0.004Kpa，方向与平面的内法向一致；

（3）主动力 主动力为 （KPa）其中 表示纤维的伸缩率，实验中选择的纤维方向 ；

图3给出了主动力与纤维伸缩率之间的关系。从图3中可以看出两者之间具有线性关系，这与表达式 的对应效果是一致的，由此表明仿真结果是正确的。图4给出了三种测试条件下，体积比率（即当前体积与初始体积之比）在牛顿迭代过程中的变化情况，由图4可知在迭代过程趋向于收敛的情况下，体积比率逐渐趋向于1，这与不可压缩限制完全相符。

4 结束语

本文基于混合有限元方法建立了心肌纤维弹性力学模型，并对心脏组织在体力、牵引力和主动力作用下的形态学变化进行了仿真，实验结果表明，混合有限元方法克服了基于位移有限元方法在求解不可压缩问题出现的锁定问题，可以有效仿真心肌组织的不可压缩性。另外，本文对主动力的仿真则为今后进行电力耦合仿真提供了理论基础。

参考文献：

[1] NOBLE D. A modification of the Hodgkin—Huxley equations applicable to Purkinje fibre action and pacemaker potentials[J]. The Journal of Physiology， 1962， 160（2）： 317-352.

[2] HODGKIN A L， HUXLEY A F. A quantitative description of membrane current and its application to conduction and excitation in nerve[J]. The Journal of physiology， 1952， 117（4）： 500-544.

[3] Arís Sánchez R. Electromechanical large scale computational models of the ventricular myocardium[J/OL].Universitat Politecnica de Catalunya， 2014：149[2014-11-06].http：//hdl.handle.net/10803/284398.

[4] Kirk N R. An adaptive， preconditioned， electromechanical model for the simulation of cardiac arrhythmias[D]. Leeds： The University of Leeds， 2012.

[5] COSTA K D， HOLMES J W， MCCULLOCH A D. Modelling cardiac mechanical properties in three dimensions[J]. Philosophical Transactions of the Royal Society of London A： Mathematical， Physical and Engineering Sciences， 2001， 359（1783）： 1233-1250.

[6] USYK T P， LEGRICE I J， MCCULLOCH A D. Computational model of three-dimensional cardiac electromechanics[J]. Computing and visualization in science， 2002， 4（4）： 249-257.

[7] VETTER F J， MCCULLOCH A D. Three-dimensional stress and strain in passive rabbit left ventricle： a model study[J]. Annals of biomedical engineering， 2000， 28（7）： 781-792.

[8] THORVALDSEN T， OSNES H， SUNDNES J. A mixed finite element formulation for a non-linear， transversely isotropic material model for the cardiac tissue[J]. Computer methods in biomechanics and biomedical engineering， 2005， 8（6）： 369-379.

[9] GUREV V， PATHMANATHAN P， FATTEBERT J L. A high-resolution computational model of the deforming human heart[J]. Biomechanics and modeling in mechanobiology， 2015： 1-21.

[10] F. Brezzi， M. Fortin， Mixed and Hybrid Finite Element Methods[M]. Heidelberg：Springer-Verlag， 1991.

[11] USYK T P， MAZHARI R， MCCULLOCH A D. Effect of laminar orthotropic myofiber architecture on regional stress and strain in the canine left ventricle[J]. Journal of elasticity and the physical science of solids， 2000， 61（1-3）： 143-164.

[12] HOLZAPFEL G A， OGDEN R W. Constitutive modelling of passive myocardium： a structurally based framework for material characterization[J]. Philosophical Transactions of the Royal Society of London A： Mathematical， Physical and Engineering Sciences， 2009， 367（1902）： 3445-3475.

[13] SIFAKIS E， BARBIC J. FEM simulation of 3D deformable solids： a practitioner's guide to theory， discretization and model reduction[C]//ACM SIGGRAPH 2012 Courses， Los Angeles： ACM， 2012： 20.

[14] Holzapfel G A. Nonlinear solid mechanics[M]. Chichester： Wiley， 2000.

[15] HOLZAPFEL G A， NIESTRAWSKA J A， OGDEN R W， et al. Modelling non-symmetric collagen fibre dispersion in arterial walls[J]. Journal of The Royal Society Interface， 2015， 12（106）： 20150188.

[16] LINGE S， LINES G， SUNDNES J. Solving the heart mechanics equations with Newton and quasi Newton methods–a comparison[J]. Computer methods in biomechanics and biomedical engineering， 2005， 8（1）： 31-38.

[17] FORTIN M， TARDIEU N， FORTIN A. Iterative solvers for 3D linear and nonlinear elasticity problems： Displacement and mixed formulations[J]. International Journal for Numerical Methods in Engineering， 2010， 83（13）： 1780-1802.

[18] BANGERTH W， HARTMANN R， KANSCHAT G. deal. II-a general-purpose object-oriented finite element library[J]. ACM Transactions on Mathematical Software （TOMS）， 2007， 33（4）： 24.