CN113435082B - 一种用于骨折手术机器人的力学调控仿真系统 - Google Patents

Abstract

本发明公开了一种用于骨折手术机器人的力学调控仿真系统，包括以下几个模块：(1)模型处理模块：确定患骨的具体尺寸并形成二维模型；(2)数值分析模块：对二维模型进一步处理并进行有限元分析；(3)力学刺激效果模块：对模型施加力学刺激，仿真模拟骨折康复机器人所施加的机械刺激对患骨愈合过程的影响；(4)联合仿真处理模块：修正了二维模型上所施加的力学刺激，并使模型处理模块、数值分析模块、力学刺激效果模拟模块以及联合仿真模块的功能连接起来形成一个迭代周期。仿真模拟骨折手术机器人对患骨进行一个周期的康复治疗；本发明系统通过在患骨有限元模型上进行康复刺激的仿真得到了最佳的康复策略。

Description

一种用于骨折手术机器人的力学调控仿真系统

技术领域

本发明涉及力学调控仿真系统，尤其涉及一种用于骨折手术机器人的力学调控仿真系统。

背景技术

如何实现骨折自然快速愈合是研究骨折术后康复的重要和核心内容，但骨折愈合是一个极其复杂的生物学再生过程，受到机体内环境中内分泌、局部血液循环、生长因子等众多生化信号的影响，内环境的影响因素与骨折愈合的关系只停留在定性实验研究层面难以提出定量的康复方法。沃尔夫定律在十九世纪被提出，它表示人体骨骼的愈合生长与力学环境有关：高压力促进骨骼生长，低压力促使骨骼吸收，自此力学调控模型开始被广泛研究。到二十世纪五十年代Ilizarov通过大量动物实验提出可控张力会刺激骨外膜和骨内膜的骨痂生长。力学调控模型的研究被用来应用到二期骨折愈合的分析和模拟上来，也基于此类模型提出了各种各样的促进二期骨折术后康复的被动力学刺激手段。然而大多数存在骨组织分化对应的力学刺激范围不完整和力与组织分化的时间特性相对模糊的问题。

尽管理论的发展时间已久，骨折愈合力学刺激的研究由于其问题的独特性。难以实现精准的康复模拟和定量的衡量和评价。现有研究只能为经验丰富的临床医师提供一定的指导性建议，不能满足现代临床医用的实际要求。

发明内容

本发明的目的在于克服已有技术的缺点，提供一种高精度、低误差、效果好的用于骨折手术机器人的力学调控仿真系统。

本发明解决上述技术问题采用的技术方案为：

本发明的一种用于骨折手术机器人的力学调控仿真系统，包括：

模型处理模块，所述的模型处理模块用于读取骨折区域的术后CT影像，将CT影像导入mimics软件进行阈值分割，把感兴趣的模型部分分割成一个蒙版，然后将该蒙版重建成蒙版三维模型，在蒙版三维模型上采用mimic软件内部的测量功能进行包括骨折间隙、血肿区域的高度、胫骨的长度、髓内管的半径、皮质骨的半径、血肿区域的半径在内的变量值测量，以测量得到的变量值作为模型关键字参数建立二维模型，并在二维模型中皮质骨的顶端施加轴向力载荷，再利用Python语言编写Abaqus的二维模型cae分析文件；

数值分析模块，数值分析模块首先利用python语言为模型处理模块输出的二维模型cae文件中的模型进行网格划分，将划分好的网格区域通过人工识别划分为三个区域：中央骨痂模型区域、骨髓模型区域、皮质骨模型区域，然后为中央骨痂模型区域的网格单元设定纤维材料属性，骨髓模型区域的网格单元设定骨髓材料属性，皮质骨模型区域的网格单元设定皮质骨材料属性，得到一个新的二维模型cae文件，再采用Abaqus中土力学渗流耦合方法对该新的二维模型cae分析文件进行处理得到新的二维模型的odb有限元分析结果文件，该odb有限元分析结果文件包含二维模型所有网格单元的主应变、流体流速和所有节点的位移和支反力；所述的材料属性包括：杨氏模量、渗透率、泊松比、孔隙度、固体和流体压缩模量；

