CN106777582A - 一种基于组织分化的长骨骨折愈合仿真系统 - Google Patents

Abstract

一种基于组织分化的长骨骨折愈合仿真系统，设计生物医学工程领域。本发明用来预测骨折愈合的复杂过程，探寻最佳的骨折愈合方案。所述系统包括骨折区域几何建模模块、骨折区域生物力学有限元分析模块、骨痂单元组织分化模块和程序终止判断模块。骨折区域几何建模模块用于建立骨折区域三维几何模型；骨折区域生物力学有限元建模模块用于对建立的三维几何模型进行有限元分析，得到单元力学刺激；单元组织分化模块用于仿真组织分化，使单元内各组织含量得到更新；程序终止判断模块用于判断程序是否终止。本发明将骨折区域看作双相多孔弹性模型，能够更加准确的模拟骨折愈合过程，为探寻最佳的骨折愈合方案提供有益帮助。

Description

一种基于组织分化的长骨骨折愈合仿真系统

技术领域

本发明涉及生物医学工程领域，特别涉及一种基于组织分化的长骨骨折愈合仿真系统。

背景技术

骨折是一种常见的创伤，骨折的高发性使得对骨折机理及促进愈合的研究尤为迫切，一旦骨折发生，与其它组织损伤修复不同的是，骨折不是靠纤维结缔组织连接，而是骨组织的完全再生。然而，并不是所有的骨折都可以完成愈合，有时会发生延迟愈合甚至是不愈合。骨折延迟愈合或者不愈合会引起患肢疼痛，功能障碍，导致患者失业，由此造成很大的社会经济负担。因此，尽管关于骨折愈合的研究一直备受关注，但仍有5％～10％的骨折因各种原因发生延迟愈合甚至是不愈合。

骨折愈合受到特定的几何因素、力学因素、生物学因素影响，良好的几何因素、力学因素、生物学因素得到良好愈合效果。反之则会导致骨折的延迟愈合甚至是不愈合。目前缺少能够精确表达骨折愈合这一复杂过程的计算机仿真系统。在骨折愈合仿真系统中存在以下缺陷：

1.没有建立专门针对患者的个体化模型；

2.力学因素与骨折愈合过程没有一个确定性关系；

3.骨折愈合区域的模型和生物力学材料设置过于简化；

4.没有在同一个仿真系统中综合考虑力学因素和生物学因素对骨折愈合的影响。

发明内容

本发明的目的是为了解决现有的骨折愈合仿真中不能综合模拟力学因素和生物学因素对骨折愈合过程的影响，骨折愈合生物力学模型材料设置过于简化的缺点，而提出的一种基于组织分化的长骨骨折愈合仿真系统。

本发明的目的通过下述技术方案实现：一种基于组织分化的长骨骨折愈合仿真系统，其特征在于，所述系统包括：

骨折区域几何建模模块、骨折区域生物力学有限元分析模块、骨痂单元组织分化模块和程序终止判断模块；

骨折区域几何建模模块用于根据导入的二维断层扫描图像数据，经过图像预处理后进行骨折部位的三维表面几何模型的建立；

骨折区域生物力学有限元分析模块用于对建立好的骨折区域模型进行网格划分，施加外部载荷和设置边界条件；

骨折区域生物力学有限元分析模块还用于初始骨折区域环境的设置；

骨折区域生物力学有限元分析模块还用于计算单元力学刺激；

骨痂单元组织分化模块用于对单元内组织分化进行仿真，使单元内各组织含量得到更新，从而使单元材料属性得到更新，从而得到下一迭代步中所需要的单元力学刺激；

程序终止判断模块用于判断程序是否终止，若不满足终止条件，程序进行下一迭代步；若满足程序终止条件，则程序结束并输出愈合时间。

本发明的有益效果为：

1.本发明提出的一种基于组织分化的长骨骨折愈合仿真系统是基于windows开发语言平台来开发软件，通过自主编程实现骨折愈合过程的动态模拟，基于对话框的形式，易于操作，培训周期短；

2.将骨折区域设置为双相多孔弹性模型，相比单相模型，更加符合骨折愈合区域的生物特性，使仿真结果更加精确；

3.将血供作为变量引入骨折愈合仿真系统中，能够更加准确的模拟骨折愈合的过程以及血液对骨折愈合的影响；

4.利用模糊逻辑的方法对组织分化进行仿真，相较于传统的采用偏微分方程组进行建模的方法，更加便于理解。减少了在建模过程中偏微分方程组的数量，便于建模且减少了仿真时间；

5.通过构建骨折愈合仿真系统，可以对医生制定最优的手术方案提供指导，进而提高手术成功率、提高骨折愈质量，减少骨折不愈合和延迟愈合的情况；

6.通过构建骨折愈合仿真系统，可以对建立的仿真模型进行多次重复实验研究，减少真实的生物实验，节省时间，提高效率，节省费用，避免人道主义的争议。

综上，本发明的仿真平台克服了现有技术的缺点与不足。

附图说明

图1为基于组织分化的骨折愈合仿真系统流程图；

图2为骨折区域几何模型建立流程图；

图3为骨折区域生物力学有限元分析流程图；

