CN106227993A - 一种基于生物学机理的骨折愈合动态过程仿真方法 - Google Patents

Abstract

一种基于生物学机理的骨折愈合动态过程仿真方法，本发明涉及基于生物学机理的骨折愈合动态过程仿真方法。本发明的目的是为了解决现有的骨折愈合动态过程模型不能综合模拟骨折愈合过程中骨痂形状的动态变化、力学环境变化、生物学环境变化间的复杂关系的缺点。一、三维几何模型的建立；二、进行网格的划分；三、建立骨及骨痂的生物力学模型；四、确定仿真初始参数及施加负载与边界条件；五、分为正常血供单元和非正常血供单元，如果是正常血供区域则执行六，如果是非正常血供区域则执行七；六、建立正常血供区域的确定性数学模型；七、建立非正常血供区域的模糊数学模型；八、根据六和七建立骨折愈合的仿真过程。本发明用于生物医学工程领域。

Description

一种基于生物学机理的骨折愈合动态过程仿真方法

技术领域

本发明涉及基于生物学机理的骨折愈合动态过程仿真方法。

背景技术

事实上，不是所有的骨折都可以被修复，有时会有不修复或拖延修复，骨折延迟愈合或不愈合引起患肢疼痛，功能障碍，导致患者失业，由于骨折事故每年发生的总数较大，骨折延迟愈合或不愈合的人数就很可观，由此将带来很大的社会经济负担。骨折延迟愈合或不愈合是由于受到特定的几何因素、力学因素、生物学因素影响，因此关于骨折愈合过程及其影响骨折愈合速度和质量的研究一直备受关注，也取得了一些收获。但是受到研究手段和骨折过程复杂又不能直接观察的限制，尽管在该领域的研究一直在进步，但仍有约5％～10％的骨折因各种原因发生延迟愈合或不愈合。

目前缺少能够精确表达骨折愈合这一复杂过程的计算机仿真模型，首先是没有建立起力学因素与骨折过程间的确定性关系式；没有在同一个仿真模型中充分考虑到力学环境和生物学环境的双重影响；没有从建立针对专门患者的个体化模型角度解决仿真建模问题，几何模型和生物力学材料设置过于简化，影响模型的仿真结果；没有综合考虑骨折愈合过程中骨痂形状的动态变化、力学环境变化、生物学环境变化间的复杂关系在计算机中的动态仿真。

发明内容

本发明的目的是为了解决现有的骨折愈合动态过程模型不能综合模拟骨折愈合过程中骨痂形状的动态变化、力学环境变化、生物学环境变化间的复杂关系的缺点，而提出一种基于生物学机理的骨折愈合动态过程仿真方法。

一种基于生物学机理的骨折愈合动态过程仿真方法按以下步骤实现：

步骤一、三维几何模型的建立；

步骤二、将得到的三维几何模型导入到网格划分软件中进行网格的划分；

步骤三、在网格的划分的基础上建立骨及骨痂的生物力学模型；

步骤四、在骨及骨痂的生物力学模型基础上，确定仿真初始参数及施加负载与边界条件；

步骤五、在步骤四的基础上进行骨痂形状动态变化的仿真设计，对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；

血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域，如果是正常血供区域则执行步骤六，如果是非正常血供区域则执行步骤七；

步骤六、建立正常血供区域在不同初始骨内应力环境下，骨密度、软骨密度随时间变化的确定性数学模型；

步骤七、建立非正常血供区域的模糊数学模型，非正常血供区域的模糊数学模型包括隶属度和模糊控制规则；

步骤八、根据步骤六和步骤七建立骨折愈合的仿真过程。

本发明的有益效果为：

本发明研究目标是基于生物学机理实现骨折愈合动态过程建模与仿真，仿真***具备复杂性、动态特性、个体化特性、可靠性。通过建立初始骨折内力学参数与骨折愈合过程刚度值之间的确定数学关系式，实现能反映力学环境影响因素的仿真模型；通过建立力学参数、血管血供参数、骨组织特性参数三者间关系的模糊准则，实现能反映生物学环境影响的仿真模型；通过有限元方法实现动态力学计算，通过模糊逻辑控制器实现非正常血供区域骨愈合过程的动态更新，通过应力分配的计算和临界应力的选择实现骨痂形状的动态更新。

发明的用途：

