髋臼周围截骨术中球面软骨层的动步态接触数值模拟

李 飞,徐澍民,陈荷娟

(南京理工大学机械工程学院,南京 210094)

股骨头的前部和外侧部分通常被发育不良的髋臼不完全覆盖[1]。发育不良的髋关节接触覆盖不足会导致日常活动中接触压力增加[2-3],因此不仅会导致髋关节不稳定和脱位,而且也会影响骨关节炎的发展和其他一些并发症[4-5]。为了纠正这种异常情况,全髋关节置换术(total hip replacement,THA)是缓解晚期骨关节炎患者疼痛和恢复髋关节功能的一种选择。尽管THA的寿命已在技术上得到了改善,但对于某些年轻患者而言仍然不确定[6-7]。另一方面,髋臼周围截骨术(periacetabular osteotomy,PAO)是预防骨关节炎发展并保留髋关节原始功能的有效治疗方法,该方法更适用于年轻患者[8-9]。对于患者具有足够的剩余软骨层,PAO是治疗晚期骨关节炎有前途的替代方法[9]。

可通过数值模拟深入了解异常解剖结构与软骨退化的机械原因之间的关系[10-11]。软骨层的精确表示对于获得接触力学的真实有限元模拟至关重要[12-13]。文献[14]进行了有限元模型设计分析,比较形状记忆合金智能膝盖垫片以增强膝盖功能的性能。另外,通过在数值模拟分析中假设软骨层具有恒定的厚度,引入了一些简化的软骨层模型[15]。根据PAO的最佳手术规划,研究了软骨层厚度分布和压缩特性的影响[16]。文献[17]开发了另一个基于患者软骨层特定几何形状的三维有限元模型,以模拟基于形态学的前后软骨层的接触压力变化。大多数软骨层模型已被视为同心球窝关节,并表示为一系列具有恒定长度和刚度的离散弹簧[18-19]。由人体测量学研究表明,股骨头和髋臼像一个球形,因此髋关节的球形软骨层是合适的[20]。然而,这些研究集中在正常髋关节的单一步态上,而不是研究不同步态(动步态)时的髋臼周围截骨术软骨层接触压力分布。

基于以上考虑,现阐明在髋臼周围截骨术中在动步态载荷下数值模拟球面软骨层接触压力。通过侧向旋转截骨段来模拟髋臼周围截骨术。比较动态步态下的接触压力分布、峰值接触压力和接触面积,评估球面软骨层模型在髋臼周围截骨术中的性能,以期为髋臼周围截骨术提供生物力学上的参考。

1 原理方法

1.1 髋臼周围截骨术数值模拟几何模型

为了保证真实性和个体性,髋关节几何模型的建立依据电子计算机断层扫描(computed tomography,CT)拍摄的髋关节病人的影像数据资料。将CT图像转换成DICOM数据格式,然后再转换成BMP文件格式,分辨率为512×512像素,切片间隔为1.25 mm,整个人体髋关节由117张二维CT图像组成。基于二维切片之间的线性插值生成3D体积,并渲染STL格式存储的三角化表面。应用MESHLAB将表面表示分为股骨和髋臼两部分。图1所示为建立的球面软骨层的几何模型。

图1 球面软骨层髋关节几何模型Fig.1 Geometric model of spherical cartilage layer of hip joint

根据矢状面的对称性,仅对骨盆的一半进行建模以降低计算成本。如图2(a)所示,为了模拟PAO,使用半径为40 mm的球体切开髋臼周围的截骨段。根据骨科医生的建议决定球的位置,该位置超过软骨顶部间隙上方20 mm的位置。将截骨段逆时针侧向旋转多个角度以模拟截骨手术,其中φ是截骨段的侧向旋转角度,如图2(b)所示。

图2 髋臼周围截骨术切球侧向旋转示意图Fig.2 A lateral rotation diagram of the bone-cutting ballon PAO

1.2 球面软骨层数值模拟建模方法

软骨形态与其接触力学之间关系的准确量化对于PAO的规划至关重要。当前对PAO的了解是基于临床观察和数值计算模型,这有助于改进PAO的手术规划。如图3所示,进行几何建模后,由美国ANSYS的(ICEM CFD)在几何模型上进行有限元网格划分,其中四节点四面体单元(用于骨骼)和棱镜层单元(用于软骨层)捕获其内部应力变化。图4所示为软骨层网格截面,为了确保精度并提高计算效率,将加密的网格应用于靠近软骨层的区域,而在远离软骨层看作为刚体的次要骨头的区域生成较粗的网格。