力学刺激效果模拟模块，所述的力学刺激效果模拟模块首先利用Python语言编写的脚本提取数值分析模块得到的odb有限元分析结果文件中中央骨痂模型区域所有网格单元的主应变和流体流速，并计算得到中央骨痂模型区域每个网格单元的成骨指数s；再利用Python语言编写的脚本打开数值分析模块中得到的新的二维模型cae文件，根据每个网格单元的成骨指数s修改该网格单元的材料属性。若0<s<0.0106，则将该网格单元的材料属性改为骨吸收材料属性，若0.0106≤s<0.02667则将该网格单元的材料属性改为成熟骨材料属性，若0.02667≤s<1，则将该网格单元的材料属性改为未成熟骨材料属性，若1≤s<3，则将该网格单元的材料属性改为软骨材料属性，若3≤s，则将该网格单元的材料属性改为纤维材料属性，得到一个修改后的cae文件；

联合仿真处理模块，所述的联合仿真处理模块包括力学刺激矫正子模块和自动仿真流程控制子模块；

所述的力学刺激矫正子模块首先利用Python语言编写的脚本提取数值分析模块得到的odb文件中皮质骨底端所有节点的支反力和位移，并取所有节点中的支反力和位移的最大值，保存成一个excel文件；然后利用位移的最大值计算中央骨痂模型区域的微应变，最后利用Python语言编写的脚本打开力学刺激效果模拟模块得到的新的二维模型cae文件，修改其中施加在皮质骨顶端的载荷值为新载荷值；

所述的自动仿真流程控制子模块利用matlab编写程序设定循环次数n，然后依次连接模型处理模块、数值分析模块、力学刺激效果模拟模块以及联合仿真模块形成一个迭代周期并且进行n次循环，在每一次仿真循环过程中联合仿真模块中的力学刺激矫正子模块都会生成一个保存本次迭代皮质骨底端所有节点中最大支反力和最大位移的excel文件；

当循环结束时，利用Python语言编写的脚本打开n次循环中生成的所有excel文件，利用所有excel文件的保存结果绘图，第一幅图为皮质骨最大支反力随迭代次数的变化图，其横轴为迭代的次数，纵轴为本次迭代对应的支反力最大值；第二幅图为皮质骨最大位移值随迭代次数的变化图，横轴为迭代的次数，纵轴为本次迭代对应的位移最大值。

本发明的有益效果是：

1.通过双相多孔弹性材料的本构模拟了生长组织材料，比使用其他例如线弹性材料进行骨愈合仿真更具有准确性；

2.基于力学调控模型，骨愈合进程和效果受力学环境的影响是公认的普遍规律，利用力学调控模型，排除了在体外和仿真过程中无法实现的内环境生物因素，提供更加清晰的刺激-响应-愈合的逻辑关系。

3.利用有限元仿真技术精准量化了每个阶段的所要施加的力学刺激，解释了力学刺激与模型变化之间的准确关系。解决了目前时间特性与力学特性难以匹配的难题。

附图说明

图1为一种用于骨折手术机器人的力学调控仿真系统流程图；

图2为骨折有限元模型及各项参数值；

图3为组织分化的机械调节算法图；

图4为力学调控仿真系统软件工作流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如附图1所示本发明的一种用于骨折手术机器人的力学调控仿真系统，包括：

模型处理模块，所述的模型处理模块用于读取骨折区域的术后CT影像，将CT影像导入mimics软件进行阈值分割，把感兴趣的模型部分分割成一个蒙版，然后将该蒙版重建成蒙版三维模型。在蒙版三维模型上采用mimic软件内部的测量功能进行包括骨折间隙、血肿区域的高度、胫骨的长度、髓内管的半径、皮质骨的半径、血肿区域的半径在内的变量值测量，以测量得到的变量值作为模型关键字参数建立二维模型，并在二维模型中皮质骨的顶端施加轴向力载荷，载荷大小如可以采用300N，再利用Python语言编写Abaqus的二维模型cae分析文件。二维模型及参数标注见图2，其中骨折间隙为2a、血肿区域高度为b、胫骨长骨为2c、髓内管半径为d、皮质骨半径为e、血肿区域半径为f。

二维模型是在三维模型的基础上基于二维轴对称思想将骨愈合的三维问题简化为二维平面问题进行研究，使用二维模型的优点包括：缩减分析速度、提高工作效率、降低工作成本。二维模型中骨折间隙为骨折术后骨折断端的间隙不是术前的间隙。血肿区域的高度取两侧血肿区域依附在皮质骨上高度的平均值。血肿区域的半径取自中央骨痂血肿区的最大直径。

