不同数值方法模拟皮质骨结构弯曲断裂的准确性

doi:10.12307/2024.035

中国组织工程研究 ›› 2024, Vol. 28 ›› Issue (12): 1895-1900.doi: 10.12307/2024.035

• 数字化骨科 digital orthopedics • 上一篇下一篇

不同数值方法模拟皮质骨结构弯曲断裂的准确性

范若寻1，2，刘杰2，贾政斌3

1扬州工业职业技术学院，交通工程学院，江苏省扬州市 225000；2吉林化工学院，航空工程学院，吉林省吉林市 132022；3北京航空航天大学，生物与医学工程学院，北京市 100191

收稿日期:2023-01-30 接受日期:2023-04-07 出版日期:2024-04-28 发布日期:2023-08-22
通讯作者: 贾政斌，北京航空航天大学，生物与医学工程学院，北京市 100191
作者简介:范若寻，男，1988年生，黑龙江省大庆市人，汉族，2016吉林大学毕业，博士，副教授，主要从事汽车碰撞安全与人体损伤生物力学相关研究。
基金资助:
吉林省自然科学基金项目(YDZJ202301ZYTS250)，项目负责人：范若寻

Fracture accuracy on cortical bone structure under bending load using different numerical methods

Fan Ruoxun1, 2, Liu Jie2, Jia Zhengbin3

1School of Traffic Engineering, Yangzhou Polytechnic Institute, Yangzhou 225000, Jiangsu Province, China; 2School of Aerospace Engineering, Jilin Institute of Chemical Technology, Jilin 132022, Jilin Province, China; 3School of Biomedical Engineering, Beihang University, Beijing 100191, China

Received:2023-01-30 Accepted:2023-04-07 Online:2024-04-28 Published:2023-08-22
Contact: Jia Zhengbin, School of Biomedical Engineering, Beihang University, Beijing 100191, China
About author:Fan Ruoxun, MD, Associate professor, School of Traffic Engineering, Yangzhou Polytechnic Institute, Yangzhou 225000, Jiangsu Province, China; School of Aerospace Engineering, Jilin Institute of Chemical Technology, Jilin 132022, Jilin Province, China
Supported by:
Natural Science Foundation of Jilin Province, No. YDZJ202301ZYTS250 (to FRX)

摘要/Abstract

摘要：

文题释义：

数值方法：针对某种结构进行断裂模拟的仿真分析方法。在此文中主要指能够针对大鼠股骨皮质骨结构进行断裂模拟的仿真方法，具体包括基于单元瞬时失效理论、连续损伤力学理论以及扩展有限元方法下的3种断裂模拟方法。
弯曲断裂：指大鼠股骨结构在三点弯曲载荷作用下发生断裂失效的一种破坏模式，在该失效模式下，股骨结构发生准脆性断裂，裂纹起始于股骨中部下侧区域，在弯曲载荷作用下从下向上逐渐扩展，直至裂纹完全贯穿整个中部皮质骨结构，导致结构发生整体断裂失效。

背景：当前针对皮质骨结构的断裂模拟主要有基于单元瞬时失效、连续损伤力学以及扩展有限元理论3种方法。虽然断裂模拟文献众多，但却鲜有研究比较3种断裂方法的模拟精度。

目的：探究3种数值方法模拟皮质骨结构弯曲断裂的准确性与适用性。
方法：首先针对大鼠股骨样本进行三点弯曲实验，然后依据样本微观扫描影像建立大鼠股骨有限元模型，分别采用3种数值方法模拟皮质骨结构在三点弯曲载荷下的断裂，通过将预测断裂载荷及断裂模式与实验结果进行对比，以判定各方法模拟皮质骨结构弯曲断裂的准确性。

结果与结论：①由于单元失效策略不同，3种数值方法针对同一股骨有限元模型预测所得断裂进程存在明显差异；②结果显示基于连续损伤力学理论的断裂模拟与实验结果一致性更好；③通过与实验对比可以找到更适合模拟皮质骨结构弯曲断裂的数值方法，通过观测不同数值方法预测断裂参数差异并解析差异存在的原因，能够为各数值方法确定断裂模拟的适用范围，以提高仿真精度。

https://orcid.org/0000-0002-1604-2489 (范若寻)

中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程

关键词: 皮质骨, 断裂模拟, 单元瞬时失效, 连续损伤力学, 扩展有限元

