Numerical Modelling of Crack Propagation in Cement PMMA: Comparison of Different Criteria

aMaterials and Reactive Systems Laboratory, Mechanical Engineering Department, Djillali Liabes University of Sidi Bel-Abbes, Larbi Ben Mhidi City, Algeria bMechanics and Physics of Materials Laboratory, Mechanical Engineering Department, Djillali Liabes University of Sidi Bel-Abbes, Larbi Ben Mhidi City, Algeria cMechanics Laboratory of Lille, Ecole Polytech’Lille, University of Lille, 59655, Villeneuve d’Ascq, France


Introduction
Total Hip Replacement (THR) is one of the most successful surgical procedures ever developed where a ballsocket structure is used to replace a diseased or damaged hip joint.The replacement cup socket is usually attached to the pelvis by acrylic bone cement, which consists of polymethylmethacrylate (PMMA) powder and a liquid component of methylmethacrylate monomer (MMA).When mixing together, polymerization takes place and within a few minutes (10-20 minutes depending on the formulation) of application to the bone cavity, the mixture becomes solid 1,2 .Figure 1 show a schematic of a cemented THR where the acetabular cup is fixed by the cement to the pelvic bone.The stability of the fixation critically depends on the integrity of the bone cement under typical physiological loading conditions.At least one million loading cycles per year would be experienced by the hip joint with the maximum hip contact force up to three times of body weight during walking 3 .
The presence of defects in the cement during mixing can locally lead to regions of stress concentrations producing a possible fracture of the cement and consequently loosening of the prosthetic cup.There are three major types of defects 4 : porosities, inclusions and cracks.
It is known that cracks are the most dangerous type of defect because of the presence of stress intensity on their front.Most cracks identified in orthopedic cement are 5 : • Cracks initiated at porosities (Figure 2).

•
Cracks initiated during cement withdraw.

•
Cracks initiated at the junction between bone and cement or between cement and cup.
In the literature, there has been a little research carried out into the crack's growth path in the orthopedic cement, there have been some studies dealing with the fatigue life and/or fracture of orthopedic cement, but not following the crack's growth path 6,7,8,9 .Benbarek et al. 10 used the finite element method to analyze the propagation path of the crack in orthopedic cement around a total hip replacement.Results show that the crack propagation's path can be influenced by human body posture.Benouis et al. 11 presented a numerical modeling of crack propagation trajectory in cement of reconstructed acetabulum.The direction crack is evaluated as a function of the displacement extrapolation technique and the strain energy density theory.Kim et al. 12 determine the specific fracture mechanics response of cracks that initiate at the stem-cement interface and propagate into the cement mantle.
For a crack under the mixed-mode I/II loading conditions, a number of fracture criteria have been developed through a concerted effort by many researchers in the past decades.Erdogan et al. 13 proposed the maximum circumferential stress (MCS) criterion, which assumed that fracture occurs in the direction where the circumferential stress surrounding the crack tip is the maximum.Smith et al. 14 developed a generalized maximum tensile stress criterion by taking T-stress into consideration.Alternatively, the minimum strain energy density (MSED) criterion introduced by Sih 15 assumed that fracture occurs in the direction where the strain energy density is the minimum.Hussain et al. 16 , Palaniswamy and Knauss 17 , Nuismer 18 and Wu 19 proposed the maximum energy release rate (MERR) criterion based on Griffith's theory 20,21 proposed a new  general mixed-mode brittle fracture criterion based on the concept of maximum potential energy release rate (MPERR).
This paper presents a finite element analysis for the modeling of the crack growth problems in cement of reconstructed acetabulum using the displacement extrapolation technique (DET).This modeling is based on the maximum tangential stress criterion (MTS), the minimum strain energy density (MSED) criterion and the general fracture criterion based on the approximate expression of energy release rate G(θ), to present a comparison of the crack propagation paths.In this investigation, the effect of the inclusions and cavities on the crack trajectories in cement orthopedic was examined.

Crack Growth Criteria
In order to simulate crack propagation under linear elastic condition, the crack path direction must be determined.There are several methods used to predict the direction of crack trajectory such as the maximum tangential stress theory (or the maximum circumferential stress theory) 13 , the minimum strain energy density theory 15 and the maximum energy release rate criterion based on Griffith's theory 20 .

