Study on validation method for femur finite element model under multiple loading conditions

Acquisition of accurate and reliable constitutive parameters related to bio-tissue materials was beneficial to improve biological fidelity of a Finite Element (FE) model and predict impact damages more effectively. In this paper, a femur FE model was established under multiple loading conditions with diverse impact positions. Then, based on sequential response surface method and genetic algorithms, the material parameters identification was transformed to a multi-response optimization problem. Finally, the simulation results successfully coincided with force-displacement curves obtained by numerous experiments. Thus, computational accuracy and efficiency of the entire inverse calculation process were enhanced. This method was able to effectively reduce the computation time in the inverse process of material parameters. Meanwhile, the material parameters obtained by the proposed method achieved higher accuracy.


Introduction
At present, biomechanics of bone is one of the most popular research topics in the field of injury biomechanics. Femoral fracture is the major injury mechanisms of human lower limb. Such severe injuries can lead to long-term disability for victims [1]. Therefore, lots of femoral experiments have been conducted to investigate femoral injury mechanism in different loading conditions [2][3][4]. And many femur finite element (FE) models have been developed and validated [5][6][7][8][9][10][11][12], including commercial FE model of Total Human Model for Safety (THUMS) [13,14] and Global Human Body Models Consortium (GHBMC) [15], and the parametric human femur model [16]. By virtue of a femur FE model, its biomechanical responses are simulated under various loading conditions to help research on related tissue injury mechanism and development of protection countermeasures.
However, accurate material parameters are important prerequisites for FE models to arrive at both the reliable analysis and predictive results. Without any doubt, experiment is an effective means to provide data for estimating material constitutive parameters of the femur. Nevertheless, mechanical parameters obtained by material testing based on limited samples are rather random and the corresponding experimental results can not entirely represent the material characteristics of the population. The femur FE model established by trial-and-error method is difficult to match the specific biomechanical test results, particularly difficult to conform to the biomechanical test results under a variety of loads. In order to circumvent the disadvantages of the trial-and-error methods, an inverse engineering method is adopted to further characterize mechanical response of the material [17][18][19]. In this way, the predicted behaviors of femur such as deformation, fracture, etc. can be more consistent with those observed in the experiment, thus more effectively solving the problem that biomechanical theories are divorced from engineering practices. The inverse method provides a useful alternate approach to reasonably identify the material parameters of femur.
Quasi-static and dynamic 3-point bending mechanical experiments have been carried out to study the femur biomechanical response under various loading conditions in the literatures [20][21][22]. Moreover, numerous femur FE models have also been developed and reported published [5-16, 23, 24]. However, the inverse engineering method employed for determining the material parameters of femur is rarely studied and the inverse procedures of them are also indefinite. Guan et al. (2011) have used optimization methodology and specimen-specific FE models for material property identification of rat skull, each experimental force-displacement curve was matched to the simulated curve using the optimal material parameters [19].
Using this approach human femur constitutive material parameters were selected as design variables, and further optimized in an optimization process. Using the obtained material parameters, the femur FE model was validated under multiple loading conditions (different quasi-static loading directions and dynamic loading locations). In this study, the FE simulations were performed using LS-DYNA non-linear code (LSTC, Livermore, CA, USA).

Three-point bending experiments for femur
One of the major experiments for determining the biomechanical behaviour of femur is 3-point bending testing. According to different load rates, experiments are divided into the quasi-static and the dynamic loadings. For both test conditions, the proximal and the distal ends of femur specimen are fixed. For the quasi-static tests, the loading can be applied along the A-P (from the anterior to the posterior) or the L-M (from the lateral to the medial) directions at the mid-shaft location. In the dynamic 3-point bending experiments, the load is applied at 1/3 of the proximal part of femur, the mid-shaft, and 1/3 of its distal part separately.  summarized the femur injury tolerance characteristics under static loading conditions, without force-displacement curves in the literature. Yamada (1970) conducted the quasi-static 3-point bending biomechanical experiment of femur by placing both ends of the femur onto a rigid platform and gradually applying load at the mid-span by an impactor. Both loading conditions along the A-P and L-M directions were carried out [1]. The time histories of the applied loads and displacement of the impactor were recorded. The detailed experimental data were also published in the literature, which are used extensively in femur FE model validation. Kerrigan et al (2003) [25,26] performed dynamic 3-point bending experiments on femurs where the loading locations were at the middle shaft of femur (Mid-shaft.), at the 1/3 of femur proximal part (Prox. 1/3) and 1/3 of distal part (Distal 1/3) respectively. The associated experimental facilities, loading and boundary conditions were all illustrated and documented in detail. Therefore, their test data are used as the basis of material parameters identification and model verification in this study.

Dynamic 3-point bending experiment
The experimental test setup is shown in reference [25,26]. In the experiment, two square metal boxes were used to fix both ends of the femur specimen, and curved plates were installed at the bottom of these boxes and placed on the flat. Between the flat and the curved plates, lubricating oil was applied to allow free rotation of the square metal boxes. In this manner, influences of boundary conditions on biomechanical responses of the femur could be minimized. A solid cylinder of a diameter of 12.7 mm was used as an impactor at a speed of 1.2 m/s. The femur specimen was positioned with the lateral side on the top and the medial side on the bottom. The load was applied on the femur by the impactor on the middle shaft of the femur specimen, 1/3 of femur proximal part and 1/3 of femur distal part, respectively. Moreover, a load cell was installed at the top of the impactor to measure the force.

