# Fatigue crack propagation: In situ visualization using X-ray microtomography and 3D simulation using the extended finite element method

Emilie Ferrié, Jean-Yves Buffière, Wolfgang Ludwig, Anthony Gravouil, Lyndon Edwards

► To cite this version:

Emilie Ferrié, Jean-Yves Buffière, Wolfgang Ludwig, Anthony Gravouil, Lyndon Edwards. Fatigue crack propagation: In situ visualization using X-ray microtomography and 3D simulation using the extended finite element method. Acta Materialia, 2006, 54 (4), pp.1111-1122. 10.1016/j.actamat.2005.10.053. hal-00436215

HAL Id: hal-00436215

https://hal.science/hal-00436215v1

Submitted on 5 Apr 2023

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 CC BY 4.0 - Attribution - International License
Fatigue crack propagation: In situ visualization using
X-ray microtomography and 3D simulation using the
extended finite element method

Emilie Ferrié<sup>a</sup>, Jean-Yves Buffière<sup>a,*</sup>, Wolfgang Ludwig<sup>a</sup>,
Anthony Gravouil<sup>b</sup>, Lyndon Edwards<sup>c</sup>

<sup>a</sup> Groupe d'Etude de Métallurgie Physique et de Physique des Matériaux, INSA de Lyon, 25 av. Jean Capelle, 69621 Villeurbanne Cedex, France

<sup>b</sup> Laboratoire de Mécanique des Contact et des Solides, INSA de Lyon, 69621 Villeurbanne, France

<sup>c</sup> Department of Materials Engineering, The Open University, Milton Keynes, UK

The propagation of a semi-elliptical crack in the bulk of an ultrafine-grained Al–Li alloy has been investigated using synchrotron radiation X-ray microtomography. In this material, the studied crack, despite its small dimension, can be considered as “microstructurally long” and described in the frame of the linear elastic fracture mechanics. The extended finite element method is used to calculate the stress intensity factors along the crack front taking into account the three-dimensional geometry extracted from the tomographic images. For the same nominal value of the stress intensity factor range, crack propagation is faster in the bulk than at the surface. The observed anisotropy is attributed to the variation of the closure stress along the crack front between surface and bulk. The experimentally observed fatigue crack propagation is compared to numerical simulations. Good agreement is found when a linear variation of closure stress along the crack front is taken into account in the “3D crack propagation law” used for the simulation.

Keywords: Fatigue crack propagation; Extended finite element method; Synchrotron radiation; Microtomography; Crack closure

# 1. Introduction

In the last 20 years, efforts have been made to develop
experimental techniques for investigating the behaviour
of fatigue cracks *within the bulk* of opaque materials. For
instance, indirect techniques such as direct (or alternating)
current potential drop measurement [1–3], Rayleigh wave
resonance [4] or compliance measurements [5] have all been
used to estimate the size of fatigue cracks in bulk samples.
However, no information is obtained on the three-dimen-
sional (3D) shape of the crack inside these samples and
the estimates given by these techniques are sensitive to sev-
eral factors such as the crack/probe distance [2] or the pres-

ence of residual stress [3]. More quantitative information
on crack front shape evolution can be obtained by means
of post-mortem fractographic techniques such as heat tint-
ing [6], the “marker technique” [7] or the more widely used
“beach marking technique” [8]. This latter technique con-
sists of increasing, during a certain number of fatigue
cycles, the minimum cycling stress while the value of the
maximum stress is kept constant. By doing so, a distin-
guishable striation-free area is created on the fracture sur-
face which allows the crack shape to be later determined.
However, this technique assumes that changing the stress
amplitude during the load cycle, without changing the
maximum load value, does not alter the crack shape. The
validity of this assumption remains questionable. Further-
more, in Al alloys, fracture surface markings obtained close
to the sample surface are blurred so that the crack shape

<footnote>* Corresponding author. Tel.: +33 4 72 73 88 54; fax: +33 4 72 43 85 39.
E-mail address: Jean-Yves.Buffière@insa-lyon.fr (J.-Y. Buffière).</footnote>
can only be determined using optical profile measurement
techniques [9]. These experimental methods have mainly
been used to study the shape evolution of part-through
cracks, which are, for most cases, more representative of
cracks in real structures than through cracks. Fracto-
graphic observations of samples subjected to mode I (or
mode II) loading have shown that the shape of part-
through cracks is usually near to semi-elliptical [9,10].

Alongside these characterisation methods, analytical
models predicting 3D fatigue crack growth of part-through
cracks in specimens with simple geometries (usually plates)
have also been developed, for example in Ref. [11]. They
are normally based on the overarching assumption that
the crack has a semi-elliptical shape so that crack extension
may be calculated at the surface and at the depth using
analytical stress intensity factor equations established for
tensile or bending loading [12] or polynomial stress distri-
butions [13]. Three-dimensional finite element (FE) calcula-
tions have also been used intensively to simulate crack
propagation in more complex crack/specimen/loading con-
figurations. Load spectrum effects have been investigated
[14] and propagation of non-planar cracks under mixed-
mode conditions has been predicted [15,16].

However, direct comparisons between the crack front
shape observed experimentally and that simulated by the
FE method are rather limited: final crack shape [15,16] or
beach marks [17] observed on the fracture surface have
been compared to 3D FE results and discrepancies were
found at the surface. Crack sizes obtained by 3D FE sim-
ulations for different configurations have also been com-
pared to experimental results [18] and an error of 30%
was found in the predicted crack shape evolution. Which-
ever method is used for crack propagation simulation, ana-
lytical or numerical, the local crack extension at the surface
and also in the bulk is calculated using Paris law equations
based on fatigue crack growth measurements at the speci-
men surface. This means that simulations of the crack front
shape evolution, and, indirectly, the prediction of the fati-
gue lifetime of the component or structure, are based on
the assumption that the crack growth behaviour measured
at the surface can account for crack growth behaviour in
the bulk. This assumption is necessary due to the lack of
reliable quantitative experimental data describing long
crack propagation in the bulk of structural materials.

The aim of the study presented here is to observe exper-
imentally long crack propagation in the bulk of metallic
alloys and to predict the crack shape evolution using linear
elastic fracture mechanics (LEFM). Synchrotron radiation
X-ray microtomography is used to investigate in three
dimensions the development of a semi-elliptical crack in
an ultrafine-grained Al alloy. This high-resolution 3D
imaging technique allows the 3D structure of objects in
the bulk of a non-transparent material to be visualised with
a resolution close to 1 µm. At such a resolution, the maxi-
mum volume that can be imaged is 1 mm³. Due to this size
limitation, the observation of crack propagation in
conventional metallic alloys, with a grain size in the range