(1)通过构建骨折愈合动态过程计算机仿真***，该***能作为一个仿真试验平台，用于各种适应人体生物力学要求的内固定器材的评价和优化设计；用于非生物学因素对骨组织愈合速度和质量影响的测试分析；进行寻找有效避免骨折发生方法的研究；由于该仿真模型的可扩展性，该平台亦可实现对部分影响骨愈合进程的生物学因素进行分析和研究。

(2)通过构建骨折愈合动态过程计算机仿真***，能够实现复杂骨折愈合过程的预测，对医生制订正确的手术方案提供指导，进而提高手术成功率、提高骨折愈合质量，可以有效减少骨折不愈合和延迟愈合的情况，帮助骨折愈合者通过有效的治疗手段恢复健康生活，减轻由此带来的社会经济负担。

(3)可以利用计算机中建立的仿真模型进行多次重复试验研究，不需要真实的生物学试验，节省时间，提高效率，节省费用，避免人道主义的争议。

发明的特色之处：

(1)传统关于力学环境对骨折愈合进程的影响方面的研究主要集中于通过力学试验测试不同的力学环境对骨折愈合速度和最后骨痂刚度的影响，大部分是定性研究，没有从定量的角度以取得确定数学模型为目的探寻这个问题，本发明依据生物数学思想，建立了骨折内应力与骨折刚度的确定性数学模型。

(2)传统骨愈合过程仿真模型中很少考虑血供状态的影响，更没有把血供重建参数作为动态变量加入到仿真模型中的研究，在血供环境不好的区域，骨痂形成将会受到影响，力学环境不稳定会限制血管的重新生成，骨愈合过程是受到力学环境和生物学环境等多种因素影响的复杂过程，本发明通过模糊控制规则描述了骨折愈合的复杂过程，将血供重建作为动态变量融入仿真模型中。

(3)传统骨折愈合仿真的研究工作保持在一个静态或者只体现较少动态变量的研究水平，本发明通过引入临界应力的概念，将应力特性与骨痂生成的生物学特性联系起来，通过骨折内应力与骨折刚度的确定性数学模型和血供重建的模糊控制规则，最终实现能够同时表现骨痂形态变化、骨痂材料特性变化与血供重建过程的动态骨愈合过程计算机仿真。

通过对比仿真结果和实验数据可以得出：基于生物学机理的骨折愈合动态过程仿真***能通过有效计算不断更新骨折愈合过程中各参数的变化，准确模拟骨折愈合。

附图说明

图1a为膨胀应变隶属度函数示意图；

图1b为畸变应变隶属度函数示意图；

图1c为血液、软骨、骨密度隶属度函数示意图；

图1d为相邻单元隶属度函数示意图；

图1e为血液改变量隶属度函数示意图；

图1f为软骨密度隶属度函数示意图；

图1g为骨密度隶属度函数示意图；

图2为本发明骨折愈合动态过程仿真路线图；

图3a为本发明组织转变过程图；

图3b为本发明血管组织转变过程图；

图4为本发明外加载荷和边界条件示意图；

图5a随时间变化的折断间位移(即骨间动度))2mm稳定组仿真结果与实验数据的对此示意图；

图5b随时间变化的折断间位移(即骨间动度))3mm不稳定组仿真结果与实验数据的对此示意图。

具体实施方式

具体实施方式一：本实施方式的一种基于生物学机理的骨折愈合动态过程仿真方法具体是按照以下步骤进行的：

步骤一、三维几何模型的建立；

步骤二、将得到的三维几何模型导入到网格划分软件中进行网格的划分；

步骤三、在网格的划分的基础上建立骨及骨痂的生物力学模型；

步骤四、在骨及骨痂的生物力学模型基础上，确定仿真初始参数及施加负载与边界条件；

步骤五、在步骤四的基础上进行骨痂形状动态变化的仿真设计；

步骤六、建立正常血供区域在不同初始骨内应力环境下，骨密度、软骨密度随时间变化的确定性数学模型；

步骤七、建立非正常血供区域的模糊数学模型，非正常血供区域的模糊数学模型包括隶属度和模糊控制规则；

步骤八、根据步骤六和步骤七建立骨折愈合的仿真过程。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中三维几何模型的建立；具体过程为：

采用基于分割的三维医学图像表面重建算法对图像进行三维表面重构，得到三维几何模型；

所述图像由影像设备CT得到，数据存储格式为DICOM。

