首页> 中国专利> 基于超弹性固体相特性关节软骨两相模型的建立方法

基于超弹性固体相特性关节软骨两相模型的建立方法

页面导航

摘要
著录项
法律信息
说明书
相似文献

摘要

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

著录项

公开/公告号CN106202738A

专利类型发明专利
公开/公告日2016-12-07

原文格式PDF
申请/专利权人哈尔滨理工大学;
展开▼

申请/专利号CN201610555178.1
发明设计人王沫楠;
展开▼

申请日2016-07-14
分类号G06F17/50;
代理机构哈尔滨市松花江专利商标事务所;
代理人杨立超
地址 150080 黑龙江省哈尔滨市南岗区学府路52号
入库时间 2023-06-19 01:07:21

法律信息

法律状态公告日

法律状态信息

法律状态
2020-07-17

未缴年费专利权终止 IPC(主分类):G06F17/50 授权公告日:20171219 终止日期:20190714 申请日:20160714

专利权的终止
2017-12-19

授权

授权
2017-01-04

实质审查的生效 IPC(主分类):G06F17/50 申请日:20160714

实质审查的生效
2016-12-07

公开

公开

说明书

技术领域

本发明涉及关节软骨两相模型的建立方法。

背景技术

关节软骨主要由水和复合有机质组成。复合有机质主要包括胶原纤维和蛋白多糖，间隙流体中主要含有水，这些流体的流动不仅对软骨的力学性能有着重要影响，还与无血管组织的营养物质传输有着密切关系。同时，关节软骨使关节具有极低的摩擦系数，起到良好的润滑作用。关节软骨的润滑作用对动关节运动和承受载荷的能力具有极其重要的意义。因此，研究关节软骨的力学特性，认识其在运动过程中的应力、应变状况及其在运动中各力学量的变化规律，对医学研究、临床诊断和新型生物材料的开发，具有非常重要的理论意义和实用价值。

在关节软骨力学建模的研究初期，人们将关节软骨看作是各相同性线弹性单相材料。在一段时间里，研究人员接受了粘弹性的单相关节软骨模型，然而，通过这种本构关系得到的仿真结果总是与软骨的实际情况和实验数据有很大的差异。目前，关节软骨被视为两相模型。但是，这些研究中大都将软骨的固相看作是线弹性材料，线弹性模型只对小变形有效，针对关节软骨在力的施加过程中会产生较大变形的问题，将软骨的固相看作是线弹性材料的模型就会与实际情况有很大的差异，甚至失效。

发明内容

本发明为了解决将软骨的固相视作线弹性材料的模型针对关节软骨产生较大变形的时候与实际情况有很大的差异、甚至失效的问题。

基于超弹性固体相特性关节软骨两相模型的建立方法，包括以下步骤：

步骤1、几何模型和网格划分：

以关节为研究对象，通过获得关节的CT数据，将得到的数据以DICOM格式输出并储存在计算机中；通过MIMICS中图像分割功能将CT数据中关节软骨的区域从其他组织中分离出来，通过三维重建得到需要研究的人体组织的几何模型；

从DICOM数据到生成网格模型并获取模型的单元编号和节点坐标作为形变计算输入的原始数据，获取过程如图1所示；

步骤2、建立基于混合物理论的关节软骨两相模型及基于v-p变量的控制方程，包括以下步骤：

步骤2.1、建立关节软骨的两相模型：

将关节软骨视为由超弹性固体和理想流体组成的两相介质混合物，且两相具有独立的运动规律，用s表示固体相，f表示液体相，那么φ_s表示固体相体积分数，φ_f表示液体相体积分数；在初始构型中，两相的初始体积分数为φ_s0、φ_f0，两相的初始密度为且均匀分布；

由饱和下的体积分数式得到：

φ_s+φ_f＝1(1)

质量平衡方程为：

$> ▿ \cdot (φ_{s} v_{s} + φ_{f} v_{f}) = 0 - - - (2)$ >

其中，表示矢量微分算子符号；v_s是固体相速度，v_f是液体相速度；

在拟静态问题中，动量平衡方程转化为静力平衡方程，加速度a等于零；此研究中忽略外部体积力的作用，在两相混合物中，由间隙液体流过多孔固体所产生的摩擦阻力形成的动量交换为P_s、P_f，且有

P_s＝-P_f＝K′(v_f-v_s)(3)

式中，K′为扩散阻力系数；

