一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法技术

一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法，涉及数值模拟仿真技术领域。该发明专利技术考虑新骨形成中骨痂的粘弹塑性力学行为，编写UMAT子程序，计算应变历史，利用模糊逻辑实现基于应变调控的组织分化并进行数值分析，再现牵张成骨术中复杂的骨再生动态过程。所述方法包括：(1)对截骨术后截骨区域进行几何和有限元建模；(2)建立截骨区域包含骨痂粘弹塑性特征的生物力学模型，模拟牵张过程并计算骨痂单元应变；(3)利用模糊逻辑实现基于应变调控的组织分化，计算骨痂单元各成分含量并更新其材料属性；(4)对牵张大变形后骨痂区域有限元网格重划分并进行状态数据映射。

A numerical simulation method of distraction osteogenesis based on viscoelastic and plastic mechanical properties of tissue

【技术实现步骤摘要】
一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法
：本专利技术属于数值模拟仿真
，涉及一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法。
技术介绍
：大段骨缺损的修复一直是骨科临床中需要面临的难题，也是近些年的研究焦点之一。其修复方法包括，自体或异体骨移植，基于3D打印组织工程支架，以及牵张成骨术。一直以来牵张成骨术是临床上治疗大段骨缺损的金标准之一。牵张成骨术是通过截骨术将骨骼切开，利用牵张器对骨组织施加稳定和缓慢的牵张力以刺激组织细胞的再生和生长，促进牵张间隙中新骨形成和矿化，从而达到延长骨骼的目的。牵张成骨的骨再生过程受益于骨痂区域生物力学环境，但是目前有关该环境如何调控骨痂区域的组织的形成并不了解。因此，开发一种能够准确再现牵张成骨骨再生过程的数值仿真方法，对于探究牵张成骨中力学刺激和骨再生之间复杂的关系具有重要意义，对临床上术前探寻实施最佳牵张成骨手术方案具有一定的价值。目前，牵张成骨方面的数值模拟研究较少，缺少能够精确再现牵张成骨这一复杂骨再生过程的仿真模型。主要问题在于：(1)用于模拟牵张成骨的数值模型相对简单，通常为二维的简化模型，不能准确呈现截骨区域的几何形状和所处的力学环境；(2)在牵张成骨过程中，骨痂组织被拉伸并长时间保持在某一位置，表现出粘弹塑性力学行为。以往的模型通常将材料视为多孔弹性性质，未能真实描述出在牵张期骨痂区域的粘弹塑性行为；(3)牵张仿真时未考虑骨痂中骨吸收情况，因此无法再现骨痂的真实大小和形状。
技术实现思路
：本专利技术提供了一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法，主要目的是为了解决现有牵张成骨模型仿真过程中，仿真模型过于简化，往往只将模型材料视为多孔弹性性质，不能描述真实的骨痂区域的力学状态。此外，先前的牵张成骨模型无法再现骨痂的真实大小和形状。为此，本专利技术采用对象截骨区域个体化的几何和有限元模型，赋予牵张成骨过程中骨痂区域的粘弹塑性力学行为，再现骨痂区域力学状态，利用模糊逻辑实现基于应变调控的组织分化，开发了一种更准确的牵张成骨复杂骨再生过程的仿真方法。为了达到上述目的，本专利技术通过以下技术方案来实现：一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法，该方法包括以下步骤：步骤A1：对截骨术后截骨区域进行几何和有限元建模：获取对象目标骨骼的CT图像，对CT图像进行自动分割并进行三维重建和网格划分，得到对象截骨区域个体化的几何和有限元单元模型；截骨区域分为两部分：皮质骨和骨痂，则皮质骨和骨痂分别也被分成多个有限元单元模型；步骤A2：建立截骨区域包含骨痂粘弹塑性特征的生物力学模型，设置材料属性并进行边界条件和初始状态设置，模拟牵张并计算骨痂中各单元应变；步骤A3：利用模糊逻辑实现基于应变调控的组织分化，计算骨痂中各单元中各成分含量并更新骨痂中各单元的材料属性；步骤A4：对牵张大变形后骨痂区域的有限元网格进行重划分并进行状态数据映射，并进行循环即从步骤A1开始进行循环。作为本专利技术进一步的技术方案，步骤A2中所述特征，对步骤A1中重建的对象截骨区域的有限元模型上施加外加载荷，并设置材料属性和初始状态值，设置材料属性包括：在ABAQUS的用户材料设置皮质骨杨氏模量、泊松比和骨痂杨氏模量、泊松比、屈服应力、粘度系数(通过包含粘弹塑性生物力学模型的UMAT子程序进行定义材料类型)，骨痂区域初始阶段被假定充满结缔组织。对截骨区域有限元模型进行初始状态设置。初始状态设置包括初始皮质骨骨含量为100％、软骨含量为0％和血供为100％，初始骨痂区域骨含量为0％、软骨含量为0％和血供为0％，并将其保存在Excel文件中。在牵张过程中，骨痂组织被牵拉后会表现出粘弹塑性力学行为：由牵张位移引起的弹性应变随时间发生衰减(粘弹性变形行为)，并部分转化为固体的永久变形(塑性变形行为)；采用与时间有关的粘弹塑性模型来描述骨痂区域单元的应力-应变本构关系，包括以下步骤：步骤B1：骨痂区域粘弹塑性生物力学建模，包括：利用宾汉-麦克斯韦粘弹塑性模型建立骨痂区域的生物力学模型，该粘弹塑性模型由线性弹簧、牛顿粘壶和一个摩擦件组成。通过计算模型本构关系，编写UMAT子程序，计算骨痂区域的粘弹塑性行为，具体为；若摩擦件中应力σf小于屈服应力σyield，此时粘弹塑性模型为线弹性性质，可表示为：σ＝Esεs(1)式中，σ为粘弹塑性模型总应力，εs为线性弹簧应变，Es线性弹簧弹性模量；若摩擦件中应力σf超过临界应力σyield，此时粘弹塑性模型为粘弹塑性质，可表示为：其中，式中，η为粘度系数，t为时间。将上述方程输入到ABAQUS的UMAT子程序中，用于计算骨痂区域各单元的应变。步骤B2：应变结果数据的整理，包括：在ABAQUS中进行粘弹塑性有限元分析得到有限元模型各单元应变随时间变化曲线。对于每个牵张时间周期，使用ndiff等距分割每个牵张时间周期的应变样本，按照每个样本最大峰值刺激进行采样，得到ndiff个用于组织分化算法的畸变应变γ0和膨胀应变ε0；式中，ε1,ε2,ε3分别为骨痂区域中各单元的三个主应变。作为本专利技术进一步的技术方案，步骤A3中所述特征，利用应变调控组织分化理论，计算骨痂区域中各单元各成分含量并更新骨痂区域各单元材料属性的方法包括以下步骤：步骤C1：建立输入变量隶属度函数；模糊控制器将模型中截骨区域有限元模型(即所有皮质骨和骨痂)中每个单元的七个变量作为输入，包括：(1)单元中骨含量的百分比，(2)单元中软骨含量的百分比，(3)单元的膨胀应变，(4)单元的畸变应变，(5)单元的血供，(6)相邻单元中骨含量的百分比，(7)相邻单元血供，其中，有限元模型各单元中初始骨、软骨和血供含量通过步骤A2获得，即若单元为皮质骨中的单元则采用皮质骨中单元对应的参数，若单元为骨痂中的单元则采用骨痂中单元对应的参数；相邻单元的影响通过对每个有限元单元质心处进行采样，利用切比雪夫距离获得每个单元的相邻区域，再通过高斯核函数判断相邻区域每个相邻单元对其影响的权重，最后通过加权平均判定相邻区域内相邻单元的影响。在模糊逻辑控制的模糊输入中，输入变量的精确值通过隶属度函数转变为相应的输入变量的模糊值。输入变量的隶属度函数用模糊逻辑语言值表示；分别对应着：低，中，高等；步骤C2：建立描述在力学刺激(畸变应变和膨胀应变)下血管生成、膜内骨化、软骨生成、软骨内骨化和组织破坏等组织分化的模糊控制规则(如建立由17条if-then语言组成的模糊控制规则)，如果存在或符合规则条件，则发生相对应的组织分化，这些规则描述了在力学刺激(畸变应变和膨胀应变)下的血管生成、膜内骨化、软骨生成、软骨内骨化和组织破坏的过程。除了上述建立的用于描述组织分化的规则之外，在本专利技术方法中，还进一步结合骨痂中骨成熟和骨吸收的生物学过程，从而能够更好地预测骨痂的整体大小和形状；步骤C3:建立输出变量隶属度函数；模糊逻辑本文档来自技高网...

