骨重建的有限元计算模拟：小位移弹性损伤理论

摘要

背景

骨骼最重要的特性是通过自我修复和适应环境条件来抵抗损伤。这些适应性变化受通常称为骨重塑的生理过程调节。更好地理解这一过程需要我们将小位移假设下的弹性损伤理论应用于骨骼结构，并观察其力学行为。

结果

本研究的目的是通过考虑弹性损伤和机械刺激来模拟股骨近端的二维模型。在这里，我们提出了一个基于非线性常微分方程组的数学模型，并发展了力学问题的变分公式。然后，我们将我们的数学模型实现到有限元算法中，以研究损伤的影响。

结论

结果与现有文献一致，现有文献表明，在机械载荷作用下，受损骨结构的骨刚度下降。

介绍

骨骼是骨骼系统的主要组成部分，能够基本保持身体的形状；保护内脏；储存矿物质和油脂；参与血细胞生产；并通过将肌肉收缩力从一个部位传递到另一个部位来帮助身体运动[1].

作为一种活组织，骨骼能够在外力的影响下通过重新分配密度来优化其结构。自从Wolff发表这篇文章以来，人们提出了许多描述骨密度重新分布的理论[2–4].

这一过程称为骨重塑，后来由Huiskes等人正式开发，使用的概念是，骨重塑是由激活调节细胞并导致局部骨适应的局部机械信号诱导的[5,6].

然而，当外力高于临界水平时，骨骼可能会因增加通常修复的损伤形成而更容易骨折[7,8].

因此，更好地了解骨重塑过程有助于预防骨折和其他类型的疾病。已经开展了多项工作，将骨重塑过程与数学观点联系起来，然后进行一些计算模拟[6,7].

本研究的总体目标是使用小位移的弹性损伤理论对股骨近端进行数值模拟。首先，我们描述了力学问题，并推导了它的变分公式。接下来，我们提出了一种骨重建算法，并使用有限元方法解决了二维股骨问题。

最后，给出了一些二维计算模拟，结果与文献报道的结果完全一致。

背景

几何图形和材料特性

我们将股骨近端的二维模型视为文献中的先前研究[9–11]. 建议的几何形状如图所示1.

图1

股骨近端模型。

股骨的材料性质假定为线弹性、可损伤、各向同性和均匀。泊松比指定给υ = 0.3参考文献[12,13]; 杨氏模量（E）由骨密度ρ计算得出，表达式如下：

其中M和γ是正常数。

弹性损伤骨重建理论

股骨可以根据外部载荷进行重塑，以在承受较高应力的区域提供更大的骨基质。这些外部载荷会引起力学场的变化和股骨的损伤，例如宏观裂纹或微观裂纹[14,15].

这种损害概念是由科学界，特别是卡恰诺夫在20世纪90年代提出的[16]. 本文中使用的弹性损伤模型最初是由法国学派开发的，尤其是Jean Lemaitre的模型[17,18]. 它通过引入量化微裂纹影响的标量变量D来描述材料的本构行为。

在这一弹性损伤力学理论中，单元的弹性模量可能随着损伤的发展而逐渐降低。在小位移假设下，受损材料的弹性模量定义如下：

哪里：

D是0≤D≤1的损伤程度

E是未受损弹性的杨氏模量

是损伤弹性的实际模量。然后将损伤表示为刚度损失。

我们将弹性损伤力学理论应用于骨重建过程，并使用前面的方程得出：

在哪里？是考虑到Abdali建议的损伤的骨密度[19]，可表示为：

我们假设此函数有界为：

哪里：

ρ_最小值与再吸收骨相对应的最小密度

ρ_最大值皮质骨的最大密度

此外，设ρ₀表示初始骨密度

然后，我们介绍了Martin提出的损伤演化规律[7,20]通过方程式：

哪里：

D类₀是最初的损坏

t是时间

（f）_d日骨骼的疲劳寿命是否没有重塑

控制方程

让，是中的非空开放有界域具有Lipschitz连续边界Γ=⏴Ω。边界被分割成两个不相交的部分Γ₁和Γ₂其中Γ₁是股骨固定的边界的固定部分，Γ₂也是施加力F（F1，F2）的边界的一部分。

我们在下面的任何地方使用表示二阶对称张量的空间，并用n表示Γ的单位外法向量。为了简化，我们认为占据集合的身体不受密度f的体积力作用。建模后，定义了荷载和边界条件。

图2显示了作用在股骨上的不同载荷位置。

图2

不同的荷载位置，（a）荷载工况1，（b）荷载工况2。

最后，让我们是位移场，应力场，应变张量，骨密度函数和T>0的持续时间。

我们记得ϵ（u）是由[9,21,22]:

在这里∇代表渐变操作符。

本构关系与应力-应变关系相关[9]写为：

哪里：

Div表示散度运算符

l表示中的恒等运算符

是拉梅系数，假定该系数取决于以下所示的骨密度

表示为

表示为

具有：

杨氏模量

泊松比

世界上已经发表了许多骨重塑定律，主要针对骨密度的演变[5,6,20]. 我们采用Huiskes等人提出的定律，他们使用应变能密度作为刺激信号来控制骨重塑过程[5,6,20,22].

骨密度函数的演变[6,22]由以下非线性一阶常微分方程得到：

