腕关节有限元骨性建模及力学分析

张浩，朱建民，马南，解品亮，钟务学，徐灵军，赵海洋

(1.江苏大学临床医学院，江苏 镇江212001;2.上海市徐汇区中心医院骨科，上海200031)

手是由骨、关节、韧带、肌肉和肌腱组成的，发挥手的任何功能都必须依赖于这些结构共同而复杂的配合来完成。只有深入地探索和研究其生物力学原理及各种疾病的发生、发展原理，手外科的整体理论深度和治疗水平才有可能得到进一步发展［1]，从而更好地应用于临床。本研究探索在个人计算机平台上依据活体CT 图像结合解剖资料，构建腕关节三维有限元模型的可行性及其应用价值。

1 材料与方法

1.1 影像设备与图像处理软件

仪器:荷兰Philips64 排螺旋CT;计算机硬件:Dell latitude D630 工作站;CPU:Intel(R)Core(TM)i3504@ 3. 07 GHz;硬盘:320G;内存:12G;显卡:NVIDIA GeForce310;相关软件:Mimics10. 01，ANSYS10.0，GEOMAGIC。

1.2 方法

1.2.1 模型的建立 CT 图像采集于1 名29 岁健康男性志愿者(图1)，无腕部疾病史及外伤史，使用64层螺旋CT 扫描机对被试者上至尺桡骨中段，下至近端指骨进行薄层扫描。扫描条件:扫描电压120 kV，层厚0.67 mm，层间隔0.67 mm，扫描图像以DICOM格式保存并刻盘，并导入三维医学图像建模软件Mimics 进行腕部骨骼三维立体模型重建。

图1 腕部64 排128 层螺旋CT 扫描图像Fig 1 Scanning images by 64 row 128 layer spiral CT

本研究主要包括骨和软骨组织三维有限元的研究。从CT 扫描图像中可以很清楚得到骨与软骨的轮廓分界，人工自由选择比较合适的单位拟定为软骨数据，其余部分拟定为骨骼数据(图2)。测量CT 扫描图像中软骨轮廓，并参考文献［2]中软骨厚度的数据，设置软骨层厚度为1.2 mm。将CT 扫描得到的图像以DICOM 3.0 格式刻盘，导入到三维医学图像建模软件Mimics 10.01，根据骨与软骨轮廓分界将骨与软骨分割，然后建立实体三维模型，并以IGES 格式导出。并再次导入到ANSYS 软件中，对模型进行网格化，并对网格化后的模型进行网格优化(图3)获得不包括软骨及其附属韧带等的完整腕部骨骼有限元模型(图4a)，并在此模型中根据软骨层厚度进行手绘各软骨面(图4b)，同时记录各骨模型的单元位置及单元编号，用于其后的材料属性的赋值，以及接触分析中接触对的定义。最终得到完整的腕关节整体三维有限元模型，包括近端尺桡骨、腕骨、掌骨、近节指骨及各软骨面。

图2 参照CT 扫描图像选取节点Fig 2 Selecting node according to the CT scan image

图3 舟骨模型的网格划分及网格优化Fig 3 Griddivision and grid optimization of scaphoid model

图4 腕部整体有限元模型Fig 4 The whole finite element model

1.2.2 模型的材料属性 生物组织属于各向异性的非线性体，而且腕部结构复杂，影响因素众多。本研究着重讨论静态载荷下的腕舟骨的受力情况(主要为舟、月骨与桡骨接触面的接触应力分布情况)，同时参考文献方法，模型中仅考虑了骨与软骨材料［3-5]，以减少计算量。

1.2.3载荷设置与边界条件 将模型中的桡骨和尺骨近端在X/Y/Z 3 个方向完全约束［3]，对近节指骨各施加轴向共100 N 的压力载荷(图5)［6]。为了更好地模拟腕管载荷传递情况，本模型将相邻腕骨、桡骨与尺骨、桡骨与腕骨、尺骨与腕骨、腕骨与掌骨各关节处的软骨单元设置为接触单元。根据解剖实验数据得到载荷作用下各骨周围可能与之接触的骨骼分布情况，设置各软骨接触对。考虑软骨之间的摩擦系数较小，且关节中的组织滑液有减小摩擦的作用，因此本文没有考虑摩擦因素的影响。

图5 有限元模型约束Fig 5 The finite element model been restricted of the wrist joint with prosthesis in ANYSYS

2 结果

建立了包括桡、尺骨远端，腕骨，掌骨，近端指骨，各关节软骨面，三角纤维软骨复合体等结构的三维有限元模型，成功模拟了掌骨受轴向压力情况，并得到各腕骨应力分布图。见图6。尤其真实有效模拟轴向压力加载下舟状骨应力分布最为集中地分布于腰部(图7)，这明确解释了临床工作中舟状骨骨折多见于腰部，并且有效验证了本模型的真实性及有效性。

图6 模型整体和局部应力云图Fig 6 the whole and local stress nephogram of model

图7 舟骨应力云图Fig 7 The stress nephoqram of scaphoid

3 讨论

有限元素法是一种有效的离散化数值计算方法，Courant 于1943年首创了有限元单元法，随着医学与工程力学的结合，有限元分析作为一种新的生物力学研究手段，逐渐被应用到医学领域，成为目前生物力学研究中先进有效的方法。Brekelmans等［7]和Rybicki 等［8]在1972年首次将有限元法应用于骨科生物力学的研究。

