# Micromechanical modeling of ductile fracture of human humerus

Jamila Rahmoun, Hakim Naceur, Pascal Drazetic, Christian Fontaine

► To cite this version:

Jamila Rahmoun, Hakim Naceur, Pascal Drazetic, Christian Fontaine. Micromechanical modeling of ductile fracture of human humerus. Journal of Mechanical Engineering and Sciences, 2020, 14 (2), pp.6952-6960. 10.15282/jmes.14.2.2020.32.0544. hal-03442148

HAL Id: hal-03442148

https://uphf.hal.science/hal-03442148

Submitted on 25 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ORIGINAL ARTICLE

# Micromechanical modeling of ductile fracture of human humerus

J. Rahmoun<fnref n="1" />, H. Naceur<fnref n="1" />, P. Drazetic<fnref n="1" /> and C. Fontaine<fnref n="2" />

Phone: +33327511412

**ABSTRACT** – This paper deals with the formulation, development and validation of a newly developed micromechanical-based model for the modeling of the nonlinear ductile fracture of human humerus. The originality of the present works concerns the coupling between the micromechanical formulation based on the Mori-Tanaka homogenization scheme for cylindrical voids and the Marigo nonlinear ductile damage model based on the porosity growth. The proposed model was implemented as a User Material UMAT within the explicit dynamic software LS-DYNA and validated by numerical and experimental analysis conducted by a drop tower impact of human humerus. The outcome of the proposed multi-scale model appears to correctly predict the general trends observed experimentally via the good estimation of the ultimate impact load and the fracture patterns of the human humerus.

ARTICLE HISTORY

Revised: 22nd Apr 2020

Accepted: 30th Apr 2020

KEYWORDS

Micromechanic;

biomechanic;

multiscale modeling;

fracture;

impact;

finite element.

INTRODUCTION

Understanding the physical mechanisms of bone fracture represents a major challenge in biomechanics, since it allows the enhancement of injury criteria commonly used by European New Car Assessment Programme (Euro NCAP) for the safety of vehicle passengers and pedestrians [1]. It can also deliver a follow up of athlete's safety during their trainers avoiding risk zones of injury especially in contact sports [2]. This knowledge is essentially based on the use of the numerical models, whose prediction is assessed through the development of high resolution medical imaging and simulation software [3]. Among these models, the anthropometric test device (ATD) commonly used for crash-test or other more detailed local models simulating the interaction between bone tissue and clinical equipment such as prostheses. Their advantage lies mainly in the diversity of configurations and loading conditions and hence the optimization of time and the total benchmark cost. However, one can observe that the material constitutive laws used are often derived from the experimental characterizations carried out at the macroscopic scale ignoring the bone microarchitecture [4]. A micromechanical approach is more interesting to improve the limitations of these models [5]. Concerning the studies investigating the effects of geometrical parameters and their combinations, we can cite the study of Sam Daliri et al. [6, 7] who combined numerical and statistical analysis for the prediction of stress distribution and buckling of the cylindrical shell structures [8], analogous to humerus bone.

The present investigation shows the theoretical formulation, the numerical development and the experimental validation of a ductile damage and fracture model applied to the human humerus bone in the thermodynamics framework [9]. As the bone architecture has a matrix-inclusion morphology, the Mori-Tanaka homogenization scheme [10] was proposed for the modeling of the elastic behavior [11, 12]. For the non linear behavior, the macroscopic tangent operator was formulated. The strain rate effects on the humerus failure was considered by means of the standard model of Johnson-Cook as a preliminary trial [13]. The resulting micromechanical damage model was then implemented as User Material subroutine (UMAT) within the explicit dynamics code LS-DYNA [14].

The numerical model was validated by a drop tower test on human humerus for the determination of the mechanical response and the damage growth until complete failure of humerus. The proposed approach revealed to be appropriate for fracture analysis of biomechanical hard tissues.

METHODS AND MATERIALS

Micromechanical Formulation of the Damage model

The representative elementary volume of human humerus bone is composed of an elastic solid matrix and cylindrical voids. In the thermodynamic framework, we consider an elastic damage model caused by the void growth in which the porosity represents the damage variable d. As classically, the thermodynamical potential of the solid matrix reads:

$$ \psi = \frac{1}{2} \varepsilon : C : \varepsilon \qquad (1) $$

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Laboratory LAMIH, Polytechnic University Hauts-de-France, 59313 Valenciennes cedx 9, France</footnote>