MTS criterion
Erdogan et al. 13 were the first to propose a crack initiation criterion using stress as a critical parameter.This criterion states that direction of crack initiation coincides with the direction of the maximum tangential stress along a constant radius around the crack tip so it is called the maximum tangential stress (MTS) criterion.It can be stated mathematically as 26 : Where: K I and K II are respectively the SIFs corresponding to mode I and mode II loading.Eq. ( 2) can be solved for θ such that θ= θ 0 , which is the predicted crack initiation angle.

MSED criterion
Sih 27 used strain as a critical parameter in order to propose the minimum strain energy density (MSED) criterion.It states that the direction of crack initiation coincides with the direction of minimum strain energy density along a constant radius around the crack tip.The MSED criterion showed a good agreement with the experimental results obtained earlier by Erdogan et al. 13 .In addition, this criterion is the only one that shows the dependence of the initiation angle on material property represented by Poisson's ratio ν.In mathematical form, this approach can be stated as: Where S is the strain energy density factor, defined as: Where is the strain energy density function per unit volume, and r 0 is a finite distance from the point of failure initiation.The angle of crack propagation θ can be determined by solving the following equations 26 : 1

G(θ) criterion
In order to predict the fracture direction of cracked materials under the general mixed-mode state, Chang et al. 21resented a new general mixed-mode brittle fracture criterion based on the concept of maximum potential energy release rate (MPERR).This criterion can be easily degraded to the pure mode fracture criterion, and can also be reduced to the commonly accepted fracture criteria for the mixed-mode I/ II state.The criterion can be described as:

Numerical Calculation of Stress Intensity Factors
Fracture mechanics is based on the determination of stress intensity factors.It is therefore important to develop a numerical model capable of calculating these factors for different geometries of cracked structures under different boundary conditions.In this paper, the displacement extrapolation method 28,29 is used to calculate the stress intensity factors K I and K II as follows: , ( )  In order to obtain a better approximation of the field near the crack tip, special quarter point finite elements proposed by Barsoum 30 are used where the mid-side node of the element in the crack tip is moved to 1/4 of the length of the element L.

Methodology of Crack Modeling
In this paper, the APDL code has been employed to create the program to simulate the mixed-mode crack propagation in the cement.The displacement extrapolation technique is used, to determine the SIFs at each increment of the crack length.
The crack propagation is characterized by successive propagation steps performed without user interaction.Each step consists of: 1. Setting the geometrical with initial crack and input material properties data of the problem.10.Return to Step 3: The algorithm is repeated until ultimate failure of material or by using another criterion for termination of the simulation process.Figure 4 shows the Flow-chart of the prepared APDL code based on the combination of the finite element analysis and the crack growth criteria.a circular surface that was mated with congruent spherical acetabular socket.The mechanical properties of cement, cup and all subregions of the acetabulum bone are given in Table 1.

Geometric, Material's Definition and Loading Conditions
The model was generated from a roentgenogram of a 4mm slice normal to the acetabulum through the pubic symphysis and ilium.The inner diameter of the UHMWPE cup was 54mm and the cement thickness was taken as 2mm 24,25 .Figure 5a shows the geometrical model used in this investigation.The initial crack of length of a= 50μm emanating from the interface cup/cement is presumed to exist in the cement layer (Figure 5b).The model was divided into seven regions (Figure 5a) of different elastic constants with isotropic material properties assumed in each region.The femoral head was modelled as In order to approach the reality, force F=2400N is applied to the centre of the femoral head (Figure 6).It corresponds to a three times of average weight of man (80kg) 5 .These conditions of loading are considered as the extreme case of an effort which the prosthesis can support.It is the force which man can exert on only one support at time of rise of staircase 33 .The loading cases analysed in this study correspond to 0° of the metallic implant position.Figure 6 shows boundary conditions acting on the model.To simulate the connections with the remainder of the basin, we blocked the nodes corresponding to the junction with the iliac bone 34 .The boundary conditions imposed on the geometrical model are taken from previous work 10,11,32,35,36 .
u x and u y indicate the displacements in the x and y directions, respectively.

