Fast calculation of vehicle-road coupled response based on moving frequency response function

Abstract

Vehicle–road coupled system is inherently time–varying, and its responses are traditionally calculated using time–domain methods which involves significant computational effort. Aiming to improve the efficiency of response calculation for the coupled system, this paper proposes a fast calculation method in frequency domain, based on the newly developed moving frequency response function (FRF). Firstly, considering the vibration characteristics of an infinitely long road, the road response is straightforwardly expressed using the road impulse response function (IRF). Subsequently, the concept of the road moving IRF is proposed and derived with respect to the moving observation points. The moving FRF is then obtained by applying Fourier transform, which allows the responses of the road moving observation points to be established in frequency domain for fast calculation under moving loads. Furthermore, by analyzing the vehicle–road coupled vibrations, based on the vehicle FRF and road moving FRF, a formula for the vehicle–road coupling force is derived in frequency domain, along with an expression for the responses at the vehicle–road contact points. Finally, the approach is illustrated in numerical simulations of vehicle–road coupled systems, and its computational efficiency and accuracy are verified through comparison with currently popular methods.

Keywords

vehicle-road coupled vibration frequency domain frequency response function impulse response function numerical simulations

Introduction

Vibration responses of vehicle–road coupled systems can provide valuable evidence for vehicle and road design (Liu et al., 2024; Zhang et al., 2022), as well as road fatigue life prediction (Pillai and Talukdar, 2023; Zhang and Cai, 2012). Study of such systems has been an important topic in road engineering for decades. It also plays a significant role in the analysis of vehicle driving quality and safety (Wang et al., 2020).

Currently, for the dynamic analysis of the vehicle–road coupled system, the equations of motion are usually constructed and employed in time domain. Common methods include the separate iterative method (SIM) (Krishnanunni and Rao, 2019; Wang et al., 2021) and the coupled degree of freedom method (CDFM) (Deng et al., 2016; Zhai et al., 2013). Yang et al. (2010) proposed a three-dimensional vehicle–pavement–foundation coupled model, where the system response was calculated through coupled equations of motion. Ding et al. (2014) simplified the pavement as a Timoshenko beam on the Pasternak foundation and the vehicle as a spring–mass–damper oscillator. The vibration equations for each subsystem were established using the mode superposition method (MSM). Krishnanunni and Rao (2019) decoupled the vehicle–pavement system into two subsystems and established separate equations of motion for the vehicle and the pavement. Ma et al. (2022) separated the vehicle–pavement coupled system into a quarter vehicle subsystem with two degrees of freedom (DOFs) and a viscoelastic multilayered asphalt pavement subsystem. The system response was solved by combining the finite-difference method and the integral transform technique. In order to improve the calculation efficiency and accuracy, Koh et al. (2003) proposed the moving element method (MEM) by discretizing the track model into elements that move together with the vehicle, eliminating the need to continuously track vehicle position. Xu et al. (2023) proposed a reduced-plate model transmission method for fast dynamic analysis of vehicle–road systems based on the MEM. All these approaches first establish the system equations of motion and then use a step-by-step time integration for solving them (Krishnanunni and Rao, 2019; Xu et al., 2023). When the vehicle or road subsystem is more complex and has a larger number of DOFs, the efficiency of solving the coupled vehicle–road response will be affected.

If the impulse response function (IRF) of the system is known, the Duhamel’s integral can be utilized to calculate the response of the system by the convolution method, which avoids solving the equations of motion and improves the calculation efficiency (Mazilu, 2007, 2010a; Mazilu et al., 2011; Sheng et al., 2020). In order to use the Green’s function for train-track coupling analysis, Metrikine and Popp (1999) used the concept of track equivalent stiffness to simplify the calculation. Mazilu (2010b) simplified the rail into an infinite-length Timoshenko beam and the vehicle into a two-mass oscillator, and Green’s function (also known as IRF) of the track was applied then to simulate the wheel–rail interaction. Sheng et al. (2016) derived the time-domain moving Green’s function of the track and established the connection with the stationary Green’s function. Zhang et al. (2020) constructed the time-domain moving Green’s function of the entire track by linear interpolation, utilizing the periodicity and symmetry of the track Green’s function. The wheel–rail force and the vibration response were obtained iteratively by establishing a vibration balance equation, which significantly improved computational efficiency.

Methods based on the frequency response function (FRF) can convert the convolution-based calculation used in the IRF method into a simple product in frequency domain, which improves the calculation efficiency. Sun and Chen (2021) proposed a time-frequency hybrid Green’s function method for dynamic analysis of periodic slab tracks. In this approach, the slab track dynamic response is calculated using the convolution of the time-domain Green’s function and the product of the frequency-domain Green’s function. Zhao and Wang (2020) utilized the Fourier transform to convert the time-domain convolution integrals into simple multiplicative calculations in the frequency domain, which reduced the difficulty of the road response calculation. Sheng et al. (2004) assumed that the stiffness and mass of the track are continuous and combined with the Fourier transform to perform train-rail coupled response calculations and spectral analysis in the frequency domain.

However, since the vehicle–road coupled system is an inherently time-varying system, existing frequency-domain methods often need to be combined with time-domain methods for response calculation, and the computational efficiency is generally low. In order to improve the calculation efficiency of the vehicle-road coupled system, this paper proposes a fast calculation method in frequency domain, based on the newly developed moving FRF, which can decouple the coupled system by directly establishing a frequency-domain product equation linking moving loads to responses at moving observation points. Furthermore, by combining the frequency response of the vehicle, the equations for calculating the coupled vehicle–road response are derived, and the vehicle-road coupled system is solved in frequency domain. It has significant advantages in the fast calculation and analysis of the vehicle-road (or train-track) coupled vibration. In practice, via the methods proposed in this study, the interaction loads of vehicle-road (train-track) can be estimated quickly, and this benefits road health monitoring and vehicle vibration control.

The remainder of this paper is structured as follows. First, a fast calculation method of road response under moving load is derived and expressed based on the proposed moving FRF. Then, this method is expanded upon to enable fast calculation of the coupled vehicle–road response. Next, the road response under a single moving load is calculated to verify and illustrate the proposed moving FRF methods. Finally, the efficient calculation of the vehicle–road coupled system is verified through numerical simulations.

Road response calculation under moving load based on moving FRF

Impulse response function (IRF) of infinitely long road with homogeneous cross-section