of 10–100 µm, is limited to the short crack regime. Micro-
structurally short cracks have already been investigated in
coarse-grained cast Al alloys [19] and in cast iron [20].
Short cracks are strongly influenced by microstructural fea-
tures such as grain boundaries and exhibit multiple branch-
ing so that no modelling of the 3D propagation has been
attempted yet.

In the fine-grained alloy studied here, the crack size
(which is in the range 100–500 µm) is at least 100 times
the grain size, and can hence be considered as microstruc-
turally long as defined by Miller [21]. The extended FE
method (X-FEM) is used to calculate the stress intensity
factors along the crack front taking into account the real
3D geometry of the crack extracted from the tomographic
images. Based on experimental observations and quantita-
tive analysis of crack propagation in the bulk, a “3D crack
growth law” is established to predict the crack front shape
evolution in the bulk of the material.

## 2. Experimental

### 2.1. Material

The material used in this study is a powder metallurgy
Al–Li alloy based on AA5091 alloy obtained from Aero-
space Metal Composites Ltd, UK. The alloy composition
was Al–4.0Mg–1.2Li–1.0C–0.5O and was produced from
powders using a mechanical alloying process, followed by
hot isostatic pressing (hipping) [22]. The billet was supplied
as an upset forged 42 mm thick plate in the as-forged (T1)
condition. The material exhibited a yield stress of 450 MPa
and a tensile strength of 505 MPa, values similar to those
previously reported for 5091 in the T1 condition [22–24].
The tensile performance of Al–Mg–Li–C alloys is derived
from significant Hall–Petch strengthening, due to their very
fine-grain structure, solution and dispersion hardening.
They typically possess a high Mg concentration in compar-
ison to heat treatable Al–Li alloys, but the Li concentration
is deliberately kept below 1.5%, to prevent formation of the
age hardening δ' phase (Al₃Li). During the mechanical
alloying, surface layers of oxides on the powder particles
are repeatedly grown, then fractured and dispersed within
the particles. The grains, which are typically smaller than
1 µm, contain significant distributions of fine (50–100 nm)
precipitates, which are believed to be Al₄C₃, Al₂O₃ and
MgO. The Al₂O₃ and MgO phases derive from fragmenta-
tion of the surface oxides from the as-received materials.
The Al₄C₃ phase, by contrast, derives only from the reac-
tion between the aluminium and the process control agent
used to aid the sintering process [25].

The main advantage of using this alloy for the present
study is that it exhibits an exceptionally linear crack path
compared to ingot metallurgy Al alloys. The ultrafine
grains promote homogeneous deformation and prevent
crystallographic cracking so that the crack shape is not dis-
turbed by microstructural features, at least not at the level
of the spatial resolution employed in this study. This fine
microstructure also endows the alloy with low levels of
crack closure [23,25] and yield stress-insensitive fatigue
crack propagation behaviour so that the properties of the
as-forged T1 materials [23,27] are indistinguishable from
those of the fully heat treated condition [26].

## 2.2. X-ray microtomography: a 3D imaging technique

In situ fatigue tests monitored by X-ray microtomogra-
phy were carried out on beamline ID 19¹ of the European
synchrotron radiation facility (ESRF) in Grenoble, France.
This experimental station is dedicated to high-resolution
X-ray imaging and features a highly coherent X-ray beam,
a precision mechanics sample stage and a high-resolution
detector system. The latter consists of a 2024² CCD camera
coupled via microscope optics to a scintillator that
transforms X-ray photons into visible light. In the work pre-
sented here, a multilayer crystal is used to monochromatise
the incoming white beam to an energy of 20.5 keV. This
value allows imaging of Al samples of millimetric dimen-
sions. During a tomographic scan, the sample is rotated over
180° in order to acquire 1500 two-dimensional (2D) projec-
tion radiographs. From this set of projections, a quantita-
tive 3D map of the attenuation coefficient distribution
within the sample is produced by means of a standard fil-
tered backprojection tomographic reconstruction algo-
rithm. The spatial resolution obtained in the reconstructed
images is of the order of 1 µm, a value comparable to the res-
olution of an optical microscope. Due to the coherence of
the X-ray beam, the tomographic images contain, in addi-
tion to conventional “absorption contrast”, a “phase con-
trast” contribution due to propagation and Fresnel
diffraction of the highly coherent wavefield over the finite
specimen/detector distance. This phase contrast results in
an edge enhancement and is particularly useful for visualis-
ing small details for which the “absorption contrast” is very
low [28]. In our case, it improves crack-tip detection [19].
However, a too strong phase contrast will tend to spoil the
resolution and lead to reconstruction artefacts. For this rea-
son, the phase contrast contribution, is controlled by adjust-
ing the specimen/detector distance so that the diffraction
fringes only appear locally at the edges of heterogeneities.

## 2.3. *Fatigue machine and specimen geometry*

The sample was imaged in a dedicated fatigue machine
designed at INSA Lyon to perform in situ fatigue test at
the ESRF. Its specific features are the following: the load
frame consists of a poly(methyl methacrylate) tube nearly
transparent to X-rays such that the sample, under load,
is visible over 180° during the scan acquisition. The
mechanical design of the cyclic tension loading mechanism
enables operation at cycling frequencies of up to 80 Hz,
minimising thereby the cycling time required for a fatigue

test. This is a crucial point as, for a given experiment, the
average beamline availability is 72 h. The compact design
of the device and the low vibration level enabled us to
directly mount it on the sample stage of the tomograph
of beamline ID19. (A more detailed description can be
found in Ref. [29].)

At the spatial resolution employed (0.7 µm pixel size),
and in parallel beam configuration, the maximum gauge vol-
ume is limited by the detector area and amounts to approx-
imately 1.4 mm³. Small fatigue specimens with a minimum
cross-section of 1 mm × 1 mm were spark cut from a piece
of as-forged 5091 material. The geometry of the sample used
is represented in Fig. 1(a). A thin (2 µm) rectangular notch,
100 µm wide and 20 µm deep, was produced in the sample
using focused ion beam machining and is represented in
Fig. 1(b). This notch is located at the centre of one of the
specimen faces and acts as a crack initiation site.

## 2.4. *Fatigue experiment*

An in situ fatigue test was performed in air, at constant
stress amplitude with σmax = 220 MPa, a stress ratio
R = 0.1, a frequency of 40 Hz and at room temperature.
The value of σmax was chosen such that the stress intensity
factor calculated at the surface of a semi-circular crack
with a 100 µm radius is 2.86 MPa m<fnref n="1" />/² which is close to
the long crack threshold value of 2.7 MPa m<fnref n="1" />/² given in
Ref. [26]. Crack initiation was monitored by simple 2D
radiographs acquired at regular intervals. Once a crack
was detected, tomographic scans were recorded regularly
(every 1000 or 2000 cycles) to monitor crack propagation
in the bulk of the material. The specimen was imaged under

