# 2D fatigue crack propagation in rails taking into account actual plastic stresses

Benoît Trollé, Marie-Christine Baietto, Anthony Gravouil, Si Hai Mai, Benoit Prabel

► To cite this version:

Benoît Trollé, Marie-Christine Baietto, Anthony Gravouil, Si Hai Mai, Benoit Prabel. 2D fatigue crack propagation in rails taking into account actual plastic stresses. Engineering Fracture Mechanics, 2014, 123 (SI), pp.163-181. 10.1016/j.engfracmech.2014.03.020. hal-01951565

HAL Id: hal-01951565

https://hal.science/hal-01951565

Submitted on 7 Jul 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
2D fatigue crack propagation in rails taking into account actual
plastic stresses

B. Trollé<sup>a,b,*</sup>, M.-C. Baietto<sup>a</sup>, A. Gravouil<sup>a,c</sup>, S.H. Mai<sup>b</sup>, B. Prabel<sup>d</sup>

Due to increase axle loads, repeated passages of the wheels, rolling contact fatigue cracks initiate in the surface or subsurface of the rails. These defects can propagate and lead to the rail failure. A two-scale frictional contact fatigue crack model developed within the X-FEM framework is used to address the cracked rail problem, calculate the stress intensity factors, perform crack propagation and fatigue life prediction. Actual plastic stresses, quantified through a dedicated software developed by SNCF within the consortium IDR2, are taken into account in the propagation simulation via projection of the asymptotic mechanical fields. The effects of those actual plastic stresses are investigated on the crack growth path and rate. Interactions between multiple cracks of a network are also analyzed.

# 1. Introduction

The repeated train traffic leads to crack initiation at the rail surface or within the subsurface, up to the development of very complex 3D crack network, like squats and head-checks. Rail failure may occur and potentially derailment too. Costly maintenance operations are deployed to avoid such a situation.

Rolling contact loading is a time-dependent, multi-axial and nonproportional loading. This cyclic loading leads to an asymptotic stress state within the rail that influences the crack behavior. The cracks are also submitted to large compressive stresses inducing crack face closure and frictional contact.

Previous works about fatigue crack growth in the rails are available in the literature. The role of liquid entrapment using finite element analysis [1–3] or boundary element method [4,5] was analyzed on the crack growth mechanism. The influence of different parameters, such as elastic foundation [3,6], the initial crack geometry or the crack face friction coefficient [3,5,7,8] was also investigated on the stress intensity factors (SIFs) or on the crack propagation path [9]. It is well known that the distribution of the residual stress field affects the crack propagation behavior [10–15], but most of previous studies use simplified techniques. This residual stress field is accounted as the addition of uniform tensile or compressive conditions at the boundary of the considered domain [8,16,17].

To improve the understanding of the crack initiation and propagation mechanisms, a global methodology has been developed thanks to a long-term collaboration between French railway organizations (SNCF, RATP, RFF), rail producer (Tata Steel)

<footnote><sup>a</sup>Université de Lyon, CNRS INSA-Lyon, LaMCoS, UMR5259, 20 Avenue Albert Einstein, F69621 Villeurbanne Cedex, France</footnote>

<footnote><sup>b</sup>SNCF, Direction Innovation et Recherche, 40 Avenue des Terroirs de France, 75611 Paris Cedex 12, France</footnote>

<footnote><sup>c</sup>Institut Universitaire de France, France</footnote>

<footnote><sup>d</sup>CEA, DEN, DM2S, SEMT, DYN, CEA-Saclay, 91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote>* Corresponding author at: Université de Lyon, CNRS INSA-Lyon, LaMCoS, UMR5259, 20 Avenue Albert Einstein, F69621 Villeurbanne Cedex, France. Tel.: +33 171323324.</footnote>

<footnote>E-mail address: benoit.trolle@sncf.fr (B. Trollé).</footnote>
Nomenclature

<table><tr><td>ε(0)</td><td>initial elastic strain for the cracked structure problem</td></tr><tr><td>ε<sub>p</sub></td><td>initial plastic strain for the cracked structure problem</td></tr><tr><td>σ(0)</td><td>initial elastic stresses for the cracked structure problem</td></tr><tr><td>σ<sub>e</sub></td><td>elastic stress tensor</td></tr><tr><td>σ<sub>plastic</sub></td><td>initial plastic stresses for the cracked structure problem</td></tr><tr><td>ε<sub>e</sub></td><td>elastic strain tensor</td></tr><tr><td>da</td><td>crack growth rate</td></tr><tr><td>μ<sub>crack</sub></td><td>friction coefficient between the crack faces</td></tr><tr><td>μ<sub>wheel-rail</sub></td><td>wheel-rail friction coefficient</td></tr><tr><td>σ</td><td>Cauchy stress tensor</td></tr><tr><td>K</td><td>hook tensor</td></tr><tr><td>t</td><td>local traction field</td></tr><tr><td>u</td><td>global displacement field in the structure</td></tr><tr><td>w</td><td>local displacement field</td></tr><tr><td>θ</td><td>initial angle between the crack and the running surface</td></tr><tr><td>a</td><td>Hertzian contact semi-length</td></tr><tr><td>da</td><td>crack extension length</td></tr><tr><td>dN</td><td>cycle jump</td></tr><tr><td>k<sub>1</sub></td><td>mode I stress intensity factor at the tip of a virtual crack extension</td></tr><tr><td>k<sub>2</sub></td><td>mode II stress intensity factor at the tip of a virtual crack extension</td></tr><tr><td>K<sub>eq</sub></td><td>equivalent stress intensity factor at the actual crack tip</td></tr><tr><td>K<sub>II</sub></td><td>mode II stress intensity factor at the actual crack tip</td></tr><tr><td>K<sub>l</sub></td><td>mode I stress intensity factor at the actual crack tip</td></tr><tr><td>l</td><td>initial crack length</td></tr><tr><td>P<sub>max</sub></td><td>maximal Hertzian pressure</td></tr><tr><td>x<sub>c</sub></td><td>distance between the crack mouth and the maximal Hertzian pressure</td></tr></table>

and research institutes and universities (LaMCoS/INSA Lyon, LMS/Polytechnique, MECAMIX, IFSTTAR) within the consortium IDR2 (Initiative for Development and for Research on Rails) based on a theoretical, numerical and experimental approach. Concerning the numerical simulation, it starts from the train traffic to a fatigue life assessment [18].

* Step 1: a railway multi-body dynamics is performed, accounting for the actual rail and wheel geometries and load conditions, to determine contact conditions at the wheel-rail interface.

* Step 2: the stabilized cyclic mechanical state of the rail is calculated using a 3D elasto-plastic finite element simulation and an original, time-cost efficient direct stationary algorithm.

* Step 3: a fatigue analysis of the rail is performed according to Dang Van criterion to define the critical sites for crack initiation within the rail.

* Step 4: the crack (s) are considered and their behavior and growth are analyzed.

This paper focuses on this latter point. A two dimensional Extended Finite element model accounted for actual residual stresses and frictional contact between crack faces is combined with a mixed mode nonproportional crack propagation model previously developed. In Section 2, the theoretical basis of the two-scale XFEM strategy for cracked body problem with frictional contact at the interface accounting for actual residual stresses is presented. Then, in Section 3 the fatigue crack growth procedure is briefly introduced. Section 4 illustrates the reliability of the results. Parametric studies are then performed to analyze the influence of the tangential loading, the friction coefficient between the crack faces and the initial crack orientation on the crack growth path and rate. In the last part of this section, the actual residual stress influence is studied and the role of neighboring cracks in the crack propagation mechanism is emphasized.