<footnote n="2">²Laboratory of Anatomy, Hospital University of Lille, 59000 Lille, France</footnote>
where $\underline{\epsilon}$ is the macroscopic strain tensor and where $\mathbb{C}$ is the macroscopic stiffness tensor obtained from the Mori Tanaka model [10, 15]:

$$ \mathbb{C} = \mathbb{C}_s + (\mathbb{P} - \mathbb{C}_s^{-1})^{-1} : [(1-d)\mathbb{I} + d(\mathbb{I} - \mathbb{P} : \mathbb{C}_s)^{-1}]^{-1} \quad (2) $$

where $\mathbb{I}$ is the identity tensor and where $\mathbb{P}$ is the Hill tensor [16] which depends on the shape and orientation of voids and on the solid matrix elasticity tensor $\mathbb{C}_s$ (see [11, 17]).

From the Marigo [9] formulation, we can obtain the following yield function $g$ as a function of the damage variable $d$ and the thermodynamical force $F^d$:

$$ g = F^d - R(d) \le 0 \quad (3) $$

with $R(d) = k + \xi d$. The two damage variables $k$ and $\xi$ can be identified from experimental investigations. The standard normality rule reads:

$$ \dot{d} = \dot{\lambda} \frac{\partial g}{\partial F^d} = \dot{\lambda} \quad (4) $$

where $\dot{\lambda}$ is the damage multiplier calculated by means of the consistency condition $\dot{g} = 0$. Then, by including the evolution law in the stress tensor, we obtain the tangent formulation of the elastic damage law:

$$ \underline{\dot{\epsilon}} = \mathbb{L} : \underline{\dot{\epsilon}} \quad (5) $$

where the tangent operator $\mathbb{L}$ writes:

$$ \mathbb{L} = \mathbb{C} - \frac{1}{H} \left( \frac{\partial \mathbb{C}}{\partial d} : \underline{\underline{\epsilon}} \right) \otimes \left( \frac{\partial \mathbb{C}}{\partial d} : \underline{\underline{\epsilon}} \right) \quad (6) $$

with $H = -\frac{\partial g}{\partial d} = \xi + \frac{1}{2} \underline{\underline{\epsilon}} : \frac{\partial^2 \mathbb{C}}{\partial d^2} : \underline{\underline{\epsilon}}$.

## Strain Rate Effect

The strain rate $\dot{\epsilon}$ can have a significant influence in the case of dynamic loading. Then, to consider the strain rate effect on the humerus bone failure strain, we adopt the simplified Johnson-Cook rupture model [18]:

$$ \epsilon_r^{dyn} = \epsilon_r^{sta} \left( 1 + D_4 \ln \frac{\dot{\epsilon}}{\dot{\epsilon}_0} \right) \quad (7) $$

where $\epsilon_r^{sta}$ is the yield strain at the failure in quasi-static and $\epsilon_r^{dyn}$ the one in dynamic. $\dot{\epsilon}_0$ represents the initial strain rate and the parameter $D_4$ can be determined from experiments. Then, we can obtain the critical damage in dynamic by the following equation:

$$ d_r^{dyn} = d_0 + (d_r^{sta} - d_0) \left( 1 + D_4 \ln \frac{\dot{\epsilon}}{\dot{\epsilon}_0} \right)^2 \quad (8) $$

where $d_0$ and $d_r^{sta}$ are the initial and the static critical damage respectively.

## Three-Point Bending Experiments

Three-point bending experiments of three fresh human humeri were conducted on a drop tower apparatus of 2.30 m height. A drop mass fixed to a carriage was allowed to fall freely while guided on vertical rails. The carriage height can be adjusted to adapt the impact velocity. A laser sensor M7L for accurate measurement against all surfaces was fixed to the drop tower bed in order to record the successive positions of the impactor. The two epiphyses of each humerus are fixed by an epoxy resin in two aluminum prismatic brackets (Figure 1) and the diaphysis was equipped with five strain gauges.
<figure><img src="image_3.png" /><figcaption>Figure 1. Three-point bending experiment.</figcaption></figure>

## Multi-Scale Finite Element Model

The finite element model of the humerus mechanical response under impact was elaborated from the digitalized STL format obtained by medical scans of the system humerus-brackets. First, the 3D CAD geometry of the diaphysis was generated using the CATIA © v6.0 software [19] and the inner and outer surfaces of the humerus diaphysis were reconstituted. Next, a first mapped finite element mesh was generated using only quads. Then, a 3D mesh was extruded using 7920 hexahedrons by means of the HYPERMESH © v2017 software [20].

