## 2 Methods
### 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.  
Figure 7: A single edge notch (SEN) specimen at: (**a**) 3,099,999 cycles; (**b**) 13,275,999 cycles. Displacements are increased 10 times.  
The vertical force density, in N/m\({}^{3}\), 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.006 386^{2}}}, \tag{27}\]  
where \(p\) is the force density (N/m\({}^{3}\)), \(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 \((x,z)\) has been divided by 3 and applied to each of three nodes under this position.  
Shear traction forces are taken from Figures 5 and 6 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 Figures 8 and 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.  
Figure 8: Half of the load area divided into four parts with a tri-linear function for each part. Functions describe the shear traction stress values in the load area. The other half of the load area is a mirror image. The axis directions and node size are the same as in the rail head’s model.  
Figure 9: Surface shear traction data from [48] (shown with symbols) and the tri-linear functions used to describe the shear traction values in the load area.