【技术保护点】
1.一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法，其特征在于，该方法包括以下步骤：/n步骤A1：对截骨术后截骨区域进行几何和有限元建模：获取对象目标骨骼的CT图像，对CT图像进行自动分割并进行三维重建和网格划分，得到对象截骨区域个体化的几何和有限元单元模型；截骨区域分为两部分：皮质骨和骨痂，则皮质骨和骨痂分别也被分成多个有限元单元模型；/n步骤A2：建立截骨区域包含骨痂粘弹塑性特征的生物力学模型，设置材料属性并进行边界条件和初始状态设置，模拟牵张并计算骨痂中各单元应变；/n步骤A3：利用模糊逻辑实现基于应变调控的组织分化，计算骨痂中各单元中各成分含量并更新骨痂中各单元的材料属性；/n步骤A4：对牵张大变形后骨痂区域的有限元网格进行重划分并进行状态数据映射，并进行循环即从步骤A1开始进行循环。/n

【技术特征摘要】
1.一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法，其特征在于，该方法包括以下步骤：
步骤A1：对截骨术后截骨区域进行几何和有限元建模：获取对象目标骨骼的CT图像，对CT图像进行自动分割并进行三维重建和网格划分，得到对象截骨区域个体化的几何和有限元单元模型；截骨区域分为两部分：皮质骨和骨痂，则皮质骨和骨痂分别也被分成多个有限元单元模型；
步骤A2：建立截骨区域包含骨痂粘弹塑性特征的生物力学模型，设置材料属性并进行边界条件和初始状态设置，模拟牵张并计算骨痂中各单元应变；
步骤A3：利用模糊逻辑实现基于应变调控的组织分化，计算骨痂中各单元中各成分含量并更新骨痂中各单元的材料属性；
步骤A4：对牵张大变形后骨痂区域的有限元网格进行重划分并进行状态数据映射，并进行循环即从步骤A1开始进行循环。