## 2. Two-scale XFEM strategy for cracked body problem with frictional contact at the interface accounting for actual residual stresses

### 2.1. Two-scale strategy for cracked body problem with frictional contact at the interface

Different authors have proposed mixed weak formulation to take into account the contact and the friction between the crack faces (Fig. 1) within the X-FEM framework [19–22]. We here use the approach proposed in [20] and consider a cracked body $\Omega \in \mathbb{R}^3$. Contact and friction can occur along the crack faces $\Gamma^+$ and $\Gamma^-$ (Fig. 2(a)). Under small displacement and small strain assumptions, the interface $\Gamma$ ($\Gamma = \Gamma^+ \cup \Gamma^-$) is assumed as an autonomous entity with its own behavior possibly
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of the phenomena occurring at the crack scale.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Cracked body problem (a), divided into a global structure problem (b) and a local interface problem (c).</figcaption></figure>

nonlinear. This fracture problem is divided into a global ($\Omega \setminus \Gamma$ in Fig. 2(b)) and a local (crack scale in Fig. 2(c)) problem. The global problem is defined with its own primal and dual variables, **u** the displacement field and **σ** the Cauchy stress tensor respectively. The local problem is defined with its own primal and dual variables, **w** is the interface displacement field and **t** is the interface traction field respectively. Let **n** be the outward unit normal to $\Omega \setminus \Gamma$ and $\mathbf{n}_\Gamma$ is the outward unit normal to $\Gamma$.

We assume quasi-static formulation and write the governing equations as follows for the global problem at a given time $t \in [0; T]$:

$$ \mathbf{div} \, \boldsymbol{\sigma} + \mathbf{f} = 0 \quad \text{in} \quad \Omega \setminus \Gamma_c \qquad (1) $$

$$ \boldsymbol{\sigma}\mathbf{n} = \mathbf{s} \quad \text{on} \quad \partial_2 \Omega \qquad (2) $$

$$ \mathbf{u} = \mathbf{v} \quad \text{on} \quad \partial_1 \Omega \qquad (3) $$

where **s** and **v** are the prescribed load on $\partial_2 \Omega$ and displacement on $\partial_1 \Omega$ respectively, and **f** the body force vector. Both global and local problems obey a constitutive law possibly nonlinear in the bulk and on the interface respectively:

$$ \boldsymbol{\sigma} \equiv \boldsymbol{\sigma}(\mathbf{u}, u_i, t) \quad \text{in} \quad \Omega \setminus \Gamma_c \qquad (4) $$

$$ \mathbf{t} \equiv \mathbf{t}(\mathbf{w}, w_i, t) \quad \text{on} \quad \Gamma_c \qquad (5) $$

where $u_i$ and $w_i$ are the internal variables of the global and local problem respectively and $t$ a given time. Furthermore, we add the coupling between the global and the local problems both on the primal and dual variables:

$$ \mathbf{u} = \mathbf{w} \quad \text{on} \quad \Gamma_c \qquad (6) $$

$$ \boldsymbol{\sigma}\mathbf{n} = \mathbf{t} \quad \text{on} \quad \Gamma_c \qquad (7) $$

Eqs. ((1)-(3), (6), (7)) consist in the strong formulation of the crack model with contact and friction between the crack faces. As already shown in [23], an equivalent three field weak formulation (the global displacement field **u**, the local displacement field **w** and the local tractions field **t**) at a given time $t \in [0; T]$ can be introduced (see [20,24]).
$$0 = - \int_{\Omega} \boldsymbol{\sigma} : \boldsymbol{\epsilon}(\mathbf{u}^*) d\Omega + \int_{\Omega} \mathbf{f} \cdot \mathbf{u}^* d\Omega + \int_{\partial\Omega} \mathbf{s} \cdot \mathbf{u}^* dS + \int_{\Gamma_c} \mathbf{t} \cdot \mathbf{u}^* dS + \int_{\Gamma_c} (\mathbf{p} - \mathbf{t}) \cdot \mathbf{w}^* dS \\ + \int_{\Gamma_c} (\mathbf{u} - \mathbf{w}) \cdot \mathbf{t}^* dS \quad \forall \mathbf{u}^* \in \mathcal{U}_0, \forall \mathbf{w}^* \in \mathcal{W}, \forall \mathbf{t}^* \in \Lambda, \forall t \in [0; T] \quad (8)$$

with

$$\mathbf{u} \in \mathcal{U}, \mathcal{U} = \{\mathbf{u} \in H^1(\Omega \setminus \Gamma_c) / \mathbf{u} = \mathbf{v} \text{ on } \partial_1 \Omega\} \quad (9)$$

$$\mathbf{u}^* \in \mathcal{U}_0, \mathcal{U}_0 = \{\mathbf{u} \in H^1(\Omega \setminus \Gamma_c) / \mathbf{u} = 0 \text{ on } \partial_1 \Omega\} \quad (10)$$

$$\mathbf{w} \in \mathcal{W}, \mathbf{w}^* \in \mathcal{W}, \mathcal{W} = \{\mathbf{w} \in H^1(\Gamma_c)\} \quad (11)$$

$$\mathbf{t} \in \Lambda, \mathbf{t}^* \in \Lambda, \Lambda = \{\mathbf{t} \in L^2(\Gamma_c)\} \quad (12)$$

and is equivalent to the previous strong formulation (1, 2, 3, 6, 7) (see [23]). This weak formulation (8) allows a global-local X-FEM formulation handling successfully the structure scale (global), the crack and the interfacial nonlinear local scales. Here, an homogeneous, isotropic material with an elastic linear behavior is considered. The stress strain law is then written according to the fourth order Hooke tensor K:

$$\boldsymbol{\sigma} = \mathbf{K} \boldsymbol{\epsilon}(\mathbf{u}) \quad \text{in} \quad \Omega \setminus \Gamma_c \quad (13)$$

Furthermore, at the local scale, interface quantities (w, t) obey contact unilateral law and Coulomb's friction law for the normal and tangential problem respectively. w and t interface fields on the crack faces $\Gamma_c^+$ and $\Gamma_c^-$ are expressed in the local basis (n_c, t_c) relative to the crack

$$\mathbf{w} = w_N \mathbf{n}_c + w_T \mathbf{t}_c \quad \text{and} \quad \mathbf{t} = t_N \mathbf{n}_c + t_T \mathbf{t}_c \quad (14)$$

where ($w_N$, $t_N$) and ($w_T$, $t_T$) are normal and tangential scalar interface fields, respectively. Further relative interface displacements (opening and sliding) between the crack faces at a given interface point x are given by:

$$[w_N(\mathbf{x}, t)] = w_N^+(\mathbf{x}, t) - w_N^-(\mathbf{x}, t) \quad \text{and} \quad [w_T(\mathbf{x}, t)] = w_T^+(\mathbf{x}, t) - w_T^-(\mathbf{x}, t) \quad (15)$$

This yields for the local equations of the interface behavior:

$$\begin{aligned} \bullet \text{ opening } [w_N(\mathbf{x}, t)] > 0 & \rightarrow \mathbf{t}^+(\mathbf{x}, t) = \mathbf{t}^-(\mathbf{x}, t) = 0 \\ \bullet \text{ contact } [w_N(\mathbf{x}, t)] = 0 & \rightarrow \mathbf{t}^+(\mathbf{x}, t) = -\mathbf{t}^-(\mathbf{x}, t) \\ \bullet \text{ sticking } |t_T(\mathbf{x}, t)| < \mu_c |t_N(\mathbf{x}, t)| & \rightarrow \Delta[w_T(\mathbf{x}, t)] = 0 \\ \bullet \text{ sliding } |t_T(\mathbf{x}, t)| = \mu_c |t_N(\mathbf{x}, t)| & \rightarrow \exists \gamma > 0 / \Delta[w_T(\mathbf{x}, t)] = -\gamma t_T^+(\mathbf{x}, t) \end{aligned} \quad (16)$$