由表面模型构建实体模型的过程，就是由表面模型的三角片序列及上下底面构建实体模型的面链表，由表面模型的顶点序列构建实体模型的顶点链表，同时建立起体、环、边、半边链表及各链表中结点的指向关系。以边界模型表达的实体构建过程由一系列欧拉操作实现。基本的欧拉操作包括如下互逆的5对：MVFS，MEV，MEF，MEKR，KFMRH；KVFS，KEV，KEF，KEMR，MFKRH。其中M表示构造，K表示删除，S、E、V、F、R、H分别表示体、边、顶点、面、环、孔。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中将得到的三维几何模型导入到网格划分软件中进行网格的划分；具体过程为：

由于网格划分软件会生成许多数据，而本发明只需要节点坐标和单元编号，将得到的三维几何模型导入到MATLAB中进行预处理，只提取目标数据，根据目标数据生成后续有限元计算所需要的单元编号和节点坐标两个文件；

所述单元编号和节点坐标两个文件为txt文本格式的文件；

节点坐标文件包含三列数据，三列数据分别代表每个节点的空间坐标值；

单元编号文件包含四列数据，四列数据分别为每个单元的四个节点的节点序号。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三中在网格的划分的基础上建立骨及骨痂的生物力学模型；具体过程为：

针对骨和骨痂的力学特性，采用连续介质力学方法建立线弹性材料、非线性材料的本构方程，采用格林方法，通过骨及骨痂的势能函数来描述生物组织力学特性；

骨为线弹性材料，弹性势能W(X)是无限小应变张量的二次函数：

$W\left(X\right)=\frac{\mathrm{\λ}}{2}{\left({\mathrm{trE}}_L\right)}^2+{\mathrm{\μtrE}}_L^2$

式中，λ,μ是拉梅(Lame)常数，E_L为线性应变张量；

线弹性材料单元应力计算公式为：

{ε}＝[S]·{σ}，

式中，{σ}为线弹性材料单元应力，{ε}为线弹性材料单元应变，[S]为柔度矩阵；

$\[S\]=\left[\begin{array}{cccccc}\frac{1}{E_{12}}& -\frac{v_{12}}{E_{12}}& -\frac{v_{31}}{E_3}& & & \\ -\frac{v_{12}}{E_{12}}& \frac{1}{E_{12}}& -\frac{v_{31}}{E_3}& & 0& \\ -\frac{v_{13}}{E_{13}}& -\frac{v_{13}}{E_{12}}& \frac{1}{E_3}& & & \\ & & & \frac{1}{{\mathrm{\μ}}_{12}}& & \\ & 0& & & \frac{1}{{\mathrm{\μ}}_{13}}& \\ & & & & & \frac{1}{{\mathrm{\μ}}_{13}}\end{array}\right]$

式中，E₁₂、E₃为弹性模量，ν₁₃、ν₁₂为泊松比，μ₁₃为剪切模量，为材料工程常数，1,2,3分别代表骨骼的径向、切向和轴向；

骨痂为非线性弹性材料，弹性势能W表示为：

$W=C_1{e}^{\mathrm{\β}(I_1-3)}-\frac{C_1\mathrm{\β}}{2}(I_2-3)$

式中：I₁、I₂分别是应变张量的第一、第二主不变量，β、C₁分别是材料常数；

非线性弹性材料单元应力计算公式为：

dσ＝D_Tdε，

式中，dσ为非线性弹性材料单元应力，dε为非线性弹性材料单元应变，D_T为切线弹性矩阵；

切线弹性矩阵D_T为：

$D_T=\frac{{\∂}^{2}W}{\∂E^2}$

式中，E为格林应变张量。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤四中在骨及骨痂的生物力学模型基础上，确定仿真初始参数及施加负载与边界条件；具体过程为：

根据患者实际骨折状况设置各单元初始血供参数、初始组织特性参数，根据确定的骨折固定方式，添加负载和边界条件；

添加负载和边界条件有两种施加的方法，一是施加位移，二是施加力。

材料是刚性的或者边界材料比我们研究的材料更硬，可以通过在一系列顶点施加给定的位移来模拟这样的接触碰撞；如果边界材料的刚度小于或等于被研究材料刚度，可以通过施加线性弹簧力来建模。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述步骤五中在步骤四的基础上进行骨痂形状动态变化的仿真设计，对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域，如果是正常血供区域则执行步骤六，如果是非正常血供区域则执行步骤七；；具体过程为：