图4为骨痂单元组织分化模糊控制示意图；

图5为组织分化与骨痂单元力学刺激之间的关系示意图；

图6为骨痂单元力学刺激隶属度函数；

图7为软骨组织含量隶属度函数；

图8为骨组织含量隶属函数；

图9为血供隶属度函数。

具体实施方式

具体实施方式一：如图1所示，本实施方式所述的一种基于组织分化的长骨骨折愈合仿真系统包括：

骨折区域几何建模模块1、骨折区域生物力学有限元分析模块2、骨痂单元组织分化模块3、程序终止判断模块4；

骨折区域几何建模模块1用于根据导入的二维断层扫描图像数据，经过图像预处理后进行骨折部位的三维表面几何模型的建立；

骨折区域生物力学有限元分析模块2用于对建立好的骨折区域模型进行网格划分，施加外部载荷和设置边界条件；

骨折区域生物力学有限元分析模块2还用于初始骨折区域环境的设置；

骨折区域生物力学有限元分析模块2还用于计算单元力学刺激；

骨痂单元组织分化模块3用于对单元内组织分化进行仿真，使单元内各组织含量得到更新，从而使单元材料属性得到更新，从而得到下一迭代步中所需要的单元力学刺激；

程序终止判断模块4用于判断程序是否终止，若不满足终止条件，程序进行下一迭代步；若满足程序终止条件，则程序结束并输出愈合时间。

具体实施方式二：如图1～9所示，本实施方式中，所述的骨折区域几何建模模块1实现其功能的具体过程为：

采用基于分割的三维医学影像表面重建算法对图像进行三维表面重构，通过阈值筛选、交互式分割和三维重建过程得到三维表面几何模型；

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

其他组成及连接与具体实施方式一相同。

具体实施方式三：如图1～9所示，本实施方式中，所述的骨折区域生物力学有限元分析模块2实现其功能的具体过程为：

1)将骨折区域三维表面几何模型进行网格化分，使连续的几何模型离散化，得到骨折区域有限元模型；

所述的网格划分包括面网格划分和体网格划分两个步骤；面网格划分过程用于将三维表面模型进行优化，包括：表面模型优化，平滑处理，修补漏洞；表面模型的优化通过减小表面模型的三角面片来实现，该过程只需将相邻的两个顶点合并到一个新的顶点上，并延续原有的拓扑关系；平滑处理的过程中，对三维的面网格模型进行去噪；修补漏洞的过程中，通过将模型当中的空洞提取成空间多边形，然后对空洞多边形进行三角化的方法实现；体网格划分的过程是将面网格模型进行拉伸、旋转步骤来实现的；

通过网格划分得到的骨折区域有限元模型包括单元编号和节点坐标两部分；

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

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

2)在有限元模型上施加外加载荷，并设置边界条件。载荷的大小由骨所承受的力的大小决定，实验对象不同，所受的力也不同；

3)对骨折区域有限元模型进行骨折区域初始环境设置。骨折区域由皮质骨和骨痂区域两部分组成。初始骨折区域环境设置包括皮质骨材料属性赋值，皮质骨血供赋值，初始骨痂材料属性赋值，初始骨痂区域血供赋值；

骨折之后，骨折断端处的皮质骨收到损害，血供也遭到破坏，所以将距离骨折断端5mm内的皮质骨血供设置为0％，其余部分皮质骨结构完整，血供条件良好，所以将其余部分的皮质骨血供设置为100％；

骨痂外周区域可以接受其周围组织提供的血供，所以将骨痂外周3mm内血供设置为30％，骨痂内部血供设置为0％。

4)将骨折区域看作双相多孔弹性模型，由多孔弹性理论得到骨痂单元的本构方程，平衡方程和几何方程，并通过有限单元法计算骨痂单元应力刺激S，具体过程为：

a.本构方程

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；α，α'分别为各向同性弹性面的Biot系数和轴向Biot系数；p为骨痂单元中的流体压力；M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅为脱水的弹性模量矩阵分量；

其中，M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅脱水的弹性模量矩阵分量表达式如下所示：

M₄₄＝E_r/2(1+ν_r) (6)

M₅₅＝G' (7)

式中，E_r，ν_r分别是各同性弹性层的弹性模量和泊松比；E_z，ν_z分别是轴向弹性模量和泊松比；G'为剪切模量；

b.平衡方程

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；r为径向半径；

c.几何方程

式中，ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；u_r，u_θ，u_z分别为三个方向上的位移；r为径向半径；

通过上述方程的求解得到骨痂单元的正应变σ_rr，σ_θθ，σ_zz，由正应变可得到骨痂单元受到的畸变应变：

式中，D为骨痂单元受到的畸变应变；σ_rr，σ_θθ，σ_zz分别为各个方向上的正应变；

骨痂单元中液体的流速V为：

其中，k为骨痂中液体的达西渗透系数；u液体粘度；p为液体压力；

由此可得到骨痂单元所受的力学刺激S为：

其中，D为骨痂单元受到的畸变应变；V为骨痂单元中液体流速；a，b分别为经验常数。

其他组成及连接与具体实施方式一或二之一相同。

