摘要
目的 探究适合模拟中空皮质骨结构在拉伸与压缩断裂模式下的失效准则。方法 根据前期压缩与弯曲实验数据,结合断裂仿真,对比不同失效准则模拟结果,确定模拟准确度。结果 压缩载荷下,应用主应变与不变量应变失效准则模拟断裂载荷与实验数值差异小于5%,同时断裂模式也与实验相符,说明两种失效准则可准确预测皮质骨压缩断裂;弯曲载荷下,应用等效应变与不变量应变失效准则模拟皮质骨断裂载荷与实验数值差异小于5%,且断裂模式也与实验相符,说明两种失效准则可准确预测弯曲断裂。结论 不同失效准则模拟的准确程度主要取决于应变增速是否符合当前载荷环境下骨单元变形实际情况,应变增速过快或过缓会导致结构提前或滞后断裂。断裂模拟方法能够适应多数部位与形状的皮质骨结构,可广泛应用于探究不同载荷环境下的适合失效准则,从而协助获取各部位皮质骨强度极限,为临床中提高骨断裂预测精度与掌握骨折发生条件提供数据支撑。
关键词:
失效准则
皮质骨
压缩
三点弯曲
断裂模拟精度
骨折发生率逐年递增,如何预防骨折已成为重要研究内容 。预防骨折发生前提是掌握主要承载部位断裂条件并进行预判,皮质骨作为主要承载单位,多数骨折均由皮质骨失去承载能力引起 。因此,获取皮质骨在不同载荷环境下的强度极限对于骨折预防至关重要。
由于实验样本限制,常采用断裂仿真预测皮质骨强度极限,目前模拟技术已相对成熟,各种仿真方法的模拟准确度也已得到验证 。但不同研究针对相似结构所预测强度极限往往存在差异,如有研究认为方形牛股骨皮质骨的弯曲失效载荷约为160 N,但也有仿真预测该结构的弯曲失效载荷分别为120、80 N 。由此可见,即使同一结构在相同边界条件下的预测断裂载荷也存在差异,这将影响该部位骨折条件的掌握与判断。
多数模拟均将应变作为判定单元力学状态的指标,区别在于所采取失效准则不同 。皮质骨主要呈准脆性断裂,断裂载荷主要取决于表观失效应变,而表观失效应变又与单元失效应变相关。因此,判定单元力学状态的失效准则可在一定程度上决定强度极限 。目前主要有4种失效准则被用来判定单元损伤与失效,分别为单一方向应变、主应变、等效应变,以及不变量应变,其中单一方向应变为加载方向应变,主应变为未考虑剪切应变的正应变,等效应变由主应变计算得到,不变量应变由正应变与剪应变组合而成 [11-14] 。由于计算公式不同,即使针对相同结构,4种应变在同一时刻的数值也存在差异 。由此可知,不同载荷下的皮质骨断裂模拟所适合的失效准则也应有所不同,但却鲜有研究讨论不同失效准则对于预测精度的影响。
皮质骨断裂以压缩或拉伸失效模式为主,由于表现出拉压不对称特性,故总应变中的正应变与剪应变占比随载荷环境变化 [16-17] 。因此,采用适合失效准则对于预测结构在压缩或拉伸失效模式下的断裂至关重要。基于前期已获得针对皮质骨的三点弯曲与压缩实验数据,本文采用不同失效准则进行中空皮质骨结构断裂模拟,通过对比模拟结果,确定断裂仿真精度,找到适合失效准则。本文断裂仿真方法能够适应多数部位与形状的皮质骨结构,可广泛应用于确定不同载荷作用下的适合失效准则,从而协助获取各部位皮质骨强度极限,为临床中掌握骨折发生条件提供数据支撑。
本文只进行模拟分析,但需借助前期实验数据,即应用4根股骨样本进行三点弯曲测试,加载位置为股骨干中段皮质骨区域,支撑跨距20 mm;应用4根完整股骨,在中段截取5 mm长的皮质骨结构进行轴向压缩实验[见 图1 (a)] 。
将实验所得微观影像导入Mimics 17.0生成几何模型,随后导入ABAQUS 6.14,应用C3D8单元生成三维有限元模型。本文仿真对象为中空皮质骨结构,在建模过程中将图像分割灰度阈值范围设置为800~2 000 HU,以排除股骨中部皮质骨腔内存在的微量松质骨 。
弯曲工况中,依据试验机中3根圆柱实际尺寸与相对位置,在股骨上方建立1个刚性圆柱作为压头,在下方建立两个刚性圆柱作为支撑,以模拟实验边界条件。前期实验显示,3根刚性圆柱表面均为特制,表面粗糙度较大,故弯曲过程中圆柱均未与股骨发生横向滑动 。基于此,本文将下方两根刚性圆柱与股骨接触区域设置为绑定接触关系[见 图1 (b)]。
压缩工况下,在皮质骨上方建立1个参考点,将其与皮质骨上表面耦合。轴向位移载荷作用于该参考点,同时约束参考点除轴向以外的其他自由度,并将皮质骨模型下表面的自由度完全约束(见 图2 )。
连续损伤力学模型较适合模拟准脆性断裂。虽然该方法并不新颖,但已证明适合模拟皮质骨结构断裂 。该方法最重要部分为定义损伤变量表达式,以达到刚度矩阵逐渐减小效果。分别建立单元在拉伸与压缩方向的损伤变量表达式 :
式中: D 为损伤变量; ε C 为单元压缩应变; ε T 为单元拉伸应变; ε fc 为皮质骨材料的压缩失效应变阈值; ε ft 为皮质骨材料的拉伸失效应变阈值; L c 为三维单元特征长度值,定义为单元体积开立方。
由前期实验观察得知,皮质骨结构在三点弯曲载荷下发生拉伸失效,在压缩载荷下发生压缩失效。因此,单元应变在弯曲工况下与材料拉伸失效应变阈值比较,在压缩工况下与材料压缩失效应变阈值比较。当单元应变超过失效应变时,刚度矩阵开始衰减,衰减规律依据式(1)或(2)。当损伤变量 D 接近1时,单元完全失效,无法承载,直至失效单元达到一定数量,皮质骨无法承载,发生表观失效。
皮质骨作为横向同性材料,需分别赋予横向与纵向弹性模量。依据前期实验,本批次皮质骨纵向与横向弹性模量均值分别为19.571、15.692 GPa,皮质骨纵向与横向方向泊松比为0.3 。同时,皮质骨材料的拉伸与压缩失效应变阈值分别为2.61%和4.35% 。
单一方向应变为单元在加载方向应变,可直接与材料失效应变比较,其中正值与拉伸失效应变比较,负值与压缩失效应变比较。主应变为不考虑剪应变下的正应变,将最大主应变与拉伸失效应变比较,将最小主应变与压缩失效应变比较。
式中: ε eq 为等效应变; ε 1 、 ε 2 、 ε 3 分别最大、中间、最小主应变。
分别将等效应变与拉伸和压缩失效应变阈值进行比较,由于等效应变永为正值,故将压缩失效应变阈值也变为绝对值进行比较。
不变量应变为单元在各个方向应变的综合值 :
式中:ε iv 为不变量应变; ε [
{
"name": "text",
"data": "xx"
}
] 、 ε [
{
"name": "text",
"data": "yy"
}
] 、 ε [
{
"name": "text",
"data": "zz"
}
] 、 γ [
{
"name": "text",
"data": "xy"
}
] 、 γ [
{
"name": "text",
"data": "yz"
}
] 、 γ [
{
"name": "text",
"data": "xz"
}
] 为单元应变的6个分量。
分别将不变量应变与拉伸和压缩失效应变阈值进行比较,由于不变量应变永为正值,故将压缩失效应变阈值也变为绝对值进行比较。
网格敏感性分析以本文压缩载荷下的断裂模拟为例,网格尺寸以50 μm为起始,每次减小5 μm,连续进行模拟,直至所得载荷-位移曲线达到收敛,即两种网格尺寸模型所得的断裂载荷差异小于5% 。结果显示,不同网格尺寸模型预测所得曲线形状相似,说明网格尺寸对预测精度有影响,但对于失效模式影响不大。当尺寸细化至20 μm以下时,曲线已基本达到收敛,可以看出5~15 μm网格模型预测所得断裂载荷差异均较小(见 图3 )。考虑到本文所采用断裂模拟方法不允许裂纹穿过单元内部,在保证仿真精度前提下,应尽可能细化网格。基于此,本文所有皮质骨有限元模型网格尺寸均选取为5 μm。
针对准脆性断裂,载荷-位移曲线主要取决于表观刚度与载荷下落时刻 。由不同工况下采用4种失效准则模拟皮质骨断裂所得载荷-位移曲线可见,当材料输入参数相同,仅采用不同失效准则进行断裂模拟对表观刚度未造成影响,不同失效准则所引起差异主要在于结构断裂载荷。压缩工况下,应用单向应变失效准则模拟所得断裂载荷明显大于其他模拟及实验数据,应用等效应变失效准则模拟所得断裂载荷明显小于其他模拟及实验数据,应用主应变与不变量应变失效准则模拟所得断裂载荷差异较小,并与实验数据接近,二者与实验断裂载荷差异均小于5%。三点弯曲工况下,应用单向应变与不变量应变失效准则模拟所得结果与压缩工况中相似,应用主应变失效准则模拟所得断裂载荷明显小于其他模拟及实验数据,应用等效应变失效准则模拟所得断裂载荷差异较小,而且也与实验数据也较为吻合(见 图4 )。通过结果对比可以看出,应用主应变与等效应变失效准则在两种工况中的断裂模拟结果与实验吻合情况相反。
图4 压缩与弯曲工况下4种皮质骨模型应用不同失效准则预测载荷-位移曲线对比
本文结果显示,当应用不同失效准则判定单元力学状态时,皮质骨失效单元出现位置不同,而且失效时刻也不同。压缩工况下,单向应变失效呈现较分散的失效单元布局,主要裂纹方向与轴向约呈45°夹角;等效应变失效得到大面积的损伤与失效单元,主要呈现横、纵向两种裂纹,但这两种断裂模式均与实验不符;主应变失效与不变量应变失效均呈现与加载方向成一定角度的轴向裂纹扩展,这与实验断裂模式相符(见 图5 )。
图5 压缩工况下同一皮质骨模型在不同失效准则下模拟所得断裂模式
弯曲工况下,单向应变、等效应变,以及不变量应变断裂模式与实验断裂模式符合,裂纹均从股骨下半部分产生并沿横向扩展,直至发生拉伸失效;而主应变失效与其他模拟不同,其裂纹主要出现在压头与股骨接触位置,随着上表面单元逐渐失效而发生压缩破坏,与实验结果不符(见 图6 )。
图6 弯曲工况下同一股骨皮质骨模型在不同失效准则下模拟所得断裂模式
基于载荷-位移曲线以及断裂模式对比得出,压缩载荷下,主应变与不变量应变失效能够准确模拟皮质骨断裂过程;弯曲载荷下,等效应变与不变量应变失效能够准确模拟皮质骨断裂过程。
皮质骨结构在不同载荷环境下主要呈现拉压两种失效模式,压缩失效主要由于压缩载荷作用,例如腰椎压缩性骨折;拉伸失效则可能由弯曲载荷作用形成,例如交通碰撞中的长骨弯曲断裂 。基于此,本文选取压缩与三点弯曲载荷工况进行断裂模拟。对于皮质骨这种准脆性材料,相比应力,应变判据能够更加及时准确地反映骨单元力学状态的改变 [6,11] 。因此,本文主要针对不同应变失效准则进行断裂模拟,目的是提高仿真精度,协助在临床研究中准确预测皮质骨结构断裂载荷。
针对单向应变,设定在加载方向应变与失效应变比较,这种设置只关注单元在加载方向的力学状态,误差较小,但同时也拒绝了在其他方向发生失效的可能,导致结构完全失效所需加载过大,从而引起断裂载荷增大,该解释与MacNeil等 研究相符。本文发现,失效单元较为分散[见 图5 (a)],这是由于只允许单元在加载方向失效,而无法在其他方向失效所致。应用单向应变模拟弯曲断裂与实验断裂模式相近[见 图6 (a)],这主要是由于弯曲工况下,压头进行轴向加载,而股骨上半部分为横向压缩,下半部分为横向拉伸,故难以确定单向应变失效准则的加载方向。基于此,本文设定单元在 X 、 Y 、 Z 方向的应变根据自身正负值分别与拉伸或压缩失效应变阈值比较,由此避免了在压缩工况中出现的单元唯一方向失效的问题。因此,断裂模式与实验相符,但所预测断裂载荷仍高于其他模拟及实验数据,推测原因是某一方向的单向应变在加载过程中增速缓慢,导致单元失效不及时所造成 [10,12] 。
由于主应变与等效应变均不包含剪应变,故二者在模拟过程中本质区别不大,在皮质骨断裂模拟过程中也均有一定应用 [14-15] 。压缩工况下,主应变模拟更加准确;弯曲工况下,等效应变失模拟则更加准确。具体而言,压缩工况下,等效应变增速快,单元较早达到失效条件,结构提前断裂,所预测断裂载荷小于实验数值,同时等效应变仿真所得失效单元数量较大[见 图5 (b)],也证明其增速过快,单元过早失效;弯曲工况下,则变成主应变增速过快,较早达到失效条件,导致结构提前断裂,所预测断裂载荷小于实验测试数值,同时失效单元主要出现在压头与股骨接触的上半部分,发生压缩失效,与实验拉伸失效相反[见 图6 (b)],说明股骨结构上侧单元的最小主应变增长过快,率先满足压缩失效条件,从而造成压缩破坏模式。
针对不变量应变,目前在皮质骨断裂模拟中应用较少,但在其他脆性材料中则经常应用进行断裂模拟,例如纤维复合材料,与皮质骨构造较为相似,均是在基体中布置增强结构 。因此,在纤维复合材料中适用,可能也同样适用于皮质骨断裂的判定。本文结果显示,应用不变量应变失效模拟得到的断裂载荷与实验接近,而且断裂模式也与实验吻合,说明该失效准则同时适合拉伸与压缩失效模式下的皮质骨断裂模拟,适合原因可能与该应变构成相关。不变量应变计算公式涉及各个方向的应变分量,剪应变在计算过程中占据一定比重,这与失效准则只关注正应变或单一方向应变有明显区别。因此,在应变迭代过程中能够及时反映单元在剪切层面的变形,使其更加符合骨单元的实际受力与变形情况。同时,仿真与实验结果接近还可能由于不变量应变的增速与实际情况较为符合。压缩工况下,由于加载方向与骨单元方向基本为平行状态,故变形过程中单元剪应变不大,即不变量应变此时主要由正应变确定。在不考虑剪应变前提下,不变量应变增长速度明显小于等效应变[见式(4)],故单元未过早发生失效,所预测断裂载荷符合实验测试。弯曲工况下,剪应变则比较明显,故不变量应变增速较快,计算结果与等效应变相近。与实验对比说明二者的应变变化趋势均较为符合实际情况,故预测贴合度较高。
综上所述,应用各种失效准则进行模拟的准确程度主要取决于应变增速是否符合当前载荷环境下的骨单元变形实际情况,单元应变增速过快或过缓会导致皮质骨结构提前或滞后断裂。不同载荷环境下,骨单元应变增速不同,因此,在弯曲与压缩载荷下的适合失效准则也有所区别。总体来说,采用不同失效准则对于皮质骨断裂模拟精度确实会产生影响,而何种应变增速符合骨单元实际情况,则需要进行仿真尝试,以找出最符合的失效准则。
本文的局限性如下:① 仿真分析样本量不足。限于前期实验样本限制,在两种工况下各建立4个皮质骨模型进行断裂模拟,尚未达到统计学数量,但本文目的是对比各种失效准则模拟准确性,并不是标定力学参数,故样本数量对结论不会造成大的影响;② 只建立了宏观皮质骨结构,并未针对微观皮质骨组织进行断裂模拟,故无法预测骨单元在不同载荷环境下的强度极限。下一步拟进行皮质骨不同层次的建模,以从微观层次讨论结构力学性能的变化特点。
作者贡献声明: 范若寻负责论文设计、撰写与修订工作;王溢童负责有限元仿真分析;苗前程负责有限元模型建立;胡辰负责仿真数据整理与分析;贾政斌负责断裂仿真程序代码调试。