Abstract: BACKGROUND: Current fracture simulation for cortical bone structure is mainly based on three numerical methods: the element instantaneous failure, continuum damage mechanics, and extended finite element methods. Although many studies focus on cortical bone fracture simulation, few have compared the differences in prediction accuracy using the three numerical methods.
OBJECTIVE: To probe the accuracy and applicability of the three numerical methods in simulating cortical bone fracture under bending load.
METHODS: The rat femur samples were primarily used to perform the three-point bending experiment. The rat femoral finite element models were established based on the micro-CT images of the femur samples and the three numerical methods were used to conduct the fracture simulations under three-point bending loads. The predicted fracture loads and fracture patterns were compared with the experimental data to determine the accuracy of various numerical methods in simulating cortical bone fracture.
RESULTS AND CONCLUSION: (1) The discernible differences in the failure processes could be observed in the same finite element model under the three numerical simulations due to different element failure strategies. (2) The simulation results showed that the fracture simulation using the continuum damage mechanics method was in better agreement with the experimental results. (3) The numerical method that was suitable for simulating cortical bone fracture under bending load could be determined by comparing it with experimental results. The variations in the fracture parameters were observed, and the reason for the differences in the predicted results using different numerical methods was also discussed, which aided in determining the range of applicability of structural fracture simulation for each numerical method and then improving the simulation accuracy.

Key words: cortical bone, fracture simulation, element instantaneous failure, continuum damage mechanics, extended finite element

中图分类号:

范若寻, 刘杰, 贾政斌. 不同数值方法模拟皮质骨结构弯曲断裂的准确性[J]. 中国组织工程研究, 2024, 28(12): 1895-1900.

Fan Ruoxun, Liu Jie, Jia Zhengbin. Fracture accuracy on cortical bone structure under bending load using different numerical methods[J]. Chinese Journal of Tissue Engineering Research, 2024, 28(12): 1895-1900.

图/表（结果） 2

应用3种仿真策略针对同一大鼠股骨有限元模型进行断裂模拟，均模拟出结构的表观断裂。图3显示了4个股骨样本与有限元模型所得载荷-位移曲线对比。可以看出虽然针对同一模型进行断裂模拟，但不同数值方法预测所得载荷-位移曲线形状却存在一定差异。EIF与XFEM仿真曲线在弹性阶段基本重合，而CDM仿真曲线在弹性阶段的表观刚度呈逐渐下降趋势，并且随弯曲载荷增大，其表观刚度与前两个仿真模拟所得数值差距表现出逐渐增大的趋势。弹性阶段过后，CDM仿真曲线在表现出较短暂的屈服特性后随即发生整体断裂失效，EIF与XFEM仿真曲线则基本未表现出屈服现象，而是从弹性阶段直接进入断裂阶段。

通过对比断裂参数可以看出，EIF断裂模拟下，股骨有限元模型发生完全断裂时刻最早，但由于在弹性阶段表观刚度未见下降，其预测所得断裂载荷仍高于CDM仿真所得数值，而且该方法模拟下的断裂时间也明显早于实验结果，但断裂载荷与实验结果差异不大。XFEM模拟得到的断裂时间与实验较为吻合，但断裂载荷却明显大于实验结果。CDM仿真预测结果无论断裂时刻还是断裂载荷均与实验结果一致性更好。
图4展示了应用3种不同数值方法模拟下的断裂模式对比以及与实验断裂模式的对比。图4A显示了应用EIF模拟所得断裂模式，其主裂纹出现在股骨中部上方区域，并从上方开始随弯曲增加逐渐向下方扩展，直至贯穿皮质骨外径导致表观断裂发生。应用XFEM与CDM模拟所预测断裂模式较为相似，与应用EIF模拟所得断裂模式存在一定差异。图4B，C显示主裂纹均出现在股骨中部下方区域，并且随弯曲增加，裂纹从下至上逐渐扩展，直至贯穿皮质骨外径，导致结构发生表观断裂，该断裂模式与图4D所示的实验断裂模式较为相似。应用XFEM与CDM模拟所得断裂模式的不同之处则在于损伤单元数量差异。应用CDM进行模拟，只要达到损伤条件，单元便会发生损伤，因此从图4B中可以看到在约束与加载区域，均出现一定数量的绿色损伤单元，而应用XFEM模拟，由于整个股骨有限元模型均设置为富集区域，根据XFEM仿真特点，损伤单元只出现在最先满足损伤条件的区域，因此从图4C中可以看到只有极少数损伤单元最初出现在上方刚性圆柱附近，但由于这些损伤单元无法达到失效条件，因此在股骨中部下方区域又重新出现损伤单元，并且随着载荷增加，正是这些损伤单元逐渐演化为失效单元，最终导致结构发生表观断裂。