图3 球面软骨层髋关节有限元模型Fig.3 Finite element model of spherical cartilage

图4 球面软骨层网格截面Fig.4 Mesh cross-section of spherical cartilage

1.3 材料特性、载荷及动步态边界条件

基于人体结构组织的复杂性,为简化模型降低模拟成本,不区分股骨与髋骨盆的皮质骨和海绵骨,股骨头和髋骨盆为各向同性弹性材料,并统一定义弹性模量为17 GPa,泊松比为0.3[21-22]。基于日常生活的步态,所以弹性模量定为短时负荷状态时的11.25 MPa,泊松比定为0.45[23-24]。

如图5所示,体内产生的髋关节数据应用于股骨头中心，在髋骨盆模型的上部限制x、y、z3个方向的自由度,这一约束固定了上部的位移。根据对称性,在对称面限制x方向的自由度，如图5中圈出部分,模型中对耻骨韧带部分也做了相应的简化,也同样限制x方向的位移。

图5 边界约束条件Fig.5 Boundary bound condition

在正常的步行周期中,平均合力及其分量的变化如图6(a)所示,数据来自Bergmann等[25-26]的研究。由于软骨层之间的摩擦系数较低(0.01～0.02),因此忽略了摩擦力[27]。如图6(b)所示,髋关节合力的大小、方向以及大腿骨角θ都随步态而变化。选取其中具有代表性的6个步态进行有数值模拟接触分析。将这6个步态数据列于表1中,并加以说明描述。大腿骨和地面垂线之间的夹角和 6个选定的分析步态如下。

表1 6个步态时的髋关节力Table 1 Hip strength in six gaits

图6 正常步态合力及其各个分量变化Fig.6 Normal gait force and its various component change stuitals

步态1最初状态,开始迈脚的时刻,大腿骨与垂线之间角度为24°。

步态2最大载荷发生在脚后跟蹬地时刻,大腿骨与垂线间角度为15°。

步态3大腿骨与地面垂直状态,大腿骨与垂线间角度为0°。

步态4步态载荷的第2个峰值状态,大腿骨与垂线间角度为-10°。

步态5第2个大腿骨与地面垂直状态,大腿骨与垂线间角度为0°。

步态6最小步态载荷状态,大腿骨与垂线间角度为20°。

1.4 软骨层数值模拟接触策略

软骨层接触分析中采取单面接触,具体为股骨头软骨层作为从接触体、髋臼软骨层作为主接触体来进行接触计算,这样既提高了收敛性，也有利于应力分布的连续性。并且由于初始状态的不确定性使得计算收敛变得困难,所以在计算过程中适当放宽两个接触条件参数:接触容许穿透量和节点分离反力。在参数调整过程中首先必须保证收敛,再适当地牺牲一定精度,但又不能使精度太低而背离实际情况太多。

在利用MSC.MARC对髋骨节模型计算之前,还必须对模型做一些检查和优化。优化是重新给模型的节点和单元编号,使得数值分析的总体刚度矩阵的元素分布合理,计算时占用尽可能少的CPU时间、内存和磁盘空间。求解器可以利用刚度矩阵的对称性、带状分布、稀疏等特性,提高计算速度,缩短计算时间。对于MSC.MARC,采用Gbbs-Pool-Stk的带宽优化算法来对节点编号进行优化。

2 数值模拟结果

在发育不良的髋关节上用Bergmann等[25-26]髋关节合力来数值模拟计算软骨层接触峰值压力和接触面积,研究PAO术后的接触力学环境变化情况。将选定的正常行走期间6个动步态动作的接触峰值压力和接触面积记录在表2中。根据接触峰值压力,找出在动步态中具有最高风险的步态为第4步态,即脚尖着地时的步态4。

表2 6个步态时的接触峰值压力和接触面积Table 2 Contact peak pressure and contact area in six gaits

2.1 球面软骨层的接触压力分布和接触面积

图7～图12分别为前述6个动态步态动作的球面软骨层间的接触压力分布和接触面积数值模拟结果。

图7 步态1球面软骨层间接触压力分布和接触面积Fig.7 Gait 1 contact pressure distribution and contact area between spherical cartilage layers

图8 步态2球面软骨层间接触压力分布和接触面积Fig.8 Gait 2 contact pressure distribution and contact area between spherical cartilage layers

