基于Abaqus二次开发用户自定义场的功能梯度材料人体骨骼建模.doc

基于Abaqus二次开发用户自定义场的功能梯度材料人体骨骼建模.doc 基于Abaqus二次开发用户自定义场的功能梯度材料 人体骨骼建模 摘要: 用Abaqus提供的接口进行二次开发，提出功能梯度材料人 体骨骼的建模方法.提出使用用户自定义场子程序(USDFLD)处理大变形 和大位移等复杂情况的算法以及数据传递的办法，从而避免使用复杂的用 户自定义材料子程序(UMAT);研究网格划分密度对结果和网格收敛 性的影响，以确定有限元离散人体功能梯度材料的合理性. 关键词: 功能梯度材料; 二次开发; 用户子程序; 用户自定义场; USDFLD 中图分类号: R318.01; TB34 文献标志码: B Modeling of functionally graded material human skeleton based on secondary development of user defined field subroutine in Abaqus FENG Mina， HU Muyuanb， ZHENG Bailinb， HE Pengfeib (a.Department of Physical Education; b.College of Aerospace Engineering and Mechanics， Tongji University， Shanghai 200092， China) Abstract: The modeling method of functionally graded material human bones is proposed by the secondary development on the interface in Abaqus. The algorithm that uses user defined field subroutine(USDFLD) to solve large deformation and large displacement and the data transfer method are proposed， which avoids the complex user defined material subroutine(UMAT). The effect of mesh density on the analysis results and mesh convergence is researched to ensure the the rationality using finite element to discretize functionally graded material human bones. Key words: functionally graded material; secondary development; subroutine; user defined field; USDFLD 收稿日期: 2013-06-24 修回日期: 2013-07-17 作者简介: 冯敏(1976—)，女，江西景德镇人，讲师，硕士，研究 方向为运动人体科学，(E-mail)fengmin0726@163.com 通信作者: 胡牧原(1989—)，男，浙江温州人，硕士研究生，研究 方向为固体力学，(E-mail)charlie.yaha@gmail.com 0 引 言 人类骨骼是一个非常复杂的系统，外层有较为致密的皮质骨，内层有 多孔的松质骨以及中心的骨髓腔.密质骨富含有机质，是基于哈佛氏管的 板状骨胶原纤维复合材料;松质骨由大量骨小梁构成，其间的空隙由外到 内逐渐增大. 目前，研究骨骼的力学性质主要有仿真和试验2种方法.试验一般采 用动物骨骼或者医用尸体;仿真一般从细观的骨单元结构入手，如胞元法 等，但该法计算建模的成本过高，仅局限于医学范畴内应用.对于体育生 物力学这种更关注骨骼整体力学性能、传力机制的学科，应找出适当的工 程简化方法. 1 骨骼的功能梯度属性 BARON[1]提出利用超声波探测骨骼孔隙率的梯度，并发现年龄与其密切相关，并建立骨骼的功能梯度理论模型.骨骼的力学参数直接与骨质的孔隙率相关，文献[2]中有胫骨皮质骨截面径向外侧、中部、内侧3点骨质孔隙率的统计数据，见表1和2. 表 1 各年龄段胫骨皮质骨质孔隙率和梯度 Tab.1 Cortex tibia porosity and gradient of different ages 表 2 胫骨松质骨质孔隙率 Tab.2 Cancellous tibia porosity % 将容易测定的骨质孔隙率转化为与弹性模量有关的函数，这是骨骼均匀化、建立功能梯度模型的关键.考虑到骨骼明显的各向异性，文献[3]给出刚度矩阵的几个独立项与孔隙率的关系，c11=ET(1-νTLνLT)Δ c12=ET(νTT+νTLνLT)Δ c13=ET(νLT+νTTνLT)Δ c33=EL(1-νTTνTT)Δ c44=c55=GLT c66=GTT Δ=ν2TT+2νLTνTL+2νLTνTLνTT (1)式中:L为皮质骨轴向;T为截面横向;ET=20.40-0.250×p，EL=25.98-0.274×p;泊松比νLT=0.3+0.001×p，νTT=0.28-0.0002×p，νTL=0.3;剪切模量GLT=8.51-0.116×p，GTT=7.97-0.097×p.p为骨质孔隙率，由细观有限元 模型及对比试验求得.[3] 利用骨骼具有功能梯度的属性，建立连续、均匀化的模型，以代替复杂的细观骨单元模型，减小计算规模，以适合体育生物力学研究需要. 三维功能梯度问题基本不可能求得解析解，必须借助有限元;Marc虽然可以在图形界面直接设置求解功能梯度问题，但其前处理模块相对较弱，对于骨骼这种具有复杂几何、复杂边界条件的问题，处理难度较大.Abaqus比Marc拥有更强大的前处理模块，能更加快捷地建立参数化的模型，为优化提供较大的便利.同时，Abaqus也具有值得信赖的非线性计算能力，但直到6.12版本，Abaqus/CAE模块仍未提供可以直接设置功能梯度的方法，需通过用户子程序，利用Abaqus的接口进行二次开发. 本文利用Abaqus平台编写通用且兼顾计算效率的二次开发用户子程序，对人体骨骼建立功能梯度模型，为后续复杂骨骼几何结构建模以及生物力学应用和仿真建立基础. 2 骨骼功能梯度模型的建立 2.1 几何参数和边界条件 建模对象为胫骨，忽略骨骼的不外形.几何参数基于华东地区成人骨骼标本统计值[4]，算例模型几何参数见图1，其中，人为规定骨髓腔不明显的内径为2.5 mm. 图 1 算例模型几何参数，mm Fig.1 Geometric parameters of model，mm 模型下端固支，上端施加一个大小为80 N?m沿z轴方向的弯矩，根据膝关节韧带所施加的弯矩确定.[5] 2.2 骨骼力学参数 根据表1和2的数据差值，即可得到孔隙率关于空间坐标的函数，代入式(1)即为功能梯度控制方程.骨骼被认为是横观各向同性的材料[6]，本文为控制变量个数，直观地反应骨骼模型变形与弹性模量之间的关系，优化二次开发程序调试效率，因此，暂将骨骼模型简化为各向同性.选择表1中20,29年龄段和表2中数据，以二次多项式形式拟合，E(R)=-0.1061R2+3.4165R-6.35 (2) 在有限元法中，功能梯度材料一般将控制方程离散到每个单元;而在Abaqus/CAE模块中，一般选择用户材料子程序(UMAT)或者用户自定义场子程序(USDFLD).UMAT需要用户使用FORTRAN语言或者C/C++语言编程，在UMAT中提供自定义材料本构模型的雅可比矩阵，即应力增量对应变增量的变化率，相对比较复杂;USDFLD主要应用在设置随场变化的量，恰好符合功能梯度材料的特性，弹性模量随坐标变化且程序相对UMAT简单，也与Marc处理功能梯度的思路一致.因此，本文采用USDFLD作为二次开发的接口. 2.3 用户自定义场子程序USDFLD 大多数材料属性都可通过子程序USDFLD定义.该子程序将会返回一个场变量fi，然后对不同的fi给出相应的材料属性，其中fi=f(R). 在Abaqus中定义材料时，首先在Abaqus/CAE材料编辑模块中定义Depvar以及User Defined Field选项，前者为子程序返回变量个数，后者用于在inp文件中写入USER DEFINED FIELD，激活调用子程序USDFLD;然后在Elastic面板中增加场变量个数，并定义场变量与弹性模量的关系.为方便子程序表达式的构造，根据式(2)给出的功能梯度材料模量控制 方程，结合模型的几何参数，给出材料属性参数，见表3，数值由线性插值求得. 表 3 子程序返回的场量与弹性模量的关系 Tab.3 Relation between field variable returned by subroutine and modulus 子程序一般有固定定义格式，需写入外部文件，在提交计算时利用外部编译器进行编译.子程序USDFLD格式为: #定义子程序，定义、传递变量; SUBROUTINE USDFLD (FIELD， STATEV， PNEWDT， DIRECT， T， CELENT， TIME， DTIME， CMNAME， ORNAME， NFIELD， NSTATV， NOEL， NPT， LAYER， KSPT， KSTEP， KINC， NDI， NSHR， COORD， JMAC， JMATYP， MATLAYO， LACCFLA) #导入Abaqus/Standard迭代增量步 INCLUDE ‘ABA_PARAM.INC’ #子程序执行算法部分，按照式(2)给出 RETURN END 其中，FIELD即为目标场量fi;COORD为单元积分点坐标，包含3个分量，由求解器在每个增量步开始时导入子程序.根据式(2)和表3，子程序中算法部分为 FIELD(1) = COORD(1)**2+COORD(3)**2 子程序被Abaqus/Standard求解器调用的过程见图2.值得注意是，单 元对应的fi在每个增量步的每次迭代都会被更新，因此，弹性模量也会 在每次迭代前被更新. 图 2 USDFLD子程序运行流程 Fig.2 Calculation process of subroutine USDFLD 3 结果分析 3.1 非线性大变形对结果的影响 若只以上述算法执行计算，其结果不是十分理想，截面正应力不对称，见图3(a);对比减小载荷，只发生小变形的情况，见图3(b). 图 3 大、小变形截面正应力分布对比 Fig.3 Comparison of sectional normal stress distribution under large deformation and small deformation 造成上述现象的原因是调用COORD(*)返回的是单元积分点变形后的 坐标.单元在计算过程中始终保持更新fi和弹性模量，若模型只发生小变 形，最终结果影响不大;大变形或者大位移会导致物体偏离最初定义的功 能梯度场，导致正应力结果不对称.y向正应力分布对比见图4. 图 4 y向正应力分布对比 Fig.4 Comparison of normal stress distribution in y direction 3.2 子程序修正 USDFLD子程序虽然形式比UMAT简单得多，但在大变形、大位移条件 下存在缺陷，必须修改程序，使得单元在发生位移后仍可维持单元初始的 fi和弹性模量，并在迭代更新过程中保持恒定. USDFLD仅支持通过GETVRM语句调用前一增量步储存在积分点的单元 信息，包括积分点现时坐标，而不能提供如Marc的单元初始坐标.其他涉及到位移的几何坐标都储存于节点，不能被调取，限制程序的算法. 子程序编程不能定义全局变量，子程序内的变量会在遍历单元的过程中被覆盖，调用储存函数会因为函数的生命周期不能控制而丢失储存的值，因此，一般子程序只能对前一增量步的数据进行操作. STATEV数组是用于储存单元状态数据的一个数组，也是唯一一个可以储存于单元积分点的自定义变量，不会在遍历更新单元过程中被覆盖.该变量只能被下一步迭代读取，若不重新定义，将会被清零. STATEV数组可被GETVRM语句通过‘SDV’(Solution Depend Values)调用.通过语句 GETVRM(‘SDV’， ARRAY， JARRAY， FLGRAY， JRCD， JMAC， JMATYP， MATLAYO， LACCFLA) STATEV(1)=ARRAY(1) FIELD(1)=STATEV(1) 反复继承前一步定义的STATEV数组，实现将初始状态的fi传递到后续的计算，使得每个单元的弹性模量在整个计算过程中保持恒定.该方法具有支持并行、只读写内存、无明显瓶颈的优点.程序修正后截面正应力分布对比见图5. 图 5 程序修正后截面正应力分布对比 Fig.5 Comparison of sectional normal stress distributions calculated by modified subroutine 综合上述设定功能梯度控制方程，修正大变形、大位移影响，得到简 单、通用、高效的功能梯度材料建模二次开发程序.由图5可知，功能梯度材料与普通材料截面正应力分布不同，功能梯度材料外侧模量较大，因此正应力倾向有外侧材料承担. 4 网格收敛性 截面径向功能梯度的模拟是基于单元离散的，即单元内部并没有功能梯度，本文分别从模型轴向、截面周向和截面径向等方面进行研究. 4.1 截面径向和周向单元数 截面径向的网格密度至关重要，因为功能梯度沿径向变化，径向网格足够密时才能拟合连续变化的情况.模型各截面y轴方向正应力最大值形成的曲线见图6，曲线左右两端对应模型上下两端，由于约束，曲线存在不稳定的波动，因此主要考察中部稳定段.图6(a)为径向单元密度分别为8，16，24，32和40时的几条曲线;随着单元密度增大，轴向的正应力增大，超过30后趋于收敛，见图6(b);超过30后，单元的长细比已经非常大.增加截面周向的单元个数，见图6(c)，中段截面的正应力最大值将会减小，各曲线分别为轴向划分20，30，40和50时的情况.当周向划分超过50段时，y=135 mm处截面的最大正应力趋于收敛，见图6(d). 4.2 轴向单元数量 轴向单元数量主要影响单元的长细比.在静力学问题中，对曲线两端的结果影响较大，但中段影响较小，轴向单元数量-轴向正应力曲线见图7.可以结合具体问题，兼顾计算效率，确定轴向单元数量. (a)径向不同网格密度截面内最大正应力与截面位置关系 (b)梁中点截面内最大正应力与径向网格密度关系 (c)周向不同网格密度截面内最大正应力与截面位置关系 (d)梁中点截面内最大正应力与周向网格密度关系 图 6 不同网格密度对截面正应力分布的影响 Fig.6 Curves of maximum sectional normal stress under different mesh densities 图 7 轴向不同网格密度截面内最大正应力与 截面位置关系 Fig.7 Curves of axial element number vs axial normal stress 5 结束语 研究Abaqus建立骨骼功能梯度材料模型的方法，通过用户自定义场 子程序USDFLD实现按功能梯度控制方程以及单元坐标，赋予单元不同的 弹性模量，离散地模拟功能梯度材料.开发的子程序能够修正大变形、大 位移带来的影响，通用性、效率都能保证;能充分利用Abaqus的非线性 计算能力.建模计算的结果可以进一步与现有功能梯度梁解析解进行比对 验证.[7-8] 本文用最简单的算例提出人体骨骼功能梯度的有限元法，并考虑大变 形的算法，实现Abaqus在功能梯度问题上的应用，为下一步骨骼实际外 形建模以及各向异性参数的设置带来极大的方便.目前，人体骨骼功能梯 度建模依赖全局坐标，不能使用Abaqus方便的局部坐标;控制方程需经 过几步转换，不如Marc便利，需要进一步优化.参考文献: [1] BARON C. Using the gradient of human cortical bone properties to determine age-related bone changes via ultrasonic guided waves[J]. Ultrasound Med & Biol， 2012， 38(6): 972-981. [2] BOUSSON V， MEUNIER A， BERGOT C， et al. Distribution of intracortical porosity in human midfemoral cortex by age and gender[J]. J Bone & Miner Res， 2001， 16(7): 1308-1317. [3] GRIMAL Q， RAUM K， GERISCH A， et al. A determination of the minimum sizes of representative volume elements for the prediction of cortical bone elastic properties[J]. Biomech & Modeling Mechanobiology， 2011， 10(6): 925-937. [4] 刘仪， 蒋青， 路通. 华东地区成人骨骼标本胫骨平台几何参数 的测量及其临床意义[J]. 苏州大学学报: 医学版[J]. 2007， 27(2): 214-215. LIU Yi， JIANG Qing， LU Tong. The geometric parameters of the tibia platform of East China adult bone specimens and its clinical significance[J]. Suzhou Univ J Med Sci， 2007， 27(2): 214-215. [5] 郑诚功， 黄昌弘， 魏鸿文， 等. 膝关节的生物力学性能简介 [J]. 中华骨科杂志， 2006， 26(12): 862-864. ZHENG Chenggong， HUANG Changhong， WEI Hongwen， et al. The biomechanics properties of knee joint[J]. Chin J Orthop， 2006， 26(12): 862-864. [6] 陈秉智. 计算骨力学若干问题研究[D]. 大连: 大连理工大学， 2002. [7] 杨青， 郑百林， 张锴， 等. 功能梯度双材料复合悬臂梁受集 中剪力作用的弹性力学解[J]. 兰州交通大学学报， 2013， 32(1): 20-24. YANG Qing， ZHENG Bailin， ZHANG Kai， et al. Elastic solutions for a bilayer functionally graded cantilever beam under concentrated shear loads[J]. J Lanzhou Jiaotong Univ， 2013， 32 (1): 20-24. [8] 李娟， 郑百林. 岩石流变对不同锚杆应力影响的数值分析[J]. 计算机辅助工程， 2012， 12(6): 39-43. LI Juan， ZHENG Bailin. Numerical analysis on effect of rock rheology on different anchor stress[J]. Comput Aided Eng， 2012， 12(6): 39-43. (编辑 陈锋杰)