The rigid brackets have been modeled by solid-shell elements [21] and the humerus epiphyses were embedded into the brackets using a prismatic solid having the resin material parameters. The system humerus-brackets was positioned on a plate which has been modeled by rigid elements. The impactor has been modeled using shell elements and an elastic material (see Figure 2). The simulation was carried out using the Explicit algorithm of LS-DYNA © R11 software [14].

<figure><img src="image_4.png" /><figcaption>Figure 2. Finite Element model</figcaption></figure>

## RESULTS

### Experimental Local Identification of the Model Parameters

The first step was the experimental identification of the local elastic and damage behaviors. As shown in the Figure 3, small specimens were dissected from the same humerus and tested mechanically. First, nanoindentation tests were performed on four samples (I1 to I4) of 8 × 5 × 2 mm size using a nanoindenter Nano XP from MTS factory equipped with a Berkovich tip in order to identify the elastic properties of the bone solid matrix. This nanoindenter has a load resolution of 10 nN, a displacement resolution of 0.01 nm and a maximum indentation load of 500 mN.

We obtain an average value of $E_s = 17.37 \pm 1.97$ GPa for the elastic modulus which is in agreement with the results available in the literature [22]. The value of the Poisson ratio $v_s$ was fixed to 0.3 as commonly proposed for long bones in the literature [23, 24].

<figure><img src="image_5.png" /><figcaption>Figure 3. Specimens localization on the humerus</figcaption></figure>
A second test was performed by incrementing the applied displacement at a constant velocity of 0.05 mm/min in order to evaluate, until failure, the specimens behavior as a function of the damage evolution. Cyclic three-points bending experiments on two specimens F₁ and F₂ of 58 × 5 × 3 mm size were performed using an Hounsfield H5KT machine directly connected to a computer recording data of load and displacement in real time. The H5KT testing machine is also equipped cell of 5 kN and a servomotor for displacement. The results showed an average value of 100 ± 38 N/m for the initial bending stiffness and of 92 ± 30 N/m for the final one.

Finally, cyclic compression tests were performed on six specimens (C₁ to C₆) of 2 × 3 × 2.8 mm using an ElectroPuls E3000 test machine from Instron company equipped with two compression plates. This machine designed for dynamic and static testing on a wide range of materials and components possesses a linear stroke of 60 mm and a dynamic load capacity of 3 kN. The compression stiffness was 9119 ± 783 N/mm in the longitudinal direction and 9568 ± 425 N/m in the transverse one.

## Numerical Local Identification of the Model Parameters

The next step was the numerical identification of the local parameters of the proposed micromechanical model by using the basic equations of the damage for bar or beam structures and exploiting the results of the previous cyclic compression and bending tests.

For the identification of the damage variables (k and ξ) under compression loading, a one dimension analytical model using bar theory in large strains was proposed. The one dimension hypothesis (Figure 4) leads to the following expression of the damaged elastic modulus as function of the elastic modulus of the solid matrix E_s and the damage variable d:

$$E = E_s(1 - d) \qquad (9)$$

By means of our micromechanical model and the results of cyclic compression experiments, we can obtain this explicit expression of the critical damage as a function of the critical strain:

$$d_r = d_0 + \frac{E_s}{2\xi} \epsilon_r^2 \qquad (10)$$

<figure><img src="image_6.png" /><figcaption>Figure 4. 1D model of humerus specimen under compression loading</figcaption></figure>

The evolution of the damage in this case with respect to the compression strain of the specimen C₁ is given in the Figure 5.
<figure><img src="image_7.png" /><figcaption>Figure 5. Damage evolution with respect to the strain for compression specimen $C_1$.</figcaption></figure>

This procedure was applied to the six specimens ($C_1$ to $C_6$). The identification results of compression damage parameters are summarized in Table 1.

<table><caption>Table 1. Summary of damage parameters (case 1: compression specimens).</caption><thead><tr><th>Specimen</th><th>ξ</th><th>ε<sub>r</sub></th></tr></thead><tbody><tr><td>C<sub>1</sub></td><td>8.8</td><td>0.022</td></tr><tr><td>C<sub>2</sub></td><td>9.1</td><td>0.022</td></tr><tr><td>C<sub>3</sub></td><td>8.3</td><td>0.020</td></tr><tr><td>C<sub>4</sub></td><td>8.3</td><td>0.020</td></tr><tr><td>C<sub>5</sub></td><td>8.1</td><td>0.023</td></tr><tr><td>C<sub>6</sub></td><td>7.1</td><td>0.023</td></tr><tr><td>Mean</td><td>8.28</td><td>0.0216</td></tr></tbody></table>