where $[w_N(\mathbf{x}, t)]$ and $[w_T(\mathbf{x}, t)]$ are the normal and tangential jump at the interface, $\mu_c$ is the friction coefficient at the interface

## 2.2. X-FEM discretization and nonlinear solver dedicated to interfacial frictional contact

The eXtended Finite Element Method [25] enriches the approximation space so that it allows to account for discontinuities like cracks. It addresses nonlinear problems, no explicit crack mesh and no remeshing during the crack growth process. Compared to the FEM, X-FEM consists in introducing additional functions to the discretization basis following the partition of unity framework [26]. In this respect, due to the presence of a crack, one defines an appropriate approximation displacement field composed of three terms: (i) the first term or the standard finite element approximation, characterizing the continuum behavior of the body, (ii) the second term corresponding to an arbitrary displacement discontinuity along the crack faces, (iii) and the latter capturing the asymptotic behavior of the stresses near the crack tip. One assumes that the state vector $\mathbf{X}_n = (\mathbf{u}_n, \boldsymbol{\sigma}_n, \mathbf{w}_n, \mathbf{t}_n)$ is known at time $t_n$. Within a quasi-static incremental framework, the next stage consists of calculating the unknown state vector $\mathbf{X}_{n+1}$ at time step $t_{n+1}$. The LATIN nonlinear solver consists of an iterative strategy between a local stage from $\mathbf{X}_n^i$ to $\mathbf{X}_n^{i-1/2}$ and a global stage from $\mathbf{X}_n^{i+1/2}$ to $\mathbf{X}_n^{i+1}$. The local stage corresponds to a set of local equations, possibly nonlinear (Eq. 16), and the global stage to a set of global linear equations (Eq. 8). This two-steps approach requires search directions $E^{G-L}$ and $E^{L-G}$ between the set of equations (L) and (G). The process is repeated until convergence is reached.

The introduction of the X-FEM discretization of the fields and the conditions introduced by the stabilization process [24] lead to the discretized linear system (25). A priori formula were recently proposed to set up the parameters controlling the convergence rate (search direction and stabilization parameter) to obtain the best convergence rate possible [27] and improve the method efficiency by reducing the computing time.

$$\begin{bmatrix} \mathbf{K} & 0 & -\mathbf{K}_{ut} \\ 0 & \mathbf{K}_{ww} & \mathbf{K}_{wt} \\ -\mathbf{K}_{ut}^T & \mathbf{K}_{wt}^T & \mathbf{K}_{ft} \end{bmatrix} \begin{pmatrix} \mathbf{U}_{i+1} \\ \mathbf{W}_{i+1} \\ \mathbf{T}_{i+1} \end{pmatrix} = \begin{pmatrix} \mathbf{F} \\ \mathbf{K}_{wt} \cdot \mathbf{T}_{i+1/2} + \mathbf{K}_{ww} \cdot \mathbf{W}_{i+1/2} \\ \mathbf{K}_{ft} \cdot \mathbf{T}_{i+1/2} \end{pmatrix} \quad (17)$$
The evaluation of the mechanical state in the rail is crucial for the prediction of its durability. The repeated contact stresses at the wheel/rail interface during rolling motion are responsible for plastic deformations within the contact zone. The residual stress field in the rail section were studied in a lot of works [10–15,28,29]. Nevertheless when dealing with crack propagation simulation, residual stresses are basically introduced in the model applying forces [8,16,17] at the boundaries of the considered domain. The current study focuses on the analysis of the influence of actual computed residual stresses on crack propagation.

## 2.3. Introduction of actual plastic stresses in the X-FEM strategy

The residual stresses considered here correspond to the plastic strain accumulation due to the repeated passages of the wheels on the rails and to the manufacturing process [10]. They result from the rail cyclic stress and strain computation under the action of a great number of these contact loads. The difficulty lays on the fact that plastic deformations and high stress gradients can occur in the area near the rail-wheel contact zone due to the high level of the rolling-sliding contact loads. Therefore, very refined finite element meshes are used to capture accurately the stress gradients. Simulating a great number of passages by running classical finite element simulations based on classical integration schemes would be very lengthy and require considerable computing ressources. To avoid these drawbacks, alternative methods were proposed like the direct stationary method [29–31], going straight to the steady state in the case of a repeatedly moving load and applied to rail modeling. It is based on the following ideas:

* the simulation is performed in a reference moving with the load, i.e., the wheel reference, that allows to take advantage of the time independance of mechanical quantities;

* the regime is supposed to be stationary in this reference;

* the rail material flow is modeled using Eulerian formulations;

* elastoplastic behavior such as elastic or plastic shakedown can be considered;

* refined meshing is only necessary in the vicinity of the contact zone;

This alternative method allows proposing a very time cost-efficient computing algorithm and easy to use. An elastoplastic constitutive law with a linear kinematic hardening is used for the steel of the rail.

Using this algorithm, two kinds of asymptotic states may be obtained. The former is elastic shakedown (Fig. 3(a)) where the limit cycle is elastic or plastic shakedown (Fig. 3(b)) with an elasto-plastic limit cycle. The latter is the ratcheting for which plastic strains increase continuously without reaching an asymptotic state and thus a limit cycle. The attention is focused here on elastic shakedown. Plastic strains are plotted in Fig. 4(a)–(c) for elastic and plastic shakedown respectively. Elastic shakedown is characterized by all components of the plastic deformation tensor constant along all the streamlines, here the horizontal lines of the mesh whereas it is not the case for plastic shakedown.

The problem is expressed as follows. In the initial state ($t = -\infty$), the rail is considered as a virgin structure. At $t = 0$, it is assumed that the rail is in the asymptotic state and that shakedown has occured before crack initiation. We then perform fatigue crack growth simulation (see Fig. 5).

For a given time $t$, the equilibrium equations, the constitutive law and the boundary equations are summarized hereafter:

Equilibrium: $div \sigma_e(t) = 0$ in $\Omega$

Constitutive law: $\sigma_e(t) = K\epsilon_e(u(t))$ in $\Omega$

Dirichlet condition: $\sigma_e(t)n = F(t)$ on $\partial_2 \Omega$

Neumann condition: $u(t) = u_d(t)$ on $\partial_1 \Omega$

(18)

<figure><img src="image_4.png" /><figcaption>Fig. 3. Rail elastoplastic behavior: asymptotic state.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. $\epsilon_{\infty}$ plastic under loading in the case of elastic shakedown (b) and plastic shakedown (c) in the area of interest (a).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Notation for the introduction of residual stresses.</figcaption></figure>

with

$$ \sigma_e(t) = \sigma(0) + \Delta\sigma(t) \qquad (19) $$

$$ \epsilon(t) = \epsilon(0) + \Delta\epsilon(t) \qquad (20) $$

where $\sigma(0)$ are the elastic stresses related to $\epsilon(0)$ (Fig. 6). If no loading is assumed, we only have plastic strains $\epsilon_p$ and the residual stresses are the plastic stesses defined by $\sigma_{plastic} = K\epsilon_p$

Moreover

$$ \epsilon(t) = \epsilon_p + \epsilon_e(t) \qquad (21) $$

Multiplying (21) by $K$ the left and right hand sides and introducing (20) in the left hand side, we obtain:

$$ K(\epsilon(0) + \Delta\epsilon(t)) = K\epsilon_p + \sigma(t) \qquad (22) $$

<figure><img src="image_7.png" /><figcaption>Fig. 6. Output of the stationary algorithm.</figcaption></figure>
We now introduce (22) into (19) to express $\Delta\sigma(t)$ as:

$$ \Delta\sigma(t) = \mathbf{K}(\Delta\epsilon(t) - \epsilon_p) + \mathbf{K}\epsilon(0) - \sigma(0) \quad (23) $$

which is equivalent to:

$$ \mathbf{K}\Delta\epsilon(t) - \Delta\sigma(t) = \mathbf{K}\epsilon_p - (\mathbf{K}\epsilon(0) - \sigma(0)) \quad (24) $$

This equation is introduced into the Eq. (8). The linear system (Eq. (25)) is obtained with $\Delta\mathbf{U}_{i+1}$ instead of $\mathbf{U}_{i+1}$ as unknown vector. Before post-processing the results to compute the stress intensity factors, an update step should be performed in the same way than in Eqs. (19) and (20) to obtain the field $\mathbf{U}$.

$$ \begin{bmatrix} \mathbf{K} & 0 & -\mathbf{K}_{ut} \\ 0 & \mathbf{K}_{ww} & \mathbf{K}_{wt} \\ -\mathbf{K}_{ut}^{T} & \mathbf{K}_{wt}^{T} & \mathbf{K}_{tt} \end{bmatrix} \begin{pmatrix} \Delta\mathbf{U}_{i+1} \\ \mathbf{W}_{i+1} \\ \mathbf{T}_{i+1} \end{pmatrix} = \begin{pmatrix} \mathbf{F} - \mathbf{F}(0) \\ \mathbf{K}_{wt} \cdot \mathbf{T}_{i+1/2} + \mathbf{K}_{ww} \cdot \mathbf{W}_{i+1/2} \\ \mathbf{K}_{tt} \cdot \mathbf{T}_{i+1/2} \end{pmatrix} \quad (25) $$

with:

$$ \mathbf{F}(0) = \langle \mathbf{F}_l(0) \rangle \quad (26) $$

$$ \mathbf{F}_l(0) = \int_{\Omega} \mathrm{Tr}[(\mathbf{K}\epsilon_p - (\mathbf{K}\epsilon(0) - \sigma(0)))\epsilon(\phi_l)]d\Omega \quad (27) $$

The internal generalized force $\mathbf{F}_l(0)$ has to be added to the right hand side of the linear system. Its estimation is based on fields computed on the stationary method mesh. The method to compute the residual stresses assumes no gradient in the moving load direction. The mesh is hence designed by a domain translation according to this direction. Different requirements are needed for the mesh created for the crack propagation simulation. Indeed, a fine mesh under the moving wheel-rail contact is needed to capture the stress and strain gradients, but an even finer mesh in the crack propagation area is necessary to determine accurately the K-dominance zone and the stress intensity factors at the crack tip. To reduce the number of degrees of freedom and optimize the computing resources, coarse and fine mesh domains are designed versus either the residual stress computation or the crack propagation simulation. Therefore residual fields estimated according to the stationary method with a coarser mesh are transferred onto the finer one defined for the propagation.

## 2.4. Projection step to introduce the generalized forces on the discretization used for the crack propagation simulation

The third step is concerned with the projection of those residual fields onto the mesh specially designed for the crack propagation. Two main categories of methods exist to project mechanical fields between meshes with different refinement. Very briefly, either the fields are projected from mesh 1 to mesh 2 and a balance step is introduced to ensure the respect of the equilibrium on mesh 2 or only part of the fields are projected on mesh 2, computed according to the constitutive law and the local variables are updated. This is this method which is used in this work. Under the assumptions of elastic shakedown and considering an unloaded body we therefore have $\mathbf{u}(0) = \epsilon(0) = \sigma(0) = \mathbf{0}$ and we only consider the plastic stresses:

$$ \mathbf{F}_l(0) = \int_{\Omega} \mathrm{Tr}[(\mathbf{K}\epsilon_p)\epsilon(\phi_l)]d\Omega \quad (28) $$

The plastic stresses are defined by $\mathbf{K}\epsilon_p = \sigma_{plastic}$ and

$$ \begin{pmatrix} \Delta\mathbf{U}_{i+1} \\ \mathbf{W}_{i+1} \\ \mathbf{T}_{i+1} \end{pmatrix} = \begin{pmatrix} \mathbf{U}_{i+1} \\ \mathbf{W}_{i+1} \\ \mathbf{T}_{i+1} \end{pmatrix} \quad (29) $$

The plastic stress field projection onto the finer mesh followed by the computation of the equivalent generalized forces on this mesh allows to enforce the equilibrium of the structure afterwards. Initially, the field $\sigma_{plastic}$ is defined at the Gauss points of the coarse discretization. The projection is performed in two steps. The element field is first computed on the nodes of the coarse mesh elements using an interpolation based on the shape functions and a least square to compute the value at a given node. Then the projection is performed element by element from the coarse mesh $m_1$ to the fine mesh $m_2$ using the linear shape functions as interpolation functions. The interpolations could have been performed using different functions than the shape functions. This choice has been done for practical reasons. The last step consists in computing the generalized forces corresponding to the projected plastic stresses on all the nodes of the mesh $m_2$ (Fig. 7).

As the material is assumed to behave elastically when dealing with the crack propagation simulation, an initial elastic resolution is performed accounting only for generalized forces corresponding to the projected plastic stresses. The constitutive law is automatically verified and the stress fields are balanced on the fine discretization. The asymptotic mechanical fields are then projected on the fine mesh designed for crack propagation analysis (see Fig. 8). Those fields are the initial state of the propagation simulation. This state is considered as permanent and nonuniform. No redistribution of plastic stresses is performed throughout the crack growth.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Projection procedure of the plastic stresses from a mesh $m_1$ to a mesh $m_2$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. 2D projection of the asymptotic mechanical fields computed on mesh $m_1$ on the mesh used for the crack propagation simulation (mesh $m_2$).</figcaption></figure>

### 3. Fatigue crack growth procedure

Rolling fatigue cracks experience nonproportional mixed mode loading during the load traveling corresponding to a characteristic sequence of tensile stresses followed by a shear cycle, with complex sequences of opening/contact and stick/slip distributions along the crack faces. These distributions govern the stress intensity factors and the growth behavior, both rate and path. Some studies [32] on nonproportional loading have shown the inability of classical criteria like the Maximum Tangential Stress MTS, the Maximum Shear Stress MSS, the minimum strain energy density $S_{min}$, the maximum energy release rate $\dot{G}_{max}$, etc. to determine the crack propagation direction under out of phase loading histories. The crack growth path $\theta$ direction and rate $da/dN$ are here predicted according to Hourlier and Pineau's criterion [33]. This criterion was successfully used previously under both rolling [8] and fretting fatigue [34–36]. The crack path predictions were in very good agreement with experimental observations.

The mixed-mode fatigue crack propagation model is briefly recalled here. A more detailed explanation can be found in [8] and more recently in [36]. The first step consists in computing the stress intensity factors (SIFs) histories at crack tip during the load traveling over the rail (one wheel passage) defined as one load cycle. This load cycle is divided into time steps corresponding to successive wheel positions with respect to the crack location. The crack body problem is solved at each time step using the two scale strategy (2.1). The SIFs are computed using integral methods.

