人体生物力学研究中的有限元仿真计算

djh2003

虽然是几年前的文章了，但是基本技术思路没有变化，仍有参考价值
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
原文见：丁峻宏, 丁海. 生物力学研究中的数值仿真计算. 高性能计算发展与应用, 2010, (2): 31-35.
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

人体生物力学研究中的有限元仿真计算
        丁峻宏 上海超级计算中心     上海 201203 jhding@ssc.net.cn
        丁  海 蚌埠医学院附属医院   安徽 233000

摘  要: 本文首先对有限元方法及其在生物力学中的应用进行了简要回顾，指出了目前有限元法在生物力学研究中的分类，给出了具有普遍特性的生物力学有限元仿真实现流程。在大量文献基础上，整理分类出了有限元方法是如何应用在人体结构组织器官研究诸多方面并分别给出了具体引证实例。最后介绍了我中心在开展生物力学数值仿真合作研究方面所具有的资源优势，以及已经取得的一些计算成果。
关键词：ANSYS；接触；装配；结构分析



概  述  
有限元法是适应电子计算机的使用而发展起来的一种比较新颖和有效的数值计算方法。该方法起源于20世纪50年代航空工程中飞机结构的矩阵分析方法，1956年波音飞机工程师Turner，Clough等人在首次将有限元法用于飞机机翼的结构分析通过对复杂连续系统的离散化，将问题转而利用计算机进行处理，有限元这一名称在1960年被正式提出。经过数十年的发展，目前有限元法几乎可以用来解决所有的连续介质和场问题，包括静力、动力以及与时间有关的变化问题，等等。可以说，有限元方法已经成为工程设计中不可或缺的一种重要方法，在大型结构强度分析、变形分析、疲劳分析、传热分析、电磁场分析、流体流动分析等方面扮演越来越重要的角色。
人类经过长期的劳动进化后,人体骨骼已形成了一个几乎完美的力学结构。然而在对人体力学结构进行力学研究时,力学实验几乎无法直接进行,这时用有限元数值模拟力学实验的方法恰成为一种有效手段。经过长期发展，有限元法在生命科学研究中的应用,也取得了很多成绩，尤其在人体生物力学研究中,更显示出其极大优势。
20世纪60年代,在心血管系统的力学问题研究中,有限元法得到了初步应用。从70年代起,开始应用于骨科生物力学研究,最初应用于脊柱。80年代后,应用范围逐步扩展到颅面骨、颌骨、股骨、牙齿、关节、颈椎、腰椎及其附属结构等生物力学研究中。
    本文综合了最近几年来国内外诸多从事生物力学数值仿真研究的多篇文献，着重描述了有限元数值仿真方法是如何以及应用到哪些人体研究方面的，最后介绍了我中心所拥有的一些可用于人体生物力学研究的工程计算资源，并给出应用实例。
1  有限元法在生物力学研究中的应用分类
在生物力学领域，根据研究工作的性质和对象，有限单元法的应用研究工作大致可分为有限元理论算法的改进应用、医疗器械的力学性能优化测评、以及生物力学实验仿真等几个方面。
首先，在开展生物力学有限元工作的过程中,许多学者作了有限元应用理论的基础工作。针对各种问题需要，新的分析模型、算法和设计被提出来。如针对人体组织特性，利用等参变换方法、非线性以及复合材料分析方法和Galerkin 加权残值法等不同数值仿真的调整和改进。
其次，在牙科、骨科以及关节等人体部位的整形治疗中，常常会用到辅助性的医疗器械。力学性能的好坏往往决定了其临床应用价值的大小，因此对其力学性能的评价很重要。除了传统实验手段外，用有限元法进行的力学模拟具有实验时间短、费用少、可模拟复杂条件、力学性能测试全面及可重复性好等优点。另外，不少有限元软件具有优化设计模块，可以根据需要根据力学性能指标进行结构优化，改进原有设计。
最后也是应用最广泛的，就是结合计算结构力学和计算流体力学现有的数值软件，针对人体不同部位进行力学仿真。利用现有有限元软件日趋强大和完善的建模功能及其接口工具，可以拟实建立三维人体骨骼、肌肉、血管等器官组织，并模拟其生物力学材料特性。可以模拟各种类型的边界条件和载荷约束（几何约束、固定载荷、冲击载荷、温度特性等)，进行结构静力学、动力学、疲劳、流体力学等各种类型的仿真模拟，从而获得在不同虚拟实验条件下任意部位的变形、应力/ 应变分布、内部能量变化、流动特性以及极限破坏预测等特性。
2  生物力学有限元仿真路线
人体内部非常复杂，不同的组织器官的力学特性差别很大，因此具体建模计算时，也各具特点。图1给出了一个典型的实现生物力学有限元仿真的实现方法及顺序。