On the other hand, the identification of the damage parameters under three points bending loading is made in the same way as in compression case. Here, the exact beam bending theory is used.

The 1D hypothesis (Figure 6) makes it possible to arrive at an explicit relationship of the critical damage as a function of the curvature $\chi$:

$$d_r = d_0 + \frac{E_s}{2\xi} e^2 \chi^2 \frac{L_0}{2} \qquad (11)$$

<figure><img src="image_8.png" /><figcaption>Figure 6. 1D model of humerus specimen under bending</figcaption></figure>
Then, the evolution of the damage of the specimen F1 under bending loading is given in the Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7. Damage evolution with respect to the strain for bending specimen F₁</figcaption></figure>

This procedure was applied again to the second specimen F₂. The identification results of bending damage parameters are summarized in the Table 2.

<table><caption>Table 2. Summary of damage parameters (case 2: bending specimens).</caption><thead><tr><th>Specimen</th><th>ξ</th><th>ω<sub>r</sub></th><th>χ<sub>r</sub></th><th>ε<sub>r</sub></th><th>d<sub>r</sub></th></tr></thead><tbody><tr><td>F<sub>1</sub></td><td>9.5</td><td>2.1</td><td>0.0109</td><td>0.0202</td><td>0.2919</td></tr><tr><td>F<sub>2</sub></td><td>9.1</td><td>2.2</td><td>0.0109</td><td>0.0164</td><td>0.2747</td></tr><tr><td>Mean</td><td>9.3</td><td>2.15</td><td>0.0109</td><td>0.0183</td><td>0.2833</td></tr></tbody></table>

Finally, by means of the cyclic three point bending tests, the parameters of the Johnson-Cook model are determined and listed in the Table 3.

<table><caption>Table 3. Summary of the Johnson-Cook parameters (case 2: bending specimens).</caption><thead><tr><th>Specimen</th><th>d<sub>0</sub></th><th>ξ</th><th>ε<sup>sta</sup><sub>r</sub></th><th>D<sub>4</sub></th><th>d<sup>sta</sup><sub>r</sub></th><th>d<sup>dyn</sup><sub>r</sub></th></tr></thead><tbody><tr><td>F<sub>1</sub></td><td>0.029</td><td>9.5</td><td>0.0202</td><td>-0.3315</td><td>0.2919</td><td>0.0553</td></tr><tr><td>F<sub>2</sub></td><td>0.045</td><td>9.1</td><td>0.0164</td><td>-0.2959</td><td>0.2747</td><td>0.0799</td></tr><tr><td>Mean</td><td>0.037</td><td>9.3</td><td>0.0183</td><td>-0.3137</td><td>0.2833</td><td>0.0676</td></tr></tbody></table>

According to the results of the table 4, it can be seen that the values obtained for the two specimens F₁ et F₂ are very close. Thus, the found critical damage value at the failure was $d_r \approx 0.02$. This maximum damage value will be retained and used for the validation of the proposed model of the human humerus, as a criterion for eliminating elements that will reach this damage value.

## DISCUSSION

The aim of this section is to evaluate the predictive capabilities of the proposed micromechanical model and in particular its ability to reproduce the main features of mechanical behavior of the human humerus under dynamic loading.

The Figure 8(a) represents the impact load as a function of time for the three humerus 485, 514 and 582. As we can observe, the load curve of the FE elastic model without damage, performs increasingly fast oscillations over the period of 10 ms. Because the humerus bone is assumed to be elastic and undamageable (no degradation), only the first oscillation is to be retained in this study and compared with the physical measurements. Indeed, the rest of the curve has no physical meaning, because in practice the humerus breaks at a critical damage value. Moreover, it can be seen from the three
figures that the maximum impact load is correctly estimated by the FE elastic model without damage compared to the experimental results.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Mechanical response of humerus under impact: (a) Experiment vs numerical (b) Predicted fracture patterns.</figcaption></figure>

The introduction of damage in the FE model was made by means of the parameters $d_r^{dyn}$ and $\xi$. The filter SAE 1000 was, then, applied in order to smooth the dynamic oscillations and to obtain a mean curve of the FE model with damage.
This model had allowed to reach the first peak of the elastic curve and to estimate more precisely the maximum impact load with a minimum error of 7%.