The second step is concerned with the crack growth path and rate determination. The crack growth path is predicted according to Hourlier and Pineau's criterion [36]. It assumes that a crack follows the easiest path available that is the path along the direction $\beta$ that maximizes the growth rate at the tip of an infinitesimal virtual part $s$ of the kinked crack. Values of $k_1$ and $k_2$, the local stress intensity factors at the tip of the virtual infinitesimal part $s$ inclined at an angle $\beta$ of the original
<figure><img src="image_10.png" /><figcaption>Fig. 9. Parameters introduced in the Hourlier and Pineau's criterion to compute the propagation angle $\beta$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Flowchart of the propagation procedure.</figcaption></figure>

crack, are calculated from $K_I$ and $K_{II}$ versus the angle $\beta$ (Fig. 9). Corresponding $\Delta k_1$ and $\Delta k_2$ range are then calculated for each angle $\beta$. A dedicated mixed-mode propagation law is then used to predict the crack growth rate [8], based on an equivalent stress intensity factor range.

$$ \frac{da}{dN} = 2.10^{-9} (\Delta K_{eq})^{3.33} \qquad (30) $$

with

$$ \Delta K_{eq} = (\Delta K_{I}^{2} + 0.772 * \Delta K_{II}^{2})^{0.5} \qquad (31) $$

A crack extension is added to the current crack and the corresponding increment of number of cycles is estimated according to the crack propagation law. This procedure is repeated until either a critical crack size or a given number of cycles is reached (see Fig. 10).

## 4. Validation: X-FEM modeling of a surface breaking crack under a moving Hertzian loading

In this section, results obtained according to the X-FEM model previously described here are compared with results from the literature. The influence of different factors on the SIFs $K_I$ and $K_{II}$ and the crack growth behavior are analyzed. Factors like the wheel-rail friction coefficient, the crack face friction coefficient, the crack inclination with respect to the free surface, the plastic stresses and the crack network interactions are investigated.

Sets of wheel/rail contact forces, wheel diameters and crack morphologies used here for the comparison are those defined in [8]. The fatigue crack model is based on the crack analogy with continuous distributions of dislocations combined with unilateral contact analysis. Detailed calculations of mode I and mode II effective SIFs ranges under nonproportional sequential loading cycles can be performed at crack tips, for any crack geometry (straight or kinked), crack number, accounting for interactions between cracks [37,38]. A 6 mm long crack inclined at 15° with respect to a half-plane (the rail) is considered.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Definition of the parameters.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Multi-scale mesh to capture the different scales of the problem.</figcaption></figure>

The crack friction coefficient is assumed to be equal to 0.1. A fully sliding Hertzian contact load is considered with peak pressure, contact length and friction coefficient equal respectively to P<sub>max</sub> = 845 MPa, 2a = 13.5 mm, μ<sub>wheel-rail</sub> = 0.4 (see Fig. 11). A rectangular domain, with dimensions L = 300 mm times H = 172 mm, is considered as a half-plane. The finite element mesh used is presented in the Fig. 12. The refined zone is a rectangular domain containing the crack which size is automatically adapted. For the reference case and in the initial configuration, the dimensions are l = 10 mm and h = 4 mm. The ratio between the biggest and the smallest elements is higher than 500. The fine mesh domain is obviously located under the contact zone to capture accurately the stress and strain gradients. Further, in the area where the crack will propagate, the mesh is even finer to ensure a reliable modeling of the K-dominance zone. The SIF obtained according to both models are presented Fig. 13 versus x<sub>c</sub>/a. The leading edge of the load distribution is located at the surface breaking tip of the crack at x<sub>c</sub>/a = -1. Note that in [8] the SIF computation is based on the displacement jump between the crack faces whereas an energetic approach and domain integral at the crack tip vicinity are used here. Therefore, different SIFs are obtained according to both

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison of (a) K<sub>I</sub> and (b) K<sub>II</sub> values obtained using CASTEM and [8].</figcaption></figure>
methods in the particular case of a crack under full sticking conditions during several steps of the moving load. Constant $K_{II}$ values are obtained according to [8] whereas slightly varying ones are obtained here due the stress variations during the load traveling. This effect can be observed in the Fig. 13 for the maximun $K_{II}$ values using our model (CAST3M [39] in Fig. 13). Despite those differences the results are in very good agreement.

# 5. Parametric study

The influence of different parameters on crack growth path and rate are investigated. Parameters corresponding to the so-called reference case are listed in Table 1.

## 5.1. Influence of the friction coefficient between the wheel and the rail

The crack inclination is imposed as well as the direction of the moving load, from negative to positive x values. The influence of the friction coefficient between the wheel and the rail is analyzed and the effects of both positive and negative values are considered. The convention used is defined in Fig. 11 as well as the direction of the moving load. Values of 0.025, 0.1 and 0.4 are considered. For $x_c/a < -1$, the load distribution has not moved over the crack whereas for $x_c/a > 1$, the load is located over the crack before traveling away. The tangential load and the normal load effects on the relative sliding and opening along the crack faces can either be added or subtracted depending on the sign of the friction coefficient and on the position of the load with respect to the crack. This is illustrated in Figs. 14(b) and 15(b) for $K_{II}$ and Figs. 14(a) and 15(a) for $K_I$.

The numerical crack propagation is then analyzed. The procedure described in Fig. 10 is used. The crack path is determined at each numerical cycle, the crack growth is simulated with a 4% crack total length increment and the corresponding increment of number of cycles is quantified. 25 numerical load cycles are simulated, except if a crack self-arrests or branches upwards out of the domain. For the crack propagation, a higher friction coefficient between the wheel and the rail causes an earlier branching. $\Delta K_{II}$ is increased and the crack branches. Depending on the orientation of the tangential load, the cracks branch earlier upwards to the rail surface (Fig. 16(a)) or downwards in the rail web (Fig. 16(b)).

## 5.2. Influence of the crack friction coefficient

The influence of the friction coefficient between the crack faces is studied, all the other parameters being unchanged (see Table 1). An increasing in the crack friction coefficient decreases the relative sliding between the crack faces and consequently reduces $\Delta K_{II}$ (see Fig. 17(b)) whereas it has obviously no influence on the crack opening and on $\Delta K_I$ (see Fig. 17(a)). The crack friction coefficient and the wheel-rail friction coefficient have thus opposite influences on the crack growth rate [8,17].

<table><caption>**Table 1**
Reference parameters for the parametric study.</caption><thead><tr><td>l (mm)</td><td>θ (°)</td><td>μ<sub>crack</sub></td><td>μ<sub>wheel-rail</sub></td><td>P<sub>max</sub> (MPa)</td><td>2a (mm)</td></tr></thead><tbody><tr><td>6</td><td>15</td><td>0.5</td><td>0.025</td><td>845</td><td>13.5</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** (a) $K_I$ and (b) $K_{II}$ computed for different positive friction coefficients $\mu_{wheel-rail} > 0$ (braking situation) ($l = 6$ mm, $\theta = 15^{\circ}$, $\mu_{crack} = 0.5$).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15. (a) $K_I$ and (b) $K_{II}$ computed for different negative friction coefficients $\mu_{wheel-rail} < 0$ (speeding up condition) ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.5$).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Crack propagation path for different friction coefficient between the wheel and the rail in braking situation (a) and speeding up condition (b).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. (a) $K_I$ and (b) $K_{II}$ computed for different crack friction coefficient $\mu_{crack}$ ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{wheel-rail} = 0.025$).</figcaption></figure>
The SIF range $\Delta K_{II}$ is thus dependent on the value of the crack friction coefficient. The crack growth rate and direction are also modified as shown in Fig. 18. The crack propagation is divided into two periods. During the first one, a co-planar crack extension is obtained whatever the value of $\mu_{crack}$. Then, differences in the crack path are predicted. Increasing $\mu_{crack}$ value causes (i) an earlier crack branching and (ii) the crack branches downwards instead of upwards to the upper surface. This is more dangerous as the former may lead to rail failure whereas the latter may simply lead to a spall.

