Energy release rate of small cracks under finite multiaxial straining

The energy release rate of small cracks governs fatigue crack nucleation. A method is presented here to efficiently and accurately evaluate the energy release rate of such cracks, arbitrarily oriented, under general conditions of finite multiaxial loading. As a motivation, the dependence on crack length is then investigated. It is demonstrated that the energy release rate of small cracks is proportional to the crack length and that the proportionality factor is a function of the far-field parameters only. An attempt is then made to search for a general expression of this proportionality factor under simple loading conditions.


INTRODUCTION
Studies carried out in the recent years have brought to light the physical phenomena governing fatigue life of rubber (Cam et al. 2004;Le Gorju 2007). It turns out to be driven by the growth of small cavities transforming then into small cracks and propagating throughout the material up to a critical size that leads to a major loss of material properties. The propagation of these small cracks represents the main stage of fatigue life.
Regarding small crack growth prediction two distinct approaches are usually considered. The crack initiation approach is based on the evaluation of the mechanical fields of a crack-free material in order to study how a small flaw would propagate when subjected to these mechanical conditions; this approach lies thus within the framework of continuum mechanics. Some predictors have already been developed to predict crack initiation (Mars 2002;Verron and Andriyana 2008). To the contrary, the crack propagation approach studies the propagation of an existing small crack embedded in the material (Gent et al. 1964;Lake and Lindley 1965).
Actually, in some simple cases it is possible to reconcile these two approaches. Indeed, considering a small crack of length c under plane stress uniaxial tension, Rivlin and Thomas have been able to factorize the energy release rate T (Rivlin and Thomas 1953): (1) where W stands for the strain energy density of the far-field region and k is a factor that depends only on the far-field loading conditions. This result has later been extended to both pure shear and equibiaxial tension (Yeoh 2002). Thus, it turns out that for a small crack subjected to simple loading conditions, the energy release rate is proportional to the crack length and the proportionality factor is a function of the far-field parameters only. In this paper we aim at developping a method to easily estimate the energy relase rate of small cracks under arbitrary plane stress loading conditions and then at investigating the factorization of the energy release rate with respect to the crack length.

A simple tool to model a small crack under arbitrary loading conditions
A method is needed first to model a small crack, arbitrarily oriented, under general conditions of far-field multiaxial loading.

Present definition of a small crack
In order to define the concept of small crack we introduce first the notion of "boundary" as a circular region drawn around the crack, on which we compare solutions for the "body" and the "crack neighbourhood" (see Figure 1): • "body" refers to the solution obtained for a deformed body without considering the effects of a crack. As we move away form the edge, the gradients induced by edge effects get smaller and imply an upper limit on the size of the "boundary". • "crack neighbourhood" refers to the solution obtained for a crack embedded in an infinite medium, with specified far-field loading conditions. As we move away from the crack, the gradients induced by the crack get smaller and imply a lower limit on the size of the "boundary".
A crack is then said to be small when one can draw a boundary around it such that the solution on the boundary, in both the "body" and in the "crack neighbourhood", is constant to within a given tolerance.
Note here that this definition is very general because edge effect is considered. However, in the present study, the body is made up of an isotropic hyperelastic material and is subjected to multiaxial homogeneous loading conditions; thus gradients of the mechanical fields in the body are nul. The body-crack neighbourhood boundary is then chosen at a location where gradients tend to zero.

