组织液-纤维环流固耦合的椎间盘软组织损伤演变预测方法技术

本发明专利技术公开了一种组织液

【技术实现步骤摘要】
组织液
‑
纤维环流固耦合的椎间盘软组织损伤演变预测方法


[0001]本专利技术属于骨科生物力学领域，涉及考虑椎间盘组织液机械作用的基础上的组织液
‑
纤维环流固耦合的椎间盘软组织损伤演变预测方法。

技术介绍

[0002]椎间盘由中央髓核、纤维环和上下软骨终板组成，可以吸收缓冲腰椎活动时的机械荷载。其中纤维环绕髓核连接上下软骨终板，主要由蛋白聚糖和多层的定向胶原纤维网格构成，为椎间盘提供支撑并实现高拉伸强度，是椎间盘的主要承重部件。
[0003]随着年龄的增长，椎间盘内部会发生微结构和化学成份的变化，这些变化会影响其承受日常载荷的能力。椎间盘发生损伤退变是引起下腰痛的主要原因，由于体外尸体试验并不能完全反映椎间盘在人体内的机械环境，有限元被用于研究人体腰椎间盘损伤的一种新方法，通过数值计算得到椎间盘应力、应变和孔隙压力的分布，并与体外试验互相验证。然而现有的椎间盘损伤模型都存在一定局限性：
[0004]首先是模型忽略了组织液的流入流出，因此难以捕捉到椎间盘的蠕变行为；
[0005]其次，采用不可压缩的单一固相模型不能反应出液相(组织液)的机械贡献。

技术实现思路

[0006]针对以上模型的局限性，本专利技术给出了一个新的解决方案。在椎间盘多孔弹性模型的基础上，提出纤维环的多孔
‑
损伤数值计算模型，即将椎间盘纤维环视为由不可压缩的液体和可压缩的各向异性固相组成。新的建模思路将多孔模型和超弹性各向异性损伤模型相结合，多孔模型部分采用传统的渗流
‑
应力耦合分析，得到椎间盘纤维环在固相有效应力下的真实损伤演化，从而可以更准确的描述椎间盘纤维环生物力学损伤响应。利用这种方法可以捕捉到椎间盘组织的蠕变行为，并在空间上和时间上对特定工况下的椎间盘纤维环的胶原纤维损伤进程进行预测。
[0007]为了实现上述技术目的，本专利技术采用如下技术手段：
[0008]一种组织液
‑
纤维环流固耦合的椎间盘软组织损伤演变预测方法，结合多孔模型和超弹性各向异性损伤模型对椎间盘的力学响应进行数值模拟，实际模拟计算过程通过ABAQUS有限元软件实现，该方法包括以下步骤：
[0009]1)在每一计算步开始时，根据渗流
‑
应力耦合分析方法，计算液相压力p以及液相流动带来的单元体积比J；将有限元积分点处的变形梯度F解耦为体积变化部分F
vol
和等容变化部分用解耦后的变形梯度分别计算固相体积部分应变能密度Ψ
vol
和各向同性等容应变能密度Ψ
iso
；
[0010]2)用子程序UEXTERNALDB导入外部纤维数据，储存在纤维方向数组中，在单元计算时，与纤维方向数组中的单元编号进行匹配，读取当前单元的两族胶原纤维方向并计算未损伤的各向异性应变能密度函数
[0011]3)根据当前未损伤的各向异性应变能密度函数对损伤准则进行判断，计算纤维环中的各向异性损伤系数D
i
，并基于考虑损伤的固相应变能密度函数，在空间参考下计算考虑损伤的固相柯西应力σ
e
以及刚度张量
[0012]4)根据应力和刚度张量公式编写UMAT子程序来表征纤维环的材料损伤行为，对特定工况下某一时间点的椎间盘损伤位置和损伤值进行预测。
[0013]所述步骤1)中渗流
‑
应力耦合分析的计算利用ABAQUS有限元软件中的Soil分析步，将总柯西应力σ看成是固相有效应力σ
e
和液相压力pI之和,其中 I为二阶单位张量。并根据Darcy定理由孔隙压力梯度计算组织液的等效净流速场 v，从而得到每点处液体流出导致的单元体积变化：
[0014]σ＝σ
e
+pI
[0015]其中，μ
f
为动力粘度系数；v是等效净流速场，表示液体相对固相的流动速度，k为应变依赖的渗透率系数，由初始孔隙比e0和初始渗透率k0确定：
[0016][0017]其中M为修正参数，根据试验数据拟合确定，e是现时孔隙比，与单元体积比的关系为：J＝(1+e)/(1+e0)。
[0018]所述步骤1)中，高斯积分点处的变形梯度解耦公式为其中F
vol
＝J
1/3
I，应变能密度的各向同性部分采用Neo
‑
Hooken本构模型：
[0019][0020]其中D1和C
10
是材料参数，是等容部分右柯西
‑
格林变形张量的第一不变量，J是体积比，表示现时体元大小与初始体元大小之比，同时也是变形梯度F的行列式值。
[0021]所述步骤2)中，未损伤纤维环的各向异性应变能密度函数采用 Gasser
‑
Holzapfel
‑
Ogend模型进行描述，两族胶原纤维的初始方向分别由两个纤维方向向量m0和n0表示，未损伤的各项异性部分应变能密度函数由结构张量表示，未损伤的各项异性部分应变能密度函数由结构张量和确定：
[0022][0023]其中，分别是等容部分的右柯西
‑
格林变形张量的第四、六不变量；
[0024]k1表示纤维模量，
[0025]k2是指数系数，为无量纲量；
[0026]损伤折减后的固相各向异性应变能密度函数为：
[0027][0028]所述步骤2)中，确定两族胶原纤维方向的具体流程为：
[0029]21)首先将划分的10层环状基质网格的外表面提取出来，得到共10层二维网格面，输出二维网格面以及对应层基质单元的inp文件，为后续胶原纤维方向的确定作准备；
[0030]22)处理inp文件得到每个单元在对应面上的两个对角线方向，将该单元的对角线方向作为纤维环在该单元处的胶原纤维方向，作为外部数据储存在新的inp文件中，第一列代表单元编号，后续3列是该单元处胶原纤维方向在空间上的x,y,z分量；
[0031]23)在UMAT子程序中通过UEXTERNAL子程序对纤维方向进行读取，作为全局变量储存在数组中，在计算流程中，对单元编号进行识别，提取数组中的纤维方向并进行归一化。
[0032]步骤3)中所述纤维环中的各向异性损伤系数D
i
的计算具体流程为：
[0033]31)根据等容变形梯度计算当前步下各向异性部分的
[0034]32)计算损伤判断准则中的指标并与历史最大指标值进行比较，若损伤面发生移动，根据下式计算损伤折减系数以及纤维环中的各向异性损伤系数D
i
[0035][0036][0037]其中，β
f
是胶原纤维的损伤演化规律控制参数；和分别表示损伤初始时刻和完全损伤时刻的对应损伤指标。
[0038]所述步骤3)中固相有效应力是基于应变能密度函数计算的，根据有限变形理论，在材料参考系下，纤维环的PK2固相应力S...

