*2.2. Methods*

In the course of analysis of the dynamic impact of a rail vehicle on the railway surface, an original computational model was created, allowing determination of the impact of differentiated rail support on the dynamic response of the entire structure. The results of numerical calculations were verified in field tests carried out on the actual railway line in places where there was a change in the stiffness of the elements of the railway line, loaded with vehicles of different weights, moving at different speeds. The adopted research method is a theoretical model and uses experimental verification.

#### 2.2.1. Theoretical Model

As a starting point for development of a model of a dynamically loaded railway surface, the Bernoulli-Euler beam, located on the elastic Winkler substrate, was adopted. Since the railway track has a symmetrical structure, half of its construction was considered, i.e., one rail track and the load falling on it, as well as the load-bearing system. All the supporting layers of the structure below the rail were parameterized by substitute coefficients of elasticity and damping. The model takes into account the vertical displacement of the rail caused by a multiaxis rail vehicle with a constant axle load moving uniformly. A simplification of the mapping of a rail vehicle to a moving concentrated force with constant pressure was adopted. The validity of this approach has been proven, inter alia, in [24] in which a number of variants of rail vehicle mapping with varying levels of complexity were analysed and it was found that the method of load mapping does not have a significant impact on the obtained displacement values of the railway engineering object. The same representations of a moving rail vehicle were also used in [25,26], each time giving correct, real results.

The model in this paper also takes into account the fact that in individual cross-sections, the ground layer and the railway surface can be characterized by variable values, resulting from the use of different materials and different design solutions [10,27]. The model under consideration is presented in Figure 2.

After carrying out the appropriate analysis and after taking into account the data depending on position and time, the differential equation of the deflection line of the rail course (3) was obtained. This is a differential equation of the fourth order for the spatial variable and the second order for time.

$$\frac{\partial^2 \mathbf{w}(\mathbf{x}, \mathbf{t})}{\partial \mathbf{t}^2} \cdot \mathbf{m}\_\mathbf{s} + \frac{\partial \mathbf{w}(\mathbf{x}, \mathbf{t})}{\partial \mathbf{t}} \cdot \mathbf{c}(\mathbf{x}) = -\frac{\partial^4 \mathbf{w}(\mathbf{x}, \mathbf{t})}{\partial \mathbf{x}^4} \cdot (\mathbf{E} \cdot \mathbf{J})(\mathbf{x}) - \mathbf{w}(\mathbf{x}, \mathbf{t}) \cdot \mathbf{k}(\mathbf{x}) + \frac{\mathbf{P}(\mathbf{x}, \mathbf{t})}{\mathbf{d} \mathbf{x}}.\tag{3}$$

#### **Figure 2.** Rail line model.

The vertical displacement of the rail will depend on the place where the load is applied and on the time that has elapsed since the load was applied at a given point. The coefficients of the Equation (3) are described earlier.

The relation presented above was the basis for further numerical analyses. The finite differences method was used to solve it. It makes it possible to solve differential equations by replacing differential operators with difference operators. They are defined on a set of points, called a grid. Grid elements are nodes. The considered rail element was divided into sections of a length " Δx" and then, the values of individual quantities in the grid nodes were calculated. With a sufficiently high density of such a division, it is possible to achieve very accurate results, converging to those provided by the analytical solution. Thanks to this assumption, the initial extended differential equation is replaced by a system of algebraic equations [28]. To derive formulas for difference operators, the assumption was used that the desired values of the function between individual nodes are connected by fragments of the parabola. Difference schemes occurring in the finite differences method refer to the time–space grid with the parameters " Δt" and " Δx". The solutions in nodes, defined by coordinates (x, t), are calculated.

In order to solve the differential Equation (3) by the finite differences method, taking into account dynamic loads, it was necessary to determine the boundary conditions and the initial conditions of the issue under consideration. The following assumptions have been introduced:


In addition, to start the numerical calculation process, one artificial time layer was introduced, and for each subsequent time layer of the calculation process, two artificial nodes for the right and left boundary conditions were introduced. The correctness of such an assumption was confirmed in the calculation process. The time-space grid with computational nodes along with the initial and boundary conditions is shown in Figure 3.

The following boundary conditions were assumed: (1) the vertical displacement in the first "fictitious" node marked by the number "1" is equal to the displacement in the first node on the rail course marked by the number "3" and taken with a minus sign, (2) the vertical displacement in the first support node marked by the number "2" is equal to 0, 3), the vertical displacement in the last support node marked by the number "n + 2" is equal to 0, 4), and the vertical displacement in the last "fictitious" node marked by the number "n+3" is equal to the displacement in the last node on the rail course marked by the number "n + 1" and taken with a minus sign.

To determine the initial conditions, the following assumptions were made: (1) the value of the vertical displacement of a given node at the initial moment is equal to 0, (2) at the initial moment, the first node on the rail course marked by the number "3" begins to vibrate at a certain low speed and low acceleration, the values of which depend on the value of the load.


**Figure 3.** Differential time–space grid—authors.

Thus, in order to start the numerical calculation procedure, it was assumed that the initial value of the vertical displacement of each of the nodes, resulting from the assumed initial conditions, is equal to the initial value at the time of time two steps earlier than the moment when the load is applied. This value depends on: (1) load value, (2) geometric and material characteristics of the rail, (3) surface damping value, and (4) adopted values of the time–space grid of the finite differences method, i.e., (a) spatial step and (b) temporal step.

It is worth emphasizing here that the assumed numerical value of the initial condition is used only to start the calculation process and has a negligible impact on further results and vertical displacements of the rail course calculated only using this value without rolling stock load would be negligibly small, which was confirmed in the calculations carried out.

In the calculations, there are a number of coefficients that ensure the compliance of the material parameters of individual elements of the rail surface and the position of the loading forces with the current position of the loading vehicle on the analyzed rail section. The introduced coefficients are:


Such a set of coefficients made it possible to describe each cross-section at any time with a precisely defined set of data. These coefficients and the others that occur in the calculation process are shown in Table 2.

Particular attention should be paid to the coefficient of effective stiffening of the rail, which reflects the method of fastening the rail in various structural solutions of the railway surface. It was introduced after defining the phenomenon that the effective bending stiffness of the rail is different for ballast and ballastless surfaces [9]. In places where there is a ballast surface, its value is assumed as 1.00; in places where there is a ballastless surface, the value of 1.30 and for sheath rail systems (for example, the ERS system), where the rail is immersed in an elastic mass that limits its freedom of bending, the value of this coefficient is equal to 1.50. It should be noted that these are proprietary values, verified during the in situ research, but may require clarification and detailing with more measurements and analyses.

An analysis of the convergence of spatial and temporal steps was also performed. On its basis, the value of the spatial step Δx at the level of 0.05 m and the time step at a level lower than the critical time, which was defined as:

$$
\Delta \mathbf{t}\_{\mathbf{kr}} = \frac{2\Delta \mathbf{x}^2}{\pi} \cdot \sqrt{\frac{\mathbf{m}\_s}{1.50 \cdot \mathbf{E} \cdot \mathbf{J}}} \,. \tag{4}
$$


**Table 2.** A list of parameters that must be defined before starting the calculation.

As a result of the performed analyses and appropriate transformations, a relation was obtained, which determines the value of the vertical displacement of a given node "i" in the next considered time moment "j + 1":

