CN107391784B - 一种基于拓扑优化技术的松质骨多孔结构建模方法 - Google Patents

Abstract

本发明涉及生物医疗技术领域，一种基于拓扑优化技术的松质骨多孔结构建模方法，包括以下步骤：(1)建立松质骨模型，(2)求解松质骨模型的体积分数，(3)求解松质骨模型的力学性能，(4)提出本方法的拓扑优化模型，(5)求解多孔结构，(6)还原CAD模型。本发明考虑到真实松质骨的性能，设计得到的模型更加符合真实松质骨结构的要求，并且大大减少了模型的数据量，满足加工和后处理的要求，为进一步的骨科临床手术应用奠定了基础。

Description

一种基于拓扑优化技术的松质骨多孔结构建模方法

技术领域

本发明涉及一种基于拓扑优化技术的松质骨多孔结构建模方法，属于生物医疗技术领域。

背景技术

股骨是人体最大的长管状骨，它起着支撑人体的重要作用。股骨头位于股骨的近端，它是股骨中一个十分重要的部位。近年来，股骨头坏死越来越多的出现在病人的生活中，通过股骨头自身的恢复机制很难使其复原，因此需要通过外部植入体替换的方法使其恢复正常。传统的植入体多是实心植入体，这种植入体虽然具有良好的刚度和强度，但是它的生物相容性差，主要原因是股骨头内部完全是松质骨结构，松质骨是由网状结构的骨板排列而成，具有多孔的特性，除了起到支撑的作用之外，同时也为营养物质的传送以及骨髓的流动提供一个固定的通道。

为了解决实心植入体存在的问题，研究人员提出了采用多孔结构植入体代替实心植入体。多孔结构植入体是由一系列周期性排列的多孔结构组成，它具有多孔的特性，有利于骨组织的内生长，且骨组织与多孔植入体的结合界面会更加牢固。与传统实心植入体相比，它不仅生物相容性更好，而且可以通过改变体积分数以及结构形式的方式对植入体的整体的性能进行调控。因此，多孔结构植入体在很长一段时间内受到了人们的广泛关注。

传统的多孔结构设计方法通常可以分为以下两种，分别是构造实体几何法和基于图像的建模方法。构造实体几何法是一种以CAD技术为基础的程序化建模方法。它是以简单体元为基础，利用现有的工程软件UG、CATIA、PRO-E、SolidWorks等，通过布尔运算得到复杂的曲面以及空间结构的方法。通常用到的体元有圆柱、圆锥、长方体、球体等，将其进行求交集或者并集等操作，可以得到复杂的空间结构。构造实体几何法进行多孔结构重建所得到的模型具有结构简单、便于加工以及后处理的特点，但是无法对结构的整体性能加以控制。基于图像的建模方法是以现有的医学图像CT以及MRI等为基础，直接利用图像进行模型建立的方法，这种方法最大的优势就是可以保证模型的真实性。图像建模方法得到的多孔结构模型具有真实结构的性能，应用于数值分析尚可，但是给加工以及后期处理带来困难。

发明内容

为了克服现有技术中存在的不足，本发明目的是提供一种基于拓扑优化技术的松质骨多孔结构建模方法。该方法将松质骨的力学性能和体积分数应用于多孔结构的设计中，使得到的结果既能满足真实松质骨的力学性能要求，又能满足工业加工的要求。拓扑优化技术可以根据给定的目标函数和约束条件对特定性能的结构进行设计，结果具有模型数据量小，便于加工等特点，为进一步的骨科临床手术应用奠定了基础。

为了实现上述发明目的，解决己有技术中存在的问题，本发明采取的技术方案是：一种基于拓扑优化技术的松质骨多孔结构建模方法，包括以下步骤：

步骤1、建立松质骨模型，为了保证松质骨模型的真实性，本发明基于医学图像建立松质骨模型，具体包括以下子步骤：

(a)根据设计要求选取股骨头并对其进行物理处理，利用Micro-CT仪器扫描得到松质骨结构的医学图像；

(b)利用Mimics软件对Micro-CT仪器扫描得到的医学图像进行处理，包括窗宽和窗位调整、阈值分割、区域增长，得到具有清晰骨小梁结构的松质骨图像，并建立松质骨结构模型；