图1  典型生物力学仿真实现方法及路线
从一般从事有限元数值仿真人士的角度来看，上述方法和路线的大部分在本质上和处理机械设计等工业领域仿真问题的处理手段基本一致，不同的就是医学计算模型的前处理建模环节更具特点。目前传统设计领域的模型构建会借助PRO/E、UG、CATIA等著名的三维CAD设计软件直接创建模型，而人体组织构型多是曲面轮廓且没有统一规律，因此在生物力学仿真过程中，往往首先执行并读取CT断层扫描数据，进行三维轮廓造型时，需要借助一些软件比如比利时Materialise公司的交互式的医学影像控制系统MIMICS，以及美国Able Software公司开发的医用三维图形建模系统3D-DOCTOR。对模型建模精度要求更高时，可将骨骼、血管等三维重建模型转入GEOMAGIC(Raindrop 公司，美国)等软件中对数据中的图像进行修补并转化为NURB 曲面， 最后再转入不同的有限元分析软件中进行后续处理步骤。
3  有限元在生物力学中的典型应用
鉴于人体本身复杂性和可认知程度，有限元软件的发展水平和应用面不同，以及专业和跨学科领域研究人员的数量和能力的差异，从国内外来看，人体骨骼生物力学是目前该领域应用最为广泛和成熟的一个方向，它也可以归为计算结构力学理论和方法的一类应用。
3.1 　颅面骨、颌骨有限元研究
头颅及颞下关节是有限元在生物力学中研究的重点之一。通过建立包括鼻上颌复合体、下颌骨及牙齿在内的颅颌面硬组织形态的三维有限元模型，用来研究颅颌面硬组织在正中矢状面上的形态特征及因生长、正畸矫治、正颌手术引起的形态变化。也可以对下颌骨坚强内固定下的应力遮挡作用进行有限元力学分析，对下颌骨体部、角部骨折在骨愈合的不同时期、不同的咬合形式及不同的内固定方法时的应力遮挡率进行了计算分析。图2摘自文献，描述的是通过建立的下颚-颌骨-牙齿咬合分析模型以及LS-DYNA有限元软件，可以研究在牙齿不同咬合形式下颞颌关节中的应力分布情况。

图2  下颚-颌骨-牙齿咬合分析模型
3.2  脊柱有限元研究
脊柱生物力学仿真是有限元法在生物力学中研究的最早、分析的最多、也是临床上应用最广泛的领域。现今研究脊柱的工作使模型不仅能逼真地模拟椎骨、椎间盘,还能将周围的韧带、肌肉直接或间接地加入模型，使模拟更加真实与完善。这些工作不仅要求建立逼真的脊椎模型，而且要求测试椎间盘、周围韧带、肌肉的各种力学性能。
有限元在颈椎生物力学中的研究对象又可以细分为椎体、椎间盘、后部结构以及肌肉韧带等软组织。此外，内固定器械的生物力学研究，也有助于选择正确手术方法，以取得最佳矫形和固定效果。图3所描述的头颅－颈椎模型，曾用于进行人体颈部位置冲击损伤方面的研究。
有限元分析在脊柱腰椎段的应用，也覆盖了生理负载及外来负载下腰椎各部分应力分布，手术内固定及人工假体，脊柱内固定对邻近脊柱结构影响，骨质疏松椎体压缩性骨折，以及肌肉和韧带在有限元模型中的应用等多个方面，有力促进了脊柱动力学(载荷下的脊柱运动)、运动学(椎体间运动)和脊椎及椎间盘内部的应力应变等各种研究。图4描述的文献中的ABAQUS腰椎模型，曾用于对不同程度骨质疏松病症所进行的研究。

图3  颅骨－颈椎冲击分析模型

