基于混合有限元的心肌纖維動力學仿真方法

徐宇　袁永峰　王寬全

摘 要：心臟的節律性收縮是其泵血功能的動力學基礎，研究其收縮行為可以為心臟循環系統運轉提供科學指導。心臟的收縮是心臟被動力和主動力共同作用的結果。被動力行為與心肌組織的材料性質有關，主動力行為心臟受到電生理調控。本文基于混合有限元方法建立心肌纖維彈性力學模型，仿真實驗表明該模型可以有效仿真心臟組織在體力、牽引力和主動力作用下的形態學變化。綜上本文提出的方法可以為今后建立電力耦合模型和心臟動力學研究提供理論指導。

關鍵詞：心肌纖維力學建模；有限彈性力學；混合有限元

中圖分類號：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.