图9 步态3球面软骨层间接触压力分布和接触面积Fig.9 Gait 3 contact pressure distribution and contact area between spherical cartilage layers

图10 步态4球面软骨层间接触压力分布和接触面积Fig.10 Gait 4 contact pressure distribution and contact area between spherical cartilage layers

图11 步态5球面软骨层间接触压力分布和接触面积Fig.11 Gait 5 contact pressure distribution and contact area between spherical cartilage layers

图12 步态6球面软骨层间接触压力分布和接触面积Fig.12 Gait 6 contact pressure distribution and contact area between spherical cartilage layers

2.2 较高风险步态下髋臼周围截骨术中软骨层接触压力

在较高风险的步态4,通过侧向旋转髋臼截骨段10°～30°模拟PAO手术。图13所示为球面软骨层的PAO模型。

图13 球面软骨层PAO模型Fig.13 Spherical cartilage pao model

根据MSC.MARC的数值模拟接触计算,将较高风险的步态4时的接触峰值压力和接触面积记录在表3中。

表3 步态4不同转角的接触峰值压力和接触面积Table 3 Gait 4 contact pressure and contact area at different angles

比较图14中PAO术前和术后的球面软骨层接触压力变化,PAO术前位于软骨层前外侧的应力比较集中,在PAO术后(截骨段旋转10°～30°) 得到了明显释放,接触压力分布趋向均匀。比较表2和表3,接触峰值压力从术前的3.73 MPa降到PAO旋转截骨段20°时的1.82 MPa,下降约51%。

图14 球面软骨层模型接触压力分布Fig.14 Contact pressure distribution of spherical cartilage model

2.3 较高风险步态下髋臼周围截骨术中软骨层接触面积

图15所示为球面软骨层PAO软骨层模型在较高风险的步态4下,旋转髋臼截骨段10°～30°模拟PAO术后的接触面积。其中,左上角注明准确计算的接触面积。从PAO术前和术后的接触面积来看,接触面积在PAO术后显著增加约40%,从而说明髋臼对股骨头的覆盖增大。

图15 球面软骨层模型接触面积Fig.15 Contact area of spherical cartilage model

3 讨论

为了更好地评价球面软骨层数值模拟模型,反映软骨层的接触压力分布,分析了一例样本患者的髋关节,并研究了基于提出的PAO软骨层模型在动步态下的接触力学行为。通过对模型的分析可以发现,步态4是这一患者的最危险步态。原因是髋臼的前外侧边承受了更多的载荷,而且在步态4时的载荷方向也朝向前外侧,这共同导致了前外侧的应力集中现象。通过旋转截骨段模拟PAO可以看到,旋转截骨段来矫正髋臼的覆盖取得了增加接触面积的效果,从而相应地改善了接触压力分布,最终增加了股骨头的稳定性。提供了PAO术前规划的参考依据,以确定合理的PAO手术规化。

对髋关节来说接触压力分布和接触压力的大小同样重要,而目前为止对于髋关节接触压力分布的人体实验数据还相对较少,文献[28]在这方面取得了比较重要的实验数据,将传感器置于先天覆盖不全的人工髋关节的10个部位,进行接触压力测量。步行动作中接触压力的分布如图16所示,最大接触压力位于软骨层3号位置,其次是5号位置,相似于本研究髋臼周围截骨手术前球面软骨层接触压力分布(图7～图12),这就意味着实验结果和接触数值模拟结果得出了同样的结论。也就是,髋臼覆盖不全的髋骨节软骨层的前外侧部位承受较大的应力,也表明了软骨层应力的集中区域受到大腿骨头合力与髋臼软骨层形状的重要影响。

图16 髋关节接触压力分布[28]Fig.16 Hip joint contact pressure distribution[28]

4 结论

旨在评估髋关节的软骨建模方法以模拟髋臼周围截骨术。基于球面软骨层的建模方法,进行了PAO术前和术后动步态数值模拟接触分析。数值模拟了接触压力分布、峰值接触压力和接触面积。并对球面软骨层的建模方法进行验证。

由于接触的软骨表面一致,因此接触面积较大,在球面软骨建模方法中获得了相对连续的接触压力分布和较低的峰值接触压力。对于定量评估,球面软骨层建模方法适合于对具有近似球面股骨头和髋臼的患者进行长期预测模拟；而球面软骨层模型不可能反映骨的真实形状,所以它不适合模拟退化的平状股骨头。