参考文献

[1] 宫赫,张萌,邹珊珊.肌骨系统生物力学建模2021年研究进展[J].医用生物力学,2022,37(1):18-26.
[2] 李木生,林梓凌,郑利钦,等.基于断裂力学的骨质疏松性大鼠股骨颈骨折有限元仿真[J].医学研究生学报,2019,32(9):915-919.
[3] BOWMAN L, LOUCKS AB. In vivo assessment of cortical bone fragility. Curr Osteoporos Rep. 2020;18(1):13-22.
[4] BALA Y, ZEBAZE R, SEEMAN E. Role of cortical bone in bone fragility. Curr Opin Rheumatol. 2016;27(4):406-413.
[5] 安兵兵,李凯,张东升.皮质骨断裂力学行为的实验研究[J].力学学报,2010,42(6):1164-1171.
[6] 王伟军,刘杰,刘军,等.基于不同应变判定准则模拟皮质骨断裂的准确性[J].中国组织工程研究,2022,24(24):3828-3833.
[7] MACNEIL JA, BOYD SK. Bone strength at the distal radius can be estimated from high-resolution peripheral quantitative computed tomography and the finite element method. Bone. 2008;42(6):1203-1213.
[8] KRAIEM T, BARKAOUI A, MERZOUKI T, et al. Computational approach of the cortical bone mechanical behavior based on an elastic viscoplastic damageable constitutive model. Int J Appl Mech. 2020;12(7):2050081.
[9] CRONIN DS, WATSON B, KHOR F, et al. Cortical bone continuum damage mechanics constitutive model with stress triaxiality criterion to predict fracture initiation and pattern. Front Bioeng Biotechnol. 2022;10:1022506.
[10] SALEM M, WESTOVER L, ADEEB S, et al. Prediction of fracture initiation and propagation in pelvic bones. Comput Method Biomec. 2022;25(7):808-820.
[11] YADAV RN, UNIYAL P, SIHOTA P, et al. Effect of aging on microstructure and fracture behavior of cortical bone as determined by experiment and Extended Finite Element Method (XFEM). Med Eng Phys. 2021;93: 100-112.
[12] KUMAR A, GHOSH R. A review on experimental and numerical investigations of cortical bone fracture. Proc Inst Mech Eng H. 2022; 236(3):297-319.
[13] LI JW, GONG H. Fatigue behavior of cortical bone: a review. Acta Mech Sinica. 2020;37(3):516-526.
[14] FANG J, GAO JZ, GONG H, et al. Multiscale experimental study on the effects of different weight-bearing levels during moderate treadmill exercise on bone quality in growing female rats. Biomed Eng Online.2019;18:33.
[15] DAI XF, FAN RX, WU HJ, et al. Investigation on the differences in the failure processes of the cortical bone under different loading conditions. Appl Bionics Biomech. 2022;2022:3406984.
[16] GINER E, BELDA R, ARANGO C, et al. Calculation of the critical energy release rate Gc of the cement line in cortical bone combining experimental tests and finite element models. Eng Fract Mech. 2017; 184:168-182.
[17] 张睿.老龄骨退化及骨质疏松防治的多尺度试验与数值仿真研究[D].长春:吉林大学,2014.
[18] JIA ZB, GONG H, HU SM, et al. Influence of design features of tibial stems in total knee arthroplasty on tibial bone remodeling behaviors. Med Eng Phys. 2017;48:103-113.
[19] FAN RX, LIU J, JIA ZB, et al. Determination of a tissue-level failure evaluation standard for rat femoral cortical bone utilizing a hybrid computational-experimental method. Proc Inst Mech Eng H. 2018; 232(1):80-89.
[20] KUMAR A, SHITOLE P, GHOSH R, et al. Experimental and numerical comparisons between finite element method, element-free Galerkin method, and extended finite element method predicted stress intensity factor and energy release rate of cortical bone considering anisotropic bone modelling. Proc Inst Mech Eng H. 2022;233(8):823-838.
[21] HAMBLI R, ALLAOUI S. A robust 3D finite element simulation of human proximal femur progressive fracture under stance load with experimental validation. Ann Biomed Eng. 2013;41(12):2515-2527.
[22] GAZIANO P, FALCINELLI C, VAIRO G. A computational insight on damage-based constitutive modelling in femur mechanics. Eur J Mech A-Solids. 2022;93:104538.
[23] 李兴国,安兵兵,张东升,等.非自相似分层次结构对骨组织断裂力学性能分析[J].医用生物力学,2021,36(1):30-35.
[24] FAN RX, GONG H, ZHANG XB, et al. Modeling the mechanical consequences of age-related trabecular bone loss by XFEM simulation. Comput Math Method M. 2016;2016:3495152.
[25] GAO X, CHEN MH, YANG XG, et al. Simulating damage onset and evolution in fully bio-resorbable composite under three-point bending. J Mech Behav Biomed Mater. 2018;81:72-82.
[26] CHEN YG, WANG WS, LI X. Fracture analysis of cortical bone under the condition of cement line debonding and osteon pullout. Int J Biomath. 2018;11(2):1850023.
[27] KHOR F, CRONIN DS, WASTON B. Importance of asymmetry and anisotropy in predicting cortical bone response and fracture using human body model femur in three-point bending and axial rotation. J Mech Behav Biomed Mater. 2018;87:213-229.
[28] LAURENTC, BOHME B, MENGONI M, et al. Prediction of the mechanical response of canine humerus to three-point bending using subject-specific finite element modelling. Proc Inst Mech Eng H. 2016;230(7):639-649.