$$\begin{split} \mathbf{w}\_{\mathrm{i}}\mathbf{w}\_{\mathrm{i}}^{\mathrm{j}+1} &= \frac{1}{\mathbf{S}\_{\mathrm{i}}} \cdot \left\{ \boldsymbol{\Delta}\mathbf{t}^{2} \cdot (\mathbf{S}\_{\mathrm{4}} \cdot \mathbf{w}\_{\mathrm{i}-2}\mathbf{}^{\mathrm{j}} + (-\mathbf{4} \cdot \mathbf{S}\_{\mathrm{4}}) \cdot \mathbf{w}\_{\mathrm{i}-1}\mathbf{i}^{\mathrm{j}} + (\mathbf{6} \cdot \mathbf{S}\_{\mathrm{4}} - \mathbf{z}\_{\mathrm{i}} \cdot \mathbf{k}\_{\mathrm{i}}) \cdot \mathbf{w}\_{\mathrm{i}}^{\mathrm{j}} + (-\mathbf{4} \cdot \mathbf{S}\_{\mathrm{4}}) \cdot \mathbf{w}\_{\mathrm{i}+1}\mathbf{i}^{\mathrm{j}} \\ &+ \mathbf{S}\_{\mathrm{4}} \cdot \mathbf{w}\_{\mathrm{i}+2}\mathbf{i}^{\mathrm{j}} + \mathbf{H}\_{\mathrm{i}}^{\mathrm{j}} \cdot \mathbf{P}\_{\mathrm{i}}) - \mathbf{S}\_{\mathrm{2}} \cdot \mathbf{w}\_{\mathrm{i}}^{\mathrm{j}-1} - \mathbf{S}\_{\mathrm{3}} \cdot \mathbf{w}\_{\mathrm{i}}^{\mathrm{j}} \right\} \end{split} (5)$$

$$\mathbf{S}\_1 = \mathbf{m}\_{\text{e}} + \frac{3}{2} \cdot \mathbf{n}\_2 \cdot \mathbb{H}\_{\text{i}}^{\dagger} \cdot \mathbf{c}\_{\text{i}} \cdot \Delta \mathbf{t} \tag{6}$$

$$\mathbf{S}\_2 = \mathbf{m}\_8 + \frac{1}{2} \cdot \mathbf{n}\_2 \cdot \mathbf{H}\_i^\dagger \cdot \mathbf{c}\_i \cdot \Delta t \tag{7}$$

$$\mathbf{S}\_3 = -2 \cdot (\mathbf{m}\_\\$ + \mathbf{n}\_2 \cdot \mathbf{H}\_i^{\dot{j}} \cdot \mathbf{c}\_i \cdot \Delta t) \tag{8}$$

$$\mathbf{S}\_4 = -\frac{\mathbf{n}\_1 \cdot \mathbf{E} \cdot \mathbf{J}}{\Delta \mathbf{x}^4} \tag{9}$$

$$\mathbf{P}\_1 = \frac{\mathbf{P}}{\Delta \mathbf{x}}\tag{10}$$

The scheme of operation of the created algorithm in relation to individual nodes "i" and time moments "j" is presented in the differential time—space grid shown in Figure 3. An explicit method of solving differential equations was used [29], in which it is important to properly give the initial conditions of the calculation process.

An explicit method of searching for solutions on the time–space grid of the finite differences method was used [29]. The new value wij+<sup>1</sup> (orange in Figure 3) is based on the values already known, calculated in the previous steps (green). The first unknown to be calculated is w31. This is the vertical displacement of the first node located on the rail course at the first considered time moment (orange color). For the first calculation, therefore, five values from the time moment (j = 0, green color) and one fictitious value from the time moment (j = −1, green color) will be used. This shows how important it is to

correctly give the initial condition, because otherwise any data taken into account for the first calculation would be equal to 0.

The horizontal dimension of the grid is equal to "n + 3". The vertical dimension of the grid is marked as the value "m" and corresponds to the number of time steps by which the analysis of the vertical displacement of the rail line in individual nodes will be carried out. It is important that the numerical duration of the analysis allows the value of the maximum displacement of each point to which the dynamic load is applied to be determined. It is also important to determine the number of time steps " Δt" of the whole calculation process. This quotient should be rounded up to the whole to obtain an integer and also increased by 2 to take into account the time moments marked in Figure 3 as " −1" and "0".