### 5.3. Influence of the initial crack orientation

The influence of the crack inclination on the crack growth is analyzed in this section. The loading conditions and the crack parameters are defined in Table 1. Different values of $\theta$ are considered (15°, 30°, 45°, 60° and 75°) from shallow angled to nearly perpendicular crack to the upper surface. $K_I$ and $K_{II}$ variations at the $\theta$ angled crack tip are plotted versus the load location during a cycle (see Fig. 19(a) and (b)). A shallow inclination reduces the relative slip between the crack faces but eases the crack opening, for the configurations considered. Hence, $\Delta K_{II}$ increases with increasing $\theta$ while $\Delta K_I$ decreases. The crack growth driving force under mixed mode loading being a shear-tensile mechanism, the crack rate and branching direction are influenced. The crack branches downwards whatever $\theta$, except for $\theta = 75°$ for which a co-planar extension is predicted (see Fig. 20).

### 5.4. Influence of the plastic stresses

The influence of plastic stresses on the crack growth behavior is now investigated. The procedure described in Section 2.4 is used. This plastic stress field has been computed on a different mesh using a dedicated numerical tool briefly introduced in Section 2.3. The plastic field introduced is a combination of the residual stresses coming from the manufacturing process and the ones deriving from the plastic strain accumulation due to the train traffic.

#### 5.4.1. Plastic stress field tensor

The components of plastic stresses $\sigma_{plastic_{xx}}$, $\sigma_{plastic_{yy}}$, $\sigma_{plastic_{xy}}$ and $\sigma_{plastic_{zz}}$, in the center of the considered domain, are plotted in Fig. 21. This figure shows that a complex stress state exists within the rail subsurface, with successive layers of different thicknesses and different magnitude of tensile and compressive stresses: first, tensile stresses (about 25 MPa) exist at the surface of the rail, followed just below by high compressive stresses (about -110 MPa) within a 1 or 2 mm thick layer of

<figure><img src="image_19.png" /><figcaption>Fig. 18. Crack propagation path for different local friction coefficient $\mu_{crack}$ ($l = 6$ mm, $\theta = 15°$, $\mu_{wheel-rail} = 0.025$).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. (a) $K_I$ and (b) $K_{II}$ computed different initial crack orientation $\theta$ ($l = 6$ mm, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.025$).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 20. Crack propagation path versus the initial crack orientation $\theta$ ($l = 6$ mm, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.025$).</figcaption></figure>

the rail subsurface, then tensile stresses again (50 MPa) and finally again compressive stresses (maximum –50 MPa) for 10 mm depth. We also notice that shear stresses are localized just beneath the running surface and vanish in the rail web. This plastic stress distribution results from the accumulation of plastic strain caused by repeated wheel rail contact in the near surface vicinity and deeper in the rail subsurface from the manufacturing process, characterized by the C-shape of residual stress field. The $\sigma_{xx}$ stress component is shown after projection on the mesh used for the simulation of the crack propagation and the initial elastic re-balanced step (Fig. 21). As already mentioned, this initial state is considered as a non-uniform and permanent state for the fracture problem.

### 5.4.2. Influence of the plastic stress field

SIF and the crack growth behavior are now analyzed and compared with and without accounting for the plastic stresses, for the reference case (data in Table 1). The SIF are plotted in Fig. 22(a) and (b) for mode I and II respectively. As the plastic stresses vary from tensile to compressive stresses from the crack mouth to the tip, they inhibit the crack opening during the wheel traveling and $K_I = 0$. $\Delta K_{II}$ is further reduced as the compressive stresses decrease the relative sliding between the crack faces.

Concerning the crack growth, the plastic stresses have a significant influence. First, the crack rate is reduced and secondly, after a longer co-planar extension, the crack branches slightly upwards to the free surface towards the zone submitted to tensile stresses (Fig. 23), instead of branching downwards. Two major consequences are highlighted through this example: the plastic stresses may increase the fatigue life and may be responsible for a change in the damage mechanism. From failure to spall for the considered case. Accounting for actual plastic stresses is thus of utmost importance for the fatigue prediction.

The next step of this parametric study is concerned with multiple cracks and their interactions. In rails, cracks developed often in networks, either on the rail head or at the gauge corner [40]. It is important to quantify the modifications in the SIFs

<figure><img src="image_22.png" /><figcaption>Fig. 21. Plastic stress field after projection on the mesh used for the propagation simulation on a vertical line at the center of the domain.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 22. (a) $K_i$ and (b) $K_{ii}$ computed with the plastic stresses ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.025$).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 23. Crack propagation path with and without plastic stresses ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.025$).</figcaption></figure>

due to crack interactions with respect to the values obtained for a single crack and to expand this analysis to crack propagation.

## 5.5. Crack network

The crack growth algorithm is modified in case of multiple cracks. The growth increment is imposed equal to $\Delta a_{max} = 4\%$ of the current crack length in the single crack analysis. In the multiple crack analysis, the maximum growth increment is still imposed equal to 4% and the crack growth increment for each crack $i$ $\Delta a_i$ is determined proportionally to the equivalent stress intensity factor range to the power of the Paris law, such as:

<figure><img src="image_25.png" /><figcaption>Fig. 24. (a) $K_i$ and (b) $K_{ii}$ computed for a network of 2 cracks in a squat configuration ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.01$, $\mu_{wheel-rail} = 0.025$). The reference case corresponds to same configuration with only a single crack.</figcaption></figure>
$$ \frac{\Delta a_i}{\Delta a_{max}} = \left( \frac{\Delta K_{eq_i}}{\Delta K_{eq_{max}}} \right)^{3.33} \qquad (32) $$

with

$$ \Delta K_{eq} = (\Delta K_{I}^{2} + 0.772 * \Delta K_{II}^{2})^{0.5} \qquad (33) $$

<figure><img src="image_26.png" /><figcaption>Fig. 25. Crack propagation path for 2 cracks in a squat configuration ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.01$, $\mu_{wheel-rail} = 0.025$).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 26. Schematic representation of the initial configuration with 9 parallel cracks ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.1$).</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 27. (a) $K_I$ and (b) $K_{II}$ computed for a network of 9 parallel cracks ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.1$). Reference stands for the SIFs obtained for a single crack in the same configuration.</figcaption></figure>
### 5.5.1. Two symmetrical cracks

The first case is concerned with two cracks labeled 1 and 2 (see Fig. 25). They have the same length 6 mm, their initial inclination is $\theta = 15^\circ$, the crack friction coefficient is $\mu_{crack1} = \mu_{crack2} = 0.5$. The spacing is equal to 1 mm. Crack 1 surface breaking tip is located at the origin of the reference axis. The load position is defined by $x_c/a$. The loading conditions are those defined in Table 1. $K_I$ and $K_{II}$ variations versus the load location are presented in Fig. 24 for cracks 1 and 2, accounting for crack interactions. Two reference curves are added, corresponding to the SIF variations computed at crack 1 tip and crack 2 tip in the single crack configuration.

The effects of the interactions can thus be quantified. Differences in mode II SIF are not significant. The major difference concerns the crack opening of both cracks and thus $\Delta K_I$. Crack 1 opening is enhanced by the relative displacements along crack 2 for $x_c/a < -2$ while crack 2 opening is enhanced by the relative displacements along crack 1 for $x_c/a > 2$.