具体实施方式四：如图1～9所示，本实施方式中，所述的骨痂单元组织分化模块3实现其功能的具体过程为：

1)设置输入变量隶属度函数

骨痂单元组织分化过程由六个输入变量协同决定，分别为单元力学刺激，骨痂单元骨含量，骨痂单元软骨含量，骨痂单元血供，周围骨痂单元骨含量，周围骨痂单元血供。输入变量的精确值通过隶属度函数转变为相应输入变量的隶属度。输入变量隶属度函数由梯形函数表示；

2)输入隶属度函数模糊化

将输入变量的隶属度函数用模糊逻辑语言值来表示：

单元力学刺激可分为五个等级，分别为：低，中，高；

单元血供可分为三个等级，分别为：低，中，高；

单元骨含量可分为三个等级，分别为：低，中，高；

单元软骨含量可分为三个等级，分别为：低，中，高。

3)设置输出变量隶属度函数

输出变量共有三个，分别为骨痂单元血供改变量，骨痂单元骨含量改变量和骨痂单元软骨含量改变量。设置输出变量隶属度函数，输出隶属度函数由梯形函数表示，且将输出隶属度函数用如下的模糊逻辑语言值表示：

骨痂单元血供改变量分为三个等级，分别为：低，中，高；

骨痂单元骨含量改变量分为三个等级，分别为：低，中，高；

骨痂单元软骨含量改变量分为三个等级，分别为：低，中，高。

4)模糊规则的编写

骨痂内组织分化主要包括：血管再生，膜内骨化，软骨生成，软骨钙化，软骨骨化五个过程。通过模糊规则分别表示这五个组织分化过程，由此可得到组织分化后，骨痂单元中骨的模糊逻辑语言值，软骨的模糊逻辑语言值以及血供的模糊逻辑语言值。具体过程如下：

过程一、血管再生

如果单元应力刺激为非高且骨痂单元血供为低且周围骨痂单元血供为非低，那么骨痂单元血供增加；

如果单元应力刺激为非高且骨痂单元血供为中且周围骨痂单元血供为高，那么骨痂单元血供增加；

如果单元应力刺激为非高且骨痂单元血供为高，那么骨痂单元血供增加；

过程二、膜内骨化

如果骨痂单元软骨含量为低且单元应力刺激为中且骨痂单元血供为高且周围骨痂单元血供为高，那么骨痂单元骨含量增加；

过程三、软骨生成

如果骨痂单元骨含量为非高且骨痂单元软骨含量为低且单元应力刺激为中，那么骨痂单元软骨含量增加；

如果骨痂单元软骨含量为非低且单元应力刺激为高，那么骨痂单元软骨含量增加；

如果骨痂单元软骨含量为高且单元应力刺激为中，那么骨痂单元软骨含量增加；

过程四、软骨钙化

如果骨痂单元软骨含量为非低且单元应力刺激为中且骨痂单元血供为非低且周围骨痂单元骨含量为非低，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

如果骨痂单元软骨含量为非低且单元应力刺激为高且骨痂单元血供为非低且周围骨痂单元骨含量为非低，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

过程五、软骨骨化

如果骨痂单元骨含量为高且骨痂单元软骨含量为低且单元应力刺激为中且骨痂单元血供为非低且周围骨痂单元骨含量为高，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

如果骨痂单元骨含量为高且骨痂单元软骨含量为低且单元应力刺激为低且骨痂单元血供为非低且周围骨痂单元骨含量为高，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

5)输出变量反模糊化

通过模糊规则的作用，得到输出变量的模糊逻辑语言值，利用面积中心法得到输出变量的隶属度；

6)更新输出变量含量

将得到的输出变量隶属度经过转化得到输出变量的改变量，从而得到经过组织分化后，骨痂单元骨组织，软骨组织，结缔组织含量以及骨痂单元血供。具体过程如下：

a.计算经组织分化后骨痂单元骨组织含量

式中，为第n个骨痂单元在第k+1迭代步中骨组织的含量；为第n个骨痂单元在第k迭代步中骨的含量；为单位时间第n个骨痂单元在第k迭代步中骨的生成量；Δt为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中骨的改变量由下式得到：

式中，为单位时间第n个骨痂单元在第k迭代步中骨的生成量；为骨痂单元中骨含量的隶属度；T_Bone在骨痂单元骨的转化率；

b.计算经组织分化后骨痂单元软骨骨组织含量

式中，为第n个骨痂单元在第k+1迭代步中软骨的含量；为第n个骨痂单元在第k迭代步中软骨的含量；为单位时间第n个骨痂单元在第k迭代步中软骨的生成量；Δt为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中软骨的改变量由下式得到：

式中，为单位时间第n个骨痂单元在第k迭代步中软骨的生成量；为骨痂单元中软骨含量的隶属度；T_Cartilage在骨痂单元软骨的转化率；

c.计算经组织分化后骨痂单元血供

式中，为第n个骨痂单元在第k+1迭代步中结缔组织的含量；为第n个骨痂单元在第k迭代步中结缔组织的含量；为单位时间第n个骨痂单元在第k迭代步中结缔组织的生成量；Δt为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中软骨的改变量由下式得到：

