## 4 Discussion  
This paper aimed to evaluate the prediction accuracy and applicability of cortical bone fracture in the three different numerical methods. The simulation results predicted by different numerical methods were compared with the experimental data to determine the respective prediction accuracy on the basis of the same material parameters and FE model. The fracture parameters and fracture patterns obtained from the three simulations were compared and analyzed to reveal the cause of the variations in the predicted results and to find a suitable numerical method for simulating cortical bone fracture. The simulation results showed visible differences in the failure processes, fracture parameters, and fracture patterns predicted by the three simulation methods for the same FE model. First, the cause of the variations in the elastic stage of the load-displacement curves was discussed. The slope of the curve in the elastic stage in the EIF and XFEM simulations was the same and almost did not drop, but the CDM simulation displayed a decreasing trend. This phenomenon was mainly related to the different element failure strategies in the numerical simulation methods. In the EIF simulation, the element included two states, i.e., the elastic and failure states. As the principal strain in the element reached the material  
Figure 5: Comparison of the fracture patterns predicted in the three numerical method simulations and performed in the experiment: (**a**) fracture simulation based on the EIF method; (**b**) fracture simulation based on CDM method; (**c**) fracture simulation based on XFEM; (**d**) experiment. failure strain, the element stiffness immediately dropped to 0. The element did not go through the damage evolution process and failure occurred directly. Thus, the apparent stiffness in the elastic stage was not decreased. Though the element would experience the damage initiation and propagation stages in the XFEM simulation, the traction-separation behavior followed by the crack extension is customized in the ABAQUS software to describe the intrinsic matrix of the normal and tangential stresses at the crack interface associated with the normal and tangential separation by elastic relations, resulting in the linear stiffness degradation behavior during the damage evolution stage [26; 27]. Thus, the apparent stiffness in the curve was also not decreased in the elastic phase. In the CDM simulation, the element also included the damage initiation and propagation stages. However, the elements follow the nonlinear stiffness degradation criterion during the damage evolution, unlike XFEM. In the elastic phase, a few damaged elements appeared, and the element stiffness gradually decreased with the increasing damage variable, resulting in a slightly dropping trend in the slope of the curve. Because the bending moment was not great at the beginning of the loading, the number of damaged elements was small, and the damaged elements had not yet entered the failure stage and still had the bearing capacity, the dropping trend in the slope was not obvious. The number of the damaged and failed elements gradually increased with bending load, so the slope of the curve in the later elastic stage fell more visibly, which was different from the results obtained by the other two numerical methods. Additionally, Figure 4 demonstrated that the slope displayed a tendency of slowly dropping at the beginning and obviously decreasing at the middle and end of loading in the experimental load-displacement curves. This phenomenon indicated that despite the mechanical behavior in the cortical bone structure and general lack of a visible yielding phase, the osteon in cortical bone may nevertheless exhibit plastic mechanical properties, which is why the quasi-brittle fracture occurred in the cortical bone structure under three-point bending load [28; 29]. Second, the causes of the variations in fracture time predicted by the three numerical methods were discussed.