图4  腰椎骨质疏松分析模型
3.3 关节有限元分析
人体关节尤其是大腿骨两端的髋关节以及膝关节，一直以来也是病症多发部位，利用限元数值模拟是一种对人体关节力学结构进行研究的有效方法。上肢的肘关节、腕关节的研究常常与骨折以及其他骨骼创伤性疾病的应力分析联系在一起。而在髋关节方面，有限元分析最为广泛地应用于全髋关节置换的研究，分析全髋关节置换术前术后髋关节应力的分布情况，而且还可对骨水泥残余应力的细致分析和假体设计进行研究。如图5所示，对于膝关节分析来说，建立一个完整的三维ANSYS计算模型，不仅可以了解各部位的应力分布和工作原理，还有助于人工膝关节置换的合理设计。

图5  膝关节接触分析模型
3.4 足部有限元分析
当FEM应用于足部生物力学研究时，复杂的骨胳几何结构、边界条件和材料的不均匀性等问题便找到了可能的解决途径，人们也尝试对足部骨骼、软组织等结构内部的应力传递机理进行力学解释。如图6所示，此类计算分析模型不仅可以分析Lisfranc损伤和Midfoot融合等足部疾病，研究例如Hansen氏病和糖尿病人发生的足骨变形等病理现象的力学成因，还可以就鞋垫的舒适性、高跟鞋的致病性等日常问题进行分析。

图6  足部骨骼-韧带-软组织模型
3.5 人体软组织和流质的模拟
除了骨骼以外，人体软组织的研究也在不断深入。人体软组织研究主要针对人体运动系统皮肤以下骨骼之外的肌肉、韧带、筋膜、肌腱、滑膜、脂肪、关节囊等组织以及周围神经、血管。一般在以骨骼为主要研究对象的同时，如果需要考虑软组织的存在和影响，处于使用方便灵活的角度出发，通常都尽量在结构分析软件里面寻求模拟方案。如图7所示，这时候寻找合适的非线性材料本构模型来模拟对应的脑内多种不同物质特性就非常重要。图8则利用有限元方法来研究眼球周围直肌和眼窝脂肪对眼球悬挂所起的作用。

图7  颅脑内各组成部分有限元模型

图8  眼球-直肌-脂肪有限元分析模型
有时血液等流质因素也不可忽视，如开展心血管等疾病如主动脉瘤的研究时，血液动力学参数如壁面切应力、压力和血流速度等与动脉瘤的生长及破裂有着重要联系。这时就要借助CFX、FLUENT等流体力学分析软件来进行胸主动脉瘤的血流动力学分析，获得血液流场的流线、速度矢量、血管壁面压力等关注对象的分布和变化情况，见图9。

图9  血液流动流线云图
4  中心资源和研究基础
4.1 可为生物力学研究提供的资源
诚如相关学者所指出的，与国外相比，目前国内的生物力学有限元工作存在不少问题，比如研究课题多是追随跟踪国外研究进展，建立的有限元模型粗糙简单，材料本构和参数不够准确，许多研究仍停留在结构静力分析状态，对各类结构动力学、流体力学、热力学及力热耦合场研究较少，尤其是不同的研究工作者从不同的实际问题出发，重复建模浪费大量人力、物力，没有国内的公共资源。对此，上海超级计算中心目前所拥有的硬件、软件和人力资源正可以为国内从事生物力学研究的学者和科研项目提供很好的计算平台和技术帮助。
上海超级计算中心目前拥有的魔方（“曙光5000A”超级计算机）采用145台曙光TC2600-CB85刀片服务器以及82台8路A950服务器构建计算节点，提供总额达180Tflops的主频峰值计算能力，按照应用领域不同分为A/B/C三个区。
其中A区专门用来服务来自工程领域的用户，它由82个机架式服务器构成，主要部署各类工程计算商业软件，每个服务器含8个AMD Barcelona 1.9GHz四核处理器，128GB共享内存，节点之间采用InfiniBand光纤网络全线速互联，总计算能力达到20TB，配置磁盘存储100TB。
目前在A区中配备了许多通用商业有限元软件，如针对结构力学分析的ANSYS、ABAQUS、NASTRAN、MARC、LS-DYNA等软件，以及针对流体力学分析CFX、FLUENT、STAR-CD等业内公认的当前最好最新的数值仿真软件。对于有些大模型计算的前后处理工作，我部门还配置了多台工作站，双路四核HP工作站，最大内存达16GB，处理上百乃至上千万网格模型都不成问题。
更重要的的是中心的一支高学历研究团队，在结构力学和流体力学方面具有扎实的理论基础和丰富的专业有限元课题项目经验，不仅可以服务和保障用户以正确高效地使用上述计算软硬件资源，还可以联合进行课题研发或提供专业CAE项目咨询。
4.2 已开展的一些合作性研究
我中心从事工程计算的技术团队，利用单位的软硬件资源曾和来自上海医学领域的不少用户合作，完成了一些生物力学方面的研究工作。
以髋关节研究为例，图10所示的ABAQUS计算模型，主要通过有限元分析研究不同程度的股骨头缺损对全髋关节表面置换后生物力学的影响。图11所示的计算模型，可以将计算区域缩小后，建立更加精细的有限元全接触非线性计算模型，从而研究手术后不同强度载荷对置换关节和髋臼的性能的影响。而图12建立的现为超精细髋臼模型的单元和节点数都超过了300万，利用LS-DYNA软件，可以观察在人工股骨头植入髋臼过程中，髋臼的接触部位和磨损状态的变化情况。