引言

皮质骨结构的断裂特性将直接决定整骨结构的力学性能与骨折风险[1]。因此，无论是探究预防还是治疗骨折的有效方法，均需掌握皮质骨结构在不同载荷环境下的断裂进程[2]。基于此，已有较多实验通过观测皮质骨结构从受力变形到断裂的全过程，描绘载荷-位移曲线，计算断裂参数，极大促进了皮质骨结构的断裂特性研究[3-5]。
随着断裂力学与有限元仿真技术的发展，针对皮质骨结构的断裂模拟研究也逐渐增多。当前断裂模拟方法主要集中在以下3种：首先为基于单元瞬时失效(element instantaneous failure，EIF)理论的断裂模拟，该方法设置单元达到失效条件时刚度直降为0，通过失效单元积累模拟结构整体断裂，模拟过程相对容易实现且收敛性较好[6-7]；第2种为基于连续损伤力学(continuum damage mechanics，CDM)理论的断裂模拟，该方法设置单元损伤时刚度按一定规律逐渐降低，通过失效单元积累模拟结构整体断裂，由于单元刚度逐渐降低的特性，该方法模拟适用范围较广，但单元损伤演化规则较难确定[8-9]；第3种为基于扩展有限元方法(extended finite element method，XFEM)的断裂模拟，该方法通过使用特别位移函数扩展节点自由度，并允许单元内不连续场存在，同时设置单元力学状态存在损伤起始与演化两个阶段，是较为有效的断裂模拟方法[10-11]。
上述数值模拟方法各有利弊，目前针对不同皮质骨结构的断裂模拟也各有应用，但却鲜有研究对比探究这3种数值方法的断裂模拟精度，以至于即使针对同一皮质骨结构，采用不同方法进行断裂模拟所得结果也可能存在差异[12-13]。此次研究拟探究应用这3种数值方法模拟皮质骨结构弯曲断裂的准确性与适用性。首先获取大鼠股骨样本并进行三点弯曲实验；然后，基于微观影像数据建立大鼠股骨有限元模型，分别采用3种数值方法进行断裂模拟；通过将预测结果与实验数据对比，判定各模拟方法的准确性，同时针对所预测断裂参数进行对比分析，揭示仿真结果差异的原因，并以此确定适合皮质骨结构弯曲断裂的模拟方法，为提高皮质骨结构断裂模拟精度提供理论依据。中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程

材料方法

1.1 设计大鼠股骨样本三点弯曲实验与大鼠股骨有限元模型在3种数值仿真方法下的断裂模拟。
1.2 时间及地点三点弯曲实验于2022年7月在北京航空航天大学生物与医学工程学院力学实验室完成，断裂模拟于2022年9-12月在吉林化工学院航空工程学院虚拟仿真实验室完成。
1.3 材料仿真模拟对象为4个大鼠股骨有限元模型。实验对象为4根大鼠股骨，来自4只3月龄雌性Wistar大鼠，由吉林大学实验动物中心统一提供，许可证号：2018-238，实验开始前测得4只大鼠体质量为355.61-377.15 g。大鼠经麻醉后处死，剔除股骨附近肌肉与软组织，取其左侧股骨。该操作实施符合《赫尔辛基宣言》和吉林大学第一医院对研究的相关伦理要求，批准号为2018-238。实验过程遵循了国际兽医学编辑协会《关于动物伦理与福利的作者指南共识》和本地及国家法规。
1.4 方法
1.4.1 三点弯曲实验描述首先针对大鼠股骨样本进行微观影像扫描(Skyscan 1076，Skyscan，Belgium)，而后将其依次置于电子试验机进行三点弯曲实验。实验弯曲跨距设置为
20 mm，试验机压头以0.5 mm/min的速度匀速向下，实行准静态加载[14]。加载直至股骨中段皮质骨区域发生完全断裂为止，分别记录下各样本的载荷-位移曲线，样本加载示意图如图1所示。

