CN109101686A - 一种基于骨小梁结构形态和力学性能的多孔结构设计方法 - Google Patents

Abstract

本发明涉及生物医疗工程技术领域，一种基于骨小梁结构形态和力学性能的多孔结构设计方法，包括以下步骤：(1)建立骨小梁的三维模型，(2)求解骨小梁的形态学参数，(3)导出骨小梁体网格模型，(4)计算骨小梁的力学性能，(5)优化求解多孔模型，(6)还原实体结构。本发明考虑真实骨小梁的结构形态特点和力学性能，设计的多孔结构实现了结构和功能的统一，极大地减少了设计中的计算量和复杂度，易于后期的加工处理，为骨病的植入物置换治疗提供依据和方法。

Description

一种基于骨小梁结构形态和力学性能的多孔结构设计方法

技术领域

本发明涉及一种基于骨小梁结构形态和力学性能的多孔结构设计方法，属于生物医疗工程技术领域。

背景技术

骨小梁是机体的椎骨及股骨等长骨端内的网状结构，起支撑造血组织和保证血液等物质流通的作用。研究发现，当骨小梁结构产生一定程度的损坏时，整体骨结构的性能明显降低。此外，骨小梁的结构破损不能通过机体自身的恢复机制复原，只能靠药物或者去除破损进行植入物固定等方式治疗。由于破坏方式和受损程度不一致，植入物无法应用于所有损伤治疗中。目前针对骨折等情况的治疗方法，多是通过植入钛合金等生物相容性较好的实心金属棒进行支撑或固定。然而，金属植入物多为实心植入体，密度远大于真实骨组织结构，植入后应力抵抗作用较大，导致骨细胞不易生长并易引起机体的不适，且在骨连接部位会造成一定程度的骨结构损坏使固定失败。

骨小梁结构具有较高的孔隙率，这是骨组织支撑作用的重要影响因素。因此，为了减少实心植入物对机体的损坏，研究人员提出设计仿骨小梁的多孔结构的方法。仿生多孔结构与实心植入物相比，不仅具有良好的生物相容性，而且可以通过对特定性能的控制实现结构的个性化设计。

现有的仿生多孔植入物设计方法主要包含结构和性能设计两大类，结构设计主要是利用数学拟合的方法拟合不同形状的柱状或者板状梁结构，然后通过布尔运算构造不同形态的类骨小梁结构。这种方法可以实现结构的模拟，使之具有一定的弯曲度，但是未考虑力学性能的要求。性能设计主要是以力学性能为目标进行多孔结构设计，然后利用工程软件UG、CAD等重建规则多孔结构。该设计法可以对结构的力学性能加以控制，结构简单，更利于后期的加工实现，但是由于应力分布问题易导致骨细胞生长不连续。

发明内容

为了克服现有技术中存在的不足，本发明目的是提供一种基于骨小梁结构形态和力学性能的多孔结构设计方法。该方法将骨小梁的结构形态和力学性能应用到结构设计中，同时利用拓扑优化方法对结构设计的目标函数和约束条件加以控制，使设计可以降低计算量和复杂度及实现结构和功能的相互统一，更加利于后期的实现。

为了实现上述发明目的，解决已有技术中所存在的问题，本发明采取的技术方案是：一种基于骨小梁结构形态和力学性能的多孔结构设计方法，包括以下步骤：

步骤1、建立骨小梁的三维模型，为了保证正确表征骨小梁的结构形态，重建骨小梁的三维模型作为结构设计的基础，具体包括以下子步骤：

(a)根据实验设计要求截取股骨头内骨小梁结构进行物理处理，利用Micro-CT仪器扫描获得骨小梁的DICOM格式图像；

(b)利用Mimics软件对骨小梁图像进行预处理，包括调整窗宽和窗位，阈值分割和区域增长，去除孤立结构单元，得到边界清晰的骨小梁结构图像，并重建骨小梁三维结构模型；

步骤2、求解骨小梁的形态学参数，为了准确衡量骨小梁的结构形态特点，利用形态学参数来进行定量测量，具体包括以下子步骤：

(a)导出Mimics中建立的三维模型的STL格式文件，并将其导入到Magics软件中得到骨小梁的体积信息V_m，结合包含孔洞和骨小梁的总体积V求解结构的体积分数V_m/V；

(b)利用MATLAB程序对Mimics重建的目标结构的图像进行截取操作，保证每一个目标区域具有相同的图像尺寸和结构体积，同时对结构的整体结构图像不产生影响；

(c)将截取的骨小梁图像导入至图像处理软件ImageJ中，利用BoneJ插件求解骨小梁的形态学参数，包括各向异性即Anisotropy、连通性即Connectivity、结构模型指数即Structure Model Index、厚度即Trabecular Thickness、数量即Trabecular Number和间隙即Trabecular Separation，然后利用测量各向异性中获得的拟合椭球的长轴计算结构方向θ；