步骤2、求解松质骨模型的体积分数，将在Mimics中建立得到的松质骨模型导出为STL格式，并将其导入到Magics软件中得到松质骨模型的体积信息，并基于此体积信息对松质骨模型的体积分数V_m/V进行求解，其中V_m是松质骨实心部分的体积，V是松质骨实心部分和空心部分的总体积；

步骤3、求解松质骨模型的力学性能，选用代表体元法对模型施加约束条件和外载荷，并基于ANSYS软件对松质骨模型的9个力学性能参数进行求解，具体包括以下子步骤：

(a)对松质骨模型施加约束，可通过公式(1)进行描述，

式中，1,3,5行分别代表x＝0，y＝0和z＝0的面分别在x方向，y方向和z方向的位移为0；u，v，w分别是x，y，z三个方向的位移；a，b，c分别是代表体元三个边的长度；2,4,6行分别代表x＝a，y＝b，z＝c的面的位移分别为δ_x，δ_y，δ_z，它们的值由ANSYS求解得到；

(b)求解松质骨模型的等效弹性模量，包括三个等效弹性模量，分别是x方向的弹性模量E_x，y方向的弹性模量E_y，z方向的弹性模量E_z，可通过公式(2)、(3)、(4)进行描述，

式中，p为施加的外载荷，它的大小不影响结果；a为模型x方向的边长，b为模型y方向的边长，c为模型z方向的边长；

(c)求解松质骨模型的等效泊松比，包括三个等效泊松比，分别对应于xy面的泊松比v_xy，yz面的泊松比v_yz和zx面的泊松比v_zx，可通过公式(5)、(6)、(7)进行描述，

式中，δ_xy为在x＝a的面施加载荷p后y＝b的面的位移，δ_yz为在y＝b的面施加载荷p后z＝c的面的位移，δ_zx为在z＝c的面施加载荷p后x＝a的面的位移；

(d)求解松质骨模型的等效剪切模量，包括三个等效剪切模量，分别是与xy面平行的剪切模量G_xy，与yz面平行的剪切模量G_yz，与zx面平行的剪切模量G_zx，可通过公式(8)、(9)、(10)进行描述，

步骤4、提出本方法的拓扑优化模型，其中包括目标函数、约束条件以及控制方程，具体包括以下子步骤：

(a)本发明引入拓扑优化中的材料模型，选用Solid Isotropic Material WithPenalization(SIMP)方法，可通过公式(11)进行描述，

E(x)＝ρ(x)^pE⁰ (11)

式中，E(x)为材料点x处的单元材料的弹性模量，E⁰为材料本身的弹性模量，ρ(x)为单元材料的密度，p是惩罚因子；

(b)为了能够实时监测迭代过程中的材料模型的等效性能，采用均匀化方法通过公式(12)、(13)对迭代过程中的材料模型性能进行求解，

式中，Y为单胞域，

是符合上式的广义位移函数；v_i表示单胞中定义的连续函数；E_ijkl表示单胞中材料的弹性常数；i，j,k，l表示坐标编号，二维取1，2，三维取1，2，3，E_ijmn表示y_n处的四阶弹性张量，m，n同样是坐标编号，二维取1，2，三维取1，2，3，y_n，y_j表示不同维度的坐标，V_Ω表示设计域，通过公式(12)求解得到广义位移函数

然后通过公式(13)求解得到等效的弹性模量；

(c)基于SIMP材料插值模型和均匀化方法，本发明以松质骨结构的九个性能参数和体积分数为设计目标提出了结构建模优化列式，可通过公式(14)进行描述，

式中，ρ_e表示单元密度，具体可以参考SIMP方法；e表示单元编号，N表示单元的总数目；f(ρ_e)表示目标函数，E_x、E_y、E_z、G_xy、G_yz、G_zx、v_xy、v_yz、v_zx为多孔结构的九个设计参数，同时也是要设计要达到的性能，

表示设计迭代过程中九个等效性能实时的值，V^upp表示设计的体积上限，V表示单元材料分布的总体积，E(x)＝ρ(x)^pE⁰表示引入的密度插值函数，Homogenization表示均匀化方法的控制方程；

步骤5、求解多孔结构，首先确定设计域，给定一个r×r×r的正方体网格设计域，建议给定敏度过滤半径r_min为1.5，建议最大迭代步数为100步，体积约束，给定实心材料无孔隙条件下的弹性模量和泊松比，控制某一特定阈值以上的单元在ANSYS中显示，所述阈值特定值由具体情况进行选取，范围在0-1之间；