2.按照权利要求1所述的一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法，其特征在于，
步骤A2中所述特征，对步骤A1中重建的对象截骨区域的有限元模型上施加外加载荷，并设置材料属性和初始状态值，设置材料属性包括：在ABAQUS的用户材料设置皮质骨杨氏模量、泊松比和骨痂杨氏模量、泊松比、屈服应力、粘度系数，骨痂区域初始阶段被假定充满结缔组织。对截骨区域有限元模型进行初始状态设置。初始状态设置包括初始皮质骨骨含量为100％、软骨含量为0％和血供为100％，初始骨痂区域骨含量为0％、软骨含量为0％和血供为0％，并将其保存在Excel文件中；
在牵张过程中，骨痂组织被牵拉后会表现出粘弹塑性力学行为：由牵张位移引起的弹性应变随时间发生衰减，并部分转化为固体的永久变形；采用与时间有关的粘弹塑性模型来描述骨痂区域单元的应力-应变本构关系，包括以下步骤：
步骤B1：骨痂区域粘弹塑性生物力学建模，包括：
利用宾汉-麦克斯韦粘弹塑性模型建立骨痂区域的生物力学模型，该粘弹塑性模型由线性弹簧、牛顿粘壶和一个摩擦件组成。通过计算模型本构关系，编写UMAT子程序，计算骨痂区域的粘弹塑性行为，具体为；
若摩擦件中应力σf小于屈服应力σyield，此时粘弹塑性模型为线弹性性质，可表示为：
σ＝Esεs(1)
式中，σ为粘弹塑性模型总应力，εs为线性弹簧应变，Es线性弹簧弹性模量；
若摩擦件中应力σf超过临界应力σyield，此时粘弹塑性模型为粘弹塑性质，可表示为：



其中，
式中，η为粘度系数，t为时间；
将上述方程用于计算骨痂区域各单元的应变；
步骤B2：应变结果数据的整理，包括：
在ABAQUS中进行粘弹塑性有限元分析得到有限元模型各单元应变随时间变化曲线。对于每个牵张时间周期，使用ndiff等距分割每个牵张时间周期的应变样本，按照每个样本最大峰值刺激进行采样，得到ndiff个用于组织分化算法的畸变应变γ0和膨胀应变ε0；






式中，ε1,ε2,ε3分别为骨痂区域中各单元的三个主应变。


3.按照权利要求1所述的一种基于组织粘弹塑性力学特性的牵张成骨数值仿真方法，其特征在于，步骤A3中所述特征，利用应变调控组织分化理论，计算骨痂区域中各单元各成分含量并更新骨痂区域各单元材料属性的方法包括以下步骤：
步骤C1：建立输入变量隶属度函数；
模糊控制器将截骨区域有限元模型(即所有皮质骨和骨痂)中每个单元的七个变量作为输入，包括：(1)单元中骨含量的百分比，(2)单元中软骨含量的百分比，(3)单元的膨胀应变，(4)单元的畸变应变，(5)单元的血供，(6)相邻单元中骨含量的百分比，(7)相邻单元血供。其中，有限元模型各单元初始骨、软骨和血供含量通过步骤A2获得；即若单元为皮质骨中的单元则采用皮质骨中单元对应的参数，若单元为骨痂中的单元则采用骨痂中单元对应的参数；相邻单元的影响通过对每个有限元单元质心处进行采样，利用切比雪夫距离获得每个单元的相邻区域，再通过高斯核函数判断相邻区域每个相邻单元对其影响的权重，最后通过加权平均判定相邻区域内相邻单元的影响；在模糊逻辑控制的模糊输入中，输入变量的精确值通过隶属度函数转变为相应的输入变量的模糊值。输入变量的隶属度函数用模糊逻辑语言值表示；分别对应着：低，中，高等；
步骤C2：建立描述在力学刺激下血管生成、膜内骨化、软骨生成、软骨内骨化和组织破坏等组织分化的模糊控制规则，如果存在或符合规则条件，则发生相对应的组织分化，这些规则描述了在力学刺激下的血管生成、膜内骨化、软骨生成、软骨内骨化和组织破坏的过程；
除了上述建立的用于描述组织分化的规则之外，在本发明方法中，还进一步结合骨痂中骨成熟和骨吸收的生物学过程，从而能够更好地预测骨痂的整体大小和形状；
步骤C3：建立输出变量隶属度函数；
模糊逻辑控制最终输出骨痂区域中各单元血供改变量、骨痂区域各单元骨含量改变...

【专利技术属性】
技术研发人员：杨海胜，付瑞森，刘有军，冯懿俐，
申请(专利权)人：北京工业大学，
类型：发明
国别省市：北京;11