图10  全髋关节表面置换中的股骨头分析模型

图11  人工股骨头和髋臼受力的接触分析

图12  髋臼精细模型的动力学仿真
5  结论
随着计算机软硬件的不断发展，以及跨学科研究人才的不断涌现，有限元分析方法在人体生物力学研究领域中作为新兴研究手段的作用将会更加明显和重要，极具发展潜力。在拓宽高性能计算走进医学领域的同时，我中心愿与国内医学类专业人士一道，将人体生物力学研究工作进一步推向深入。
参考文献
[1]        李人宪. 有限元法基础. 国防工业出版社, 2002.
[2]        张美超, 钟世镇. 国内生物力学中有限元的应用研究进展.  解剖科学进展, 2003, 9(1): 53-56.
[3]        崔红新, 程方荣, 王健智. 有限元法及其在生物力学中的应用. 中医正骨, 2005, 17(1): 53-55.
[4]        胡辉莹, 钟世镇, 聂晨阳. 人体骨骼生物力学中有限元分析的研究进展. 广东医学, 2007, 28(9): 1532-1534.
[5]        关哲, 马迅. 有限元法在颈椎生物力学研究中的应用及进展. 临床医药实践杂志, 2008, 17(8): 641-643.
[6]        赵文志, 苏晋, 陈秉智等. 有限元分析在腰椎生物力学研究中的应用. 中国组织工程研究与临床康复, 2009, 13(30): 5927-5930.
[7]        顾耀东, 李建设. 有限元分析在步态生物力学研究中的应用进展. 北京体育大学学报, 2007, 30(8): 1080-1082.
[8]        付文宇, 乔爱科. 基于CT图像的胸主动脉瘤模型数值模拟. 科技导报, 2009, 27(20): 27-31.
[9]        王涛, 王丙强, 陈小龙等. 魔方(曙光5000A)超级计算机应用环境. 高性能计算发展应用, 2009, (3): 27-41.
[10]        Gaivile Pileicikiene, Algimantas Surna, Rimantas Barauskas. Finite element analysis of stresses in the maxillary and mandibular dental arches and TMJ articular discs during clenching into maximum intercuspation, anterior and unilateral posterior occlusion. Stomatologija, Baltic Dental and Maxillofacial Journal, 2007, 9(4): 121-128.
[11]        Ee Chon Teo, Qing Hang Zhang, Tian Xia Qiu. Finite Element Analysis of Head-neck Kinematics Under Rear-end Impact Conditions. ICBPE 2006.
[12]        MA.Tyndyk, V.Barron1, PE.McHugh. Effect of Osteoporosis on the Biomechanics of the Thoracolumbar Spine: Finite Element Study. European Cells and Materials, 2005, 10 (Suppl. 3):73.
[13]        Raúl Lesso-Arroyo, Julio Cesar Sánchez Jiménez, Ramón Rodríguez Castro. Biomechanical Behavior of the Knee Joint Using ANSYS. 2004 International ANSYS Conference Proceedings.
[14]        J. Cheung, G. de Vrles, B. Nigg. Biomechanical Effects of Midfoot Fusion - A Finite Element Study. Journal of Biomechanics, 2007, 40, Issue null: S326-S326.
[15]        Tamer El Sayed, Alejandro Mota, Fernando Fraternali. Biomechanics of traumatic brain injury. Comput. Methods Appl. Mech. Engrg. 197 (2008) 4692–4701.
[16]        Sander Schutte, Sven P.W. van den Bedemb, Fred van Keulen, et al. A Wnite-element analysis model of orbital biomechanics. Vision Research, 2006, (46): 1724–1731.