数值分析模块，数值分析模块首先利用python语言为模型处理模块输出的二维模型cae文件中的模型进行网格划分，将划分好的网格区域通过人工识别划分为三个区域：中央骨痂模型区域、骨髓模型区域、皮质骨模型区域。然后为中央骨痂模型区域的网格单元设定纤维材料属性，骨髓模型区域的网格单元设定骨髓材料属性，皮质骨模型区域的网格单元设定皮质骨材料属性。得到一个新的二维模型cae文件，再采用Abaqus中土力学渗流耦合方法对该新的二维模型cae分析文件进行处理得到新的二维模型的odb有限元分析结果文件，该odb有限元分析结果文件包含二维模型所有网格单元的主应变、流体流速和所有节点的位移和支反力；所述的材料属性包括：杨氏模量、渗透率、泊松比、孔隙度、固体和流体压缩模量。

所述的每个网格单元的杨氏模量、渗透率、泊松比、孔隙度、固体和流体压缩模量的值可以参照下表选取：

表1 材料属性表

所述的土力学渗流耦合方法原理如下：

骨痂断端区域在生长中为软组织的生长，软组织中含有大量的流动流体，相对于可变形的固体材料组件移动。这种游离组织液可用于多种目的，包括用于运输、加热或冷却的对流路径，以及有助于结构的机械功能。因此，了解和量化软组织的机械反应，包括游离组织液的作用是很重要的，所以软组织必须考虑为双相多孔弹性材料不能忽略其流体性质。骨痂断端区域的软组织在生长中随着愈合时间的延长会分化为纤维组织、软骨组织、未成熟的骨组织、成熟骨组织等4种软组织。这4种软组织组织的材料属性均要视为包含饱和流体且可变形的双相多孔弹性材料。但是双相多孔弹性材料的混合物理论模拟极其复杂且并未有现有的商用软件能够解决此类问题。而土力学的研究与生物力学有关，其中孔隙弹性理论和混合体理论在该领域的理论环节频繁出现。在生物力学方向已经有学者使用混合体理论进行讨论和分析，但是后期研究又基于孔隙弹性理论。多孔弹性理论以连续体的位移和平均相对流体位移为场变量，已经证明Abaqus的土力学渗流耦合仿真结果可以实现多孔弹性理论的本构模型仿真，具体可以参见：参考文献：Wu J Z,Herzog W,Epstein M.有限元软件ABAQUS在双相组织生物力学建模中的应用[J].Journal of Biomechanics,1998,31(2):165-169.)

以下是针对利用土力学渗流耦合本构模型分析多孔弹性软组织材料需要的输入量(每个网格单元的杨氏模量、渗透率、泊松比、孔隙度、固体和流体压缩模量)的推导过程。根据双相多孔弹性理论，各种软组织可以模拟为存在于每个网格单元中(有限元分析中划分的网格单元)的固体和流体成分的混合物。由于这种连续性，体积分数之和必须等于1，并且混合物的表观总密度ρ可表示如下：

ρ＝(1-n)ρ^s+nρ^f (1)

其中n是每个网格单元的孔隙率(等于流体体积分数)，和ρ^s、ρ^f分别是每个网格单元的固相的和液相的表观密度。n_R表示孔隙度，参考构型中的孔隙度如下：

其中J是当前每个网格单元体积dV与每个网格单元参考体积dV_R的比率，或者也是相对于参考构型的变形梯度的决定因素：

其中

是哈密顿算子，标量场通过哈密顿算子运算形成一个矢量场。x是每个网格单元的变形向量。

作用在固体和孔隙总面积上的总应力σ(拉伸正值))被分离为晶间应力(固相应力)σ′和对应于流体无粘性时的流体应力的孔隙压p(压缩正值)、I是单位矩阵。如公式(4)

σ＝σ′-pI (4)

假设忽略重力场、对流体相和固液两相之间的相对加速度，动量守恒得出

div(σ)＝ρü_s (5)

其中是ü_s是每个网格单元中固体的加速度，使用公式(4)，此公式为

如果假设每个网格单元材料属性是各向同性的，如果假设无限小的应变位移关系，可以写出如下本构定律

σ′＝λeI+2με (7)

其中ε表示每个网格单元的应变张量，e表示每个网格单元的扩张应变，λ和μ是与每种材料杨氏模量和泊松比有关的拉梅常数。有限应变可以通过使用总拉格朗日公式或更新拉格朗日公式来计算。有两种可能的方式可以改变组织的体积。第一个是从组织中排出液体。这就产生了根据Darcy′s的相对流体速度

其中，

是每个网格单元中相对于固体速度的相对流体速度矢量，k是每个网格单元的渗透率矢量，v是每个网格单元的运动粘度。组织体积变化的第二种方式是通过压缩固体和液体相本身。如果假设多孔固体的压缩模量远小于非多孔固体的固有压缩模量，则一般质量守恒定律得出：