Finite Element Model
The Finite element standard code ANSYS 37 has been employed for modeling the problem.A typical FE model is shown in Figure 7a.The special quarter point singular elements proposed by Barsoum 30 are investigated for modeling the singular field near the crack tip in cement layer (Figure 7b).For a=50μm, the mesh discretization consists of 66583 elements and 136529 nodes.The numerical simulation is performed in plane stress conditions.

Evolution of SIFs
In order to simulate the behavior of a crack under mixed mode loading (mode I+II); we considered an example of a crack of length a localized in the cement layer.This crack is inclined to the horizontal x-axis from -35° to 90°.
Figures 8a and 8b show, respectively, the variations of the SIFs K I and K II according to the crack inclination α, for various crack length.The results obtained show that: 1.The SIF K I is maximum when the angle α=40° (Figure 8a), then it decreases gradually with the increase of the positive angle and tends to a very low value as one approaches the interface cup/cement.2. The SIF K I decreases with the angle α, up to a minimal value corresponds to an angle α=-18° from which factor K I takes values negative with the decrease of the angle α.These negative values of K I indicate that the crack lips are closed due to stress compression around the crack tip. 3. The SIF K II is null when the inclination angle α=40° and increases with this angle (Figure 8b).
The factor K II decreases with the angle α up to a minimal value corresponds to an angle α≈-18° from which the curve takes on a downward look with the decrease of angle α.Through these results, we can conclude that the initial crack angle of 40° represents the initial direction of the crack propagation according to the opening mode (mode-I) with: K I =K Imax and K II =0.
Figure 9 shows the evolution of SIFs K I and K II during crack propagation extension obtained by MTS criterion.These results are compared with those obtained by MSED and G(θ) criterion.The plotted curves show a good agreement between these three approaches.
A good correlation between these criterion are verified for another initial crack direction (α=0° and 60°).The crack paths obtained are shows in Figure 10.

Crack propagation simulation
The values of SIFs are used to determine in Figure 11, the final crack path by calculating the crack propagation direction at each time step.The three crack trajectories are obtained for initial crack propagation α=0°.This comparison shows a good correlation between the three simulations.
In order to determine the effect of a geometry defect on the crack propagation, we chose a circular cavity of radius r = 100 μm located in cement layer, at a vertical distance y = 0.25mm from the crack tip.Foucat 34 found that the size of the cavities existing in the cement is between a few micromillimeters to 1 mm.The geometrical model is meshed by triangular elements in cement layer and by special quarter point singular elements around the crack-tip (Figure 12).
Figure 13 presents three steps of crack propagation path using the GF criterion.We find that the crack is moving towards the cavity.This is because the cavity creates a stress drop that will change the maximum principal stress in the cement layer and draw the crack.Once the cavity has passed, this crack shifts under mode-I loading, and it departs slightly from the cavity.The comments we get correspond with those obtained by Benouis et al. 11 and by Boulenouar et al. 36 to study the crack growth path in the cement mantle using MSED criterion.
Figure 14a and 14b show a comparison between the crack trajectories obtained in cement layer with and without  a default, respectively.The crack trajectories obtained by MTS and GF criteria are compared with the results obtained by Benouis et al. 11 , using the MSED approach.The three criteria give good results on the crack propagation path and the results between them are very close.This behavior of crack propagation confirms the comments obtained in the case of a cracked path with a hole 38,39,40 .Without a defect, the crack changes the direction of propagation with an angle of 40°, and is propagated according to the opening mode (mode-I).In what follows, we propose to studied the crack propagation in our model containing an inclusion of radius r = 100 μm. Figure 15 shows the typical mesh model near the inclusion and the crack.We recall that the mechanical properties of the cement are: Young's modulus E 1 =2300 MPa and Poisson's ratio ν 1 =0.3 and the inclusion considered is characterized by its Young's modulus E 2 and Poisson's ratio ν 2 =ν 1 .
To demonstrate the inclusion effect on the propagation path, we evaluated for three ratios E 2 /E 1 = 1, 0.1 and 10, the SIFs and the crack direction for each increment Δa to predict then the propagation path using MTS approach.Figure 16 show respectively the steps of crack trajectory for three ratios E 2 /E 1 = 1, 0.1 and 10.A calculation made by our integrated model shows that: a) When the matrix (cement layer) and inclusion have the same mechanical properties (E 2 /E 1 =1), the crack is propagated in its possible direction as a single homogeneous material (Figure 16  a).b) In the case of the inclusion less rigid than the matrix (E 2 /E 1 =0.1), we find a result similar to that which had been obtained for the cracked cement with a cavity (see Figure 14b); this inclusion changes the stress distribution and attracts the crack.Once the inclusion exceeded the crack in its possible direction (Figure 16b).c) In the case of the inclusion is more rigid than the cement (E 2 /E 1 =10), the crack is this time slightly delayed (Figure 16c).Figure 17 shows the effect of inclusion on the final crack trajectories obtained by MTS and GF criteria (with E 2 /E 1 = 1, 0.1 and 10).Practically, the selected criteria give the same behavior of crack propagation explained for Fig. 16, using criterion MTS.
For better presenting the comparison between the crack path obtained by the three criteria, Figure 18 illustrates the final crack trajectories predicted by MTS and GF criteria, for each ratio E 2 /E 1 .Under the same loading conditions, these results are compared with the results obtained in reference 11,36 , using the MSED theory.This comparison gives a good correlation between the three calculations obtained for each ratio.

