Multiple scale modeling for cortical bone fracture in tension using X-FEM

Élisa Budyn∗,1 – and Thierry Hoc∗∗ ∗

Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, 842 West Taylor Street, Chicago IL. 60607, USA. ∗∗

Department of Material Science, LMSSMat, UMR 8579, Ecole Centrale Paris, Grande Voie des Vignes, Chatenay Malabry, F-92295.

We present a multiple scale method for modeling multiple crack growth in cortical bone under tension. The four phase composite Haversian microstructure is discretized by a Finite Element Method. The geometrical and mechanical bone parameters obtained by experiments mimic the heterogeneity of bone at the micro scale. The cracks are initiated at the micro scale where a critical elastic-damage strain driven criterion is met and are grown until complete failure in heterogeneous linear elastic media when a critical stress intensity factor criterion is reached. The cracks are modeled by the eXtended Finite Element Method. The simulations provide the global response at the macroscopic level and stress and strain fields at the microscopic level. The model emphasizes the importance of the microstructure on bone failure in assessing the fracture risk.

ABSTRACT.

RÉSUMÉ. Nous proposons une méthode multi-échelle pour modéliser la propagation de fissures multiples dans l’os cortical sous traction. La microstructure Haversienne composite à quatre phases est discrétisée par une méthode éléments finis. Les paramètres du tissu osseux géométriques et mécaniques obtenus par l’expérience décrivent l’hétérogénéité à l’échelle microscopique. Les fissures sont initiées à l’échelle micro où un critère élastic-endommagement piloté en déformation est atteint et sont propagées jusqu’à rupture complète dans des milieux linéaires élastique hétérogènes. Les fissures sont modélisées par la méthode des éléments finis enrichis (X-FEM). Les simulations montrent la réponse globale à l’échelle macroscopique et les champs de contraintes et de déformations à l’échelle microscopique. Le modèle montre l’importance de la microstructure dans le mécanisme de rupture en vue de déterminer un risque fracturaire. KEYWORDS:

finite element, X-FEM, cortical bone, fracture, composite, multi-scale method.

MOTS-CLÉS :

element finis, X-FEM, os cortical, rupture, composite, méthode multi-échelle.

REMN – 16/2007, X-FEM, pages 215 à 238

216

REMN – 16/2007, X-FEM

1. Introduction We present a multiple scale method for modeling multiple crack growth in cortical bone under tension. The model is decomposed into two parts : the modeling of random Haversian bone microstructures as a multi phase composite and the modeling of multiple crack growth inside this heterogeneous microstructure. Haversian cortical bone microstructure has been extensively studied as a factor influencing the mechanical behavior of the failure process of bone (Katz, 1976; Aoubiza et al., 1996; Gottesman et al., 1980; Hogan, 1992). At the micro scale, cortical bone is composed of densely packed concentric lamellar structures (osteons) that are embedded in an interstitial matrix (Figure 1). The cement line is an interface between the osteon and the interstitial matrix. Bone porosity is principally due to Haversian canals. Moreover, bone is a living tissue that can adapt both its structure and its architecture to its mechanical environment (Currey, 2003). This process produces a heterogeneous material with a high variability in osteon bone-mineral density that is related to the level of maturation of the newly-formed tissue. A mechanistic understanding of the role of this complex microstructure is strongly relevant to predicting the risk of fracture associated with age and disease. From a statistical point of view, several studies have demonstrated that cortical bone fracture toughness can be related to osteon density and porosity (Wright et al., 1977; Hui et al., 1988; Schaffler et al., 1988; O.Brien et al., 2005). At the microstructural scale, considerable work has been developed to study the initiation and growth of microcracks that have been shown to act as a stimulus for bone remodeling. In in vivo conditions, cement lines act as a microstructural barrier for the majority of cracks. In contrast, in transverse loading the underlying microstructure does not have any significant influence on the crack propagation (Nalla et al., 2003). Due to the complex microstructural hierarchy, the mechanical behavior and the fracture mechanisms have been very difficult to predict. Haversian cortical bone has often been compared to a composite where discontinuities within the material may generate mechanical stress concentration sites for crack initiation. Several models using homogenization theory have been developed to predict macroscopic behavior (Gottesman et al., 1980; Katz, 1971; Aoubiza et al., 1996), but they do not resolve the spatial distribution of local strain. Dong et al. (2005) developed a self-consistent model for multi-coated cylinders to obtain more precise bounds for cortical bone. Some attempts to model larger unit cells have been developed by Prendergast et al. (1996), however the geometrical and mechanical local anisotropy were not included. The objective of this work is to develop a model for Haversian cortical bone based on a multiscale approach in order to study the influence of the microstructure features on the strain localization pattern, which is strongly relevant in the mechanism of bone fracture. To solve this problem the eXtended Finite Element Method (X-FEM) which is a numerical method for treating arbitrary discontinuities without remeshing (Belytschko et al., 1999; Belytschko et al., 2001) is an efficient tool. In this method, cracks are

Multiple scale modeling for cortical bone

217

initiated at the micro scale at osteonal locations where a critical elastoplastic strain driven criterion is met. The cracks can therefore be randomly located within the microstructure. In the present paper, the Haversian microstructure that is a four phase composite: Haversian canal, osteons, cement line and matrix, is discretized by a finite element method. The geometrical and mechanical bone parameters obtained by nanoindentation measurements are randomly distributed to mimic the nature of bone at the micro scale. The discontinuities are represented by enriching the standard finite element approximation space. For cracks in the framework of linear elastic fracture mechanics, suitable enrichments are the Heaviside function for the crack and the Westergaard field for the crack tip. X-FEM is an application of the partition of unity (Melenk et al., 1996) and was introduced in (Belytschko et al., 1999; Dolbow et al., 2001; Moës et al., 1999) and later adapted for multiple cracks in Budyn et al. (2004). The crack topology is represented by vector level sets, a method that was developed by Osher et al. (1988) and Sethian (1999) for problems of interface tracking and later improved by Burchard et al. (2001) and Osher et al. (2002) for the evolution of curves. Here we use a narrow banded vector level set method developed by Ventura et al. (2002) and Ventura et al. (2003), because this method simplifies the process of freezing the existing level sets for cracks (Belytschko et al., 2001). The cracks are then grown in heterogeneous linear elastic media when a critical stress intensity factor criterion is met. The cracks with the maximum stress intensity factors are grown so that they remain approximately at the critical stress intensity factor by adjusting the load parameter and solving for the corresponding displacement. The solution method follows a "crack length control" algorithm. In the case of competing crack tips, a stability analysis is performed to determine the crack configuration path that leads to the maximum decrease in the potential energy. Stress intensity factors are computed by means of an interaction integral (Moran et al., 1987). The outline of this paper is as follows. In a first part, the description of the Finite Element Method to create the random Haversian bone microstructure is given (Section 2). The eXtended Finite Element Method for initiation and multiple crack growth problems until coalescence of cracks and percolation is presented (Section 3.3). Section 4 presents the results and compared with experimental results. Conclusions are given in Section 5.