Figure 4 shows a flowchart of the described solution algorithm. In an orderly way, it illustrates the subsequent steps necessary to obtain the final result, which are the values of vertical displacement of the rail in individual nodes "i" of the model and in individual time moments "j", caused by the dynamic impact of the rail vehicle on the railway surface.

**Figure 4.** Block diagram of the algorithm.

### 2.2.2. Experimental Verification

The theoretical model should reliably describe reality so that in similar cases it is not necessary to perform field research, but only to use the theoretical modeling. In order to verify and validate the created algorithm, in situ studies of vertical displacements of dynamically loaded railway rail were carried out. The research was carried out using the laser scanning technology. The advantages of measurements performed in this technology include: high accuracy and automation, as well as the speed of measurements and the lack of the need to destroy the tested object or exclude it from operation. Among the disadvantages are: large volume of data and their long processing time, and the inability to conduct measurements during bad weather conditions, as well as the high price of the scanner and its sensitivity to mechanical damage [30]. The choice of research technology was preceded by an analysis of the literature, which gave grounds for obtaining satisfactory accuracy of the results.

Studies of vertical rail displacement caused by a passing train using the laser scanning technology were carried out in cooperation with PKP (PKP Polskie Linie Kolejowe S.A. Railway Lines Plant in Ostrów Wielkopolski, Poland) at the Railway Lines Plant in Ostrów Wielkopolski. The measurements were carried out in tracks located within two engineering structures and their access zones. In order to analyze the impact of the threshold effect on the dynamic impact of the rail vehicle on the railway surface, the following were selected: (1) place with a uniform type of surface on the object and on the trail and (2) place where there is a jumping change in the type of surface in front of and behind the object.

