基于两种不同建模方式的胫骨有限元对比分析

唐勇涛，魏思奇，程刚，赵建华，蔡繁繁，欧阳汉斌



基于两种不同建模方式的胫骨有限元对比分析

唐勇涛，魏思奇，程刚，赵建华，蔡繁繁，欧阳汉斌

【摘要】目的探讨两种不同建模方式对胫骨有限元模型生物力学特性的影响。方法采集1例正常志愿者胫骨CT断层扫描数据，对其胫骨模型进行三维重建，分别重建出单纯胫骨外轮廓模型以及带松质骨和髓腔结构的腔体模型，采用相应的灰度公式赋值法和均一赋值法完成模型材料参数赋值，同时对两组设定相同的载荷和边界条件，模拟站立位状态下胫骨上段受到的轴向压缩和扭转载荷，对比两组模型的胫骨应力和位移分布云图、轴向压缩刚度以及扭转刚度。结果两组不同方式构建的胫骨模型总体位移分布基本一致；灰度赋值组的胫骨有限元模型在两种载荷作用下产生的应力峰值和位移幅值均小于均一赋值组，其轴向压缩刚度和扭转刚度均高于均一赋值组；均一赋值组模型在胫骨中下段位置的应力分布相较于灰度赋值组更为集中，相对符合临床实际情况。结论两种建模方式均可满足胫骨有限元分析需求，均一赋值法建模具有更高的准确性，模型应用的灵活性更高。

【关键词】胫骨；骨重建；模型，解剖学；有限元分析；材料赋值；应力，物理；生物力学；计算机模拟

作者单位：518106广东，深圳市光明新区人民医院骨一科（唐勇涛，魏思奇，程刚，赵建华，蔡繁繁）；510515广州，南方医科大学人体解剖学教研室（欧阳汉斌）

E-mail：1241195475@qq.com

近年来，以生物力学为基础的有限元分析方法在人体骨骼结构力学分析中发挥着巨大作用。相对于传统生物力学实验手段，有限元仿真方法能够更为有效地预测骨骼在生理性和异常冲击受力状态下的应力分布以及多个方向的应变情况[1]，预测准确性高[2-3]。随着螺旋CT的广泛应用和医学影像三维重建技术的日趋成熟，有限元分析方法也在不断完善，学者们提出各种类型的建模方法，以获得更为准确的预测模型[4-5]。胫骨作为人体下肢骨的一个重要承载结构，有关其尸体标本结构特性的生物力学研究已有文献报道，但运用有限元分析方法对胫骨进行生物力学研究，各家报道不一，特别是在胫骨模型材料属性赋值方法的选用上，尚未得出一致的结论[6-8]。本文就胫骨有限元分析建模的材料属性赋值问题展开初步研究，旨在阐明不同建模方法对模型最终分析结果的影响，为胫骨生物力学研究提供可靠的参考依据。

1　材料与方法

1.1胫骨模型数据采集与三维重建

选取一名正常志愿者右侧胫骨作为数据采集来源，排除受检者受试肢体相关骨骼结构先天异常、损伤及退变等病理变化。采用64排螺旋CT机（西门子公司，德国）沿胫骨干长轴方向进行扫描，范围覆盖膝关节及踝关节，扫描层距1 mm，共获得连续横断面1 000张512×512像素的CT图像，并将获得的CT图像以通用DICOM标准格式刻录至光盘导出影像工作站。进一步将胫骨影像数据导入Mimics 14.0软件（Materialise公司，比利时）进行胫骨皮质外轮廓、胫骨髓腔、胫骨远近端的图像分割处理，将分割的各个结构蒙板提取处理并计算出胫骨皮质骨、松质骨以及髓腔的三维模型，导出后以STL文件形式保存（图1）。

1.2模型表面处理与胫骨腔体结构的实体建模

图1　胫骨髓腔内部复合结构的二维图像分割及三维重建