2. Cortical bone microstructure We consider three dimensional unit cells of cortical bone of 0.5 millimeter width and 0.1 millimeter height. Cortical bone microstructures are represented by a four phase composite composed by an interstitial matrix and osteonal fibers (hollow discs in two dimensions, hollow tubes in three dimensions) coated by viscous cement lines and hollowed by Haversian canals considered as free boundaries where normally cytoplasmic fluid flows freely, see Figure1b.

218

REMN – 16/2007, X-FEM

Γt t = λ tο

Γcr1

2 l2 Γcr

l1

Ω Γcr3

y

l3

l4

Γcr4

x

B

A

(a)

Γu

(b)

Figure 1. (a) Light microscope bone observation. (b) Schematic model of the cortical unit cell

2.1. Geometrical description All geometrical parameters are experimentally measured from light microscope pictures of bovine bone samples as in Figure 1a. In our measurements the volumetric fraction of the interstitial matrix is about 34 percent and the average width of the cement lines is close to 5 µm. The distributions of the diameters of the osteons and the Haversian canals are statistically measured from the microscopic observations and regularized by a Gaussian fit (see Figure 2) denoted Ho for the osteons and Hc for the Haversian canals. The average diameters are 35 µm for the Haversian canals and 140 µm for the osteons. The random microstructure is constructed using a hard sphere scheme for which details can be found in Torquato (2000). In this scheme each osteon i is represented in two dimensions as a disc and defined by a position Pi (x, y) and a diameter Di . We assume the osteons are impenetrable. The construction process can be outlined as follows: the position Pi for one osteon is randomly chosen inside the unit cell. A diameter Di is randomly chosen within the Ho distribution. If the osteon volume is included in the unit cell without overlapping another osteon, it is accepted; otherwise, it is rejected. The procedure is repeated until the surface fraction of osteon covers more than 65 percent. Finally, a diameter for the Haversian canal for each osteon is randomly chosen within the Hc distribution in Figure 2. We check that the obtained porosity is between 5 to 9 percent according to our microscopic observations. The algorithm to construct the microstructure is programed in Fortran 90 and then en-

Multiple scale modeling for cortical bone

219

100 90 80

percent

70 60 50 40 30 experiment osteon regularization osteon experiment Havers regularization Havers

20 10 0 0

50

100 150 diameter (µ m)

200

250

Figure 2. Experimental and regularized distribution of the size of the osteon diameters and the Haversian canal diameters

coded in Python as a script file. This geometry is first discretized using a commercial three dimensional FEM code (ABAQUS, 2004) and linear tetrahedron elements for their efficiency and adaptability to complicated geometries. A preliminary local strain field analysis is performed to determined the locations where the cracks are initiating. The same geometry is then discretized in a two dimensional X-FEM code using Maltab (MATLAB7, 2006) and Gmsh (GMSH, 2006) and quadratic triangles. The longitudinal axis of the osteons is parallel to the z-direction.

2.2. Mechanical properties In this study the material properties are assigned in each phase based on experimental observations. The four phases of the interstitial matrix, the Haversian canals, the cement lines and the osteons can be described by four different constitutive laws. However, the bone remodeling process produces a heterogeneous material with a high variability in osteon bone-mineral densities and geometries.Therefore, the local mechanical properties of the osteons and the cement lines will be different in each osteon. The osteons are chosen to be transversally isotropic in the three dimensional FEM Abaqus model used for preliminary strain field analysis and isotropic in the two dimensional XFEM Matlab/Gmsh final model. Hoc et al. (2006) have performed nanoindentation measurements of the local Young’s modulus in wet conditions for over a hundred osteons in the longitudinal osteonal axis (Fan et al., 2002). It is shown that the local Young’s modulus and the bone-mineral content is reasonably correlated (r2 =0.75).We assumed the Young’s moduli to be directly related to the mineral content in our model. Humid nanoindentation measurement was also per-

220

REMN – 16/2007, X-FEM 100 90 80

percent

70

E exp. T ET regul. E exp. L E regul. L

60 50 40 30 20 10 0 5

10

15 E (GPa)

20

25

Figure 3. Experimental and regularized distributions of the Young’s moduli of the osteons measured by nanoindentation in the transversal plane and along the axis of the osteons

formed in the transversal direction. The obtained experimental distributions of the Young’s moduli were regularized by a Gaussian fit (see Figure 3) denoted HT in the transversal direction and HL in the longitudinal direction. Gaussian distributions for local longitudinal and transversal Poisson’s ratio HνL and HνT are constructed from micro-extensometry measurements on millimetric bone samples. The transverse shear modulus GT is described by a Gaussian distribution HG of mean value and standard deviation taken from Reilly et al. (1975). The attribution of the material properties for each osteon i can be outlined as follows: a probability pi is randomly chosen to which corresponds a longitudinal Young’s modulus within HL distribution, a transversal Young’s modulus within HT and the shear modulus within HG . For the Poisson’s ratio, the probability is fixed to 1 − pi in order to obtain the longitudinal and transversal Poisson’s ratio νLi and νTi . In our model, we imposed a larger Poisson’s ratio to osteons with lower Young’s moduli considering that less mineralized osteons contain more collagen (Katz et al., 1971). The critical stress intensity factors of the cement lines are linearly correlated to its Young’s modulus by a factor 10−4 to fit in the range of 0.7 to 2 MPa (Shelton et al., 2003; Hazenberg et al., 2006). Table 1 summarizes the material properties of the osteons. Lakes has shown cement lines present a specific chemical composition which gives an isotropic viscoelastic behavior (Lakes et al., 1979). This composition confers different elastic properties compared to the osteons that are encircled by the cement lines. The Young’s modulus is taken 25 percent lower than the Young’s modulus of the osteon it encircles. The Poisson’s ratio is chosen to be 0.49 (Sabelman et al., 1997 Jul 20-23). Although only the elastic properties of the cement lines are used in this study, The critical stress intensity factors of the cement lines are linearly correlated

Multiple scale modeling for cortical bone

221

Table 1. Microscopic elastic properties Osteons Mean value Standard deviation Extremum min Extremum max Matrix

ET (GPa) 11.51 1.84 7.79 15.55 12.66

EL (GPa) 20.71 3.10 14.33 24.71 22.78

νT 0.17 0.06 0 0.343 0.153

νL 0.18 0.06 0 0.353 0.162

GT (GPa) 3.57 0.25 2.82 4.3 3.93