The crack growth is then predicted. The results are presented in Fig. 25. Crack 1 branches and extends parallel to the rail surface whereas crack 2 experiences a co-planar growth, as the one obtained in the one-crack configuration (Fig. 18). After several propagations steps, $\Delta K_{eq1}$ is decreasing causing a reduction in the growth rate. These predictions are coherent with observations made in [41–43] for a similar crack configuration: the second crack propagates deeper in the rail than the first one.

### 5.5.2. Nine parallel cracks

The test case is presented in the Fig. 26. 9 initial parallel cracks, labeled from 1 to 9, of the same length ($l = 6$ mm, $\theta = 15^\circ$) with equal spacing of 3.375 mm (equal to $a/2$, a quarter of the Hertzian contact zone) are considered. The crack coefficient of friction $\mu_{crack} = 0.5$ is the same for all cracks. $\mu_{wheel-rail} = 0.1$ (braking situation). Cracks labeled 3 to 6 are located at the center of the network. $K_I$ and $K_{II}$ are plotted in Figs. 18(b) and 27(b) respectively, for the nine cracks. The reference SIF variations in mode I and II, computed for the one-crack configuration, are also plotted for comparison and estimation of the interaction effects. The influence of the interactions between cracks are visible. Note that the crack inclination combined with the load direction and the sign of the tangential load favors crack opening, enlarged for all cracks when compared to the one-crack reference case, as the domain compliance is increased. A maximum $\Delta K_I$ increase of 150% is obtained. Concerning the relative sliding between the crack faces, the load sharing between the cracks leads to a shielding effect. $\Delta K_{II}$ is thus significantly reduced (between 0 and 25 % depending on the crack).

The computations for the fatigue crack network propagation are then carried out. The fatigue crack growth is co-planar before branching downward in the case of a single crack (Fig. 18(a)). At each load cycle, a growth increment is added to each crack with respect to the value of the corresponding equivalent SIF range along the direction determined angle $\theta_i$. The calculation is rapidly stopped as the crack merging is not analyzed in this paper.

<figure><img src="image_29.png" /><figcaption>Fig. 28. Crack path after 5 propagation steps for a network of 9 parallel cracks $\theta$ ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.1$).</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 29. (a) $K_I$ and (b) $K_{II}$ computed for a network of 9 parallel cracks with plastic stresses ($l = 6$ mm, $\theta = 15^\circ$, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.1$).</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. 30.** Propagation path for a network of 9 parallel cracks without and with plastic stresses $\theta$ ($l = 6$ mm, $\theta = 15^{\circ}$, $\mu_{crack} = 0.5$, $\mu_{wheel-rail} = 0.1$).</figcaption></figure>

The crack path is represented for each crack of the network in Fig. 28. No significant change is obtained for cracks 1, 2, 8 and 9 whereas cracks 3, 4, 5 and 7 branch parallel to the surface. The effects of the crack interactions are not the same for all the cracks and are very complex. They depend for each crack on both the equivalent stress intensity factor range that governs the crack rate and on the maximum equivalent stress intensity factor range at the tip of an infinitesimal segment that determines the growth angle. It is interesting to notice that a shape of spalls is predicted for a few crack and propagation in the depth of the rail for others.

### 5.5.3. Introduction of plastic stresses in the crack network

Finally, the simulation is carried out again but accounting for plastic stresses. Mode I and mode II SIFs are plotted in Fig. 29(a) and (b) respectively. The compressive plastic stresses extend up to the crack tips and inhibit crack opening for all cracks, except for crack 9 due to its particular location at the edge of the crack network. The relative sliding along all the crack faces is also reduced and a significant decrease in mode II stress intensity factor range is noted for all cracks.

The crack network propagation is also different. All cracks branch and extend in parallel to the free surface in the direction of the moving load (see Fig. 30). Crack 9 behaves differently, being less influenced by the crack interactions.

## 6. Conclusion

This paper aims at predicting fatigue crack growth and branching conditions under rolling contact fatigue taking into account actual plastic stresses and not uniform ones. The model we present, combining the X-FEM discretization and the LATIN nonlinear solver, is a robust tool capable of taking into account contact and friction at the interface and to simulate efficiently fatigue crack growth. Furthermore, the plastic fields are computed with a dedicated tool based on a direct stationary algorithm. Elastic shakedown is likely to occur for the stabilized rail mechanical state. Using this assumption and thanks to the LATIN nonlinear solver, the introduction of plastic stresses is performed easily by only adding a generalized internal force term in the right hand side of the equation describing the cracked body equilibrium. The developed numerical tool within CAST3M [39] is formed by the proposed two-scale strategy taking into account the fields acting at the rail scale (plastic stresses), the contact and friction between the crack faces at the crack scale, the fatigue Hourlier and Pineau's crack branching criterion and Paris type fatigue crack growth laws. The X-FEM discretization of the structure problem allows very easy propagation parametric studies since the crack mesh is not required to be conform with the structure one. Furthermore, it is very interesting to keep the same mesh for a parametric study. Indeed, all the simulations can then been compared at a fixed level of accuracy. The role of liquid is here assumed to simply reduce friction between the crack faces and to enhance the relative movements.

The crack growth is influenced by multiple parameters, in particular by plastic stresses, contact and interfacial crack frictional effects and neighboring cracks. Hence, the variation of one of the mentioned parameters results in changes of the crack propagation path. These parametric studies also give us a better understanding of the major damage observed under various contact conditions. We have seen that the friction coefficient between the crack faces and the wheel and the rail have an antagonist effect on the growth rate. This point has been largely pointed out in the literature [5,7,8]. But in this paper propagation simulation are led under all those conditions. This emphasizes an earlier crack branching when the friction coefficient between the wheel and the rail is high. Depending on the orientation of the tangential loading, rolling contact fatigue cracks are expected to propagate upwards for accelerating train situations and downwards when the train is braking. The work on the initial crack orientation has demonstrated that once cracks have branched downwards, for a given set of parameters, a single crack follows the same propagation direction. This paper also demonstrates the influence of plastic stresses on the propagation path for fatigue cracks in rail. In this work we chose to place the crack tip in the zone of high compressive stresses which lead to a decrease of the SIFs. The consequence was a longer coplanar propagation in the horizontal strip of compressive stresses.

Those results rely on the study of a single crack in a 2D domain. The study of the squat configuration seems to show good agreements with observations that were reported in some works [41,42]. The crack network investigation showed that crack propagation behavior are modified and for the considering case that few cracks branch earlier in the horizontal direction whereas others propagate downwards in the rail. A shielding effect has been observed for the $K_{II}$ values of the center cracks. It is interesting to highlight that they are subjected to higher opening forces when dealing with a parallel crack network due to the orientation of the tangential load and the crack initial inclination. The introduction of plastic stresses in the crack network illustrates once again the high role of the plastic stress fields in the fatigue crack propagation behavior.
Further improvements of the modeling include the rail bending or thermal loading and 3D crack propagation modeling.
The quantitative results require to be compared with the measurements coming from the railway network. This task implies
a statistical work to be able to compare the results of this deterministic model with the variations of very local effects in
practical situations.

References

[1] Bogdanski S, Lewicki P, Szymaniak M. Experimental and theoretical investigation of the phenomenon of filling the rcf crack with liquid. Wear 2005;258:1280-7.

[2] Bogdanski S, Lewicki P. 3D model of liquid entrapment mechanism for rolling contact fatigue cracks in rails. Wear 2008;265:1356-62.

[3] Farjoo M, Pal S, Daniel W, Meehan PA. Stress intensity factors around a 3d squat form crack and prediction of crack growth direction considering water entrapment and elastic foundation. Engng Fract Mech 2012;94:37-55.