1.4.2 股骨有限元模型建立首先将大鼠股骨微观扫描影像导入至MIMICS软件，以重构股骨三维几何模型，经过优化整理后将三维几何模型导入至ABAQUS软件，应用C3D8单元建立三维大鼠股骨有限元模型[15]。由于股骨在三点弯曲载荷作用下发生表观断裂位置为中段皮质骨区域，所以此文断裂模拟对象实际为大鼠股骨皮质骨结构。因此所建立股骨有限元模型，只包含外部皮质骨结构，并未建立内部松质骨结构。为模拟三点弯曲实验边界条件，在股骨模型上方与下方分别建立刚性圆柱，其中3根刚性圆柱与股骨模型的相对位置与实验一致，上方刚性圆柱与股骨模型之间的相互作用关系设置为无摩擦接触，下方刚性圆柱与股骨模型之间的相互作用关系设置为绑定连接[16]。在上方刚性圆柱参考点处施加轴向位移载荷，同时约束下方刚性圆柱所有自由度，以模拟弯曲工况，具体边界条件示意图如图2所示。

经课题组前期针对3月龄Wistar大鼠股骨皮质骨的研究，通过对皮质骨样本进行纳米压痕实验测得皮质骨组织的纵向和横向弹性模量。基于此，对所建立股骨有限元模型的纵向与横向弹性模量分别赋予10 383 MPa与8 769 MPa，泊松比设置为0.3[17-18]。同时3月龄Wistar大鼠皮质骨组织的临界失效应变以及临界能量释放率也已在前期研究中得到，即皮质骨材料的拉伸失效应变阈值为2.6%，压缩失效应变阈值为4.25%，在弯曲载荷下失效的临界能量释放率为0.16 N/mm[19-20]。基于此，3种数值方法进行断裂模拟所需材料输入参数均为已知。

1.4.3 基于EIF理论的断裂模拟基于EIF的断裂模拟通过ABAQUS软件中的子程序USDFLD实现。在程序中设置一个与单元刚度相关联的变量，此变量初始值为0，随着载荷增加，当单元受弯矩产生的最大或最小主应变达到皮质骨材料临界拉伸或压缩失效应变阈值时，该变量立刻变为1，致使单元刚度骤降为0，并发生失效，随着弯曲载荷增大，当失效单元达到一定比例时，结构发生整体断裂失效[6-7]。
1.4.4 基于CDM理论的断裂模拟基于CDM的断裂模拟通过ABAQUS软件中的子程序UMAT实现。在程序中设定一个损伤变量，当单元发生损伤时，即单元受弯矩产生的最大或最小主应变大于皮质骨材料临界拉伸或压缩失效应变阈值时，损伤变量开始按照一定规律增加，当其增加到0.999时，单元失效。随着弯曲载荷增大，失效单元逐渐增多，当其积累到一定程度时，结构无法有效承载，发生整体断裂。皮质骨结构在准静态载荷作用下基于损伤变量的应力-应变关系如下[21]：
σ = Cdε (1)
Cd = (1-D)C (2)
式(1)中σ为单元应力张量，Cd为单元损伤弹性矩阵张量，ε为单元应变张量；式(2)中D为单元损伤变量，C为单元完整弹性矩阵张量。
应用CDM理论模拟骨结构断裂时，关键问题是描述单元的渐进损伤演化规律。由于皮质骨结构的破坏过程主要基于应变准则控制，因此损伤演化规律依据文献选取为[22]：
D = 0 (ε≤εf) ；D = 1-e[1-(ε/εf)] (ε > εf) (3)
式(3)中ε为单元最大或最小主应变，εf为皮质骨材料临界拉伸或压缩失效应变阈值。
1.4.5 基于XFEM的断裂模拟基于XFEM的断裂模拟中，设置有限元模型中所有单元均为富集单元，通过在ABAQUS软件中编制UDMIGINI子程序判定单元在拉伸或压缩方向的损伤起始，当单元受弯矩产生的最大或最小主应变超过皮质骨材料临界拉伸或压缩失效应变阈值时，单元即进入损伤状态。单元损伤演化规律基于软件自定义的牵引-分离力学行为。随着裂纹开始扩展，裂纹界面处单元胶粘刚度逐渐降低，致使牵引力下降，当胶粘刚度退化为0时，界面损伤单元彻底失去牵引力，导致结构发生整体断裂[23]。
1.5 主要观察指标针对大鼠股骨有限元模型进行弯曲断裂模拟，预测其在EIF、CDM以及XFEM三种数值模拟方法下的断裂载荷与断裂模式，通过与三点弯曲实验数据进行对比，判定各断裂模拟方法的准确性。