通过设置临界应力值，按照力学分布的状态，当单元应力大于等于临界应力值时，单元为需要保留的单元，当单元应力大于等于临界应力值时，单元为被组织吸收的单元；

对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；

血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域，如果是正常血供区域则执行步骤六，如果是非正常血供区域则执行步骤七；

临界应力用最大第三主应变ε_3,max与一个系数e相乘来计算，临界应力＝eε_3,max。

整个仿真过程的力学解算也是通过有限元思想来进行解算模型的设计；利用肉芽组织的低应力矢量定义骨痂的空间分配，由于软组织的存在，可以产生较大范围的骨折内运动，由此可能产生更大的骨痂，压缩应变强度最大允许量的比例被用来确定骨痂形状和尺寸。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：所述系数e的取值在0-1之间。

其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：所述步骤六中建立正常血供区域在不同初始骨内应力环境下，骨密度、软骨密度随时间变化的确定性数学模型；具体过程为：

初步分析骨折内应力与骨愈合速度之间的内在联系；统计分析中，选取符合要求的材料数据，建立愈合时间与骨密度、软骨密度变化的关系曲线；则正常血供区域在不同初始骨内应力环境下，骨密度、软骨密度随时间变化的确定性数学模型为：

f(t)＝a(t^m-1)+1

f′(t)＝a′(t^m′-1)+1

式中，a、a′为血供区域在不同初始骨内应力下的系数；m、m′为血供区域在不同初始骨内应力下的指数；f(t)为骨密度，f′(t)为软骨密度，t为时间。

材料数据的整理、统计及对照研究；

选择满足下列条件的研究资料进行分析：血管损伤小(充血为100％)的骨折或相应区域；试验过程中进行了骨密度、软骨密度和骨折内运动情况的记录；初始骨折内应力已知或通过步骤三计算获取；实验研究具有连续性，至少有7周的结果被统计；骨折方式和部位相近。对试验进行分类分析，对已有的试验数据或新设计的实验数据通过应力大小、固定方式、测量方法、被测量、缝隙进行分类；

实验材料一般是记录下来8周内位移即骨折内运动变化情况，通过建立骨折内运动与骨密度、软骨密度之间的变化关系，将试验材料统一起来，同时确定好统一的单位和衡量标准，这些是进行数据对比分析的基础，对各组实验数据分别分析找出规律和特点。

建立经验公式的主要步骤为：a.根据整理的数据表，选取合适的坐标，绘出数据折线或散点图；b.根据折线或散点图的形状判断函数关系，建立回归方程；c.根据回归方程，求出待定常数，并进行相关性检验。

其它步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式与具体实施方式一至八之一不同的是：所述步骤七中建立非正常血供区域的模糊数学模型，非正常血供区域的模糊数学模型包括隶属度和模糊控制规则；具体过程为：

步骤七一、建立隶属度函数；

用充血、软骨密度、骨密度三个组织浓度变量，相邻单元的充血、相邻单元的骨密度，膨胀应变、畸变应变两个力学刺激作为模糊输入，经过21条模糊控制规则描述组织分化过程；用充血、骨密度、软骨密度的改变量作为输出；首先建立7个输入3个输出语言变量的语言值，然后以组织学实验结果和细胞培养实验结果为基础建立如图1a、图1b、图1c、图1d、图1e、图1f、图1g所示的隶属度函数；

设置7个输入3个输出语言变量的语言值，建立隶属度函数；

7个输入为充血、软骨密度、骨密度、相邻单元的充血、相邻单元的骨密度、膨胀应变、畸变应变；

3个输出为充血的改变量、骨密度的改变量、软骨密度的改变量；

步骤七二、建立模糊控制规则；

模糊控制规则主要描述愈合进程中组织转变过程，这些规则是模糊模型的基础，本研究采用if A and B then C的语句描述，其主要形式如表1所示。

规则1-4为血管重建，具体过称为：

血肿的刺激使得骨折处的毛细血管再生，毛细血管再生是从已存在的血管中产生新血管的过程，并在力学的调控下产生生物学效应而诱导血管重建；

规则5和规则6为膜内骨化，具体过称为：

