CN117993254A - 一种基于非线性大变形力学原理的全心脏模拟仿真方法 - Google Patents

Abstract

本发明公开了一种基于非线性大变形力学原理的全心脏模拟仿真方法，属于计算力学与心脏医学的交叉领域。针对心脏实体模型，第一步：先对全心脏模型划分有限元网格，得到心脏有限元模型，在此基础上，在心房内部生成血液粒子；利用基于规则的方法，给心房心室上的网格节点赋予相应的心肌纤维方向，并在心内膜上生成导电网络系统；第二步：针对建立的全心脏模型，基于非线性大变形方法，利用心肌力电耦合模型、流固耦合模型、接触力学有限元模型三个主要的理论模型，通过大变形有限元编程，模拟电信号沿着导电网络传播，激发心脏发生宏观变形，并驱动血液运动的基本功能。本发明可用于心血管病的科研和教学展示。

Description

一种基于非线性大变形力学原理的全心脏模拟仿真方法

技术领域

本发明属于计算力学与心脏医学的交叉研究领域，具体涉及一种基于非线性大变形力学原理的全心脏模拟仿真方法，即心脏在电驱动下的除极、复极以及引发的心脏周期性变形、血液流动和瓣膜运动的全心脏计算机模拟仿真方法。

背景技术

人的心脏大概拳头那么大，一天跳动约10万次，泵出约5000升血液，是一个非常辛苦的器官。心脏任何一个部分出现异常，其他部分都会受到牵连。因此，在临床上，心血管病患者往往同时患有多种类型的心血管病，医学上用一个专门的词“合并”来描述这一现象，比如患者A患有主动脉瓣狭窄合并房颤和心室肥大。由于这一基本事实，心脏医学领域客观上需要一个全心脏模型来研究心血管病的发病机理，优化临床上的诊断治疗方案。这样的全心脏模型不但需要包括心房、心室、瓣膜、血液等心脏的主要结构，还要能够模拟心脏变形-电传导的力电耦合作用、血液和心肌的流固耦合作用以及心脏不同部分之间相互接触的接触力学作用等主要作用过程。基于核磁共振成像以及超声影像资料的三维模型重构技术为建立全心脏实体模型奠定了基础。同时，有限元网格划分技术以及非线性大变形有限元等一系列计算力学方法为研究心脏力电耦合、流固耦合、接触等方面的现象提供了技术手段，并为最终对心脏主要过程的高精度仿真和演示奠定了基础。

心脏跳动是一个多尺度多物理的复杂过程，涉及纳米尺度离子通道的周期性开合、微米尺度心肌纤维的收缩和松弛以及宏观尺度心脏的跳动，同时包括力学变形、电传导扩散以及化学反应。由于心脏过程的高度复杂性，为了降低模拟计算的难度，现有的心脏模型往往只考虑心脏的一部分结构或者一部分功能，比如左心室模型只考虑心脏左心室的电传导和变形行为；双心室电传导模型虽然同时考虑左右两个心室，但是该类模型将心脏视为刚体，只研究心电传导扩散行为；研究二尖瓣的模型常常将心房、心室理想化为圆柱腔，只考虑血液与瓣膜的相互作用。这些只模拟单一功能的局部心脏模型虽然具有相对较高的计算精度，但由于模型太简单太理想化，无法模拟真实复杂的心脏过程，无法研究临床上遇到的复杂现象，也无法全面展示心脏运行和发病的机制。

发明内容

针对现有技术中存在的问题，本发明提供一种基于非线性大变形力学原理的全心脏模拟仿真方法，将全心脏几何建模、全心脏有限元大变形力电耦合模拟、血液心肌流固耦合模拟以及相邻部分的接触有限元模拟等一系列基于力学原理的计算力学技术结合在一起，从而能够模拟并展示心脏的所有主要结构和主要功能，并与临床上观察到的主要心脏现象进行对比分析。

为达到上述目的，本发明采用如下技术方案：

一种基于非线性大变形力学原理的全心脏模拟仿真方法，包括以下步骤：

