基于混合有限元的心肌纖維動(dòng)力學(xué)仿真方法

徐宇　袁永峰　王寬全

摘 要：心臟的節(jié)律性收縮是其泵血功能的動(dòng)力學(xué)基礎(chǔ)，研究其收縮行為可以為心臟循環(huán)系統(tǒng)運(yùn)轉(zhuǎn)提供科學(xué)指導(dǎo)。心臟的收縮是心臟被動(dòng)力和主動(dòng)力共同作用的結(jié)果。被動(dòng)力行為與心肌組織的材料性質(zhì)有關(guān)，主動(dòng)力行為心臟受到電生理調(diào)控。本文基于混合有限元方法建立心肌纖維彈性力學(xué)模型，仿真實(shí)驗(yàn)表明該模型可以有效仿真心臟組織在體力、牽引力和主動(dòng)力作用下的形態(tài)學(xué)變化。綜上本文提出的方法可以為今后建立電力耦合模型和心臟動(dòng)力學(xué)研究提供理論指導(dǎo)。

關(guān)鍵詞：心肌纖維力學(xué)建模；有限彈性力學(xué)；混合有限元

中圖分類號(hào)：TP391.41 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：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引 言

心臟的計(jì)算模型可以分為電生理模型和力學(xué)模型兩大類。心臟電生理模型的研究始于Noble[1]，其中是將Hodgkin[2]關(guān)于電興奮在槍烏賊神經(jīng)上傳導(dǎo)的模型運(yùn)用到心臟的浦氏纖維中。此后，心臟電生理建模取得了長(zhǎng)足的發(fā)展。然而，心臟的力學(xué)模型研究滯后于電生理模型[3-4]，在心臟力學(xué)建模領(lǐng)域存在著諸多尚未解決的難題。目前，國(guó)外的文獻(xiàn)雖然提出了心臟力學(xué)的建模方法，但很多文獻(xiàn)存在不一致和不準(zhǔn)確的地方，比如在處理不可壓縮約束時(shí)，某些文獻(xiàn)在應(yīng)力表達(dá)式加入了項(xiàng) ，有些文獻(xiàn)加入了項(xiàng) ；有些文獻(xiàn)將弱形式和基于最小能量泛函的變分形式混用。另外，在模型中考慮心臟的不可壓縮性[5]會(huì)給模型的數(shù)值求解帶來難度，研究人員主要通過罰函數(shù)法[6-8]和多域有限元法[9，10]來解決該問題。罰函數(shù)法主要用于基于位移的有限元框架中，其形式簡(jiǎn)單，只需要計(jì)算位移變量，然而在材料趨向于不可壓縮時(shí)，該方法會(huì)出現(xiàn)鎖定現(xiàn)象[11]。多域有限元方法能夠很好地避免鎖定現(xiàn)象，對(duì)心臟不可壓縮性質(zhì)進(jìn)行良好建模?；谝陨峡紤]，本文首先嚴(yán)格按照有限彈性力學(xué)的方法建立了心臟的力學(xué)模型，然后給出了基于混合有限元方法的模型求解流程，最后通過實(shí)驗(yàn)對(duì)方法進(jìn)行了測(cè)試。

1基于有限彈性力學(xué)的心臟力學(xué)模型

心臟由左右心房、左右心室四個(gè)腔室構(gòu)成，左心室是四個(gè)腔室中體積最大的，左心室收縮時(shí)產(chǎn)生的高壓將血液泵送至左心房[12]。為了能承受高壓，左心室的心室壁比較厚，分別由心內(nèi)層，心肌層和心外層構(gòu)成。本文的力學(xué)模型主要針對(duì)左心室，并在左心室的一個(gè)六面體形狀的切塊上進(jìn)行仿真分析。

圖1為心臟切塊形變示意圖，用 表示未發(fā)生形變的心臟切塊，其上的點(diǎn)記為 ，稱該坐標(biāo)為材料坐標(biāo)[13]；用 表示發(fā)生形變后的心臟切塊，其上的點(diǎn)記為 ，稱該坐標(biāo)為空間坐標(biāo)[13]。將初始構(gòu)形中的質(zhì)點(diǎn) 映射到當(dāng)前構(gòu)形中的質(zhì)點(diǎn) 的函數(shù)稱為形變函數(shù)：

接著需要考慮如何度量心臟形變程度，本文使用有限彈性力學(xué)中的Cauthy應(yīng)變張量 來衡量心臟形變的大小，Cauthy應(yīng)變張量 的表達(dá)式為：

公式（6）中的 為Cauthy應(yīng)變張量， 為二階單位張量。