<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Geometry of the sample used for the in situ fatigue test monitored by X-ray microtomography. (b) 3D illustration of the notch obtained by focused ion beam machining and located at the centre of one of the specimen faces.</figcaption></figure>

<footnote n="1">¹ Beamline ID19 website: http://www.esrf.fr/exp_facilities/ID19/homepage/id19 homepage.html.</footnote>
load such that the crack was fully open on the 3D images. After 18,000 fatigue cycles, a 100% tensile overload was applied. After this overload, the load level was set to its initial value ($\sigma_{max} = 220$ MPa with $R = 0.1$) and the fatigue experiment was continued until failure.

## 2.5. Processing of the tomographic images

Each tomographic 3D image acquired during the fatigue test shows the crack shape “inside” the aluminium alloy. Amira 3D imaging software was used to visualise the 3D raw images and also to perform the segmentation operation required to extract the crack from the surrounding microstructure [30]. Reconstruction artefacts due to phase contrast produce spurious images so that a manual segmentation operation had to be performed. Once extracted, the crack was defined as an independent object and its full 3D geometry was visualised on 3D renditions. As an illustration, the 3D rendition of the crack at 16,000 cycles is displayed in Fig. 2. In this figure, the two orthogonal images that intersect the crack are slices extracted from the reconstructed 3D tomographic images. The rectangular focused ion beam notch is indicated in the $(X-Y)$ plane of the crack by black dotted lines.

## 2.6. Calculation of the local stress intensity factor values along the crack front

In this section, we present briefly the method used to calculate mixed-mode stress intensity factors $(K_I, K_{II}$ and $K_{III})$ along the crack front obtained from the 3D tomographic images. The numerical technique presented here couples the X-FEM [31] with the level sets method [32], which is used to define the crack geometry.

### 2.6.1. Crack modelling by enrichment functions

In FE methods, the presence of cracks in a structure must be taken into account in mesh generation: the mesh

<figure><img src="image_3.png" /><figcaption>Fig. 2. 3D rendition of the crack after 16,000 cycles: the rectangle in black dotted lines represents the notch in the $(X-Y)$ plane of the crack. The two orthogonal grey-scale images are slices along the $(X-Z)$ and $(Y-Z)$ planes extracted from the raw tomographic dataset.</figcaption></figure>

must adequately define the crack geometry. Special elements and considerable mesh refinement near the crack tip are necessary to capture accurately the asymptotic displacement fields. The X-FEM alleviates the shortcomings associated with meshing the crack surface by representing the crack geometry using additional functions called enrichment functions. The enrichment method, as described in Ref. [33], can be summarised as follows. The displacement field calculated in the structure “without crack” is locally modified by adding, at specific nodes, the nodal values of the enrichment functions. These functions define the crack geometry by “modelling” the discontinuity introduced by the crack in the displacement field. As the discontinuities at the crack front and along the crack surface are different, two enrichment functions, named $H$ and $F_J$, are necessary to model the entire crack.

The function $H$, defined as a generalised Heaviside function, models the displacement field along the crack surface. It is used to enrich the nodes of the element that are cut by the crack surface and can be defined as follows:

$$H(x) = \begin{cases} +1 & \text{for } (x - x^*) \cdot \vec{n} \ge 0, \\ -1 & \text{otherwise.} \end{cases} \qquad (1)$$

In Eq. (1), $x$ is a given point of the $(X, Y, Z)$ space, $x^*$ is a point of the crack surface and $\vec{n}$ is a vector in $x^*$ normal to the crack surface (see Fig. 3(a)).

The enrichment $F_J$ is used to enrich the nodes of the elements that contain the crack front. $F_J$ consists of a span of functions which incorporate the radial and angular behaviour of the asymptotic crack-tip displacement fields. $F_J$ is given by the following equation where $r$ and $\theta$ are expressed in the local polar coordinate system of the crack front [34] (see Fig. 3(a)):

$$\{F_J(\mathbf{x})\} = \left\{ \begin{array}{l} \sqrt{r} \sin\left(\frac{\theta}{2}\right), \sqrt{r} \cos\left(\frac{\theta}{2}\right), \sqrt{r} \sin\left(\frac{\theta}{2}\right) \sin(\theta), \\ \sqrt{r} \cos\left(\frac{\theta}{2}\right) \sin(\theta) \end{array} \right\}. \qquad (2)$$

Fig. 3(b) illustrates the enrichment strategy in two dimensions when the crack is not aligned with the element edges. The key issue is now to define the 3D crack geometry independently of the mesh of the structure in order to select the nodes that must be enriched by $H$ and $F_J$.

### 2.6.2. Implicit representation of the crack using the level sets method

The crack front is defined as the intersection between a surface $\Gamma_{surface}$ that defines the crack surface and a surface $\Gamma_{front}$ that defines the crack front. These two surfaces, defined for a semi-elliptical crack, are displayed in Fig. 4.

In the general case, these surfaces can have arbitrary 3D shape, and are defined by means of level sets functions. The concept of level sets functions was initially introduced by Osher et al. to track moving interfaces [32]. A level sets function $\gamma(x)$ (where $x$ is a point of the $(X, Y, Z)$ space), used to represent an arbitrary surface $\Gamma$, can be defined as a scalar
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) Part of the crack front and its local polar coordinate system used to express the enrichment function $F_J$ and $H$ as well as the displacement fields. The vector $n$ is the normal in $x^*$ to the crack surface, and $x$ is a point of the $(X, Y, Z)$ space inside the structure. (b) 2D schematic representation of the crack (black curve) and the mesh of the structure that contains the crack (squared grid). The nodes of elements cut by the crack are enriched by the $H$ function and are represented by open circles. The nodes of the elements that contain the crack tip are enriched by the $F_J$ functions span and are represented by open rectangles.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. The surface in the $(X-Y)$ plane is the iso-zero surface of the level set function $\phi(x)$ and represents the crack surface. The “curved” surface normal to the $(X-Y)$ plane is the iso-zero surface of $\psi(x)$ and represents the crack front, which corresponds to the intersection of these two iso-zero surfaces.</figcaption></figure>

field, where each scalar is the signed distance to the surface $\Gamma$. Two level sets functions, $\phi$ and $\psi$, define the surfaces $\Gamma_{\text{surface}}$ and $\Gamma_{\text{front}}$, respectively, which give the following equations:

$$ \begin{aligned} \phi(x) &= \pm \|x - x_{\Gamma_{\text{surface}}}\| \quad \text{where } x_{\Gamma_{\text{surface}}} \in \Gamma_{\text{surface}} \quad \text{and} \\ \psi(x) &= \pm \|x - x_{\Gamma_{\text{front}}}\| \quad \text{where } x_{\Gamma_{\text{front}}} \in \Gamma_{\text{front}}, \end{aligned} \quad (3) $$