哪里：

B、 s和K是实验常数

U（σ（U），\1013；（U））是应变能密度[9]调节重塑过程：

设：是两个张量的双收缩，产生标量。

最后，问题的方程组定义如下。

问题：查找这样的话

在哪里？

变分公式

该模型的变分公式包含位移场的变分方程。

（等式1、3和4）的弱形式如下：我们求位移场u∈ V={φ∈ 【H】¹(Ω)]²; Γ上φ=0₁}这样的话

其中v是测试函数。

时间离散化

龙格-库塔二阶方法是一种用于求解常微分方程的数值技术2:

使用：

∆t是时间步长

是初始数据

提出的算法

采用不同的算法将材料属性分配到网格模型的元素中[23–25]. 其中一些算法是基于无损伤的假设。本研究的目的是用有限元方法模拟股骨的弹性损伤。分层算法的示意图如图所示三.

图3

提出的骨重塑算法的示意图。

提出的算法可以概括为以下七个步骤：

• 步骤1。定义全局模型：几何体、载荷条件和初始骨密度分布。考虑了具有均匀密度分布ρ的初始模型的重塑₀=0.8克/厘米^三.

• 第2步。确定杨氏模量、泊松比和拉梅系数。

• 步骤3。通过求解位移场的线性变分方程来计算位移。

• 步骤4。使用有限元方法评估每个离散位置的应变能密度U。

• 步骤5。证明机械刺激是否会导致骨并置、骨吸收或平衡。然后，更新骨骼密度。

• 步骤6。更新损坏情况。

• 第7步。检查收敛情况。收敛准则是根据迭代过程中质量的变化而设定的。当满足收敛准则时，得到最终拓扑；否则，迭代过程将从步骤2继续。

结果和讨论

众所周知，研究骨的生物力学行为具有重要意义。本工作旨在模拟弹性损伤股骨，以提供有关骨骼几何拓扑和材料特性的有用信息。

以下数据[7,9,19,20,22,24]受雇于：

ρ_最小值 = 0.01克/厘米^三, ρ_最大值 = 1.74克/厘米^三，K=0.004焦耳/克，s=0.1，ρ₀ = 0.8克/厘米^三，γ=3，B=1（克/厘米^三)²（MPa.UT）^- 1，D₀ = 0.8，f_d日 = 3年，f=0 N/m²，Δt=10^- 5美国犹他州，υ = 0.3，M=3790（兆帕）（克/厘米^三)^- 3，（F₁（情况1）=1000 N，F₂（情况1）=1200 N），（F₁（情况2）=1200 N，F₂（情况2）=1500 N）。

假设域中存在弹性损伤，就机械应力期间股骨近端的骨重塑进行了几个计算模拟。同时，我们将稳定力和固定约束作为该模型的边界条件。前面描述的拟议算法在FreeFEM++中实现（参见[26])，针对不同的负载情况执行多个仿真。考虑了两种载荷情况来评估它们对骨密度分布的影响。

图4对荷载工况1中无损伤和有损伤的骨密度分布进行了比较。

图4

载荷工况1下从初始时间（左侧）到最终时间（右侧）的骨密度，（a）未受损股骨，（b）受损股骨。

如图所示4（a） 是从最初到最后时间的未受损骨密度分布。可以看出，机械载荷触发了骨重塑过程，尤其是在靠近载荷的区域。与其他区域相比，该区域的骨密度较高。

图4（b） 绘制了从最初到最后一段时间内受损骨骼密度分布的变化。它导致密度在初始时间内增加，然后逐渐降低。但它永远不会达到未受损骨骼的密度。

在图中5，示出了未受损和受损骨骼在载荷情况2下的密度分布。也可以观察到与上图相同的趋势。

图5

载荷工况2下从初始时间（左侧）到最终时间（右侧）的骨密度，（a）未受损股骨，（b）受损股骨。

研究结果表明，在接近机械负荷的区域，骨密度高于正常值；在远离机械负荷的区域，骨密度低于其他区域。结果与Li等人获得的结果相似[23]，Sharma等人[24]和Li等人[27].

比较数字4和5可以观察到，骨密度的降低进一步导致骨刚度的杨氏模量降低，从而增加骨折的风险。在Tomaszewski等人的工作中也观察到了相同的结果[28].

从这一比较中，我们发现，在骨骼中，受损结构的密度在最初的时候降低，然后在最后的时候可以在一定程度上修复损伤本身。但只有当负载不太高，以至于自我修复机制能够跟上不断增加的损伤时，才会发生这种情况[7,14,29,30].

在这种情况下，尽管存在许多生物因素，如免疫、激素和血流动力学刺激，但仅考虑机械刺激的自我修复机制；从而因个人而异[24,27,31].

本文介绍的工作可以应用于不同的模型以及骨科生物材料的研究，并有助于进一步的研究。

结论

本研究采用二维数学模型和基于有限元的数值技术，建立了一个模拟弹性损伤骨重建的模型。研究了应变和损伤对骨结构的影响。

本文的结果表明，在请求区域，骨密度很重要；并使我们能够观察到与文献研究结果的良好一致性。

因此，从生物力学角度来看，最好使用有限元方法模拟三维股骨，以便更好地了解骨的行为。