步骤一：利用有限元网格划分方法对全心脏模型划分有限元网格，构建全心脏模型，所述全心脏模型包括心脏肌肉的几何有限元模型、血液模型、心电传导网络模型和心肌纤维分布模型；在全心脏有限元模型基础上，在左右心房内部生成血液粒子，模拟回流到左、右心房的动脉、静脉血液；利用基于规则的方法，给心房心室上的有限元节点赋予相应的心肌纤维方向，并在心房、心室内膜上生成导电网络系统；

步骤二：针对第一步建立的全心脏模型，基于非线性大变形方法，利用心肌力电耦合模型、流固耦合模型、接触力学有限元模型，对全心脏模型进行大规模并行计算，所述计算包括：心脏肌肉的力电耦合计算、心脏肌肉与血液的流固耦合计算以及瓣膜相邻部分的接触有限元计算；通过大变形有限元编程，模拟电信号从窦房结发出，沿着导电网络逐步传播到心房、心室，在传播过程中，心肌被电信号激活，发生局部收缩变形，随着被激活区域越来越多，心脏发生宏观收缩、舒张变形，并将血液压入心室或者主动脉、肺动脉中，从而模拟心脏跳动驱动血液运动的基本功能；

步骤三：利用第二步计算过程中产生的有限元节点的位置和电势数据以及血液速度、压强统计数据，借助可视化软件展示电信号在心脏导电网络以及心肌中的传播、心脏肌肉被电信号激活之后的变形以及血液在心脏中流动的基本心脏过程。模拟结果包括电信号沿导电网络的传播、电信号激发心脏肌肉发生的心脏收缩、舒张变形、心脏肌肉电势演化，即除极、复极过程，以及血液在心脏变形过程中的流动。

进一步地，所述步骤一中，心脏肌肉的几何有限元模型通过对心脏实体模型划分有限元网格获得，心脏实体模型包括心房、心室、瓣膜、腱索以及与心脏相连的主动脉、肺动脉。

进一步地，所述步骤一中的血液模型的构建方法为：基于流体力学中的光滑粒子动力学方法，根据心脏肌肉的几何有限元模型给定的心房内表面范围，在左、右心房内各生成一定量(5万)的光滑粒子。

进一步地，所述步骤一中的心电传导网络模型由全心脏模型内表面的一系列有限元网格节点组成，用来模拟窦房结、房间束、房室结、左右束支和浦肯野导电网络；电信号以设定好的速度，从右心房代表窦房结的有限元网格节点出发，在右心房和左心房内部的房间束传播，到达代表房室结的有限元网格节点之后再沿着代表左右束支和浦肯野导电网络的有限元网格节点在左右心室中传播；电信号传播的位置与该位置处的有限元节点到代表窦房结的有限元节点之间的距离相关。

进一步地，所述步骤一中的心肌纤维分布模型是在心房和心室所有有限元节点处给定一个单位方向矢量，该单位方向矢量指向该处心肌纤维的分布方向，不同位置处的心肌有限元节点具有不同的单位方向矢量。

进一步地，所述步骤二中，心脏肌肉的力电耦合计算包括：通过牛顿迭代方法求解弱形式的心肌力平衡方程和电传导方程，在每一个迭代步中，求解每一个迭代步的刚度矩阵和余量，待余量收敛到一个特定的小量(如0.000001)之后，停止迭代，更新所有有限元节点的位置坐标和电势值，进入下一个时间步。

进一步地，所述步骤二中，心脏肌肉与血液的流固耦合计算包括：步骤一中生成的代表血液的光滑粒子按照流体力学原理运动，光滑粒子与步骤一中生成的心脏肌肉的几何有限元模型的网格之间通过法向排斥力避免光滑粒子侵入到心脏肌肉的几何有限元模型的网格内部，光滑粒子对心脏肌肉的几何有限元网格的表面压力作为边界条件，应用到刚度矩阵和余量的计算中。

