## 4 Discussion
The fracture time in the XFEM and CDM simulations was consistent with the experimental results, but in the EIF simulation was earlier than the experimental fracture time. This phenomenon was mainly caused by the element in the EIF simulation lacking a damage evolution process, namely the stiffness changed instantly from its initial value to zero when it reached the damage initiation criterion. In actuality, the osteon still has a certain load-bearing capacity even damage occurs in cortical bone structure [30]. The EIF method directly specified the damage element as failure, which weakened the structural load-bearing capacity and caused the apparent failure to occur earlier than expected. Meanwhile, the causes of the variations in the fracture load predicted by the three numerical methods were also discussed. The fracture load predicted by the EIF and CDM simulations was comparable to the experimental results but considerably less than the value in the XFEM simulation. This phenomenon was caused by the fact that the slope of the load-displacement curve in the XFEM simulation was not dropped before the apparent failure occurred. The predicted fracture load in EIF simulation should express an excessive value because the slope of the curve in the elastic stage was also not decreased; however, the fracture load was lower than the experimental result. This phenomenon was caused by the fact that the load-bearing ability in the cortical bone structure was weakened due to the absence of damage evolution stage in the element in the EIF simulation process, which led to the premature occurrence of the apparent failure. The differences in the fracture patterns predicted by the three numerical methods can be seen in Figure 5. The reason for the differences in the fracture patterns was mainly related to the various criteria in the element damage and failure processes. The principal strain was high near the indenter at the initial loading, so all three simulations demonstrated that a limited number of damaged elements appeared in the upper side of the middle femur near the indenter. Because the damaged element lost its load-bearing capacity in the EIF simulation, the complete elements near the failed elements must bear more load, which could cause these complete elements to go into failure state quickly due to overload. Therefore, the crack propagated from the upper to lower sides of the middle femur in the EIF simulation. In the CDM and XFEM simulations, at the beginning of loading, the damaged elements in the upper side of the middle femur did not fail and still retained carrying capacity, so the complete elements close to the damaged elements were not overloaded and did not fail. In the middle of loading, the bending moment was transferred from the constraint location to the lower side of the middle femur, and the bending was sufficient to lead to failure in the damaged elements. Therefore, the crack propagated from the lower to upper sides of the middle femur. Furthermore, it has been reported that the femur fracture under three-point bending load is mainly because the continuous propagation of the lower surface crack in the transverse and longitudinal directions cut off the central section of the femoral shaft, which is consistent with the failure processes in the CDM and XFEM simulations [31; 32]. The above discussion concluded that the EIF numerical method was more suitable for simulating the brittle fracture because the damage evolution process in the element was not set up and the apparent stiffness was not decreased in the load-displacement curve. Though the damage evolution process was set in the XFEM simulation, the available traction-separation model in the ABAQUS assumes the linear stiffness degradation behavior during the damage evolution in the element, so the apparent stiffness was also not decreased in the elastic stage, which caused the excessive fracture load. Therefore, the XFEM simulation may also be more suitable for simulating the brittle fracture.