式中，为单位时间第n个骨痂单元在第k迭代步中结缔组织的生成量；为骨痂单元中结缔组织含量的隶属度；T_Perfusio_n在骨痂单元结缔组织的转化率；

d.计算经组织分化后骨痂单元结缔组织含量

骨痂单元由骨组织，软骨组织，结缔组织三部分组成，三者之间的关系如下所示：

μ_Bone+μ_Cartilage+μ_ConnTissue＝1 (26)

式中，μ_Bone为骨痂单元骨组织含量；μ_Cartilage为骨痂单元软骨组织含量；μ_ConnTissue为骨痂单元结缔组织含量；

由此可得组织分化后骨痂单元结缔组织含量。

其他组成及连接关系与具体实施方式一至三之一相同。

具体实施方式五：如图1～9所示，本实施方式中，所述的程序终止判断模块4实现其功能的具体过程为：

1)判断骨痂单元材料属性

判断当前骨痂单元材料属性是否与骨的材料属性相同，若不同，则程序执行下面的2)步骤；若不等，程序执行3)步骤；

2)更新骨痂单元材料属性

若骨痂单元材料属性不等于骨的材料属性，则骨痂单元材料属性进行更新，进入下一个迭代步，骨痂材料属性更新公式如下：

式中，为第n个骨痂单元在第k+1迭代步中的弹性模量；E_Bone为骨的弹性模量；为第n个骨痂单元在第k+1迭代步中骨的含量；E_Cartilage为软骨的弹性模量；为第n个骨痂单元在k+1迭代步中的软骨的含量；E_ConnTissue为结缔组织的弹性模量；为第n个骨痂单元在k+1迭代步中结缔组织的含量；

式中，为第n个骨痂单元在第k+1迭代步中的泊松比；ν_Bone为骨的泊松比；为第n个骨痂单元在第k+1迭代步中骨的含量；ν_Cartilage为软骨的泊松比；为第n个骨痂单元在k+1迭代步中的软骨的含量；ν_ConnTissue为结缔组织的泊松比；为第n个骨痂单元在k+1迭代步中结缔组织的含量；

3)程序结束

若所有骨痂单元材料属性等于骨的材料属性，程序终止并输出愈合时间。

其他组成及连接关系与具体实施方式一至四之一相同。

实施例：

为了说明本系统的使用方法，下面具体举一个例子阐述本发明的操作过程。

模拟羊跖骨骨折愈合过程

1.羊跖骨骨折区域有限元模型的建立

1)建立几何模型

把CT图像数据利用三维医学图像表面重建算法对图像进行三维表面重构，得到羊跖骨骨折区域三维几何模型。

2)网格划分

上述的三维几何模型导入网格划分软件中进行网格化分，将得到的有限元模型导入到MTLAB中进行预处理，只提取目标数据，根据目标数据生成后续有限元计算所需要的单元编号和节点坐标文件。单元编号和节点坐标两个文件为txt文本格式的文件。单元编号文件包含四列数据，四列数据分别为每个单元的四个节点序号，节点坐标文件包含散三列数据，三列数据分别为每个节点的空间坐标值。

2.骨折区域初始环境设置

骨折区域初始环境设置包括皮质骨材料属性赋值，初始骨痂单元材料属性赋值。将骨折区域看成是基于混合物的双相多孔弹性模型，皮质骨和骨痂单元材料属性如表1所示。

骨折区域初始环境设置还包括皮质骨血供赋值和骨痂区域血供赋值。骨折之后，骨折断端处的皮质骨收到损害，血供也遭到破坏，所以将距离骨折断端5mm内的皮质骨血供设置为0％，其余部分皮质骨结构完整，血供条件良好，所以将其余部分的皮质骨血供设置为100％；骨痂外周区域可以接受其周围组织提供的血供，所以将骨痂外周3mm内血供设置为30％，骨痂内部血供设置为0％。

表1各组织材料属性

3.骨折区域生物力学有限元分析

对建立好的有限元模型施加边界条件和载荷。设定骨折区域下端固定约束长度，即将羊跖骨下端约束的那段所有节点的自由度均赋值为0，包括3个方向的位移和3个方向的旋转。由于羊在正常行走时，跖骨所受力的大小为500N，所以在骨折区域施加大小为500N的边界载荷。与求解羊跖骨相关的经验常数a和b分别为a＝0.0375μm/s，b＝3μm/s。

接下来对有限元模型进行生物力学有限元分析，得到骨折区域单元的畸变应变和流体速度，从而得到骨折区域单元收到的单元力学刺激。

4.组织分化过程

将由上一步得到的单元力学刺激，以及骨痂单元骨组织含量、周围骨痂单元骨组织含量、骨痂单元软骨组织含量、骨痂单元血供和周围骨痂单元血供作为输入变量，导入到由模糊规则中，得到经组织分化后骨痂单元骨组织含量、骨痂单元软骨组织含量和骨痂单元血供，如图2所示。其中，单元力学刺激和组织分化之间的关系如图3所示，骨组织、软骨组织和结缔组织之间的组织分化关系如图4所示，输入变量及输出变量的隶属度函数如图5、6、7、8所示。