【技术保护点】

【技术特征摘要】
1.一种组织液
‑
纤维环流固耦合的椎间盘软组织损伤演变预测方法，其特征在于，结合多孔模型和超弹性各向异性损伤模型对椎间盘的力学响应进行数值模拟，实际模拟计算过程通过ABAQUS有限元软件实现，该方法包括以下步骤：1)在每一计算步开始时，根据渗流
‑
应力耦合分析方法，计算液相压力p以及液相流动带来的单元体积比J；将有限元积分点处的变形梯度F解耦为体积变化部分F
vol
和等容变化部分用解耦后的变形梯度分别计算固相体积部分应变能密度Ψ
vol
和各向同性等容应变能密度Ψ
iso
；2)用子程序UEXTERNALDB导入外部纤维数据，储存在纤维方向数组中，在单元计算时，与纤维方向数组中的单元编号进行匹配，读取当前单元的两族胶原纤维方向并计算未损伤的各向异性应变能密度函数3)根据当前未损伤的各向异性应变能密度函数对损伤准则进行判断，计算纤维环中的各向异性损伤系数D
i
，并基于考虑损伤的固相应变能密度函数，在空间参考下计算考虑损伤的固相柯西应力σ
e
以及刚度张量4)根据应力和刚度张量公式编写UMAT子程序来表征纤维环的材料损伤行为，对特定工况下某一时间点的椎间盘损伤位置和损伤值进行预测。2.根据权利要求1所述组织液
‑
纤维环流固耦合的椎间盘软组织损伤演变预测方法，其特征在于，所述步骤1)中渗流
‑
应力耦合分析的计算利用ABAQUS有限元软件中的Soil分析步，将总柯西应力σ看成是固相有效应力σ
e
和液相压力pI之和，其中I为二阶单位张量，并根据Darcy定理由孔隙压力梯度计算组织液的等效净流速场v，从而得到每点处液体流出导致的单元体积变化：σ＝σ
e
+pI其中，μ
f
为动力粘度系数；v是等效净流速场，表示液体相对固相的流动速度，k为应变依赖的渗透率系数，由初始孔隙比e0和初始渗透率k0确定：其中M为修正参数，根据试验数据拟合确定，e是现时孔隙比，与体积比的关系为：J＝(1+e)/(1+e0)。3.根据权利要求1所述组织液
‑
纤维环流固耦合的椎间盘软组织损伤演变预测方法，其特征在于，所述步骤1)中，高斯积分点处的变形梯度解耦公式为其中F
vol
＝J
1/3
I，应变能密度的各向同性部分采用Neo
‑
Hooken本构模型：其中，D1和C
10
是材料参数；是等容部分右柯西
‑
格林变形张量的第一不变量；J是体积比，表示现时体元大小与初始体元大小之比，同时也是变形梯度F的行列式值。4.根据权利要求1所述组织液
‑
纤维环流固耦合的椎间盘软组织损伤演变预测方法，其特征在于，所述步骤2)中，未损伤纤维环的各向异性应变能密度函数采用Gasser
‑
Holzapfel
‑
Ogend模型进行描述，两族胶原纤维的初始方向分别由两个纤维方向向量m0和n0表...

【专利技术属性】
技术研发人员：郭小明，许万玟，
申请(专利权)人：东南大学，
类型：发明
国别省市：