to its Young’s modulus by a factor 10−4 and increased by 40 percent. The Haversian canals are modeled by free boundaries and assume no viscous contact between the cytoplasmic fluid and the canal walls. The interstitial matrix is considered homogeneous and transversally isotropic with Young’s moduli about 10 percent stronger than the mean values of the osteonal moduli. This assumption is consistent with nanoindentation tests performed by Rho et al. (2002) that reveal higher Young’s moduli in the matrix than in the osteons. The Poisson’s ratio are chosen 10 percent lower than the mean value of the osteonal Poisson’s ratio.

3. Multiple crack growth by X-FEM 3.1. Elastic-damage strain driven criterion to initiate cracks A linear elastic analysis using a FEM discretization with a commercial code Abaqus is applied to cortical bone microstructures. The transversal Young modulus E2 was determined for 4 different microstructure sizes (0.2, 0.3, 0.4 and 0.5 mm) and was found to be constant (11.5 GP a) for microstructures of size greater than 0.3 mm. Therefore we choose a microstructure of 0.5 mm height and width and 0.1 mm thickness that is statistically significant . As many biological materials, failure in bone can be described by an elastic-damage strain driven criterion. In this study we determine the initial crack locations in the microstructure by using the critical maximum principal strain. When a critical value of 0.4 % (Pattin et al., 1996; Bayraktar et al., 2004) has been reached in the elastic calculation, cracks are initialized for the XFEM calculation. Tension tests are performed on the microstructure (0.5 mm) given in Figure 4 in the y-direction to determine regions where the critical yield strain is reached. Note that the strain localization zone direction depends on the loading direction (Budyn et al., 2005). The FEM analysis shows that cracks can be initiated perpendicularly to the direction of maximum principal stress, mode I in tension. Under vertical tension loading, initial cracks are placed horizontally starting from the edge of Haversian

222

REMN – 16/2007, X-FEM

Figure 4. 0.5 mm unit cell under 0.3% strain tension on the top face where the local region over 0.4% are in light grey near the Haversian canals

canals. Note that we progressively increases the displacement at the top of the face until local zones with a strain above 0.4% appear; Figure 4 shows the local maximum principal strain fields when the top face of the bone cell is submitted to 0.3%.

3.2. Description of the problem For our X-FEM multiple crack growth problem, we consider a two dimensional elastic unit cell Ω with boundary Γ, nc cracks with surfaces Γcr = {Γα cr , α = 1 to nc }, and nt crack tips as shown in Figure 1b. The normal to a surface is denoted by n. The cracks are assumed to be traction free. Prescribed tractions t are imposed on the boundary Γt and prescribed displacements are imposed on Γu . According to 3.1 for the boundary conditions we can prescribe a uniform traction along the top edge of the cell and fix the displacement in the y-direction along the bottom edge and the displacement in the x-direction of node A in Figure 1b. The implementation is limited to linear elastic fracture mechanics in np multi phase unit cells. The equilibrium equation and boundary conditions are given by: ∇ · σ = 0 in Ω

[1]

σ · n = t on Γt

[2]

¯ on Γu u=u

[3]

¯ is the where σ is the stress, u(x) is the displacement field, uT = {ux , uy }; u prescribed displacement on Γu .

Multiple scale modeling for cortical bone

223

Each phase of material is considered elastic governed by C the tensor of elastic moduli, and in small deformation. The length of each crack is denoted ℓi ; the set of crack lengths is represented by a matrix ℓ = {ℓi }, i = 1 to nt . The imposed traction t depends linearly on a scalar parameter called the load factor λ: t = λto , where to is a reference traction field. For a preliminary study of multiple crack growth in cortical bone modeling, the crack propagation is described by linear elastic fracture mechanics for brittle multi phase body by the following Lagrangian form Budyn et al. (2004): Z nt X

L(ℓ, u) = W(ℓ, u) +

i=1

Γicr

Gic dℓi

[4]

where W(ℓ, u) is the potential energy of the system, Gic is the critical energy release rate at crack tip i. Gic is a constant parameter characteristic of each phase and a function of (x, y). If we define nact as the number of active crack tips, the second term in the right hand side of Equation [4] is the energy dissipated during the growth of the nact active crack tips. The potential energy W(ℓ, u) can be decomposed into the strain energy W int and the load potential W ext : W(ℓ, u) = W int (ℓ, u) − W ext (u) where W int (ℓ, u) =

1 2

Z

[5]

ǫ(ℓ, u) : C : ǫ(ℓ, u) dΩ

[6]

Z

[7]

Ω\Γcr

W ext (u) = λ

Γt

to · u dΓ

The equilibrium states of the body Ω correspond to the stationary points of [4] or points on the boundary of the feasible domain, so the solution u ∈ U corresponds to: δL = δu W(ℓ, u)δu +

nt  X ∂W(ℓ, u) i=1

∂ℓi

 + Gic δℓi > 0,

∀δu ∈ Uo , ∀δℓi > 0

[8] where δu W(ℓ, u) is the variation of W with respect to u, δu is the variation of the displacement and dℓi is a crack length differential and: U = {uu ∈ H1 (Ω \ Γcr ) , [[u · n ¯]] > 0 on Γα ¯ on Γu } cr , u = u

[9]

Uo = U ∩ {uu = 0 on Γu }

[10]

where n ¯ is the normal to the crack and H1 is the Hilbert space of functions with square integrable derivatives. Note that the displacement field is discontinuous across the crack and that the jump [[]] in the normal displacement is required to be non-negative.

224

REMN – 16/2007, X-FEM

The above gives: δu W(ℓ, u) = 0

∂W(ℓ, u) + Gic = −Gi + Gic > 0, ∀i ∈ {1, .., nt } ∂ℓi

[11] [12]

where the energy release rate is defined as: Gi = −

∂W(ℓ, u) ∂ℓi

[13]

and [11] is the equilibrium Equation and [12] is the Griffith criterion for each crack tip i. As it is often common to express the crack growth law in terms of stress intensity √ factors, the energy release rate at tip i can be rewritten Ki = E ′ Gi where E ′ is the effective Young’s modulus. At each tip i, we compute the equivalent stress intensity i factor Keq (Paris et al., 1965) expressed as the norm of the contribution in mode I and mode II of the interaction integrals (Yau et al., 1980; Belytschko et al., 1999) as a surface integral (Budyn et al., 2004): q
2
i i 2 Keq = KIi + KII [14] The crack growth law in linear elastic fracture mechanics of Equation [12] can be written:  i if 0 < Keq < Kci ∆ℓi = 0 ( no growth ) [15] i if Keq = Kci ∆ℓi ≥ 0 ( growth )