[4] Fletcher D, Hyde P, Kapoor A. Modelling and full-scale trials to investigate fluid pressurisation of rolling contact fatigue cracks. Wear 2008;265:1317-24.

[5] Fletcher D, Smith L, Kapoor A. Rail rolling contact fatigue dependence on friction, predicted using fracture mechanics with a three-dimensional boundary element model. Engng Fract Mech 2009;76:2612-25.

[6] Farjoo M, Daniel W, Meehan PA. Modelling a squat form crack on a rail laid on an elastic foundation. Engng Fract Mech 2012;85:47-58.

[7] Seo J, Kwon S, Jun H, Lee D. Fatigue crack growth behavior of surface crack in rails. Proc Engng 2010;2:865-72.

[8] Dubourg M-C, Lamacq V. A predictive rolling contact fatigue crack growth model: onset of branching, direction, and growth: role of dry and lubricated conditions on crack patterns. J Tribol 2002;124:680-8.

[9] Brouzoulis J, Ekh M. Crack propagation in rails under rolling contact fatigue loading conditions based on material forces. Int J Fatigue 2012;45:98-105.

[10] Ringsberg JW, Lindåck T. Rolling contact fatigue analysis of rails culminating numerical simulations of the rail manufacturing process and repeated wheel-rail contact loads. Int J Fatigue 2003;25(6):547-58.

[11] Kapoor A, Fletcher DI. Post hatfield rolling contact fatigue: the effect of residual stress on contact stress driven growth in rail. Part 1: The model. Tech. rep., Newrail Report No. WR061106-2, Newcastle University; 2006.

[12] Kapoor A, Fletcher DI. Post hatfield rolling contact fatigue: the effect of residual stress on contact stress driven growth in rail. Part 2: Data. Tech. rep., Newrail Report No. WR061106-3, Newcastle University; 2006.

[13] Kapoor A, Fletcher DI. Post hatfield rolling contact fatigue: the effect of residual stress on contact stress driven growth in rail. Part 3: Further data. Tech. rep., Newrail Report No. WR061106-4, Newcastle University; 2006.

[14] Kapoor A, Fletcher DI, Franklin FJ, Smith L, Hyde P. Post hatfield rolling contact fatigue, the effect of residual stress on contact stress driven growth in rail: comparison of the hatfield and alternative uk rails using models to assess the effect of residual stress on crack growth from rolling contact fatigue. Tech. rep., Newrail Report No. WR061106-6, Newcastle University; 2006.

[15] Seo JW, Goo BC, Choi JB, Kim JK. Effects of metal removal and residual stress on the contact fatigue life of railway wheels. Int J Fatigue 2008;30:2021-9.

[16] Bogdanski S, Brown MW. Modelling the three-dimensional behaviour of shallow rolling contact fatigue cracks in rails. Wear 2002;253:17-25.

[17] Mellings S, Baynham J, Adey R. Automatic crack growth prediction in rails with bem. Engng Fract Mech 2005;72:309-18.

[18] Nguyen-Tajan M, Funfschilling C. A numerical modeling strategy for the rolling contact fatigue analysis of rails. In: WCRR; 2011.

[19] Dolbow J, Moës N, Belytschko T. An extended finite element method for modelling cack growth with frictional contact. Comput Methods Appl Mecha Engng 2001;190:6825-46.

[20] Ribeaucourt R, Baietto-Dubourg M-C, Gravouil A. A new fatigue frictional contact crack propagation model with the coupled x-fem/latin method. Comput Methods Appl Mech Engng 2007;196:3230-47.

[21] Géniaut S, Massin P, Moës N. A stable 3d contact formulation for cracks using x-fem. Eur J Comput Mech 2007;16(2):259-75.

[22] Liu F, Borja RI. Stabilized low-order finite elements for frictional contact with the extended finite element method. Comput Methods Appl Mech Engng 2010;199:2456-71.

[23] Pierrès E, Baietto M-C, Gravouil A, Morales-Espejel G. 3d two scale x-fem crack model with interfacial frictional contact: application to fretting fatigue. Tribol Int 2010;43:1831-41.

[24] Gravouil A, Pierrès E, Baietto M-C. Stabilized global-local x-fem for 3d non-planar frictional crack using relevant meshes. Int J Numer Methods Engng 2011;88:1449-75.

[25] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Engng 1999;46:131-50.

[26] Melenk J, Babuska I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Engng 1996;39:289-314.

[27] Trollé B, Gravouil A, Baietto M-C, Nguyen-Tajan TML. Optimization of a stabilized x-fem formulation for frictional cracks. Finite Elem Anal Des 2012;59:18-27.

[28] Skyttebol A, Josefson B, Ringsberg J. Fatigue crack growth in a welded rail under the influence of residual stresses. Engng Fract Mech 2005;72(2):271-85.

[29] Dang Van K, Maitournam M, Moum Z, Roger F. A comprehensive approach for modeling fatigue and fracture of rails. Engng Fract Mech 2009;76:2626-36.

[30] Nguyen Q. On the elastic plastic initial-value problem and its numerical integration. Int J Numer Methods Engng 1977;11:817-32.

[31] Dang Van K, Maitournam M. Steady-state flow in classical elastoplasticity: application to repeated rolling and sliding contact. J Mech Phys Solids 1993;41(11):1691-710.

[32] Bold PE, Brown M, Allen R. A review of fatigue crack growth in steels under mixed-mode i and ii loading. Fatigue Fract Engng Mater Struct 1992;15:965-77.

[33] Hourlier F, Dhondt H, Truchon M, Pineau A. Etude du branchement des fissures de fatigue sous sollicitations bimodales non proportionnelles: influence de la cinétique de propagation du matériau. Tech. rep., Institut de reccherche de la sidérurgie française IRSID; 1982.

[34] Lamacq V, Dubourg M-C. Modelling of initial fatigue crack growth and crack branching conditions under fretting conditions. Fatigue Fract Engng Mater Struct 1999;22:535-42.

[35] Dubourg M-C, Lamacq V. Stage ii crack propagation direction under fretting fatigue loading: a new approach in accordance with experimental observations. In: Hoeppner DW, Chandrasekaran V, Elliot CB, editors. Fretting fatigue: current technology and practices, ASTM, 1367, American Society for Testing and Materials, Baltimore, MD; 2000. p. 436-50.

[36] Baietto M-C, Pierrès E, Gravouil A, Berthel B, Fouvry S, Trollé B. Fretting fatigue crack growth simulation based on a combined experimental and xfem strategy. Int J Fatigue 2013;47:31-43.

[37] Dubourg M-C, Villechaise B. Analysis of multiple fatigue cracks. Part i: Theory. ASME J Tribol 1992;114:455-61.

[38] Dubourg M-C, Godet M, Villechaise B. Analysis of multiple fatigue cracks. Part ii: Results. ASME J Tribol 1992;114:462-8.

[39] http://www.cast3m.cea.fr/.

[40] Tillberg J, Larsson F, Runesson K. A study of multiple crack interaction at rolling contact fatigue loading of rails. Proc Inst Mech Eng, Part F: J Rail Rapid Transit 2009;22(4):319-30.

[41] Cannon D, Pradier H. Rail rolling contact fatigue research by the european rail research institute. Wear 1996;191:1-13.

[42] Steenbergen M, Dollevoet R. On the mechanism of squat formation on train rails. Part i: Origination. Int J Fatigue 2013;47:361-72.

[43] Grassie SL, Fletcher DI, Hernandez EAG, Summers P. Studs: a squat-type defect in rails. Proc Inst Mech Eng, Part F: J Rail Rapid Transit 2012;243-56.