心肌組織發(fā)生變形后，組織內(nèi)各部分之間產(chǎn)生了相互作用的內(nèi)力，本文將引進(jìn)有限彈性力學(xué)中的應(yīng)力來描述心肌組織間的內(nèi)力。心肌組織在內(nèi)力和外力的共同作用下達(dá)到動(dòng)態(tài)平衡，在本文的模型中，則只是考慮具有兩個(gè)狀態(tài)的擬靜態(tài)平衡。平衡方程是聯(lián)系心臟組織的運(yùn)動(dòng)學(xué)行為與心臟組織的材料性質(zhì)的橋梁，根據(jù)動(dòng)量守恒定律可以建立心臟組織的應(yīng)力平衡偏微分方程：

公式（7）中， 表示散度算子， 為形變梯度張量， 為第二Piola-Kirchhoff應(yīng)力張量[14]， 為心臟單位體積受到的外力。另外，根據(jù)角動(dòng)量守恒定律可以知道，第二Piola-Kirchhoff應(yīng)力張量 為對(duì)稱張量，利用該性質(zhì)可以在模型的求解中避免不必要的張量運(yùn)算，節(jié)省時(shí)間開銷。

本文在心臟的材料坐標(biāo)系下進(jìn)行模型求解，所以問題的求解區(qū)域?yàn)樾呐K的初始構(gòu)形 。研究將 的邊界劃分為互不相交的兩部分 和 ，其中 表示狄利克雷邊界，該邊界上的質(zhì)點(diǎn)位移值為一給定的函數(shù)； 表示紐依曼邊界，該邊界上的質(zhì)點(diǎn)受到沿著邊界法向的牽引力，牽引力與心臟組織的形變無關(guān)。邊界條件的數(shù)學(xué)描述如下：

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

2 混合有限元求解方法

混合有限元方法通過在應(yīng)變能函數(shù)中引進(jìn)拉格朗日乘數(shù)項(xiàng)來滿足不可壓縮的限制，并且把靜壓力作為與位移域獨(dú)立的未知變量，這樣可以有效地避免“鎖定現(xiàn)象”[8]的發(fā)生。該方法對(duì)應(yīng)的應(yīng)變能函數(shù)為：

公式（19）中， 表示關(guān)于初始位置 的梯度， 表示張量的對(duì)稱化運(yùn)算， 為彈性張量[17]，其它符號(hào)含義同上文。對(duì)于線性系統(tǒng)（17），本文使用迭代法[17]求解。

3 數(shù)值實(shí)驗(yàn)

本文的數(shù)值實(shí)驗(yàn)使用C++語言和有限元庫deal.II[15]實(shí)現(xiàn)。在數(shù)值實(shí)驗(yàn)中，選用 （cm）的心臟切塊，切塊被劃分成 的網(wǎng)格，如圖2（a）所示，總的自由度為2 443，切塊的 邊界為狄利克雷邊界，規(guī)定該邊界上的位移為0，除此之外沒有其它邊界條件。選用文獻(xiàn)[11]中的指數(shù)應(yīng)變能函數(shù)進(jìn)行仿真，其數(shù)學(xué)形式可表述為：

本文設(shè)計(jì)了如下三組實(shí)驗(yàn)：

（1）體力 心臟切塊受到的體力 N/cm3

（2）牽引力 牽引力作用于 平面，大小為0.004Kpa，方向與平面的內(nèi)法向一致；

（3）主動(dòng)力 主動(dòng)力為 （KPa）其中 表示纖維的伸縮率，實(shí)驗(yàn)中選擇的纖維方向 ；

圖3給出了主動(dòng)力與纖維伸縮率之間的關(guān)系。從圖3中可以看出兩者之間具有線性關(guān)系，這與表達(dá)式 的對(duì)應(yīng)效果是一致的，由此表明仿真結(jié)果是正確的。圖4給出了三種測(cè)試條件下，體積比率（即當(dāng)前體積與初始體積之比）在牛頓迭代過程中的變化情況，由圖4可知在迭代過程趨向于收斂的情況下，體積比率逐漸趨向于1，這與不可壓縮限制完全相符。

4 結(jié)束語

本文基于混合有限元方法建立了心肌纖維彈性力學(xué)模型，并對(duì)心臟組織在體力、牽引力和主動(dòng)力作用下的形態(tài)學(xué)變化進(jìn)行了仿真，實(shí)驗(yàn)結(jié)果表明，混合有限元方法克服了基于位移有限元方法在求解不可壓縮問題出現(xiàn)的鎖定問題，可以有效仿真心肌組織的不可壓縮性。另外，本文對(duì)主動(dòng)力的仿真則為今后進(jìn)行電力耦合仿真提供了理論基礎(chǔ)。

參考文獻(xiàn)：

[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.