生命力分享 http://blog.sciencenet.cn/u/bio

博文

生物力学三维有限元建模中非实体单元的应用

已有 7880 次阅读 2010-7-26 09:42 |个人分类:生物力学|系统分类:论文交流

这篇文章去年被收在在高等教育出版社出版的《数字医学的现状与未来》(工程前沿 第11卷)。近日与朋友讨论医学建模的问题,而这篇文章是我两年前写的一份在建模方面的经验小结,部分体现了我当时的一些经验和思路。由于大多数朋友没有这本书,网上又没有这本书的电子版,所以我就暂时把它贴在博客里。文章中有些图片不好粘贴,也就省略了,一方面也是为了适当保护高等教育出版社的版权。
由于这篇文章是个经验总结,所以当时写的更像博客,而不像严格意义上的论文。部分不好在期刊上或者不够资格在期刊上发表的想法可以自由挥洒。
如果有朋友对此感兴趣,可以私下向我索要全文,毕竟没有图看起来不是太爽。如果部分观点能得到认可,不胜欣慰,当然,由于不是太严谨,其中自然也有纰漏,欢迎探讨。
为保护知识产权,引用请注明出处:牛文鑫,樊瑜波*.生物力学三维有限元建模中非实体单元的应用.张绍祥,付征(主编).数字医学的现状与未来,工程前沿(11).北京:高等教育出版社.2009414-428.
生物力学三维有限元建模中非实体单元的应用
牛文鑫 樊瑜波*
北京航空航天大学生物与医学工程学院 100191 北京
 
   有限元法是生物力学研究的重要工具,有限元模型的建立在其生物力学应用中占举足轻重的地位。本文总结了作者部分建模经验,通过几个典型例子阐述了非实体单元在生物力学模型中的应用,主要包括不同单元在模拟韧带组织上的比较,关于以面代体问题的验证思考,自适应封闭壳-杆结构的进一步阐释和扩展应用,双截棍型LARS人工韧带模型和宫灯型非融合固定系统等效力学模型的建立。恰当地使用非实体单元,可以大大简化有限元模型,并可以更接近所分析问题的本质。
关键词:生物力学;有限元方法;有限元建模;单元类型
 