$\Gamma_{\text{surface}}$ is the iso-zero surface of the scalar field $\phi: \phi(x) = 0$. Similarly, $\Gamma_{\text{front}}$ is the iso-zero surface of $\psi: \psi(x) = 0$ see (Fig. 4).

The values of $\psi$ and $\phi$ are stored at each node of the elements of the structure. The nodal values of $\psi$ and $\phi$ give the precise location of the crack and control whether a node has to be enriched or not. In the study presented here, the level sets functions $\psi$ and $\phi$ are built from the real crack shape obtained from the 3D tomographic image by calculating the signed distance function to $\Gamma_{\text{surface}}$ and $\Gamma_{\text{front}}$ using a C++ subroutine.

The values of the stress intensity factors, $K_I$, $K_{II}$ and $K_{III}$, along 3D curved cracks are determined using a general crack-tip interaction integral introduced by Gosz and Moran [35]. Details concerning the calculation of the mixed-mode stress intensity factors can be found, for instance, in Refs. [33–35].

### 2.6.3. Advantage of the X-FEM

The numerical method used in this study is very general and allows the calculation of the stress intensity factors $K_I$, $K_{II}$ and $K_{III}$ for a crack with an arbitrary 3D geometry. Because of the addition of local enrichment functions, accurate values of the stress intensity factors can be obtained with a relatively coarse mesh compared to the mesh used at the crack tip in conventional FE methods. Moreover, the crack geometry and the FE mesh of the structure are defined separately, and no remeshing of the entire structure is necessary when the crack size/shape changes.

## 3. Results

After 14,000 fatigue cycles, a crack initiated at the notch was detected on a 2D radiograph. A tomographic scan of the part of the specimen containing the crack was acquired and corresponds to the “first scan”. During crack propagation, nine scans were recorded, reconstructed and segmented to obtain 3D renditions of the crack at different stages of its evolution. Projections in the $(X-Y)$ plane of seven 3D renditions (chosen among the nine) are displayed in Fig. 5. The dotted lines on the 3D crack renditions represent the location of the crack front at the previous step. The position along the crack front is defined by the angle $\omega$ as shown in Fig. 5(a). Qualitative and quantitative information on the evolution of the 3D crack geometry obtained from the 3D renditions of the crack is presented in this section.

### 3.1. Crack shape evolution

The crack plane, as shown in Fig. 2, appears almost perfectly flat and is normal to the $z$ axis which is parallel to the stress direction. The crack front shape remains close to a
<figure><img src="image_7.png" /><figcaption>Fig. 5. Projections in the (X-Y) plane of 3D renditions of the crack at different stages of its development. (a) The angle $\omega$ is used to determine the location of a given point along the crack front. (b) The black dotted lines on the 3D rendition represent the location of the crack front at the previous stage (a), and similarly for stages (c)–(g).</figcaption></figure>

semi-ellipse during crack propagation. The aspect ratio $a/c$, defined as the ratio of the crack depth $a$ and the half surface length $c$, as defined in Fig. 5(a), is increasing steadily from 1.1 to 1.4 during crack growth (see Table 1). The effect of the 100% tensile overload on the crack front shape (see Fig. 5(d)) is discussed in Section 3.3.

<table><caption>Table 1
Evolution of the aspect ratio during crack propagation</caption><thead><tr><th>N</th><th>14,000</th><th>16,000</th><th>18,000</th><th>19,000</th><th>20,000</th><th>22,000</th><th>24,000</th></tr></thead><tbody><tr><td>a/c</td><td>1.19</td><td>1.23</td><td>1.29</td><td>1.35</td><td>1.37</td><td>1.39</td><td>1.44</td></tr><tr><td colspan="8">↑<br>100% tensile overload</td></tr></tbody></table>

## 3.2. Fatigue crack growth rate at the surface and in the bulk

The crack size at the surface (2c), for $\omega = 0^{\circ}$, and in the bulk (a), for $\omega = 90^{\circ}$, are measured on the 3D renditions of the crack (see Fig. 5). The fatigue crack growth rates at the surface, $dc/dN$, and in the bulk, $da/dN$, as a function of the stress intensity factor range $\Delta K = (K_{max} - K_{min})$ are shown in Fig. 6. The values of $K_{max}$ and $K_{min}$ are calculated for $\omega = 0^{\circ}$ and for $\omega = 90^{\circ}$ using the X-FEM presented in Section 2.6. Note that $K_{max}$ and $K_{min}$ are the stress intensity factor under mode I, the values of $K_{II}$ and $K_{III}$ are taken equal to 0 as the crack, subjected to uniaxial tension loading, is normal to the stress axis. For comparison, the

<figure><img src="image_8.png" /><figcaption>Fig. 6. The fatigue crack growth rates measured at the surface ($\omega = 0^{\circ}$) and at the deepest point ($\omega = 90^{\circ}$) of the crack are represented by plain black points and plain triangles, respectively. The sudden drops between 18,000 and 20,000 cycles correspond to post-overload retardation. The fatigue crack growth law measured for a CT specimen in Ref. [19] is given for comparison.</figcaption></figure>
fatigue crack growth law obtained on CT specimens by Venkateswara Rao and Ritchie [26] for the same material is also represented by open squared symbols in Fig. 6.

A discrepancy between the fatigue crack growth rate measured at the surface and in the bulk can be observed in Fig. 6: the slopes of the two curves $da/dN = f(\Delta K)$ and $dc/dN = f(\Delta K)$ are similar but the former is shifted towards lower values of $\Delta K$. Thus for the same $\Delta K$, the crack propagates faster in the bulk than at the surface implying that, at least for the specimen geometry used here, the crack growth behaviour is anisotropic. Thus, a fatigue crack growth law determined at the surface does not account for the crack behaviour in the bulk. Furthermore, if used for fatigue life calculation, the surface experimental law would lead to a non-conservative prediction.

### 3.3. Overload

After the acquisition of the tomographic scan at 18,000 cycles, a 100% tensile overload was applied to the sample. In Figs. 5(d) and (e), it can be seen that, during the 2000 cycles following the overload, the crack growth is considerably reduced, especially at the surface where the crack propagates three time less than in the bulk.

