摘要 目的 针对传统有限元法(finite element method,FEM)分析全髋关节置换(total hip arthroplasty,THA)后压电股骨重建时精度低的问题,采用边光滑有限元法(edge-based smoothed finite element method,ES-FEM)对植入假体后压电股骨近端的骨重建进行仿真分析。方法 根据自适应骨重建理论,建立假体-压电股骨模型。基于模型的背景网格构建光滑域,引入梯度光滑技术,求解出光滑的重建刺激,进而得到术后压电股骨近端的密度分布。结果 植入假体后,受力点由股骨头转移到假体,出现应力屏蔽现象,股骨内部表观密度的分布发生明显变化。相比于FEM,ES-FEM在一定程度上能软化数值模型,提高仿真精度。在相同的网格下,电势和密度的求解精度分别提高27%和30%左右。结论 采用ES-FEM能够更精确地模拟出THA术后压电股骨近端的骨重建进程,为THA临床研究提供有效的理论依据。
关键词:
股骨近端
全髋关节置换
自适应骨重建
压电效应
边光滑有限元法
人体的骨骼系统由骨骼、关节和韧带构成,其中关节的作用是实现相邻骨骼的相对运动并传递载荷 [1-2] 。当关节的功能因疾病或创伤受损时,不仅整个骨骼的功能受到破坏,患者的正常生活也会受到重大影响 [3] 。多数情况下,当关节表面发生不可逆转的变化或关节功能严重受损时,更换关节是使患者恢复正常且疼痛感较低的有效方式之一 [4] 。在人工关节置换术中,全髋关节置换十分常见。然而在手术后,假体和股骨同时承担了原本只有股骨承担的外部载荷,导致骨骼承受的应力减小,因而出现了应力屏蔽现象。根据Wolff定律,应力屏蔽将会导致骨骼产生自适应骨重建,进而引起假体的松动 [5] 。因此,术后股骨重建过程的精确模拟,在临床研究中具有重要意义。目前常用的研究方法为数值仿真法,这是一种低成本、高效率,并且对患者友好的方法 [6-7] 。
有限元法(finite element method, FEM)作为一种应用广泛的数值方法,常被用于分析假体植入对骨重建的影响 [8] 。井野等 [9] 模拟健康股骨和THA术后股骨的骨重建过程,通过对比股骨应力分布与密度分布的变化,研究假体植入对股骨重建的影响。Gubaua等 [10] 建立三维人体股骨重建模型,基于细胞动力学,模拟假体植入后与应力屏蔽相关的骨骼适应现象。Levadnyi等 [11] 采用传统FEM,评估不同假体模型对骨密度分布的影响。虽然FEM在THA术后骨重建的模拟中得到广泛应用,但是其本身存在一些缺陷:① FEM本身的刚度“过硬”,导致其求解精度并不高;② FEM求解精度严重依赖于网格尺寸,尺寸越小的网格得到的求解精度越高,但这也极大地增加了工作量;③ 上述研究在进行骨重建分析时,并没有考虑骨具有压电效应这一特性,而由于上述缺陷,FEM在压电问题上的计算精度也不高 [12-14] 。针对上述问题,许多新的方法被相继提出。有研究团队在FEM的基础上引入梯度光滑技术,提出基于边光滑有限元法(edge-based smoothed finite element method,ES-FEM),使模型刚度得到有效软化,从而提高计算精度 [15-16] 。基于ES-FEM,有研究团队探讨静态下的二维压电结构,获得理想的系统刚度,并取得了良好的计算结果 [17] 。
基于上述研究,本文建立了假体-压电股骨模型,利用ES-FEM求解该模型的离散系统方程,从而得到光滑的应变能密度与电能密度。随后基于自适应骨重建理论,模拟植入假体后股骨近端骨重建的过程。结果表明,相比于FEM,ES-FEM可以得到更精确的仿真结果,为全髋关节置换术后压电股骨近端骨重建过程的研究提供有效参考。
骨骼是一个动态结构,当骨骼中的机械环境发生改变(例如出现植入物破坏内部骨形成和吸收的平衡时),骨骼就会经历适应性重建过程,导致骨表观密度的分布发生变化 [18] 。根据自适应骨重建理论,以能量密度 U 作为刺激信号控制重建过程,得到股骨内部表观密度的变化规律 [19-20] :
式中: ρ 为骨的表观密度; g 为实验常数; B 为重构系数; ω 为阈值系数,即重建激励必须偏离参考状态超过一定范围时才会触发自适应过程。
式中: σ 为应力; ε 为应变; D 和 E 分别为电位移和电场强度。
股骨的压电效应通常被描述为力场和电场在本构关系中的耦合,这种耦合关系由密度函数 ρ/ρ * 进行调节,具体可以表示为:
式中: ρ * 为表观密度的参考值; e 表示压电常数矩阵; θ 为介电常数矩阵; C E 为弹性矩阵,可由泊松比 ν 和弹性模量 E 求得。在压电骨中, ν =0 . 3是常数,与密度无关; E 与密度有关,随着密度的更新而更新,则
骨的压电效应主要源于内部的胶原纤维,其等价于六边形对称的晶体 [21] 。因此,股骨的压电常数矩阵 e 和介电常数矩阵 θ 具体表示为:
假设重建过程在[0, T ]的时间范围内进行,时间步长为Δ t 。利用一阶Adams-Bashforth法求解式(1),可得到第 n 个迭代步中密度:
对于传统的FEM,由于其本身固有缺陷导致求解式(2)时计算精度较低,而ES-FEM对系统的梯度场进行光滑处理,可以得到更精确的重建刺激,进而得到更接近实际情况的股骨表观密度分布。
对于ES-FEM,其插值在离散单元内进行,而积分则是在重新构建的光滑域内完成。因此,在进行THA术后股骨的骨重建模拟时,问题域 Ω 首先被离散为一系列3节点三角形单元。然后基于离散单元的边,问题域被进一步划分为一系列光滑域 满足 且 其中 N SDs 代表光滑域数目,即背景网格中边的数目。 图1 展示了二维问题中ES-FEM光滑域的构造方式。对于内部单元的边 k ,顺次连接边的两个端点与相邻两个单元的中心点构成光滑域。对于问题域边界处的边,其所在单元中点和边 k 两个端点连线构成的三角形区域即为该边的光滑域。
在ES-FEM中,任意单元内部的位移 u ={ u , v } T 和电势 φ 可以分别由单元内各个节点的位移 d [
{
"name": "text",
"data": "i"
}
] 和电势 φ [
{
"name": "text",
"data": "i"
}
] 线性插值得到:
式中: N [
{
"name": "text",
"data": "i"
}
] 为单元节点 i 的形函数值。
引入梯度光滑技术,在光滑域内对应变和电场强度进行光滑操作可得到:
式中: 为Heaviside型分段函数,具体可写为
将式(12)代入式(10)和(11)中,压电股骨模型的光滑梯度场可以进一步表示为:
式中: S [
{
"name": "text",
"data": "k"
}
] 为光滑域内支持节点的总数; 和 分别为光滑应变矩阵和光滑电场强度矩阵
n [
{
"name": "text",
"data": "j"
}
] 代表边界 外法向向量的分量。
将式(13)和(14)代入式(2)中,可得到光滑域内的骨重建激励,具体表示为:
在假体-压电股骨模型中,利用式(13)和(14),可以得到问题域内的能量总和,即泛函
根据哈密顿变分原理,能量泛函 L 的变分形式具体表示为:
式中: F t 为系统受到的外力; Q t 为系统受到的电激励,具体可以表示为:
应变刚度矩阵、应变-电势刚度矩阵、电势刚度矩阵的计算公式如下:
1 . 4 . 1 骨重建的模拟流程 本文采用MATLAB软件开发了THA术后压电股骨重建的ES-FEM代码程序,具体步骤如下:① 利用3节点三角形单元离散假体-压电股骨模型,基于背景网格创建光滑域;② 输入迭代天数 T 、迭代步长 Δt 以及材料参数包括刚植入假体时的密度 ρ 0 、压电矩阵 e 以及介电矩阵 θ ;③ 根据骨表观密度,更新弹性模量 E 和密度函数 ρ [
{
"name": "text",
"data": "n"
}
] /ρ * ;④ 完成光滑域内刚度矩阵、外部载荷向量的求解与组装;⑤ 施加边界条件并求解线性方程组,获得各个节点的位移值和电势值,进而得到各个光滑域的应变和电场强度;⑥ 利用密度变化率函数求解股骨内部的表观密度分布 ρ [
{
"name": "text",
"data": "n"
}
] ,然后重复步骤③～⑥ ,直至迭代天数为 T 或骨密度不再发生变化(见 图2 )。
1 . 4 . 2 材料属性与重建参数 表1 列出了在模拟过程中用到的重建系数和材料参数,其中假体材料选用钛合金。
1 . 4 . 3 载荷及约束条件 进行THA术后骨重建模拟时,首先需要得到刚植入假体时股骨的初始密度,完整的股骨模型见 图3 (a)。采用2.5 mm的3节点三角形单元对问题域进行离散,节点总数为1 463,单元总数为2 699。为模拟THA操作,首先从股骨小转子上方约1.5 cm处至股骨颈外侧大转子根部进行切除,然后将选好的假体插入股骨中,建立THA术后假体-压电股骨模型[见 图3 (b)]。
该模型的离散方式与完整股骨模型保持一致,去掉股骨头后模型的节点数为1 197,单元数为2 172。为了更好地验证ES-FEM在THA术后骨重建仿真过程中的精度,将ANSYS中6节点三角形单元得到的计算结果作为参考解,单元总数为8 665,节点总数为17 748,远大于编程所用的节点数[见 图3 (c)]。因此,将其作为参考解是合理且可信的。
通过对比可知,植入假体后,接触力的受力点由股骨头变为假体,受力情况发生变化。约束条件保持不变,即股骨底部最左端完全固定且电势为零,底部其余部分仅限制垂直方向上的运动。股骨所承受载荷由3种情况所替代,即单腿站立(case 1)、髋外展(case 2)和髋内收(case 3)。每种载荷情况都包括作用在股骨头处的接触力和转子处肌肉产生的反力,3种载荷情况见 表2 [23] 。
本文首先利用ANSYS对完整压电股骨的骨重建过程进行模拟,得到完整股骨的密度分布。由 图4 (a)可知,100 d时股骨平均密度不再发生变化,重建过程趋于稳定[见 图4 (a)]。此时,股骨内部表观密度的分布如下:① 股骨的皮质骨已经形成;② 股骨近端松质骨的分布符合实际情况;③ Ward’s三角区的低密度区域已经清晰地显露出来[见 图4 (b)]。所得结果符合股骨近端密度分布的基本特征,故可以作为置换术开始时股骨内部的初始密度 [24] 。置换术后,股骨头被切除,从而得到刚植入假体时股骨的初始密度[见 图4 (c)]。
术后第100天,观察ES-FEM、FEM以及参考解得到的股骨内部总位移分布发现,FEM和ES-FEM得到的位移分布趋势与参考解相同。但是与FEM相比,尤其是在骨干上半部分和转子区域,ES-FEM的位移分布与参考解更加接近[见 图5 (a)]。
由于在THA术后的重建过程中考虑了股骨的压电效应,故可以得到股骨内部第100天的电势分布[见 图5 (b)]。结果表明,在电势的最大值区域,ES-FEM得到的结果与参考解更接近。为了定量描述ES-FEM的优势,电势最大值误差 e p 定义如下:
式中: φ num 和 φ ref 分别为数值解和参考解的电势最大值。
表3 列举了0～100 d每隔20 d由ES-FEM和FEM得到的电势最大值误差。可以看到,ES-FEM得到的误差比FEM更小,以第100天的电势为例,ES-FEM将求解误差降低了约27%,说明ES-FEM的精度更高。
根据自适应骨重建理论,可以获得ES-FEM、FEM以及参考解得到的术后100 d股骨密度分布[见 图6 (a)]。为了更加直观地展示THA后股骨内部密度分布的变化情况, 图6 (b)展示了植入假体100 d后与刚植入假体时的密度差值图。通过对比可知:① 植入假体以后,与假体接触部分的骨密度明显降低;近端皮质骨层明显变薄;② FEM得到的术后100 d的密度分布与参考解的偏差较大,主要体现在骨干上半部分;③ ES-FEM得到的股骨内部密度分布与参考解吻合程度较高,误差较小。为了定量展示两种方法在精度上的差别, 图6 (c)绘制了术后重建过程中股骨平均密度的变化。通过对比可以发现,FEM、ES-FEM与参考解密度变化的趋势相同,但是ES-FEM得到的结果与参考解之间的误差更小,大约减小30%,尤其在迭代20 d后对比更加明显。
针对在考虑骨压电效应的情况下传统FEM进行THA术后股骨重建模拟时精度较低的缺陷,本文构建了假体-压电股骨模型,利用ES-FEM对THA术后股骨近端的骨重建进行数值模拟,得到植入假体后股骨内部的位移、电势以及密度分布。通过与传统FEM得到的结果进行对比,突出ES-FEM在精度方面的优势。
(1) 植入假体后,接触力的受力位置发生改变,由原来的股骨头变为了假体头部。由于在模拟过程中考虑了股骨的压电效应,故与假体-压电股骨模型的边界条件相对应,得到了重建后的位移与电势分布(见 图5 )。因为ES-FEM中引入了梯度光滑技术,适当地软化了刚度矩阵,因此,与FEM相比,ES-FEM得到的位移分布与电势分布更精确,明显与参考解更接近。
(2) THA术后,由于股骨和假体同时承担了原本只有股骨单独承担的载荷,故压电股骨内部出现了明显的应力屏蔽现象。与刚植入假体时相比,压电股骨的皮质骨区域明显变薄,骨密度显著降低。植入假体后,外侧皮质骨层由近端至远端骨密度下降程度逐渐减小[见 图6 (b)]。
(3) 虽然FEM和ES-FEM都可以模拟THA术后压电股骨的重建过程,但是在重建过程中,ES-FEM在精度方面具有明显优势。在THA术后骨重建的过程中,FEM得到的股骨平均密度与参考解之间始终保持很大的差距,而ES-FEM则表现得十分稳定,其解与参考解始终非常接近[见 图6 (c)]。
(4) ES-FEM仅需要在FEM背景网格上进行操作,不需要额外划分新的网格,并没有增加前处理的工作量。而且由于引入了梯度光滑技术,使得ES-FEM具有与真实模型更接近的刚度。因此,在相同网格情况下,ES-FEM可以得到比FEM更精确的THA术后股骨重建仿真结果。
由于植入假体后真实股骨内部的变化需要通过CT扫描的方式得到,因此在后续研究中,为了进一步提高模型的精确度,提升在临床研究中的应用价值,可以借助CT图对此模型进行进一步的佐证。
理论分析和模拟结果表明,假体植入后,股骨的受力情况发生变化,影响了股骨近端骨重建的进程。在考虑骨压电效应的情况下,本文采用的ES-FEM能够更精确模拟全髋关节置换后股骨重建的过程,为关节置换的临床研究提供有效参考。
(1)
(2)
(3)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
(24)
(25)
(26)
(27)
微信
QQ空间
微博
QQ