其中

是每个网格单元中固体的速度矢量，

是每个网格单元中相对于固体速度的相对流体速度矢量，Q是与流体和固体的可压缩性有关的Biot材料参数，该数值的计算方式如下

其中K_f和K_s分别是每个网格单元的流相和固相的固有压缩模量。

基本的有限元假设是根据节点的值对每个单元的坐标、位移和孔压势进行插值。在上述这些方程的基础上，采用标准Galerkin方法进行有限元空间离散。使用欧拉反向时间积分在几个时间步长处评估空间离散化。因此，需要知道六个参数(杨氏模量、渗透率、泊松比、孔隙度、固体和流体压缩模量)才能完整地定义多孔弹性可压缩各向同性模型即完整分析该骨折的有限元模型。

力学刺激效果模拟模块，所述的力学刺激效果模拟模块用于分析在模型处理模块中建立的二维模型模型上施加力学刺激后对骨折愈合的影响。这个过程可以仿真模拟骨折手术机器人所施加的机械刺激对患骨愈合过程的影响。

该模块首先利用Python语言编写的脚本提取数值分析模块得到的odb有限元分析结果文件中中央骨痂模型区域所有网格单元的主应变和流体流速，并按照公式(11)、(12)计算中央骨痂模型区域每个网格单元的成骨指数s。

其中γ为每个网格单元的偏应变，计算公式如下

其中ε₁是中央骨痂模型区域该网格单元的第一主应变、ε₂为中央骨痂模型区域该网格单元的第二主应变、ε₃是中央骨痂模型区域该网格单元的第三主应变、v是每个网格单元的流体流速。这四个量均由脚本从odb结果文件中提取得到，a的值为3.75％，b的值为3μms^-1。

再利用Python语言编写的脚本打开数值分析模块中得到的新的二维模型cae文件，根据每个网格单元的成骨指数s修改该网格单元的材料属性。若0<s<0.0106，则将该网格单元的材料属性改为骨吸收材料属性，若0.0106≤s<0.02667则将该网格单元的材料属性改为成熟骨材料属性，若0.02667≤s<1，则将该网格单元的材料属性改为未成熟骨材料属性，若1≤s<3，则将该网格单元的材料属性改为软骨材料属性，若3≤s，则将该网格单元的材料属性改为纤维材料属性，得到一个修改后的cae文件。其中各种材料属性的具体值参见上述表1。

联合仿真处理模块，所述的联合仿真处理模块包括力学刺激矫正子模块和自动仿真流程控制子模块；

所述的力学刺激矫正子模块首先利用Python语言编写的脚本提取数值分析模块得到的odb文件中皮质骨底端所有节点的支反力和位移，并取所有节点中的支反力和位移的最大值，保存成一个excel文件；然后利用位移的最大值计算中央骨痂模型区域的微应变，最后利用Python语言编写的脚本打开力学刺激效果模拟模块得到的新的二维模型cae文件，修改其中施加在皮质骨顶端的载荷值为新载荷值。

若其2％≤微应变＜10％则保持施加在二维模型中皮质骨顶端的当前轴向力载荷不变；

若微应变≥10％，则

若微应变≤2％，则

若计算的新载荷大于500N则取新载荷为500N。

所述的自动仿真流程控制子模块利用matlab编写程序设定循环次数n，然后依次连接模型处理模块、数值分析模块、力学刺激效果模拟模块以及联合仿真模块形成一个迭代周期并且进行n次循环，在每一次仿真循环过程中联合仿真模块中的力学刺激矫正子模块都会生成一个保存本次迭代皮质骨底端所有节点中最大支反力和最大位移的excel文件；

当循环结束时，利用Python语言编写的脚本打开n次循环中生成的所有excel文件，利用所有excel文件的保存结果绘图，第一幅图为皮质骨最大支反力随迭代次数的变化图，其横轴为迭代的次数，纵轴为本次迭代对应的支反力最大值；第二幅图为皮质骨最大位移值随迭代次数的变化图，横轴为迭代的次数，纵轴为本次迭代对应的位移最大值。这样通过皮质骨最大支反力随迭代次数的变化图可以直观看出模型在周期内是否存在超载情况，通过皮质骨最大位移值随迭代次数的变化图可以得到骨折手术机器人在一个康复治疗周期内的路径。