The research was possible thanks to cooperation with P.P.H.WObit E.K.J.Ober s.c. (https://wobit.com.pl/en) (accessed on 11 March 2022). This company made the laser scanner of the scanCONTROL LLT2610-50 series available for testing free of charge under a loan agreemen<sup>t</sup> concluded with the Faculty of Civil Engineering and Geodesy of the Military University of Technology in Warsaw. The measuring set included the following elements: (1) laser scanner of the scanCONTROL LLT2610-50 series, (2) a tripod with mounting screws enabling permanent fixation of the measuring device, (3) a power set, and (4) a cable enabling the scanner to be connected to a computer. The scanner allows to measure up to 2,560,000 points per second at a frequency of up to 4000 Hz. Measurement accuracy is 2 μm. The scanner is enclosed in a metal casing, which ensures an IP65 degree of protection. At the bottom of the scanner there is a light source and a receiver. There are two connectors on the housing: (1) Ethernet (used to communicate with a PC or Plc controller) and (2) multifunctional (used, among others, for power supply). To communicate the device with a computer, a cable with an RJ45 connector on one side and a connector dedicated to the scanner on the other [31] is used. For the proper operation of the entire measuring system, one needs a computer with the appropriate software installed and a power generator providing a constant power source. Before starting measurements, it is necessary to check the correctness of the scanner's connection with the computer, select the appropriate software configuration and assign it to the device, as well as define the scanner's work settings by setting parameters such as exposure time and the number of measured profiles per second. In selection of the correct parameters, a preview of the measurement, which is displayed in real time, in one of the software modules, is helpful. In addition, one must specify a name and location to save the result file. The recording is made automatically after the measurements are completed. The scanner measuring set is shown in Figure 5.

The scanner was placed transversely to the track axis, so that the infrared radiation it sent was directed at the rail foot. The holder with the scanner was each time mounted on a metal rod about 0.65 m long driven into the railway crushed stone in a way that prevented any movement of the measuring system. The measuring stations are shown in Figures 6–9.

**Figure 5.** The scanner measuring set.

**Figure 6.** Measuring station—railway viaduct—km 83,808 LK No. 272.

**Figure 7.** Measuring station—jumping change in the stiffness of the surface—access to the railway viaduct at km 88,882 LK No. 272.

**Figure 8.** Railway viaduct at km 83,808 LK No. 272—measurement of the vertical displacement of a rail dynamically loaded with a train of Koleje Wielkopolskie.

**Figure 9.** Railway viaduct at km 88,882 LK No. 272—measurement of the vertical displacement of a rail dynamically loaded with a freight train.

#### **3. Results and Discussion**

#### *3.1. General Comparisons of Calculation Results to In Situ Test Results*

The results of the measurements were compared to the results of calculations obtained from the developed theoretical model for given surfaces and vehicles that occurred during the field tests. The results are summarized in Table 3.

Differences of up to 3.5% were obtained between the results. A satisfactory convergence of the results obtained in in situ measurements with the results calculated using the developed algorithm and the finite differences method was observed. The characters of the charts are in all the cases the same. It was noticed that the curve mapping the measured vertical displacements of the rail line is more irregular, while the calculation curve is smoothed and much more symmetrical. In the case of calculations for low speeds

of the order of 30 km/h (measurement point 4), a higher density of results was noted due to the need to adopt a longer analysis time and a denser time step. In the case of measurements for high speeds, of the order of more than 100 km/h (measuring point 5) an overlap of displacements caused by adjacent axes was noted. This phenomenon is caused by the relatively short wheelbase in the bogie of the rail vehicle and the short time passing between successive loads. This issue, however, did not affect the recorded maximum value of the vertical displacement of the rail.

**Table 3.** Comparison of results obtained in in situ measurements to the theoretical results received from the developed algorithm using the finite differences method.


At the same time, factors that may have a direct impact on the results of measurements and calculations must be indicated. The accuracy of the measurements carried out is decisively influenced by the adoption of the appropriate parameters of the laser scanner, such as: length and frequency of "exposure". It is also important to fix the device in such a way that it remains in a constant position throughout the measurements.

In the context of the accuracy of calculations, the parameters of the computer used for the analysis are of grea<sup>t</sup> importance. The better they are, the more it is possible to assume a denser division of the model into nodes and time moments, which directly translates into the accuracy of calculations. In addition, it is important to adopt appropriate numerical values of data and input parameters for the railway surface. While in the case of such attributes as the mass of the rail, the modulus of longitudinal elasticity of rail steel or the geometric moment of inertia of the rail cross-section, choosing the right values does not raise doubts, in the case of the modulus of elasticity of the substrate or the damping coefficient, it is much more difficult to assume the right quantities. In this work, theoretical values determined based on the analysis of the literature on the subject were used, which, however, was confirmed by the high convergence of the results of calculations and measurements.

The influence of the dynamic impact of a rail vehicle on the railway surface within the transition zones in front of and after the engineering object, depending on its design, was analyzed.

As a starting point, a combination of a ballast and ballastless surface was adopted, where there is a step change in the parameters of the structure. A ballast surface with concrete sleepers and a good subtrack surface, as well as a ballastless surface in the rail system in the ERS type sheath, were adopted. Based on the obtained results, it is clear that within the transition zone, the dynamic effects on the surface are much greater and weaken rapidly when changing the type of structure.

For comparison with the above results, a solution was analyzed in which the transition zone in front of the engineering object was strengthened, which was reflected in the form of a 50% higher value of the elastic modulus of the substrate. Based on the obtained results, smaller differences in dynamic effects on the surface were noticed at the tested points.

The effect of the gradual change in the elasticity of the rail support within the transition zone on the magnitude of the dynamic impact of the rail vehicle on the surface was also analyzed. It was assumed that the modulus of elasticity of the support increases evenly every one meter over a length of Ln1 = 10 m in front of the object from the value characteristic for the ballast surface equal to k1.0 = 45 MPa to the value characteristic for the ballastless surface equal to k2.0 = 99 MPa. In this case, a gradual and mild decrease in the magnitude of dynamic effects on the surface in subsequent fragments of the model was noticed.

In order to compare the impact of the reinforcement of transition zones and the gradual change in the elasticity of the rail support on the reduction of the magnitude of dynamic interactions on the railway surface, the calculated results were tabularly compiled. Table 4 contains the maximum values of vertical displacement of the rail at critical points within the place of change of the type of surface (from ballast to ballastless) i.e., (1) just before the point of change and (2) just after the point of change.

**Table 4.** Comparison of three variants of the construction of the railway surface within the zone in front of the engineering facility.



By strengthening the transition zones, a 14% reduction in the magnitude of the dynamic impact on the railway surface within the site of the change of the type of surface in front of the engineering site was achieved. By applying a gradual change in the elasticity of the rail support, a 25% reduction in the dynamic effects on the railway surface within the site of the change of the type of surface in front of the engineering object was achieved.

From the analysis, it was concluded that a gradual change in the elasticity of the rail support within the transition zones of the engineering object reduces the negative impact of the dynamic impact of the rail vehicle on the railway surface and this solution is better compared to a step change in the elasticity of the support.

#### *3.2. Detailed Analysis*

Detailed results of numerical analysis are given in Figures 10–14.

## *3.3. Discussion*

The paper analyzed the dynamic impact of a rail vehicle on various solutions of the railway surface structure, with particular emphasis on the phenomenon of the threshold effect, which occurs within the transition zones of the engineering facility. The problem of locally variable stiffness of the railway surface has been identified, which in turn may lead to accelerated degradation of the structure and the need to incur increased expenditures on maintaining the infrastructure in proper condition.

Based on the analytical and numerical considerations, a computational model was created, thanks to which it is possible to determine the impact of various variants of the rail support on the dynamic response of the entire structure. The starting point for the deliberations was the Bernoulli–Euler beam. The dynamic load caused by the passage of a multiaxle rail vehicle and different vibration-damping properties of different types of surfaces are taken into account. As a consequence, the fourth-order differential equation

was obtained. It was solved by the finite differences method. A script in MATLAB was developed for a numerical solution of the problem.

**Figure 10.** The dependence of the maximum vertical displacement of the rail caused by the dynamic load, depending on the position of the node in the model in the case of a step change in the parameters of the structure.

**Figure 11.** Vertical displacement caused by dynamic load, depending on time in the case of a step change in the parameters of the structure for three different nodes: 1. node on the transition zone (ballast surface), 2. node in the place of a jumping change in the parameters of the structure (connection of the surface), and 3. node on the engineering object (ballastless surface).

**Figure 12.** Dependence of the maximum vertical displacement of the rail caused by dynamic load depending on the position of the node in the model in the case of reinforcement of transition zones.

**Figure 13.** Vertical displacement caused by dynamic load, depending on time in the case of reinforcement of transition zones for three different nodes: (1) node on the transition zone (ballast surface), (2) node in the place of change of structural parameters (surface connection), and (3) node on the engineering object (ballastless surface).

In order to verify and validate the created algorithm, in situ studies of vertical displacements of dynamically loaded railway rail were carried out. For this purpose, laser scanning technology was used.

Thus:

(1) The effectiveness of the finite differences method in the context of solving multiparameter differential equations describing the dynamics of the surface loaded by a passing railway vehicle has been confirmed. The algorithm used allowed for precise determination of the impact of technical and operational parameters on the magnitude

of dynamic interactions of a rail vehicle on the railway surface. At the same time, it should be emphasized that development of an algorithm using the finite differences method does not involve the need to use complicated and expensive computer software.


**Figure 14.** The dependence of the maximum vertical displacement of the rail caused by the dynamic load depending on the position of the node in the model in the case of a gradual change in the elasticity of the rail support within the transition zone.
