## 2 Methods
### 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 \(\times\) 0.024 \(\times\) 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 \(\times\) 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). 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\({}^{3}\), [PERSON]'s ratio 0.3, and [PERSON]'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 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 \(\pm 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\times 10^{-7}\) m. It should have no effect on the PD model's accuracy.