(d)将求解的形态学参数结果进行换算，按照结构体素单元数目与结构体积的比例关系，将骨小梁形态学参数计算结果转换为以毫米为单位的参数表示结果；

步骤3、导出骨小梁体网格模型，为了减少模型网格划分时的尺寸差异，在Mimics软件中三维重建后进行体网格划分，选用Hexahedral 8-point单元类型，smoothing模块迭代次数为4，平滑因子为0.5，并选中Close small holes、Filter small parts和improveConnectivity选项，去除骨小梁结构中的孤立结构；

步骤4、计算骨小梁的力学性能，结合代表体元法和有限元法，将力学性能等效为x、y、z方向的弹性模量E，x、y、z方向的泊松比v和x、y、z方向的剪切模量G，基于ANSYS和MATLAB程序利用代表体元的边长和受力后形变长度进行求解，x方向的弹性模量通过公式(1)进行描述，

式中，E_x是x方向的弹性模量，P为代表体元的外载荷，a、b和c分别为x、y和z方向的边长，δ是施加载荷后形变长度，代表体元在x方向施加外载荷后，利用x方向和y方向产生的形变计算泊松比值，通过公式(2)进行描述，

式中，v_xy是x方向施加载荷后的泊松比，δ_xy是x方向载荷y方向形变，δ_xx是x方向载荷x方向的形变，利用求解的弹性模量和泊松比求解剪切模量，通过公式(3)进行描述，

式中，G_xy是x方向施加载荷的剪切模量；

步骤5、优化求解多孔模型，利用变密度法和均匀化法设计目标函数、约束条件及控制方程，进行仿骨小梁结构的模型优化，具体包括以下子步骤：

(a)引入材料模型设计结构，选用Solid Isotropic Material withPenalization即SIMP模型进行材料设计，通过公式(4)进行描述，

E(i)＝ρ(i)^pE⁰ (4)

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

(b)采用均匀化法求解力学性能等效参数，通过式(5)、(6)进行描述，

式中，Y为单胞域，是符合上式的广义位移函数；v_i表示单胞中定义的连续函数；E_ijkl表示单胞中材料的弹性常数；i，j，k，l，m，n表示坐标编号，二维取1，2，三维取1，2，3，E_ijmn表示y_n处的四阶弹性张量，y_n，y_j表示不同维度的坐标，V_Ω为设计域，示表示等效弹性模量；

(c)以骨小梁的结构方向θ，体积分数V_m/V和力学性能等效参数E、v、G为目标函数，以结构体积、单元尺寸、结构内球半径、网格参数、惩罚系数、敏度过滤半径、迭代次数、阈值为约束条件，利用MATLAB程序进行优化计算，得到灰度单元集合模型，随后将模型导入ANSYS中进行结构后处理工作，当单元值大于设定阈值时，标记为实体单元，当单元值小于阈值时，规定为孔，进而得到拓扑优化计算后的清晰连通的实体和孔洞单元模型；

步骤6、还原实体结构，利用UG软件将ANSYS中显示的孔洞模型还原为实体模型。

本发明有益效果是：一种基于骨小梁结构形态和力学性能的多孔结构设计方法，包括以下步骤：(1)建立骨小梁的三维模型，(2)求解骨小梁的形态学参数，(3)导出骨小梁体网格模型，(4)计算骨小梁的力学性能，(5)优化求解多孔模型，(6)还原实体结构。与已有技术相比，本发明考虑真实骨小梁的结构形态特点和力学性能，设计的多孔结构实现了结构和功能的统一，极大地减少了设计中的计算量和复杂度，易于后期的加工处理，为骨病的植入物置换治疗提供依据和方法。

附图说明

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

图2是本发明结合代表体元法求解等效力学性能约束载荷施加图。

图3是本发明优化求解多孔模型实体图。

具体实施方式

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

如图1所示，一种基于骨小梁结构形态和力学性能的多孔结构设计方法，包括以下步骤：

步骤1、建立骨小梁的三维模型，为了保证正确表征骨小梁的结构形态，重建骨小梁的三维模型作为结构设计的基础，具体包括以下子步骤：

(a)根据实验设计要求截取股骨头内骨小梁结构进行磨削、蒸煮等物理处理，利用Micro-CT仪器扫描获得骨小梁的DICOM格式图像；

(b)利用Mimics软件对骨小梁图像进行预处理，包括调整窗宽和窗位，阈值分割和区域增长，去除孤立结构单元，得到边界清晰的骨小梁结构图像，并重建骨小梁三维结构模型；