进一步地，所述步骤二中，瓣膜相邻部分的接触有限元计算包括计算全心脏模型相邻部分之间的接触作用，包括四个瓣膜的瓣叶之间的接触的计算，计算过程中，考虑瓣膜瓣叶之间的相互排斥力，避免一部分瓣叶侵入另一部分瓣叶的有限元网格中，产生非物理的现象，并将瓣膜相邻区域的有限元网格之间的相互作用力应用到刚度矩阵和余量的计算中。

进一步地，所述步骤三中，演示模拟结果包括心脏肌肉的力电耦合计算、心脏肌肉与血液的流固耦合计算以及瓣膜相邻部分的接触有限元计算的计算结果，即每一计算时间步心脏每一个有限元网格节点的位置和电势、光滑粒子的位置、速度和压强，将每一个时间步心脏的变形和电势变化、血液流动和瓣膜开合状态进行演示。

进一步地，所述步骤一包括：基于核磁共振成像或者超声技术，可以得到患者心脏的一系列二维影像资料，通过三维重构技术和有限元网格划分技术，建立患者心脏的有限元模型，包括心房、心室、瓣膜、乳头肌、腱索以及与心脏相连的主要动脉、静脉血管；在获得的有限元模型基础上，提取心房、心室初始状态下内表面的构型数据，以心房心室内表面为边界，在其内部，通过计算机编程生成一定量的模拟血液行为的光滑粒子。

进一步地，所述步骤二包括：在非线性大变形有限元的牛顿迭代计算框架下，在每一个增量步，利用力电耦合模型、流固耦合模型、界面接触模型计算系统整体刚度矩阵和余量，通过牛顿迭代方法使余量趋于一个小量(0.000001)，从而得到心脏肌肉有限元节点的位置信息以及电势演化信息，并利用光滑粒子动力学方法，得到血液流体光滑粒子的位置、速度等相关力学量。

进一步地，所述步骤三包括：将步骤中每一个时间步的计算结果保存到相应的文件中，利用开源可视化软件展示电信号在心脏中的传播并激活心脏肌肉收缩，引发心电演化，触发心脏跳动、并驱动血液流动的心脏主要过程。

与现有技术相比，本发明的有益效果在于：

(1)本发明建立包括血液在内的心脏所有主要结构的全心脏模型，这是通过计算模拟方法全面系统研究心脏行为以及从全局角度研究心脏发病机理的前提；

(2)本发明可以同时模拟心脏的力电耦合、流固耦合、接触三个主要心脏过程，即可以模拟心脏在电信号的传播扩散过程中激活心脏肌肉并引发心脏收缩、舒张的变形过程；可以模拟血液在心脏收缩压力驱动下流向并穿越瓣膜进入动脉静脉血管或者从心房进入心室；可以模拟瓣膜在血液冲击之下打开并在血液流过之后瓣叶重新贴合的过程；

(3)本发明通过同时模拟这些重要的心脏过程，可以研究各种心脏因素之间的相互联系，通过调节相关参数，可以研究各种病理心脏的运动规律，以及心脏发病的机理等一系列具有临床意义的重要问题；

(4)本发明利用计算过程中生成的数据，用开源可视化软件非常直观的展示复杂的心脏物理、力学过程，给患者和心脏医学相关从业人员解释相关心脏病的发病机理。

综上所述，本发明可以定量的给出心脏在任意心动周期(正常人约600毫秒)内每一时刻心脏每个区域的变形和电传导状态，即心肌电势激活和传导状态，同时，可以通过调控心肌力学、电学性能参数来模拟不同病理情况(如心梗、心律失常)下的心脏变形和电传导行为；利用本发明不但可以研究一些典型心脏病，如结构性心脏病、电紊乱引发的心脏病的发病机理，同时可以更加形象直观的展示正常心脏以及病理心脏的运行机理，为心血管病的科研和教学提供研究和展示方法。

附图说明

图1为本发明的一种基于非线性大变形力学原理的全心脏模拟仿真方法的流程图。

具体实施方式

