基于计算流体力学的骨骼肌自激活变形的数值仿真方法与流程

本发明属于生命健康医疗，涉及基于计算流体力学的骨骼肌自激活变形的数值仿真方法。

背景技术：

1、骨骼肌是人体的重要构成部分，起着保护人体内部组织、牵引人体骨骼运动等功能。包含骨骼肌的人体动力学数值模型在人体运动学、医疗器械开发、临床疾病诊疗等领域中都发挥着关键作用。因而，开发合适的数值方法来准确的模拟骨骼肌的真实变形行为对人体动力学模型的发展及应用都起着重要作用。

2、在骨骼肌自激活力学变形行为的数值仿真中，关键的一步是准确提取并模拟骨骼肌纤维在三维空间中的排布。针对这一挑战性问题，目前主要有3种方法。第一种方法是从骨骼肌核磁图像数据(mri)中计算得到骨骼肌纤维的排布信息(damon et al.,diffusion-tensor mri-based skeletal muscle fiber tracking.)，即纤维各空间点的方向矢量，然后将这些信息传入到骨骼肌的有限元数值分析模型中，进行后续的力学变形分析。该方法的优点是纤维方向矢量是基于临床影像计算得到的，即信息来源于实际，但从mri影像中计算各纤维的方向矢量是个挑战性工作，目前还没有软件能自动实现该功能，不利于方法的广泛应用。第二种数值仿真方法是基于nurbs方法重建三维骨骼肌的几何，然后从nurbs实体中计算出各纤维的方向矢量，随后将该信息输入给骨骼肌的数值模型中(lu et al.,modelling skeletal muscle fibre orientation arrangement.)。相比于第一种方法，该方法可以更高效地计算出三维骨骼肌的纤维排布方向矢量信息，但构建骨骼肌三维nurbs实体模型仍是一较繁琐的过程，致使该方法并没有被广泛采纳。第三种获取三维骨骼肌纤维排布的方法是采用计算流体动力学分析得到三维骨骼肌内部流场的分布，并以此近似为骨骼肌的纤维分布(inouye et al.,fiber tractography for finite element modelingof transversely isotropic biological tissues of artibrary shape usingcomputational fluid dynamics.)。已有研究证明，该方法可以很好地近似具有贯穿式长纤维的骨骼肌纤维分布，并且，相较于前两种方法，该方法可以较容易地计算出骨骼肌纤维分布信息，但是，目前开发的方法都需要混用多种商业分析软件，操作较繁琐且不同软件之间的信息传递较复杂。因而，针对骨骼肌纤维排布的数值模拟目前还需要一种简便且高效的方法。

技术实现思路

1、为解决上述问题，本发明提供一种基于计算流体力学及材料本构子程序的骨骼肌自激活力学变形行为的高效仿真模型方法。

2、为了达到上述目的，本发明采用的技术方案是：

3、基于计算流体力学的骨骼肌自激活变形的数值仿真方法(简称cfd-fe)，包括如下步骤：

4、第一步：采用图像处理软件，从骨骼肌核磁图像(mri)中重建骨骼肌的三维几何模型；

5、第二步：采用有限元网格将第一步得到的骨骼肌三维几何模型进行离散化，获得骨骼肌有限元数值网格模型。

6、第三步：进行三维骨骼肌的计算流体动力学分析；

7、基于第二步得到的骨骼肌有限元数值网格模型，进行流体动力学分析来得到骨骼肌纤维分布信息，具体如下：将骨骼肌的上附着端定义为流体流入端，骨骼肌的下附着端定义为流体流出端，骨骼肌壁面定义为无摩擦滑动接触，进而在有限元分析软件的cfd板块中进行稳态流动分析；分析结果输出的流场分布线即假定为骨骼肌的纤维分布，通过调整有限元网格的尺寸大小和密度来改变cfd结果文件中的流场密度，进而得到一种满足计算精度要求的最适宜的网格密度，最终获得骨骼肌三维流场数据。

8、第四步：骨骼肌纤维分布信息的传递；

9、本步骤的目的是将cfd分析得到的骨骼肌三维流场数据传递给有限元分析软件的静力学板块进行骨骼肌力学变形分析，该步骤使用有限元分析软件结合编程软件进行。首先用编程软件从cfd结果文件中提取出各节点上描述流场分布的方向矢量信息。随后，由于骨骼肌固体本构模型中应力、应变的计算是在骨骼肌高斯点上进行，因而，从骨骼肌固体本构模型节点方向矢量数据库中搜寻出离各高斯点坐标最近的节点，并将该节点的流场方向矢量视为对应最近高斯点上的骨骼肌纤维方向。最后，将骨骼肌纤维方向信息传递给静力学板块中骨骼肌本构模型，进行后续的力学变形分析。

10、第五步：纤维分布信息应用于骨骼肌本构模型；

11、基于有限元分析软件的材料用户子程序功能开发了骨骼肌本构模型：将骨骼肌看作是一纤维增强的复合材料，其中基体是由不可压缩的软组织构成，因而其力学行为为不可压缩的超弹性力学行为；骨肌中的纤维是其主要构成部分，负责自激活骨肌收缩行为，使用经典的hill三元件模型来模拟。为了描述这些力学行为，基于超弹性材料力学本构理论，从骨骼肌的能量函数出发推导并建立其力学本构关系，具体如下：

12、骨骼肌的总能量u写为：

13、u＝ui+uf+ui (1)

14、其中，ui为骨骼肌基体中的能量；uf为骨骼肌纤维中的能量；uj为骨骼肌体积变化产生的能量。

15、骨骼肌基体中的能量ui和体积变化产生的能量uj为：

16、

17、

18、其中，是消除了体积变化的cauchy-green应变张量的第一不变量，b和c是材料参数，j是变形梯度的雅可比矩阵，d是压缩常数。

19、骨骼肌纤维中的能量uf为：

20、

21、其中，σs为串联弹性元件(see)的应力，σp为平行元件(pe)的应力，两者组成纤维中的总应力；是纤维拉伸比，λs是see中的拉伸比。在求解出三个能量函数ui、uj和uf后，再对能量函数进行求导便可以得到本构关系。

22、同时施加于收缩元件(ce)的肌肉激活函数表示如下：

23、

24、ht(ti,tb)＝{1-exp[-s·(ti,tb)]} (6)

25、其中，n1是激活前后的激活水平，n2是激活期间的激活水平，t0是激活时间，t1是失活时间，s是指数因子。

26、针对具体要分析的骨骼肌，其有限元网格采用第三步中cfd分析时使用的网格，仅需要将网格类型从流体网格改为固体网格，骨骼肌的材料模型采用上述骨骼肌本构模型。在骨肌自激活行为的仿真中，将骨骼肌的上附着端进行固支约束，使用本构关系中的肌肉激活函数来控制骨骼肌的收缩行为。

27、本发明的有益效果：本发明方法首先通过计算流体动力学数值仿真分析得到三维骨骼肌内部流场的分布，利用骨骼肌内部流场来近似骨骼肌纤维分布，然后利用程序代码将流场信息输入到骨骼肌有限元数值模型中的各高斯点，进而进行后续的力学变形仿真分析。本发明简单易行，全部流程可在常见有限元分析软件的计算流体动力学(cfd)及静力学模块完成，因而，一套有限元网格模型就能完成所需的数值计算，大大简化了前处理中纤维分布信息传入cae模型的步骤并降低了复杂度。