where ∆ℓi is the crack growth increment of tip i. We do not consider any closure of the cracks, though we note that ∆ℓi < 0 is not a valid solution. If the crack closes, the inequality [[u · n ¯ ]] > 0 must be enforced. The cracks are grown in the direction of the maximum hoop stress, θi :    q
2 1 i i i i i KI /KII ± KI /KII + 8 [16] θ = 2 arctan 4 The above gives two directions; we choose the angle that corresponds to the posii tive maximum hoop stress. The stress intensity factors KIi and KII are computed by an interaction integral. The equivalent stress intensity factor is computed by (14). The expressions for the J-integral and the interaction integral as domain integrals can be found in Moran et al. (1987). An explicit algorithm is used to satisfy both Equations [11] and [12]. The Equilibrium [11] is discretized and solved for a prescribed traction to and then the load parameter λ is adjusted to satisfy [12], so that the crack with the maximum stress intensity factor met its critical value. The crack growth increments are set at the beginning of each step based on a “crack length control” scheme so that the evolution is controlled by a monotonically increasing function of the total crack length ∆ℓtot , i.e. the sum of all active crack lengths (Bocca et al., 1991). In some cases, several cracks are close to their critical stress intensity factors, the number of nact active

Multiple scale modeling for cortical bone

225

crack to grow is determined by stability analysis and a set of competitive crack tips Ncomp is determined as in Budyn et al. (2004). A stability analysis is conducted by applying a second variation of the Lagrangian form [4] to determine the most stable crack configuration evolution corresponding to that with the minimum energy dissipation (Budyn et al., 2004). The derivatives of the energy release rate with respect to the crack length (Nguyen et al., 1985; Nguyen et al., 1990), are computed by the generalized X-FEM formulation (Budyn et al., 2004) of the FEM formulation developed by Suo et al. (1989), Suo et al. (1992) and Suo et al. (1998), All subdeterminants of this matrix [∂Gi /∂ℓj ] are computed at time tn−1 and the maximum subdeterminant gives the set of tips Nact that will grow at time step tn determined by:       2 ∂ (Gi (σo )) ∂ L(ℓ, u) = max det ≥0 Nact = i ∈ Ncomp  max det − ∂ℓi ∂ℓj ∂ℓj {i,j}∈Ncomp {i,j}∈Ncomp [17] The variation of Gic in Equation [17] vanishes as we considered strong discontinuities in bone material properties and Gic constant in each material phase. Note that when multiple cracks grow at the same time, the Griffith criterion [12] is exactly satisfied at the crack with the maximum stress intensity factor met its value, and approximately satisfied at the other cracks.

3.3. X-FEM discretization The eXtended Finite Element Method (X-FEM) is applied to approximate the displacement field and incorporate the crack representation with a step function enrichment for the discontinuity of the interior of a crack and the asymptotic near-tip displacement field enrichment at the tips of the cracks (Belytschko et al., 1999; Moës et al., 1999; Budyn et al., 2004). Figure 5 displays the enrichment scheme. In a mesh with a set of nodes I, all corner nodes of elements crossed by a crack will be enriched. J n contains the set of corner nodes of the elements cut by the crack n enriched by the step function of crack n (circled nodes in Figure 5). Km contains the set of corner nodes of the element that contains the crack tip m enriched by the branch function (squared nodes in Figure 5). Nc is the set of cracks and Nt the set of crack tips in the entire model. The crack geometry of crack n is described by narrow banded level set functions interpolated by signed distance function f n (x) through (Ventura et al., 2003). The displacement field is based on modified functions (Stazi et al., 2002) adapted for multiple crack problems (Budyn et al., 2004) is as follows: ! nc X nt X 4 X X X X h n ¯n m ¯m ˜ ˜ u (x) = NI (x)uI + NJ (x)aJ HJ (x)+ NK (x) blK FlK (x) I∈I

n=1J∈J n

m=1K∈Km

l=1

[18] ˜J (x) where NI (x) are the shape functions for the continuous displacement field; N are the shape functions applied to the enrichment field. This choice of shape functions

226

REMN – 16/2007, X-FEM

(a)

(b)

Figure 5. X-FEM representation of a crack. The circled nodes are enriched by the step function of crack n and the square nodes are enriched with the tip enrichment of tip m. The function f gives the distance of a point to the crack

for the enrichment field is explained in (Stazi et al., 2002; Chessa et al., 2003). aJn ¯ n (x) of crack n and are the additional unknowns for the modified step enrichment H J m blK are the additional unknowns for the tip enrichment of tip m for the modified lth m branch function F¯lK (x); and xJ is the position of node J. The modified enrichment n ¯ and F¯ m are given in Budyn et al. (2004). As the enrichment for each functions H J lK crack is local, the size of the global stiffness matrix of the system is not much larger than that of the standard stiffness (Budyn et al., 2004).

3.4. Crack tip reaching a free boundary or another crack Within each step, the length of each active crack tip is increased in the direction of the maximum hoop stress by: ∆ℓtot ∆ℓi = [19] nact where nact is the number of active cracks that grow at the beginning of step n; nact is determined by a stability analysis at the end of step n − 1 described in (17). Three possibilities can arise before percolation: a crack can reach the external free boundary, reach the free boundary of an Haverse canal or coaslesce with another crack. When this occurs, the crack tip is annihilate or “kill” the crack tip and its near enrichment removed and replaced by either a simple step function for free boundaries or a “junction” enrichment for the linked discontinuities derived in Budyn et al. (2004). The decision to join boundaries or an another cracks is based on heuristic considerations based on the idea that a very small amount of brittle material subjected to a

Multiple scale modeling for cortical bone

227

tensile stress will break and therefore join the boundary or another crack (Kienzler et al., 2002; Rubinstein, 1996; Demir et al., 2001; Deng et al., 1992; Nemat-Nasser et al., 1987). Note that a special enrichment is developed in this model when a crack joins another crack in the same element where the tip of the connected crack is located. In this special case of junction, the near tip enrichment of the connected crack should be removed as well and changed into a step enrichment to ensure a perfect strong discontinuity within this entire element that contains the junction. For each crack tip i that is identified to grow, we consider a virtual increment ∆ℓivirt that is either set, when approaching free boundaries, to ∆ℓivirt = ∆ℓtot , where ∆ℓtot is the total crack growth per step, or more restrictively when approaching cracks, to ∆ℓivirt = max(∆ℓi , ri ), where ∆ℓi is the increment of growth of [19] and ri is the radius of the domain of computation of the interaction integral at tip i; the radius ri being about twice the element length. When the mesh is fine, we have observed that when the elements of the annular domain on which the J-integral is computed are partially truncated by the internal (Haverse canals) or external boundary of the model or contains another crack or another phase, the stress intensity factors still gives the correct direction of propagation. When another phase encounters the domain of calculation of the J-integral, the correct direction of crack propagation is still obtained because the gradient between the moduli of the phases is not too large. However the restrictive approach criterion is considered when many cracks interact. The length of the added increment ∆ℓi is also adjusted when the growing crack encounters a free boundary or another crack so that the tip dies on the boundary edge.

