Vehicle Response to Kinematic Excitation, Numerical Simulation Versus Experiment

: The article is devoted to the numerical simulation and experimental veriﬁcation of a vehicle’s response to kinematic excitation caused by driving along an asphalt road. The source of kinematic excitation was road unevenness, which was mapped by geodetic methods. Vertical unevenness was measured in 0.25 m increments in two longitudinal proﬁles of the road spaced two meters apart with precise leveling realized by geodetic digital levels. A space multi-body computational model of a Tatra 815 heavy truck was adopted. The model had 15 degrees of freedom. Nine degrees of freedom were tangible and six degrees of freedom were intangible. The equations of motion were derived in the form of second-order ordinary differential equations and were solved numerically by the Runge–Kutta method. A custom computer program in MATLAB was created for numerical simulation of vehicle movement (eps = 2 − 52 ). The program allowed simulation of quantities such as deﬂections, speeds, accelerations at characteristic points of the vehicle, and static or dynamic components of contact forces arising between the wheel and the road. The response of the vehicle (acceleration at characteristic points) at different speeds was experimentally tested. The experiment was numerically simulated and the results were mutually compared. The basic statistical characteristics of experimentally obtained and numerically simulated signals and their power spectral densities were compared.


Introduction
A moving load's effect on transport structures has been observed in literature since 1849. After the collapse of the Chester Rail Bridge in England in the year 1947, civil engineer R. Willis and mathematician G. G. Stokes published work [1,2] in which they tried to clarify the cause of the accident and a method for solving equations of motion describing the problem. Subsequently, a flood of articles appeared that tried to solve the problems of moving loads analytically, with varying degrees of applicability. The development of computer technologies has stimulated the transition from analytical to numerical methods of a solution. Initially, the effects of moving loads were addressed in the context of dynamic analysis of railway bridges, and gradually attention was paid to road bridges. A complete overview of the procedures up to 1959 is published in [3]. An overview of the moving load effects on road bridges up to 1975 is published in [4]. The Czech and Slovak Republics are known for their high-level theoretical approach to the solution of moving load effects on transport structures. The work of V. Koloušek and L. Frýba are known all over the world [5][6][7]. A generation of younger authors, who already use numerical methods in their solutions, followed up on the work of their teachers [8][9][10][11]. Over time, attention has shifted from bridges to roads and vehicles [12][13][14]. Computational models of vehicles with different degrees of aptness, quarter-1D models, half-2D models, and spatial-3D models began to be developed. Single-layer or multi-layer models of asphalt or concrete pavements began to be created. Real unevenness of the road surface began to be taken into account. A contribution to the development was the appearance of FEM (Finite Element Method) on the scene. Finite methods have significantly expanded the possibilities of the solution. Matrix-oriented programming languages of a higher level than MATLAB [15] contributed to the development of numerical methods and the programming activities of engineers. The current situation in engineering practice is very inclined to solve interaction problems that have a theoretical and experimental foundation [16]. In general, it is stated that the research is aimed at solving increasingly complicated mathematical methods, which are used as a mathematical apparatus for solving these tasks. The final results of a theoretical task contain the large number of factors that are included in the calculation. If the problem is complicated, a detailed solution is needed to monitor and evaluate all aspects [17][18][19][20]. Monitoring, recording, and evaluating experimental measurements also form a very complicated part of the task. The process is time consuming and requires very precise work. However, it represents a possibility of verifying numerically obtained results. Therefore, the penetration of these two tasks is very valuable, and the results connected with it have a great weight in engineering and scientific practice [21][22][23].
There are many articles paying attention to numerical studies. Fewer papers are devoted to experimental studies. The contribution sought here is, as far as possible, to compare the results of numerical simulations of vehicle movement along the road with the results of experimental tests. This comparison will help engineers realize the current possibilities of applying numerical simulations in practice. It can also motivate them to look for new effective solutions. In the present state of knowledge, this work works with a combination of numerical and experimental analyses; with such a large amount of data study creates an original contribution to solving this engineering problem.