步骤6、还原CAD模型，将ANSYS中显示的多孔结构利用CAD软件建立得到符合实际的CAD模型。

本发明有益效果是：一种基于拓扑优化技术的松质骨多孔结构建模方法，包括以下步骤：(1)建立松质骨模型，(2)求解松质骨模型的体积分数，(3)求解松质骨模型的力学性能，(4)提出本方法的拓扑优化模型，(5)求解多孔结构，(6)还原CAD模型。与已有技术相比，本发明考虑到真实松质骨的性能，设计得到的模型更加符合真实松质骨结构的要求，并且大大减少了模型的数据量，满足加工和后处理的要求，为进一步的骨科临床手术应用奠定了基础。

附图说明

图1是本发明方法步骤流程图。

图2是本发明基于代表体元法求解复合材料等效性能约束施加图。

其中：(a)是主视图，(b)是左视图。

图3是本发明得到的多孔结构模型图。

具体实施方式

下面结合附图对本发明作进一步说明。

如图1所示，一种基于拓扑优化技术的松质骨多孔结构建模方法，包括以下步骤：

步骤1、建立松质骨模型，为了保证松质骨模型的真实性，本发明基于医学图像建立松质骨模型，具体包括以下子步骤：

(a)根据设计要求选取股骨头并对其进行物理处理，利用Micro-CT仪器扫描得到松质骨结构的医学图像；

(b)利用Mimics软件对Micro-CT仪器扫描得到的医学图像进行处理，包括窗宽和窗位调整、阈值分割、区域增长，得到具有清晰骨小梁结构的松质骨图像，并建立松质骨结构模型；

步骤2、求解松质骨模型的体积分数，将在Mimics中建立得到的松质骨模型导出为STL格式，并将其导入到Magics软件中得到松质骨模型的体积信息，并基于此体积信息对松质骨模型的体积分数V_m/V进行求解，其中V_m是松质骨实心部分的体积，V是松质骨实心部分和空心部分的总体积；

步骤3、求解松质骨模型的力学性能，选用代表体元法对模型施加约束条件和外载荷，并基于ANSYS软件对松质骨模型的9个力学性能参数进行求解，具体包括以下子步骤：

(a)以图2为例，对复合材料模型施加约束，可通过公式(1)进行描述，

式中，1,3,5行分别代表x＝0，y＝0和z＝0的面分别在x方向，y方向和z方向的位移为0；u，v，w分别是x，y，z三个方向的位移；a，b，c分别是代表体元三个边的长度；2,4,6行分别代表x＝a，y＝b，z＝c的面的位移分别为δ_x，δ_y，δ_z，它们的值由ANSYS求解得到；

(b)求解松质骨模型的等效弹性模量，包括三个等效弹性模量，分别是x方向的弹性模量E_x，y方向的弹性模量E_y，z方向的弹性模量E_z，可通过公式(2)、(3)、(4)进行描述，

式中，p为施加的外载荷，它的大小不影响结果；a为模型x方向的边长，b为模型y方向的边长，c为模型z方向的边长；

(c)求解松质骨模型的等效泊松比，包括三个等效泊松比，分别对应于xy面的泊松比v_xy，yz面的泊松比v_yz和zx面的泊松比v_zx，可通过公式(5)、(6)、(7)进行描述，

式中，δ_xy为在x＝a的面施加载荷p后y＝b的面的位移，δ_yz为在y＝b的面施加载荷p后z＝c的面的位移，δ_zx为在z＝c的面施加载荷p后x＝a的面的位移；

(d)求解松质骨模型的等效剪切模量，包括三个等效剪切模量，分别是与xy面平行的剪切模量G_xy，与yz面平行的剪切模量G_yz，与zx面平行的剪切模量G_zx，可通过公式(8)、(9)、(10)进行描述，

步骤4、提出本方法的拓扑优化模型，其中包括目标函数、约束条件以及控制方程，具体包括以下子步骤：

(a)本发明引入拓扑优化中的材料模型，选用Solid Isotropic Material WithPenalization(SIMP)方法，可通过公式(11)进行描述，

E(x)＝ρ(x)^pE⁰ (11)

