Peridynamic Analysis of Rail Squats

Featured Application

This study highlights a novel application of peridynamics to rail surface defects by improving the computational method for fracture mechanics. Rail squats is a rail surface defect that can consume over 70% of track maintenance budget. This novel application of peridynamics will enable a better preventative and predictive track maintenance strategy, enhancing public safety while saving hundreds of million euros annually.

Abstract

Rail surface defects are a serious concern for railway infrastructure managers all around the world. They lead to poor ride quality due to excess vibration and noise; in rare cases, they can result in a broken rail and a train derailment. Defects are typically classified as ‘rail studs’ when they initiate from the white etching layer, and ‘rail squats’ when they initiate from rolling contact fatigue. This paper presents a novel investigation into rail squat initiation and growth simulations using peridynamic theory. To the best of the authors’ knowledge, no other comprehensive study of rail squats has been carried out using this approach. Peridynamics are well-suited for fracture problems, because, contrary to continuum mechanics, they do not use partial-differential equations. Instead, peridynamics use integral equations that are defined even when discontinuities (cracks, etc.) are present in the displacement field. In this study, a novel application of peridynamics to rail squats is verified against a finite element solution, and the obtained simulation results are compared with in situ rail squat measurements. Some new insights can be drawn from the results. The outcome exhibits that the simulated cracks initiate and grow unsymmetrically, as expected and reported in the field. Based on this new insight, it is apparent that peridynamic modelling is well-applicable to fatigue crack modeling in rails. Surprisingly, limitations to the peridynamic analysis code have also been discovered. Future work requires finding an adequate solution to the matter-interpenetration problem.

1. Introduction

Rail surface defects are a critical safety concern for railway infrastructure owners and operators all over the world. They undermine the safety and operational reliability of both moderate- and high-speed trains in passenger suburban, metro, urban, mixed-traffic, and freight rail systems. Furthermore, the cost of rail replacements due to such defects has become a significant portion of the whole track maintenance costs, especially in European countries, e.g., Austria, Germany, and France [1].

Traditionally, two different defects are classified: rail studs and rail squats [2]. Rail studs initiate from the white etching layer (WEL) due to wheel slides or excessive traction and grow horizontally 3–6 mm below the rail surface. Rail squats propagate from surface cracks initiated by rolling contact fatigue (RCF), and grow at similar depth of 3–6 mm below the rail surface. Both defects are shown in Figure 1. As a result, the rail surface becomes depressed and passing wheels create excess vibration, noise, and impact loads. This leads to uncomfortable rides for passengers [3], and in cases where impact forces exceed acceptable limits the safety of track components can be compromised [3,4,5,6,7]. Rail squats and studs have been observed in all arrays of track geometries and gradients, in all types of track structures, and in all operational rail traffics. Squats are often found in tangent tracks, in high rails of moderate-radius curves, and in turnouts with vertical, unground rails. Due to the high potential damage caused by rail squats and studs, several research and development projects have been initiated around the world to investigate the causes of, and feasible solutions to, these defects.

Computational rail squat and stud modeling has been the topic of several studies. A finite-element (FE) analysis with a two-dimensional (2D) elastic-plastic model under the assumption of plane strain was used to investigate crack growth from the WEL in [8]. The researchers found that the crack growth direction in the interface between the base material and the WEL is determined by the discontinuity of a material rather than the stress state, and that cracks tend to grow along the interface between the WEL and the rail material, because it is comparatively hard for a crack to propagate into the rail material. Field observations and a numerical analysis in [9] showed that squats initiate as a result of differential wear and differential plastic deformation. Numerical simulations in [10] have also shown that the growth of squats is related to some eigenmodes of the wheel–track interaction system and the high-frequency vibration at wheel–rail contact plays an important role. The probability of rail squat initiation from surface defects based on a transient stress analysis was studied in [11] using an FE model of the vehicle–track interaction. The results showed that when a defect is smaller than 6 mm, its chance to grow into a squat is very small, and when it is larger than 8 mm and in the middle of the running band, the chance is large. RCF occurring on Chinese high-speed rails and wheels was investigated in [12]. Based on field observations and a numerical simulation, it was concluded that indentations seem to be the main cause of RCF. If relatively small but deep indentations exist, then peak von Mises stress can occur both on the surface and at the bottom of the crack, but stress at the bottom is likelier to create RCF cracks [13,14].