This discrepancy between retardation at the surface and in the bulk can be quantified by measuring the amplitude of crack retardation on the fatigue crack growth rate curves displayed in Fig. 6. Both curves exhibit a sudden drop between 18,000 and 20,000 cycles: the $dc/dN = f(\Delta K)$ curve decreases from $9 \times 10^{-9}$ m/cycle down to $4 \times 10^{-9}$ m/cycle whereas the $da/dN = f(\Delta K)$ curve only decreases from $1.4 \times 10^{-8}$ m/cycle to $1.26 \times 10^{-8}$ m/cycle. After cycle 20,000, the crack growth rate recovers its original trend. From this it can be concluded that post-overload retardation appears to be more important at the surface than in the bulk.

## 4. Simulation of the 3D crack propagation

In the previous section it has been shown that crack propagation at the point ($\omega = 90^\circ$) furthest into the bulk is faster than at the surface ($\omega = 0^\circ$). The objective of this section is to investigate the possibility of modelling crack growth anisotropy in order to predict more precisely the evolution of the crack front shape in the bulk.

### 4.1. Analytical and numerical simulation of crack propagation with a single crack growth law

The discrepancy between the fatigue crack growth rates at the surface and in the bulk allows us to conclude that a unique crack growth law (e.g., with unique coefficients C and m for the Paris law) cannot be used to predict the evolution of the crack front shape in 3D. It will indeed be necessary to establish a “3D crack propagation law” to simulate anisotropic crack propagation.

However, as a preliminary study, we will simulate crack propagation by integrating the same Paris law over the

crack front, as is usually done in crack simulation. By doing so, it will be possible to determine if the discrepancy between simulated and real crack shape is localised at the sample surface or if it affects the whole crack front. The law used for simulation is the fatigue crack growth law measured for $\omega = 90^\circ$. This law, defined by Eq. (4) and represented by a dotted line in Fig. 6, has been chosen for the basis of the modelling, rather than the “surface” law, because numerical calculations of the stress intensity factors are more reliable in the bulk than at the surface. (This point is discussed in Section 5.)

$$ \frac{da(\omega)}{dN} = C \cdot [K_{\max}(\omega) - K_{\min}(\omega)]^m \quad (4) $$
with $C = 10^{-9.2}$ and $m = 3.5$.

In Eq. (4), $K_{\max}(\omega)$ and $K_{\min}(\omega)$ are the values of the local stress intensity factor calculated at a given point of the crack front localised by the angle $\omega$. The constants C and m have been determined from the $da/dN = f(\Delta K)$ curve in Fig. 6 but without considering the points between cycles 18,000 and 20,000 which are characteristic of post-overload retardation. In the following section, the evolution of the crack front shape between cycles 14,000 and 18,000 will be simulated using both X-FEM (cf. Section 2.6) and an analytical method based on the assumption that the crack front shape is exactly semi-elliptical.

In the numerical simulation, the stress intensity factors are calculated along the crack front using the X-FEM, taking into account the measured 3D crack geometry. The evolution of $K_{\max}(\omega)$ along the crack front is obtained by fitting the “raw” numerical values with a polynomial function. Note that the values of $K_{\max}(\omega)$ calculated at the first points located just under the sample surface are not taken into account in the polynomial interpolation. The local crack extension $da(\omega)$ is normal to the crack front and its value is calculated at each node by integrating the crack growth law defined by Eq. (4). The new crack front location is calculated for a difference in number of cycles $dN = 4000$. The initial crack is the experimental crack at cycle 14,000 obtained from the first tomographic scan.

In the analytical simulation, the simulated crack front is assumed to be a perfect semi-ellipse: the crack extension $da(\omega)$ is calculated by integrating Eq. (4) only at the two points located at $\omega = 0^\circ$ and $\omega = 90^\circ$. The values of $K_{\max}(0^\circ)$ and $K_{\max}(90^\circ)$ are obtained from the Newman and Raju equations for a semi-elliptical crack in tension [36]. The initial crack is a semi-ellipse where a and c are the crack size at $\omega = 0^\circ$ and $\omega = 90^\circ$, respectively, as measured on the experimental crack at cycle 14,000.

### 4.2. Comparison between the simulated and experimental crack front shapes

In the following analysis the simulated crack front obtained using analytical values of K is called the “analytical crack front” and the crack front simulated with the numerical method is called the “numerical crack front”.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Simulation of crack propagation using a single Paris law given by Eq. (4) between cycles 14,000 and 18,000. The crack represented by open squared symbols corresponds to the numerical simulation using the stress intensity factors given by the X-FEM. The crack represented by black points is obtained by “analytical simulation” using the Raju and Newman equations to calculate the stress intensity factors.</figcaption></figure>

The experimental crack fronts extracted from the tomographic images (solid lines), the numerical crack front (open squared symbols) and the analytical crack front (filled circles) are compared in Fig. 7. This figure shows that the shape of the numerical crack front, obtained without any shape assumption, is comparable to the perfect semi-ellipse obtained with the analytical method (except for $0° < ω < 15°$). The discrepancy between the “experimental” and “simulated” crack fronts, both analytical and numerical, is a maximum at the surface and tends to decrease towards the bulk. The crack size for $ω = 90°$ is overestimated only by 3% in the numerical simulation and by 10% in the analytical simulation. The crack size at the surface is overestimated by 35% in the analytical simulation and by 25% in the numerical simulation. From this it can be concluded that the discrepancy between the simulated and experimental crack shapes, when integrating a unique law over the whole crack front, is not local but affects the shape of the entire crack between $ω = 0°$ and $ω = 180°$. This discrepancy is more important when the crack is described as a perfect semi-ellipse than when the real crack shape is taken into account.

### 4.3. Numerical simulation using a 3D effective crack growth law

One possible reason for the observed crack growth anisotropy between the surface and the bulk is the variation of the closure stress along the crack front: the crack will then propagate faster in the bulk because the closure stress in this region of the crack front is smaller and hence the effective driving force is higher than at the surface. The validity of this assumption as well as the 3D nature of the closure phenomena are discussed in Section 5.

The objective of this section is to establish a 3D crack propagation law that accounts for the crack growth anisotropy induced by closure. The local crack extension $da(ω)$ will be correlated with an effective stress intensity factor $ΔK_{eff}(ω)$ that takes into account the 3D variation of the closure stress along the crack front.

Elber [37] proposed a modified Paris equation to account for the effect of closure on the crack growth rate: $da/dN = C(ΔK_{eff})^m$ where $ΔK_{eff} = (K_{max} - K_{op})$ is the effective stress intensity factor range and $K_{op}$ is the stress intensity factor at which the crack opens. To account for the variation of the closure stress along the crack, we will modify Elber’s equation by introducing a new variable $ω$ such that $da/dN$ and $ΔK_{eff}$ vary along the crack front: $da(ω)/dN = C(ΔK_{eff}(ω))^m$.