引 言
有限元法是生物力学研究的一个重要工具。由于生物组织在几何、材质、力学行为等方面的复杂特性,模型的建立长期以来在有限元法生物力学应用中作为主要研究内容。除非对特定问题的关注,生物组织基本上都不宜用二维模型模拟。但由于建模困难,早期的生物力学研究大量应用了二维有限元模型[1]。随着各种建模技术发展[2],目前生物力学研究已经整体跨入到三维模型时代。
有限元的基本原理就是将一个由无限个质点构成的连续体离散成有限个小单元,单元之间通过节点相连,单元之间的信息无损地通过节点传递。若干个单元的刚度矩阵集合成连续体的整体刚度矩阵,并通过数学形式表达出来。从另一个角度来讲,也可以看作是将一个无限自由度的连续体简化成有限的节点所遍历的空间,每个节点又具有有限的自由度,节点之间的信息通过单元传递,单元在传递节点信息过程中对信息进行加工和改造,模型的特性主要由单元的类型和材质决定。所以有限元分析的结果很大程度上取决于所采用的单元类型。在应用过程中,有限元法发展了很多单元形式。在生物力学三维有限元模型的构建过程中,为合理简化模型,降低模型的自由度,或更方便地模拟不易使用三维实体单元划分的结构,需要经常采用一些非实体单元构建模型。本文将结合具体实例,讨论这一技巧性和应用性较强的问题。
二、应用举例
(一) 关于韧带建模的讨论
韧带是维持正常生理活动的重要组织,也是生物力学研究所关心的重要结构之一。由于韧带的特殊几何和材料特性,以及三维重建常用的断层医学影像对韧带的显影都不太理想,所以用有限元方法模拟韧带存在很多困难[3]。在有限元建模中通常使用实体单元、膜单元或杆单元来模拟韧带。Bandak[4]应用膜单元在有限元模型中模拟足踝韧带。Cheng等[5]使用实体单元来模拟跖筋膜。其他绝大多数有限元模型均使用仅受拉的杆单元来模拟人体各种韧带[6-16]。但是对这些模拟方法的合理性和局限性,以及对分析结果的影响缺乏深入讨论。
本文使用两个相对运动的半圆形球体,一端加载,另一端分别以一定厚度的实体单元、与实体单元面积一致的膜单元和仅受拉杆单元来在120度扇面上模拟韧带结构。两个半球韧带附着扇面区域网格被细化,为考察杆单元的简化处理对模型的影响,分别建立不同数目韧带的模型共7个,韧带数目分别为126122754109,各模型中总截面积平均分配给所有单元,同时109是这个区域相对的所有节点对数,即韧带模拟的最大数目。模型如图1所示。
1                 2                3                4
图1 韧带模型比较(1.实体单元;2.膜单元;3.一个杆单元;4.109个杆单元.)
上述9个模型在相同的边界条件和加载条件下计算,模拟关节运动。除膜单元模型边界体现较大误差外,其他模型的最大应力都出现在韧带与半球体结合处,最大位移出现在加载端。但是由于实体单元具有的剪切强度,实体单元的最大位移(27.892长度单位)明显小于其他模型最大位移,最大von Mises应力(61.4应力单位)也有差异。膜单元模型的开放边界部分产生极大的计算误差,导致了不切实际的最大位移(461873长度单位)和最大von Mises应力(123应力单位),而半球体最大位移和最大von Mises应力只有207.038长度单位和55.5应力单位。如图2,把膜单元双侧边界去掉,可以看到膜单元在位移和应力两个结果上都表现出明显的类似瀑布的效应,证明可以通过杆单元来模拟该结构。从图3可见,在本文模型中韧带数目大于6时,可以认为韧带数目对位移结果影响不大,在韧带数目大于27时,可以认为韧带数目对应力结果影响不大。相对于109个最大匹配,这个密度非常小,所以由于使用杆单元简化而引入的误差在一般问题分析中可以忽略。
1                                          2
图2 膜单元计算结果的瀑布效应(1. 位移;2.应力)
图3 模拟韧带数目对结果的影响
使用杆单元模拟韧带有效利用了模型骨结构中划分的节点,未增加新的自由度,无论从操作上还是实际计算上,都节约了时间。杆单元可以通过截面和材料设置充分利用解剖和实验数据,便于参数化研究和模型设计[17]。另外,杆单元也便于输出计算韧带承受拉力,便于根据量化的结果发现规律。但是韧带在生理条件下同时具有牵张和包绕的功能,通过韧带的包绕可以有效防止骨脱位,杆单元在这方面存在不足。同时需要指出的是,如果关注韧带局部承受应力应变,或研究韧带损伤机理,使用实体单元建模应该作为首选[18]
(二) 关于以面代体问题的讨论
壳单元用来模拟厚度方向尺寸远小于另外两维尺寸,且垂直于厚度方向的应力可以忽略的结构。杨济匡与方海峰[19]建立的人体下肢三维有限元模型使用壳单元划分骨组织实体,大大减少了网格数量,但并未对该方法给出深入讨论和分析。我们之前建立的全股骨三维有限元模型生成了9340个二阶四面体单元,共14909个节点[20],而使用与之对应的二阶四边形单元来对构成体的表面进行壳单元划分,仅形成1101个单元和3139个节点,且划分过程更简便。在设置该壳单元厚度为1 mm的条件下,与实体模型一致的材料参数和边界条件、加载条件不能保证与实体模型相同的输出结果。其中输出位移与应力为实体模型相应输出结果的约4倍以上,而整体分布规律基本一致。通过对弹性模量的调节,可以控制位移结果的输出,在设置弹性模量为50~51 MPa时(实体模型弹性模量为12 MPa),可以保证所有位移输出与实体模型输出结果的误差在10%之内,基本满足生物力学计算误差要求。但是无论改变弹性模量还是泊松比,对应力输出结果的影响都不显著。同时,可以对非主要关注的硬组织设置为刚体,对其以壳单元划分,可以起到类似实体单元划分的效果[21]。总之,可以根据实验结果通过反向参数设计,简化部分不感兴趣的内容,应用该方法建立旨在分析运动学结果的特定问题的局部等效模型,并且该模型的应力结果同样具有定性分析的比较意义。
(三) 自适应封闭壳-杆结构的应用
自适应封闭壳-杆结构由牛文鑫等[6, 7, 9-11]提出并最初应用于足底软组织的三维建模。足部骨结构构成异常复杂,对足底软组织进行实体有限元划分非常困难,极易形成不良网格,影响足骨的分析精度。而且足底软组织一般不是有限元分析的重点部位,且构成复杂,材料参数很难确定,实体单元的应用将浪费大量的计算资源在不必要的位置上。自适应封闭壳-杆结构是壳单元与杆单元的统一构成体,其中壳单元是用来划分所描述实体部分边界,是该结构的骨架部分。杆单元由两部分组成,一部分为连接壳单元与其依附单元,这部分杆单元是该结构的实质部分,数量众多且不增加任何节点。另一部分杆单元为连接壳边缘的自适应约束杆,是该结构的灵魂,变开放的结构为封闭结构,避免在计算中开放的边界如同上一节中的膜单元一样被压溃和撕裂,更好地模拟有着广泛力学联系的复杂生理系统。
该方法可以扩展应用在生物力学其他模型的建立过程中,例如在关节软骨的的建模中已被成功应用[12-15]。关节软骨是附着在骨组织上,可以通过MR图像确定其相对位置,但由于关节软骨相对骨组织一般比较薄,导致有限元建模一定困难,使用自适应封闭壳-杆结构可以巧妙地解决该问题。如图4为膝关节股骨和胫骨平台软骨建模中的壳单元和封闭杆单元。
1                                         2
图4 自适应封闭壳-杆结构在关节软骨建模中的应用(1.股骨远端;2.胫骨平台)
(四) 内植物建模的两个实例
对内植医疗器械进行有效的生物力学分析,可以评价其应用效果,改进其设计,对临床实践具有重要的指导意义。但现代医学工程的迅速发展使得内植物往往具有独特的材料性质和结构设计,给生物力学建模带来一定困难。但对这些结构进行巧妙地简化,往往能起到很好的效果。
1双截棍型人工韧带模型
LARS人工韧带是法国生产的新一代产品[22],为保证重建后韧带与正常生理各种应力相适应,避免扭转剪切应力引起韧带断裂和减少滑膜炎的发病,该韧带将关节内部分设计为自由纤维,区别骨隧道内编织结构。完全使用实体单元对该韧带进行建模存在很多困难,而使用弹簧单元结合实体单元划分的方式可以顺利解决该问题。我们针对该结构,在模拟前交叉韧带重建术时提出了双截棍型有限元模型,隧道内结构以等直径圆柱模拟,并以实体单元划分,圆柱对侧对应的18对节点可以联接成18个弹簧单元。LARS韧带的材料参数来自法国ITF实验室测试结果,确定为弹性模量170 MPa所建立的LARS人工韧带和膝关节整体模型如图5所示。
应用所建立的前交叉韧带重建膝关节三维有限元模型,执行与Segawa[23, 24]相同的边界条件和加载条件,给股四头肌施加40N的载荷,股二头肌及半腱肌分别施加10 N的载荷。分别在屈膝角度为0°、30°、60°、90°及120°时,测定股骨隧道前、后、内、外壁应力。
1                                         2
图5 (1).LARS人工韧带及其双截棍型有限元模型;(2).ACL重建整体膝关节有限元模型
在整个膝关节活动范围内, 隧道壁的接触应力在前后壁明显较大,且有一定的规律性,而内外侧壁接触应力非常小,无明显规律。在屈膝0°-120°范围内,移植物对股骨隧道的前壁有持续性接触应力。膝关节完全伸直时,前壁所受的力最大,而后壁几乎不受力。此后随屈膝角度增大,前壁的接触应力逐渐减小,而后壁的接触应力逐渐增大在屈膝120°时后壁所受的力达到最大。
  10°隧道后壁的接触应力明显高于其他4个隧道,40°隧道前壁的接触应力明显高于另外4个隧道。在整个屈膝范围内, 20°隧道和25°隧道前壁的受力无差别,但在屈膝30°和90°时, 20°隧道的后壁受力明显高于25°隧道。比较25°隧道和30°隧道,就后壁受力而言,在屈膝60°和90°时, 25°隧道略高于30°隧道( 25°隧道分别为0.97 MPa1.49 MPa, 30°隧道分别为0.84 MPa和1.39 MPa) ,而在屈膝120°时, 30°隧道后壁所受的力( 2.39 MPa)则略高于25°隧道后壁所受的力(2.15 MPa) ; 而在前壁, 30°隧道所受的力(屈膝0°时为3.73 MPa,屈膝60°时为1.07 MPa,屈膝90°时为0.65 MPa)明显高于25°隧道所受的力(屈膝0°时为3.22 MPa,屈膝60°时为0.72 MPa,屈膝90°时为0.33 MPa)