Computational Model of Vehicle
A multi-body computational model of vehicles can be created on three different qualitative levels: 1D-quarter model, 2D-plane model, and 3D-spatial model. Every model has its advantages and disadvantages, and under certain assumptions, it can be used as a solution for real practical problems [5]. The behavior of the multi-body computational model is described by a system of ordinary differential equations. The equations of motion can be derived in the sense of the component element method [24]. The relation between the components of displacements {r(t)}, corresponding to individual degrees of freedom, and deformations of joining elements {d(t)} can be expressed by the transposed static matrix [A] T : Dependence between elastic forces in joining elements (in the sense of the action of mass objects on joining elements) and its deformations is described by the equation where [k] is the diagonal stiffness matrix of joining elements. Dependence of damping forces on the velocity of deformations { . d(t)} is described by the equation where [b] is a diagonal damping matrix. The derivation according to time t is denoted by a dot above the symbol. Friction forces are considered as Resulting forces in joining elements in action on mass objects are Sign (-) is due to the principle of action and reaction. From the forces in joining elements {F JE (t)}, the static equivalents corresponding to individual degrees of freedom The gravity forces {F G } and reactions in supports {F RS (t)} must be added to the forces corresponding to individual degrees of freedom {F DF (t)}. In this manner we obtain the complete vector of forces {F R (t)} acting on the computing model of vehicle: The system of equations of motion describing the vibration of the computational model of vehicle is then expressed by the relation where [m] is a diagonal mass matrix. For this paper the spatial model of the vehicle was adopted, Figure 1. The model has 15 degrees of freedom and 10 joining elements. Tangible degrees of freedom (in number 9) correspond to movements of the mass points of the model, and intangible degrees of freedom (in number 6) correspond to movements of the contact points with the road.  {r} = [r 1 (t), r 2 (t), r 3 (t), r 4 (t), r 5 (t), r 6 (t), r 7 (t), r 8 (t), r 9 (t), u 5 (t), u 6 (t), u 7 (t), u 8 (t), u 9 (t), u 10 (t)] T (9) By mechanical multiplication of the matrices according to Equations (2) to (8), we obtain the needed Equation (11) The mass matrix [m] in Equation (11) is diagonal. The elements of the main diagonal are presented in the row matrix [m] D (12). The last six elements are zero, since the computational model has also massless degrees of freedom u 5 (t), u 6 (t), u 7 (t), u 8 (t), u 9 (t), u 10 (t).
The first nine equations in the system (11) are equations of motion and the last six equations give the expressions for the reactions of the road: ..

Road Unevenness
Road unevenness is the source of vehicle kinematic excitation. It is an integral part of the numerical solution of the vehicle response. It is possible to generate a random road profile of a certain category A-H according to the international standard (ISO 8606) [25] or to map unevenness on the selected section of the road. For this article, the geodetic method of exact leveling was used to identify road unevenness. Vertical road unevenness was measured in 0.25 m increments on a 110 m long section in two profiles, left and right, 2 m apart. The road has a small slope, 0.48%. In numerical calculations, the slope of the road is neglected. The left and right road profiles are shown in Figure 2.
Road unevenness is random. Unevenness in the monitored section on the left and right longitudinal road profile represents only one implementation of the random process. The basic statistical characteristics that characterize the longitudinal profile from a statistical point of view are given in Table 1.  It is generally accepted that vertical road unevenness satisfies Gauss's law of probability distribution. A visual view of this assumption can be made from Figure 3, where the probability distribution histograms for vertical unevenness in the left and right longitudinal road profiles are compared to Gauss's law. It should be noted that the length of the 110 m section is too short for such a comparison to be made seriously. The quality of a road surface is assessed by an international standard [7] that divides roads into 8 categories A to H according to the values of power spectral densities (PSD). The smoothed PSD S u (Ω) can be approximated in a log-log scale by a straight line. It is valid that or also where Ω is path angular frequency in [1/m], Ω 0 = 1 [1/m] is reference angular frequency, and k = 2. Figures 4 and 5 show the PSD of the left and right road profile on the log-log scale. The road can be included in category B.