式中，E(x)为材料点x处的单元材料的弹性模量，E⁰为材料本身的弹性模量，ρ(x)为单元材料的密度，p是惩罚因子；

(b)为了能够实时监测迭代过程中的材料模型的等效性能，采用均匀化方法通过公式(12)、(13)对迭代过程中的材料模型性能进行求解，

式中，Y为单胞域，

是符合上式的广义位移函数；v_i表示单胞中定义的连续函数；E_ijkl表示单胞中材料的弹性常数；i，j,k，l表示坐标编号，二维取1，2，三维取1，2，3，E_ijmn表示y_n处的四阶弹性张量，m，n同样是坐标编号，二维取1，2，三维取1，2，3，y_n，y_j表示不同维度的坐标，V_Ω表示设计域，通过公式(12)求解得到广义位移函数

然后通过公式(13)求解得到等效的弹性模量；

(c)基于SIMP材料插值模型和均匀化方法，本发明以松质骨结构的九个性能参数和体积分数为设计目标提出了结构建模优化列式，可通过公式(14)进行描述，

式中，ρ_e表示单元密度，具体可以参考SIMP方法；e表示单元编号，N表示单元的总数目；f(ρ_e)表示目标函数，E_x、E_y、E_z、G_xy、G_yz、G_zx、v_xy、v_yz、v_zx为多孔结构的九个设计参数，同时也是要设计要达到的性能，

表示设计迭代过程中九个等效性能实时的值，V^upp表示设计的体积上限，V表示单元材料分布的总体积，E(x)＝ρ(x)^pE⁰表示引入的密度插值函数，Homogenization表示均匀化方法的控制方程；

步骤5、求解多孔结构，首先确定设计域，给定一个r×r×r的正方体网格设计域，建议给定敏度过滤半径r_min为1.5，建议最大迭代步数为100步，九个目标性能参数，体积约束，给定实心材料无孔隙条件下的弹性模量和泊松比，控制某一特定阈值以上的单元在ANSYS中显示，所述阈值特定值由具体情况进行选取，范围在0-1之间；综合考虑，本发明给定一个30×30×30的正方体网格设计域，共含有27000个单元，引入SIMP材料插值模型，选用惩罚因子为3.，给定敏度过滤半径r_min＝1.5，最大迭代步数为100步，目标性能参数为E_x＝2.041，E_y＝1.609，E_z＝2.519，v_xy＝0.28，v_yz＝0.24，v_zx＝0.20，G_xy＝0.797，G_yz＝0.649，G_zx＝1.049，体积约束为松质骨的体积分数V_m/V＝0.3274，控制阈值为0.45；

步骤6、还原CAD模型，将ANSYS中显示的多孔结构利用CAD软件建立得到符合实际的CAD模型，最终呈现一个完整的多孔结构模型如图3所示。

本发明优点在于：一种基于拓扑优化技术的松质骨多孔结构建模方法，考虑到真实松质骨的性能，设计得到的模型更加符合真实松质骨结构的要求，并且大大减少了模型的数据量，满足加工和后处理的要求，为进一步的骨科临床手术应用奠定了基础。

Claims (1)

1.一种基于拓扑优化技术的松质骨多孔结构建模方法，其特征在于，包括以下步骤：

步骤1、建立松质骨模型，为了保证松质骨模型的真实性，基于医学图像建立松质骨模型，具体包括以下子步骤：

(a)根据设计要求选取股骨头并对其进行物理处理，利用Micro-CT仪器扫描得到松质骨结构的医学图像；

(b)利用Mimics软件对Micro-CT仪器扫描得到的医学图像进行处理，包括窗宽和窗位调整、阈值分割、区域增长，得到具有清晰骨小梁结构的松质骨图像，并建立松质骨结构模型；

步骤2、求解松质骨模型的体积分数，将在Mimics中建立得到的松质骨模型导出为STL格式，并将其导入到Magics软件中得到松质骨模型的体积信息，并基于此体积信息对松质骨模型的体积分数BV/TV进行求解，其中BV是松质骨实心部分的体积，TV是松质骨实心部分和空心部分的总体积；

步骤3、求解松质骨模型的力学性能，选用代表体元法对模型施加约束条件和外载荷，并基于ANSYS软件对松质骨模型的9个力学性能参数进行求解，具体包括以下子步骤：