膜内骨化是一种骨生长的成骨方式，尤其是在成骨过程中没有软骨生成时，膜内骨化是***转变成骨组织的直接转化形式；初期，骨髓间充质干细胞开始增殖聚集，血管生成活性增加，毛细血管开始增长；骨髓间充质干细胞进而发展成为成骨细胞，开始产生类骨质；同时，类骨质开始以相同的速率矿化；在骨表面，由成骨细胞促使生成新的类骨质；而成骨细胞逐渐变成骨细胞。这种成骨方式仅仅在稳定的力学环境和充足的血液供应条件下，并且相邻单元具有较高骨密度时才会发生。

规则7-9为软骨的形成，具体过称为：

在固定不良，骨折断端活动较大，同时伴随血液供给不足的情况下，断端周围组织内的间充质细胞才会分化成软骨细胞，继而产生软骨。

规则10-13为软骨钙化过程，具体过称为：

软骨骨痂形成后，成骨细胞在高应力刺激下合成并分泌出溶胶原、蛋白多糖和糖蛋白，构成骨的有机质；

溶胶原逐渐聚合成胶原纤维，并和成骨细胞一起被基质包围，逐渐发生钙化；这个过程的发生与血液供给无关。

规则14、15为软骨骨化过程，具体过称为：

在软骨钙化区内，营养物质逐渐弥散发生障碍，软骨细胞逐渐凋亡，钙化基质逐渐降解，释放血管生长因子，此时，大量的毛细血管和成骨细胞侵入，矿化的软骨基质逐渐产生类骨质，发生骨化现象；因此，在血液供给良好的条件下，完全钙化的软骨才能发生骨化；在本模型中，完全钙化的软骨是以高的骨密度和低的软骨密度表达的。规则不仅预测了骨密度的增加而且也预测了相同数量软骨密度的减少。

除了正常状态下的骨痂分化过程，规则16-21为在膨胀应变为4％～8％，-4％～-8％，畸变应变14％～20％的加载条件下引起的骨组织萎缩或者没产生诱发骨组织分化条件时保持的血供、骨密度、软骨密度不变；

模糊规则主要描述的组织转变过程如图2所示；

其它步骤及参数与具体实施方式一至八之一相同。

具体实施方式十：本实施方式与具体实施方式一至九之一不同的是：所述步骤八中根据步骤六和步骤七建立骨折愈合的仿真过程；具体过称为：

有限元模型被建立以后，对有限元模型进行模型的解算，得出单元应力，当单元应力大于等于临界应力值时，单元为需要保留的单元，当单元应力大于等于临界应力值时，单元为被组织吸收的单元；

对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；

血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域；

如果是正常血供单元，通过正常血供区域骨密度和软骨密度随时间变化的确定性数学模型更新其状态；

如果是非正常血供单元，通过模糊数学模型确定血供信息的变化和骨组织特性的更新；

完成后再重复步骤一至步骤八；直到最后的骨痂生成结束不再有可吸收的单元；(通过删除相应网格单元的边及多余的节点)；具体实施步骤如图3a、图3b所示：

所述有限元模型包括三维几何模型、网格划分、生物力学模型、负载和边界条件、仿真初始参数。

采用以下实施例验证本发明的有益效果：

实施例一：

本实施例一种基于生物学机理的骨折愈合动态过程仿真方法具体是按照以下步骤制备的：

模拟绵羊跖骨骨折的愈合，在模型顶端施加载荷大小为F＝500N，假设这是在绵羊正常行走中足够满足跖骨可能承受的最大的力。底部为固定约束部分，其自由度均为0，包括三个方向的位移和三个方向的旋转。外部固定架的安放使得模型只能发生轴向位移，其它方向自由度均为0。图4为三维外固定横向截骨模型的外加载荷和固定约束的示意图。骨的弹性模量为4000，泊松比为0.36，软骨的弹性模量为200，泊松比为0.45，表2为初始血供参数和材料性能参数选取；图5a和图5b分别为2mm间隙稳定情况和3mm间隙不稳定情况的模拟结果。

表2初始血供参数及材料性能参数选取

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (10)

1.一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：一种基于生物学机理的骨折愈合动态过程仿真方法具体是按照以下步骤进行的：

步骤一、三维几何模型的建立；

步骤二、将得到的三维几何模型导入到网格划分软件中进行网格的划分；

步骤三、在网格的划分的基础上建立骨及骨痂的生物力学模型；

步骤四、在骨及骨痂的生物力学模型基础上，确定仿真初始参数及施加负载与边界条件；

步骤五、在步骤四的基础上进行骨痂形状动态变化的仿真设计，对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；

