摘要 目的 针对骨重建压电效应仿真分析时传统有限元法(finite element method,FEM)数值解精度差的问题,提出了一种基于边的光滑有限元法(edged-based smoothed FEM,ES-FEM)模型。方法 基于三角形背景网格构建光滑域,依托梯度光滑技术获得光滑应变梯度与光滑电场梯度,在光滑伽辽金弱形式的框架下构造系统离散方程。结果 采用上述模型能反映压电效应下骨重建过程中骨密度变化和电势分布,相比于FEM,ES-FEM能够在一定程度上提高骨重建仿真结果的精度。结论 提出的边光滑有限元模型能够更准确地模拟出骨重建过程,该方法对骨重建压电效应问题的准确预测为骨类疾病临床研究提供有效的理论依据。
关键词:
骨重建
边光滑有限元法
梯度光滑技术
压电效应
数值算法
骨组织作为一种活性生物器官,其随着载荷变化不断地进行生长、加强和再吸收的过程,被称为骨重建 [1-3] 。Wolff [4] 提出了关于骨的演化规律:骨骼生长会受到力学刺激的影响而改变其结构,用之则强,废用则弱,即Wolff定律。Fukada等 [5] 通过实验发现,骨骼演化过程中伴随有局部电场的电子转移,并把这一现象归结为压电效应。
对骨重建进行研究,分析压电效应对骨组织密度分布的影响,是预测骨生长的关键。目前常用的研究方法有实验法、数值仿真法。在实验方面,黄克勤等 [6] 通过对动物小腿骨施加电流研究压电效应对骨重建的影响。深田荣一等 [7] 分析骨的压电特性和压电膜对骨生长的作用机制。刘今朝等 [8] 采用生物电压放大器测量湿骨试样的电信号,并探究压电效应对骨重建的影响。富东慧 [9] 通过分析梯形载荷作用下电压的波形特点,发现切应力是影响骨压电的主要因素。然而,在研究骨重建的压电效应时,由于观察指标较多,难以处理数据结果。数值仿真通过构建计算模型进行仿真分析,从而能够直接观察骨密度的演化,弥补了实验法的缺陷。边界元法(boundary element method,BEM)是骨重建研究中具有影响力的一种数值仿真方法 [10] 。Vannessa等 [11] 运用压电边界积分方程研究骨组织的演化行为。Gonzalez等 [12] 采用孔隙弹性介质分析骨折愈合过程。Miguel 等 [13] 基于BEM模拟骨组织在剪切作用下的压电效应,并探究电刺激对骨重建的影响,然而边界元法因离散方程的满秩矩阵限制了计算骨密度的效率。在骨重建的数值仿真中,有限元方法(finite element method,FEM)是另外一种常用的模拟方法。安梅岩等 [14] 通过ANSYS有限元软件,利用内部骨重建方程分析外载荷下骨密度的演化过程。Fernandez等 [15-16] 通过建立数值模型探究骨压电对骨重建的影响。陈秉智等 [17-18] 基于应变能最优化准则建立数值模拟方法,研究骨骼内部重建的机理和规律。刘海等 [19] 应用骨量与给定载荷之间的关系,模拟计算骨重建后骨密度的变化。虽然FEM有着广泛的工程应用,但在模拟仿真骨重建时,不同的网格会导致计算结果差距较大,原因是FEM本身过于依赖网格划分,网格质量直接影响数值结果的精度和稳定性。
提高数值仿真的精度是提升骨重建问题计算结果有效性的关键。综合有限元法和无网络法的共同优势,有研究团队提出了光滑有限元方法(smoothed finite element method,S-FEM) [20-23] 。该方法通过梯度光滑操作来“软化”离散系统的刚度,从而实现对计算精度的大幅提升 [24-26] 。其中,基于边的光滑有限元法(edged-based smoothed FEM,ES-FEM)在二维压电结构的静态和频率分析中取得了较好的效果 [27] 。鉴于此,本文将ES-FEM引入骨重建的压电效应数值模拟。采用三角形单元对问题域进行离散,基于背景网格各边构建光滑域,依托梯度光滑技术获得光滑应变和光滑电场梯度,借助广义光滑Galerkin弱形式得到系统离散方程。基于上述模型,分析骨重建过程中骨密度的演化规律,进而对骨模型表面施加电荷,分析其电势的分布。结果表明,相比于FEM,ES-FEM具有更高的精度,可提供更加有效的骨重建压电效应仿真分析结果。
骨重建的压电效应表现为:在机械刺激的作用下,骨组织内部的局部电场发生电子转移,内部应力和介电位移发生改变,从而影响骨密度。应力和介电位移可表示为:
式中: σ 和 ε 分别为机械应力和机械应变; W 和 E 分别为介电位移和电场强度; c E 为弹性矩阵; e 为压电矩阵; θ ε 为介电矩阵,即
式中: e [
{
"name": "text",
"data": "ij"
}
] 为压电常数( i =1,2,3; j =1,2,…,6), β [
{
"name": "text",
"data": "ij"
}
] 为介电常数( i , j =1,2,3)。
对于本文考虑的骨重建压电效应问题,边界条件应满足 Γ = Γ t ∪ Γ u ∪ Γ D ∪ Γ φ ,其中, Γ t 为力边界; Γ D 为电位移边界; Γ u 为位移边界; Γ φ 为电势边界。
σn=t [
{
"name": "text",
"data": "N"
}
] , 在 Γ t 上
Wn=q [
{
"name": "text",
"data": "N"
}
] , 在 Γ D 上
φ=φ [
{
"name": "text",
"data": "D"
}
] , 在 Γ φ 上
对于骨重建问题,由于三维骨模型较为复杂,相关边界也并不明确,相关参考数据亦较为匮乏,而二维骨模型也能较好地反映骨重建问题。为了验证此方法在该问题上的适用性,可以近似地将三维结构简化为二维平面问题,本文采用Workbench建立二维骨模型。采用3节点三角形网格将问题域 Ω 离散为 N e 个单元,共有 N n 个节点,它们满足: 和 每个单元内的机械位移 u ( x )和电势 φ ( x )可通过场函数 d [
{
"name": "text",
"data": "i"
}
] 、 φ [
{
"name": "text",
"data": "i"
}
] 由线性插值得到:
式中: N [
{
"name": "text",
"data": "i"
}
] ( x )为节点 i 的形函数。
为进行梯度光滑操作,需在背景网格的基础上进一步构造基于边的光滑域。如 图1 所示,基于三角形单元的各边,顺次连接边的两个端点与对应三角形单元的中点创建光滑域。问题域内的光滑域个数等于单元的边数。当边在问题域的边界上时,光滑域由3条边组成;反之,当边在问题域的内部时,光滑域由4条边组成。
依托梯度光滑技术 [] ,光滑域 内光滑应变和光滑电场强度可表示为:
在本文建立的数值模型中,应变和电场强度需要通过位移和电势得到,而位移和电势在域边界上是连续的,且光滑函数在光滑域内是可微的,因此依据高斯散度定理,光滑域内的积分可变为域边界上的积分:
将式( 5 )、( 6 )分别代入( 9 )、( 10 ),在光滑域 内,可得如下形式的光滑梯度项:
在光滑域内,二维压电问题的广义光滑伽辽金弱形式如下:
将式(11)、(12)代入式(13),经整理可得如下系统离散方程:
式中: 为光滑应变刚度矩阵, 为光滑电势刚度矩阵, 为光滑应变-电势刚度矩阵。
根据人体股骨组织的模型及其边界条件,载荷施加在股骨和转子处[见 图2 (a)]。通常采用3种载荷工况表示日加载情况,每种工况都包含作用于股骨头的分布载荷和外展肌所产生的反作用力。 图2 (b)所示为日常状态下日载荷变化规律。其中,每种载荷工况的最小值为0,最大值见 表1 [13]
Fyhrie等 [28] 用数学公式定量描述骨重建,其方程可表示为:
式中: ρ 为骨密度,d ρ/ d t 为密度随时间的变化率; B 和 k 均为试验常数,总能量密度 S = S [
{
"name": "text",
"data": "u"
}
] + S φ 。其中, S u 为应变能密度, S φ 为电能密度,计算公式为:
在骨重建的压电效应过程中,各个参数随骨密度的改变而变化。骨组织弹性模量为 E ( ρ )= Mρ [
{
"name": "text",
"data": "γ"
}
] 。其中, M =3 790 Pa/(kg·m -2 ) 2 , γ =3。此外, e 31 =1.507 65×10 -9 (C/mm 2 )× ρ 3 , e 33 =1.872 09×10 -9 (C/mm 2 )× ρ 3 , e 15 =3.576 43×10 -9 (C/mm 2 )× ρ 3 , β 11 =88.54×10 -12 (F/mm)× ρ 3 , β 33 =106.248×10 -12 (F/mm)× ρ 3 , B =1(g·cm -3 ) 2 (MPa·d) -1 , k =4 mJ/g, κ =0.3 [16] 。
基于ES-FEM不依赖于网格的特性,此处采用简单三角形单元网格对骨模型进行划分,节点总数为789,单元总数为1 444。依据人体骨质代谢情况,设定骨组织最小密度 ρ min =0.01 g/cm 3 ,最大密度 ρ max =1.740 g/cm 3 ,初始密度 ρ 0 =0.8 g/cm 3 。骨重建过程均通过MATLAB程序编译并进行仿真计算。
图3 (a)和(b)分别给出采用FEM-T3和ES-FEM在同一网格节点数下计算得到的骨密度分布云图, 图3 (c)为与FEM-T3相同节点数的4节点四边形单元(FEM-Q4)的计算结果。考虑到该问题无解析解,采用节点间距极小的标准有限元解(111 49个节点)作为参考解[见 图3 (d)]。结果表明,采用FEM-T3和FEM-Q4得到的解均与参考解有较大的差别,而采用ES-FEM得出的云图与参考解吻合较好,精度较高。
为了更方便地观察ES-FEM在精度方面的优势, 表2 列出了采用不同方法计算得到的 图2 (a)中 AB 段相关点骨密度值与相对误差,如选取 X 坐标为3.321 cm的点,FEM-T3和FEM-Q4的相对误差分别可达到12.00%和9.13%,而ES-FEM的相对误差仅为0.96%。该结果更明确地验证了ES-FEM具有比FEM较高的精度。
为了定量地分析ES-FEM与FEM在精度方面的差异,选取骨密度演化差距较大的 图2 (a)中路径 AB , 图4 (a)给出了使用不同数值算法的计算结果,可以直观地发现:① 在骨模型 AB 段中间部分,FEM-T3和FEM-Q4的计算结果与参考解的误差较大,是因为此处受到载荷方向的作用较大;② ES-FEM的计算结果整体上接近参考解,误差较小。该结果表明,ES-FEM在研究二维骨重建问题时具有较高的精度,同时也验证了本文所构造模型的正确性。
由于日载荷的周期循环作用使总能量密度不断发生改变,最终骨重建达到平衡状态。 图4 (b)为日载荷作用第1天骨模型中 A 点[见 图2 (a)]总能量密度的变化规律。可以看出:① 总能量密度随日载荷的变化而变化。例如,2.4～7.2 h和14.4～19.2 h载荷不变,总能量密度也是平稳的;② 相比于FEM-T3和FEM-Q4,采用ES-FEM得到的总能量密度更接近参考解。
在骨重建的压电效应中,可以通过在骨组织表面施加电荷来改变骨密度 [13] 。为了分析施加电荷对骨重建的影响,重建期后,在不施加载荷的一段时间内,分别在转子和骨干上施加表面电荷(见 图5 )。
首先,考虑在骨模型的转子处施加2×10 -9 C·mm -2 的表面电荷,分别给出采用节点总数为789、单元总数为1 444的FEM-T3和ES-FEM计算得到的电势分布图,再采用和FEM-T3相同节点数的FEM-Q4得到电势分布。由于该问题无解析解,故采用非常密集的网格(11 149节点)通过标准FEM得到的解作为参考解[见 图6 (a)]。可以看出:① 施加表面电荷能够改变骨模型的电势,特别是在着力点附近,通过改变总能量密度来影响骨密度的变化;② FEM-T3和FEM-Q4的解在施加电荷的转子处均与参考解有较大的差别,随着与着力点距离的增加,差别逐渐减小;③ 采用ES-FEM得到的云图与参考解在整个问题域吻合较好,尤其是在着力点附近较为明显,如骨模型中 A 点电势参考解为-45.26 V,采用ES-FEM计算结果为-46.37 V,而FEM-T3和FEM-Q4得到的结果分别为-47.48、-46.93 V。由此可见,相比于传统FEM,ES-FEM精度较高。
为了验证表面电荷对骨重建的影响,同样地,在骨干处施加5×10 -8 C·mm -2 表面电荷,使用FEM-T3、ES-FEM、FEM-Q4和参考解得到电势分布[见 图6 (b)]。可以发现:① 在骨干处施加表面电荷可以改变骨干附近的电势,证实了表面电荷能够有效地影响骨重建;② FEM-T3和FEM-Q4的计算结果在骨干处与参考解有较大的偏差,随着与骨干距离的增加,FEM-T3和FEM-Q4的结果与参考解的误差出现减小的趋势,而ES-FEM在整体上能给出与参考解基本一致的计算结果。该结果进一步证实了ES-FEM能够提高数值解的精度,故ES-FEM具有很高的实用价值。
针对传统FEM在人体股骨组织重建的压电效应数值分析时精度较低的缺陷,本文构建了人体股骨重建压电效应仿真的ES-FEM模型,研究压电效应下骨重建过程及表面电荷对股骨电势的影响。理论分析和数值算例表明:
(1) ES-FEM采用3节点三角形单元即可完成对复杂骨模型的网格划分,有效地减少了网格部分的负担。
(2) 在骨重建过程中,梯度光滑技术使ES-FEM的结果不仅能够准确地模拟出骨重建过程,而且比采用相同网格FEM-T3和相同节点的FEM-Q4模型更加精确。
(3) 骨组织受到力学因素后作出内部重建等反应,受力大小和位置直接影响重建过程,在受力载荷大的区域骨密度不断增加;反之,骨密度则不断减小。该结果证明了骨组织在受载荷作用的地方不断再生;反之,则逐渐吸收。
(4) 在模型转子处施加电荷,能够明显地改变着力点附近的电势值,采用ES-FEM得到的电势值优于FEM的计算结果,通过在骨干上施加电荷验证了这一结果。同时也说明电刺激可以有效地改善骨重建,证明电疗可作为骨重建过程中的一个额外刺激。
本文基于边光滑FEM的骨重建压电效应方法对骨重建问题是一种有效的计算方法,为生物组织的宏观力学性能仿真分析提供了行之有效的途径,下一步准备将其拓展到其他更加复杂的场景中。
,
刘易 3
(1)
(2)
(5)
(6)
(
(
(
(
(
(
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
微信
QQ空间
微博
QQ