Femur FE Model
Femoral bone is composed of the outer cortical bone and the inner cancellous bone. In femur FE models, elastic-plastic materials and viscoelastoplastic materials are selected to simulate behaviors of the cortical bone and the cancellous bone, respectively. The geometric modeling of the femur was generated from the femur CT imaging. The femur was divided into one shaft with proximal and distal ends which was shown in Figure 1, and different thicknesses and materials parameters were set. Quadrilateral shell element and hexahedral element were used to simulate cortical bone and cancellous bone, respectively. The cortical bone has higher strength than the cancellous bone. Thus, the cancellous bone at femoral diaphysis has a slight influence on the impact response. Therefore, the elastic-plastic material parameters of the cortical bone at femoral diaphysis were studied.

Figure 1. Femur model
According to the 3-point bending mechanical experiment setup, five FE simulation models were developed, including quasi-static test with loading along the A-P direction, quasi-static test with loading along L-M direction, dynamic test with loading at the mid-shaft location, dynamic test with loading at the proximal 1/3 location, and dynamic test with loading at the distal 1/3 location.
The FE models were generated based on experimental conditions. A rigid cylinder with a diameter of 25mm was used to apply load at the mid-shaft of femur at a constant speed of 0.01m/s, while releasing only the degree of freedom of its movement along the vertical direction, to simulate the quasi-static experiments. For the dynamic experiment simulations, a rigid cylinder with a diameter of 12.7mm moving at a speed of 1.2m/s and impacting at proximal 1/3, mid-shaft, and distal 1/3 of the femur, respectively. The time history curves of contact force and the displacement of the impactor were outputted in all simulations.

Validation Method
There were five loading conditions for femoral 3-point bending biomechanical tests (quasi-static with A-P and L-M loading direction, dynamic L-M loading direction at proximal 1/3, mid-shaft, and distal 1/3 of the femur). Tens of thousands of FE simulation calculations should be executed if the genetic algorithm was adopted to optimize the inverse problem model directly, which would lead to low efficiency of inverse calculation. Therefore, a method of the femur material parameters identification based on the sequential response surface was adopted in this study.

Sequential Response Surface Method
Traditional polynomial response surface that constructs approximation model inside the entire design space frequently give rise to dissatisfactory precisions. Here the Sequential Response Surface Method (SRSM) was employed to select sampling point in the initial design subdomain first, followed by polynomial response surface model approximation specific to the mathematical model (e.g., FE model) of these sampling points. In this way, optimization design points of such an initial subdomain could be acquired and then used as the centre of the next subdomain. Through movement and scaling, etc. (Figure 2), new subdomains were updated, and polynomial approximation and optimization was conducted. Successively and iteratively, results approaching the globally optimal solution were obtained. Based on continuous movement and scaling of the design subdomains, the optimization problem finally became a problem to solve an extremum inside a smaller domain. As a result, a response surface approximation model with a higher precision and the optimal design points could be achieved. Such a new approach was featured with fast convergence and high approximation precision.  Figure 2. The subdomain center ranking k+1 is assumed to be the optimal design point of a subdomain at k; and, the size of k+1 subdomain space remains unchanged or reduces correspondingly. Relation of the area changes between rankings k and k+1 is given as follow: where, ( ) x , and is usually taken to be 0.5; γ, ranging from 0 to 1, is the value of λ i when the approximate optimal design point * ( ) k i x is located at upper/lower bound of the subdomain at k. To ensure that the scaling of all variables remains consistent in renewal processes of subdomains of every generation, λ is taken as the partial factor, which is also known as the variation factor, of all variables and its expression is: Now, the upper and lower bounds of a design variable ranking i in a design subdomain at k+1 are denoted as: Therefore, the optimization is ultimately changed into a problem of solving the extremum within a small area by means of constant design subdomain movement and scaling. Furthermore, both response surface approximation model of high precision and optimal design points are acquired. In addition to effectively eliminate noises, this method can also improve the convergence performance of the modified algorithm.

Inverse Strategy for Material Parameters Identification Based on FEM-SRSM
The response surface model of the objective function was constructed by using finite element methodsequential response surface method (FEM-SRSM) combined with experimental data of different loading conditions in different design subspaces. Then genetic algorithm is used to optimize the polynomial response surface model to improve the computational efficiency of the material parameters inverse process. The inverse process, that is, the process of approximating the simulation curve of the finite element model to the experimental force-displacement curve measured.
Mean square error of experimental curves and simulated curves was used as the objective function during this inverse process. The value of the objective function in each iteration was calculated, so as to determine whether the optimization requirements were met. If the optimization requirements were not reached, the response surface model would continually update by the movement and scaling of the interest domain until the convergence was reached and the material parameters were output. where, K refers to the number of experimental curves, k = 1, 2, ··· K. m = 1, 2, ··· P k , P k is the number of calculating points on the curve. X is the design variable. W m is the weight coefficient. f m (X) is the calculated value of response surface approximation model. G m is the value of experimental test point. S m is the residual proportion coefficient.
Initially, the experimental force-deformation curves from quasi-static and dynamic tests used as the matching target were selected, and then the sample points were taken in the entire design space by using the D-optimal experimental design to construct an approximate model of femur responses under each loading condition. Then the objective functions were formed according to Eq. (4). Finally, the optimization was performed by genetic algorithm and the approximate model were continually updated by moving and scaling the domain of interest until the optimization process converged.