针对一般心脏实体模型，本发明首先利用有限元网格划分方法对全心脏模型划分有限元网格，从而得到全心脏有限元模型；心脏组织学专家利用解剖和磁共振成像方法研究发现心房心室上的心肌纤维分布具有一定的规律性，例如，心室的肌肉纤维按照螺旋的方式排布，按照这种心肌纤维分布规律，通过计算机编程在心房心室模型的有限元节点上生成纤维分布矢量的方法被称为“基于规则(rule-based)”的方法。利用“基于规则”的方法，给心房心室上的节点赋予相应的心肌纤维方向，并在心房、心室内膜上生成导电网络系统；然后在第一步生成的全心脏有限元模型基础上，在左右心房内部生成血液粒子，模拟回流到左、右心房的动脉、静脉血液；最后针对第一步建立的包括血液、心肌纤维分布、导电网格在内的全心脏有限元模型，基于非线性大变形方法，利用心肌力电耦合模型、流固耦合模型、接触力学有限元模型三个主要的理论模型，通过大变形有限元编程，模拟电信号从窦房结发出，沿着导电网络逐步传播到心房、心室，在传播过程中，心肌被电信号激活，发生局部收缩变形，随着被激活区域越来越多，心脏发生宏观收缩、舒张变形，并将血液压入心室或者主(肺)动脉中，从而模拟心脏跳动驱动血液运动的基本功能。

下面按照图1所给的流程，从构建全心脏模型、大规模并行计算和模拟结果可视化三个方面对本发明的技术方案进行更加详细、完整地描述。

如图1所示，本发明的一种基于非线性大变形力学原理的全心脏模拟仿真方法包括如下步骤：

步骤1、构建全心脏模型，全心脏模型包括几何有限元模型、血液模型、心电传导网络模型和心肌纤维分布模型四个部分。

几何有限元模型是需要首先获得的最基本部分。如果要研究心脏的一般运动规律，只要购买一般的心脏实体模型即可，如美国Zygote公司出售的全心脏模型就具有非常精细的内部结构。如果要研究具体患者的心脏，则需要用患者心脏的核磁共振或超声影像资料，借助三维重构技术，将心脏二维断层扫描图重构为三维心脏实体。目前该过程仍然需要深度人工参与，借助建模软件来实现，是一个耗时耗力的工作。建立的全心脏实体模型具有心房、心室、瓣膜、乳头肌以及主要的动脉静脉血管等心脏主要结构。然后利用有限元网格划分软件给全心脏实体模型划分四面体有限元网格，即可获得全心脏的几何有限元模型，该模型包括约30万个四面体单元、12万左右个节点以及处于模型内表面和外表面上的三角形面单元。通过几何结构的分析，从四面体几何有限元模型中分割出心房、心室的内表面三角单元集合，即确定出心房和心室内表面所包围的封闭空间。

血液模型的建立是在左心房和右心房的内表面空间中分别随机生成一定数量(5万)的光滑粒子用以在下一步计算中模拟血液的行为。

心电传导网络模型是建立在心房心室内表面三角有限元节点基础之上。首先在右心房内表面的节点中按照窦房结和房室结的大致位置，选择两个相应的内表面节点代表窦房结和房室结节点，即模拟电信号从代表窦房结的有限元节点发出，在心房内表面的导电网络中传播，当电信号到达代表房室结的有限元节点之后，从该有限元节点沿着左右束支分别传导到左右心室的导电网络中。在右心房的内表面上，在代表窦房结的有限元节点和代表房室结的有限元节点之间随机生成三条互不相交的传导线路，并将三条传导线路上的内表面有限元节点标记出来。按照同样的做法，在左心房内表面上，也随机生成三条导电线路，并标记出线路上的有限元节点。在左心室和右心室的室间隔内表面上，分别生成两条导电网络，连接代表房室结的有限元节点和心室底部的内表面节点，并从底部选择的节点出发，按照分叉树的形式分别在两个心室的内表面生成浦肯野导电网络线路，并且标记线路上的有限元节点，作为导电网络节点。这些导电网络节点将在电信号抵达时被激活，电信号进一步从这些有限元节点传导到心脏的其他部位。