4. Results We consider a square of cortical bone of width 0.5 mm containing 58 cracks initiated after the analysis of the local strain field shown in Figure 4. The init cell is loaded by tension in the x-direction; the force is prescribed at the top edges and the displacement in the y-direction is prescribed equal to zero at the bottom edge. The displacement of the bottom left corner node in the x-direction is prescribed equal to zero in Figure 1b. The material properties are random and representative of cortical bone microstructure as described previously and Figure 6 shows the distribution of the local transverse Young’s moduli. The force-deflection curves, expressed as a nominal stress versus nominal strain curve, are shown for a osteons/cement lines SIF’s contrast of about 6 and 3 respectively in Figure 7a and c. The evolution of the nominal stress at each time step are shown for a osteons/cement lines SIF’s contrast of about 6 and 3 respectively in Figure 7b and d.In the present study, analysis is static and we consider two cases when the average Stress intensity factor (SIF) of the cement lines are six times or three times tougher as the osteons. In the first case the load deflection curve (Figure 7a) shows three phases: from step 1 to 36 the cracks are growing inside the osteons, which produces a strain hardening, from step 36 to 68, one dominant crack is able to cross through the cement line of the osteon it originated

228

REMN – 16/2007, X-FEM

4

x 10

0.25

1.9 0.2 1.8 0.15

58

57 38 5037 4926

12

46 25 3445

33 22

1.7 0.1

4

1.6

0.05

1.5

0

1.4

−0.05

1.3 1.2

−0.1

1.1 −0.15

24 14

−0.2

0.9

23

−0.1

0

(a)

0.1

0.2

51

30

29

43

20

19

13 36

48 47 41 6

40 5 54

−0.2

44

52 1

9

1

−0.25 −0.3

3 2

10

42

21

11

8 39 16 53

0.3

35

7 28

15 32

5631

18 55

27 17

(b)

Figure 6. (a) Heterogeneous E2 Young’s moduli distribution (MPa). (b) Initial microstructure that is implemented to start the XFEM simulations; the locations of the initial 58 cracks are determined by the ABAQUS pre-calculation

Table 2. Case when the cement lines are three times tougher than the osteons. The active crack’s names are given in Figure 6b. Several events can happen: A: reaches the outer boundary, B: connects to crack 7, C: reaches an Haversian canal, D: reaches crack 35 but cannot connect 35 again, E: reaches an Haversian canal, F: reaches the outer boundary, G: reaches the outer boundary, H: reaches the outer boundary step 2-5 6-10 11 12 13-15 16-18 19 active crack 7 8 7 36 4 3 4 event step 21 22-25 26-29 30 31-41 42 43-44 active crack 35 36 8 36 8 7 8 event A B C step 47-49 50-51 52-53 54 55 active crack 48 42 10 51 58 event E F G H

20 3 45-46 7 D

200

200

150

150

nominal stress (MPa)

nominal stress (MPa)

Multiple scale modeling for cortical bone

100

50

0 0

2

4 6 8 nominal strain (%)

10

100

50

0 0

12

20

40

80

100

(b)

200

200

150

nominal stress (MPa)

nominal stress (MPa)

60 step

(a)

100

50

0 0

229

2

4 6 8 nominal strain (%)

(c)

10

12

150

100

50

0 0

20

40

60

80

100

step

(d)

Figure 7. (a) Load deflection curve (SIF’s contrast of 6), (b) evolution of the load versus time (SIF’s contrast of 6). (c) load deflection curve (SIF’s contrast of 3), (d) evolution of the load versus time (SIF’s contrast of 3). The ratio of cement line average SIF’s to osteonal average SIF’s is about 6 (upper curves, the work of separation is 1.48 × 10−3 µN.mm ) and about 3 (lower curves, the work of separation is 1.73 × 10−4 µN.mm)

from, and damage the matrix. During this period a strain softening is observed as the matrix looses progressively its structural integrity. Finally from 68 to 88 the cracks grow until complete percolation, but the microstructure had already lost most of its integrity. In the second case when the cement line are three times tougher than the osteons, The same patterns are observed in the global load deflection curve: from step 1 to 20 the cracks are growing inside the osteons, which produces a strain hardening, from step 20 to 44, one dominant crack crosses through the

230

REMN – 16/2007, X-FEM

(a)

(b)

(c)

(d)

Figure 8. (a) Initial crack paths, (b) Crack paths ate step 36 (hardening phase), (c) Crack paths at step 68 (softening phase), (d) Crack paths at step 88 (broken microstructure). The ratio of cement line average SIF’s to osteonal average SIF’s is about 6

cement line of the osteon it originated from, and damages the matrix, which results in a strain softening. Finally from 44 to 55 the cracks grow until complete percolation. We also note the effect of the cement lines over the strain field distribution in Figure 4b and the crack paths in Figure 8 and Table 2 and Figure 9d and Table 3. In Budyn et al. (2005), the structural function of the cement lines was studied in transversal compression by comparing the local strain field distributions in the same osteonal geometry in a microstructure including cement lines around the osteons and another microstructure without cement lines. It was noticed in the absence of cement lines the osteons tend to transfer more deformation to the matrix; with cement lines,

Multiple scale modeling for cortical bone

(a)

(b)

(c)

(d)

231

Figure 9. (a) Initial crack paths, (b) Crack paths ate step 20 (hardening phase), (c) Crack paths at step 44 (softening phase), (d) Crack paths at step 55 (broken microstructure). Then the ratio of cement line average SIF’s to osteonal average SIF’s is about 3

the osteons remains confined and exhibit a wider range of deformations. This phenomenon can partially explain the hardening phase observed in the global behavior due to primary growth of cracks inside as many osteons as possible before damaging the matrix. Despite very similar crack paths in Figures 8 and 9, by comparing the two load deflection curves in Figure 7a, we observe that more inhomogeneity makes the response richer (Belytschko, 2006). Note that the hardening is more developed when the cement lines are tougher and more energy is needed to break the microstructure: the work of separation is 1.48 × 10−3 µN.mm when the cement lines are six times tougher than the osteons and 1.73 × 10−4 µN.mm when the cement lines are three times tougher than the osteons. Therefore the cement lines appear as critical elements

232

REMN – 16/2007, X-FEM