首先将上述模型（皮质骨、近端松质骨、远端松质骨、髓腔）的STL文件分别导入Geomagic 12软件（Geomagic公司，美国）中进行模型表面平滑处理和精确曲面构建；然后在三维建模软件UG NX 8.5（西门子公司，德国）中进行含髓腔结构的胫骨建模相关布尔操作，最后得到完整实体模型（图2），并以STP文件格式导出皮质骨和远近端松质骨部件。

1.3灰度赋值法与均一赋值法胫骨有限元模型建立

将胫骨实体模型导入有限元分析软件Abaqus 6.10（达索公司，法国）中划分四面体网格，并以inp文件形式导入Mimics赋予材料属性，Mimics将根据材料属性公式[9]，在原位二维影像的相应位置自动将灰度转换为局部有限元网格单元的弹性模量和泊松比（图3）；对均一赋值法模型采用均一材料属性在Abaqus软件中完成赋值，将两种骨组织结构均视为各向同性材料，其中皮质骨弹性模量为16 700 MPa，松质骨弹性模量为155 MPa，泊松比均取值0.3[10]。

1.4轴向压缩及扭转载荷与边界条件设定

图2　胫骨三维模型的表面处理与实体建模

图3　灰度赋值法构建的胫骨模型材料属性分布图

根据文献报道，模拟人体站立位状态下胫骨平台关节面受到的轴向压缩载荷为700 N、扭转载荷为20 Nm，载荷作用区域范围位于胫骨平台上方以及胫骨的干骺端上部[11]。另外，将胫骨下端踝关节关节面上方2.5 cm以下范围设置为6个自由度完全约束的边界条件，然后建立2个分析步以先后施加轴向压缩和扭转载荷。完成设置后，利用Abaqus软件的网格划分功能对模型进行四面体网格划分，采用C3D4实体单元，获得两组网格化模型。其中灰度赋值模型单元数为217 006个，均一赋值模型单元数为549 535个。最后将模型提交到有限元求解器中进行运算。

2　结果

2.1基于灰度赋值法和均一赋值法的胫骨模型应力、位移分布云图

最终提取出两组胫骨模型的应力分布云图见图4。两种受力状态下灰度赋值组胫骨的应力幅值要显著低于均一赋值组：灰度赋值组轴向压缩的最大应力峰值为70.57 MPa，扭矩作用下最大应力峰值为17.21 MPa；而均一赋值组轴向压缩的最大应力峰值为108.40 MPa，扭矩作用下最大应力峰值为27.42 MPa。此外，灰度赋值组胫骨干整体应力分布较为均匀，高应力区主要位于骨干中段以下；均一赋值组胫骨干的整体应力分布主要集中在胫骨中下段交界处，呈现出明显的应力集中现象。

两组胫骨模型的位移分布如图5所示。两种受力状态下灰度赋值组胫骨发生的位移幅值略低于均一赋值组：灰度赋值组轴向压缩最大位移为13.97 mm，扭矩作用下最大位移为3.25 mm；而均一赋值组轴向压缩最大位移为15.75 mm，扭矩作用下最大位移为3.87 mm。此外，两组模型无论是在轴向压缩还是在扭转作用下，胫骨总体产生的位移分布均较一致，总体趋势符合生理状态下胫骨受力后发生变形的规律。

2.2两组胫骨模型轴向压缩刚度和扭转刚度对比结果

通过计算载荷-位移曲线的线性斜率，得到两组模型的整体刚度，灰度赋值组胫骨模型的轴向压缩刚度和扭转刚度（50.11 N/mm，13.36 Nm/°）均高于均一赋值组（44.44 N/mm，11.23 Nm/°）。

图4　两组胫骨模型在轴向压缩与扭转作用下的应力分布云图（单位：MPa）

3　讨论