讨论

当前针对皮质骨结构的断裂模拟方法主要有基于EIF、CDM理论以及XFEM三种。此文针对同一股骨有限元模型，在赋予相同材料参数的基础上，通过将3种数值方法断裂模拟结果与实验数据进行对比，以确定各方法的模拟准确性，同时针对3组断裂力学参数及断裂模式进行对比分析，揭示仿真结果存在差异原因，探讨3种数值方法的模拟适用范围，并以此找到适合模拟皮质骨结构弯曲断裂的数值方法。
仿真结果显示针对同一有限元模型，3种模拟方法预测所得断裂进程、断裂参数以及断裂模式均存在一定差异，这里拟探讨并解释差异产生原因。首先可以看出3种仿真曲线在弹性阶段存在明显差异，其中EIF与XFEM仿真曲线在弹性阶段的斜率基本一致，CDM仿真曲线的斜率呈下降趋势，说明其表观刚度随加载逐渐降低。造成这种现象的原因主要与各模拟方法设置的单元失效策略不同有关。对于EIF断裂模拟，单元在模拟过程中只存在两种状态，即处于弹性状态或失效状态，单元主应变一旦达到材料失效应变，其刚度则立刻降为0，即从弹性转为失效状态，无法体现出屈服状态，所以曲线在弹性阶段的表观刚度未见降低[6-7]。针对XFEM模拟，虽然其在断裂模拟中设置模型单元包含损伤起始与扩展两个阶段，但软件中同时也自定义裂纹扩展所遵循的牵引-分离行为用弹性关系描述裂纹界面的法向应力、切应力与法向分离、切向分离相关联的本构矩阵，致使单元在损伤演化过程中仍遵循线性刚度退化行为[11，24]。因此，结构所呈现载荷-位移曲线在弹性阶段的表观刚度未见明显降低，同时在加载后期也未表现出明显的塑性屈服。针对CDM模拟，在断裂模拟中也同样设置单元包含损伤起始与扩展两个阶段，但与XFEM不同的是，其单元在损伤演化过程中遵循非线性刚度退化准则。当结构处于弹性阶段时，就已有少量单元发生了损伤，并且单元刚度随损伤变量增加而降低，导致载荷-位移曲线斜率下降，但由于在加载初期弯矩不大，发生损伤单元数量不多，且损伤单元也未失效，仍具备承载能力，所以曲线斜率下降不明显。随着弯矩增大，损伤与失效单元逐渐增多，所以在加载后期曲线斜率下降较为明显，这与另外两种仿真方法所得结果存在一定差异。同时，从图3中也可看出，针对4个股骨样本实验测得载荷-位移曲线斜率在弹性阶段初期缓慢下降，在后期有明显下降趋势。该现象说明皮质骨结构虽然在总体上未呈现出明显的塑性屈服阶段，但内部骨单位可能仍具有一定塑性力学特性，所以股骨样本在三点弯曲载荷作用下呈现出准脆性断裂行为[4，25]。
从图3中还可以看出不同数值方法预测得出的断裂时刻也存在差异，其中EIF模拟发生完全断裂时间明显早于实验数据，CDM与XFEM模拟下的断裂时刻与实验结果较为接近。产生这种现象的原因主要是由于EIF模拟断裂过程中缺乏单元损伤演化阶段，即单元刚度直接从初始值变为0，而在实际情况中，对于皮质骨材料，即使达到材料失效应变阈值发生损伤，内部骨单位也仍具备一定承载能力，而单元瞬时失效的模拟策略，直接将损伤单元定义为失效，显然弱化了结构承载能力，导致断裂提前发生[26]。与此同时，图3显示3种数值方法模拟得到的断裂载荷也存在异同。EIF与CDM模拟所得断裂载荷均明显小于XFEM模拟数值。产生这种现象的原因主要是由于XFEM模拟中结构在发生断裂失效前，载荷-位移曲线斜率始终未变，结构表观刚度未见降低，致使弹性阶段后的断裂载荷处于高值状态，与实验结果不符。EIF断裂模拟理应也由于结构表观刚度未降低而导致断裂载荷高于实验结果，但由于仿真策略弱化了单元承载特性，致使断裂提前发生，所以断裂载荷并未上升到XFEM的模拟结果。
图4显示3种数值方法针对同一有限元模型预测得到的断裂模式存在差异，该现象主要与单元损伤及失效准则设置不同相关。加载初期，模型在上方刚性圆柱附近应变较高，因此3种模拟均显示股骨中段上方区域附近率先出现少量损伤单元。对于EIF模拟，由于单元损伤即失效，立刻失去承载能力，所以失效单元附近，即上方刚性圆柱附近的完好单元需承载更多载荷，由于过载导致这些单元很快也进入失效状态。因此，EIF模拟显示裂纹从股骨中段上方区域逐渐扩展至下方。对于CDM与XFEM模拟，加载初期上方损伤单元并未失效，仍具备承载能力，所以上方损伤单元附近的单元并未过度承载，因此股骨中段上方区域基本未出现失效单元。在加载中期，弯矩主要作用于股骨中部下侧，足够大的弯曲使股骨中段下方单元发生损伤与失效。因此，裂纹从下侧起始逐渐扩展至上侧。同时，根据文献报道，股骨在三点弯曲载荷作用下也是从中段下侧区域开始产生裂纹，并逐渐向上扩展直至断裂，文献所述断裂模式与CDM及XFEM模拟的断裂模式相近[27-28]。
通过上述分析可知，EIF断裂模拟由于未设置单元损伤演化过程，所以更适合模拟脆性断裂，由于此文股骨样本展示出微弱的塑性屈服特性，所以无法准确模拟其断裂进程。XFEM模拟虽然设置了单元损伤与演化过程，但由于软件自定义的牵引-分离行为设定单元遵循线性刚度退化行为，致使曲线表观刚度未见明显下降，导致断裂载荷过大。因此，XFEM相对也更适合模拟脆性断裂。相比之下，CDM则更适合模拟皮质骨结构的弯曲断裂，由于其断裂策略设置了单元损伤演化过程，所以在模拟过程中能够呈现出较微弱的塑性屈服特征，这与皮质骨结构在弯曲载荷作用下的准脆性断裂特性更加符合。因此，从结果对比看出，CDM仿真预测所得载荷-位移曲线以及断裂模式与实验结果一致性更好。
此文在进行断裂模拟的同时也存在一些缺点。首先所研究样本量过小，由于仅选取4个股骨模型进行模拟，无法进行数据统计；然后，此文仅进行了大鼠股骨结构在弯曲载荷下的断裂模拟，还缺乏在其他力学环境如压缩与冲击载荷下的模拟与实验研究，无法对结构断裂模拟方法提出一个综合性、全面性的建议。虽然存在上述缺点，但此文主要目的为探究不同数值方法的模拟精度，以此为提高皮质骨结构弯曲断裂模拟准确性提供理论依据。
结论：此文针对同一皮质骨结构进行不同仿真策略下的弯曲断裂模拟。经与实验结果对比得知，相比于EIF与XFEM，CDM模拟与实验数值更加符合。分析发现仿真与实验差异原因主要与各模拟方法所设置单元失效策略不同相关。EIF断裂模拟中由于设置单元直接从弹性进入失效阶段，无法体现出结构的塑性屈服特性，导致结构断裂起始位置及裂纹扩展路径与实验不符；XFEM断裂模拟虽然设置了单元的损伤起始与演化，但由于其所遵循的牵引-分离模型默认单元遵循线性刚度退化行为，导致断裂载荷过大，与实验结果不符；CDM断裂模拟由于设置单元的非线性连续损伤，其断裂参数与实验结果相对一致。因此，基于CDM理论的数值方法较适合模拟皮质骨结构在三点弯曲载荷下的断裂。中国组织工程研究杂志出版内容重点：人工关节；骨植入物；脊柱；骨折；内固定；数字化骨科；组织工程