Table 3. Case when the cement lines are six times tougher than the osteons. The active crack’s names are given in Figure 6b. Several events can happen: A: reaches the outer boundary, B: reaches the outer boundary, C: reaches the outer boundary, D: connects to crack 7, E: reaches an Haversian canal, F: reaches the outer boundary, G: reaches the outer boundary, H: reaches the outer boundary, I: reaches the outer boundary, J: connects to crack 6 step 2-5 6-10 11 12 13-15 16-18 19 20 active crack 7 8 7 36 4 3 4 3 event step 21 22 23 24-26 27-29 30-32 33 34 active crack 35 22 21 5 6 28 27 15 event A B step 35 36 37-39 40-51 52-53 54-55 56-57 58-62 active crack 16 32 7 8 6 5 6 8 event C step 63 64-66 67 68 69-70 71 72-73 74-76 active crack 7 36 7 10 7 42 10 21 event D E F G step 77-78 79-81 82 83-88 active crack 14 13 58 8 event H I J

in the protection of bone against fracture and relevant structural elements to study their evolution in osteoporosis. When the cement lines are less tough, other physical results changes: a snap-back is less pronounced, the softening phase is smoothened. Some experimental observation of the effect of the cement lines as osteonal barriers can be seen in Mohsin et al. (2006). However more experimental data are necessary to investigate the local fracture properties such as the stress intensity factors. Once the matrix is damaged, softening appears at the macroscopic level as the structural integrity of bone is lost.

5. Conclusion The results of this preliminary study present a multiple scale approach to describe multiple crack growth in cortical bone under tension. The Haversian microstructure is discretized by a statistical finite element model of cortical bone material to investigate the influence of the local geometrical and mechanical parameters of the microstructure on the mechanical behavior of bone material. Even if in our model majority of geometrical and material parameters were based on physical data, some limitations can be given. For the geometrical description the osteons are modeled by disjoint circu-

Multiple scale modeling for cortical bone

233

lar tubes, which is an idealized geometry and does not allow them to connect to each other through Volkman canals (Jones et al., 2004). We study bone at the osteonal level and do not include osteocyte lacunae (Qiu et al., 2005; Prendergast et al., 1996). We also modeled interstitial matrix as homogeneous material. In our study we focused our investigations on the effect of the toughness heterogeneities due to microstructure. Other heterogeneities such as the spatial distribution and shape of the osteons, cement line thickness, the initial defect orientation and distribution ... influence the crack paths and will be studied in further work. To implement crack growth the eXtended Finite Element Method is particularly suitable for multiple crack growth because remeshing is avoided. Higher order elements have been chosen, they are quadratic for the standard displacement field and linear for the enrichment to give an accurate crack solution. In contrast to boundary element methods, the method can easily handle the heterogeneities of the microstructure. The method is explicit, applied to static crack growth and is developed to satisfies exactly both the equilibrium [11] and the brittle crack growth law [12] at each time step of the load deflection curve. A stability analysis based on energy consideration enable to solve the case of competitive crack growth. The response of the unit cell is tracked until complete failure when the cracks have percolated the cell, joined and reached the free boundaries. Some improvement in precision in the load deflection curve will be obtained by using smaller steps and finer meshes to refine the heuristic criterion for crack bridging and crack connecting to the free boundaries. In the present study we considered traction loading applied to the top edge of the cell. We are currently investigating consistency between the average local stress field when a displacement is applied on the top edge and the average local strain field when the corresponding traction loading to the first displacement boundary case is applied to the top of cells containing defects. We have already studied these issues for cells of uniform materials containing defects and the global responses was extremely similar in both cases when load BC or corresponding displacement BC were imposed to the top of the cells. For heterogeneous cells of 0.5 mm without defect it is also possible to recover the same macro field by applying displacement BC or the traction BC corresponding to the first load BC case. The model accesses two different scales: a macro scale at the unit cell level and a micro scale for the strain field inside the osteons, in a mimetic osteonal microstructure. At macroscopic scale, the model is able to provide the overall stress and strain response. At the microscopic level, even if the matrix is homogeneous its strain field is not. Moreover, our results show that cement lines have an important role in isolating the osteons from the matrix. This explains why cement lines are implicated with energy in fracture processes in the deflection of crack propagation by slowing it down (Advani et al., 1987) or debonding the osteon (Guo et al., 1995). Cement lines have also been linked to energy absorption (Akkus et al., 2001) and viscous damping (Lakes, 2001; Dong et al., 2004) in allowing osteon movement inside the interstitial matrix.

234

REMN – 16/2007, X-FEM

In summary, this model presents preliminary results that characterize the local properties of cortical bone and emphasizes the importance of the cement lines as an efficient structural component to ensure internal movement at the material level and to prevent the progression of localization damage zones and cracks. Our results suggest how critical circular porosity (Haversian canals) are sources of localization and fracture nucleation for transversal tension loading. More validations of experiments will be developed for this type of numerical approach, where a dialog between numerical modeling and experimental investigation can be established. Further investigations of the physiological and mechanical properties of the cement lines and the matrix will also help to refine the modeling of damage and fracture in cortical bone. In the future, introduction in this model of more experimental damage and fracture parameters need to be determined to make the model more physiological and extended to compression loading. This will lead towards a multiscale study of bone fracture in an entire bone.

Acknowledgments The authors are grateful for the research support of the University of Illinois at Chicago. The first author is also grateful to the CNRS for a Post-doctoral fellowship. We are also grateful to Professor Ted Belytschko for his comments and suggestions.

6. References ABAQUS, User’s Manual, Hibbit, Version 6.3, Karlsson & Sorensen, Providence, 2004. Advani S. H., Lee T. S., Martin R. B., “ Analysis of crack arrest by cement lines in osteonal bone”, Proc. ASME Annual Winter Meeting, vol. 3, p. 57-581, 1987. Akkus O., Rimnac C. M., “ Cortical bone tissue resists fatigue fracture by deceleration and arrest of microcrack growth”, Journal of Biomechanics, vol. 34, p. 757-764, 2001. Aoubiza B., Crolet J. M., Meunier A., “ On the mechanical characterization of compact bone structure using the homogenization theory”, Journal of Biomechanics, vol. 29, n◦ 12, p. 1539-1547, 1996. Bayraktar H. H., Morgan E. F., Niebur G. L., Morris G. E., Wong E. K., Keaveny T. M., “ Comparison of the elastic and yield properties of human femoral trabecular and cortical bone tissue”, Journal of biomechanics, vol. 37, p. 27-35, 2004. Belytschko T., “ Private communications”, 2006. Belytschko T., Black T., “ Elastic Crack Growth in Finite Elements With Minimal Remeshing”, International Journal for Numerical Methods in Engineering, vol. 45, n◦ 5, p. 610-620, 1999. Belytschko T., Moës N., Usui S., Parimi C., “ Arbitrary discontinuities in finite elements”, International Journal for Numerical Methods in Engineering, vol. 50, n◦ 4, p. 993-1013, 2001. Bocca P., Carpinteri A., Valente S., “ Mixed mode fracture of concrete”, International Journal of Solids and Structures, vol. 27, p. 1139-1153, 1991.