心肌纤维分布模型的建立：心房心室上的心肌纤维对心脏的力学变形和电传导具有决定性的影响。按照心脏组织学专家通过解剖等手段研究发现的心房心室上的心肌纤维分布规律，通过在计算心脏领域内广泛使用的称为“基于规则”的算法，在第一步中生成的心脏有限元模型中，在每一个心房心室的有限元网格节点上生成一个单位矢量，该矢量与心脏组织学家在心脏相应位置处观察到的心肌纤维方向一致。

至此，就得到了包含血液、心肌纤维分布以及导电网络的全心脏有限元模型。

步骤2、大规模并行计算，包括力电耦合计算、流固耦合计算、接触有限元计算三部分。

力电耦合计算、流固耦合计算、接触有限元计算即心脏肌肉在电信号激发下发生的力电耦合计算、血液与心脏肌肉和瓣膜部分的计算、心脏相连部分(如瓣膜瓣叶)之间的接触有限元计算三个部分。由于整个计算问题中存在几何非线性、材料非线性和接触非线性三类非线性问题，求解过程采用牛顿迭代方法(Newton-Raphson)。

在力电耦合计算中，心脏肌肉每一处局部应力包含主动应力和被动应力两个部分，主动应力是由电驱动，当节点电势从-80毫伏跃迁到20毫伏时，心肌纤维会产生相应的主动应力，除了主动应力，心肌纤维还存在由力学领域常用的非线性大变形本构模型Neo-Hooken模型描述的被动应力部分。应力满足平衡方程。节点电势会随着时间推移按照电传导方程不断向周围心肌扩散。将力平衡方程和电传导方程写成相应的弱形式，在每一个时间步，写出系统的整体刚度矩阵和余量，通过牛顿迭代，使得余量逐步趋近于一个小量(0.000001)，从而求得所有有限元网格节点的位置矢量和电势标量值。

血液心肌的流固耦合计算也是在这一框架下进行。血液的光滑粒子对心肌或者瓣膜的作用计入到表面力，血液的光滑粒子按照光滑粒子动力学规律运动，心脏肌肉或者瓣膜对光滑粒子的作用计入边界光滑粒子的合外力，并将作用逐步传递到其他光滑粒子。心脏心房的一些部分常常与心室相互接触，心脏的四个瓣膜的瓣叶之间也相互接触，心脏相连部分之间的接触问题对正确模拟心脏基本功能非常重要。相互接触的两个部分都是柔性体，本发明采用罚函数法考虑两部分之间的接触作用力，在每一个迭代时间步内，当一部分的节点侵入到另一部分的三角形面单元内部后，利用罚函数法，计算法向接触力和接触刚度，并计入系统整体刚度和余量中，选择合适的罚参数，通过多次牛顿迭代，更新接触部分之间的接触深度，避免相邻部分之间出现非物理的穿插现象。

步骤3、模拟结果可视化，包括电信号沿导电网络的传播、心脏变形、电势演化、血液流动四部分。

从步骤2的大规模并行计算中，可以得到每一个时间步心脏所有有限元节点的位置和电势以及血液的光滑粒子的位置、速度、压强等力学量。借助开源可视化软件的显示功能，只要把并行计算得到的每一个时刻的数据写成可视化软件要求的数据格式，就可以直观的展示心脏在每一个时刻的变形、电势演化以及血液流体在心脏内部的流动行为。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，包括以下步骤：

步骤一：利用有限元网格划分方法对全心脏模型划分有限元网格，构建全心脏模型，所述全心脏模型包括心脏肌肉的几何有限元模型、血液模型、心电传导网络模型和心肌纤维分布模型；在全心脏有限元模型基础上，在左右心房内部生成血液粒子，模拟回流到左、右心房的动脉、静脉血液；利用基于规则的方法，给心房心室上的有限元节点赋予相应的心肌纤维方向，并在心房、心室内膜上生成导电网络系统；