$> K^{'} = \frac{{(φ_{f})}^{2}}{κ} - - - (4)$ >

式中，κ为软骨渗透率；固体相动量为P_s，液体相动量P_f；

令σ_s和σ_f分别为固相和液相的柯西应力张量，关节软骨在当前构形的动量平衡方程为：

$> ▿ \cdot σ_{s} + P_{s} = 0 - - - (5)$ >

$> ▿ \cdot σ_{f} + P_{f} = 0 - - - (6)$ >

由软骨的不可压缩性和力学特性得到关节软骨弹性固体相和非粘性流体相的力学本构方程为：

σ_s＝-φ_spI+σ_e(7)

σ_f＝-φ_fpI(8)

其中，I是单位矩阵，p为压力；σ_e为固体相的有效应力也叫做与固相的变形量相一致的弹性应力张量；

经过以上推导得出由质量平衡方程(2)、动量平衡方程(5)(6)和两相混合物的本构方程(7)(8)组成了关节软骨的两相模型；

根据第二Piola-Kirchhoff应力张量的定义可知第二Piola-Kirchhoff应力张量S与弹性应力张量σ_e的关系为：

σ_e＝J^-1FSF^T＝FSF^T(9)其中，F为变形梯度矩阵；J为弹性体积比，是变形梯度矩阵F的雅可比行列式；S为第二Piola-Kirchhoff应力张量，该参数在步骤2.3中求解，通过S的联结而得到关节软骨的超弹性两相介质模型；

在有限变形中，考虑到物体受力前后空间位置的改变，令t＝0时刻物体在空间所占据的区域V₀为初始构形；在当前时刻t，物体所占据的空间区域V表示为当前构形；

以初始构形为参考构型，利用物质描述法可得出体元、面元在有限变形中的变化，即当前构形中的单元物质的体积与参考构型中的体积满足以下关系：

dV＝JdV₀(10)

步骤2.2、建立超弹性模型：

关节软骨的超弹性模型由Helmholtz应变能函数来表示，即：

$> Ψ = α_{0} \frac{e^{α_{1} (I_{1} - 3) + α_{2} (I_{2} - 3)}}{{I_{3}}^{β}} - - - (14)$ >

α₀、α₁、α₂、β是材料的力学参数，来源于数学模型与实验数据的拟合结果；I₁、I₂、I₃分别是C的第一、第二、第三主不变量；

对于不可压缩材料而言，I₃＝det(C)＝1，简化后的超弹性模型为：

$> Ψ = α_{0} e^{α_{1} (I_{1} - 3) + α_{2} (I_{2} - 3)} - - - (15)$ >

步骤2.3、求解第二Piola-Kirchhoff应力张量S：

关节软骨的超弹性固体相，应力应变不再满足线性关系，弹性矩阵与应变张量有关，第二Piola-Kirchhoff应力张量S和Green应变张量E都是以初始构形为参考构型的应力和应变张量，在一个瞬时变形过程中，本构关系为：dS＝DdE；

Green应变张量E为：

$> E = (\begin{matrix} E_{11} & E_{12} & E_{13} \\ E_{21} & E_{22} & E_{23} \\ E_{31} & E_{32} & E_{33} \end{matrix}) - - - (16)$ >

$> S = \frac{\partial Ψ}{\partial E} = Ψ (\begin{matrix} 2 α_{1} + 4 α_{2} (1 + E_{22} + E_{33}) & - 8 α_{2} E_{12} & - 8 α_{2} E_{13} \\ - 8 α_{2} E_{12} & 2 α_{1} + 4 α_{2} (1 + E_{11} + E_{33}) & - 8 α_{2} E_{23} \\ - 8 α_{2} E_{13} & - 8 α_{2} E_{23} & 2 α_{1} + 4 α_{2} (1 + E_{11} + E_{22}) \end{matrix}) - - - (17)$ >

弹性张量D表示为：

$> D = (\begin{matrix} D_{1111} & D_{1122} & D_{1133} & D_{1123} & D_{1113} & D_{1112} \\ D_{2211} & D_{2222} & D_{2233} & D_{2223} & D_{2213} & D_{2212} \\ D_{3311} & D_{3322} & D_{3333} & D_{3323} & D_{3313} & D_{3312} \\ D_{2311} & D_{2322} & D_{2333} & D_{2323} & D_{2313} & D_{2312} \\ D_{1311} & D_{1322} & D_{1333} & D_{1323} & D_{1313} & D_{1312} \\ D_{1211} & D_{1222} & D_{1233} & D_{1223} & D_{1213} & D_{1212} \end{matrix}) - - - (18)$ >

式中，D中的元素S_i″j″表示S中的第i″行第j″列的元素，E_k″l″表示E中的第k″行第l″列的元素；

令：

2α₁+4α₂(1+E₂₂+E₃₃)＝G

2α₁+4α₂(1+E₁₁+E₂₂)＝H

2α₁+4α₂(1+E₁₁+E₃₃)＝B

得到对称矩阵D，该矩阵中的对称部分用sym表示：

$> D = Ψ (\begin{matrix} G^{2} & 4 α_{2} + G B & 4 α_{2} + G H & - 8 α_{2} E_{23} G & - 8 α_{2} E_{13} G & - 8 α_{2} E_{12} G \\ B^{2} & 4 α_{2} + H B & - 8 α_{2} E_{23} B & - 8 α_{2} E_{13} B & - 8 α_{2} E_{12} B \\ H^{2} & - 8 α_{2} E_{23} H & - 8 α_{2} E_{13} H & - 8 α_{2} E_{12} H \\ 64 α_{2}^{2} E_{23}^{2} - 8 α_{2} & 64 α_{2}^{2} E_{13} E_{23} & 64 α_{2}^{2} E_{12} E_{23} \\ s y m & 64 α_{2}^{2} E_{13}^{2} - 8 α_{2} & 64 α_{2}^{2} E_{13} E_{12} \\ 64 α_{2}^{2} E_{12}^{2} - 8 α_{2} \end{matrix}) - - - (19)$ >

步骤2.4、建立基于v-p变量的控制方程

本研究中选取速度v、压力p为未知量建立基于v-p变量的有限元方程，进而表达关节软骨基于v-p变量的控制方程；

为了消除流体相的速度，先将式(3)和式(8)带到式(6)中，通过式(4)得到：

$> φ_{f} (v_{f} - v_{s}) = - κ ▿ p - - - (20)$ >

在饱和条件式(1)下，将式(20)带入式(2)中，得到只包含固体相速度的关系式：

$> ▿ \cdot (v_{s} - κ ▿ p) = 0 - - - (21)$ >

接下来由动量守恒方程(5)(6)和本构方程(7)(8)联立可得：

$> ▿ \cdot (σ_{e} - p I) = 0 - - - (22)$ >

由式(21)和式(22)就构成了基于v-p变量的控制方程；

步骤3、建立关节软骨力学平衡方程并进行有限元计算：

步骤3.1、建立边界条件：

流体通量定义为则基于v-p变量的控制方程的边界条件为：

$> u_{s} (x, t) = {\overline{u}}_{s}, x \in Γ_{u_{s}} - - - (23)$ >

$> σ (x, t) \cdot n^{*} = \overline{T}, x \in Γ_{T} - - - (24)$ >

$> p (x, t) = \overline{p}, x \in Γ_{p} - - - (25)$ >

$> - κ ▿ p \cdot n^{*} = \overline{Q}, x \in Γ_{Q} - - - (26)$ >

式中，x指在某一时刻t时当前构型中的空间点；为已知的固相位移，为总的载荷力，为已知的压力，n^*是边界Γ的外法线向量；Γ_p和Γ_Q分别代表液体相体积域Ω_f的边界，和Γ_T分别代表固体相体积域Ω_s的边界，且边界满足条件：

步骤3.2、构造形函数：

定义有限元四面体单元的形函数为：

N＝[N_i′·I_3×3>j′·I_3×3>m′·I_3×3>p′·I_3×3](44)

其中，

$> N_{i^{'}} = \frac{1}{6 V_{e}} (a_{i^{'}} + b_{i^{'}} x^{'} + c_{i^{'}} y^{'} + d_{i^{'}} z^{'}) - - - (45)$ >

式中，V_e是四面体单元的体积，a_i′、b_i′、c_i′、d_i′为与N_i′对应的系数；

N_j′、N_m′、N_p′与N_i′类似，有

$> \begin{matrix} N_{j^{'}} = \frac{1}{6 V_{e}} (a_{j^{'}} + b_{j^{'}} x^{'} + c_{j^{'}} y^{'} + d_{j^{'}} z^{'}) \\ N_{m^{'}} = \frac{1}{6 V_{e}} (a_{m^{'}} + b_{m^{'}} x^{'} + c_{m^{'}} y^{'} + d_{m^{'}} z^{'}) \\ N_{p^{'}} = \frac{1}{6 V_{e}} (a_{p^{'}} + b_{p^{'}} x^{'} + c_{p^{'}} y^{'} + d_{p^{'}} z^{'}) \end{matrix} - - - (46)$ >

步骤3.3、建立关节软骨力学平衡方程：

基于有限元求解的单元类型为四面体单元，构造插值函数；通过基于v-p变量的控制方程(21)、(22)以及边界条件(23)-(26)，采用伽辽金加权余量法建立关节软骨力学平衡方程；

步骤3.4、求解关节软骨力学平衡方程，实现关节软骨两相模型的仿真。

本发明具有以下优点：

1、目前的关节软骨被视为两项模型。但是，这些研究中大都将软骨的固相看作是线弹性材料，线弹性模型只对小变形有效。针对关节软骨在力的施加过程中会产生较大变形的问题，本发明提出基于超弹性固体相特性建立关节软骨两相模型，实现了关节软骨两相模型的建立和仿真。

2、本发明针对关节软骨的结构特点，将关节软骨描述为由不可压缩的非粘性液体相和具有超弹性、横观各向同性的固体相所组成的二相混合物，并且在变形的条件下还会引起软骨组织渗透性的改变。

3、依据Helmholtz应变能函数定义固体相的超弹性特性，液体相定义为理想流体，该模型能够描述关节软骨的非线性、不可压缩性和渗透性。通过伽辽金加权余量法获得基于v-p变量的有限元系统平衡方程，通过有限差分法对平衡方程进行数值计算。

4、本发明针对关节软骨在力的施加过程中产生较大变形的情况，本发明建立的关节软骨两相模型也能够准确的反应实际关节软骨的情况，与实际关节软骨的情况的相比，误差低于13％。

附图说明

图1为从DICOM数据到生成网格模型并获取模型的单元编号和节点坐标的流程示意图。

图2为单轴压缩的应力-速度实验曲线和有限元仿真结果对比图。

具体实施方式

具体实施方式一：

基于超弹性固体相特性关节软骨两相模型的建立方法，包括以下步骤：

步骤1、几何模型和网格划分：

从DICOM数据到生成网格模型并获取模型的单元编号和节点坐标作为形变计算输入的原始数据，获取过程如图1所示；

步骤2、建立基于混合物理论的关节软骨两相模型及基于v-p变量的控制方程，包括以下步骤：

步骤2.1、建立关节软骨的两相模型：

由饱和下的体积分数式得到：

φ_s+φ_f＝1(1)

质量平衡方程为：

$> ▿ \cdot (φ_{s} v_{s} + φ_{f} v_{f}) = 0 - - - (2)$ >

其中，表示矢量微分算子符号；v_s是固体相速度，v_f是液体相速度；

P_s＝-P_f＝K′(v_f-v_s)(3)

式中，K′为扩散阻力系数；

$> K^{'} = \frac{{(φ_{f})}^{2}}{κ} - - - (4)$ >

式中，κ为软骨渗透率；固体相动量为P_s，液体相动量P_f；

令σ_s和σ_f分别为固相和液相的柯西应力张量，关节软骨在当前构形的动量平衡方程为：

$> ▿ \cdot σ_{s} + P_{s} = 0 - - - (5)$ >

$> ▿ \cdot σ_{f} + P_{f} = 0 - - - (6)$ >

由软骨的不可压缩性和力学特性得到关节软骨弹性固体相和非粘性流体相的力学本构方程为：

σ_s＝-φ_spI+σ_e(7)

σ_f＝-φ_fpI(8)

其中，I是单位矩阵，p为压力；σ_e为固体相的有效应力也叫做与固相的变形量相一致的弹性应力张量；

经过以上推导得出由质量平衡方程(2)、动量平衡方程(5)(6)和两相混合物的本构方程(7)(8)组成了关节软骨的两相模型；

根据第二Piola-Kirchhoff应力张量的定义可知第二Piola-Kirchhoff应力张量S与弹性应力张量σ_e的关系为：

σ_e＝J^-1FSF^T＝FSF^T(9)

其中，F为变形梯度矩阵；J为弹性体积比，是变形梯度矩阵F的雅可比行列式；S为第二Piola-Kirchhoff应力张量，该参数在2.3中求解，通过S的联结而得到关节软骨的超弹性两相介质模型；

以初始构形为参考构型，利用物质描述法可得出体元、面元在有限变形中的变化，即当前构形中的单元物质的体积与参考构型中的体积满足以下关系：

dV＝JdV₀(10)

步骤2.2、建立超弹性模型：

关节软骨的超弹性模型由Helmholtz应变能函数来表示，即：

$> Ψ = α_{0} \frac{e^{α_{1} (I_{1} - 3) + α_{2} (I_{2} - 3)}}{{I_{3}}^{β}} - - - (14)$ >

α₀、α₁、α₂、β是材料的力学参数，来源于数学模型与实验数据的拟合结果；I₁、I₂、I₃分别是C的第一、第二、第三主不变量；

对于不可压缩材料而言，I₃＝det(C)＝1，简化后的超弹性模型为：

$> Ψ = α_{0} e^{α_{1} (I_{1} - 3) + α_{2} (I_{2} - 3)} - - - (15)$ >

步骤2.3、求解第二Piola-Kirchhoff应力张量S：

Green应变张量E为：

$> E = (\begin{matrix} E_{11} & E_{12} & E_{13} \\ E_{21} & E_{22} & E_{23} \\ E_{31} & E_{32} & E_{33} \end{matrix}) - - - (16)$ >

弹性张量D表示为：

式中，D中的元素S_i″j″表示S中的第i″行第j″列的元素，E_k″l″表示E中的第k″行第l″列的元素；

令：

2α₁+4α₂(1+E₂₂+E₃₃)＝G

2α₁+4α₂(1+E₁₁+E₂₂)＝H

2α₁+4α₂(1+E₁₁+E₃₃)＝B

得到对称矩阵D，该矩阵中的对称部分用sym表示：

步骤2.4、建立基于v-p变量的控制方程

本研究中选取速度v、压力p为未知量建立基于v-p变量的有限元方程，进而表达关节软骨基于v-p变量的控制方程；

为了消除流体相的速度，先将式(3)和式(8)带到式(6)中，通过式(4)得到：

$> φ_{f} (v_{f} - v_{s}) = - κ ▿ p - - - (20)$ >

在饱和条件式(1)下，将式(20)带入式(2)中，得到只包含固体相速度的关系式：

$> ▿ \cdot (v_{s} - κ ▿ p) = 0 - - - (21)$ >

接下来由动量守恒方程(5)(6)和本构方程(7)(8)联立可得：

$> ▿ \cdot (σ_{e} - p I) = 0 - - - (22)$ >

由式(21)和式(22)就构成了基于v-p变量的控制方程；

步骤3、建立关节软骨力学平衡方程并进行有限元计算：

步骤3.1、建立边界条件：

流体通量定义为则基于v-p变量的控制方程的边界条件为：

$> u_{s} (x, t) = {\overline{u}}_{s}, x \in Γ_{u_{s}} - - - (23)$ >

$> σ (x, t) \cdot n^{*} = \overline{T}, x \in Γ_{T} - - - (24)$ >

$> p (x, t) = \overline{p}, x \in Γ_{p} - - - (25)$ >

$> - κ ▿ p \cdot n^{*} = \overline{Q}, x \in Γ_{Q} - - - (26)$ >

步骤3.2、构造形函数：

定义有限元四面体单元的形函数为：

N＝[N_i′·I_3×3>j′·I_3×3>m′·I_3×3>p′·I_3×3](44)

其中，

$> N_{i^{'}} = \frac{1}{6 V_{e}} (a_{i^{'}} + b_{i^{'}} x^{'} + c_{i^{'}} y^{'} + d_{i^{'}} z^{'}) - - - (45)$ >

式中，V_e是四面体单元的体积，a_i′、b_i′、c_i′、d_i′为与N_i′对应的系数；

N_j′、N_m′、N_p′与N_i′类似，有

步骤3.3、建立关节软骨力学平衡方程：

步骤3.4、求解关节软骨力学平衡方程，实现关节软骨两相模型的仿真。

具体实施方式二：

本实施方式所述步骤3.3建立关节软骨力学平衡方程的具体过程如下：

有限元求解的单元类型为四面体单元，选择形函数作为微分方程和边界条件的权函数，令加权余量积分表达式中的权函数分别为N_J和-N_J，得到相应的加权余量表达式，通过高斯散度定理和总牵引力的定义得到简化后的加权余量积分表达式为：

$> \begin{matrix} \int_{Ω} {▿ {N_{J}}^{T} : σ_{e} - {N_{J}}^{T} ▿ \cdot {N_{J}}^{T} - κ (▿ p \cdot ▿ {N_{J}}^{T})} d Ω \\ = \int_{Γ_{T}} {N_{J}}^{T} \cdot \overline{T} d Γ + \int_{Γ_{Q}} {N_{J}}^{T} \cdot \overline{Q} d Γ \end{matrix} - - - (29)$ >

将问题域离散化成单元形式，对速度和压力进行插值：

v＝Nv_e(30)

p＝Np_e(31)

其中，ν_e代表四面体单元速度，p_e代表四面体单元压力；

则在一个有限元四面体单元网格上的加权余量表达式的矩阵形式为：

$> (\begin{matrix} e_{n} & e_{n} \end{matrix}) {(\begin{matrix} 0 & {G_{e}}^{T} \\ G_{e} & H_{e} \end{matrix}) (\begin{matrix} v_{e} \\ p_{e} \end{matrix}) + (\begin{matrix} M_{e} \\ 0 \end{matrix})} = (\begin{matrix} e_{n} & e_{n} \end{matrix}) (\begin{matrix} T \\ Q \end{matrix}) - - - (32)$ >

式中，e_n为一个单位向量，由于e_n是非零常向量，因此可以直接消掉；

$> G_{e} = \int_{Ω_{0}} - {N_{J}}^{T} ▿ \cdot {NdΩ}_{0}, {G_{e}}^{T} = \int_{Ω_{0}} - (▿ \cdot N) {N_{J}}^{T} {dΩ}_{0}, H_{e} = \int_{Ω_{0}} - ▿ \cdot {N_{J}}^{T} κ ▿ \cdot {NdΩ}_{0},$ >

$> M_{e} = \int_{Ω_{0}} t r [{(▿_{0} N_{J})}^{T} F S] {dΩ}_{0}, T = \int_{Γ_{T_{0}}} {N_{J}}^{T} \cdot \overline{T} J_{Γ} {dΓ}_{0}, Q = \int_{Γ_{Q 0}} {N_{J}}^{T} \cdot \overline{Q} J_{Γ} {dΓ}_{0};$ >

$> ▿ {N_{J}}^{T} : {FSF}^{T} = t r [{(▿_{0} N_{J})}^{T} F S], ▿_{0} = ▿ \cdot F;$ >

令：

$> Y = (\begin{matrix} 0 & {G_{e}}^{T} \\ G_{e} & H_{e} \end{matrix}), v = (\begin{matrix} v_{e} \\ p_{e} \end{matrix}), M = (\begin{matrix} M_{e} \\ 0 \end{matrix}), f = (\begin{matrix} T \\ Q \end{matrix}) - - - (33)$ >

得到最终的关节软骨力学平衡方程：

Yv+M＝f(34)

式中，Y是一个容量矩阵，M是在固体相中与弹性应力张量相关的节点力的非线性弹性向量，f是外加载荷力向量。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：

本实施方式所述步骤3.4所述的求解关节软骨力学平衡方程的具体过程如下：

通过有限差分法对关节软骨力学平衡方程(34)进行数值计算：

首先在时域上采用增量法，将时间t离散化为若干个时间点，即：

0＝t₀<t₁<…<t_n<t_n+1<…<t_N＝t(35)

时间增量表示为：

△t_n+1＝t_n+1-t_n(36)

采用完全拉格朗日法，令初始时刻t₀＝0的构形作为参考构形，并且该参考构形不随时间的变化而发生变化；在t_n+1时刻式(34)写成：

Y_n+1v_n+1+M_n+1＝f_n+1(37)

通过Newton-Raphson方法，M的线性化形式由以下递推形式获得：

$> \begin{matrix} M_{n + 1}^{i} = (\begin{matrix} M_{n + 1}^{i} \\ 0 \end{matrix}), \\ M_{n + 1}^{i} = M_{n + 1}^{i - 1} + K_{n + 1}^{i - 1} \cdot u_{n + 1}^{i} \end{matrix} - - - (38)$ >

式中，K_n+1是固相的切线刚度矩阵，K_n+1是的变分；i代表迭代次数；

由梯形法确定迭代位移的递推关系为：

$> {Δu}_{n + 1}^{i} = u_{n + 1}^{i} - u_{n + 1}^{i - 1} = {ωΔtΔv}_{n + 1}^{i} - - - (39)$ >

式中，u为节点位移，为t_n+1时刻第i次迭代的节点位移；ω是指定的时间积分参数，0≤ω≤1，且当ω≥1/2时，为隐式积分，方程无条件稳定；

速度的增量表示为：

$> v_{n + 1}^{i} = v_{n + 1}^{i - 1} + {Δv}_{n + 1}^{i} - - - (40)$ >

在速度的第i次迭代中，方程(34)可写成：

$> Y_{n + 1}^{i - 1} v_{n + 1}^{i} + M_{n + 1}^{i} = f_{n + 1}^{i - 1} - - - (41)$ >

由(38)(39)(40)(41)能够得到：

$> (Y_{n + 1}^{i - 1} + {ωΔtK}_{n + 1}^{i - 1}) {Δv}_{n + 1}^{i} = f_{n + 1}^{i - 1} - Y_{n + 1}^{i - 1} v_{n + 1}^{i - 1} - M_{n + 1}^{i - 1} - - - (42)$ >

在△t_n+1时间步长内，对式(42)反复进行迭代，直到速度增量满足收敛准则；在一个时间步长内，迭代的收敛准则为：

$> \frac{| {Δv}^{i} |}{| v^{i} |} < T O L, T O L = 10^{- 4} - - - (43)$ >

在当前步长内迭代收敛发生时，将得到的解作为初始值更新到下一个时间步长内，并且进行重复迭代直到收敛再次发生；反复迭代并进行数值计算，最终实现关节软骨两相模型的仿真。反复迭代并进行数值计算的过程在Visual Studio平台C语言环境下编程实现。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：

本实施方式步骤3.4中所述的反复迭代并进行数值计算的迭代求解流程如下：

(a)在△t＝t_n+1-t_n增量上定义已知初始速度v_n和初始位移u_n，且有

$> v_{n + 1}^{0} = v_{n}, u_{n + 1}^{0} = u_{n} + {Δtv}_{n}$ >

(b)设置迭代次数：i＝1

(c)使用来计算方程(42)中的矩阵和其余未知量，解算这个线性代数方程得到第一个速度增量

$> (Y_{n + 1}^{0} + {ωΔtK}_{n + 1}^{0}) {Δv}_{n + 1}^{1} = f_{n + 1}^{0} - Y_{n + 1}^{0} v_{n + 1}^{0} - M_{n + 1}^{0}$ >

(d)结束第一次迭代，更新位移和速度向量：

$> v_{n + 1}^{1} = v_{n + 1}^{0} + {Δv}_{n + 1}^{1}, u_{n + 1}^{1} = u_{n + 1}^{0} + {ωΔtΔv}_{n + 1}^{1}$ >

(e)更新迭代次数到i

(f)使用更新矩阵和向量，解线性方程得到速度增量

$> (Y_{n + 1}^{i - 1} + {ωΔtK}_{n + 1}^{i - 1}) {Δv}_{n + 1}^{i} = f_{n + 1}^{i - 1} - Y_{n + 1}^{i - 1} v_{n + 1}^{i - 1} - M_{n + 1}^{i - 1}$ >

(g)结束当前迭代，更新速度和位移向量：

$> v_{n + 1}^{i} = v_{n + 1}^{i - 1} + {Δv}_{n + 1}^{i}, u_{n + 1}^{i} = u_{n + 1}^{i - 1} + {ωΔtΔv}_{n + 1}^{i}$ >

回到步骤(e)继续进行新的迭代，直到速度的增量满足收敛准则。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：

本实施方式步骤3.2中所述的V_e、a_i′、b_i′、c_i′、d_i′的具体形式如下

$> V_{e} = \frac{1}{6} (\begin{matrix} 1 & x_{i^{'}}^{'} & y_{i^{'}}^{'} & z_{i^{'}}^{'} \\ 1 & x_{j^{'}}^{'} & y_{j^{'}}^{'} & z_{j^{'}}^{'} \\ 1 & x_{m^{'}}^{'} & y_{m^{'}}^{'} & z_{m^{'}}^{'} \\ 1 & x_{p^{'}}^{'} & y_{p^{'}}^{'} & z_{p^{'}}^{'} \end{matrix}) - - - (47)$ >

$> \begin{matrix} a_{i} = (\begin{matrix} x_{j^{'}}^{'} & y_{j^{'}}^{'} & z_{j^{'}}^{'} \\ x_{m^{'}}^{'} & y_{m^{'}}^{'} & z_{m^{'}}^{'} \\ x_{p^{'}}^{'} & y_{p^{'}}^{'} & z_{p^{'}}^{'} \end{matrix}), b_{i} = - (\begin{matrix} 1 & y_{j^{'}}^{'} & z_{j^{'}}^{'} \\ 1 & y_{m^{'}}^{'} & z_{m^{'}}^{'} \\ 1 & y_{p^{'}}^{'} & z_{p^{'}}^{'} \end{matrix}) \\ c_{i} = - (\begin{matrix} x_{j^{'}}^{'} & 1 & z_{j^{'}}^{'} \\ x_{m^{'}}^{'} & 1 & z_{m^{'}}^{'} \\ x_{p^{'}}^{'} & 1 & z_{p^{'}}^{'} \end{matrix}), d_{i} = - (\begin{matrix} x_{j^{'}}^{'} & y_{j^{'}}^{'} & 1 \\ x_{m^{'}}^{'} & y_{m^{'}}^{'} & 1 \\ x_{p^{'}}^{'} & y_{p^{'}}^{'} & 1 \end{matrix}) \end{matrix} - - - (48)$ >

式中，(x′_i′，y′_i′，z′_i′)、(x′_j′，y′_j′，z′_j′)、(x′_m′，y′_m′，z′_m′)、(x′_p′，y′_p′，z′_p′)分别为四面体单元的四个顶点坐标。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：

本实施方式步骤2.1所述的弹性体积比J＝1，由于关节软骨被视为不可压缩成型材料，所以J＝1。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：

本实施方式步骤2.1中所述的软骨渗透率其中L是常数，下标“0”代表未变形前的初始配置。关节软骨不仅具有超弹性，而且在变形的条件下还会引起组织渗透性的改变。

其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：

本实施方式所述中的φ_s0＝0.2，κ₀＝2.519×10^-15m⁴/(Ns)，L＝0.0848。

其它步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：

本实施方式步骤2.2中所述的α₀＝0.1084Mpa，α₁＝0.592Mpa，α₂＝0.0846Mpa。

其它步骤及参数与具体实施方式一至八之一相同。

本发明针对关节软骨在力的施加过程中产生较大变形的情况，本发明建立的关节软骨两相模型也能够准确的反应实际关节软骨的情况。如图2所示，通过与相同条件下的实验结果对比，有限元分析结果在数值和趋势上都与实验数据曲线吻合，图中“实验结果”为实际关节软骨的实验数据，“仿真结果”为本发明的建立的基于超弹性固体相特性关节软骨两相模型的仿真数据。通过仿真结果与实验曲线的对比发现：基于超弹性固体相特性关节软骨两相模型能够精确描述关节软骨的力学特性，与实际关节软骨的情况的相比，误差低于13％，验证了模型的精确性和有限元程序的有效性。

去获取专利，查看全文>

相似文献

专利
中文文献
外文文献

1. 基于超弹性固体相特性关节软骨两相模型的建立方法 [P] . 中国专利： CN106202738B . 2017.12.19
2. 基于超弹性固体相特性关节软骨两相模型的建立方法 [P] . 中国专利： CN106202738A . 2016-12-07
3. DEVICE FOR MEASURING THE PHASE CONCENTRATION OF A TWO-PHASE FLOW, A SYSTEM AND A METHOD FOR MEASURING FLOW PARAMETERS OF A TWO-PHASE FLOW, CAPABLE OF DETECTING THE VOLUME CONCENTRATION OF A SOLID PHASE, A FLOW RATE OF A TWO PHASE FLOW, AND A MASS FLOW RATE OF WORKING MEDIUM [P] . 韩国专利： KR20100051037A . 2010-05-14

机译：一种用于测量两相流的相密度的装置，一种用于测量两相流的相量的系统和方法，其能够检测固体相的体积浓度，两相流的流速以及工作介质的质量流量
4. Direct production process of a length of continuous thin two- phase stainless steel strip having excellent superplasticity and surface properties [P] . 美国专利： US5259443A . 1993-11-09

机译：一段具有超塑性和表面特性的连续薄两相不锈钢薄带的直接生产工艺
5. Three-phase regulator calculating method for controlling asynchronous motor, involves generating single input single output model from two-phase Park model, and calculating three-phase regulator from single input single output regulator [P] . 法国专利： FR2922061A1 . 2009-04-10

机译：用于控制异步电动机的三相调节器计算方法，包括从两相Park模型生成单输入单输出模型，以及从单输入单输出调节器计算三相调节器