Firstly, as shown in Figure 1(a), let $h (x_{2}, t_{2}, x_{1}, t_{1})$ denote the road impulse response function (IRF), which represents the vertical response of road position $x_{2}$ at time $t_{2}$ to an impulse excitation $δ$ ( ${t - t}_{1}$ ) applied at road position $x_{1}$ at time $t_{1}$ . Similarly, as shown in Figure 1(b), $h (x_{2} - z, t_{2} - τ, x_{1} - z, t_{1} - τ)$ represents the vertical response of road position $(x_{2} - z)$ at time $(t_{2} - τ)$ to the impulse excitation applied at road point $(x_{1} - z)$ at time $(t_{1} - τ)$ .

The road is assumed to be infinitely long with a constant homogeneous cross-section. The structural dynamic characteristics at the two excitation points, i.e., $x_{1} and x_{1} - z$ in Figure 1(a) and (b), are the same. Considering that the respective distances between the observed points ( $x_{2} and x_{2} - z$ in Figure 1) and the excitation points ( $x_{1} and x_{1} - z$ ) are the same, the impulse responses functions $h (x_{2} - z, t_{2} - τ, x_{1} - z, t_{1} - τ)$ and $h (x_{2}, t_{2}, x_{1}, t_{1})$ equal to each other, since the same time interval $(t_{2} - t_{1})$ extends between the observation time and the excitation time:

h (x_{2} - z, t_{2} - τ, x_{1} - z, t_{1} - τ) = h (x_{2}, t_{2}, x_{1}, t_{1})

(1)

As shown in Figure 1(c), let the impulse excitation $δ (t)$ apply at the coordinate origin at the origin time. The corresponding response at the point $(x_{2} - x_{1})$ at time $(t_{2} - t_{1})$ is denoted as $h_{0} (x_{2} - x_{1}, t_{2} - t_{1})$ , that is, $h_{0} (x_{2} - x_{1}, t_{2} - t_{1}) = h (x_{2} - x_{1}, t_{2} - t_{1}, 0, 0)$ . In equation (1), let $z = x_{1}$ , $τ = t_{1}$ . Equation (1) can be then simplified:

h_{0} (x_{2} - x_{1}, t_{2} - t_{1}) = h (x_{2}, t_{2}, x_{1}, t_{1})

(2)

Therefore, for an infinitely long road with a homogeneous cross-section, the IRF can be simplified from a four-argument function to a two-argument function. In other words, the IRF $h (x_{2}, t_{2}, x_{1}, t_{1})$ of any point on the road can be calculated by the IRF at a specific point $h_{0} (x_{2} - x_{1}, t_{2} - t_{1})$ . This observation can be used to improve the calculation efficiency of the IRF.

Response calculation for moving observation point based on moving impulse response

Let $y (x, t)$ denote the road response at the observation point $x$ to a load $f (t)$ that moves with a constant velocity $v$ , as shown in Figure 2. According to Duhamel’s integral, the response $y (x, t)$ can be calculated.

y (x, t) = \int_{0}^{t} h (x, t, v τ, τ) f (τ) d τ

(3)

where

h (x, t, v τ, τ)

is the response of point

x

at time

t

to an impulse excitation

δ (t - τ)

applied at point

z_{p} (τ) = v τ

at time

τ

. By combining equations (2) and (3), the expression of

y (x, t)

can be further simplified as:

y (x, t) = \int_{0}^{t} h_{0} (x - v τ, t - τ) f (τ) d τ

(4)

Assume a set of loads move on a road with a constant velocity $v$ . Let $f_{v i} (t)$ denote the $i$ -th moving load with its initial coordinate $x_{i}$ along the road. Via equation (2), the IRF can be obtained as $h_{0} (x - x_{i} - v τ, t - τ) = h (x, t, x_{i} + v τ, τ)$ . The road response $y (x, t)$ of point $x$ to all moving loads can be expressed as:

y (x, t) = \sum_{i = 1}^{n} \int_{0}^{t} h_{0} (x - x_{i} - v τ, t - τ) f_{v i} (τ) d τ

(5)

where

n

is the number of the moving loads.

In practice, road responses at positions of moving loads are often concerned. Assume that the response observation point moves with the same velocity $v$ as the loads. Let $l$ denote the initial position of the observation point, and substitute the coordinate of the moving observation point $x (t) = l + v t$ into equation (5). This yields the following equation.

y (l + v t, t) = \sum_{i = 1}^{n} \int_{0}^{t} h_{0} (l - x_{i} + v (t - τ), t - τ) f_{v i} (τ) d τ

(6)

In order to further simplify equation (6), let ${h_{v} (l - x_{i}, t - τ) = h}_{0} (l - x_{i} + v (t - τ), t - τ)$ and $y_{v} (l, t) = y (l + v t, t)$ . The vertical response of the road point $x (t) = l + v t$ can be then rewritten as:

y_{v} (l, t) = \sum_{i = 1}^{n} \int_{0}^{t} h_{v} (l - x_{i}, t - τ) f_{v i} (τ) d τ

(7)

Since the observation points in the IRF $h_{v} (l, t)$ are moving, the $h_{v} (l, t)$ presented in this paper is called moving IRF.

Road response calculation based on moving FRF

In the moving IRF $h_{v} (l - x_{i}, t)$ in equation (7), let all the parameters (that is, $l$ , $x_{i}$ , and $v$ ) besides $t$ be constant in time. The Fourier transform can be applied on both sides of equation (7), which yields the frequency-domain road response ${\tilde{y}}_{v} (l, ω)$ as shown in the following equation.

{\tilde{y}}_{v} (l, ω) = \sum_{i = 1}^{n} {\tilde{h}}_{v} (l - x_{i}, ω) {\tilde{f}}_{v i} (ω)

(8)

As a result, the time-domain convolution of equation (7) is decoupled and represented in frequency domain as a simple product, which improves the calculation efficiency. The Fourier transform ${\tilde{h}}_{v}$ of the moving IRF $h_{v}$ is further on called the moving FRF in this paper.

According to equation (8), the response ${\tilde{y}}_{v} (x_{j}, ω)$ at the position of the $j$ -th moving load is expressed as:

{\tilde{y}}_{v} (x_{j}, ω) = \sum_{i = 1}^{n} {\tilde{h}}_{v} (x_{j} - x_{i}, ω) {\tilde{f}}_{v i} (ω)