步骤二：针对第一步建立的全心脏模型，基于非线性大变形方法，利用心肌力电耦合模型、流固耦合模型、接触力学有限元模型，对全心脏模型进行大规模并行计算，产生的有限元节点，所述计算包括：心脏肌肉的力电耦合计算、心脏肌肉与血液的流固耦合计算以及瓣膜相邻部分的接触有限元计算；通过大变形有限元编程，模拟电信号从窦房结发出，沿着导电网络逐步传播到心房、心室，在传播过程中，心肌被电信号激活，发生局部收缩变形，随着被激活区域越来越多，心脏发生宏观收缩、舒张变形，并将血液压入心室或者主动脉、肺动脉中，从而模拟心脏跳动驱动血液运动的基本功能；

步骤三：利用第二步计算过程中产生的有限元节点的位置和电势数据以及血液速度、压强统计数据，借助可视化软件展示电信号在心脏导电网络以及心肌中的传播、电信号激发心脏肌肉发生的心脏收缩、舒张变形、心脏肌肉电势演化，即除极、复极过程，以及血液在心脏变形过程中的流动。

2.如权利要求1所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤一中，心脏肌肉的几何有限元模型通过对心脏实体模型划分有限元网格获得，心脏实体模型包括心房、心室、瓣膜、腱索以及与心脏相连的主动脉、肺动脉。

3.如权利要求2所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤一中的血液模型的构建方法为：基于流体力学中的光滑粒子动力学方法，根据心脏肌肉的几何有限元模型给定的心房内表面范围，在左、右心房内各生成一定量的光滑粒子。

4.如权利要求2所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤一中的心电传导网络模型由全心脏模型内表面的一系列有限元网格节点组成，用来模拟窦房结、房间束、房室结、左右束支和浦肯野导电网络；电信号以设定好的速度，从右心房代表窦房结的有限元网格节点出发，在右心房和左心房内部的房间束传播，到达代表房室结的有限元网格节点之后再沿着代表左右束支和浦肯野导电网络的有限元网格节点在左右心室中传播；电信号传播的位置与该位置处的有限元节点到代表窦房结的有限元节点之间的距离相关。

5.如权利要求2所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤一中的心肌纤维分布模型是在心房和心室所有有限元节点处给定一个单位方向矢量，该单位方向矢量指向该处心肌纤维的分布方向，不同位置处的心肌有限元节点具有不同的单位方向矢量。

6.如权利要求1所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤二中，心脏肌肉的力电耦合计算包括：通过牛顿迭代方法求解弱形式的心肌力平衡方程和电传导方程，在每一个迭代步中，求解每一个迭代步的刚度矩阵和余量，待余量收敛到阈值之后，停止迭代，更新所有有限元节点的位置坐标和电势值，进入下一个时间步。

7.如权利要求6所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤二中，心脏肌肉与血液的流固耦合计算包括：步骤一中生成的代表血液的光滑粒子按照流体力学原理运动，光滑粒子与步骤一中生成的心脏肌肉的几何有限元模型的网格之间通过法向排斥力避免光滑粒子侵入到心脏肌肉的几何有限元模型的网格内部，光滑粒子对心脏肌肉几何有限元模型的网格表面压力作为边界条件，应用到刚度矩阵和余量的计算中。

8.如权利要求6所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤二中，瓣膜相邻部分的接触有限元计算包括计算全心脏模型相邻部分之间的接触作用，包括四个瓣膜瓣叶之间的接触计算，计算过程中，考虑瓣膜瓣叶之间的相互排斥力，避免一部分瓣叶侵入另一部分瓣叶的有限元网格中，产生非物理的现象，并将瓣膜相邻区域的有限元网格之间的作用刚度和相互作用力应用到刚度矩阵和余量的计算中。

9.如权利要求6-8之一所述的一种基于非线性大变形力学原理的全心脏模拟仿真方法，其特征在于，所述步骤三中，演示模拟结果包括心脏肌肉的力电耦合计算、心脏肌肉与血液的流固耦合计算以及瓣膜相邻部分的接触有限元计算的计算结果，即每一计算时间步心脏每一个有限元网格节点的位置和电势、光滑粒子的位置、速度和压强，将每一个时间步心脏的变形和电势变化、血液流动和瓣膜开合状态进行演示。