步骤2、求解骨小梁的形态学参数，为了准确衡量骨小梁的结构形态特点，利用形态学参数来进行定量测量，具体包括以下子步骤：

(a)导出Mimics中建立的三维模型的STL格式文件，并将其导入到Magics软件中得到骨小梁的体积信息V_m，结合包含孔洞和骨小梁的总体积V求解结构的体积分数V_m/V；

(b)利用MATLAB程序对Mimics重建的目标结构的图像进行截取操作，保证每一个目标区域具有相同的图像尺寸和结构体积，同时对结构的整体结构图像不产生影响；

(c)将截取的骨小梁图像导入至图像处理软件ImageJ中，利用BoneJ插件求解骨小梁的形态学参数，包括各向异性即Anisotropy、连通性即Connectivity、结构模型指数即Structure Model Index、厚度即Trabecular Thickness、数量即Trabecular Number和间隙即Trabecular Separation，然后利用测量各向异性中获得的拟合椭球的长轴计算结构方向θ；

(d)将求解的形态学参数结果进行换算，按照结构体素单元数目与结构体积的比例关系即1毫米对应32体素，将骨小梁形态学参数计算结果转换为以毫米为单位的参数表示结果；

步骤3、导出骨小梁体网格模型，为了减少模型网格划分时的尺寸差异，在Mimics软件中三维重建后进行体网格划分，选用Hexahedral 8-point单元类型，smoothing模块迭代次数为4，平滑因子为0.5，并选中Close small holes、Filter small parts和improveConnectivity选项，去除骨小梁结构中的孤立结构；

步骤4、计算骨小梁的力学性能，结合代表体元法和有限元法，将力学性能等效为x、y、z方向的弹性模量E，x、y、z方向的泊松比v和x、y、z方向的剪切模量G，基于ANSYS和MATLAB程序利用代表体元的边长和受力后形变长度进行求解，代表体元示意图如图2所示，代表体元由两种不同材料组成，P为x方向施加外载荷，在图中左侧和底部添加约束，a、b、c分别为代表体元x、y、z方向的边长，x方向的弹性模量通过公式(1)进行描述，

式中，E_x是x方向的弹性模量，P为代表体元的外载荷，a、b和c分别为x、y和z方向的边长，δ是施加载荷后形变长度，代表体元在x方向施加外载荷后，利用x方向和y方向产生的形变计算泊松比值，通过公式(2)进行描述，

式中，v_xy是x方向施加载荷后的泊松比，δ_xy是x方向载荷y方向形变，δ_xx是x方向载荷x方向的形变，利用求解的弹性模量和泊松比求解剪切模量，通过公式(3)进行描述，

式中，G_xy是x方向施加载荷的剪切模量；

步骤5、优化求解多孔模型，利用变密度法和均匀化法设计目标函数、约束条件及控制方程，进行仿骨小梁结构的模型优化，具体包括以下子步骤：

(a)引入材料模型设计结构，选用Solid Isotropic Material withPenalization即SIMP模型进行材料设计，通过公式(4)进行描述，

E(i)＝ρ(i)^pE⁰ (4)

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

(b)采用均匀化法求解力学性能等效参数，通过式(5)、(6)进行描述，

式中，Y为单胞域，是符合上式的广义位移函数；v_i表示单胞中定义的连续函数；E_ijkl表示单胞中材料的弹性常数；i，j，k，l，m，n表示坐标编号，二维取1，2，三维取1，2，3，E_ijmn表示y_n处的四阶弹性张量，y_n，y_j表示不同维度的坐标，V_Ω为设计域，表示等效弹性模量；

(c)以骨小梁的结构方向θ，体积分数V_m/V和力学性能等效参数E、v、G为目标函数，以结构体积(3×3×3)cm、、单元尺寸(2×2×2)cm，结构内球半径(2-2.5mm)、网格参数、惩罚系数(p＝3)、敏度过滤半径(r＝1.5)、迭代次数(建议120次)、阈值(0.5)为约束条件，利用MATLAB程序进行优化计算，得到灰度单元集合模型，随后将模型导入ANSYS中进行结构后处理工作，当单元值大于设定阈值时，标记为实体单元，当单元值小于阈值时，规定为孔，进而得到拓扑优化计算后的清晰连通的实体和孔洞单元模型；

步骤6、还原实体结构，利用UG软件将ANSYS中显示的孔洞模型还原为实体模型，如图3所示。

本发明优点在于：一种基于骨小梁结构形态和力学性能的多孔结构设计方法，考虑真实骨小梁的结构形态和力学性能，设计的仿生多孔结构能实现结构和功能的统一，并且大大减少设计中的计算量和复杂度，利于后期的加工和处理，为临床骨病的植入物置换治疗提供依据和方法。

Claims (1)