本联合仿真处理模块根据有限元分析结果修正在二维模型上所施加的力学刺激并使模型处理模块、数值分析模块、力学刺激效果模拟模块以及联合仿真模块的功能连接起来形成一个迭代周期。最终分析整个周期内的力学刺激效果。这个过程可以仿真模拟骨折手术机器人对患骨进行一个周期的康复治疗。

以上对本发明的描述仅仅是示意性的，而不是限制性的，所以，本发明的实施方式并不局限于上述的具体实施方式。如果本领域的普通技术人员受其启示，在不脱离本发明宗旨和权利要求所保护范围的情况下，做出其他变化或变型，均属于本发明的保护范围。

Claims (1)

1.一种用于骨折手术机器人的力学调控仿真系统，其特征在于包括：

模型处理模块，所述的模型处理模块用于读取骨折区域的术后CT影像，将CT影像导入mimics软件进行阈值分割，把感兴趣的模型部分分割成一个蒙版，然后将该蒙版重建成蒙版三维模型，在蒙版三维模型上采用mimic软件内部的测量功能进行包括骨折间隙、血肿区域的高度、胫骨的长度、髓内管的半径、皮质骨的半径、血肿区域的半径在内的变量值测量，以测量得到的变量值作为模型关键字参数建立二维模型，并在二维模型中皮质骨的顶端施加轴向力载荷，再利用Python语言编写Abaqus的二维模型cae分析文件；

数值分析模块，数值分析模块首先利用python语言为模型处理模块输出的二维模型cae文件中的模型进行网格划分，将划分好的网格区域通过人工识别划分为三个区域：中央骨痂模型区域、骨髓模型区域、皮质骨模型区域，然后为中央骨痂模型区域的网格单元设定纤维材料属性，骨髓模型区域的网格单元设定骨髓材料属性，皮质骨模型区域的网格单元设定皮质骨材料属性，得到一个新的二维模型cae文件，再采用Abaqus中土力学渗流耦合方法对该新的二维模型cae分析文件进行处理得到新的二维模型的odb有限元分析结果文件，该odb有限元分析结果文件包含二维模型所有网格单元的主应变、流体流速和所有节点的位移和支反力；所述的材料属性包括：杨氏模量、渗透率、泊松比、孔隙度、固体和流体压缩模量；

力学刺激效果模拟模块，所述的力学刺激效果模拟模块首先利用Python语言编写的脚本提取数值分析模块得到的odb有限元分析结果文件中中央骨痂模型区域所有网格单元的主应变和流体流速，并计算得到中央骨痂模型区域每个网格单元的成骨指数s；再利用Python语言编写的脚本打开数值分析模块中得到的新的二维模型cae文件，根据每个网格单元的成骨指数s修改该网格单元的材料属性， 若0<s<0.0106，则将该网格单元的材料属性改为骨吸收材料属性，若0.0106≤s<0.02667则将该网格单元的材料属性改为成熟骨材料属性，若0.02667≤s<1，则将该网格单元的材料属性改为未成熟骨材料属性，若1≤s<3，则将该网格单元的材料属性改为软骨材料属性，若3≤s，则将该网格单元的材料属性改为纤维材料属性，得到一个修改后的cae文件；

联合仿真处理模块，所述的联合仿真处理模块包括力学刺激矫正子模块和自动仿真流程控制子模块；

所述的力学刺激矫正子模块首先利用Python语言编写的脚本提取数值分析模块得到的odb文件中皮质骨底端所有节点的支反力和位移，并取所有节点中的支反力和位移的最大值，保存成一个excel文件；然后利用位移的最大值计算中央骨痂模型区域的微应变，最后利用Python语言编写的脚本打开力学刺激效果模拟模块得到的新的二维模型cae文件，修改其中施加在皮质骨顶端的载荷值为新载荷值；

所述的自动仿真流程控制子模块利用matlab编写程序设定循环次数n，然后依次连接模型处理模块、数值分析模块、力学刺激效果模拟模块以及联合仿真模块形成一个迭代周期并且进行n次循环，在每一次仿真循环过程中联合仿真模块中的力学刺激矫正子模块都会生成一个保存本次迭代皮质骨底端所有节点中最大支反力和最大位移的excel文件；

当循环结束时，利用Python语言编写的脚本打开n次循环中生成的所有excel文件，利用所有excel文件的保存结果绘图，第一幅图为皮质骨最大支反力随迭代次数的变化图，其横轴为迭代的次数，纵轴为本次迭代对应的支反力最大值；第二幅图为皮质骨最大位移值随迭代次数的变化图，横轴为迭代的次数，纵轴为本次迭代对应的位移最大值。

Images