这项研究借助有效的有限元建模,分析了股骨隧道角度对隧道壁接触应力的影响,认为股骨隧道与股骨纵轴在矢状面所成的角度对股骨隧道扩大确有影响,此角度在25°时能有效减轻移植物对股骨隧道的接触应力。
2宫灯型非融合固定系统力学等效模型
SCIENT’X公司发明的Semi-rigid腰椎后路固定钉棒系统代表着革命性的腰椎非融合固定新理念[25]。该系统包含一个受控微动关节,具有°的伸曲活动度和的纵向位移,作为“新铰链”起到吸震器的作用。这种独特的结构为生物力学有限元建模提出了挑战,我们抛弃了从形态上模拟该系统的观点,抓住问题的本质,建立了宫灯型力学等效模型。模型的两端设置为刚体圆柱体(r=5.73 mm),实体单元划分。中间“灯体”为连接圆柱边缘的16条等距杆三维杆单元,杆长10 mm。杆单元材料参数设置为自定义的特殊双线性弹性,在单元轴向变形为区间范围内,材料弹性模量,在轴向变形大于0.2 mm时,材料弹性模量。在极值情况下,圆柱体同一直径上两端的杆单元分别变形-0.2 mm+0.2 mm的条件下,经过计算的圆柱体直径恰好保证两圆柱体互相偏转。如此,该模型在力学上完全等效于所模拟的动态内固定系统,而且形式简洁美观,形似宫灯,具有一定的创造性和实用性。我们将这个模型添加到L3-5三维有限元模型中,并模拟轴向压缩(500 N),前屈、后伸和侧屈(10 N.m)状态,进行有限元分析,并与没有进行动态固定的情况进行比较。
6. Semi-rigid腰椎后路固定钉棒系统及其宫灯型三维有限元模型
轴向压缩状态:融合性固定L4-5后,不执行手术的情况下L3有强烈的向右后方倾斜趋势,执行动态固定后有效抑制了这种趋势。动态固定前,L3下关节突有较大的应力集中,高达19.6 MPa,动态固定后应力集中现象消失。动态固定后最大应力出现在椎体后侧边缘,为10.1 MPa,动态固定前该位置应力为14.8MPa。动态固定前后L3-4椎间盘较大应力分散分布在纤维环各区域,最大应力从固定前的后侧5.55MPa降低到固定后的前侧5.45 MPa
前屈状态:融合性固定L4-5后,L3-4不执行手术的情况下与正常模型旋转角度相近,但旋转轴略有前移,执行动态固定后,旋转角度降低为正常的约1/3,旋转轴大幅后移。动态固定前,L3下关节突有较大的应力集中,高达26.4 MPa,动态固定后应力集中现象消失。较大的应力主要分布在椎弓根,动态固定前后最高分布为8.13 MPa7.41 MPa。动态固定后两侧椎弓根应力分布趋于对称。动态固定前后L3-4椎间盘应力分布基本一致,较大应力主要分布在纤维环前部,最大应力从固定前的7.08 MPa降低到固定后的4.20 MPa
后伸状态:融合性固定L4-5后,L3-4不执行手术的情况下与正常模型旋转角度相近,但旋转轴略有上移,执行动态固定后,旋转角度降低为正常的约1/2,旋转轴后移。动态固定前,L3下关节突有较大的应力集中,高达36.0MPa,动态固定后应力集中现象消失。较大的应力主要分布在椎弓根,动态固定前后最高分布为11.3MPa8.32MPa。动态固定后两侧椎弓根应力分布趋于对称。动态固定前后L3-4椎间盘应力分布基本一致,较大应力主要分布在纤维环前部,最大应力从固定前的7.13 MPa降低到固定后的4.20 MPa
侧屈状态:融合性固定L4-5后,L3-4不执行手术的情况下比正常模型旋转角度略大,但旋转轴上移,执行动态固定后,旋转角度降低为正常的约1/3,旋转轴较固定前下移。动态固定前,L3下关节突有较大的应力集中,高达17.1 MPa,动态固定后应力集中现象消失。动态固定前L3椎体受压侧最高应力亦接近17.1 MPa,动态固定后除不考虑加载位置无意义的较高应力外,较高应力转移到了椎弓根与椎板附近区域,分布比较均匀,最高应力未超过10 MPa。动态固定前后L3-4椎间盘应力分布基本一致,较大应力集中分布在纤维环受压侧,最大应力从固定前的9.11 MPa降低到固定后的2.13 MPa
该分析应用创新性的模型设计,得出结论认为动态固定能够在保留1/3~1/2椎体自由度的条件下有效降低小关节该节段小关节的应力集中,并使得应力分布相对均匀,同时能较大程度上缓解椎间盘的压力。但是,该模型亦有一定的局限性。由于在模型的设计上仅考虑了轴向拉伸和压缩,对剪切引起的变形不能很好的模拟。因为在脊柱轴向旋转状态下,对内固定装置的作用主要为剪切力,所以本模型仅在忽略脊柱产生耦合运动的条件可以较好地分析前屈、后伸、侧屈和轴向压缩等状态下动态内固定系统对腰椎节段的影响。
 三 结论