Multiple scale modeling for cortical bone

235

Budyn E., Funfschilling C., Jeulin D., Meunier A., Hoc T., “ Effect of strain localization of cortical bone using a multiscale modeling”, Journal of Biomechanics, 2005. submitted. Budyn E., Zi G., Moës N., Belytschko T., “ A method for multiple crack growth in britle materials without remeshing”, International Journal for Numerical Methods in Engineering, vol. 61, n◦ 10, p. 1741-1770, 2004. Burchard P., Cheng L. T., Merriman B., Osher S., “ RMotion of curves in three spacial dimensions using a level set approach”, International Journal for Numerical Methods in Fluids, vol. 170, p. 720-741, 2001. Chessa J., Wang H. W., Belytschko T., “ On construction of blending elements for local partition of unity enriched finite elements”, International Journal for Numerical Methods in Engineering, vol. 57, p. 1015-1038, 2003. Currey J. D., “ The many adaptations of bone”, Journal of Biomechanics, vol. 36, p. 1485-1495, 2003. Demir I., Zbib H. M., Khaleel M., “ Microscopic analysis of crack propagation for multiple cracks, inclusions and voids”, Theoretical and applied Fracture Mechanics, vol. 36, p. 147164, 2001. Deng H., Nemat-Nasser S., “ Bridging Multi-scale Method for Localization Problems”, Mechanics of Materials, vol. 13, p. 15-36, 1992. Dolbow J., Moës N., Belytschko T., “ An extended finite element method for modeling crack growth with frictional contact”, Computer Methods in Applied Mechanics and Engineering, vol. 190, n◦ 51-52, p. 6825-6846, 2001. Dong X. N., Zhong X. N., Huang Y., Guo X. E., “ A generalized self-consisitent estimate for the effective elastic moduli for fiber-reinforced composite materials with multiple transversely isotropic inclusions”, International Journal of Mechanical Sciences, vol. 47, n◦ 6, p. 922940, 2005. Dong Z. N., Guo X. E., “ Geometric determinants to cement line debonding and osteonal lamellae failure in osteon pushout tests”, Journal of Biomechanical Engineering ASME, vol. 126, p. 387-390, 2004. Fan Z., Swadener E., Rho J. Y., Roy M. E., Pharr G. M., “ Anisotropic properties of human tibial cortical bone as measured by nanoindentation”, Journal of Orthopaedic Research, vol. 20, p. 806-810, 2002. GMSH, User’s Manual, GMSH, Version 62, Geuzaine & Remacle, Cleveland USA & Louvain Belgium, 2006. Gottesman G., Hashin Z., “ Analysis of viscoelastic behaviour of bones on the basis of microstructure”, Journal of Biomechanics, vol. 13, p. 89-96, 1980. Guo X. E., He M. Y., Goldstein S. A., “ Analysis of viscoelastic behaviour of bones on the basis of microstructure”, Bioengineering Conference - ASME, vol. BED-29, p. 112-117, 1995. Hazenberg J. G., Taylor D., Lee T. C., “ Mechanisms of short crack growth at constant stress in bone”, Biomaterials, vol. 27, p. 2114-2122, 2006. Hoc T., Henry L., Verdier M., Aubry D., Sedel L., Meunier A., “ Effect of microstructure on the mechanical properties of Haversian cortical bone”, Bone, vol. 38, p. 466-474, 2006. Hogan H. A., “ Micromechanics modeling of Haversian cortical bone properties”, Journal of Biomechanics, vol. 25, n◦ 5, p. 549-556, 1992.

236

REMN – 16/2007, X-FEM

Hui S. L., Slemenda C. M., Johnston C. C., “ Age and bone mass as predictors of fracture in a prospective study”, J. clin. Invest., vol. 81, p. 1804-1809, 1988. Jones A. C., Sheppard A. P., Sok R. M., Arns C. H., Limaye A., Averdunk H., Brandwood A., Sakellariou A., Senden T. J., Milthorpe B. K., Knackstedt M. A., “ Three-dimensional analysis of cortical bone structure using X-ray micro-computed tomography”, Physica A, vol. 339, p. 125-130, 2004. Katz J. L., “ Hard tissue as composite material - I. Bounds on the elastic behavior”, Journal of Biomechanics, vol. 4, n◦ 5, p. 455-473, 1971. Katz J. L., “ Hierarchical modeling of compact Haversian bone as a fiber reinforced material”, Advanced in Bioengineering - ASME, vol. , p. 17-18, 1976. Katz J. L., Ukraincik K., “ On the anisotropic elastic properties of hydroxyapatite”, Journal of Biomechanics, vol. 29, n◦ 4, p. 221-227, 1971. Kienzler R., Herrmann G., “ Fracture criteria based on local propertiesof the Eshelby tensor”, Mechanics Reasearch Communications, vol. 29, n◦ 6, p. 521-527, 2002. Lakes R., “ Viscoelastic Properties of Cortical Bone”, Bone mechanics, vol. 11, p. 1-15, 2001. Lakes R., Saha S., “ Cement line motion in bone”, Nature, vol. 204, p. 501-503, 1979. MATLAB7, User’s Manual, MATLAB, Version 7, The MathWorks, United States, 2006. Melenk J. M., Babuska I., “ The partition of unity finite element method: Basic theory and applications”, Computer Methods in Applied Mechanics and Engineering, vol. 39, p. 289314, 1996. Moës N., Dolbow J., Belytschko T., “ A Finite Element Method for Crack Growth without Remeshing”, International Journal for Numerical Methods in Engineering, vol. 46, n◦ 1, p. 131-150, 1999. Mohsin S., O’Brien F. J., Lee T. C., “ Osteonal crack barriers in ovine compact bone”, Journal of Anatomy, vol. 208, p. 81-89, 2006. Moran B., Shih C. F., “ A general treatment of crack tip coutour integrals”, International Journal of Fracture, vol. 35, p. 295-310, 1987. Nalla R. K., Kinney J. H., Ritchie R. O., “ Mechanistic fracture criteria for the failure of human cortical bone”, Nature Materials, vol. 2, p. 164-168, 2003. Nemat-Nasser S., Hori M., “ Toughening by partial or full bridging of cracks in ceramics and fiber reinforced composites”, Mechanics of Materials, vol. 6, p. 245-269, 1987. Nguyen Q. S., Stolz C., “ ENDOMAGEMENT, FATIGUE, RUPTUE. - Sur le problème en vitesse de propagation de fissure et de déplacement en rupture fragile ou ductile. Note de Nguyen Quoc Son et Claude Stolz, présentée par Paul Germain”, Comptes-Rendus de l’Académie des Sciences de Paris, vol. t. 301, Série II, n◦ 10, p. 661-664, 1985. Nguyen Q. S., Stolz C., Debruyne G., “ Energy methods in fracture mechanics: stability, bifurcation and second variations”, European Journal of Mechanics, A/Solids, vol. , n◦ 2, p. 157-173, 1990. O.Brien F. J., Taylor D., Lee T. C., “ The effect of bone microstructure on the initiation and growth of microcracks”, Journal of Orthopaedic Research, vol. 23, p. 475-480, 2005. Osher S., Cheng L. T., Kang M., Shim Y., Tsai Y. H., “ Geometricoptics in a phase-spacebased level set and Eulerian framework”, Journal of Computational Physics, vol. 179, n◦ 2, p. 622-648, 2002.