腕关节各腕骨形态极不规则，不能用简单的图形或公式计算出来，目前的建模方法通常是通过CT、X 线、磁共振、超声等其他医学影像资料得到研究对象的连续切面的二维图像，然后用计算机三维重建技术来得到三维模型。本研究利用64 排螺旋CT 进行薄层扫描，基于Mimics 软件首先绘制出三维模型，然后放入GEOMAGIC 软件中进行网格划分，再导入Mimics 中，然后使用布尔运算将其与相接触的骨骼相减，然后再一次导入GEOMAGIC 中进行光滑处理和重新的网格划分(此时的模型由于进行了布尔运算，与之前模型有了很大的差别，故需要重画网格)，然后再导入Mimics 中，进行整体放大1.02 倍的操作，这样就可以解决之前模型各组件之间完全没有重合的问题。通过以上操作后模型不仅有重合部分，而且重合部分很小［9-10]。

本研究首次将尺骨和桡骨远端、腕骨以及掌骨、近端指骨、各关节软骨建立成三维有限元模型，模拟轴向荷载从掌骨经腕骨传递到远端尺骨和桡骨的受力情况。本模型对掌骨受到100 N 轴向压力荷载时各腕骨接触应力的分布情况进行了三维有限元分析，并获得应力分布云图。结果表明，模型中应力集中分布区域，正是我们临床工作中骨折常见部位，如腕舟骨，应力主要分布于腰部，临床中舟骨骨折多见于腰部.这对我们进一步探讨外荷载作用下整个腕部结构的力学行为提供了一个可操作的平台和可靠的理论依据。

然而本研究在约束的施加方面发现新问题。首先，软骨最多只能在3 个节点上施加约束，不然会导致计算不收敛，具体原因暂不清楚。其次，腕骨、掌骨、指骨，初次试验于每块骨骼的下端施加约束，在计算后发现，无论施加多大的力(100 ～1000 N)，模型的最大变形和最大应力，以及应力的分布都是相同的。我们推测其原因是，模型满足本构关系，即应力应变的乘积为常数。在约束较多的情况下，模型只能有非常小的位移。本研究中位移为0.98 mm，在这样的情况下，应力值相应确定，而施加轴向压力超出此范围，则无法计算，是否可以说明在此情况下已导致骨折，还有待于进一步试验，也是我们下一步研究的一部分。

腕部功能的完整描述还应该同时包括肌腱和韧带的作用，而本模型最大的不足就是因数据提取困难，未进行肌腱及韧带建模。随着实验手段和测试技术的发展，通过建立并利用腕部结构三维有限元模型，计算和分析腕部整体结构或各组成部分如腕骨、韧带、肌肉及肌腱等组织结构在各种力荷载作用下的应力应变及位移情况，将会成为腕部生物力学的一个重要研究方向。

另外本研究在接触算法中应用MPC 算法，使得计算更具精度。bonded 接触算法使用了惩罚方法，由于接触刚度和穿透，可能会影响结果的精度，即使是线性问题，通常也需要秩代，这是因为使用了接触刚度。MPC 算法不需要操作者手工定义MPC 方程，提高了求解效率，不再需要通过多次尝试来提高和保证求解精度，求解方程不需要秩代，不仅可以约束平移自由度，而且可以约束旋转自由度，使模型后续的有限元分析更具优势。本研究建立的腕关节三维有限元模型的应用价值在于后续导入有限元软件可以进行腕部相关的生物力学分析、计算，并作为腕关节病理有限元模型的构建基础，为腕关节外科提供一定依据。

有限元的研究随着计算机技术的进步日益凸显成熟，这也将推动腕关节有限元模型的研究，本文做了这方面的初步尝试，难免存在不足，例如因数据提取困难，未进行肌腱及韧带建模，另外由于为单个样本研究，普遍的统计学意义有待于验证。

［1]顾玉东，王澍寰，侍德，等.手外科学［M].上海:上海科学技术出版社，2004:85.

［2]Streekstra GJ，Braskamp P，vander Leij C，et a1.Cartilage thickness measurement in the sub-millimeter range［C].∥Barillot C，Haynor DR，Hellier P，eds. Medical Image Computing and Computer-assisted Intervention. Berlin:Springer-Verlag，2004:950-958.

［3]Carrigan SD，Whiteside RA，Pichora DR，et al. Development of a three-dimensional finite element model for carpal load transmission in a static neutral posture［J]. Ann Biomed Eng，2003，31(6):718-725.

［4]Anderson DD，Deshpande BR，Daniel TE，et al. A threedimensional finite element model of the radiocarpal joint:distal radius fracture step-off and stress transfer［J]. Iowa Orthop J，2005，25:108-117.

［5]Dacho AK，Baumeister S，Germann G，et al. Comparison of proximal row carpectomy and midcarpal arthrodesis for the treatment of scaphoid nonunion advanced collapse(SNACwrist)and scapholunate advanced collapse(SLAC-wrist)in stage II［J]. J Plast Reconstr Aesthet Surg，2008，61(10):1210-1218.

［6]Ural A. Prediction of Colles＇fracture load in human radius using cohesive finite element modeling［J]. J Biomech，2009，42(1):22-28.

［7]Brekelmans WA，Poort HW，Slooff TJ，et al.A new method to analyse the mechanical behaviour of skeletal parts［J].Acta Orthop Scand，1972，43(5):301-317.

［8]Rybicki EF，Simonen FA，Weis EB Jr.On the mathematical analysis of stress in the human femur［J]. J Biomech，1972，5(2):203-215.

［9]方驰华，周五一，黄立伟，等. 虚拟中国人女性一号肝脏图像三维重建和虚拟手术的切割［J]. 中华外科杂志，2005，43 (11):748-752.

［10]刘文军，钟世镇. 虚拟中国人女性一号松质骨图像数据的配准与三维重建［J]. 中国临床解剖学杂志，2004，22 (4):380-383.