血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域，如果是正常血供区域则执行步骤六，如果是非正常血供区域则执行步骤七；

步骤六、建立正常血供区域在不同初始骨内应力环境下，骨密度、软骨密度随时间变化的确定性数学模型；

步骤七、建立非正常血供区域的模糊数学模型，非正常血供区域的模糊数学模型包括隶属度和模糊控制规则；

步骤八、根据步骤六和步骤七建立骨折愈合的仿真过程。

2.根据权利要求1所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤一中三维几何模型的建立；具体过程为：

采用基于分割的三维医学图像表面重建算法对图像进行三维表面重构，得到三维几何模型；

所述图像由影像设备CT得到，数据存储格式为DICOM。

3.根据权利要求2所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤二中将得到的三维几何模型导入到网格划分软件中进行网格的划分；具体过程为：

将得到的三维几何模型导入到MATLAB中进行预处理，只提取目标数据，根据目标数据生成后续有限元计算所需要的单元编号和节点坐标两个文件；

所述单元编号和节点坐标两个文件为txt文本格式的文件；

节点坐标文件包含三列数据，三列数据分别代表每个节点的空间坐标值；

单元编号文件包含四列数据，四列数据分别为每个单元的四个节点的节点序号。

4.根据权利要求3所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤三中在网格的划分的基础上建立骨及骨痂的生物力学模型；具体过程为：

采用连续介质力学方法建立线弹性材料、非线性材料的本构方程，采用格林方法，通过骨及骨痂的势能函数来描述生物组织力学特性；

骨为线弹性材料，弹性势能W(X)是无限小应变张量的二次函数：

$W\left(X\right)=\frac{\mathrm{\λ}}{2}{\left({\mathrm{trE}}_L\right)}^2+{\mathrm{\μtrE}}_L^2$

式中，λ,μ是拉梅(Lame)常数，E_L为线性应变张量；

线弹性材料单元应力计算公式为：

{ε}＝[S]·{σ}，

式中，{σ}为线弹性材料单元应力，{ε}为线弹性材料单元应变，[S]为柔度矩阵；

$\[S\]=\left[\begin{array}{cccccc}\frac{1}{E_{12}}& -\frac{v_{12}}{E_{12}}& -\frac{v_{31}}{E_3}& & & \\ -\frac{v_{12}}{E_{12}}& \frac{1}{E_{12}}& -\frac{v_{31}}{E_3}& & 0& \\ -\frac{v_{13}}{E_{12}}& -\frac{v_{13}}{E_{12}}& \frac{1}{E_3}& & & \\ & & & \frac{1}{{\mathrm{\μ}}_{12}}& & \\ & 0& & & \frac{1}{{\mathrm{\μ}}_{13}}& \\ & & & & & \frac{1}{{\mathrm{\μ}}_{13}}\end{array}\right]$

式中，E₁₂、E₃为弹性模量，ν₁₃、ν₁₂为泊松比，μ₁₃为剪切模量，为材料工程常数，1,2,3分别代表骨骼的径向、切向和轴向；

骨痂为非线性弹性材料，弹性势能W表示为：

$W=C_1{e}^{\mathrm{\β}(I_1-3)}-\frac{C_1\mathrm{\β}}{2}(I_2-3)$

式中：I₁、I₂分别是应变张量的第一、第二主不变量，β、C₁分别是材料常数；

非线性弹性材料单元应力计算公式为：

dσ＝D_Tdε，

式中，dσ为非线性弹性材料单元应力，dε为非线性弹性材料单元应变，D_T为切线弹性矩阵；

切线弹性矩阵D_T为：

$D_T=\frac{{\∂}^{2}W}{\∂E^2}$

式中，E为格林应变张量。

5.根据权利要求4所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤四中在骨及骨痂的生物力学模型基础上，确定仿真初始参数及施加负载与边界条件；具体过程为：

根据患者实际骨折状况设置各单元初始血供参数、初始组织特性参数，根据确定的骨折固定方式，添加负载和边界条件；

添加负载和边界条件有两种施加的方法，一是施加位移，二是施加力。

6.根据权利要求5所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤五中在步骤四的基础上进行骨痂形状动态变化的仿真设计，对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域，如果是正常血供区域则执行步骤六，如果是非正常血供区域则执行步骤七；；具体过程为：