不同数值方法模拟皮质骨结构弯曲断裂的准确性

Fracture accuracy on cortical bone structure under bending load using different numerical methods

PDF

可视化

摘要/Abstract

引用本文

使用本文

图/表（结果） 2

参考文献

相关文章 15

引言

材料方法

讨论

文章快阅

延伸阅读

编辑推荐

Metrics

本文评价

[1]	霞黑达·依拉尔江, 尼加提·吐尔逊, 热依拉·库尔班, 白布加甫·叶力思, 陈欣. 三维有限元分析不同特征骨组织中应力的分布规律[J]. 中国组织工程研究, 2024, 28(8): 1277-1282.
[2]	程云忠, 刘玉增, 海涌, 关立, 潘爱星, 张希诺, 陶鲁铭, 李越. 皮质骨轨迹螺钉研究现状与发展趋势的文献计量及可视化分析[J]. 中国组织工程研究, 2023, 27(4): 513-519.
[3]	马赛, 刘波, 李楠, 何达. 机器人辅助与传统徒手置入皮质骨轨迹螺钉在老年骨质疏松性腰椎手术应用中的比较[J]. 中国组织工程研究, 2023, 27(27): 4332-4336.
[4]	陈鉴权, 陈茂水, 吕洲明, 陈荣彬, 余照宇, 刘万鹏, 林新源, 林定坤. 有限元分析皮质骨轨迹螺钉联合椎弓根螺钉固定椎体运动单元的生物力学性能[J]. 中国组织工程研究, 2023, 27(22): 3457-3462.
[5]	鄢来军, 葛海雅, 周斌, 糜大国, 王鹏翔, 李楠. 斜外侧腰椎椎间融合与双侧皮质骨轨迹螺钉置入后腰椎生物力学的有限元分析[J]. 中国组织工程研究, 2023, 27(13): 1969-1974.
[6]	郑永泽, 郑利钦, 何兴鹏, 陈心敏, 李木生, 李鹏飞, 林梓凌. 基于ABAQUS软件股骨颈骨折的扩展有限元建模分析[J]. 中国组织工程研究, 2022, 26(6): 853-857.
[7]	蔡风, 余波, 曾铎, 陈钦灿, 廖琦. 皮质骨轨迹螺钉内固定治疗老年骨质疏松性腰椎疾病[J]. 中国组织工程研究, 2022, 26(3): 403-407.
[8]	周志豪, 阿拉法特·卡哈尔, 王轶希, 刘东山, 希尔艾力·麦麦提, 石文杰, 帕尔哈提·热西提. 传统椎弓根螺钉与改良皮质骨轨迹置钉技术的生物力学性能有限元分析[J]. 中国组织工程研究, 2022, 26(18): 2789-2794.
[9]	林书, 胡豇, 万仑, 唐六一, 王跃, 俞阳. 骨科机器人辅助皮质骨轨迹螺钉内固定治疗腰椎退行性疾病[J]. 中国组织工程研究, 2022, 26(15): 2356-2360.
[10]	任航宁, 居来提·买提肉孜, 帕尔哈提·热西提, 罗辉卿. 皮质骨轨迹椎弓根系统拔出力下腰椎的有限元分析[J]. 中国组织工程研究, 2021, 25(36): 5771-5776.
[11]	魏源标, 郭惠智, 张顺聪. 皮质骨轨迹螺钉固定对相邻节段影响的有限元分析[J]. 中国组织工程研究, 2021, 25(18): 2799-2804.
[12]	王伟军, 刘杰, 刘军, 贾正斌, 范若寻. 基于不同应变判定准则模拟皮质骨断裂的准确性[J]. 中国组织工程研究, 2020, 24(24): 3828-3833.
[13]	宗治国, 刘肃, 马朋朋, 张鑫, 李伟, 苏峰, 张春林, 郭建文, 文毅. 植入不同直径椎弓根螺钉后腰椎体与椎弓根结合部位骨折力学的稳定性[J]. 中国组织工程研究, 2020, 24(18): 2794-2798.
[14]	覃海阔, 罗世兴. 股骨近端不同平面皮质骨厚度值、X射线灰度值与女性髋部脆性骨折的相关性[J]. 中国组织工程研究, 2020, 24(18): 2867-2872.
[15]	王宁, 崔婷婷, 李永奇, 赵彬彬, 仲维剑, 马国武. 自体块状皮质骨修复种植区颌骨缺损[J]. 中国组织工程研究, 2020, 24(17): 2675-2679.