Numerical Simulation of Vehicle Movement along a Road
A computer program in the MATLAB was developed for numerical simulation of vehicle movement along the road. Because it is a discrete computational model, the problem is described by a system of ordinary differential equations. The system of Equation (11) breaks down into two systems. The equations corresponding to the mass degrees of freedom in Equation (9) represent the equations of motion themselves in Equation (13), which are the subject of a numerical solution. The equations corresponding to the mass-less degrees of freedom in Equation (6) represent the equations for calculating the subsoil reactions in Equation (14). The second-order equations of motion are replaced by the first-order equations with the help of suitable substitution: The subject of the solution is then the system of the first-order equations of the type The ode45 procedure was used to solve the system of differential Equations (19) and (20). Relative error tolerance RelTol = 1 × 10 −3 , absolute error tolerance AbsTol = 1 × 10 −6 . The output of the numerical integration was realized in steps of 0.0005 s, which corresponds to a sampling frequency of 2000 Hz. At the beginning of the simulation at time t = 0, the vehicle is at a rest in a state of static equilibrium and is in contact with the road. The rear axle is at the beginning of the monitored road section. At the end of the simulation, the front axle is at the end of the monitored section. The vehicle moves at a constant speed throughout the section. The total weight of the vehicle is 21,500 kg. The damping in the vehicle tires is neglected in the simulation.
The output of the simulation is the time courses of all kinematic quantities (deflections, speeds, accelerations) in the characteristic points of the vehicle and the contact forces between the wheel of the vehicle and the road. The list of monitored quantities is as follows:

Results
To assess the degree of the correctness of the numerical simulation of the vehicle response when the vehicle moves along the road, experimental tests were realized. The experiment was performed with a Tatra T815 vehicle whose parameters were used in the numerical simulation process. During the test, vertical accelerations were monitored at three points on the right side of the vehicle: at the center of gravity (CG) of the vehicle's sprung mass-sensor, A1; on the right front axle (FA)-sensor, A2; and on the right rear axle (RA)-sensor A3. See Figure 11. Acceleration sensors Brüel & Kjaer type BK4508B were used to monitor the response. The measuring string was composed of these components: sensor, amplifier with a bandpass filter, analog-digital interface, computer. The signal from the sensor was routed through a coaxial cable. A detailed view of the sensor and its location on the front and rear axles are shown in Figure 12. The measuring string was located in the vehicle. The sampling frequency was 2000 Hz. The length of the monitored section of the road was 100 m. On this section of road, the average vehicle speed was measured as follows: steel strips and accelerometers (D2 and D3) were placed on the road in the transverse direction at a distance of 100 m, Figure 13. The average speed was calculated from the known distance and travel time. The passage of the wheel through the steel strip is manifested as a sharp impulse in the accelerometer record. By comparing the records of accelerometers D2 and D3, it is possible to determine the travel time, Figure 14. Travel time t t is determined as the difference between the end time t e and the start time t s of the passage, t t = t e − t s . The average speed of the vehicle is determined as the ratio of the section length and the passage time, v = l s /t t . The time record from the sensors on the road was synchronized with the time record from the sensors on the vehicle. A total of 12 vehicle runs were performed in the experiment at different speeds, from 15.18 km/h to 52.95 km/h.