通过设置临界应力值，当单元应力大于等于临界应力值时，单元为需要保留的单元，当单元应力大于等于临界应力值时，单元为被组织吸收的单元；

对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；

血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域，如果是正常血供区域则执行步骤六，如果是非正常血供区域则执行步骤七；

临界应力＝eε_3,max

式中，ε_3,max为最大第三主应变，e为系数。

7.根据权利要求6所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述系数e的取值在0-1之间。

8.根据权利要求7所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤六中建立正常血供区域在不同初始骨内应力环境下，骨密度、软骨密度随时间变化的确定性数学模型；具体过程为：

建立愈合时间与骨密度、软骨密度变化的关系曲线；则正常血供区域在不同初始骨内应力环境下，骨密度、软骨密度随时间变化的确定性数学模型为：

f(t)＝a(t^m-1)+1

f′(t)＝a′(t^m′-1)+1

式中，a、a′为血供区域在不同初始骨内应力下的系数；m、m′为血供区域在不同初始骨内应力下的指数；f(t)为骨密度，f′(t)为软骨密度，t为时间。

9.根据权利要求8所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤七中建立非正常血供区域的模糊数学模型，非正常血供区域的模糊数学模型包括隶属度和模糊控制规则；具体过程为：

步骤七一、建立隶属度函数；

设置7个输入3个输出语言变量的语言值，建立隶属度函数；

7个输入为充血、软骨密度、骨密度、相邻单元的充血、相邻单元的骨密度、膨胀应变、畸变应变；

3个输出为充血的改变量、骨密度的改变量、软骨密度的改变量；

步骤七二、建立模糊控制规则；

规则1-4为血管重建，具体过称为：

血肿的刺激使得骨折处的毛细血管再生，毛细血管再生是从已存在的血管中产生新血管的过程，并在力学的调控下产生生物学效应而诱导血管重建；

规则5和规则6为膜内骨化，具体过称为：

膜内骨化是***转变成骨组织的直接转化形式；初期，骨髓间充质干细胞开始增殖聚集，血管生成活性增加，毛细血管开始增长；骨髓间充质干细胞进而发展成为成骨细胞，开始产生类骨质；同时，类骨质开始以相同的速率矿化；在骨表面，由成骨细胞生成新的类骨质；而成骨细胞变成骨细胞；

规则7-9为软骨的形成，具体过称为：

在骨折断端活动，血液供给不足的情况下，断端组织内的间充质细胞会分化成软骨细胞，继而产生软骨；

规则10-13为软骨钙化过程，具体过称为：

软骨骨痂形成后，成骨细胞在应力刺激下合成并分泌出溶胶原、蛋白多糖和糖蛋白，构成骨的有机质；

溶胶原聚合成胶原纤维，并和成骨细胞被基质包围，发生钙化；

规则14、15为软骨骨化过程，具体过称为：

在软骨钙化区内，营养物质发生障碍，软骨细胞凋亡，钙化基质降解，释放血管生长因子，此时，毛细血管和成骨细胞侵入，矿化的软骨基质产生类骨质，发生骨化现象；因此，在血液供给的条件下，完全钙化的软骨才能发生骨化；

规则16-21为在膨胀应变为4％～8％，-4％～-8％，畸变应变14％～20％的加载条件下引起的骨组织萎缩或者没产生诱发骨组织分化条件时保持的血供、骨密度、软骨密度不变。

10.根据权利要求9所述一种基于生物学机理的骨折愈合动态过程仿真方法，其特征在于：所述步骤八中根据步骤六和步骤七建立骨折愈合的仿真过程；具体过称为：

有限元模型被建立以后，对有限元模型进行模型的解算，得出单元应力，当单元应力大于等于临界应力值时，单元为需要保留的单元，当单元应力大于等于临界应力值时，单元为被组织吸收的单元；

对于保留的单元，依据血供参数值区分为正常血供单元和非正常血供单元；

血供参数值100％为正常血供区域，血供参数值小于100％为非正常血供区域；

如果是正常血供单元，通过正常血供区域骨密度和软骨密度随时间变化的确定性数学模型更新其状态；

如果是非正常血供单元，通过模糊数学模型确定血供信息的变化和骨组织特性的更新；

完成后再重复步骤一至步骤八；直到最后的骨痂生成结束不再有可吸收的单元；

所述有限元模型包括三维几何模型、网格划分、生物力学模型、负载和边界条件、仿真初始参数。