As a first approach, the variation of $K_{op}(ω)$ is taken to be linear with the angle $ω$ and is therefore given by $K_{op}(ω) = a · ω + b$. The value of the constants $a$ and $b$ are determined from the values of $K_{op}$ at the surface $K_{op}(0°)$ and at the deepest point $K_{op}(90°)$. Data from the literature are used to determine $K_{op}(0°)$: the closure response of an Al–Li powder metallurgy alloy, very similar to the alloy studied here, was investigated in Ref. [26]. The value of $K_{op}$ measured at the surface of CT specimens, for $R = 0.1$, is found to be equal to $0.4K_{max}$. We assume here that the ratio $K_{op}(ω = 0°)/K_{max}$ remains constant during the fatigue test and that there is no closure at the deepest point for $ω = 90°$, which gives $K_{op}(90°) = 0.1K_{max}(90°)$. The evolution of the closure stress along the crack front is given by the linear relation

$$K_{op} = -\frac{0.6K_{max}}{\pi} \cdot \omega + 0.4K_{max} \quad (5)$$

The 3D crack growth law taking into account the variation of the closure stress along the crack front is given by

$$\frac{da(\omega)}{dN} = C \left[ \left( \frac{0.6}{\pi} \omega + 0.6 \right) K_{\max} \right]^m \quad (6)$$

In Eqs. (5) and (6), the angle $ω$, expressed in radians, is defined for $[0, π/2]$ and $K_{max}(ω)$ is obtained from a polynomial interpolation of the values calculated by the X-FEM as described in the previous section.

The crack propagation is simulated, using Eq. (6), between 14,000 and 24,000 cycles by calculating the location of the simulated crack front every 2000 cycles. The initial crack front is the crack front at 14,000 cycles obtained from the tomographic images (see Fig. 8(a)). The local crack extension values $da(ω)$ are calculated along this initial crack front for $dN = 2000$ to predict the location of the crack front at 16,000 cycles (open triangles in Fig. 8(a)). The process is iterated as follows. The $da(ω)$ values are calculated along the simulated crack front at step $N$ (black triangles) to establish the location of the crack front at step $(N + dN)$ (open triangles). The obtained simulated crack fronts are compared to the experimental crack front (plain grey lines) in Figs. 8(b) and (d). However, as Eq. (6)
<figure><img src="image_10.png" /><figcaption>Fig. 8. (a)–(d) The plain grey lines represent the location of the experimental crack front obtained from the tomographic image. The open triangles represent the simulated crack fronts obtained by integrating, for 2000 cycles, the 3D crack growth law defined by Eq. (6). (a, c) The initial crack front used for propagation is the experimental crack front after 14,000 and 20,000 cycles, respectively. (b, d) The initial crack front used for propagation is the simulated crack front (black triangles) that correspond to 16,000 and 22,000 cycles, respectively.</figcaption></figure>

does not model the crack retardation induced by the over-
load, crack propagation is not simulated between 18,000
and 20,000 cycles. The experimental crack front at step
20,000 cycles is, hence, used to predict the shape of the
crack front at step 22,000 and 24,000 cycles (see Figs.
8(c) and (d)).

In Figs. 8(a)–(d), the simulated and experimental crack
front shapes appear to be very similar even in the region
close to the sample surface: in this region the average error
between simulated and experimental crack size is 6%,
whereas at the depth the error is less than 2%. These values
are far less than the error found between experimental and
simulated crack sizes calculated using a single Paris law as
described in the previous section. It can thus be concluded
that the 3D crack propagation law, based on a linear var-
iation of K<sub>op</sub> along the crack front, gives a good prediction
of the crack growth, at least for the crack/sample geometry
considered in this study.

5. Discussion

5.1. *Long crack behaviour and applicability of LEFM*

The approach presented in this study relies on the
assumption that small-scale yielding conditions apply
despite the physically short dimensions of the crack
(100–500 µm). This assumption seems justified since the
estimated monotonic plastic zone size is “small” com-
pared to the crack size. Indeed, the ratio of the monotonic
plastic zone size r<sub>p</sub> and the crack size has been calculated
at the surface r<sub>pc</sub>/c (assuming plane stress conditions) and
in the bulk r<sub>pa</sub>/a (assuming plane strain conditions). For
the crack length considered in this study the values range
between 0.09 < r<sub>pc</sub>/c < 0.13 and 0.016 < r<sub>pa</sub>/a < 0.018. At
the surface, these values are rather high but remain
acceptable. In the bulk, they seem to conform well to
the applicability limits of LEFM that can be found in
the literature [38]. It seems therefore reasonable to corre-
late the stress intensity factor range calculated along the
crack front with the fatigue crack growth rate. It should
be remembered that the values of the stress intensity fac-
tors at the free surface/crack intersection calculated with
numerical methods should be treated cautiously and are
usually obtained, as already mentioned by several
authors, by extrapolating values calculated “inside” the
sample [39,40].

5.2. Crack growth rate anisotropy between surface and bulk

Concerning the observed anisotropy between the fatigue
crack growth rates measured at the surface and in the bulk,
two remarks can be made. First, the fatigue crack growth
rate measured at the surface of the studied sample is close
to the fatigue crack growth rate obtained for long cracks in
conventional CT specimens [26]. This result suggests that
the crack propagation observed on a small sample is repre-
sentative of the long crack behaviour described by LEFM
in larger samples of standard geometry.

Second, the crack propagates faster in the bulk than at
the surface. Such crack growth anisotropy has already been
observed in 7075-T651 Al alloy for semi-elliptical cracks by
McDonald and Daniewicz [41]. To account for these
results, Solanki et al. have established a “local” Paris
law, defined by local coefficients Cᵢ and mᵢ [42]. The values
of Cᵢ and mᵢ were deduced by linear interpolation between
the experimental values of C and m, measured in the depth
and at the free surface. However, in our case, crack prop-
agation is expected to be isotropic as the grain structure
of the Al–Li alloy studied is ultrafine and equiaxed and
possesses no crystallographic texture. Thus, there is no
“physical” justification for allowing the C and m coeffi-
cients to vary along the crack front as they are intrinsic
material parameters. As mentioned in Section 4.1, it was
therefore assumed in this study that crack closure is the ori-
gin of the observed crack growth anisotropy.

5.3. 3D nature of closure: experimental results