1.一种基于骨小梁结构形态和力学性能的多孔结构设计方法，其特征在于包括以下步骤：

步骤1、建立骨小梁的三维模型，为了保证正确表征骨小梁的结构形态，重建骨小梁的三维模型作为结构设计的基础，具体包括以下子步骤：

(a)根据实验设计要求截取股骨头内骨小梁结构进行物理处理，利用Micro-CT仪器扫描获得骨小梁的DICOM格式图像；

(b)利用Mimics软件对骨小梁图像进行预处理，包括调整窗宽和窗位，阈值分割和区域增长，去除孤立结构单元，得到边界清晰的骨小梁结构图像，并重建骨小梁三维结构模型；

步骤2、求解骨小梁的形态学参数，为了准确衡量骨小梁的结构形态特点，利用形态学参数来进行定量测量，具体包括以下子步骤：

(a)导出Mimics中建立的三维模型的STL格式文件，并将其导入到Magics软件中得到骨小梁的体积信息V_m，结合包含孔洞和骨小梁的总体积V求解结构的体积分数V_m/V；

(b)利用MATLAB程序对Mimics重建的目标结构的图像进行截取操作，保证每一个目标区域具有相同的图像尺寸和结构体积，同时对结构的整体结构图像不产生影响；

(c)将截取的骨小梁图像导入至图像处理软件ImageJ中，利用BoneJ插件求解骨小梁的形态学参数，包括各向异性即Anisotropy、连通性即Connectivity、结构模型指数即Structure Model Index、厚度即Trabecular Thickness、数量即Trabecular Number和间隙即Trabecular Separation，然后利用测量各向异性中获得的拟合椭球的长轴计算结构方向θ；

(d)将求解的形态学参数结果进行换算，按照结构体素单元数目与结构体积的比例关系，将骨小梁形态学参数计算结果转换为以毫米为单位的参数表示结果；

步骤3、导出骨小梁体网格模型，为了减少模型网格划分时的尺寸差异，在Mimics软件中三维重建后进行体网格划分，选用Hexahedral 8-point单元类型，smoothing模块迭代次数为4，平滑因子为0.5，并选中Close small holes、Filter small parts和improveConnectivity选项，去除骨小梁结构中的孤立结构；

步骤4、计算骨小梁的力学性能，结合代表体元法和有限元法，将力学性能等效为x、y、z方向的弹性模量E，x、y、z方向的泊松比v和x、y、z方向的剪切模量G，基于ANSYS和MATLAB程序利用代表体元的边长和受力后形变长度进行求解，x方向的弹性模量通过公式(1)进行描述，

式中，E_x是x方向的弹性模量，P为代表体元的外载荷，a、b和c分别为x、y和z方向的边长，δ是施加载荷后形变长度，代表体元在x方向施加外载荷后，利用x方向和y方向产生的形变计算泊松比值，通过公式(2)进行描述，

式中，v_xy是x方向施加载荷后的泊松比，δ_xy是x方向载荷y方向形变，δ_xx是x方向载荷x方向的形变，利用求解的弹性模量和泊松比求解剪切模量，通过公式(3)进行描述，

式中，G_xy是x方向施加载荷的剪切模量；

步骤5、优化求解多孔模型，利用变密度法和均匀化法设计目标函数、约束条件及控制方程，进行仿骨小梁结构的模型优化，具体包括以下子步骤：

(a)引入材料模型设计结构，选用Solid Isotropic Material with Penalization即SIMP模型进行材料设计，通过公式(4)进行描述，

E(i)＝ρ(i)^pE⁰ (4)

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

(b)采用均匀化法求解力学性能等效参数，通过式(5)、(6)进行描述，

式中，Y为单胞域，是符合上式的广义位移函数；v_i表示单胞中定义的连续函数；E_ijkl表示单胞中材料的弹性常数；i，j，k，l，m，n表示坐标编号，二维取1，2，三维取1，2，3，E_ijmn表示y_n处的四阶弹性张量，y_n，y_j表示不同维度的坐标，V_Ω为设计域，表示等效弹性模量；

(c)以骨小梁的结构方向θ，体积分数V_m/V和力学性能等效参数E、v、G为目标函数，以结构体积、单元尺寸、结构内球半径、网格参数、惩罚系数、敏度过滤半径、迭代次数、阈值为约束条件，利用MATLAB程序进行优化计算，得到灰度单元集合模型，随后将模型导入ANSYS中进行结构后处理工作，当单元值大于设定阈值时，标记为实体单元，当单元值小于阈值时，规定为孔，进而得到拓扑优化计算后的清晰连通的实体和孔洞单元模型；

步骤6、还原实体结构，利用UG软件将ANSYS中显示的孔洞模型还原为实体模型。