5.判断程序终止条件

判断当前骨痂单元材料属性是否等于骨的材料属性，若不等，则利用更新当前骨痂单元材料属性，程序进入下一个迭代步；若相等，则程序结束，并输出骨折愈合时间。

Claims (5)

1.一种基于组织分化的长骨骨折愈合仿真系统，其特征在于，所述系统包括：

骨折区域几何建模模块(1)、骨折区域生物力学有限元分析模块(2)、骨痂单元组织分化模块(3)和程序终止判断模块(4)；

骨折区域几何建模模块(1)用于根据导入的二维断层扫描图像数据，经过图像预处理后进行骨折部位的三维表面几何模型的建立；

骨折区域生物力学有限元分析模块(2)用于对建立好的骨折区域模型进行网格划分，施加外部载荷和设置边界条件；

骨折区域生物力学有限元分析模块(2)还用于初始骨折区域环境的设置；

骨折区域生物力学有限元分析模块(2)还用于计算单元力学刺激；

骨痂单元组织分化模块(3)用于对单元内组织分化进行仿真，使单元内各组织含量得到更新，从而使单元材料属性得到更新，进而得到下一迭代步中所需要的单元力学刺激；

程序终止判断模块(4)用于判断程序是否终止，若不满足终止条件，程序进行下一迭代步；若满足程序终止条件，则程序结束并输出愈合时间。

2.根据权利要求书1所述的一种基于组织分化的长骨骨折愈合仿真系统，其特征在于：所述的骨折区域几何建模模块(1)实现其功能的具体过程为：

采用基于分割的三维医学影像表面重建算法对图像进行三维表面重构，通过阈值筛选、交互式分割和三维重建过程得到三维表面几何模型；

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

3.根据权利要求书1所述的一种基于组织分化的长骨骨折愈合仿真系统，其特征在于：所述的骨折区域生物力学有限元分析模块(2)实现其功能的具体过程为：

1)将骨折区域三维表面几何模型进行网格化分，使连续的几何模型离散化，得到骨折区域有限元模型；

所述的网格划分包括面网格划分和体网格划分两个步骤；面网格划分过程用于将三维表面模型进行优化，包括：表面模型优化，平滑处理，修补漏洞；表面模型的优化通过减小表面模型的三角面片来实现，该过程只需将相邻的两个顶点合并到一个新的顶点上，并延续原有的拓扑关系；平滑处理的过程中，对三维的面网格模型进行去噪；修补漏洞的过程中，通过将模型当中的空洞提取成空间多边形，然后对空洞多边形进行三角化的方法实现；体网格划分的过程是将面网格模型进行拉伸、旋转步骤来实现的；

通过网格划分得到的骨折区域有限元模型包括单元编号和节点坐标两部分；

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

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

2)在有限元模型上施加外加载荷，并设置边界条件。载荷的大小由骨所承受的力的大小决定，实验对象不同，所受的力也不同；

3)对骨折区域有限元模型进行骨折区域初始环境设置。骨折区域由皮质骨和骨痂区域两部分组成。初始骨折区域环境设置包括皮质骨材料属性赋值，皮质骨血供赋值，初始骨痂材料属性赋值，初始骨痂区域血供赋值；

距离骨折断端5mm内的皮质骨血供设置为0％，其余部分皮质骨血供设置为100％；

骨痂外周3mm内血供设置为30％，骨痂内部血供设置为0％；

4)将骨折区域看作双相多孔弹性模型，由多孔弹性理论得到骨痂单元的本构方程，平衡方程和几何方程，并通过有限单元法计算骨痂单元应力刺激S，具体过程为：

a.本构方程

$\left\{\begin{array}{c}{\mathrm{\σ}}_{rr}\\ {\mathrm{\σ}}_{\mathrm{\θ}\mathrm{\θ}}\\ {\mathrm{\σ}}_{zz}\\ {\mathrm{\τ}}_{r\mathrm{\θ}}\\ {\mathrm{\τ}}_{\mathrm{\θ}z}\\ {\mathrm{\τ}}_{rz}\end{array}\right\}=\left\{\begin{array}{cccccc}{M}_{11}& M_{12}& M_{13}& 0& 0& 0\\ M_{12}& M_{11}& M_{13}& 0& 0& 0\\ M_{13}& M_{13}& M_{33}& 0& 0& 0\\ 0& 0& 0& M_{44}& 0& 0\\ 0& 0& 0& 0& M_{55}& 0\\ 0& 0& 0& 0& 0& M_{55}\end{array}\right\}\left\{\begin{array}{c}{\mathrm{\ϵ}}_{rr}\\ {\mathrm{\ϵ}}_{\mathrm{\θ}\mathrm{\θ}}\\ {\mathrm{\ϵ}}_{zz}\\ {\mathrm{\γ}}_{r\mathrm{\θ}}\\ {\mathrm{\γ}}_{\mathrm{\θ}z}\\ {\mathrm{\γ}}_{rz}\end{array}\right\}-\left\{\begin{array}{c}\mathrm{\α}\\ \mathrm{\α}\\ {\mathrm{\α}}^{\′}\\ 0\\ 0\\ 0\end{array}\right\}p---\left(1\right)$

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；α，α'分别为各向同性弹性面的Biot系数和轴向Biot系数；p为骨痂单元中的流体压力；M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅分别为脱水的弹性模量矩阵分量；