As conventional methods of closure stress measure-
ments are all surface measurements, the experimental char-
acterisation of the 3D nature of closure phenomena is a
difficult task which requires the development of specific
methods to evaluate the contribution of surface and bulk
to the global closure response of the material. Qualitative
results have been obtained for instance in Refs. [43,44]:
specimens were sectioned normally to the crack propaga-
tion plane to observe directly the crack front through the
specimen thickness. These experiments confirm the fact
that the closure stress is higher at the surface than in the
bulk. This phenomenon has been attributed to the varia-
tion of constraint along the crack front: plane stress condi-
tions prevail at the surface which implies a larger plastic
zone size, and thus more closure than in the bulk, where
plane strain conditions prevail (a trend reflected by the val-

ues of r<sub>pa</sub>/a and r<sub>pc</sub>/c given previously). The difference
between the closure stress at the surface and in the bulk
has been determined quantitatively using a fractographic
method developed by Pelloux et al. [45]. This technique
relies on the assumption that fatigue striation spacing is
directly affected by the closure stress which can hence be
determined by comparing different striation spacings
induced by a special load sequence of varying amplitude.
This technique has been used to investigate the 3D closure
behaviour of through cracks, as for instance in Refs. [46,47]
and part-through cracks in Ref. [9]. In this latter paper,
Putra et al. measure the evolution of the opening stress,
S<sub>op</sub> = f(φ), along the crack front of semi-elliptical cracks
in 7075 T6 Al alloy (the angle φ used in Ref. [9] has the
same meaning as ω used in the present paper). However,
as fatigue striations close to the sample surface are not well
defined, this fractographic method does not allow the
determination of values of S<sub>op</sub> between φ = 0° and φ
= 15° (which corresponds to a thickness of 25–50 µm for
the crack size range observed in our experiment). In the
bulk, experimental values of S<sub>op</sub> decrease linearly between
φ = 15° and φ = 90°, which confirm our assumption, but
with a slope that is lower than the one of Eq. (5) used
for the simulation.

5.4. 3D FE simulations of plasticity-induced crack closure phenomena

For the particular ultrafine-grained Al alloy used in this
study, the main mechanism responsible for closure is plas-
ticity-induced crack closure (PICC) [23]. PICC phenomena
for semi-elliptical cracks have been analysed in 3D by the
FE method [48–52]. 3D elastic–plastic FE simulations show
that the loading part of the fatigue cycle can be described as
a “continuous unzipping process” during which the crack
opens first in the bulk. They also confirm the fact that S<sub>op</sub>
is higher at the surface than in the bulk. In Ref. [49], the evo-
lution of the closure stress has been calculated: the obtained
S<sub>op</sub> = f(ω) curve does not follow the linear trend of Eq. (5)
used in this study but decreases between ω = 0° and ω = 10°
and then remains at a constant value for 10° < ω < 90°. As
the crack/sample size ratio of the sample used in Ref. [49] is
10 times less than for our sample, this difference might be
due to the influence of the state of stress as will be explained
in the next section. It should also be emphasised that
numerical values of S<sub>op</sub> are greatly influenced by the mesh
refinement as well as the crack node release scheme used
in the FE simulation as mentioned in Ref. [50].

5.5. Influence of the state of stress on the 3D closure
behaviour

The influence of the state of stress on the variation of
S<sub>op</sub> has been studied for through cracks by comparing
elastic–plastic FE calculations performed on thin and
thick MT samples [53] (the thickness ratio between the
thick and the thin sample being 7.5). It was found that,
in thin MT samples, where plane stress conditions prevail, $S_{op}$ decreases linearly between the sample surface and the centre of the specimen. For thicker MT samples, the evolution is different: $S_{op}$ decreases quickly in a region located just under the sample surface and then remains at a constant value in the central part of the sample where plane strain conditions prevail. It can be concluded that for through cracks, the variation of $S_{op}$ along the crack front is linear only when plane stress conditions prevail in the whole sample. In our case, as the studied crack is not a through crack but a part-through crack, direct comparison is difficult. However, the results obtained in Ref. [53] seem to explain the difference, mentioned in the previous paragraph, between the $S_{op} = f(\omega)$ curve obtained by FE calculation for a large sample [49] and the linear equation (5) used in the simulation. It can thus be concluded that the influence of the state of stress on the variation of $S_{op}$ along the crack front should be taken into account, which implies that Eq. (5) should be modified whenever a different crack/sample geometry is considered.

### 5.6. Overload

It has been shown in Section 3.3 that crack retardation following the 100% overload varies along the crack front: it is maximum at the surface for ($\omega = 0°$) and decreases towards the bulk to a minimum value for ($\omega = 90°$). The fact that post-overload retardation occurs mainly at the surface has been observed on fracture surfaces, but only for through cracks in Refs. [54,55].

Retardation induced by a tensile overload can be attributed to different mechanisms such as compressive residual stress, branching and also PICC [56]. In the fine-grained alloy studied here, PICC is the main mechanism responsible for post-overload crack retardation [23]. The fact that crack retardation occurs mainly at the surface, as described above, implies that PICC is more important and, hence, the value of $K_{op}$ higher at the surface than in the bulk. This result is consistent with the assumption made in Section 4.2 that $K_{op}(0°) > K_{op}(90°)$ which was used to establish the 3D crack growth law.

## 6. Conclusions

The 3D propagation of a semi-elliptical fatigue crack in an ultrafine-grained Al–Li alloy has been visualised in situ using synchrotron radiation X-ray microtomography. This technique can be used to perform in situ characterisation of fatigue cracks in the bulk of opaque materials at a resolution close to 1 µm, which is relevant for studying crack propagation. It should be emphasised that the 3D images of the cracks have been obtained without changing the load cycles or the environmental conditions contrary to what is achieved with classic post-mortem fractography techniques. However, the stringent resolution requirements for detection of crack opening displacements at the sub-micrometre level limit the field of view to about 1 mm

and hence imply the adoption of special sample geometries with dimensions in the millimetric range.

The experimentally observed 3D crack propagation has been modelled using LEFM concepts: the local crack growth rate measured in the bulk is correlated with the stress intensity factor range $\Delta K$ calculated using the X-FEM taking into account the real 3D geometry of the crack obtained from the tomographic images. It has been shown that:

* Crack growth rates measured at the surface of the small sample used in our experiment are similar to those obtained on large conventional fatigue samples for the same alloy.

* The experimentally observed semi-elliptical crack propagates faster in the bulk than at the surface.

* Variation of the closure stress along the crack front can account for the observed crack growth anisotropy between surface and bulk. As PICC is the main mechanism responsible for closure in the alloy studied here, a larger plastic zone size at the surface (plane stress) can explain a higher closure stress than in the bulk, where plane strain conditions limit the development of plasticity and hence closure.

* Assuming that the variation of the closure stress between surface and bulk is linear, a 3D crack growth law equation has been established to simulate crack propagation. The predicted and experimental crack shapes, compared at different stages, appear to be very similar.

This means that, for the sample geometry investigated here, a single Paris law can be used to predict the observed crack growth anisotropy provided the variation of the closure stress along the crack front is taken into account.