一直以来，有限元分析方法被广泛应用于航天、土木和机械结构等工程领域，被公认为是一种有效、准确、低成本的力学结构分析方法[12-13]。在骨科生物力学领域，有限元方法可以对人体骨骼、肌肉和韧带等组织进行仿真分析，除了为骨科植入物的结构力学分析和改良设计提供可靠依据外，在人体结构生物力学基础研究方面同样具有广阔的应用前景，有效弥补了传统生物力学试验手段的不足。近年来，国内外学者围绕人体骨骼的有限元分析方法提出了不少新的理论和方法，其中材料参数模拟方法一直是领域内争论的焦点[14-15]。作为下肢承重结构的重要组成部分，胫骨结构的生物力学分析研究有了长足的发展，但对胫骨有限元建模方法方面的相关研究较少，特别是胫骨建模过程中采用何种材质赋值方法更符合临床实际情况，目前尚无定论[16]。本研究分别通过灰度赋值法和均一赋值法建立两种内部几何结构及自身材料属性参数均截然不同的胫骨有限元模型，模拟轴向压缩载荷和扭转载荷作用下胫骨结构的生物力学响应，并据此探索出一种更符合临床实际情况的建模方案。

图5　两组胫骨模型在轴向压缩与扭转作用下的位移分布云图（单位：mm）

灰度赋值建模手段比较简单，通过参考既往文献中的“灰度-弹性模量”经验公式即可完成材料属性的指派，整个建模自动化程度相对较高，理论上比较符合个体化有限元分析的需求[7]。然而，这一途径也存在局限性：首先，CT图像质量决定了材料赋值的分布，其中包括数据来源因素，如尸体标本与活体骨骼的成像差异以及成像设备的参数差异等；其次，灰度赋值是根据模型单元体素和二维影像灰度像素的位置关系进行联系的，临床应用中的一些骨骼虚拟操作如截骨、内固定等，将使骨骼相对原有的几何形态和位置发生改变，从而影响灰度赋值的准确性。

对于均一赋值法建模，既往文献也报道了支持的观点和研究结果，Vulovic等[17]利用均一赋值法分别对正常志愿者股骨三维有限元模型中的骨皮质、骨松质、ward三角区松质骨等结构进行模拟生理载荷下的有限元分析，结果显示，股骨结构的生物力学响应接近于尸体标本测试数据。本研究中采用的均一赋值法，在建模阶段和材料参数指派方面与灰度赋值法存在着不小的差异。均一赋值法要求分别建立皮质骨、松质骨和髓腔的实体模型，通过布尔运算获得仅由皮质骨和松质骨两部分组成的具有腔体结构的实体模型，在几何结构上更接近真实胫骨，并可分别对松质骨和皮质骨赋予相对恒定的弹性模量和泊松比。然而，均一赋值法建模过程同样存在一些不足，包括几何建模的准确性存在人为差异、建模路线复杂以及类似股骨近端骨小梁特征性的形态分布在均一化处理中被忽略等等。

本研究结果显示，在应力分布方面，灰度赋值组胫骨模型的应力分布相对均一赋值组更为均匀，主要原因之一可能是因为髓腔内脂肪组织同样存在灰度。尽管其赋予的弹性模量比较小，但低弹性模量单元填充整个髓腔后与骨性结构构成一个整体并共同承载传递的应力效应，这一机制可能在一定程度上造成了胫骨干应力均匀分布的现象。在应力、位移和刚度方面，灰度赋值法模型也表现出更高的抗压和抗扭转力学性能，这一结果也与大量髓腔单元填充产生的应力分载效应有关。而均一赋值组则呈现出胫骨中下段部位发生的应力集中现象，因为此处为胫骨直径由宽变窄的过渡部位，工程材料截面发生变化的位置恰恰是应力集中的位置，所以这一应力分布现象更加符合实际情况。位移分布方面，两组模型之间呈现的分布云图具有较高的一致性，通过分析位移分布方向和胫骨受力后变形的形态，均可观察到这样的一致性，表明两种建模方法获得的模型变形趋势，无论是在轴向压缩载荷还是扭转载荷，均符合正常人体胫骨受力变形规律，这也验证了本研究建模方法的有效性。而对于两组模型胫骨结构刚度的分析，如同其应力、位移幅值差异一样，也存在类似的差异，灰度赋值组胫骨模型的轴向压缩刚度和扭转刚度均大于均一赋值组。