(9)

Therefore, if the observation points move and coincide with the $n$ moving load points, their coordinates can be successively substituted into equation (9). This yields the vector of responses ${\tilde{Y}}_{v} (ω)$ as shown in the following equation.

{\tilde{Y}}_{v} (ω) = {\tilde{H}}_{v} (ω) {\tilde{F}}_{v} (ω)

(10)

The following equation presents the specific expressions for the response matrix ${\tilde{Y}}_{v} (ω)$ , the moving FRF matrix ${\tilde{H}}_{v} (ω)$ , and load matrix ${\tilde{F}}_{v} (ω)$ .

{\begin{cases} {\tilde{Y}}_{v} (ω) = {[{\tilde{y}}_{v} (x_{1}, ω), {\tilde{y}}_{v} (x_{2}, ω), \dots, {\tilde{y}}_{v} (x_{n}, ω)]}^{T} \\ {\tilde{F}}_{v} (ω) = {[{\tilde{f}}_{v 1} (ω), {\tilde{f}}_{v 2} (ω), \dots, {\tilde{f}}_{v n} (ω)]}^{T} \\ {\tilde{H}}_{v} (ω) = [\begin{array}{l} \begin{array}{l} {\tilde{h}}_{v} (0, ω) & {\tilde{h}}_{v} (x_{1} - x_{2}, ω) \\ {\tilde{h}}_{v} (x_{2} - x_{1}, ω) & {\tilde{h}}_{v} (0, ω) \end{array} & \begin{array}{l} \dots & {\tilde{h}}_{v} (x_{1} - x_{n}, ω) \\ \dots & {\tilde{h}}_{v} (x_{2} - x_{n}, ω) \end{array} \\ \begin{array}{c} ⋮ & ⋮ \\ {\tilde{h}}_{v} (x_{n} - x_{1}, ω) & {\tilde{h}}_{v} (x_{n} - x_{2}, ω) \end{array} & \begin{array}{c} ⋱ & ⋮ \\ \dots & {\tilde{h}}_{v} (0, ω) \end{array} \end{array}] \end{cases}

(11)

Calculation of vehicle–road coupled response based on moving FRF

Coupling force between vehicle and road

For vehicle–road coupling analysis, various mass–spring vehicle models with different numbers of DOFs are often employed, such as a two-DOF quarter vehicle model (Qin et al., 2019), four-DOF half vehicle model (Zhang et al., 2023) or a whole vehicle model with several DOFs (Chen et al., 2020; Jian et al., 2022; Zhu et al., 2024). Here the vehicle is considered to be a mass–spring model with $n_{c}$ DOFs, as shown in Figure 3, and with $n$ contact points between the vehicle and road. Let the matrices $M$ , $C$ , and $K$ denote respectively the mass, damping, and stiffness matrix of the vehicle model. Denote by $k_{i}$ and $c_{i}$ the tire spring stiffness and damping coefficient at the $i$ -th contact point, and let the vector $U_{c} (t)$ represent the vehicle responses along its DOFs.

Assume that the vehicle moves on the road with a constant velocity $v$ , as shown in Figure 3. Let $x_{c i}$ represent the initial coordinate of the $i$ -th contact point between the vehicle and road, so that the coordinate of the $i$ -th contact point at time $t$ is $x_{v i} (t) = x_{c i} + v t$ . Denote by $r_{v i} (t)$ , $y_{v i} (t)$ , and $f_{v i} (t)$ respectively the corresponding road roughness, the vertical road response, and the coupling force between the vehicle and road at time $t$ .

Let $B_{i}$ be the point of the vehicle that corresponds to the $i$ -th contact position between the vehicle and road, as shown in Figure 3. Then the response of the point $B_{i}$ equals the sum of the displacement and roughness of the road at the $i$ -th contact point, that is, $u_{b i} (t) = y_{v i} (t) + r_{v i} (t)$ . In dynamic analysis, the analyzed excitation applied to the structure does not include the gravity since it has been balanced at the structural equilibrium state. Let $f_{b i} (t)$ be the contact load applied on the vehicle at the contact point $B_{i}$ , and thus $f_{b i} (t) = - f_{v i} (t) - f_{g i}$ , where $f_{g i}$ is the vehicle gravity. At each contact position, the above parameters are collected and arranged as vectors, and then expressed as in the following equations.

F_{v} (t) = {- F}_{b} (t) - F_{g}

(12)

U_{b} (t) = Y_{v} (t) + R_{v} (t)

(13)

{\begin{cases} \begin{array}{l} R_{v} (t) = {[r_{1} (t), r_{2} (t), \dots, r_{n} (t)]}^{T} \\ Y_{v} (t) = {[y_{v 1} (t), y_{v 2} (t), \dots, y_{v n} (t)]}^{T} \end{array} \\ U_{b} (t) = {[u_{b 1} (t), u_{b 2} (t), \dots, u_{b n} (t)]}^{T} \\ F_{v} (t) = {[f_{v 1} (t), f_{v 2} (t), \dots, f_{v n} (t)]}^{T} \\ F_{b} (t) = {[f_{b 1} (t), f_{b 2} (t), \dots, f_{b n} (t)]}^{T} \\ F_{g} = {[f_{g 1}, f_{g 2}, \dots, f_{g n}]}^{T} \end{cases}

(14)

Vehicle equations of motion

In order to analyze the vehicle dynamics, the vehicle separated from the road system. The force acting on the contact point $A_{i}$ between the $i$ -th wheel spring and the vehicle body equals ${k_{i} u}_{b i} (t) + c_{i} {\dot{u}}_{b i} (t)$ . The vehicle equations of motion can be written as:

M {\ddot{U}}_{c} (t) + C {\dot{U}}_{c} (t) + K U_{c} (t) = I_{b} (K_{b} U_{b} (t) + C_{b} {\dot{U}}_{b} (t))

(15)

where

K_{b}

is an

n \times n

diagonal matrix with the wheel spring stiffness coefficients

k_{1}, k_{2}, \dots, k_{n}

on the diagonal,

K_{b} (i, i) = k_{i}

, and the structure of the wheel damping matrix

C_{b}

is similar to that of the matrix

K_{b}

. The matrix

I_{b}

is the position matrix of wheel loads with

I_{b} (i_{b}, i) = 1

i = 1, 2, \dots, n

, where

i_{b}

represents the number of DOF of the contact point

A_{i}

. Let

u_{{c i}_{b}} (t)

and

u_{b i} (t)

respectively denote the displacement of the points

A_{i}

and

B_{i}

. Then, the force applied at the

i

-th vehicle contact point

B_{i}

can be expressed as

f_{b i} (t) = c_{i} {\dot{u}}_{b i} (t) + k_{i} u_{b i} (t) - c_{i} {\dot{u}}_{{c i}_{b}} (t) - k_{i} u_{{c i}_{b}} (t)

. Collecting this formula for all the springs and dampers yields the matrix form of the following equation.

C_{b} {\dot{U}}_{b} (t) + K_{b} U_{b} (t) - C_{b} I_{b}^{T} {\dot{U}}_{c} (t) - K_{b} I_{b}^{T} U_{c} (t) = F_{b} (t)

(16)

Fast calculation of vehicle–road coupled response based on moving FRF

The Fourier transform applied on both sides of equations (15) and (16) yields the frequency-domain expressions shown as:

{\begin{cases} - ω^{2} {M \tilde{U}}_{c} (ω) + j ω C {\tilde{U}}_{c} (ω) + K {\tilde{U}}_{c} (ω) = I_{b} (K_{b} {+ j ω C}_{b}) {\tilde{U}}_{b} (ω) \\ j ω C_{b} {\tilde{U}}_{b} (ω) + K_{b} {\tilde{U}}_{b} (ω) - {j ω C}_{b} I_{b}^{T} {\tilde{U}}_{c} (ω) - K_{b} I_{b}^{T} {\tilde{U}}_{c} (ω) = {\tilde{F}}_{b} (ω) \end{cases}

(17)

The vector ${\tilde{U}}_{b} (ω)$ collects the responses of the vehicle contact points, and ${\tilde{U}}_{c} (ω)$ represents the responses of the vehicle body. They both can be calculated using the contact load ${\tilde{F}}_{b} (ω)$ , as expressed in the following equations.

{\tilde{U}}_{b} (ω) = {\tilde{H}}_{b} (ω) {\tilde{F}}_{b} (ω)

(18)

{\tilde{U}}_{c} (ω) = {\tilde{H}}_{c b} (ω) {\tilde{F}}_{b} (ω)

(19)

where the matrix

{\tilde{H}}_{b} (ω)

represents the FRFs of the vehicle contact points to the contact load, and the matrix

{\tilde{H}}_{c b} (ω)

represents the FRFs of the points on the vehicle body to the contact load. Via the modal analysis, they can be expressed as:

{\begin{cases} {\tilde{H}}_{c} (ω) = {(- ω^{2} M + j ω C + K)}^{- 1} \\ {\tilde{H}}_{b} (ω) = {((K_{b} + {j ω C}_{b}) - ({j ω C}_{b} + K_{b}) I_{b}^{T} {\tilde{H}}_{c} (ω) I_{b} (K_{b} + {j ω C}_{b}))}^{- 1} \\ {\tilde{H}}_{c b} (ω) = {\tilde{H}}_{c} (ω) I_{b} ({K_{b} + j ω C}_{b}) {\tilde{H}}_{b} (ω) \end{cases}

(20)

The matrix ${\tilde{H}}_{c} (ω)$ is the FRF of the vehicle body, which is also expressed in equation (20). The Fourier transform applied on the two sides of equations (12) and (13) yields the frequency-domain expressions shown as:

{\begin{cases} {\tilde{U}}_{b} (ω) = {\tilde{Y}}_{v} (ω) + {\tilde{R}}_{v} (ω) \\ {\tilde{F}}_{v} (ω) = - {\tilde{F}}_{b} (ω) - F_{g} (ω) \end{cases}

(21)

By combining equations (10), (18), and (21), the vehicle contact load ${\tilde{F}}_{b} (ω)$ can be obtained:

{\tilde{F}}_{b} (ω) = - {({\tilde{H}}_{b} (ω) + {\tilde{H}}_{v} (ω))}^{- 1} {\tilde{H}}_{v} (ω) F_{g} (ω) + {({\tilde{H}}_{b} (ω) + {\tilde{H}}_{v} (ω))}^{- 1} {\tilde{R}}_{v} (ω)

(22)

Substituting ${\tilde{F}}_{b} (ω)$ into equation (19) leads to the vehicle response ${\tilde{U}}_{c} (ω)$ . Thereupon, by substituting the expression for the road moving load ${\tilde{F}}_{v} (ω)$ shown in equation (21) into equation (10), the response ${\tilde{Y}}_{v} (ω)$ of the road contact point is obtained:

{\tilde{Y}}_{v} (ω) = - {\tilde{H}}_{v} (ω) ({\tilde{F}}_{b} (ω) + F_{g} (ω))

(23)

Finally, the time history of the moving load $F_{v} (t)$ can be recovered via the inverse Fourier transform of its frequency-domain form ${\tilde{F}}_{v} (ω)$ and used to obtain the time-domain response of the road contact point by equation (5).

The vehicle–road coupled system is inherently time-varying, and the system responses are usually calculated using time-domain methods, which is generally time consuming. This paper proposes a novel calculation method in frequency domain based on the moving FRF, where the computationally demanding time-domain deconvolution is transformed into a much simpler multiplication in frequency domain. This significantly improves the calculation efficiency of the coupled vehicle–road responses.

Numerical calculation of road responses under moving loads based on moving FRF

In order to verify the proposed methods, a numerical simulation is performed in this section to confirm the applicability and efficiency of the proposed road response calculation procedure based on the moving FRF.

Road model

As shown in Figure 4, the analyzed road is modeled as an Euler–Bernoulli beam with viscoelastic foundation and simple supports at both ends.

Let $E$ , $c_{s}$ , and $I$ denote respectively the road elastic modulus, the damping factor, and the moment of inertia of the road cross-section. The mass per meter of the road is $m$ . For the foundation, the linear elastic modulus is $k$ , and the foundation damping coefficient is $c$ . The specific values of the road model parameters are reused from Refs. (Lv et al., 2010; Yang et al., 2010, 2013) and listed in Figure 4. Regard to the parameters in Figure 4, the displacement of the road under typical moving loads practically vanishes at the distance of 10 m. Therefore, the beam is $L =$ 40 m long, and all the load and observation points are modeled in its mid part (in-between 10 m and 30 m), so that the simulated responses are equivalent to those of an infinitely long road.

Calculation of moving FRF

The road moving FRF is the frequency response of the moving observation point to an impulse excitation. It is difficult to derive analytically, and in this paper, the time-domain moving IRF is first calculated via the road finite element (FE) model, and then the Fourier transform is used to obtain the moving FRF.

The road model with a length of 40 m is established for calculation of the moving IRF when the velocity is 10 m/s. Figure 5(a) shows the moving IRF $h_{v} (0, t)$ of the road calculated using the Newmark integration method. The Fourier transform of the road IRF is the moving FRF ${\tilde{h}}_{v} (0, ω)$ , and its real and imaginary parts, as well as the absolute value are all shown in Figure 5(b).

Figure 6(a) shows the variation of the amplitude $| {\tilde{h}}_{v} (0, ω) |$ of the moving FRF under different velocities. It can be seen that the frequency spectrum of the moving FRF becomes wider and the amplitude gradually decreases as the velocity increases. Figure 6(b) shows the variation of the amplitude of the moving FRF with the distance $d$ under the moving velocity 10 m/s, where $d$ is the initial distance between the observation point and the impulse excitation. In case of $d$ < 0, the amplitude of moving FRF increases with distance $d$ increasing. In case of $d$ > 0, the amplitude gradually decreases.

Road response calculation

In this section, a single moving load is considered to illustrate and verify the applicability of the moving FRF method for road response calculation. Let both the observation point and the load move at a velocity of 10 m/s, and let the considered time interval be 10s. The sampling time step is 0.001s, i.e., the sampling frequency is $f_{s}$ = 1000 Hz.

Firstly, the road moving IRF $h_{v} (0, t)$ is calculated via the road FE model, as shown in Figure 5(a), and the corresponding moving FRF ${\tilde{h}}_{v} (0, ω)$ is obtained in frequency domain, as shown in Figure 5(b). Then, the single moving random load $f_{v} (t)$ is constructed as shown in Figure 7, with its frequency-domain amplitude $| {\tilde{f}}_{v} (ω) |$ obtained through the Fourier transform and shown in Figure 8(a). Next, the road frequency response ${\tilde{y}}_{v} (0, ω)$ is obtained via equation (10). Finally, the time-domain response $y_{v} (t)$ is obtained through the inverse Fourier transform of the frequency response ${\tilde{y}}_{v} (0, ω)$ . The frequency spectrum and the time-history curve of the calculated road response $y_{v} (0, t)$ are respectively shown in Figure 8(b) and Figure 9 (“Moving FRF”).

In order to verify the accuracy of the road response obtained using the moving FRF, the response $y_{v} (0, t)$ is additionally calculated directly in time domain using the Newmark integration method and Duhamel’s integral, and the results are shown in Figure 9 as the lines labelled “Newmark” and “Duhamel”. All three responses are in a good agreement.

Effect of frequency truncation

As seen in Figure 5(b), in the range beyond 17 Hz, the amplitude of frequency response is a strictly decreasing function of the frequency. Above a certain threshold frequency, the amplitude is small enough that its contribution to the overall response can be ignored. Denoting the truncation frequency by $\bar{f}$ , the road response can be approximated by $y_{v} (0, \bar{f}, t)$ to improve the computational efficiency.

y_{v} (0, \bar{f}, t) = \frac{1}{2 π} \int_{0}^{\bar{f}} {\tilde{y}}_{v} (0, ω) e^{j ω t} d ω

(24)

The responses computed with the truncation frequencies $\bar{f}$ equal to 17 Hz and 70 Hz are plotted in Figure 10. The third response shown in Figure 10, labeled “Theoretical solution”, is calculated using the Newmark method and defined here as $y_{v} (t)$ .

When the truncation frequency is 17 Hz, the calculated response is close to the theoretical solution, but its deviations are visible. When the truncation frequency is 70 Hz, the calculated response is indistinguishable from the theoretical solution. The effect of the truncation frequency on the calculation accuracy can be quantified by the relative error $ε (\bar{f})$ :

ε (\bar{f}) = \frac{‖ y_{v} (0, \bar{f}, t) - y_{v} (t) ‖}{‖ y_{v} (t) ‖}

(25)

The respective curve is shown in Figure 11. As expected, the error decreases as the truncation frequency increases, and it falls below 10⁻⁴ for the truncation frequencies above 70 Hz. Therefore, the frequency range of 0–70 Hz is used in the further numerical simulations to calculate the response.

Computational efficiency

The computational efficiency of the proposed methods can be analyzed as follows. For the application of the Newmark integration method, the road model length is set to 120 m to avoid the influence of the boundary support. The model is divided into 240 elements, and the total number of DOFs of the road model $N$ is 480. In the Newmark method, the number of multiplications for each time step is O $(N^{3})$ , and the number of time steps is $f_{s} \times t =$ 1000 $\times$ 10. Therefore, the total number of multiplications is O $(N^{3} f_{s} t)$ , where $N^{3} f_{s} t = 10^{12}$ .

In case of using Duhamel’s integral (Rózanski et al., 2021) via equation (7), the time length $t_{1}$ of the impulse response shown in Figure 5(a) is 2s, and the length of the impulse response is $f_{s} \times t_{1} =$ 1000 $\times$ 2, while the length of the load is $f_{s} \times t =$ 1000 $\times$ 10. According to the principle that the system response is equal to the superposition of the response under load applied at each time moment, the number of multiplications in Duhamel’s integral method is ${(f_{s} \times t_{1})}^{2} \times f_{s} \times t = {(1000 \times 2)}^{2} \times 1000 \times 10 = 4 \times 10^{10}$ .

In case of the moving FRF method, the number of multiplications (Burova and Ryapukhin, 2021) in the frequency-domain response calculation, according to equation (10), is 4 $\times \bar{f} / ∆ f$ , where $∆ f = 1 / t = 0.1$ is the frequency resolution, which amounts to 2800 multiplications. Additionally, the cost of the FFT and the IFFT should be added, which are each equal to $N / 2 \times \log_{2} N$ , where $N$ is the length of the data. The FFT is utilized to calculate ${\tilde{H}}_{v} (ω)$ from $H_{v} (t)$ , and the data length is $f_{s} \times t_{1} =$ 1000 $\times$ 2, which amounts to $1.10 \times 10^{4}$ multiplications. The IFFT is used to calculate $y_{v} (0, t)$ from ${\tilde{y}}_{v} (0, ω)$ , and the data length is 1000 $\times$ 10, which amounts to $6.64 \times 10^{4}$ multiplications. The total number of multiplications for the moving FRF method is thus $8.02 \times 10^{4}$ .

The computational costs of the three above methods are listed and compared in Table 1. It can be seen that the moving FRF method only costs 0.00,496 s to compute the road response under the moving loads, of which the computational efficiency is much higher than that of the other two methods.

Table 1.

Number of multiplication operations and computational cost for different methods.

Method	Type	Number of multiplication operations	Computational cost (s)
Newmark	Time domain	O (10¹²)	1.167
Duhamel	Time domain	$4 \times 10^{10}$	0.0877
Moving FRF	Frequency domain	$8.02 \times 10^{4}$	0.00496

Numerical simulation of coupled vehicle–road response

In this section, a vehicle–road coupled system is built and analyzed numerically to illustrate and verify the proposed fast calculation method of the system response.

Vehicle model

In this paper, a four-DOF half vehicle model is used as an example to analyze the coupled vehicle–road vibrations. Figure 12 shows the vehicle model, while the employed road model has been described in Figure 4. Table 2 lists the values of vehicle parameters which are the vehicle body mass $m_{v}$ , rotational inertia $J_{v}$ , masses of wheels $m_{1}$ and $m_{2}$ , suspension damping $c_{3}$ and $c_{4}$ , suspension stiffness $k_{3}$ and $k_{4}$ , tire damping $c_{1}$ and $c_{2}$ , tire stiffness $k_{1}$ and $k_{2}$ , as well as horizontal distances from the centroid of the vehicle to the axles $e_{1}$ and $e_{2}$ . The specific values of these parameters are borrowed from Ref. (Fu, 2015).

Table 2.

Parameters of vehicle model.

$m_{v}$ (kg)	$J_{v}$ (kg∙m²)	$c_{1}$ (N∙s/m)	$c_{2}$ (N∙s/m)	$c_{3}$ (N∙s/m)	$m_{1}$ (kg)	$m_{2}$ (kg)
25,000	1.47 × 10⁵	3.9 × 10⁴	4.3 × 10³	3 × 10⁴	1500	1000
$c_{4}$ (N∙s/m)	$k_{1}$ (N/m)	$k_{2}$ (N/m)	$k_{3}$ (N/m)	$k_{4}$ (N/m)	$e_{1}$ (m)	$e_{2}$ (m)
4 × 10⁴	3.74 × 10⁶	4.6 × 10⁶	2.47 × 10⁶	4.23 × 10⁶	2.22	2.05

In Figure 12, let $B_{1}$ and $B_{2}$ denote the contact points between the vehicle wheels and the road. Let $u_{1} (t)$ and $u_{2} (t)$ be the vertical displacement and the rotation angle of the vehicle body, respectively, and $u_{3} (t)$ and $u_{4} (t)$ denote the vertical displacement of the vehicle front axle and rear axle, respectively. Detailed formulas for the vehicle mass matrix $M$ , damping matrix $C$ , and stiffness matrix $K$ can be found in Ref. (Li et al., 2023).

Road roughness

The road roughness $r (x)$ is modeled in this paper as a superposition of a series of trigonometric functions via the power spectrum density (PSD) function, see Ref. (Zhang et al., 2021) and ISO 8086, as shown in the following equations.

r (x) = \sum_{k = 1}^{N} 2 \sqrt{G_{q} (n_{k}) ∆ n} \cos (2 π n_{k} x + θ_{k})

(26)

G_{q} (n_{k}) = G_{q} (n_{0}) {(n_{k} / n_{0})}^{- 2}

(27)

where

x

is the road coordinate and

θ_{k}

is a random phase angle uniformly distributed in [0 2

π

]. The displacement PSD is denoted by

G_{q} (n_{k})

, and

G_{q} (n_{0})

is the coefficient of roughness, where

n_{0}

is the reference spatial frequency and

n_{k}

is the spatial frequency.

The road surface grade is set to class A. The constructed roughness, shown in Figure 13, is used in the subsequent vibration analysis of the vehicle–road coupled system.

Vehicle–road coupled response

The considered road length is 240 m, and the vehicle moves with the velocity of 10 m/s. The road moving FRF is constructed via the derivation in the second Section. The moving IRFs are calculated as descripted in equation (7) by applying a unit impulse excitation on the road FE model with a length of 40 m, and then the moving FRFs are obtained by the Fourier transform. The moving FRF matrix is at last constructed according to the form shown in equation (11). Figure 14(a) shows the frequency spectrum of the moving FRFs ${\tilde{h}}_{v} (0, ω)$ , ${\tilde{h}}_{v} (4.27, ω)$ , and ${\tilde{h}}_{v} (- 4.27, ω)$ , where $4.27 m = e_{1} + e_{2}$ is the distance between the vehicle axles. The peaks of the spectra are all close to the fundamental road frequency $f_{c 1}$ . Secondly, the moving FRF matrix ${\tilde{H}}_{v} (ω)$ is calculated in frequency domain by equation (11), as well as the frequency responses ${\tilde{H}}_{b} (ω)$ of the vehicle contact points. Then, the contact loads ${\tilde{F}}_{b} (ω)$ can be obtained by equation (22). The frequency spectra of the contact loads ${\tilde{F}}_{b} (ω)$ at the vehicle contact point are plotted in Figure 14(b). At last, using the computed contact loads ${\tilde{F}}_{b} (ω)$ and the moving FRF matrix ${\tilde{H}}_{v} (ω)$ , the frequency-domain response ${\tilde{Y}}_{v} (ω)$ of the road contact points is computed by equation (23). Finally, the corresponding time histories $Y_{v} (t)$ of the road contact points are obtained by the inverse Fourier transform of ${\tilde{Y}}_{v} (ω)$ and plotted in Figure 15.

For comparison purposes, the road contact responses were calculated also by the CDFM (Deng et al., 2016; Zhai et al., 2013) to verify the accuracy of the proposed moving FRF method. The responses obtained using both methods are compared in Figure 15. Let $ε_{i}$ denote the relative error of the $i$ th road contact point response calculated by the CDFM with respect to the moving FRF method. The values of $ε_{1}$ and $ε_{2}$ , calculated according to equation (25), equal to 0.3% and 0.5%, respectively. Thus, the accuracy of the moving FRF method for calculation of the coupled vehicle–road response is initially verified.

Figure 16 plots the corresponding frequency spectra of vehicle acceleration, and Figure 17 shows the time histories of vehicle acceleration obtained via the inverse Fourier transform of u¨∼₃(ω) and u¨∼₄(ω).

Comparison and analysis

In order to further verify the accuracy and efficiency of the proposed moving FRF method, vertical acceleration responses of the wheels are additionally calculated by the CDFM (Yang et al., 2010) and SIM (Krishnanunni and Rao, 2019). The results are plotted in Figure 17 as “CDFM” and “SIM” and compared with the response obtained using the moving FRF method. All three responses are highly consistent with each other.

Furthermore, a 3D FE model of the vehicle-road coupled system, as shown in Figure 18, is established for fast computation of the system responses. The system parameters are consistent with the 2D model parameters. Specifically, the wheelbase between the left and right wheels of the 3D vehicle model is set as 2 m, while the width of the road model is 3 m. Assume that the road roughness experienced by the left and right wheels are in same condition.

The displacement $u_{1} (t)$ , velocity ${\dot{u}}_{1} (t)$ , and acceleration ${\ddot{u}}_{1} (t)$ of the vehicle body in the vertical direction using different methodologies are depicted in Figure 19, and it can be seen that the results obtained using the moving FRF method, CDFM and the SIM are essentially identical, which further validates the applicability of the proposed method in this paper.

In order to quality the efficiency improvement, the 2D model is taken as an example to calculate response of the coupled vehicle-road system. In case of using the CDFM method, the matrix dimension of the coupled system is 966 × 966. For the SIM, the dimension of the road system matrix is 960 × 960, and the dimension of the vehicle system matrix are 4 × 4. While for the moving FRF method, the dimension of the constructed FRF matrix in frequency is only 2 × 2, which depends solely on the number of contact points. For the moving FRF method, the main computational cost is the calculation of the moving FRF, as the method does not need to assemble and utilize any time-varying matrices of the coupled system. In addition, the computation is performed in frequency domain, which avoids time-consuming iterative calculations. The responses can be obtained directly based on the frequency-domain vibration information of the basic road model. In detail, for the CDFM and SIM, the total number of computational steps required is equal to

f_{s} \times t

, whereas for the moving FRF method the number of computational steps in frequency-domain is equal to

\bar{f} \times t

. If the sampling frequency

f_{s}

is set at 1000 Hz, the truncated frequency

\bar{f}

is just 70 Hz. Therefore, considering both the matrix dimensions and the number of computational steps, the computing work performed by the moving FRF method is much less than that of the CDFM and the SIM. Table 3 lists and compares the advantages and disadvantages of the proposed method, as well as the CDFM and the SIM.

Table 3.

Comparison of vehicle-road coupled calculation methods.

Method	Advantage	Disadvantage
CDFM	High computational accuracy	Large dimension of the system matrix
CDFM	Effectively application on nonlinear systems	High computational cost
SIM	Directly solving of the vehicle-road interaction forces	Repeated iterative calculations with low computational efficiency
Moving FRF	Decoupling the time-variant system into time-invariant system	Current limitation on nonlinear problems
Moving FRF	Fast calculation in frequency domain	Current limitation on nonlinear problems

To further compare the computational cost, the above three methods are performed on a same computer, and the time spent on the main computational processes of the CDFM, the SIM, and the moving FRF method is 1610 s, 686 s, and 5.9 s, respectively. It confirms that the moving FRF method has a significant advantage in fast response calculation of vehicle–road coupled systems.

Conclusions

Aiming to improve computational efficiency of vehicle–road coupled systems, this paper has proposed a fast response calculation method based on the moving FRF in frequency domain. The effectiveness of the proposed method is verified in numerical simulations, and the main conclusions are as follows:

• A vehicle–road coupled system is inherently time-varying, and the calculation of its coupled response is often a complicated process. The proposed moving FRF method establishes simple expressions for vehicle–road coupling calculations in frequency domain. It is conducive to analyses of the relationship between the vehicle response, road response, contact point response and road roughness.

• Traditional methods for calculation of vehicle–road coupled systems are often time-consuming and the computation work is often performed in time-domain. Via the proposed moving FRF, the time-variant system can be converted into time-invariant system, and the system responses can be solved and expressed in frequency domain in terms of the moving FRF and the coupling force with high computational efficiency.

• It is usually difficult to balance the computational efficiency and accuracy of vehicle–road coupled systems. By comparing with the widely used methods of CDFM and SIM, the proposed moving FRF method is verified to be simultaneously accurate and efficient.

Even though encouraging results have been obtained using the moving FRF method, it still has certain limitations in terms of nonlinear analysis of vehicle-road system. Aiming at the limitation, he method proposed in this study is going to investigate the stochastic and nonlinear characteristics of the vehicle-road coupled system based on the system motion equation. For example, take nonlinear elastic forces to equivalent the stiffness nonlinear parts in vehicles and roads, and linearize nonlinear forces through first-order Taylor expansion.

Footnotes

Acknowledgments

The authors would like to acknowledge the support of the National Natural Science Foundation of China (NSFC) (52378285), and of the China - Central and Eastern European Countries University Joint Education Project (2022206), and of the Liaoning Provincial Natural Science Foundation of China (2024-MS-020).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (NSFC) (52378285), and of the China - Central and Eastern European Countries University Joint Education Project (2022206), and of the Liaoning Provincial Natural Science Foundation of China (2024-MS-020).

ORCID iD

Jilin Hou

References

Burova

Ryapukhin

(2021) Reduction of the number of multiplication operations in digital signal processing algorithms by classical methods of discrete Fourier transform. AIP Conference Proceedings 2402: 040002.

Chen

Zhang

Wang

(2020) Rigid–flexible coupled dynamic response of steel–concrete bridges on expressways considering vehicle–road–bridge interaction. Advances in Structural Engineering 23(1): 160–173.

Deng

Cao

Wang

, et al. (2016) A multi-point tire model for studying bridge-vehicle coupled vibration. International Journal of Structural Stability and Dynamics 16(8): 1550047.

Ding

Yang

Chen

, et al. (2014) Vibration of vehicle-pavement coupled system based on a Timoshenko beam on a nonlinear foundation. Journal of Sound and Vibration 333(24): 6623–6636.

(2015) The effect of switching cracks on the vibration of a continuous beam bridge subjected to moving vehicles. Journal of Sound and Vibration 339: 157–175.

Jian

Xia

Sun

(2022) Indirect identification of bridge frequencies using a four-wheel vehicle: theory and three-dimensional simulation. Mechanical Systems and Signal Processing 177: 109155.

Koh

Ong

JSY

Chua

DKH

, et al. (2003) Moving element method for train-track dynamics. International Journal for Numerical Methods in Engineering 56(11): 1549–1567.

Krishnanunni

Rao

(2019) Decoupled technique for dynamic response of vehicle-pavement systems. Engineering Structures 191: 264–279.

Lan

Lin

(2023) Indirect damage detection for bridges using sensing and temporarily parked vehicles. Engineering Structures 291: 116459.

10.

Liu

Tao

Zhao

, et al. (2024) Dynamic response analysis for bridges subjected to moving vehicle loads by using the analytical dynamic stiffness method. Computers & Structures 292: 107240.

11.

Tian

Liu

(2010) Dynamic response solution in transient state of viscoelastic road under moving load and its application. Journal of Engineering Mechanics 136(168): 168–173.

12.

Quan

Dong

, et al. (2022) Dynamic response analysis of vehicle and asphalt pavement coupled system with the excitation of road surface unevenness. Applied Mathematical Modelling 104: 421–438.

13.

Mazilu

(2007) Green’s functions for analysis of dynamic response of wheel/rail to vertical excitation. Journal of Sound and Vibration 306(1–2): 31–58.

14.

Mazilu

(2010a) Interaction between a moving two-mass oscillator and an infinite homogeneous structure: green’s functions method. Archive of Applied Mechanics 80(8): 909–927.

15.

Mazilu

(2010b) Prediction of the interaction between a simple moving vehicle and an infinite periodically supported rail - green’s functions approach. Vehicle System Dynamics 48(9): 1021–1042.

16.

Mazilu

Dumitriu

Tudorache

(2011) On the dynamics of interaction between a moving mass and an infinite one-dimensional elastic structure at the stability limit. Journal of Sound and Vibration 330(15): 3729–3743.

17.

Metrikine

Popp

(1999) Instability of vibrations of an oscillator moving along a beam on an elastic half-space. European Journal of Mechanics - A: Solids 18(2): 331–349.

18.

Pillai

Talukdar

(2023) Fatigue life assessment of a plate girder bridge using an uncoupled iterative scheme for bridge–vehicle interaction. Journal of Bridge Engineering 28(2): 1–14.

19.

Qin

Wang

Xiang

, et al. (2019) Speed independent road classification strategy based on vehicle response: theory and experimental validation. Mechanical Systems and Signal Processing 117: 653–666.

20.

Rózanski

Sikora

Smuda

, et al. (2021) On theoretical and practical aspects of Duhamel’s integral. Archives of Control Sciences 31(4): 815–847.

21.

Sheng

Jones

CJC

Thompson

(2004) A theoretical model for ground vibration from trains generated by vertical track irregularities. Journal of Sound and Vibration 272(3–5): 937–965.

22.

Sheng

Xiao

Zhang

(2016) The time domain moving Green function of a railway track and its application to wheel–rail interactions. Journal of Sound and Vibration 377: 133–154.

23.

Sheng

Cheng

Thompson

(2020) Modelling wheel/rail rolling noise for a high-speed train running along an infinitely long periodic slab track. Journal of the Acoustical Society of America 148(1): 174–190.

24.

Sun

Chen

(2021) Modelling of periodic slab track using time-frequency hybrid Green’s function method and its application to vehicle-track dynamic interaction. Journal of Sound and Vibration 511: 116327.

25.

Wang

Zhang

Ouyang

(2020) Influence of random multi-point seismic excitations on the safety performance of a train running on a long-span bridge. International Journal of Structural Stability and Dynamics 20(4): 2050054.

26.

Wang

Zhang

Ouyang

(2021) Modeling uncertainties of vehicle-track coupled dynamic systems. Mechanics Based Design of Structures and Machines 49(7): 947–968.

27.

Yang

Zhu

, et al. (2023) A reduced-plate model transmission method for fast dynamic analysis of vehicle–pavement interaction. Journal of Sound and Vibration 548: 117554.

28.

Yang

(2010) Investigation on dynamical interaction between a heavy vehicle and road pavement. Vehicle System Dynamics 48(8): 923–944.

29.

Yang

Ding

Chen

(2013) Dynamic response to a moving load of a Timoshenko beam resting on a nonlinear viscoelastic foundation. Acta Mechanica Sinica/Lixue Xuebao 29(5): 718–727.

30.

Zhai

Xia

Cai

, et al. (2013) High-speed train–track–bridge dynamic interactions–Part I: theoretical model and numerical simulation. International Journal of Reality Therapy 1(1): 3–24.

31.

Zhang

Cai

(2012) Fatigue reliability assessment for existing bridges considering vehicle speed and road surface conditions. Journal of Bridge Engineering 17(3): 443–453.

32.

Zhang

Cheng

Sheng

, et al. (2020) Dynamic wheel-rail interaction at high speed based on time-domain moving Green’s functions. Journal of Sound and Vibration 488: 115632.

33.

Zhang

Hou

Duan

, et al. (2021) Road roughness estimation based on the vehicle frequency response function. Actuators 10(5): 1–20.

34.

Zhang

Yang

, et al. (2022) Study on crack propagation path of asphalt pavement under vehicle-road coupled vibration. Applied Mathematical Modelling 101: 481–502.

35.

Zhang

Hou

, et al. (2023) Vehicle parameter identification based on vehicle frequency response function. Journal of Sound and Vibration 542: 117375.

36.

Zhao

Wang

(2020) Dynamic pavement response analysis under moving truck loads with random amplitudes. Journal of Transportation Engineering, Part B: Pavements 146(2): 04020020.

37.

Zhu

McCrum

Keenahan

(2024) Full-scale aerodynamic study on the effects of tower wind shields and road-sign gantries on passing high-sided vehicles in the tower region at the Queensferry Crossing. Advances in Structural Engineering 27(6): 1055–1070.