Modelling boundary displacements under
plane stress multiaxial straining The small crack assumption permits us to impose the far-field state of stress directly via the displacement at the boundary. Indeed, except in a very small region around the crack, the stress and strain fields are homogeneous and are the same as in the crackfree model for which the analytical expressions of strain and stress fields along with the strain energy density are completely known.
The original model consists of a small throughcrack with orientation θ (in the undeformed configuration) under a far-field equibiaxial straining (λ 1 , λ 2 ) in the frame of reference ( , ) e e 1 2 (see Figure 2). Throughout the deformation the crack rotates and ends up with an orientation ψ with respect to its orientation in the undeformed configuration. The problem with this method is that a new geometry is required whenever the crack orientation θ is to be changed. However, from the perspective of the crack, changing the crack orientation boils down to changing the far-field loading (or equivalently the far-field straining) (see Figure 3). Thus, by performing first a change of basis and then a pull-back to maintain the crack orientation fixed throughout the deformation, we can express the associated far-field strain in the crack-based frame and capture the same range of conditions with a single FE model. Moreover, no restriction is made on the geometry of the far-field boundary which requires working out a very general expression of the displacement to be imposed to a node at the boundary. In order to conveniently study the effect of biaxiality, far-field loading parameters (λ 1 , λ 2 ) are replaced by (λ, B) with λ 2 = λ and λ 1 = λ B . Thus B = log λ 2 /log λ 1 quantifies the biaxiality and can be called the biaxiality  factor (Mars 2002). The relationship between the undeformed coordinates (X, Y) and the deformed coordinates (X′, Y′) of a point P in the crack-based frame ( , ) ′ ′ e e 1 2 (see Figure 2) under the far-field loading conditions described above is given by: and It should be pointed out that the implementation of the previous calculation (Eq. 2) and further computations carried out in this study were all obtained with the finite element software Abaqus in which the displacement of a node at the far-field boundary was imposed via a DISP subroutine.

Basic validations
First we have validated the calculations described in 2.1.2 by verifying that a crack-free model gives homogeneous loading for all B, λ and θ.
We then focused our attention on the "small crack" requirement. Indeed, far from the crack, it is desired to get a homogeneous stress field. This condition can never be met rigorously because of the finite dimension of the model and the presence of the crack. However, on using the method presented in section 2.1.1 a criterion was set to check accuracy of the approximation. The latter was considered acceptable when the relative error on the far-field stress when compared to the crack free model was less than 1%. From the experience a crack length 40 times smaller than the dimension of the model fits well into this criterion.
Note here that an additional satisfactory way to validate the calculations is to compare the energy release rate computed by Abaqus (J-integral) with the fracture mechanics solutions and from some studies carried out in finite strain (Yeoh 2002).

Energy release rates of small cracks and their factorization
It is first desired to develop a method that enables us to capture the influence of the presence of a small crack on the variation of the mechanical fields in the crack neighbourhood and that reveals how that variation makes up for the energy release rate. Using this method, the factorization of the latter with respect to the crack length is then investigated. Finally, on the basis of the numerical results obtained via finite element analysis, we examine the dependence of the energy release rate of small cracks on the farfield parameters B, λ and θ.

J-integral
The J-integral represents a way to calculate the energy release rate. It was first introduced within the framework of planar small strain as a contour path integral around the crack tip (Rice 1968). Thereafter, it was extended to planar finite strain. Indeed, on considering a contour Γ surrounding the crack tip and leaning on both faces of the crack in the undeformed configuration, the J-integral writes: where ∑ stands for the Eshelby stress tensor, q → is the crack direction vector and N → is the outward normal vector to the surface element dS in the under formed configuration.
We remind the reader that, for an elastic material, the J-integral is path-independent which permits the arbitrary choice of the contour surrounding the crack tip to calculate the energy release rate (Rice 1968). For the sake of simplicity, the far-field Eshelby stress tensor is denoted ∑ ∞ throughout the rest of this paper.
Let us focus now our attention on the contour. A rectangular contour of characteristic dimension R has been chosen (see Figure 4). On choosing R such that R/c → ∞ the evaluation gets highly simplified (see Table 1). Indeed, Γ A , Γ B and Γ C lie then in the far-field region wherein the Eshelby stress tensor is uniform and is equal to ∑ ∞ .
Using the symmetry with respect to the crack line, we deduce that: Note here that the expression above is not an approximation but the exact formula of the energy release rate. However, because the model has finite dimension, the upper bound for the integral is necessarily < + ∞. Thus, we can only compute an approximation of the energy release rate: where L is the maximum distance to the crack face in the ′ e 2 -direction.