The development of rail squats is most commonly studied using the finite-element method, which is based on the classical continuum mechanics theory. It uses spatial derivatives, which do not exist when the displacement field is discontinuous, i.e., when cracks are present. So, as a remedy, techniques of fracture mechanics must be used; however, their major drawback is that the crack path must be known a priori. Due to such limitations, independent crack branching is difficult to implement.

Peridynamic (PD) theory [15,16] was created as an alternative to continuum mechanics for problems with cracks, voids, and other discontinuities. PD uses integral, not partial-differential, equations and deformation instead of strain to compute internal forces. Since integral equations are defined even when the displacement field is discontinuous, this theory is well-suited for fracture studies. Contrary to FE analysis, in PD cracks initiate automatically and grow according to some prescribed damage law. The crack path does not have to be set at the beginning of a simulation, resulting in a more natural crack growth with branching. This theory has a large potential in fracture problems, and has been used to study damage in fiber-reinforced laminated composites [17,18,19], glass [20,21], wood [22], concrete [23,24,25], and steel [26].

In this study, rail squats are simulated using ordinary state-based peridynamic theory (PD). This technique is the recent fundamental development from the original bond-based PD theory. The state-based PD has made a significant advancement in capturing sufficient behaviors of real materials. To the best of the authors’ knowledge, and based on a critical review of the open literature, this paper is the first to present a comprehensive study of rail squat initiation and growth using peridynamic theory. We have presented some initial findings in [13], but this paper describes the full development of the model, the calibration of its parameters, and an application of coordinate-variable loads. Simulation results are compared to field measurements from [14]. The insight from this study is novel and can help further improve the technique for applications of the new theory of peridynamics to real-world problems, and help to enhance better prognostics of rail squats. Overall, the insight enhances an alternative computational method for fracture mechanics.

2. Methods

2.1. State-Based Peridynamic Theory

A brief overview of state-based peridynamic theory is presented in the following paragraphs. An extended overview can be found in [27,28,29]. A peridynamic body consists of some number of nodes each uniquely described by its volume $V_i$, density ${\rho }_i$, and position vector in the reference configuration ${\mathit{x}}_i$. An example of a 2D body is shown in Figure 2. Node ${\mathit{x}}_i$ interacts with other nodes ${\mathit{x}}_j$ through bonds (relative position vectors) ${\mathit{\xi }}_{ij}={\mathit{x}}_j-{\mathit{x}}_i$. These interactions are limited to a range called the horizon $\delta$. Nodes ${\mathit{x}}_j$ that are connected to ${\mathit{x}}_i$ are called the family of ${\mathit{x}}_i$, ${\mathit{H}}_{{\mathit{x}}_i}$. When a body deforms, node ${\mathit{x}}_i$ experiences displacement ${\mathit{u}}_i$ and moves to its deformed position ${\mathit{y}}_i={\mathit{x}}_i+{\mathit{u}}_i$. The bond in the deformed configuration is ${\mathit{y}}_j-{\mathit{y}}_i$. This deformation creates a bond force density vector ${\mathit{t}}_{ij}$ that depends on the collective deformation of all nodes in ${\mathit{H}}_{{\mathit{x}}_i}$ and an opposite bond force density vector ${\mathit{t}}_{ji}$ that depends on the collective deformation of ${\mathit{H}}_{{\mathit{x}}_j}$. Bond forces are force densities (force per volume), not stresses (force per area), because each node describes some volume. The bond deformation vectors are stored in an array called the deformation state

$${\mathit{Y}}_{{\mathit{x}}_i}=\left\{\begin{array}{c}{\mathit{y}}_1-{\mathit{y}}_i\\ ⋮\\ {\mathit{y}}_n-{\mathit{y}}_i\end{array}\right\},$$

(1)

similarly, the force density vectors are stored in an array called the force state

$${\mathit{T}}_{{\mathit{x}}_i}=\left\{\begin{array}{c}{\mathit{t}}_{i1}\\ ⋮\\ {\mathit{t}}_{in}\end{array}\right\}.$$

(2)

The bond force density vectors are computed using bond deformations:

$$\mathit{T}\left({\mathit{x}}_i\right)=\mathit{T}\left(\mathit{Y}\left({\mathit{x}}_i\right)\right),$$

(3)