其中，M₁₁，M₁₂，M₁₃，M₃₃，M₄₄，M₅₅脱水的弹性模量矩阵分量表达式如下所示：

$M_{11}=E_r(E_r-E_r{v}_z^2){(1+v_r)}^{-1}{(E_z-E_z{v}_r-2E_r{v}_z^2)}^{-1}---\left(2\right)$

$M_{12}=E_r(E_z{v}_r-E_r{v}_z^2){(1+v_r)}^{-1}{(E_z-E_z{v}_r-2E_r{v}_z^2)}^{-1}---\left(3\right)$

$M_{13}=E_r{E}_z{v}_z{(E_z-E_z{v}_r-2E_r{v}_z^2)}^{-1}---\left(4\right)$

$M_{33}=E_z^2(1-v_r){(E_z-E_z{v}_r-2E_r{v}_z^2)}^{-1}---\left(5\right)$

M₄₄＝E_r/2(1+ν_r) (6)

M₅₅＝G' (7)

式中，E_r，ν_r分别是各同性弹性层的弹性模量和泊松比；E_z，ν_z分别是轴向弹性模量和泊松比；G'为剪切模量；

b.平衡方程

$\frac{\∂{\mathrm{\σ}}_{rr}}{\∂r}+\frac{1}{r}\frac{\∂{\mathrm{\τ}}_{r\mathrm{\θ}}}{\∂\mathrm{\θ}}+\frac{1}{r}\frac{\∂{\mathrm{\τ}}_{rz}}{\∂z}+\frac{{\mathrm{\σ}}_{rr}-{\mathrm{\σ}}_{\mathrm{\θ}\mathrm{\θ}}}{r}=0---\left(8\right)$

$\frac{\∂{\mathrm{\τ}}_{r\mathrm{\θ}}}{\∂r}+\frac{1}{r}\frac{\∂{\mathrm{\σ}}_{\mathrm{\θ}\mathrm{\θ}}}{\∂\mathrm{\θ}}+\frac{1}{r}\frac{\∂{\mathrm{\τ}}_{\mathrm{\θ}z}}{\∂z}+2\frac{{\mathrm{\τ}}_{\mathrm{\θ}z}}{r}=0---\left(9\right)$

$\frac{\∂{\mathrm{\τ}}_{rz}}{\∂r}+\frac{1}{r}\frac{\∂{\mathrm{\τ}}_{\mathrm{\θ}z}}{\∂\mathrm{\θ}}+\frac{1}{r}\frac{\∂{\mathrm{\σ}}_{zz}}{\∂z}+\frac{{\mathrm{\τ}}_{rz}}{r}=0---\left(10\right)$

式中，σ_rr，σ_θθ，σ_zz为正应力，τ_rθ，τ_θz，τ_rz为剪应力；r为径向半径；

c.几何方程

${\mathrm{\ϵ}}_{rr}=\frac{\∂u_r}{\∂r}---\left(11\right)$

${\mathrm{\ϵ}}_{\mathrm{\θ}\mathrm{\θ}}=\frac{1}{r}\frac{\∂u_{\mathrm{\θ}}}{\∂\mathrm{\θ}}+\frac{u_r}{r}---\left(12\right)$

${\mathrm{\ϵ}}_{zz}=\frac{\∂u_z}{\∂z}---\left(13\right)$

${\mathrm{\γ}}_{r\mathrm{\θ}}=\frac{1}{r}\frac{\∂u_r}{\∂\mathrm{\θ}}+\frac{\∂u_{\mathrm{\θ}}}{\∂r}-\frac{u_{\mathrm{\θ}}}{r}---\left(14\right)$

${\mathrm{\γ}}_{\mathrm{\θ}z}=\frac{1}{r}\frac{\∂u_z}{\∂\mathrm{\θ}}---\left(15\right)$

${\mathrm{\γ}}_{rz}=\frac{\∂u_z}{\∂r}---\left(16\right)$

式中，ε_rr，ε_θθ，ε_zz为正应变，γ_rθ，γ_θz，γ_rz为主应变；u_r，u_θ，u_z分别为三个方向上的位移；r为径向半径；

通过上述方程的求解得到骨痂单元的正应变σ_rr，σ_θθ，σ_zz，由正应变可得到骨痂单元受到的畸变应变：

$D=\frac{1}{\sqrt{2}}\sqrt{{({\mathrm{\σ}}_{rr}-{\mathrm{\σ}}_{\mathrm{\θ}\mathrm{\θ}})}^2+{({\mathrm{\σ}}_{rr}-{\mathrm{\σ}}_{zz})}^2+{({\mathrm{\σ}}_{\mathrm{\θ}\mathrm{\θ}}-{\mathrm{\σ}}_{zz})}^2}---\left(17\right)$

式中，D为骨痂单元受到的畸变应变；σ_rr，σ_θθ，σ_zz分别为各个方向上的正应变；

骨痂单元中液体的流速V为：