从应用角度来看，均一赋值方法有限元模型应用的灵活性相对灰度赋值法建模更高，可用于复杂矫形外科手术的仿真分析及应力应变计算。李长有等[18]报道Chiari骨盆截骨术中外展肌力负荷对股骨头软骨生物力学影响的有限元分析，其采用的均一赋值骨盆有限元建模方法可实现对各种骨盆截骨术式的模拟分析，便于设置不同外展肌力等载荷条件，以获得髋关节软骨接触压应力的分布结果。戴金龙等[19]应用有限元方法比较分析传统股骨转子下封闭楔形外翻截骨术与改良“Z”形转子间外翻截骨术的术后力学差异，其采用的也是均一赋值的材料属性设定方法，在分析过程中可灵活设置股骨模型的截骨角度，以便于组间模型的对比。相比之下，灰度赋值法的有限元建模要求赋值前的网格模型空间坐标位置需和原有CT图像所覆盖的范围保持一致[20]，一旦骨骼结构的网格模型空间位置发生变化（如旋转截骨、内外翻截骨等），通过经验公式赋予材料属性将无法实现或形成错误的属性赋值，采用均一赋值法可有效解决这一问题。

综上，尽管两种胫骨有限元模型在建模和赋值操作上有明显差异，但最终通过仿真分析得到的结果在生物力学变化趋势上具有一致性，从而验证了有限元模型的有效性。而作为一种生物力学分析预测手段，均一赋值法构建的胫骨有限元模型较之灰度赋值模型能更准确地反映人体胫骨的生物力学特性，在影像数据要求和后续模型的修改应用方面也更具优势。

参考文献

[1]何旅洋,张志强,郑百林.外骨骼框架生物力学设计及强度分析[J].医用生物力学, 2014, 29(6): 504-510.

[2]Keyak JH, Rossi SA, Jones KA, et al. Prediction of femoral fracture load using automated finite element modeling [J]. J Biomech, 1998, 31(2): 125-133.

[3]Wang X, Sanyal A, Cawthon PM, et al. Prediction of new clinical vertebral fractures in elderly men using finite element analysis of CT scans [J]. J Bone Miner Res, 2012, 27(4): 808-816.

[4]汪方,石杜芳,王秋根,等.基于冰冻切片的人体骨盆有限元模型的建立与初步验证[J].复旦学报(医学版), 2010, 37(4): 384-390.

[5]Dall'Ara E, Luisier B, Schmidt R, et al. A nonlinear QCT-based finite element model validation study for the human femur tested in two configurations in vitro [J]. Bone, 2013, 52(1): 27-38.

[6]徐志才,胡广洪,黄振宇,等.胫骨模型对膝关节有限元分析结果影响的探讨[J].中国数字医学, 2014, 9(4): 69-72.

[7]张国栋,廖维靖,陶圣祥,等.股骨颈有限元分析的赋材料属性方法探讨及有效性验证[J].中国组织工程研究与临床康复, 2009, 13(52): 10263-10268.

[8]Chen G, Wu FY, Liu ZC, et al. Comparisons of node-based and element-based approaches of assigning bone material properties onto subject-specific finite element models [J]. Med Eng Phys, 2015, 37(8): 808-812.

[9]Lengsfeld M, Schmitt J, Alter P, et al. Comparison of geometry-based and CT voxel-based finite element modeling and experimental validation [J]. Med Eng Phys, 1998, 20(7): 515-522.

[10]Completo A, Duarte R, Fonseca F, et al. Biomechanical evaluation of different reconstructive techniques of proximaltibia in revision total knee arthroplasty: an in-vitro and finite element analysis [J]. Clin Biomech, 2013, 28(3): 291-298.