Quasi-static Material Parameters Identification
Quasi-static material parameters inverse model mainly aimed at acquiring the fundamental material parameters of the femur, including elastic modulus, yield stress, and tangent modulus.

Figure 3. Comparisons of simulation results and test results for quasi-static loading conditions
The optimized force-displacement curves derived from quasi-static 3-point bending simulation and experimental curves are shown in Figure 3. The optimized simulation force-displacement curves (the black solid lines) are very close to the test curves (the red dotted lines), which indicated that the optimized femoral material parameters can well characterize the biological fidelity of the femoral FE model. Figure 3 also displayed data (the green lines) related to FE model verification of femur in recent years [23,24]. According to these data, the material parameters identification method could be adopted to acquire more superior biological fidelity of a femur FE model in quasi-static loading conditions. In particular, it should be noted that the optimization verification process can achieve two kinds of loading conditions have a good simulation results. The optimized elastic modulus, yield stress and tangent modulus of the femur shaft under quasi-static 3-point bending loading condition are found to be 17.24GPa、59.93MPa, and 3.11GPa, respectively.

Dynamic Material Parameters Identification
The dynamic material parameters inverse model was primarily to identify elasticity modulus, yield stress, tangent modulus and strain rate parameters (C and P). respectively. An examination of Figure 4(b) reveals that the predicted results of the case where loading was applied at the midshaft is also compared with those obtained from simulation results by Takahashi et al. (2003) [23] and Untaroiu et al. (2006) [17]. Obviously, simulated result obtained from this study is better beyond the displacement of 8mm. Likewise, simulated result for the case of loading at the distal 1/3 location, when compared with that reported by Takahashi et al. (2003) [23] as shown in Figure 4(c), exhibits better fit the test data. Hence, the femur FE model developed in this study provided with a better capability in predicting dynamic behavior of femur under 3-point bending, thus achieving more biofidelic response. From Figure 4, the optimized elastic modulus, yield stress, tangent modulus, stain rate parameters C and P of the femur shaft are found to be 14.43GPa、87.21MPa、 0.331GPa、154.5 and 9.90, respectively.

All Loading Conditions Material Parameters Identification
During the material parameters identification under all loading conditions (combined the quasi-static and dynamic loading conditions), sensitivity analysis was carried out for femur cortical bone material parameters. Strain rate parameters C and P with rather minor influences after screening were then removed. Therefore, elasticity modulus, yield stress and tangent modulus were finally chosen as design variables. The optimized elastic modulus, yield stress and tangent modulus of the femur shaft were found to be 14.76GPa、80.38MPa and 0.438GPa, respectively.
The comparisons of simulation results and test results for combined loading conditions are presented in Figure 5. A review of these figures indicates that the femur model in simulating quasistatic and dynamic loading conditions yielded rather favorably biological fidelity results when judged from comparisons with available experimental data. Therefore, such model thus developed and validated can be readily applicable to study of biodynamics of femur.  The results for femoral cortical bone material parametric values obtained through optimization using data from quasi-static tests only, dynamic tests only, and combined quasi-static and dynamic tests are shown in Table 1. It shows that material parametric values are different for the quasi-static and dynamic bending conditions if their respective material inverse models are used. It also revealed that between the quasi-static and dynamic cases, the greatest difference occurs in tangent modulus. This may due to the strain-rate effect of C and P parameters in the dynamic condition. Further, it is noted that material parametric values for the femoral cortical bone are nearly similar between the dynamic and the combined loading conditions. Since this paper used an identical femur FE model to fit experimental curves for five loading conditions at the same time, the established model can be confidently applied to diverse loading conditions in future potential applications.

Conclusion
In this paper, a femur mathematic model validation method that combines FE model, sequential response surface model and genetic algorithm together was proposed. As a result, the model validation problem was translated into the problem of material parameters optimization. As such, this method was adopted to obtain optimal elasticity modulus, yield stress, tangent modulus, and strain rate parameters C and P. Quasi-static simulation results along A-P and L-M loading directions of femur, dynamic simulation results at the proximal 1/3, the mid-shaft and the distal 1/3, as well as the combined loading condition simulation results were all in good agreement with experimental results. Clearly, material parameters obtained by inverse engineering method presented in this paper met satisfactorily biological fidelity requirements of biomechanical FE model. In addition, it was also able to effectively minimize experiments for material parameters identification, thus, saving costs and improving efficiency.