[发明专利]基于超弹性固体相特性关节软骨两相模型的建立方法

摘要

基于超弹性固体相特性关节软骨两相模型的建立方法，本发明涉及关节软骨两相模型的建立方法。为了解决将软骨的固相视作线弹性材料的模型针对关节软骨产生较大变形的时候与实际情况有很大的差异、甚至失效的问题。本发明首先根据关节的CT数据重建得到需要研究的人体组织的几何模型，然后获取模型的单元编号和节点坐标作为形变计算输入的原始数据，建立基于混合物理论的关节软骨两相模型及基于v‑p变量的控制方程，然后建立关节软骨力学平衡方程并进行有限元计算。本发明适用于关节软骨两相模型的建立和仿真。

主权项

基于超弹性固体相特性关节软骨两相模型的建立方法，其特征在于包括以下步骤：步骤1、几何模型和网格划分：以关节为研究对象，通过获得关节的CT数据，将得到的数据以DICOM格式输出并储存在计算机中；通过MIMICS中图像分割功能将CT数据中关节软骨的区域从其他组织中分离出来，通过三维重建得到需要研究的人体组织的几何模型；从DICOM数据到生成网格模型并获取模型的单元编号和节点坐标作为形变计算输入的原始数据；步骤2、建立基于混合物理论的关节软骨两相模型及基于v‑p变量的控制方程，包括以下步骤：步骤2.1、建立关节软骨的两相模型：将关节软骨视为由超弹性固体和理想流体组成的两相介质混合物，且两相具有独立的运动规律，用s表示固体相，f表示液体相，那么φs表示固体相体积分数，φf表示液体相体积分数；在初始构型中，两相的初始体积分数为φs0、φf0，两相的初始密度为且均匀分布；由饱和下的体积分数式得到：φs+φf＝1     (1)质量平衡方程为：其中，表示矢量微分算子符号；vs是固体相速度，vf是液体相速度；在拟静态问题中，动量平衡方程转化为静力平衡方程，加速度a等于零；忽略外部体积力的作用，在两相混合物中，由间隙液体流过多孔固体所产生的摩擦阻力形成的动量交换为Ps、Pf，且有Ps＝‑Pf＝K′(vf‑vs)   (3)式中，K′为扩散阻力系数；式中，κ为软骨渗透率；固体相动量为Ps，液体相动量Pf；令σs和σf分别为固相和液相的柯西应力张量，关节软骨在当前构形的动量平衡方程为：由软骨的不可压缩性和力学特性得到关节软骨弹性固体相和非粘性流体相的力学本构方程为：σs＝‑φspI+σe    (7)σf＝‑φfpI     (8)其中，I是单位矩阵，p为压力；σe为固体相的有效应力；经过以上推导得出由质量平衡方程(2)、动量平衡方程(5)(6)和两相混合物的本构方程(7)(8)组成了关节软骨的两相模型；根据第二Piola‑Kirchhoff应力张量的定义可知第二Piola‑Kirchhoff应力张量S与弹性应力张量σe的关系为：σe＝J‑1FSFT＝FSFT      (9)其中，F为变形梯度矩阵；J为弹性体积比，是变形梯度矩阵F的雅可比行列式；S为第二Piola‑Kirchhoff应力张量；在有限变形中，考虑到物体受力前后空间位置的改变，令t＝0时刻物体在空间所占据的区域V0为初始构形；在当前时刻t，物体所占据的空间区域V表示为当前构形；以初始构形为参考构型，利用物质描述法可得出体元、面元在有限变形中的变化，即当前构形中的单元物质的体积与参考构型中的体积满足以下关系：dV＝JdV0      (10)步骤2.2、建立超弹性模型：关节软骨的超弹性模型由Helmholtz应变能函数来表示，即：α0、α1、α2、β是材料的力学参数；I1、I2、I3分别是C的第一、第二、第三主不变量；对于不可压缩材料而言，I3＝det(C)＝1，简化后的超弹性模型为：步骤2.3、求解第二Piola‑Kirchhoff应力张量S：关节软骨的超弹性固体相，应力应变不再满足线性关系，弹性矩阵与应变张量有关，第二Piola‑Kirchhoff应力张量S和Green应变张量E都是以初始构形为参考构型的应力和应变张量，在一个瞬时变形过程中，本构关系为：dS＝DdE；Green应变张量E为：弹性张量D表示为：式中，D中的元素表示S中的第i″行第j″列的元素，Ek″l″表示E中的第k″行第l″列的元素；令：2α1+4α2(1+E22+E33)＝G2α1+4α2(1+E11+E22)＝H2α1+4α2(1+E11+E33)＝B得到对称矩阵D，该矩阵中的对称部分用sym表示：步骤2.4、建立基于v‑p变量的控制方程选取速度v、压力p为未知量建立基于v‑p变量的有限元方程，进而表达关节软骨基于v‑p变量的控制方程；先将式(3)和式(8)带到式(6)中，通过式(4)得到：在饱和条件式(1)下，将式(20)带入式(2)中，得到只包含固体相速度的关系式：接下来由动量守恒方程(5)(6)和本构方程(7)(8)联立可得：由式(21)和式(22)就构成了基于v‑p变量的控制方程；步骤3、建立关节软骨力学平衡方程并进行有限元计算：步骤3.1、建立边界条件：流体通量定义为则基于v‑p变量的控制方程的边界条件为：式中，x指在某一时刻t时当前构型中的空间点；为已知的固相位移，为总的载荷力，为已知的压力，n*是边界Γ的外法线向量；Γp和ΓQ分别代表液体相体积域Ωf的边界，和ΓT分别代表固体相体积域Ωs的边界，且边界满足条件：步骤3.2、构造形函数：定义有限元四面体单元的形函数为：N＝[Ni′·I3×3 Nj′·I3×3 Nm′·I3×3 Np′·I3×3]   (44)其中，式中，Ve是四面体单元的体积，ai′、bi′、ci′、di′为与Ni′对应的系数；Nj′、Nm′、Np′与Ni′类似，有步骤3.3、建立关节软骨力学平衡方程：基于有限元求解的单元类型为四面体单元，构造插值函数；通过基于v‑p变量的控制方程(21)、(22)以及边界条件(23)‑(26)，采用伽辽金加权余量法建立关节软骨力学平衡方程；步骤3.4、求解关节软骨力学平衡方程，实现关节软骨两相模型的仿真。