Conclusion
In this paper, three crack kinking criteria and the crack paths prediction for several cases are compared in a PMMA cement layer.The maximum tangential stress criterion, the strain energy density criterion and the general fracture criterion are investigated using advanced finite element techniques.In order to obtain a better approximation of the field near the crack tip in cement of reconstructed acetabulum, special quarter point finite elements proposed by Barsoum are used.The displacement extrapolation technique is investigated to determine the SIFs under mixed-mode loading.Numerical calculations made by the finite element method show that this technique can correctly describe the stress and deformations field near the crack tip.
The initial crack angle of 40° represents the initial direction of the crack propagation according to the opening mode.
The three criteria give good results on the crack propagation path and the results between them are very close, for the cases proposed in this study.
The implementation of 2D crack propagation process in FE code can account for the influence of inclusions on the path of crack propagation in cement layer.This feature can be very interesting for propagation in multilayered or in composite parts.

Figure 1 :
Figure 1: A schematic of a cemented acetabular replacement in a total hip replacement 22,23 .
G(θ) is the energy release rate in θ direction, can be expressed as 21 : Where k=3-4n for plane strain and k= (3-n)/(1+n) for plane stress.L is the length of the element side connected to the crack tip.u i and v i (i=b, c, d and e) are the nodal displacements at nodes b, c, d and e in the x and y directions, respectively (Figure3).

Figure 3 :
Figure 3: Special elements used for displacement extrapolation method.

Figure 4 :
Figure 4: A flowchart of the main operations which make the crack propagation

Figure 7 :
Figure 7: (a) Mesh model and (b) Special elements used for DET.

Figure 8 :
Figure 8: Variation of SIFs vs. crack length a and crack inclinaison α: a) K I evolution and b) K II evolution.

Figure 9 :
Figure 9: Variation of SIFs during crack extension for α=40°: a) K I evolution and b) K II evolution.

Figure 10 :
Figure 10: Variation of SIF K I and K II during crack extension for a) α=0° and b) α=60.

Figure 12 :
Figure 12: Typical mesh model near the cavity and the crack.

Figure 13 :
Figure 13: Three steps of crack propagation predicted by GF criterion (Cavity effect).

Figure 14 :
Figure 14: Crack trajectories comparison: a) without defect and b) with defect.

Figure 15 :
Figure 15: Typical mesh model near the inclusion and the crack.

2 .
Discretization of the model by plane2 elements.3. Mesh generation; refining around the crack-tip.4. FE analysis.5. Calculation of SIFs.6. Calculation of crack propagation direction using the criteria proposed in this investigation.7. Crack arrests?If yes, go to step 8.If no, go to step 9. 8. Stop.9. Automatic delete of the crack segment and insertion of the new crack segment for new crack tip position.