Multiple scale modeling for cortical bone

237

Osher S., Sethian J. A., “ Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations”, Journal of Computational Physics, vol. 79, n◦ 1, p. 1249, 1988. Paris P. C., Sih G. C., “ Stress analysis of cracks, In: Fracture toughness and testing and its applications”, Philadelphia: Society for Testing and Materials, vol. , p. 30-83, 1965. Pattin C. A., Calet W. E., Carter D. R., “ Cyclic mechanical property degradation during fatigue loading of cortical bone”, Journal of biomechanics, vol. 29, p. 69-79, 1996. Prendergast P. J., Huiskes R., “ Microdamage and osteocyte-lacuna strain in bone : a microstructural finite element analysis”, Journal of Biomechanical Engineering - Transactions of the ASME, vol. 118, p. 240-246, 1996. Qiu S., Rao D. S., Fyhrie D. P., Palnitkar S., Parfitt A. M., “ The morphological association between microcracks and osteocyte lacunae in human cortical bone”, Bone, 2005. in press. Reilly D. T., Burstein A. H., “ The elastic and ultimate properties of compact bone tissue”, Journal of Biomechanics, vol. 8, n◦ 6, p. 393-405, 1975. Rho J. Y., Zioupos P., Currey J. D., Pharr G. M., “ Microstructural elasticity and regional heterogeneity in human femoral bone of various ages examined by nano-indentation”, Journal of Biomechanics, vol. 35, p. 189-198, 2002. Rubinstein A., “ Macrocrack-Microdefect interaction”, Journal of Applied Mechanics, vol. 53, p. 505-510, 1996. Sabelman E. E., Koran P., Diep N., Lineaweaver W. C., “ Collagen/hyaluronic acid matrices for connective tissue repair”, First Smith & Nephew international Symposium: Advances in Tissues Engineering and Biomaterials, 1997 Jul 20-23. Schaffler M. B., Burr D. B., “ Stiffness of compact bone. Effect of porosity and density”, Journal of Biomechanics, vol. 21, p. 13-16., 1988. Sethian J. A., Level sets methods & fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision and materials science, Cambridge University Press, Cambridge, U.K., 1999. Shelton D. R., Martin R. B., Stover S. M., Gibeling J. C., “ Transverse fatigue crack propagation behavior in equine cortical bone”, Journal of Materials Science, vol. 38, p. 3501-3508, 2003. Stazi F. L., Budyn E., Chessa J., Belytschko T., “ An Extended Finite Element Method with Higher-Order Elements for Curved Cracks”, Computational Mechanics, vol. 31, p. 38-48, 2002. Suo X. Z., Combescure A., “ Sur une formulation mathématique de la dérivée seconde de l’énergie potentielle en théorie de la rupture fragile”, Comptes-Rendus de l’Académie des Sciences de Paris, vol. t. 308, Série II, p. 1119-1122, 1989. Suo X. Z., Combescure A., “ Double virtual crack extension method for crack growth stability assessment”, International Journal of Fracture, vol. 57, p. 127-150, 1992. Suo X. Z., Valeta M. P., “ Second variation of energy and an associated line independent integral in fracture mechanics. II Numerical validations”, European Journal of Mechanics, A/Solids, vol. 17, n◦ 4, p. 541-565, 1998. Torquato S., Random Heterogeneous Media - Microstructure and Macroscopic Properties Interdisciplinary Applied Mathematics - Mechanics and Materials, Springer-Verlag New York, New York Oxford, 2000.

238

REMN – 16/2007, X-FEM

Ventura G., Budyn E., Belytschko T., “ Meshfree and Particle Methods and their Applications”, International Journal for Numerical Methods in Engineering, vol. 58, p. 1571-1792, 2003. Ventura G., Xu J. X., Belytschko T., “ A vector level set method and new discontinuity approximations for crack growth by EFG”, International Journal For Numerical Methods in Engineering, vol. 54, p. 923-944, 2002. Wright T. M., Hayes W. C., “ Fracture mechanics parameters for compact bone - effects of density and specimen thickness”, Journal of Biomechanics, vol. 10, p. 419-430, 1977. Yau J., Wang S., Corten H., “ A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity”, Journal of Applied Mechanics, vol. 47, p. 335-341, 1980.

ANNEXE POUR LE SERVICE FABRICATION A FOURNIR PAR LES AUTEURS AVEC UN EXEMPLAIRE PAPIER DE LEUR ARTICLE ET LE COPYRIGHT SIGNE PAR COURRIER LE FICHIER PDF CORRESPONDANT SERA ENVOYE PAR E-MAIL

1. A RTICLE

POUR LA REVUE

:

REMN – 16/2007, X-FEM 2. AUTEURS : Élisa Budyn∗,1 – and Thierry Hoc∗∗ 3. T ITRE

DE L’ ARTICLE

:

Multiple scale modeling for cortical bone fracture in tension using XFEM 4. T ITRE

ABRÉGÉ POUR LE HAUT DE PAGE MOINS DE

40 SIGNES :

Multiple scale modeling for cortical bone 5. DATE DE

CETTE VERSION

:

15th February 2007 6. C OORDONNÉES

:

DES AUTEURS

– adresse postale : ∗ Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, 842 West Taylor Street, Chicago IL. 60607, USA. ∗∗

Department of Material Science, LMSSMat, UMR 8579, Ecole Centrale Paris, Grande Voie des Vignes, Chatenay Malabry, F-92295. – téléphone : (1) 312 996 96 31 – télécopie : (33) 01 41 13 16 16 – e-mail : [email protected], [email protected] 7. L OGICIEL

UTILISÉ POUR LA PRÉPARATION DE CET ARTICLE

:

LATEX, avec le fichier de style arti le-hermes. ls, version 1.23 du 17/11/2005. 8. F ORMULAIRE

DE COPYRIGHT

:

Retourner le formulaire de copyright signé par les auteurs, téléchargé sur :

http://www.revuesonline. om

S ERVICE ÉDITORIAL – H ERMES -L AVOISIER 14 rue de Provigny, F-94236 Cachan cedex Tél. : 01-47-40-67-67 E-mail : [email protected] Serveur web : http://www.revuesonline.com