where the function $\mathit{T}\left({\mathit{x}}_i\right)$ is a material model. It is common to state $\mathit{T}\left({\mathit{x}}_i\right)\langle {\mathit{x}}_j-{\mathit{x}}_i\rangle$ or $\mathit{T}\left({\mathit{x}}_i\right)\langle {\mathit{\xi }}_{ij}\rangle$ and when referring to the force density vector ${\mathit{t}}_{ij}$ in a bond ${\mathit{\xi }}_{ij}={\mathit{x}}_j-{\mathit{x}}_i$, and similarly for deformation state and deformed bond vectors. The peridynamic equation of motion in the integral form is

$$\rho \left({\mathit{x}}_i\right)\ddot{\mathit{u}}\left({\mathit{x}}_i,\mathit{t}\right)=\underset{{\mathit{H}}_{{\mathit{x}}_i}}{{\displaystyle \int }}\left(\mathit{T}\left({\mathit{x}}_i\right)\langle {\mathit{x}}_j-{\mathit{x}}_i\rangle -\mathit{T}\left({\mathit{x}}_j\right)\langle {\mathit{x}}_i-{\mathit{x}}_j\rangle \right)d{V}_{{\mathit{x}}_j}+\mathit{b}\left({\mathit{x}}_i\right)$$

(4)

where $\rho \left({\mathit{x}}_i\right)$ is the density, $\ddot{\mathit{u}}\left({\mathit{x}}_i, t\right)$ is the acceleration, and $\mathit{b}\left({\mathit{x}}_i\right)$ is the external force density.

The contribution of a bond to the force density at a node can be weighed using an influence function $\omega \left({\mathit{x}}_i\right)$. They have been introduced in [16], and their role is further explored in [30]. The value of an influence function can depend on the length, direction, or other bond properties. It can also be used to introduce damage; remove the interaction between two nodes by setting the influence function to 0, i.e., break the bond, when some damage criterion is reached. The simplest damage criterion could be the critical stretch, in which a bond breaks when it is stretched past some critical value $s_c$:

$$\omega \left({\mathit{x}}_i\right)=\{\begin{array}{cc}1,& if{s}_{ij}<s_c\\ 0,& if{s}_{ij}\ge s_c\end{array},s_{ij}=\frac{\left|{\mathit{y}}_j-{\mathit{y}}_i\right|-\left|{\mathit{x}}_j-{\mathit{x}}_i\right|}{\left|{\mathit{x}}_j-{\mathit{x}}_i\right|}=\frac{\left|\mathbf{Y}\langle {\mathit{\xi }}_{ij}\rangle \right|-\left|{\mathit{\xi }}_{ij}\right|}{\left|{\mathit{\xi }}_{ij}\right|},$$

(5)

where $s_{ij}$ is the bond stretch. Then, the damage at a node can be defined as a ratio between the broken and the initial number of bonds [31]:

$$\varphi \left({\mathit{x}}_i\right)=1-\frac{{{\displaystyle \int }}_{H_{{\mathit{x}}_i}}\omega \left({\mathit{x}}_i\right)d{V}_{{\mathit{x}}_j}}{{{\displaystyle \int }}_{H_{{\mathit{x}}_i}}d{V}_{{\mathit{x}}_j}}.$$

(6)

The PD fatigue damage model used in this study was introduced in [32] and used in [33,34,35]. Other researchers have also developed fatigue damage models [36,37]; however, these models use bond-based peridynamic theory and simulate only the crack growth phase. A small overview of the model is given here for completeness; Equations (7) through (11) were first presented in [32].

A body undergoes some cyclic deformation between two extremes + and −, then bond strains at each extreme are $s_{ij}^{+}, s_{ij}^{-}$ and the cyclic bond strain ${\epsilon }_{ij}$ is:

$$s_{ij}^{+}=\frac{\left|{\mathbf{Y}}^{+}\langle {\mathit{\xi }}_{ij}\rangle \right|-\left|{\mathit{\xi }}_{ij}\right|}{\left|{\mathit{\xi }}_{ij}\right|},s_{ij}^{-}=\frac{\left|{\mathbf{Y}}^{-}\langle {\xi }_{ij}\rangle \right|-\left|{\mathit{\xi }}_{ij}\right|}{\left|{\mathit{\xi }}_{ij}\right|},{\mathit{\epsilon }}_{ij}=\left|s_{ij}^{+}-s_{ij}^{-}\right|.$$

(7)

For each bond, a variable called the “remaining life” ${\lambda }_{ij}\left({\mathit{x}}_i, {\mathit{\xi }}_{ij},N\right)$ is defined. It degrades at each loading cycle $N$, and a bond breaks when the remaining life is reduced to zero:

$${\lambda }_{ij}\left(N\right)\le 0.$$

(8)

At the beginning, when $N=0$:

$${\lambda }_{ij}\left(0\right)=1,$$

(9)

at each cycle in the crack nucleation phase (phase I), the change of $\lambda$ is given by

$$\frac{d{\lambda }_{ij}}{dN}\left(N\right)=\{\begin{array}{ll}-A_I{\left({\epsilon }_{ij}-{\epsilon }_{\infty }\right)}^{m_I}\hfill & ,if{\epsilon }_{ij}>{\epsilon }_{\infty }\hfill \\ 0\hfill & ,if{\epsilon }_{ij}\le {\epsilon }_{\infty }\hfill \end{array},$$

(10)

where ${\epsilon }_{\infty }$ is the fatigue limit under which no fatigue damage occurs, and $A_I, m_I$ are parameters for phase I. In phase II, the remaining life changes according to:

$$\frac{d{\lambda }_{ij}}{dN}\left(N\right)=-A_{II}{\epsilon }_{ij}^{m_{II}},$$

(11)

where $A_{II}, m_{II}$ are parameters for phase II.

The transition from phase I to phase II is handled by applying the phase I model to bonds connected to ${\mathit{x}}_i$ until there is a node ${\mathit{x}}_j$ in $H_{{\mathit{x}}_i}$ with damage

$$\varphi \left({\mathit{x}}_j\right)\ge {\varphi }_c,$$

(12)

where ${\varphi }_c$ is the damage at which phase II begins. Then, reset the remaining life of bonds connected to ${\mathit{x}}_i$ to 1 and switch to the phase II model.

2.2. Computational Model

The fatigue damage model was implemented in the open-source PD program Peridigm [38,39]. If the quasi-static analysis acceleration term in (4) is zero, then the peridynamic equation of motion in the discreet form is approximated as:

$${\displaystyle \sum }_{H_{{\mathit{x}}_i}}\left(\mathit{T}\left[{\mathit{x}}_i,t\right]\langle {\mathit{x}}_j-{\mathit{x}}_i\rangle -\mathit{T}\left[{\mathit{x}}_j,t\right]\langle {\mathit{x}}_i-{\mathit{x}}_j\rangle \right)\mathsf{\Delta }{V}_{{\mathit{x}}_j}+\mathit{b}\left({\mathit{x}}_i,t\right)=0.$$

(13)

Two techniques introduced in [32], including implicit strain simulation and time mapping, have been used to speed up the simulations. They are illustrated in Equations (14) through (17). In the case of high-cycle fatigue, the bond strains are below the elastic limit, so an elastic material model can be used. In such cases, the strain in a bond would change linearly between + and − loading conditions, so it is possible to simulate only the + loading condition and compute

$$s^{-}=R{s}^{+},R=\frac{P^{-}}{P^{+}},$$

(14)

where $R$ is the loading ratio, and $P$ is the applied load at each extreme. Then, the cyclic strain is given by:

$$\epsilon =\left|s^{+}-s^{-}\right|=\left|\left(1-R\right)s^{+}\right|.$$

(15)

The loading ratio $R=0$ was used for all simulations. Using linear time mapping ([32] also introduces exponential time mapping, but it was not used in this study), the simulation time $t$ relates to the current cycle through

$$N=\frac{t}{\tau },$$

(16)

where $\tau$ is a constant. Then, the remaining life at step $n$ in a bond ${\mathit{\xi }}_{ij}$ is given by

$$\begin{array}{l}\frac{{\lambda }^n-{\lambda }^{n-1}}{\mathsf{\Delta }t}=-A{\left({\epsilon }_{ij}^n\right)}^m\hfill \\ \frac{\mathsf{\Delta }\lambda }{\mathsf{\Delta }t}=\frac{\mathsf{\Delta }\lambda }{\mathsf{\Delta }N}\frac{\mathsf{\Delta }\mathrm{N}}{\mathsf{\Delta }t}\hfill \end{array}\to {\lambda }^n={\lambda }^{n-1}-\frac{t^n-t^{n-1}}{\tau }A{\left({\epsilon }_{ij}^n\right)}^m.$$

(17)

In PD, unlike in a fracture mechanics model with a pre-crack, it is almost impossible to develop a sharp crack surface; instead, a damaged zone is developed. For example, if damage at a node is 0.5, it could mean that 25% of the bonds on the opposite sides of a node are broken, and it could also mean that 15% are broken on one side and 35% on the other.

Crack growth in phase II is faster than that in phase I, so switching to phase II sooner would lead to faster crack growth and more conservative results. In this study, the emphasis is placed on the consideration of the most conservative case, i.e., switching to phase II at the lowest damage when it can be thought that a crack has appeared. Such a situation would happen when all of the bonds on one side of a node are broken, but other bonds remain intact. The 2D case, if the horizon is 3 times the distance between the nodes, is shown in Figure 2. An equivalent case in three dimensions (3D) would have 47 broken bonds and 75 unbroken bonds, i.e., damage of 0.385. Therefore, the fatigue model transitions from phase I to phase II when the damage at a node reaches ${\varphi }_c=0.385$.

2.3. Model of a Rail

Initially, a model of a whole UIC60 rail head had been developed. However, due to the required fine discretization, computational resources, and time constraints, only a part of the rail head has been modeled. Mesh was firstly created in the Ansys FE program using 3D eight-node solid elements. Afterwards, element centroid coordinates and their volumes were exported and converted to Peridigm’s mesh file, and both models are shown in Figure 3. The dimensions of the model were 0.03 × 0.024 × 0.03 m with a node size of $h=0.0005$ m and the horizon of $\delta =0.0015$ m. The load area, see Section 2.5. Boundary conditions, was about 0.013 × 0.013 m, which means that the distance between the edge of the model and the load area was about $8\delta$. A cartesian left-handed coordinate system was used, and the center of the model’s bottom face was located at the origin. The top face was made of $R=300$ mm and $R=80$ mm arcs, as in the specifications of the UIC60 rail.

Since the top surface was curved and mapped meshing was used, the nodes were not perfectly cubic. However, the difference between the average node volume and the volume of a cubic node is only 1.17% (see

Table 1. The cubic versus the smallest, largest, and average node volume.

Parameter	Volume, m3	% Difference
Cubic	1.25000 × 10−10	0.00%
Min	1.18960 × 10−10	−4.83%
Max	1.29750 × 10−10	3.80%
Average	1.26464 × 10−10	1.17%

). This is relevant when converting the applied loads from stress to force density, but since the difference is small, the impact is negligible.

A Linear Peridynamic Solid (LPS) [16] material model was used. The material properties were: density 7850 kg/m³, Poisson’s ratio 0.3, and Young’s modulus 189.9 GPa, obtained from [40]. The LPS model is the peridynamic equivalent to the elastic material model in continuum mechanics. It has been selected because the applied loads do not cause the material to exceed its yield strength.

To verify the peridynamic model, the displacement in the X and Y directions in an undamaged state has been compared against an FE model. The same model was used for verification of the Peridigm’s mesh creation. Movements are restricted in all directions for all nodes within $1\delta$ from the bottom. Two loadings—vertical pressure and surface shear traction—are applied to the load area at the top of the rail head. Both loadings are applied as functions in terms of node x and z coordinates; for the exact functions, please see Section 2.5 Boundary conditions. In the FE model, the pressure has been applied as distributed pressure on the top face of solid elements in the load area, and traction has been applied as horizontal forces acting on the top three layers of nodes within the load area.

Figure 4 shows the X and Y displacement in the cross-section along the centerline of the rail, and

Table 2. Maximum and minimum displacement values in the X and Y directions from the finite-element (FE) and Peridynamic (PD) simulations.

Value	X			Y
	FE, m	PD, m	Difference	FE, m	PD, m	Difference
Max	2.03 × 10−5	1.90 × 10−5	−6.95%	2.35 × 10−6	2.55 × 10−6	7.92%
Min	−6.59 × 10−7	−8.30 × 10−7	20.61%	−4.69 × 10−5	−5.07 × 10−5	7.48%

presents the maximum and minimum displacement values. It is clear that a very good agreement between models can be found. In fact, if the difference in extreme displacement values would be within ±10%, the displacements would be similar in the cross-sections. Figure 4 definitely shows a similar displacement distribution between the FE and PD models. The maximum displacement values between the FE and PD models are −6.95% and 7.92% for the X and Y directions, respectively. The difference in the minimum Y displacement is 7.48%. In the X displacement, however, it is 20.61%, which is more than what should be considered a good agreement. Though the relative difference in the X displacement is large, it should not negatively influence the simulation results, because the region with low X displacement is far from the load area, and, therefore, far from the area of interest with the growing cracks. Additionally, the absolute difference is very small, i.e., only 1.71 × 10⁻⁷ m. It should have no effect on the PD model’s accuracy.

2.4. Fatigue Damage Model Parameters

This study adopts the rail steel data from Figure 4a and Figure 5 in [40], and follows the procedure to obtain damage model parameters in [32]. Although [40] presents quite old data, it contains ε-N (strain-life), K-da/dN (Paris law), and material properties data. This is beneficial, because it assures that the data are for the same material. Other fatigue data sources have been considered [41,42,43,44,45,46], but either have only the S-N curves available, contain less data points, or do not have both ε-N and K-da/dN plots. The fatigue damage model parameters are shown in

Table 3. Fatigue damage model parameters.

	Phase I	Phase II
A	426.00	25,237.48
m	2.77	4.00
ε∞	0.00186	--

.

The parameters for phase I ($A_I, m_I,{\mathit{\epsilon }}_{\mathit{\infty }}$) are found from functions fitted to the ε-N plot (see Figure 5). The fitted power law function takes the same form, $y=a{\mathit{x}}^b$, as the phase I damage model in (10). So, parameter $m_I$ is the inverse of slope $b$:

$$b=-\frac{1}{m_I}$$

(18)

and parameter $A_I$ is calculated from the value of intercept:

$$a=\frac{-log{A}_I}{m_I}\Rightarrow A_I=\frac{1}{a^{m_I}}.$$

(19)

The fatigue limit of rail steel was determined from the function:

$${\left(\frac{\Delta \epsilon }{2}-{\epsilon }_{\infty }\right)}^{\xi }N=C,$$

(20)

where $\frac{\mathsf{\Delta }\epsilon }{2}$ is the strain amplitude, ${\epsilon }_{\infty }$ is the fatigue limit, $N$ is the number of cycles, and $\xi , C$ are constants.

A Paris law plot is required to find the parameters for phase II. In this study, R = 0.05, and moist air data from Figure 5 in [40] were used. The plot is replicated in Figure 6. The fatigue damage model in (11) has the same form as the Paris law for crack growth:

$$\begin{array}{c}\frac{da}{dN}=c\Delta K^M,\end{array}$$

(21)

where $\frac{da}{dN}$ is the crack growth speed, $c, M$ are constants, and $\mathsf{\Delta }K$ is the cyclic stress intensity factor. $\Delta K$ is proportional to the bond strain at the crack tip (in [32] called ${\epsilon }_{core}$); therefore, $m_{II}=M$, so this parameter can be obtained directly from a Paris law plot. The remaining parameter $A_{II}$, however, cannot. Instead, a simulation with some trial value $A_{II}^{\prime }$ has to be run to obtain the trial crack growth speed ${\left(\frac{da}{dN}\right)}^{\prime }$. Then, the real $A_{II}$ value can be found from [32]:

$$\begin{array}{c}{A}_{II}=A_{II}^{\prime }\frac{\frac{da}{dN}}{{\left(\frac{da}{dN}\right)}^{\prime }}\end{array}.$$

(22)

To find $A_{II}$, a single edge notch (SEN) specimen with a pre-crack in uniaxial tension is simulated. The stress intensity at a crack tip is given by:

$$K_I=\sigma \sqrt{\pi a}F\left(\frac{a}{b}\right),$$

(23)

$$F\left(\frac{a}{b}\right)=1.122-0.231\left(\frac{a}{b}\right)+10.550{\left(\frac{a}{b}\right)}^2-21.710{\left(\frac{a}{b}\right)}^3+30.382{\left(\frac{a}{b}\right)}^4,$$

(24)

where $\sigma$ is the applied stress, $a$ is the crack length, and $b$ is the specimen width. The crack tip’s location was defined as the maximum x coordinate at which all nodes through the depth of the model have damage of at least 0.385.

The model’s size is 0.05 × 0.008 × 0.003 m. It has been discretized with 150,000 nodes with a spacing of 0.0002 m using the mesh-free method described in [31]. The horizon is set to a little over 3 times the nodal spacing: 0.0006001. The model has a 0.005-m-long pre-crack on the left side to ensure that Equation (23) is applicable. A force density of 6.25 × 10¹⁰ N/m³ (equivalent to 50 MPa) has been applied to nodes within one $\delta$ of both the top and bottom, and damage is disabled for nodes within $3\delta$ from the top and bottom, to avoid unphysical behavior near the boundary conditions. Crack growth speed data only from phase II are required, so switching to phase II at low damage reduces the simulation time. The damage required for transition from phase I to phase II has been, therefore, set to 0.017. For the trial simulation, $A_{II}^{\prime }=1e6$ and $m_{II}=\mathrm{4.00.}$ An LPS material model with the same parameters as for the rail head simulation is used. The first simulation (with $A_{II}^{\prime }$) ran for 163,100 cycles, after which the crack turned upward, so Equation (23) is no longer accurate; the second simulation runs for 13,275,999 cycles until the crack splits in two. Figure 7 shows the simulation with $A_{II}$ at cycle 309,999 (top) and step 13,275,999 (bottom). The number of cycles is large because a low applied stress causes fatigue damage to increase slowly.

Since a crack grows in discrete jumps between nodes, the crack growth speed between two cycles $m$ and $n$ with such jumps has been calculated as the difference in crack length divided by the difference between the current cycle and the cycle at which the previous jump occurred:

$${\left(\frac{da}{dN}\right)}^{\prime }=\frac{a_n-a_m}{N_n-N_m},$$

(25)

where $a$ is the crack length, and $N$ is the number of cycles. Then, the ${\left(\frac{da}{dN}\right)}^{\prime }$ values are interpolated to match the $\mathsf{\Delta }K$ values from the experimental data and $A_{II}$ is calculated using Equation (22). In total, 22 $A_{II}$ values have been calculated. These values vary greatly, and the coefficient of variation is 0.7134; therefore, an average value has been used. A repeated simulation with $A_{II}$ and not $A_{II}^{\prime }$ (see Figure 6) exhibits a very good agreement with the first part of the experimental data. It was not possible to determine agreement with the latter part of the data, because the simulated crack splits into two and Equation (18) could no longer be used. A better approach (with less variance between the calculated $A_{II}$ values) might be to use the real crack growth rate $\frac{da}{dN}$ not from experimental data, but from the fitted Paris Law function. This approach will be explored in future research.

The model of a rail head uses coarser discretization than the model of an SEN specimen. Since the horizon has been kept at three node spacings for both, the actual value of the horizon is different in both simulations: 0.0015001 and 0.0006001, respectively. A change in horizon does not change the $A_I,m_I,m_{II}$ parameters (see chapter 4.3 in [32] for details), but $A_{II}$ has to be scaled with the horizon. Equation (29) in [32] provides the means to do that:

$$A_{II}\left(\delta \right)=\widehat{A_{II}}{\delta }^{\frac{m_{II}-2}{2}}\to A_{II}\left(\delta \right)=16,823,863\times {0.0015001}^{\frac{4-2}{2}}=25,237.48,$$

(26)

where $\widehat{A_{II}}$ is independent of $\delta$.

2.5. Boundary Conditions

Since nodes describe some volume, boundary conditions (BCs) must also be applied to some volumes. A BC layer thickness equal to the horizon was recommended in [47]. In [13], the researchers similarly applied loads only to a single layer of nodes on top of the rail head, and such an approach lead to poor results. BC nodes separated from the rest just after 26 thousand cycles, due to the low number of bonds over which the applied loads were distributed.

Displacement in all directions was fixed for nodes within one $\delta$ from the bottom. Additionally, damage was disabled for nodes within $3\delta$ from the bottom to avoid the concentration of unphysical damage near the BC layer.

Train wheel load data from [48] have been used in this study. The wheel–rail contact area (Figure 4f in [48]) is centered at the coordinate origin, see Figure 3b, and approximated with an ellipse with a half-axis $a=0.0066$ m, $c=0.006386$ m. Two different train wheel loadings are used: vertical pressure and surface shear traction. They are applied to a $1\delta$ thick layer on top of the rail head.

The vertical force density, in N/m³, from the elastic pressure (data from Figure 5f in [48]) can be computed from a modified ellipsoid’s formula:

$$p=\frac{1.116\times {10}^9}{h}\sqrt{1-\frac{x^2}{{0.0066}^2}-\frac{z^2}{{0.006386}^2}},$$

(27)

where $p$ is the force density (N/m³), $x, z$ are node coordinates (m), and $h$ is the node size (m). Since loads are applied to a $1\delta$ (three node spacings) thick layer, the computed value at a position $\left(x,z\right)$ has been divided by 3 and applied to each of three nodes under this position.

Shear traction forces are taken from Figure 5f and Figure 6f in [48], where they are given as a stress distribution over an area. Mesh-free discretization requires that loads are applied to discrete nodes, so the shear traction data over the whole load area had to be described by some function from which an exact value at a node could be calculated. Only half of the load area is considered, because the traction data were symmetric. Half of the load area is divided into four parts along the z axis, see Figure 8, and the shear traction values in each part are described by a tri-linear function, see Figure 8 and Figure 9. Stress values from [48] can be plotted with symbols in Figure 8 and fitted with tri-linear functions from which the exact shear traction force value at each node could be calculated. Since loads are applied to a three-node-thick layer, the calculated stress values must be converted to force density and divided by 3 before being applied to nodes.

3. Results

This problem has been simulated on a computing cluster at Riga Technical University using 4 × 36 cores. Each simulation was run for 42,884 cycles, after which the solver failed to converge. The results are shown in Figure 10 (cross-section in the longitudinal direction) and Figure 11 (cross-section in the transversal direction).

The simulation results are compared with the rail squat field measurements in [14]. Cracks have been measured at specified grid points on a tangent rail over a span of 3 years using a handheld ultrasonic testing device with the accuracy range of ±0.1 mm. The field measurement results can be seen in Figure 12 and Figure 13.

Two loadings—pressure and shear traction—have been applied to the model. The magnitude of the shear traction is not symmetric around the coordinate origin, even although the load area is. Traction is then applied in the positive x direction, and the values change as shown in Figure 8 and Figure 9. This can cause damage to develop slightly asymmetrically against the y-z plane, which is best seen in Figure 10a. Damage develops faster on the positive side of the x axis (the right side in Figure 10). Against the x-z plane, in Figure 11, damage developed symmetrically, because both the pressure and the shear traction are symmetric. The same asymmetric crack growth has been observed in the field measurements (see Figure 13). In reality, such asymmetry happens because the shear traction from a wheel rolling forward is applied in the rolling direction.

Figure 11a clearly shows that damage first develops close to the location of maximum pressure (the middle of the rail in transversal direction). In addition, the maximum damage remains under the same area (see Figure 11b–d). This is consistent with the field measurements shown in Figure 12. Cracks are deeper closer to the center of the rail head and shallower closer to the sides. This shows that they first initiated and have been growing for longer under the rolling surface.

The simulation ended unexpectedly quickly, because the fatigue resistance of a rail without any defects should definitely be above 42,884 cycles. Loads have thus been applied to a three node (one horizon) thick layer on top of the rail head. As bonds extending to nodes below this layer are broken, the applied loads are no longer transferred downwards and the loaded nodes simply moved through the layers below them; this can be seen in the c and d parts of Figure 10 and Figure 11. This is the so-called “matter-interpenetration” problem, and usually it is solved through different contact models. The simplest model—short-range force—was introduced in [31], and a better description is available in [39,49]. Other contact models and properties are presented in [50,51,52,53].

While a shot-range force contact mode is available in Peridigm, it has been implemented only for explicit and not quasi-static simulations and it does not consider contact between nodes that are bonded initially. As damage develops, it is important that the contact model reconsiders the contact between two nodes that were bonded initially but are not anymore. This limitation has been explored, and it is possible to resolve it. Future work will concentrate on how to efficiently pass data between Peridigm’s damage and contact models.

4. Discussion and Concluding Remarks

This study used a new approach to rail squat simulation: the peridynamic theory. It describes the derivation of model’s parameters, and illustrates how to apply a variable loading that is dependent on a node’s location. The simulation successfully captures the initiation of, and initial, rail squat growth. Due to limitations of the simulation, a larger crack at this stage could not be simulated. However, the simulation results are in excellent agreement with field measurements for the crack initiation phase.

Damage initiates and grows faster close to the location of maximum pressure; similar crack growth has been measured in the field. Additionally, the computational model reveals that the squat damage first grows in the direction of the applied shear traction, and the same has been shown in field measurements.

The computational model experiences a matter-interpenetration problem, where damaged nodes, no longer connected with bonds, move freely through each other, without considering possible contacts. This problem can be solved by applying a contact model; however, contact models in Peridigm do not consider the contact between nodes that were bonded initially. To solve this problem, bond damage data needs to be passed between the damage model and the contact model. Future work will focus on re-developing parts of Peridigm’s code, so that data can be passed between its damage and contact models.