$V=-\frac{k}{u}\frac{\∂p}{\∂r}---\left(18\right)$

其中，k为骨痂中液体的达西渗透系数；u液体粘度；p为液体压力；

由此可得到骨痂单元所受的力学刺激S为：

$S=\frac{D}{a}+\frac{V}{b}---\left(19\right)$

其中，D为骨痂单元受到的畸变应变；V为骨痂单元中液体流速；a，b分别为经验常数。

4.根据权利要求书1所述的一种基于组织分化的长骨骨折愈合仿真系统，其特征在于：所述的骨痂单元组织分化模块(3)实现其功能的具体过程为：

1)设置输入变量隶属度函数

骨痂单元组织分化过程由六个输入变量协同决定，分别为单元力学刺激，骨痂单元骨含量，骨痂单元软骨含量，骨痂单元血供，周围骨痂单元骨含量，周围骨痂单元血供。输入变量的精确值通过隶属度函数转变为相应输入变量的隶属度。输入变量隶属度函数由梯形函数表示；

2)输入隶属度函数模糊化

将输入变量的隶属度函数用模糊逻辑语言值来表示：

单元力学刺激可分为五个等级，分别为：低，中，高；

单元血供可分为三个等级，分别为：低，中，高；

单元骨含量可分为三个等级，分别为：低，中，高；

单元软骨含量可分为三个等级，分别为：低，中，高；

3)设置输出变量隶属度函数

输出变量共有三个，分别为骨痂单元血供改变量，骨痂单元骨含量改变量和骨痂单元软骨含量改变量。设置输出变量隶属度函数，输出隶属度函数由梯形函数表示，且将输出隶属度函数用如下的模糊逻辑语言值表示：

骨痂单元血供改变量分为三个等级，分别为：低，中，高；

骨痂单元骨含量改变量分为三个等级，分别为：低，中，高；

骨痂单元软骨含量改变量分为三个等级，分别为：低，中，高；

4)模糊规则的编写

骨痂内组织分化主要包括：血管再生，膜内骨化，软骨生成，软骨钙化，软骨骨化五个过程。通过模糊规则分别表示这五个组织分化过程，由此可得到组织分化后，骨痂单元中骨的模糊逻辑语言值，软骨的模糊逻辑语言值以及血供的模糊逻辑语言值。具体过程如下：

过程一、血管再生

如果单元应力刺激为非高且骨痂单元血供为低且周围骨痂单元血供为非低，那么骨痂单元血供增加；

如果单元应力刺激为非高且骨痂单元血供为中且周围骨痂单元血供为高，那么骨痂单元血供增加；

如果单元应力刺激为非高且骨痂单元血供为高，那么骨痂单元血供增加；

过程二、膜内骨化

如果骨痂单元软骨含量为低且单元应力刺激为中且骨痂单元血供为高且周围骨痂单元血供为高，那么骨痂单元骨含量增加；

过程三、软骨生成

如果骨痂单元骨含量为非高且骨痂单元软骨含量为低且单元应力刺激为中，那么骨痂单元软骨含量增加；

如果骨痂单元软骨含量为非低且单元应力刺激为高，那么骨痂单元软骨含量增加；

如果骨痂单元软骨含量为高且单元应力刺激为中，那么骨痂单元软骨含量增加；

过程四、软骨钙化

如果骨痂单元软骨含量为非低且单元应力刺激为中且骨痂单元血供为非低且周围骨痂单元骨含量为非低，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

如果骨痂单元软骨含量为非低且单元应力刺激为高且骨痂单元血供为非低且周围骨痂单元骨含量为非低，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

过程五、软骨骨化

如果骨痂单元骨含量为高且骨痂单元软骨含量为低且单元应力刺激为中且骨痂单元血供为非低且周围骨痂单元骨含量为高，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

如果骨痂单元骨含量为高且骨痂单元软骨含量为低且单元应力刺激为低且骨痂单元血供为非低且周围骨痂单元骨含量为高，那么骨痂单元骨含量增加，骨痂单元软骨含量减少；

5)输出变量反模糊化

通过模糊规则的作用，得到输出变量的模糊逻辑语言值，利用面积中心法得到输出变量的隶属度；

6)更新输出变量含量

将得到的输出变量隶属度经过转化得到输出变量的改变量，从而得到经过组织分化后，骨痂单元骨组织，软骨组织，结缔组织含量以及骨痂单元血供。具体过程如下：

a.计算经组织分化后骨痂单元骨组织含量

${\mathrm{\μ}}_{Bone,n}^{(k+1)}={\mathrm{\μ}}_{Bone,n}^{\left(k\right)}+\mathrm{\Δ}t\frac{d}{dt}{\mathrm{\μ}}_{Bone,n}^{\left(k\right)}---\left(20\right)$

式中，为第n个骨痂单元在第k+1迭代步中骨组织的含量；为第n个骨痂单元在第k迭代步中骨的含量；为单位时间第n个骨痂单元在第k迭代步中骨的生成量；Δt为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中骨的改变量由下式得到：

$\frac{d}{dt}{\mathrm{\μ}}_{Bone,n}^{\left(k\right)}=w_{Bone}^{\left(k\right)}{T}_{Bone}---\left(21\right)$