(a)对松质骨模型施加约束，可通过公式(1)进行描述，

式中，1,3,5行分别代表x＝0，y＝0和z＝0的面分别在x方向，y方向和z方向的位移为0；u，v，w分别是x，y，z三个方向的位移；松质骨模型代表体元在x轴方向,y轴方向和z轴方向的长度分别为a,b,c；2,4,6行分别代表x＝a，y＝b，z＝c的面的位移分别为δ_x，δ_y，δ_z，它们的值由ANSYS求解得到；

(b)求解松质骨模型的等效弹性模量，包括三个等效弹性模量，分别是x方向的弹性模量E_x，y方向的弹性模量E_y，z方向的弹性模量E_z，可通过公式(2)、(3)、(4)进行描述，

式中，p为施加的外载荷，它的大小不影响结果；a,b,c分别为松质骨模型代表体元在x轴方向,y轴方向和z轴方向的边长；

(c)求解松质骨模型的等效泊松比，包括三个等效泊松比，分别对应于xy面的泊松比v_xy，yz面的泊松比v_yz和zx面的泊松比v_zx，可通过公式(5)、(6)、(7)进行描述，

式中，δ_xy为在x＝a的面施加外载荷p后y＝b的面的位移，δ_yz为在y＝b的面施加外载荷p后z＝c的面的位移，δ_zx为在z＝c的面施加外载荷p后x＝a的面的位移；

(d)求解松质骨模型的等效剪切模量，包括三个等效剪切模量，分别是与xy面平行的剪切模量G_xy，与yz面平行的剪切模量G_yz，与zx面平行的剪切模量G_zx，可通过公式(8)、(9)、(10)进行描述，

步骤4、提出本方法的拓扑优化模型，其中包括目标函数、约束条件以及控制方程，具体包括以下子步骤：

(a)引入拓扑优化中的材料模型，选用Solid Isotropic Material WithPenalization(SIMP)方法，可通过公式(11)进行描述，

F(x)＝ρ(x)^qF⁰ (11)

式中，F(x)为材料点x处的单元材料的弹性模量，F⁰为材料本身的弹性模量，ρ(x)为单元材料的密度，q是惩罚因子；

(b)为了能够实时监测迭代过程中的材料模型的等效性能，采用均匀化方法通过公式(12)、(13)对迭代过程中的材料模型性能进行求解，

式中，Y为单胞域，

是符合上式的广义位移函数；v_i表示单胞中定义的连续函数；E_ijkl表示单胞中材料的弹性常数；i，j,k，l表示坐标编号，二维取1，2，三维取1，2，3，E_ijmn表示y_n处的四阶弹性张量，m，n同样是坐标编号，二维取1，2，三维取1，2，3，y_n，y_j表示不同维度的坐标，V_Ω表示设计域，通过公式(12)求解得到广义位移函数

然后通过公式(13)求解得到宏观等效的弹性模量

(c)基于SIMP材料插值模型和均匀化方法，以松质骨结构的九个性能参数和体积分数为设计目标提出了结构建模优化列式，可通过公式(14)进行描述，

式中，ρ_e表示单元密度，具体可以参考SIMP方法；e表示单元编号，N表示单元的总数目；f(ρ_e)表示目标函数，E_x、E_y、E_z、G_xy、G_yz、G_zx、v_xy、v_yz、v_zx为多孔结构的九个设计参数，同时也是要设计要达到的性能，

表示设计迭代过程中九个等效性能实时的值，V^upp表示设计的体积上限，V表示单元材料分布的总体积，F(x)＝ρ(x)^qF⁰，F(x)为材料点x处的单元材料的弹性模量，表示引入的密度插值函数，Homogenization表示均匀化方法的控制方程；

步骤5、求解多孔结构，首先确定设计域，给定一个r×r×r的正方体网格设计域，建议给定敏度过滤半径r_min为1.5，建议最大迭代步数为100步，九个目标性能参数，体积约束，给定实心材料无孔隙条件下的弹性模量和泊松比，控制某一特定阈值以上的单元在ANSYS中显示，所述特定阈值由具体情况进行选取，范围在0-1之间；

步骤6、还原CAD模型，将ANSYS中显示的多孔结构利用CAD软件建立得到符合实际的CAD模型。

Images