Comparison of Numerically and Experimentally Obtained Results
During the experiment, the vehicle passed along the road at a specific speed. The same phenomenon was simulated by numerical simulation. From the records of vehicle response obtained by the experiment and by numerical simulation, sections corresponding to the vehicle movement along the selected section of road in the length of 100 m were selected. These selected sections were subjected to mutual comparison. Runs at the lowest speed V = 15.18 km/h and at the highest speed V = 52.95 km/h were chosen for comparison. The values of the vertical accelerations at three points on the right side of the vehicle were compared at the center of gravity of the vehicle's sprung mass-sensor A1, on the front axle-sensor A2, and on the rear axle-sensor A3.
First, a comparison of the monitored quantities at the vehicle speed V = 15.18 km/h is presented. Time records lasted 23.71 s and contained 47,421 samples. Figure 15 shows a visual comparison of the vertical accelerations at the center of gravity of the sprung mass of the vehicle obtained experimentally and numerically. The kinematic excitation was random, so the response of the vehicle was also random. One passage of the vehicle represents only one realization of the random process. From this point of view, it did not make sense to compare the values of the two records at each point. It was more realistic to compare the basic statistical characteristics of both records. Table 2 compares the basic statistical characteristics for experimentally and numerically obtained values of vertical accelerations at the center of gravity of the vehicle sprung mass. The density of probability distribution or distribution function provide complete information about the random variable. Therefore, in the next step, these characteristics are compared with each other, Figures 16 and 17.  Frequency spectra are needed to evaluate the dynamic behavior of structures in the frequency domain. Power spectral densities (PSD) are often used. In the next step, the PSDs were also compared with each other. The Fast Fourier Transform (FFT) was used to calculate the frequency spectra. The FFT works with the number of 2 n samples, so 2 15 = 32,768 samples from the beginning of acceleration records were selected for the analysis. As we can see from Figure 18, the dominant frequencies were in the range of 0-10 Hz, so the PSD was compared in this frequency range. The comparison of the PSDs of the vertical accelerations at the center of gravity of the vehicle sprung mass for experimentally and numerically obtained records is shown in Figure 19.  Although the frequency clusters for both records are similar, the powers of individual frequencies in the simulated record were significantly lower.
The acceleration records on the right front axle were compared similarly. Figure 20 shows a visual comparison of vertical accelerations on the right front axle.  Table 3 compares the basic statistical characteristics for experimentally and numerically obtained values of vertical accelerations on the right front axle. The density of probability distribution and distribution function are compared in Figures 21 and 22.    In the third, the vertical accelerations on the right rear axle were also compared. Figure 24 shows a visual comparison of vertical accelerations on the right rear axle.    Table 4 compares the basic statistical characteristics for experimentally and numerically obtained values of vertical accelerations on the right rear axle. The density of probability distribution and distribution function are compared in Figures 25 and 26.    The comparison of PDSs in the frequency range of 0-10 Hz is made in Figure 27. All kinematic values of the vehicle response depend on the speed of movement of the vehicle. Speed is the most important parameter entering into the vehicle-road interaction. It only makes sense to talk about the influence of all other parameters in connection with a certain specific speed of the vehicle. Therefore, in the following text, a comparison is made for a higher vehicle speed V = 52.95 km h. Time records lasted 6.7985 s and contained 13,598 samples. Figure 28 shows a visual comparison of vertical accelerations in the center of gravity of the sprung mass of the vehicle.     The acceleration record had 13,598 samples. To calculate the frequency spectrum, we needed 2 n samples, that is, either 8192 or 16,384. We decided on 16,384 samples, and the missing 2786 samples were filled with zeros. The comparison of PDSs in the frequency range 0-10 Hz is shown in Figure 31.  The acceleration record had 13,598 samples. To calculate the frequency spectrum, we needed 2 n samples, that is, either 8192 or 16,384. We decided on 16,384 samples, and the missing 2786 samples were filled with zeros. The comparison of PDSs in the frequency range 0-10 Hz is shown in Figure 31. In the next step, the vertical accelerations on the right front axle were compared. Figure 32 shows a visual comparison of vertical accelerations on the right front axle.     In the last step, the vertical accelerations on the right rear axle were compared. Figure 36 shows a visual comparison of vertical accelerations on the right rear axle.    The comparison of PDSs in the frequency range 0-10 Hz is shown in Figure 39.