生物力学有限元分析方法本身就是对复杂系统的简化,在建模过程中,要善于摒弃对简单现象的理解和形似的模拟,抓住问题的主要矛盾,真正从力学关系上模拟所感兴趣的内容。在生物力学有限元建模中恰当应用非实体单元,可以起到事半功倍的作用。
 
参考文献
1.      Gefen A. Stress analysis of the standing foot following surgical plantar fascia release. J Biomech. 2002; 35(5): 629-637.
2.      牛文鑫, 丁祖泉. 三种三维有限元建模方法在跟骨模型建立中的应用和比较.医用生物力学. 2007; 22(4):345-350.
3.      Weiss JA, Gardiner JC, Ellis BJ, et al. Three-dimensional finite element modeling of ligaments: Technical aspects. Med Eng Phy. 2005; 27(10): 845-861.
4.      Bandak FA, Tannous RE, Toridis T. On the development of an osseo-ligamentous finite element model of the human ankle joint. Int J Solid Stru. 2001; 38(10-13): 1681-1697.
5.      Cheng HY, Lin CL, Wang HW, et al. Finite element analysis of plantar fascia under stretch – The relative contribution of windlass mechanism and Achilles tendon force. J Biomech. 2008; 41(9): 1937-1944.
6.      牛文鑫. 足底韧带在足弓结构维持作用中的生物力学机制研究.硕士学位论文, 2007.3, 同济大学.
7.      牛文鑫, 杨云峰, 俞光荣, 等. 人体足部三维有限元模型的有效构建方法及其合理性的实验分析研究. 生物医学工程学杂志, 2009; 26(1):80-84.
8.      张宇宸,李颉, 牛文鑫, 等. 人工颈椎间盘置换对多节段下颈椎活动影响的三维有限元分析. 中华医学杂志. 2008; 88(17):1214-1216.
9.      Yang YF, Yu GR, Zhou JQ, et al. The displacements of human foot after plantar ligaments injury: A cadaveric study and finite element analysis. The 2nd International Conference on Bioinformatics and Biomedical Engineering. 2008: 2460-2463.
10. Yang YF, Yu GR, Huang SP, et al. Effect of the plantar fasciotomy on the movement of the foot arch. The 1st International Conference on Bioinformatics and Biomedical Engineering. 2007: 486-489.
11. Yang YF, Yu GR, Niu WX, et al. Effect of the plantar ligaments injury on the longitudinal arch height of the human foot. Lecture Notes in Computer Science. v 4689 LNBI, Life System Modeling and Simulation, International Conference, 2007: 111-119.
12. Shen GS, Xu YJ, Zhou HB, et al. Effect of femoral tunnel angle on tunnel enlargement in anterior cruciate ligament reconstructions. The Seventh Asian-Pacific Conference on Medical and Biological Engineering, IFMBE Proceedings 2008; 19: 103-106.
13. 沈光思, 徐又佳, 周海斌, 等. 基于MRI的全膝关节有限元模型建立与前交叉韧带重建术模拟方案. 医用生物力学. 2008; 23(5): 353-356.
14. 沈光思, 徐又佳, 周海斌, 等. 前交叉韧带重建术股骨隧道角度与隧道壁承受应力关系的研究. 临床骨科杂志. 2008; 11(5): 469-472.
15. 姜华亮, 华锦明, 许新忠, . 正常人膝关节三维有限元模型的建立. 苏州大学学报(医学版). 2008; 28(3)421-422.
16. 杨云峰,俞光荣,牛文鑫,等. 人体足主要骨韧带结构三维有限元模型的建立及分析. 中国运动医学杂志. 2007; 26(5):542-546.
17. Cheung JTM, Zhang M, An KN. Effect of plantar fascia stiffness on the biomechanical responses of the ankle-foot complex. Clin Biomech. 2004; 19(8): 839-846.
18. 姚杰, 樊瑜波, 张明, . 前交叉韧带损伤导致膝关节继发性损伤的生物力学研究. 力学学报. 2009; 已接受.
19. 杨济匡, 方海峰. 人体下肢有限元动力学分析模型的建立和验证. 湖南大学学报(自然科学版). 2005; 32(5):31-36.
20. 刘杨, 牛文鑫, 陈正龙. 健康人体正常站立位完整股骨的生物力学有限元分析. 中国组织工程研究与临床康复. 2008; 12(52): 10251-10253.
21. 杨济匡, 许伟, 万鑫铭. 研究汽车碰撞中头部动态响应的有限元模型的建立和验证. 湖南大学学报(自然科学版). 2005; 32(2):6-12.
22. Dominkus M, Sabeti M, Toma C, et al. Reconstructing the extensor apparatus with a new polyester ligament. Clin Orthop Relat Res. 2006; (453):328-334.
23. Segawa H, Omori G, Tomita S, et al. Bone tunnel enlargement after anterior cruciate ligament reconstruction using hamstring tendons. Knee Surg Sports Traumatol Arthrosc. 2001; 9 (4): 206- 210.
24. Segawa H, Koga Y, Omori G, et al. Influence of the femoral tunnel location and angle on the contact pressure in the femoral tunnel in anterior cruciate ligament reconstruction. Am J Sports Med. 2003; 31 (3): 444 - 448.
25. Philippe G. The fate of the adjacent motion segments after lumbar fusion. J Spinal Disord Tech. 2003; 16(4):338-345.
 


https://wap.sciencenet.cn/blog-2189-347379.html

上一篇:一篇中文论文审稿意见
下一篇:半蹲式跳伞着陆踝关节生物力学性别差异实验研究
收藏 IP: .*| 热度|

0

发表评论 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 20:43

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部