Proportionality of J with respect
to the crack length Let us consider two different small cracks: one with length cand the other one with length kc where k is a real number. Both are subjected to the same farfield loading conditions. Using Eq. (7), the energy release rate of a crack of length cwrites: where ∑ 11 c (respectively ∑ 11 kc ) denotes the Eshelby stress associated with the crack of length c (respectively of length kc).
By using the substitution l = kl′, we obtain that dl = kdl′. It follows that: Let us focus now on the effect of multiplying the crack length by a factor k on the transformation of the mechanical fields with respect to the original problem of a crack of length c. We recall that the latter consists in a small crack embedded in a medium that has infinite dimension; thus changing the crack length by a factor k boils down to performing a homothetic transformation with a scale factor k and where the homothetic center is the center of the crack. And as the segment Γ + goes by the center of the crack, the mechanical fields along its length obey the homothetic transformation with respect to the crack length. Hence we have: and we deduce that: = kJ c ( ) Note here that no limits were placed on the particular plane stress loading condition. Thus, the argument we have just made is completely general and is valid for all λ, B and θ.
This result implies that it is sufficient to work out the energy release rates for one crack size, and Hence we have: that the energy release rate of all other crack sizes are then determined.

NUMERICAL RESULTS
The behaviour of the isotropic hyperelastic material was modeled by an incompressible neo-Hookean model.

Variation of the energy release rate versus crack length
In order to validate the proportionality law demon strated previously (see Eq. (15)), the dependence of the energy release rate on the crack length has been investigated for various far-field loading conditions (B ranging from -0.5 to 1, λ ranging from 1 to 5 and θ ranging from 0° to 90°). Indeed, for each and every far-field state of stress, the energy release rate has been computed twice for two different crack lengths: c and 2c. All the results obtained so far show that, for Bincreasing from -0.5 to 1 with an increment of 0.25, for λ increasing from 1.1 to 5.0 with an increment of 0.1 and θ increasing from 0° to 90° with an increment of 15°, we always verify:

Variation of the energy release rate versus crack orientation
For various biaxiality factors (B ranging from -0.5 to 1), the influence of λ on the variation of the energy release rate versus the crack orientation has been investigated. As we can observe (see Figure 5 and Figure 6), for B in [−0.5, 1] the energy release rate strictly decreases as the crack orientation increases in both small and finite strain: the maximum is reached at θ = 0° and the minimum is reached at θ = 90°. However, under equibiaxial tension (B = 1), the energy release rate is independent from the crack orientation in both small and finite strain and thus only depends on λ.

DISCUSSION: FAR-FIELD PARAMETER OF THE SCALE LAW
The proportionality of the energy release rate with respect to the crack length (section 2.3.2) permits the following factorization: where f is a function of the far-field parameters only.

Comparison with linear elastic fracture mechanics (LEFM)
From the solution of LEFM for a rubberlike material (Mars 2006), f can be factorized into: where and f 0 is a function of λ and B only. f 0 (λ, B) is basically the value of f(λ, B, θ) at θ = 0°. Thus, regarding the original problem (see Figure 2), f 0 depends only on the far-field state of stress. Note here that all the energy release rates computed in small strain for B ranging from -0.5 to 1 and θ ranging from 0° to 90° completely match this theoretical factorization.

Factorization of the energy release rate in finite strain under simple loading cases
Under equibiaxial tension, because the energy release rate is independent from the crack orientation, the previous factorization of f remains obviously valid. Under uniaxial tension, we can draw the same conclusion from the numerical results. Indeed, on comparing the factorization with the energy release rate values computed for λ ranging from 1.1 to 5, we always have: f (λ, −0.5, θ) = f (λ, −0.5, 0°) cos 2 θ (20) However, under pure shear and actually for all B in [−0.5, 1.0], a simple comparison with numerical results clearly shows that the previous factorization (Eq. 18) cannot be extended to finite strain.
As for now, all the arguments that we made in this subsection are based only on finite element results but they have the benefit to guide us towards the simplification of the expression of f under simple loading cases. From the remarks above we deduce that: J (λ, −0.5, θ, c) = (f 0 (λ, −0.5) cos 2 θ) c and J(λ, 1, θ, c) = f 0 (λ, 1) c wherein the general expression of f 0 still needs to be determined.

CONCLUSIONS
A new method for modelling small cracks under arbitrary loading states and finite straining has been presented. The latter was successfully tested against the most known cases (uniaxial tension, equibiaxial tension and pure shear). We have also investigated how the energy release rate is balanced by the distribution of configurational stresses. It has then been proved that the energy release rate of a small crack (suitably defined) always follows a linear scale law with respect to crack size, regardless of loading state.
Finally some progress has been made towards a general-purpose expression for energy release rate under arbitrary loading.