式中，为单位时间第n个骨痂单元在第k迭代步中骨的生成量；为骨痂单元中骨含量的隶属度；T_Bone在骨痂单元骨的转化率；

b.计算经组织分化后骨痂单元软骨骨组织含量

${\mathrm{\μ}}_{Cartilage,n}^{(k+1)}={\mathrm{\μ}}_{Cartilage,n}^{\left(k\right)}+\mathrm{\Δ}t\frac{d}{dt}{\mathrm{\μ}}_{Cartilage,n}^{\left(k\right)}---\left(22\right)$

式中，为第n个骨痂单元在第k+1迭代步中软骨的含量；为第n个骨痂单元在第k迭代步中软骨的含量；为单位时间第n个骨痂单元在第k迭代步中软骨的生成量；Δt为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中软骨的改变量由下式得到：

$\frac{d}{dt}{\mathrm{\μ}}_{Cartilage,n}^{\left(k\right)}=w_{Cartilage}^{\left(k\right)}{T}_{Cartilage}---\left(23\right)$

式中，为单位时间第n个骨痂单元在第k迭代步中软骨的生成量；为骨痂单元中软骨含量的隶属度；T_Cartilage在骨痂单元软骨的转化率；

c.计算经组织分化后骨痂单元血供

${\mathrm{\μ}}_{Perfusion,n}^{(k+1)}={\mathrm{\μ}}_{Perfusion,n}^{\left(k\right)}+\mathrm{\Δ}t\frac{d}{dt}{\mathrm{\μ}}_{Perfusion,n}^{\left(k\right)}---\left(24\right)$

式中，为第n个骨痂单元在第k+1迭代步中结缔组织的含量；为第n个骨痂单元在第k迭代步中结缔组织的含量；为单位时间第n个骨痂单元在第k迭代步中结缔组织的生成量；Δt为时间步长；

其中单位时间第n个骨痂单元在第k迭代步中软骨的改变量由下式得到：

$\frac{d}{dt}{\mathrm{\μ}}_{Perfusion,n}^{\left(k\right)}=w_{Perfusion}^{\left(k\right)}{T}_{Perfusion}---\left(25\right)$

式中，为单位时间第n个骨痂单元在第k迭代步中结缔组织的生成量；为骨痂单元中结缔组织含量的隶属度；T_Perfusion在骨痂单元结缔组织的转化率；

d.计算经组织分化后骨痂单元结缔组织含量

骨痂单元由骨组织，软骨组织，结缔组织三部分组成，三者之间的关系如下所示：

μ_Bone+μ_Cartilage+μ_ConnTissue＝1 (26)

式中，μ_Bone为骨痂单元骨组织含量；μ_Cartilage为骨痂单元软骨组织含量；μ_ConnTissue为骨痂单元结缔组织含量；

由此可得组织分化后骨痂单元结缔组织含量。

5.根据权利要求书1所述的一种基于组织分化的长骨骨折愈合仿真系统，其特征在于：所述的程序终止判断模块(4)实现其功能的具体过程为：

1)判断骨痂单元材料属性

判断当前骨痂单元材料属性是否与骨的材料属性相同，若不同，则程序执行下面的2)步骤；若不等，程序执行3)步骤；

2)更新骨痂单元材料属性

若骨痂单元材料属性不等于骨的材料属性，则骨痂单元材料属性进行更新，进入下一个迭代步，骨痂材料属性更新公式如下：

$E_{Callus,n}^{(k+1)}=E_{Bone}{{\mathrm{\μ}}_{Bone,n}^{(k+1)}}^3+E_{Canilage}{{\mathrm{\μ}}_{Canilage,n}^{(k+1)}}^3+E_{ConnTissue}{{\mathrm{\μ}}_{ConnTissue,n}^{(k+1)}}^3---\left(27\right)$

式中，为第n个骨痂单元在第k+1迭代步中的弹性模量；E_Bone为骨的弹性模量；为第n个骨痂单元在第k+1迭代步中骨的含量；E_Cartilage为软骨的弹性模量；为第n个骨痂单元在k+1迭代步中的软骨的含量；E_ConnTissue为结缔组织的弹性模量；为第n个骨痂单元在k+1迭代步中结缔组织的含量；

$v_{Callus,n}^{(k+1)}=v_{Bone}{\mathrm{\μ}}_{Bone,n}^{(k+1)}+v_{Cartilage}{\mathrm{\μ}}_{Cartilage,n}^{(k+1)}+v_{ConnTissue}{\mathrm{\μ}}_{ConnTissue,n}^{(k+1)}---\left(28\right)$

式中，为第n个骨痂单元在第k+1迭代步中的泊松比；ν_Bone为骨的泊松比；为第n个骨痂单元在第k+1迭代步中骨的含量；ν_Cartilage为软骨的泊松比；为第n个骨痂单元在k+1迭代步中的软骨的含量；ν_ConnTissue为结缔组织的泊松比；为第n个骨痂单元在k+1迭代步中结缔组织的含量；

3)程序结束

若所有骨痂单元材料属性等于骨的材料属性，程序终止并输出愈合时间。