Conclusions
The development of computers has created conditions for the development of simulation methods. Numerical simulation of vehicle movement along the transport structures is a real task today and can be done in real-time. This fact is confirmed by many articles dealing with this issue. Fewer articles are devoted to experimental measurements or comparisons between the results of numerical simulations and the results of experimental tests. The degree of accuracy of the various quantities obtained by numerical simulation is not the same. If differential equations are solved in the numerical simulation, where the sought functions have the meaning of deviations of certain points in time, then these quantities are calculated most accurately. The first derivatives of these functions according to time, the speed of movements, are already burdened with a larger error, and the second derivatives of these functions according to time, the acceleration of movements, are burdened with the largest error. On the contrary, when monitoring the response of a vehicle, measuring deviations is very problematic. Accelerations are measured most easily and accurately. Paradoxically, we are therefore forced to compare what we measure most accurately in the experiment with what we calculate the least accurately in the numerical simulation. When simulating the movement of a vehicle along the road, the unevenness of the road surface is a source of kinematic excitation of the vehicle, the contact between the wheel and the road is a point and changes discretely with the time step of the simulation. The source of the vehicle's kinematic excitation, in reality, is the movement of the wheel hub on the axle. This depends on the unevenness of the road, the unevenness of the tire, the pressure in the tire, and the size of the contact area between the tire and the road. The tire, therefore "smoothes out" the road profile. The real response is more settled than the response obtained by numerical simulation. It is generally accepted that road unevenness is governed by Gauss's law of distribution. If the unevenness is measured on a shorter section of the road, this may not always be evident, Figure 3. On the contrary, the response of the vehicle in all cases very well satisfies Gauss's law of distribution, Figure 40. All kinematic values of a vehicle's response depend on the speed of movement of the vehicle. Speed is the most important parameter entering into the vehicle-road interaction. It only makes sense to talk about the influence of all other parameters in connection with a certain specific speed of the vehicle. The theoretical curve (kinematic value versus vehicle speed) would have a large number of local maxima and spikes [13]. The position of the spikes is related to the discontinuous, step-changing course of the phase characteristics. The dependence of the root mean square (RMS) values for vertical accelerations in the center of gravity of the vehicle's sprung mass on the vehicle speed V used in the experiment is shown in Figure 41. It does not make sense to compare the values of the two records at each point. It is more realistic to compare the basic statistical characteristics of both records. The comparison of basic statistical characteristics looks quite satisfactory for all monitored quantities. At low vehicle speeds, it was very good. With increasing speed, the rate of agreement decreased, but it could always be considered satisfactory.
The effective value RMS is considered in dynamics as a representative quantity of a given signal. Therefore, Table 8 compares the differences (simulation-experiment) expressed in percentage (experiment represents 100%) for the observed vehicle speeds and for the three characteristic points on the vehicle. Table 8. Differences (Simulation-Experiment) of effective value RMS in (%). The biggest problem is with frequency spectra comparisons. Although the clusters of dominant frequencies in the experimental records and simulated records were similar, the powers were diametrically different. The largest disproportions related to the front axle response spectra.

Difference (Simulation-Experiment) in (%) of Experiment
Despite all the shortcomings we could mention, the methods of numerical simulation were applied in engineering practice. The advantage of numerical simulation is that it can perform interval analysis and simulate phenomena in the early stages of preparation of a structure, an analysis that could previously be obtained only experimentally on the finished structure itself; numerical simulation can also simulate phenomena that are problematic to model experimentally.
So far, an experiment is the only one way to verify the numerically obtained results. We do not consider one experiment to be sufficient. Due to the random nature of experimental measurements and the uniqueness of each passage of the vehicle on the pavement, we plan to perform another series of experimental measurements, which will provide a more extensive dataset for verification of the results obtained. This will allow the field to move the state of knowledge another step forward.