[11]Gray HA, Taddei F, Zavatsky AB, et al. Experimental validation of a finite element model of a human cadaveric tibia [J]. J Biomech Eng, 2008, 130(3): 031016.

[12]杨锋平,罗金恒,张华,等.金属延性断裂准则精度的评价[J].塑性工程学报, 2011, 18(2): 103-106.

[13]谈梅兰,武国玉,梁福祥.基于ABAQUS的连杆疲劳分析[J].中国机械工程, 2013, 24(5): 634-638.

[14]苏佳灿,张春才,陈学强,等.骨盆及髋臼三维有限元模型材料属性设定及其生物力学意义[J].中国临床康复, 2005, 9 (2): 71-73.

[15]Taddei F, Pancanti A, Viceconti M. An improved method for the automatic mapping of computed tomography numbers onto finite element models [J]. Med Eng Phys, 2004, 26(1): 61-69.

[16]Austman RL, Milner JS, Holdsworth DW, et al. The effect of the density-modulus relationship selected to apply material properties in a finite element model of long bone [J]. J Biomech, 2008, 41(15): 3171-3176.

[17]Vulovic S, Korunovic N, Trajanovic M, et al. Finite element analysis of CT based femur model using finite element program PAK [J]. J Serb Soc Computat Mech, 2011, 5(2): 160-166.

[18]李长有,王禹祥,范广宇. Chiari骨盆截骨术中外展肌张力对股骨头软骨生物力学影响的三维有限元分析[J].医用生物力学, 2007, 22(1): 79-83.

[19]戴金龙,李炫升,梁兰丽,等.改良式股骨转子间外翻截骨术的生物力学研究:三维有限元素分析[J].中华创伤骨科杂志, 2007, 9(6): 510-514.

[20]Synek A, Chevalier Y, Baumbach SF, et al. The influence of bone density and anisotropy in finite element models of distal radius fracture osteosynthesis: evaluations and comparison to experiments [J]. J Biomech, 2015, 48(15): 4116-4123.

（本文编辑：白朝晖）

基础研究

A comparison between two finite element analysis methods of human tibia via different modeling patterns

TANG Yongtao*, WEI Siqi, CHENG Gang, ZHAO Jianhua, CAI Fanfan, OUYANG Hanbin. *Orthopaedic Department I, Shenzhen Guangming New District People's Hospital, Shenzhen, Guangdong 518106, China

【Abstract】Objective To investigate the effects of two modeling approaches on biomechanical response of finite element models of human tibia. Methods A group of data of a volunteer's CT scan was obtained for 3D reconstruction of tibia, and the solid tibial model with cavity structure was then reconstructed by isolating medullary canal and cancellous bone. Material properties in 2 models were assigned via homogeneous and inhomogeneous modulus, identical load case and boundary conditions were applied to simulate the axial compression and torsion loading on the tibia during the mid-stance phase of gait circle simultaneously. Distribution of stress and displacement, stiffness of axial compression and torsion were compared between two models. Results Regarding to the displacement, two models presented similar distribution on the bone surface. Under the two load cases, inhomogeneous model resulted in lower maximum stress and maximum displacementcomparing to homogeneous model. Similarly, inhomogeneous group presented higher axial compression stiffness and torsional stiffness than homogeneous group under the same load case. A physiological stress concentration was observed on the distal part of tibia diaphysis in the homogeneous model, which relatively matched the clinical situation. Conclusions Both of the two modeling patterns could be available for finite element analysis of tibia models. However, homogeneous modeling was associated with higher accuracy for prediction of biomechanical response and flexibility in application compared to inhomogeneous modeling.

【Key words】Tibia; Bone remodeling; Models, anatomic; Finite element analysis; Material assignment; Stress, mechanical; Biomechanics; Computer simulation

中图分类号：R356.1

文献标识码：A

文章编号：1674-666X(2016)01-026-07

DOI：10.3969/j.issn.1674-666X.2016.01.005

收稿日期：（2015-12-03；修回日期：2016-01-12）