From a more general point of view, this study has shown that coupling X-ray microtomography to X-FEM provides a promising tool to assess the 3D behaviour of arbitrary shaped cracks and to perform direct comparisons of “experimental” and “simulated” crack shapes during propagation. Such data, to the best of the authors’ knowledge, are currently lacking in the literature.

## Acknowledgements

The authors are grateful to Mr. Bultreys from FEI Europe for machining the notch with the help of a FEI dual-beam FIB/SEM system and to Ian Sinclair for fruitful discussions on the 3D nature of the closure phenomenon. We also like to thank Andy Tarrant of Aerospace Metal Composites Ltd for providing the material. The ESRF and staff are acknowledged for their help and assistance during the tomographic imaging experiment.

## References

[1] Enmark M, Lucas G, Odette GR. J Nucl Mater (The Netherlands) 1992;1038–41.

[2] Makoto A, Masumi S. Eng Fract Mech 2005;72:319–34.
[3] Knight MJ, Brennan FP, Dover WD. NDT&E Int 2004;37:337-43.

[4] Zerwer A, Polak MA, Santamarina JC. J Struct Eng 2001;128:240-8.

[5] Ravichandran AS, Larsen JM. Mater Sci Eng A 1992;153:499-507.

[6] Forman RG, Mettu SR. ASTPM STP 1131 1992;1:519-46.

[7] Nadot Y, Ranganathan N, Mendez J, Béranger AS. Scr Mater 1997;37:549-53.

[8] Husset J, Lieurade HP, Maltrud F, Truchon M. Welding in the World/Soudage dans le Monde 1985;23:276-83.

[9] Putra IS, Schijve J. Fatigue Fract Eng Mater Struct 1992;15:323-38.

[10] McFadyen NB, Bell R, Vosikovsky O. Int J Fatigue 1988;10:47-9.

[11] Forman RG, Shivakumar V. ASTM STP 905 1986:59-74.

[12] Scott PM, Thorpe TW. Fatigue Eng Mater Struct 1981;4:291-309.

[13] Pommier S, Sakae C, Murakami Y. Int J Fatigue 1999;21:243-51.

[14] Fatigue crack growth computer program "Nasgro" version 3.0.
Reference manual. March; 2002.

[15] Wawrzynek PA, Martha LF, Ingraffea AR. In: Rosakis AJ et al., editors. Analytical, numerical and experimental aspects of three dimensional fracture processes, AMD-vol. 91. New York (NY): ASME; 1988. p. 321-7.

[16] Fulland M, Richard HA. Steel Res 2003;74:584-90.

[17] Wessel C, Cisilino A, Santi O, et al. Theor Appl Fracture Mech 2001;35:47-58.

[18] Riddell WT, Ingraffea AR, Wawrzynek PA. Eng Fract Mech 1997;58:293-310.

[19] Ludwig W, Buffière JY, Savelli S, Cloetens P. Acta Mater 2003;51:585-98.

[20] Marrow TJ, Buffière JY, Withers PJ, Johnson G, Engelberg D. Int J Fatigue 2003;26:717-25.

[21] Miller KJ. Fatigue Eng Mater Struct 1982;5:223-32.

[22] Pitcher PD, Bushby RS, Vine WJ, Smith AF, Tarrant AD. Mater Sci Technol 2001;17:807-14.

[23] Bray GH, Reynolds AP, Starke EA. Metall Trans A 1992;23:3055-66.

[24] Bray GH, Kaisand LR, Starke EA. Fatigue Fracture Eng Mater Struct 1995;18:551-64.

[25] Vine WJ, Pitcher PD, Tarrant AD. Mater Sci Technol 2001;17:645-50.

[26] Venkateswara Rao KT, Ritchie RO. Metall Trans A 1991;22:191-202.

[27] Sinclair I. Private communication. Unpublished work, School of Engineering University of Southampton, UK; 2005.

[28] Cloetens P, Pateyron-Salomé M, Buffière JY, Peix G, Baruchel J, Peyrin F. J Appl Phys 1997;81:58-78.

[29] Buffière JY, Ferrié E, Proudhon H, Ludwig W. Mater Sci Technol [in press].

[30] Amira 3.1. Mercury Computer System GmbH, Lepiusstr., 70 D-12163 Berlin, Germany. Available from: http://www.amiravis.com/.

[31] Belytschko T, Black T. Int J Numer Meth Eng 1999;45:601-20.

[32] Osher S, Sethian JA. J Comput Phys 1988;79:12-49.

[33] Moes N, Dolbow J, Belystchko T. Int J Numer Meth Eng 1999;46:131-50.

[34] Gravouil A, Moes N, Belytschko T. Int J Numer Meth Eng 2002;53:2569-86.

[35] Gosz M, Moran B. Eng Fract Mech 2002;69:299-319.

[36] Newman JC, Raju IS. ASTM STP 791 1983:238.

[37] Elber W. Eng Fract Mech 1970;2:37-45.

[38] Anderson TL. Fracture mechanics – fundamentand and applications. 2nd ed. Boca Raton, (FL): CRC Press; 1995.

[39] Smith RA, Cooper JF. Int J Press Ves 1989;36:326-35.

[40] Lin XB, Smith RA. Eng Fract Mech 1999;69:503-22.

[41] McDonald Jr V, Daniewicz SR. ASTM STP 1406 2001:381-96.

[42] Solanki K, Daniewicz DR, Newman Jr JC. Eng Fract Mech 2004;71:149-71.

[43] Lyndley TC, Richards CE. Mater Sci Eng 1974;14:281-93.

[44] Fleck NA, Smith RA. Int J Fatigue 1982;4:157-60.

[45] Pelloux RM, Faral M, Mc Gee WM. Fatigue Eng Mater Struct 1979;21-35.

[46] Dawicke DS, Grandt AF, Newman JC. Eng Fract Mech 1990;36:11-121.

[47] Sunder S, Dash K. Int J Fatigue 1992;3:97-105.

[48] Chermahini RG, Palmberg B, Blom AF. Int J Fatigue 1993;4:259-63.

[49] Skinner JD, Daniewicz SR. Eng Fract Mech 2002;69:1-11.

[50] Solanki K, Daniewicz SR, Newman JC. Eng Fract Mech 2004;71:149-71.

[51] Hou CY. Int J Fatigue 2004;26:1225-39.

[52] Zhang JZ, Bowen P. Eng Fract Mech 1998;60:341-60.

[53] Chermahini RG, Shivakumar KN, Newman Jr JC, Blom AF. Eng Fract Mech 1989;34:393-401.

[54] Matsuoka S, Tanaka K. Eng Fract Mech 1980;13:293-306.

[55] Bathias C, Bäilon JP. In: Hermès, editor. La fatigue des matériaux et des structures; 1997.

[56] Suresh S. Eng Fract Mech 1983;18:577-93.