As soon as the critical damage value $d_r^{dyn}$ is reached, the finite element is removed from the mesh by the software. The option $d_r^{dyn}$ was used in the software in order to obtain a first response on the humerus failure (Figure 8b). The first stalling of the FE curve with damage from the elastic curve occurs at time $t = 2$ms for the humerus 485 and 514 and at time $t = 3$ms for the humerus 582. However, the experimental response of the humerus is damped by the presence of the marrow, which is not taken into account in the FE models. This explains the differences between numerical models and experiments.

For the fracture patterns predicted by the proposed damage model, they are in agreement with the physical humerus fracture with degrees of concordance ranging from 70% to 90%. The fracture is pronounced at the bottom of the humerus slightly offset from the axis of the impactor for the humerus 485 and 514. For the humerus 582, which has a particular geometry and a greater length than the other humerus, the predicted fracture is initiated at the bottom at two symmetrical points on either side of the impactor (see the Figure 8b).

It should be noted that the fracture patterns depend highly on the mesh size according to the literature [25, 26]. Then, it is necessary to refine the mesh size in order to improve the breakthrough initiation points and to introduce the XFEM method in order to correctly predict the rupture propagation paths after initiation [27]. The proposed micromechanical model is limited only to estimating the fracture initiation points.

## CONCLUSIONS

A new theoretical formulation, development and validation of a ductile damage model applied to the human humerus bone under impact loading, has been proposed in the present investigation. The outcome of the proposed model appears to be more biofidelic than classical macroscopic-based models, since it predicts the observed physical response of several impacted human humerus with a good accuracy within the range of variability parameters inter-donors. It has to be noticed also, that the ultimate impact load that a human humerus may encounter before fracture, can be estimated with the proposed coupled model with a high accuracy.

The fracture patterns predicted by the proposed micromechanical damage model are consistent with the physical observed humerus ruptures even if this model is limited only to the fracture initiation. Further improvements will be performed to the present model by taking into account the marrow effects and fracture propagation by means of advanced Peridynamics or XFEM techniques.

## REFERENCES

[1] A. Linder and W. Svedberg, "Review of average sized male and female occupant models in European regulatory safety assessment tests and European laws: Gaps and bridging suggestions," *Accid. Anal. Prev.*, vol. 127, no. 1, pp. 156-162, 2019.

[2] S. Duma, P. Schreiber, J. McMaster, J. Crandall, C. Bass, and W. Pilkey, *Dynamic injury tolerances for long bones of the female upper extremity*. Gothenburg, Sweden: IRCOBI Conference, 1998.

[3] M. Van Eijnatten, R. Van Dijk, J. Dobbe, G. Streekstra, J. Koivisto, and J. Wolff, "Inv CT image segmentation methods for bone used in medical additive manufacturing," *Med. Eng. Phys.*, vol. 51, pp. 6-16, 2018.

[4] W.I.I Wan Iskandar Mirza, M.N Abdul Rani, M.A Yunus, R. Omar, and M.S Mohd Zin, "Alternative scheme for frequency response function measurement of experimental-analytical dynamic substructuring," *J. Mech. Eng. Sci.*, vol. 13, no. 2, pp. 4946-4957, 2019.

[5] C. Hellmich, J. Barthelemy, and L. Dormieux, "Mineral-collagen interactions in elasticity of bone ultrastructure - a continuum micromechanics approach," *Eur. J. Mech. A Solids*, vol. 23, pp. 783-810, 2004.

[6] O. Sam Daliri, M. Farahani, and M. Farhang, "A combined numerical and statistical analysis for prediction of critical buckling load of the cylindrical shell with rectangular cutout," *Eng. Solid Mech.*, vol. 7, pp. 35-46, 2019.

[7] O. Sam Daliri, and M. Farahani, "Characterization of Stress Concentration in Thin Cylindrical Shells with Rectangular Cutout Under Axial Pressure," *Int. J. Adv. Manuf. Tech.*, vol. 24, no. 3, pp. 133-141, 2017.

[8] F. Golesorkhie F, and M. Navi, "Effects of geometric variations on buckling properties of carbon nanostructures: A finite element analysis," *J. Mech. Eng. Sci.*, vol. 14, no. 1, pp. 6473-6487, 2020.

[9] A. Geromel Fischer, and J.J Marigo, "Gradient damage models applied to dynamic fragmentation of brittle materials," *Int. J. Frac.*, 2019, doi: 10.1007/s10704-019-00356-y.

[10] S. Lurie, Y. Solyaev, and K. Shramko, "Comparison between the Mori-Tanaka and generalized self-consistent methods in the framework of anti-plane strain inclusion problem in strain gradient elasticity," *Mech. Mater.*, vol. 122, pp. 133-144, 2018.

[11] J. Rahmoun, A. Auperrin, R. Delille, H. Naceur, and P. Drazetic, "Characterization and micromechanical modeling of the human cranial bone elastic properties," *Mech. Res. Commun.*, vol. 60, pp. 7-14, 2014.

[12] H. Naceur, J. Rahmoun, J. Halgrin, and P. Chabrand, "F Multiscale finite element modeling of ductile damage behavior of the human femur under dynamic loading," *Int. J. Damage Mech.*, vol. 24, no. 3, pp. 313-337, 2014.
[13] M. Murugesan, and D. W Jung, “Johnson Cook Material and Failure Model Parameters Estimation of AISI-1045 Medium Carbon Steel for Metal Forming Applications,” *Materials*, vol. 12, no. 4, pp. 609, 2014, https://doi.org/10.3390/ma120406093683.

[14] LSTC.LS-DYNA R11, “Edition: Livermore Software Technology Corporation I,” *Keyword User's manual*, 3186 pages, Livermore, California, 2018.

[15] Q.Z Zhu, J.F Shao, and D. Kondo, “A micromechanics-based thermodynamic formulation of isotropic damage with unilateral and friction effects,” *Eur. J. Mech. A Solids*, vol. 30, no. 3, pp. 316–325, 2011.

[16] R.A Hill, “consistent mechanics of composite materials,” *J. Mech. Phys. Solids*, vol. 13, pp. 213-222, 1965.

[17] J. Rahmoun, F. Chaari, E. Markiewicz, and P. Drazetic, “Micromechanical modeling of the anisotropy of elastic biological composites. Multiscale Modeling and Simulation,” *Multiscale Model Sim.*, vol. 8, no. 1, pp. 326–336 2009.

[18] A.V Sobolev M.V Radchenko, “Use of Johnson-Cook plasticity model for numerical simulations of the SNF shipping cask drop tests,” *J. Nucl. Sci. Technol.*, vol. 2, no. 4, pp. 272-276, 2016.

[19] Dassault Systemes, “CATIA V6 User’s Guide,” *3D CAD Design Software*, Vélizy-Villacoublay, France. 2013.

[20] Altair Engineering, “Hypermesh 2017 Reference Guide,” *Copyright 1986-2016 Altair Engineering Inc*, 4968 pages, 1820 East Big Beaver Road. Troy, MI 48083, United States. 2017.

[21] H. Naceur, S. Shiri, D. Coutellier, J.J Batoz, “On the modeling and design of composite multilayered structures using solid-shell finite element model,” *Finite Elem. Anal. Des.*, vol. 70, pp. 1–14, 2013.

[22] P.K Zysset, X. Edward Guo, C. Edward Hoffler, K.E Moore, and S.A Goldstein, “Elastic modulus and hardness of cortical and trabecular bone lamellae measured by nanoindentation in the human femur,” *J. Biomech.*, vol. 32, no. 10, pp. 1005-1012.

[23] V. Baca, Z. Horak, P. Mikulenka, and V. Dzupa, “Comparison of an inhomogeneous orthotropic and isotropic material models used for FE analyses,” *Med. Eng. Phys.*, vol. 30, pp. 924–930, 2008.

[24] L. Peng, J. Bai, X. Zeng, and Y. Zhou, “Comparison of isotropic and orthotropic material property assignments on femoral finite element models under two loading conditions,” *Med. Eng. Phys.*, vol. 28, no. 3, pp. 227-233, 2006.

[25] L.J Jia, W. Fujie, T. Ikai, S. Yochida, and H. Ge, “Jia L-J, Fujie W, Ikai T, Yochida S, Ge H,” *Eng. Fract. Mech.*, vol. 215, pp. 117-137, 2019.

[26] M. Korgesaar, and J. Romanoff, “Influence of mesh size, stress triaxiality and damage induced softening on ductile fracture of large-scale shell structures,” *Mar. Struct.*, vol. 38, pp. 1–17, 2014.

[27] S. Reddy Vempati, K. Brahma Raju, K. Venkata Subbaiah, “Simulation of Ti-6Al-4V cruciform welded joints subjected to fatigue load using XFEM,” *J. Mech. Eng. Sci*, vol. 13, no. 3, pp. 149–156, 2019.