System and source identification from operational vehicle responses: A novel modal model accounting for the track-vehicle interaction

TL;DR

This paper introduces a novel modal model for operational vehicle response analysis that accounts for track-vehicle interaction, enabling more accurate identification of vehicle dynamics and surface roughness during operation.

Contribution

The paper develops a specialized Track-Vehicle Interaction Modal Model that explicitly incorporates road/rail inputs, improving operational modal analysis accuracy for vehicles.

Findings

The new model accurately captures vehicle responses under realistic road conditions.

It enables cost-effective, simple tools for vehicle operational identification.

The approach allows indirect characterization of surface roughness.

Abstract

Operational Modal Analysis (OMA) is a powerful tool, widely used in the fields of structural identification and health monitoring, and certainly eligible for identifying the real in-operation behaviour of vehicle systems. Several attempts can be found in the literature, for which the usage of algorithms based on the classical OMA formulation has been strained for the identification of passenger cars and industrial trucks. The interest is mainly focused on the assessment of suspension behaviour and, thus, on the identification of the so-called vehicle rigid body modes. But issues arise when the operational identification of a vehicle system is performed, basically related to the nature of the loads induced by the roughness of rolling profiles. The forces exerted on the wheels, in fact, depending on their location, are affected by time and/or spatial correlation, and, more over, do not fit the form of white noise sequences. Thus, the nature of the excitation strongly violate the hypotheses on which the formulation of classical OMA modal model relies, leading to pronounced modelling errors and, in turn, to poorly estimated modal parameters. In this paper, we develop a specialised modal model, that we refer to as the Track-Vehicle Interaction Modal Model, able to incorporate the character of road/rail inputs acting on vehicles during operation. Since in this novel modal model the relationship between vehicle system outputs and modal parameters is given explicitly, the development of new specific curve fitting techniques, in the time-lag or frequency domain, is now possible, making available simple and cost-effective tools for vehicle operational identification. More over, a second, but not less important outcome of the proposed modal model is the usage of the resulting techniques for the indirect characterisation of rolling surface roughness, that can be used to improve comfort and safety.

Tables

Table 1. Table 1: Parameters of the full-car model [ 3 ] .

Quantity	Unit	Value	Quantity	Unit	Value
$W_1$	m	1.490	$c_r$	kN s/m	1.310
$L_{1f}$	m	1.064	$m_{u1}$	kg	57.500
$L_{1r}$	m	1.596	$m_{u2}$	kg	57.500
$m_s$	kg	1150	$m_{u3}$	kg	57.500
$j_{sx}$	kg m2	530	$m_{u4}$	kg	57.500
$j_{sy}$	kg m2	1630	$k_{ft}$	kN/m	140
$k_f$	kN/m	15.750	$k_{rt}$	kN/m	140
$k_r$	kN/m	14	$c_{ft}$	kN s/m	0.150
$c_f$	kN s/m	1.475	$c_{rt}$	kN s/m	0.150

Equations

$${\lambda }_n,{\lambda }_n^{\ast }=-{\zeta }_n{\omega }_n\pm \mathrm{i}{\omega }_n\sqrt{1-{\zeta }_n^2},$$

$${{\boldsymbol{S}}_{f}(\omega)}=\left[{\begin{array}[]{*{20}{c}}{{S_{1}}}&{}\hfil&0\\ {}\hfil&\ddots&{}\hfil\\ 0&{}\hfil&{{S_{{N}}}}\end{array}}\right].$$

$${\mathbf{R}}_q\left(\tau \right)=\sum _{n=1}^{2N}{\mathbf{\varphi }}_n{\mathbf{\psi }}_n^T{e}^{+{\lambda }_n\tau }h\left(\tau \right)+{\mathbf{\psi }}_n{\mathbf{\varphi }}_n^T{e}^{-{\lambda }_n\tau }h\left(-\tau \right),$$

$${\mathbf{S}}_q(\omega )=\sum _{n=1}^{2N}\frac{{\mathbf{\varphi }}_n{\mathbf{\psi }}_n^T}{\text{i}\omega -{\lambda }_n}+\frac{{\mathbf{\psi }}_n{\mathbf{\varphi }}_n^T}{-\text{i}\omega -{\lambda }_n},$$

$${\Gamma }_p\left(\nu \right)=\frac{\left|S_{d_i{d}_j}\left(\nu \right)\right|}{\sqrt{S_{d_i}\left(\nu \right)S_{d_j}\left(\nu \right)}},$$

$$S_{d_i{d}_j}\left(\nu \right)={\Gamma }_p\left(\nu \right)\sqrt{S_{d_i}\left(\nu \right)S_{d_j}\left(\nu \right)}{\mathrm{e}}^{-\mathrm{i}{\beta }_p\left(\nu \right)},$$

$$S_{d_i}\left(\nu \right)=S_{d_j}\left(\nu \right)=S_d\left(\nu \right);$$

$$S_{d_i{d}_j}\left(\nu \right)={\Gamma }_p\left(\nu \right)S_d\left(\nu \right),$$

$${\mathbf{S}}_r\left(\omega \right)$$

$$\displaystyle=\left[{\begin{array}[]{*{20}{c}}1&{{\it{\Gamma}_{1}}\left(\omega\right){{\rm{e}}^{+{\rm{i}}\omega{\tau_{1}}}}}&{{\it{\Gamma}_{1}}\left(\omega\right){{\rm{e}}^{-{\rm{i}}\omega{\tau_{2}}}}}&{{\it{\Gamma}_{2}}\left(\omega\right)}\\ {{\it{\Gamma}_{1}}\left(\omega\right){{\rm{e}}^{-{\rm{i}}\omega{\tau_{1}}}}}&1&{{{\rm{e}}^{-{\rm{i}}\omega{\tau_{3}}}}}&{{\it{\Gamma}_{3}}\left(\omega\right){{\rm{e}}^{-{\rm{i}}\omega{\tau_{1}}}}}\\ \vdots&{}\hfil&1&{{\it{\Gamma}_{3}}\left(\omega\right){{\rm{e}}^{{\rm{+i}}\omega{\tau_{2}}}}}\\ {{\it{\Gamma}_{2}}\left(\omega\right)}&\cdots&\cdots&1\end{array}}\right]{S_{d}}\left(\omega\right),$$

$${{\boldsymbol{G}}_{fr}}=\left[{\begin{array}[]{*{20}{c}}{{\textbf{{0}}^{\left({N-{N_{t}}}\right)\times{N_{t}}}}}&{{\textbf{{0}}^{\left({N-{N_{t}}}\right)\times{N_{t}}}}}\\ {{{\boldsymbol{K}}_{t}}}&{{{\boldsymbol{C}}_{t}}}\end{array}}\right],$$

$${{\boldsymbol{S}}_{f}}\left(\omega\right)={\boldsymbol{G}_{fr}}\left({{\boldsymbol{{D}}}\left(\omega\right)\otimes{{\boldsymbol{S}}_{r}}\left(\omega\right)}\right)\boldsymbol{G}_{fr}^{T}\quad\mbox{with}\quad{\boldsymbol{D}}\left(\omega\right)=\left[{\begin{array}[]{*{20}{c}}1&{{\rm{i}}\omega}\\ {-{\rm{i}}\omega}&{{\omega^{2}}}\end{array}}\right],$$

$$\mathbf{M}\ddot{\mathbf{q}}(t)+\mathbf{C}\dot{\mathbf{q}}(t)+\mathbf{K}\mathbf{q}(t)=\mathbf{f}(t),$$

$$\mathbf{M}\dot{\mathbf{q}}(t)-\mathbf{M}\dot{\mathbf{q}}(t)={\text{mbox}0}^{N\times N},$$

$${\boldsymbol{{x}}(t)}=\left[{\begin{array}[]{*{20}{c}}{\boldsymbol{q}(t)}\\ {{\dot{\boldsymbol{q}}}(t)}\end{array}}\right]\qquad{{\boldsymbol{\dot{x}}}}(t)=\left[{\begin{array}[]{*{20}{c}}{{\dot{\boldsymbol{q}}}(t)}\\ {{\ddot{\boldsymbol{q}}}(t)}\end{array}}\right]\qquad{\boldsymbol{{u}}(t)}=\left[{\begin{array}[]{*{20}{c}}{\boldsymbol{f}(t)}\\ {\mbox{{{0}}}}^{N\times 1}\end{array}}\right],$$

$$\mathbf{P}\dot{\mathbf{x}}(t)+\mathbf{Q}\mathbf{x}(t)=\mathbf{u}(t),$$

$${\boldsymbol{P}}=\left[{\begin{array}[]{*{20}{c}}{\boldsymbol{C}}&{\boldsymbol{M}}\\ {\boldsymbol{M}}&{\textbf{{0}}^{N\times N}}\end{array}}\right]\qquad{\boldsymbol{Q}}=\left[{\begin{array}[]{*{20}{c}}{\boldsymbol{K}}&{\textbf{{0}}^{N\times N}}\\ {\textbf{{0}}^{N\times N}}&{-{\boldsymbol{M}}}\end{array}}\right].$$

$$\mathbf{A}=\mathbf{V}\Lambda {\mathbf{V}}^{-1},$$

$$\Lambda =\mathrm{diag}\left({\lambda }_1,\dots ,{\lambda }_N,{{\lambda }_1}^{\ast },\dots ,{{\lambda }_N}^{\ast }\right),$$

$$\boldsymbol{V}=\left[{\begin{array}[]{*{20}{c}}{{{\boldsymbol{\psi}}_{1}}}&\cdots&{{{\boldsymbol{\psi}}_{N}}}&{{\boldsymbol{\psi}}_{1}^{*}}&\cdots&{{\boldsymbol{\psi}}_{N}^{*}}\\ {{\lambda_{1}}{{\boldsymbol{\psi}}_{1}}}&\cdots&{{\lambda_{N}}{{\boldsymbol{\psi}}_{N}}}&{\lambda_{1}^{*}{\boldsymbol{\psi}}_{1}^{*}}&\cdots&{\lambda_{N}^{*}{\boldsymbol{\psi}}_{N}^{*}}\end{array}}\right].$$

$${\boldsymbol{x}(t)}=\boldsymbol{V}{\boldsymbol{p}(t)}\quad\Rightarrow\quad\left[{\begin{array}[]{*{20}{c}}{\boldsymbol{q}(t)}\\ {{\dot{\boldsymbol{q}}}(t)}\end{array}}\right]=\sum\limits_{n=1}^{2N}{\left[{\begin{array}[]{*{20}{c}}{{{\boldsymbol{\psi}}_{n}}}\\ {{\lambda_{n}}{{\boldsymbol{\psi}}_{n}}}\end{array}}\right]}{p_{n}}\left(t\right),$$

$${\mathbf{M}}_a\dot{\mathbf{p}}(t)+{\mathbf{M}}_b\mathbf{p}(t)={\mathbf{V}}^T\mathbf{u}(t),$$

$${\mathbf{M}}_a={\mathbf{V}}^T\mathbf{P}\mathbf{V}=\text{diag}\left(m_{a1},\dots ,m_{aN},m_{a1}^{\ast },\dots ,m_{aN}^{\ast }\right)$$

$${\mathbf{M}}_b={\mathbf{V}}^T\mathbf{Q}\mathbf{V}=\text{diag}\left(m_{b1},\dots ,m_{bN},m_{b1}^{\ast },\dots ,m_{bN}^{\ast }\right),$$

$${{\boldsymbol{V}}^{T}}{\boldsymbol{u}}\left(t\right)={\left[{\begin{array}[]{*{20}{c}}{{{\boldsymbol{\psi}}_{1}}}&\cdots&{{{\boldsymbol{\psi}}_{N}}}&{{\boldsymbol{\psi}}_{1}^{*}}&\cdots&{{\boldsymbol{\psi}}_{N}^{*}}\end{array}}\right]^{T}}{\boldsymbol{f}}\left(t\right).$$

$$m_{an}{\dot{p}}_n\left(t\right)+m_{bn}{p}_n\left(t\right)={{\mathbf{\psi }}_n}^T\mathbf{f}(t)\text{mbox}with{m}_{bn}=-{\lambda }_n{m}_{an},$$

$$p_n\left(t\right)=\frac{{{\mathbf{\psi }}_n}^T}{m_{an}}\underset{-\infty }{\overset{t}{\int }}\mathbf{f}(ϵ){\mathrm{e}}^{{\lambda }_n\left(t-\epsilon \right)}\mathrm{d}\epsilon .$$

$$\mathbf{q}(t)=\sum _{n=1}^{2N}{\mathbf{\psi }}_n{p}_n\left(t\right)=\sum _{n=1}^{2N}\frac{{\mathbf{\psi }}_n{{\mathbf{\psi }}_n}^T}{m_{an}}\underset{-\infty }{\overset{t}{\int }}\mathbf{f}(\rho ){\mathrm{e}}^{{\lambda }_n\left(t-\rho \right)}\mathrm{d}\rho .$$

$$\mathbf{g}\left(t\right)=\sum _{n=1}^{2N}\frac{{\mathbf{\psi }}_n{{\mathbf{\psi }}_n}^T}{m_{an}}{\mathrm{e}}^{{\lambda }_{n}t}.$$

$$q_i\left(t\right)=\sum _{n=1}^{2N}\frac{{{\psi }_{in}{\mathbf{\psi }}_n}^T}{m_{an}}\underset{-\infty }{\overset{t}{\int }}\mathbf{f}\left(\rho \right){\mathrm{e}}^{{\lambda }_n\left(t-\rho \right)}\mathrm{d}\rho ,$$

$$q_j\left(t+\tau \right)=\sum _{m=1}^{2N}\frac{{{\psi }_{jm}{\mathbf{\psi }}_m}^T}{m_{am}}\underset{-\infty }{\overset{t+\tau }{\int }}\mathbf{f}\left(\sigma \right){\mathrm{e}}^{{\lambda }_m\left(t+\tau -\sigma \right)}\mathrm{d}\sigma ,$$

Full text

\floatsetup

[figure]style=plain,subcapbesideposition=top

1 Introduction

The identification of the real in-operation behaviour is of utmost importance in all the stages of the vehicle design process, from model updating, extensively used as a tool for improving the accuracy of dynamics simulations, to optimisation strategies and techniques, product life-cycle management implementations, and the verification of controllers’ performance [1, 2, 3]. In this context, a wheeled vehicle, cruising at constant speed, on a certain road profile, is a common working occurrence for designers, in which the vehicle is mainly subjected to the external loads exerted by the interaction with the rolling rough surface. From an experimental point of view, since unknown steady-state inputs are applied to the system, Operational Modal Analysis (OMA), even referred to as in-Operation or Output-only Modal Analysis, is certainly eligible to perform the system structural identification, moving from output data only [4, 5]. Specifically, in the case of linear systems subjected to operational loads, OMA allows for estimating the so-called modal parameters, that is the resonance frequencies, the damping ratios, the mode shapes, and the operational reference vectors, that can be utilised, in turn, to obtain a mathematical model of the system under test, generally referred to as the modal model [6].

Ordinary vehicle systems generally exhibit a very complex dynamics, mainly owing to nonlinearities related to shock absorbers’ operation and to the kinematics imposed by the suspension design [7, 8]. Thus, modal models synthesised under the real operational loadings actually represent an equivalent linearisation of the tested nonlinear systems around interesting and representative working points [9, 10]. Different speeds, manoeuvres, as cruising or constant radius cornering, the presence of payloads, and of different adopted designs of suspension result in different equivalent linearised models.

A novel, but even challenging application of OMA is being proposed in the framework of vehicle ride dynamics [11, 12]. Since modal parameter estimation allows for quantifying the overall damping associated to the vehicle rigid-body modes, OMA of vehicle systems can be employed as a tool for assessing the performance of different suspension systems. More over, in the case of semi-active and active suspensions, the application of OMA is extremely interesting for the validation and optimisation of controllers’ performance [12].

During the last decade, several robust methods have been developed for estimating reliable modal parameters from output-only data [13]. Basically, these methods rely on the Natural Excitation Technique (NExT) assumptions [14, 15, 16, 17]: (i) The unknown loads acting on the system have to fit the form of white noise sequences, (ii) in case of multi-point excitation, the external inputs are required to be strictly uncorrelated. The NExT assumptions are commonly satisfied in civil engineering applications, as is the case of high rise buildings and suspension bridges, where typically the structures are excited by environmental loadings due to traffic, wind and waves [18].

In this paper, we first focus on the issues that arise when an operational structural identification of a vehicle system is performed by using classical OMA based techniques. These issues are related to several aspects: (i) The presence of high modal density and closely-spaced modes; (ii) the high amount of damping owing to the presence of shock-absorbers; (iii) the nature of the loads induced by the rolling surface. Specifically, we show that the latter point implies a strong violation of the NExT assumptions. In fact, the trend of each excitation, transmitted to the vehicle through the wheel stiffness and damping, is quite far from the flat frequency shaping of white noise [19, 20]. In addition, depending on their location, the inputs acting on the wheels result affected by time and/or spatial correlation [21, 22, 23, 24, 25, 26]. As a consequence, the application of classical OMA methodologies for post-processing the vehicle output responses may lead to pronounced modelling errors and, therefore, to poorly estimated modal parameters.

Thus, by assuming that the rolling surface is an homogeneous Gaussian random field [27, 28], we formulate of a novel OMA modal model, referred to as the Track-Vehicle Interaction Modal Model (TVIMM), specialised and suitable for developing new identification procedures aimed at the estimation of the modal parameters of a vehicle system. More over, we show that the coefficients of the adopted empiric surface model can be even estimated as a further outcome of our OMA processing. Since direct measurements performed by using specific profilometers can result to be expensive, this second application of the model has relevant interest in road and railway health monitoring. In this field, in fact, the indirect characterisation of surface roughness can be employed as a tool for improving safety and comfort [29].

The rest of the paper is organised as follows. In Section 2, we obtain the representation of surface-induced forces exerted on a vehicle and discuss the issues arising when identification procedures based on the classical OMA formulation are used to process the output responses of vehicle systems. In Section 3, we propose a novel OMA formulation, specifically designed for the operational identification of vehicle systems, based on a specialised modal model accounting for the track-vehicle interaction. In Section 4, we offer a numerical demonstration of the developed formulation. Concluding remarks are summarised in Section 5. In Appendix A, we provide mathematical proof of some fundamental equations used in Section 3.

2 System and source identification from operational vehicle responses

The nature of loadings exerted on a vehicle by the rolling surface plays a fundamental role in the process of operational identification of the system. We here exploit the properties of homogeneous Gaussian random fields to derive a model of the surface-induced forces. By comparing the resulting surface inputs with those permitted by complying the NExT assumptions, we conclude that a specialised OMA formulation is needed, able to provide a correct modelling of the system and of the main excitation source, that is the surface roughness.

2.1 Classical OMA approach

By processing the output responses of a $N$ degrees of freedom (dofs) system, OMA leads to the estimation of the modal parameters, that is the poles $\lambda_{n}$, the modal vectors {{\boldsymbol{\psi}}_{n}}={\left[{\begin{array}[]{*{20}{c}}{{\psi_{1n}}}&\cdots&{{\psi_{Nn}}}\end{array}}\right]^{T}}\in{\mathbb{C}^{N\times 1}} and the operational reference vectors {{\boldsymbol{\phi}}_{n}}={\left[{\begin{array}[]{*{20}{c}}{{\phi_{1n}}}&\cdots&{{\phi_{Nn}}}\end{array}}\right]^{T}}\in{\mathbb{C}^{N\times 1}} (with $n=1,\ldots,N$), where the symbol $(.)^{T}$ indicates matrix transposition. The computation of system’s poles $\lambda_{n}$ is of fundamental interest as they contain information about the resonance frequencies $f_{n}=\omega_{n}/2\pi$ and the damping factors $\zeta_{n}$

[TABLE]

where $(.)^{*}$ indicates complex-conjugation. The application of methods of system identification to OMA has given rise to a large variety of estimation techniques [6, 30], comprising, generally, several steps. The most advanced procedures are based on the usage of the modal model, commonly utilised to complete the identification step and for validating the results of the estimation process [17]. Recently, a new approach has been proposed, allowing for the computation of modal parameters directly from the modal model [31].

In the case of classical OMA, the modal model is formulated both in the time-lag $\tau$ and in the angular frequency $\omega=2\pi f$ domains, by referring to correlation functions and power spectral densities (PSDs), respectively. Thus, by considering a system of external forces relying on the NExT assumptions, the input PSD matrix ${{\boldsymbol{S}}_{f}(\omega)}\in{\mathbb{C}^{{N}\times{N}}}$ can be written as

[TABLE]

In fact, since the inputs are required to be strictly uncorrelated white noise sequences, the matrix in Eq. (2) has only constant entries along the main diagonal. With referring to a $N$ dofs system subjected to the operational loadings Eq. (2), the formulation of the modal model can be given in terms of output correlation matrix $\boldsymbol{R}_{q}(\tau)\in{\mathbb{R}^{{N}\times{N}}}$ and output PSD matrix $\boldsymbol{S}_{q}(\omega)\in{\mathbb{C}^{{N}\times{N}}}$, as

[TABLE]

[TABLE]

where $h(.)$ indicates the Heaviside step function.

2.2 Random fields for surface profiles modelling

Generally, in case of OMA, input loads are uncontrollable and, in addition, remain unmeasured. For this reason, the analyst has to verify that the operational loads are suitable to adequately and effectively excite the system in the frequency range of interest.

A vehicle system in steady-state working conditions is subjected to both external and on-board excitations [1]. On-board forces are related to the operation of rotating parts and are specifically imputable to engine operation and to non-uniformities in assemblies and components of the driveline. These forces introduce harmonics in the spectra, whose frequencies are expected to depend on the engine rotating velocity and, in turn, on the vehicle speed. External loads are mainly related to surface roughness, that represent the only present source able to generate an effective broad band excitation. Specifically, the cut-off frequency of the input spectrum depends on roughness spatial frequency content, vehicle velocity and wheel dynamic behaviour.

Random fields for surface modelling have extensively been studied and the usage of different approaches has emerged [27, 28]. Basically, the following assumptions are commonly adopted: (i) the surface roughness is an homogeneous random field; (ii) the height of the asperities satisfies to a zero-mean Gaussian distribution; (iii) the pavement unevenness is an ergodic random process. Homogeneity implies that the statistical properties of surface roughness are independent on spatial observations. This assumption has implications similar to stationarity for one-dimensional random processes. The second assumption ensures that the output responses of a linear system subjected to pavement excitations satisfy, as well, to a zero-mean Gaussian distribution. Ergodicity guarantees that average is equal to expectation calculated over the whole ensemble.

In simple track-vehicle interaction models, the unilateral contact point hypothesis is often utilised to describe the forces transmitted to the vehicle through the wheel stiffness and damping. This assumption is made without loss of generality, by considering that the distributed contact in the wheel-pavement interface acts as a low-pass filter, whose bandwidth is governed by the contact interface itself [32]. In this case, the surface roughness can be supposed to be a single-track random process, describing the longitudinal profile along the wheel path in the travelling direction [19]. Generally, an empiric parametric model is adopted to fit the measured PSD, and the captured surface roughness is classified or employed in the simulations. The non-smoothed PSD is often approximated by a simple function, involving the usage of only few parameters. A literature survey on existing different approximations for longitudinal road profiles is presented in Ref. [20].

When more than one profile is needed, the ordinary coherence function is introduced to express the relationship between multiple tracks [33]. Considering two parallel surface profiles $d_{i}(x)$ and $d_{j}(x)$, with $x$ the longitudinal spatial variable, separated by a distance $W_{p}$, the ordinary coherence function is defined in the angular spatial frequency domain $\nu$ as

[TABLE]

where $S_{d_{i}}(\nu)$ and $S_{d_{j}}(\nu)$ are the auto-PSDs of $d_{i}(x)$ and $d_{j}(x)$, respectively, and $S_{{d_{i}}{d_{j}}}(\nu)$ is the cross-PSD between the two profiles. By definition, the coherence function is a real even function ranging from 0 to 1. As a consequence, two separated surface profiles perfectly overlap at wavelengths corresponding to amplitude values of the coherence function equal to 1. The expression of cross-PSD is obtained from Eq. (5) as

[TABLE]

where $\beta_{p}(\nu)$ is the difference between the phases related to $d_{i}(x)$ and $d_{j}(x)$.

Based on the homogeneity assumption, the following properties can be exploited for simplifying the description of multiple tracks [27, 28]:

(i) The auto-PSDs related to parallel surface profiles are coinciding and, thus, equal to the same function

[TABLE]

(ii) the cross-PSD between two parallel surface profiles is a real and even function depending on auto-PSD and coherence function

[TABLE]

in which it is specifically useful to notice the cancellation of the phase difference $\beta_{p}(\nu)$ in Eq. (8).

2.3 Problem statement

We here consider the generic $N_{t}$-wheel vehicle system represented in Fig. 1, with $N_{t}=4$, even if the following considerations can be easily extended to systems equipped with $N_{t}\neq 4$ wheels. The considered geometry includes $N_{w}=3$ different trackwidths $W_{p}$ (with $p=1,\ldots,N_{w}$) and $N_{l}=3$ different wheelbases $L_{l}$ (with $l=1,\ldots,N_{l}$). We hypothesise that the vehicle travels at constant velocity $V$ on an homogeneous Gaussian surface. By writing Eqs. from (5) to (8) in the angular frequency domain $\omega=\nu V$, we obtain the following PSD matrix of surface-induced displacements $\boldsymbol{S}_{r}(\omega)\in{\mathbb{C}^{{N_{t}}\times{N_{t}}}}$

[TABLE]

where $\tau_{l}=L_{l}/V$ indicate the time-delays between the inputs acting on different axles.

To describe the forces transmitted to the vehicle through the wheel stiffness and damping, we first introduce the following static gain matrix $\boldsymbol{G}_{fr}\in{\mathbb{R}^{{N}\times{2N_{t}}}}$

[TABLE]

where $\boldsymbol{K}_{t}$ and $\boldsymbol{C}_{t}\in{\mathbb{R}^{{N_{t}}\times{N_{t}}}}$ are diagonal matrices in which stiffness and damping terms are, respectively, collected. We, second, account for the PSDs of the derived processes, in addition to Eq. (13), in case tyre damping is included in the wheel model. We, finally, obtain the following representation for the input PSD matrix related to surface-induced forces ${{\boldsymbol{S}}_{f}(\omega)}\in{\mathbb{C}^{{N}\times{N}}}$

[TABLE]

where the matrix $\boldsymbol{D}(\omega)\in{\mathbb{C}^{{4}\times{4}}}$ allows for taking into account the contributions related to road-induced velocities and the symbol $\otimes$ denotes the Kronecker product. We stress that the sign of the imaginary part of Eq. (20) depends on the adopted definition of correlation function, where a change of definition in the time-lag domain leads to a complex-conjugate expression in the frequency domain.

By comparing Eqs. (2) and (20), we conclude that for a vehicle system subjected to operational loadings, the modal model resulting from the classical OMA approach (Eqs. (3) and (4)) is no more valid. The arisen issues are mainly related to the nature of the loads induced by the rolling surface, that violates the NExT assumptions. First, we notice that forces applied to each wheel are basically coloured excitations, with magnitude inversely proportional to the frequency raised to a certain power, and shaping depending on the auto-PSD of the road profile ($S_{d}(\omega)$ Eq. (7)) and on the wheel parameters ($\boldsymbol{G}_{fr}$ Eq. (19)). In addition, being systems equipped with more than one wheel, forces are affected by time and/or spatial correlation: (i) wheels mounted on the same axle are subjected to spatially correlated inputs ($\it{\Gamma}_{p}(\omega)$ in Eq. (13)); (ii) wheels travelling on the same path and located on different axles are subjected to time correlated inputs (${{\rm{e}}^{\pm{\rm{i}}\omega{\tau_{l}}}}$ in Eq. (13)); (iii) wheels travelling on separated paths and located on different axles are subjected to time and spatially correlated inputs ($\it{\Gamma}_{p}\left(\omega\right){{\rm{e}}^{\pm{\rm{i}}\omega{\tau_{l}}}}$ in Eq. (13)).

Owing to the aforementioned effects, the application of classical OMA methodologies to vehicle responses may lead to pronounced modelling errors and, therefore, to poorly estimated modal parameters. The formulation of a specialised modal model for correctly describing the track-vehicle interaction is required and, in turn, a specific OMA formulation is needed to understand how Eqs. (3) and (4) are modified. We stress that a modal model providing the relationship between the generic output of the system and its modal parameters in an explicit form is needed for the formulation of a whatever procedure for modal parameter estimation based on suitable curve fitting algorithms. We more over comment that since operational vehicle responses incorporate information about the surface roughness, the identification procedures based on this novel OMA formulation would allow for estimating, in addition to modal parameters, even the coefficients of the adopted empiric surface models.

3 The Track-Vehicle Interaction Modal Model

We here utilise a 7 dofs system ($N=7$), generally referred to as the full-car model, to introduce the theoretical background of the proposed Track-Vehicle Interaction Modal Model. To obtain the analytical expression of the generic system output, we solve a Duhamel integral in modal coordinates. Moving from Eq. (20), we compute the correlation matrix of surface-induced forces and achieve the TVIMM formulation by using properties of convolution integrals and Fourier transform.

3.1 Full-car model

The full-car model (Fig. 2) is the simplest mathematical description of a four-wheel vehicle ($N_{t}=4$), whose predicted output responses incorporate all the effects of the issues discussed in Section 2, which make no longer possible the usage of classical OMA methodologies for post-processing. Since full-car model offers a good trade-off between model complexity and accuracy, this linear lumped-parameters system is commonly utilised for simulating the ride dynamics of passenger and race cars [2, 3]. The model preserves the multi-input nature of road excitation and, different from half-car model, allows for evaluating the contributions to vehicle responses due to roll disturbance produced by two parallel tracks. The geometry comprises one single trackwidth $W_{1}$ and one single wheelbase $L_{1}$. In Fig. 2, we indicate the dofs of the model by using the following notation: $z_{s}$, $\varphi_{s}$, and $\theta_{s}$ represent the heave, the roll and the pitch rigid body motions of the sprung mass (the body); $z_{s1}$, $z_{s2}$, $z_{s3}$, $z_{s4}$ denote the vertical displacements at the four corners of the sprung mass (strut mounts to body); $z_{u1}$, $z_{u1}$, $z_{u3}$, $z_{u4}$ are the rattle displacements of the unsprung masses (the wheels). More over, with regards to system parameters, $m_{s}$, $j_{sx}$, $j_{sy}$ represent the mass and the moments of inertia associated with the body; $m_{u1}$, $m_{u2}$, $m_{u3}$, $m_{u4}$ denote the unsprung masses; $k_{f}$, $k_{r}$, $c_{f}$, $c_{r}$ are the stiffness and damping coefficients of the absorbers; $k_{ft}$, $k_{rt}$, $c_{ft}$, $c_{rt}$ indicate the stiffness and damping coefficients of the tyres; $L_{1f}$, $L_{1r}$ are the distances of the front and rear axle, respectively, from the center of gravity of the unsprung mass, the sum of which equals the wheelbase.

3.2 Duhamel integral in modal coordinates

Since vehicles are generally non-proportional damping systems, we adopt a representation in state-space form. Thus, we recast the motion equations of a $N$ dofs system into an equivalent set of 2$N$ first-order differential equations. Specifically, by denoting with $\boldsymbol{q}(t)\in{\mathbb{R}^{N\times 1}}$ the vector of Lagrangian coordinates and with $\boldsymbol{f}(t)\in{\mathbb{R}^{N\times 1}}$ the vector of external loads, we write the well-known set of system dynamics equations as

[TABLE]

where $t$ is the time variable, and $\boldsymbol{M}$, $\boldsymbol{C}$, and $\boldsymbol{K}\in{\mathbb{R}^{N\times N}}$ are the mass, damping, and stiffness matrices, respectively. By adding the following further set of differential equations

[TABLE]

and by making the substitutions

[TABLE]

we obtain the set of system dynamics equations in state-space form

[TABLE]

where $\boldsymbol{P}$ and $\boldsymbol{Q}\in{\mathbb{R}^{2N\times 2N}}$ are partitioned as

[TABLE]

Here, the product $-{\boldsymbol{P}^{-1}}{{\boldsymbol{Q}}}$ leads to the so-called state matrix $\boldsymbol{A}$. The eigenvalue decomposition of the state matrix $\boldsymbol{A}$ is

[TABLE]

where ${\boldsymbol{\it{\Lambda}}}\in\mathbb{C}^{\>2N\times 2N}$ is the eigenvalue matrix, containing the system poles

[TABLE]

with diag$(.)$ denoting a diagonal matrix, while ${\boldsymbol{V}}\in\mathbb{C}^{\>2N\times 2N}$ is the eigenvector matrix, having the following structure

[TABLE]

We decouple the equations of motion by using the following coordinate transformation from the physical to the modal space

[TABLE]

where {{\boldsymbol{p}}}(t)={\left[{\begin{array}[]{*{20}{c}}{{p_{1}}(t)}&\cdots&{{p_{2N}}(t)}\end{array}}\right]^{T}}\in{\mathbb{R}^{N\times 1}} is the modal state vector. By substituting Eq. (29) in Eq. (24) and pre-multiplying both sides by $\boldsymbol{V}^{T}$, we obtain 2$N$ independent differential equations, that we collect in the following compact form

[TABLE]

where $\boldsymbol{M}_{a}$ and $\boldsymbol{M}_{b}\in{\mathbb{C}^{2N\times 2N}}$ are two diagonal matrices generally referred to as modal a and modal b. Specifically,

[TABLE]

and

[TABLE]

while

[TABLE]

By considering the generic system motion equation of set Eq. (30) referred to the $n$-th mode of vibration,

[TABLE]

and by assuming zero initial conditions, we retrieve the solution in the form

[TABLE]

Eq. (35), usually written in physical coordinates, is the well-known Duhamel integral. By combining Eqs. (29) and (35), we obtain the system response in terms of the physical Lagrangian coordinates

[TABLE]

We comment that Eq. (35) can be interpreted as a convolution integral encompassing the impulse response matrix $\boldsymbol{g}(t)\in{\mathbb{R}^{N\times N}}$ of the system under study

[TABLE]

3.3 Input and output correlation functions

By considering the generic outputs $q_{i}(t)$ and $q_{j}(t)$ in Eq. (36), evaluated at the separated time instants $t$ and $t+\tau$, respectively,

[TABLE]

and

[TABLE]

where the following relations hold between the dummy integration variables

[TABLE]

we derive the resulting output cross-correlation function

[TABLE]

where the symbol $\rm{E}[.]$ indicates the expectation computed over the ensemble and ${\boldsymbol{R}}_{f}(\tau)\in{\mathbb{R}^{N\times N}}$ is the input correlation matrix related to the external forces, which, by definition, is the inverse Fourier transform of the input PSD matrix Eq. (20).

By particularising to the case of full-car model (Fig. 2) the stiffness and damping entries of the static gain matrix $\boldsymbol{G}_{fr}$ Eq. (19)

[TABLE]

we decompose the PSD matrix of surface-induced displacements $\boldsymbol{S}_{r}(\omega)$ through the following summation

[TABLE]

where $N_{w}=1$ in the case of full car model and the matrices ${\boldsymbol{\it{\Delta}}}_{0}(\omega)$ and ${\boldsymbol{\it{\Delta}}}_{1}(\omega)\in{\mathbb{C}^{N_{t}\times N_{t}}}$ are defined as

[TABLE]

[TABLE]

By substituting Eq. (43) in Eq. (20), we obtain the Fourier pair ${\boldsymbol{S}}_{f}(\omega)$ and ${\boldsymbol{R}}_{f}(\tau)$ in the case of interest

[TABLE]

and

[TABLE]

where $R_{d}(\tau)$ is the auto-correlation function associated to surface profiles and $H_{p}(\tau)$ (with $p\neq 0$) represents the inverse Fourier transform (IFT) of the ordinary coherence function $\it{\Gamma_{p}}(\omega)$ referred to two parallel tracks. The terms $\boldsymbol{d}(\tau)\in{\mathbb{R}^{{4}\times{4}}}$, ${\boldsymbol{\delta}}_{0}(\tau)$ and ${\boldsymbol{\delta}}_{1}(\tau)\in{\mathbb{R}^{{N_{t}}\times{N_{t}}}}$ correspond to a representation in the time-lag domain of the matrices included in Eq. (46), that is

[TABLE]

[TABLE]

[TABLE]

where $\delta(.)$ indicates the Dirac delta function, and the symbols $\frac{{{\text{d}}(.)}}{{{\text{d}}\tau}}$ and $\frac{{{{\text{d}}^{2}}(.)}}{{{{\text{d}}^{2}}\tau}}$ denote the first and second order derivative operators, respectively. By combining Eqs. (47) and (3.3), we achieve the final expression of the cross-correlation function ${R_{{q_{i}}{q_{j}}}}\left(\tau\right)$

[TABLE]

3.4 Modal model in the time-lag and frequency domains

Moving from Eq. (51), by means of calculations detailed in Appendix A, we obtain the relationship between the operational modal parameters ($\lambda_{n}$, $\psi_{in}$, $\alpha^{p}_{in}$, $\beta^{pl}_{in}$, $\chi^{pl}_{in}$) and the generic output cross-correlation function ${R_{{q_{i}}{q_{j}}}}(\tau)$. We here rephrase the expression of the output correlation matrix ${\boldsymbol{R}}_{q}\left(\tau\right)$, in the following compact matrix notation

[TABLE]

where the terms $\boldsymbol{\bar{\varphi}}_{n}(\tau)$ are the lag-dependent operational reference vectors, defined as

[TABLE]

By Fourier transforming Eq. (52), we obtain the expression of the output PSD matrix ${\boldsymbol{S}}_{q}\left(\omega\right)$

[TABLE]

where the terms $\boldsymbol{\varphi}_{n}(\omega)$ are the frequency-dependent operational reference vectors, defined as

[TABLE]

The modal model that we introduce through Eqs. (52) and (54), namely the TVIMM, incorporates a combination of several contributions: (i) $R_{d}(\tau)$ and $S_{d}(\omega)$ describe the statistical properties of the single surface profiles; (ii) ${{e^{\pm{\lambda_{n}}{\tau_{l}}}}}h\left({\tau\pm{\tau_{l}}}\right)$ and ${{\text{e}}^{\pm{\text{i}}\omega{\tau_{l}}}}$ account for the time correlation between the inputs applied to wheels travelling on the same path and mounted on different axles; (iii) $H_{p}(\tau)$ and ${\mathit{\Gamma}_{p}}(\omega)$ encompass the spatial correlation due to loads acting on wheels travelling on separated paths. Thus, the vehicle output responses contain information on the dynamics of the system and on its interaction with the surface profiles. Not considering this fundamental last contribution from the resulting modal model (Eqs. (52) and (54)) increase the modelling errors of the identification procedure. As a consequence, the unknown parameters result poorly estimated.

We consider useful to stress that although Eqs. (52) and (54) have been obtained by only considering the rigid-body modes of the full-car model, this novel OMA formulation is of general validity and can be utilised for real-world vehicle systems. Specifically, it is able to include and describe, without loss of generality, even the presence of deformable-body modes of vibration.

3.5 Comments on the usage of half spectra

Several robust OMA algorithms in the frequency domain rely on the estimation of the so called half spectra [18], defined as

[TABLE]

The main advantages of this representation consist in: (i) The usage of lower model orders without affecting the quality of the fitting; (ii) the adoption of well-known windowing functions for reducing the effect of leakage and the influence of samples at the higher time-lags, which are the most affected by noise.

In the case of the classical OMA formulation, half spectra $S_{q_{i}q_{j}}^{(+)}(\omega)$ are computed by Fourier transforming the correlation functions corresponding to positive time-lags

[TABLE]

We here prove that in presence of time correlation this property does no longer hold. In fact, the lag-dependent operational reference vectors $\boldsymbol{\bar{\varphi}}_{n}^{T}(-\tau)$ do not vanish at the positive time-lags

[TABLE]

implying that the Fourier transform of Eq. 52 for $\tau>0$ does not lead to the first fraction summation in Eq. 54.

4 Numerical demonstration

In this section, we offer a numerical demonstration of the effectiveness of Eqs. (52) and (54). Specifically, we first compute the responses of the full-car model to surface-induced excitations by using the following implicit analytical expression holding in the frequency domain

[TABLE]

representing the well-known input-output formula [28]. Since the full-car model utilised in this study is a linear mutliple-input mutliple-output system, the output PSD matrix $\boldsymbol{S}_{q}\left(\omega\right)$ can be calculated at each frequency through Eq. (59). We derive the frequency response function matrix $\boldsymbol{G}_{qf}(\omega)\in{\mathbb{C}^{N\times N}}$, referred to the lumped-parameters system of Fig. 2, by means of the following relation

[TABLE]

Then, we compare the obtained curves with those resulting from the TVIMM. We stress that the comparison, here provided for demonstration purposes, is made between quantities computed by following two different routes: The former is based on Eq. (59), in which the relationship between the generic output of the system and its modal parameters in not explicitly given; the latter is based on the proposed OMA for vehicles’ modal model Eqs. (52) and (54), in which the relationship is instead provided in an explicit form. Only this second formulation allows for developing specialised curve fitting techniques, needed, in turn, for modal parameter estimation. We more over comment that, based on the different approximations proposed in the literature, semi-analytical expressions of $S_{d}(\omega)$ and $\it{\Gamma}_{p}(\omega)$ have to be provided, in order to form the input PSD matrix of surface-induced forces $\boldsymbol{S}_{f}(\omega)$ (Eq. (20)).

4.1 Matrix formulation of the simulated vehicle model

For the full-car model, two alternative sets of degrees of freedom can be used to describe the dynamics of the sprung mass: (i) The heave $z_{s}$, the roll $\varphi_{s}$, and the pitch $\theta_{s}$ or (ii) the vertical translations at the right-front $z_{s1}$, left-front $z_{s2}$, and right-rear $z_{s3}$ corners. In what follows, we analyse both the two cases by introducing two different Lagrangian coordinate vectors, respectively

[TABLE]

The vector of surface-induced forces $\boldsymbol{f}(t)$ is

[TABLE]

The mass, damping and stiffness matrices of the considered vehicle model, related to the first set of Lagrangian coordinates ${\boldsymbol{q}_{1}}(t)$, are given by

[TABLE]

[TABLE]

[TABLE]

while those related to the second set ${\boldsymbol{q}_{2}}(t)$, can be expressed by

[TABLE]

where $\boldsymbol{T}\in{\mathbb{R}^{N\times N}}$ is the transformation matrix

[TABLE]

The parameters of the full-car model employed in the simulation [3] are listed in Tab. 1.

4.2 Single-profile and coherence-based models for parallel road paths

We obtain the input PSD matrix of surface-induced forces by imposing that the auto-PSD of single surface profiles $S_{d}(\omega)$ be equal to the ISO 8608 class C approximation. According to ISO 8608 [33], in fact, roughness level is classified from A to H, and the form of the fitted one-sided auto-PSD is given as follows

[TABLE]

where $\nu_{0}=1$ rad/m denotes the reference angular spatial frequency and $S_{0}$ (m3) is the amplitude of the PSD for $\nu=\nu_{0}$. A constant velocity auto-PSD can be obtained by imposing the undulation exponent $e$ equal to 2. When this occurrence is satisfied, the auto-PSD of surface-induced displacements is calculated through a simple integration of a flat-spectrum velocity signal. Rearranging Eq. (68), the expression of the PSD in the angular frequency domain becomes

[TABLE]

where $\omega_{0}=\nu_{0}V$ indicates the reference angular frequency.

Since in the classification of road profiles the angular spatial frequency usually ranges from $\nu_{a}=0.063$ rad/m to $\nu_{b}=17.77$ rad/m [33] (corresponding to wavelengths from 100 m to 0.35 m, respectively), by considering the usual range of speed values from 10 to 30 m/s, we have that the angular frequency varies in the interesting range from $\omega_{a}=1.89$ rad/s (0.30 Hz) to $\omega_{b}=177.7$ rad/s (28.3 Hz). Thus, to account for these restrictions, we introduce lower and upper limits in Eqs. (68) and (69).

The ISO 8608 class C auto-PSD and auto-correlation function are shown in Figs. 3 and 3, respectively. In general, an analytical or numerical expression of the auto-correlation function ${R_{d}}(\tau)$ can be obtained through inverse Fourier transformation of Eq. (69). Specifically, the curve depicted in Fig. 3, representing ${R_{d}}(\tau)$ for the ISO 8608 case, is build up by numerically computing the IFT of ${S_{d}}(\omega)$.

With regards to the ordinary coherence function, we hypothesize that ${\it\Gamma_{p}}(\omega)$ in Eq. (20) have the trend of the fitting model proposed by Bogsjö [24], and illustrated in Fig. 4. Experimental data acquired on very different road tracks have surprisingly shown a good agreement with this one-parameter exponentially decreasing function. At a constant velocity, the original exponential model can be rearranged in the angular frequency domain as

[TABLE]

From practical considerations, it is basically possible to conclude that ${\it\Gamma_{p}}(\omega)$ has to approach 1 for long wavelengths (small $\omega$) and 0 for short wavelengths (large $\omega$). More over, as the two tracks become closer, that is $W_{p}\rightarrow 0$, then ${\it\Gamma_{p}}(\omega)\rightarrow 1$. The curve in Fig. 4, related to $H_{p}(\tau)$, is even build up by numerically computing the IFT of ${\it\Gamma_{p}}(\omega)$.

4.3 Analysis of the results

By solving Eqs. from (25) to (28), we compute resonance frequencies ($f_{n}$), damping ratios ($\zeta_{n}$) and modal vectors ($\boldsymbol{\psi}_{n}$) of the full-car model. The obtained modes of vibration are reported in Fig. 5 for increasing values of the resonance frequency. While the computed modal parameters depend on the physical parameters of the exploited vehicle model, some considerations of general validity can be made, owing to well-known vehicle design procedures and to the fact that suspensions and wheels mounted on different axles have usually similar stiffness and damping values. Specifically, we notice the presence of high modal density and closely-spaced modes, that represent a first practical issue for the identification of vehicle systems. In the frequency band ranging from 0 to 20 Hz two main sets of rigid-body modes do exist. The three modes involving the chassis motion are all approximately located at 1 Hz, while the four modes related to rattle motions of front and rear axles are found in the range from 8 to 15 Hz depending on the vehicle parameter values (8 Hz in the case of values collected in Tab. 1). A second issue is caused by the presence of shock-absorbers with high amount of damping injected in the system, which is also a typical occurrence in vehicle systems.

In Fig. 6, we show some output cross-PSDs of the full-car model. Owing to the high damping ratios, the resonance peaks, related to modes of vibration, are not clearly evidenced in the magnitude plots. More over, time and spatial correlations of surface-induced excitations produce relevant effects in form of distortions of cross-spectra. The time delay between the loads acting on the front and rear axles introduces the most significant of these distortions. In particular, we find (i) humps appearing in the magnitude of the cross-PSDs related to chassis dofs (Figs. 6 and 6) and (ii) the phase plots referred to the cross-PSDs of front and rear axles (Figs. 6 and 6) revealing the typical saw-tooth trend caused by time-delayed signals. These distortions, well-known in the field of vehicle ride dynamics, are imputable to the so-called wheelbase filtering effect [21]. We comment that it should be expected that the presence in the output spectra of these distortion effects could hamper the correct operation of curve fitting methodologies based on the classical OMA modal model, that specifically relies on the NExT hypotheses.

By inspection of Fig. 6, we find that the curves computed by using the proposed TVIMM perfectly overlap those obtained through the usage of the input-output formula. This proves the effectiveness of the proposed Track-Vehicle Interaction Modal Model in the frequency domain (Eq. (54)). More over, some output cross-correlation functions of the full-car model are depicted in Fig. 7 and even in the case of time-lag domain (Eq. (52)), we obtain the expected overlapping of the curves.

5 Conclusions

Aim of the present paper is to introduce a novel OMA modal model, referred to as the Track-Vehicle Interaction Modal Model, needed for developing operational identification techniques, specifically designed for the case of vehicle systems and road/rail surface roughness. The effectiveness of the proposed formulation has been proved by simulating the response of a typical four-wheel vehicle model. Moving from the obtained expressions, specific estimation algorithms operating in the time-lag or in the frequency domain can be developed and implemented for practical purposes, providing new powerful tools to vehicle designers. By means of the resulting techniques, in fact, a map of vehicle system poles corresponding to the various operating conditions, as mainly at the different cruising speeds, can be achieved. The excitations produced during durability tests, in controlled conditions, could be utilised for providing the required operational loadings. Since these tests are usually already performed in scheduled test campaigns foreseen for vehicle design purposes, the usage of estimation techniques based on the modal model proposed in the paper would result simple and cost-effective. To this aim, both equipped laboratories and specific proving grounds could be successfully employed. As a further outcome of the proposed methodology, the resulting techniques could be used for the indirect characterisation of surface roughness, and, in turn, as a tool for improving safety and comfort.

Acknowledgments

The authors Giovanni De Filippis and Davide Palmieri have attended a full Doctoral Program at Politecnico di Bari, with PhD scholarships funded, respectively, through the PON01_02380 and the PON02_00576_3333585, by the Ministry of Education, University, and Research - MIUR, which is gratefully acknowledged.

Appendix A Mathematical derivation of Eqs. from (52) to (55)

Based on Eq. (19), we write the following equality

[TABLE]

where ${\boldsymbol{\psi}}_{n}^{t}\in{\mathbb{C}^{1\times N_{t}}}$ is a vector containing the modal components related to the wheels, that is, for the full-car model,

[TABLE]

By substituting in Eq. (51), we get

[TABLE]

where the terms ${\boldsymbol{J}}_{p}^{a}(\tau)$, ${\boldsymbol{J}}_{p}^{b}(\tau)$ and ${\boldsymbol{J}}_{p}^{c}(\tau)\in{\mathbb{R}^{N_{t}\times N_{t}}}$ are matrices that encompass convolution integrals involving Dirac delta functions, located at different time-lag values (Eqs. (49) and (50)), and their derivatives. We rephrase the convolution integrals referred to ${\boldsymbol{J}}_{p}^{a}(\tau)$ as

[TABLE]

where the functions $E_{1}(\tau)$, $E_{2}(\tau)$ and $E_{3}(\tau)$ represent the explicit solutions of the convolution integrals, that is

[TABLE]

[TABLE]

[TABLE]

For example, Eq. (A14) can be proved by considering the following resolution scheme

[TABLE]

Thus, by introducing the change of variable

[TABLE]

the last integral in Eq. (A) can be solved as

[TABLE]

Moving from Eq. (A), we derive the expressions of the convolution integrals referred to ${\boldsymbol{J}}_{p}^{b}(\tau)$ and ${\boldsymbol{J}}_{p}^{c}(\tau)$, that are

[TABLE]

[TABLE]

where the solutions of these convolution integrals can be retrieved from that of Eq. (A) by using the properties of Dirac delta function derivatives.

Eqs. (A), (A19) and (A20) lead to the following equality

[TABLE]

By substituting this expression in Eq. (51), we observe that the generic cross-correlation function $R_{q_{i}q_{j}}(\tau)$ can be rewritten as

[TABLE]

where the coefficients $\alpha_{in}^{pl}$, $\beta_{in}^{pl}$, $\gamma_{in}^{pl}$, $\bar{\alpha}_{in}^{pl}$, $\bar{\beta}_{in}^{pl}$, $\bar{\gamma}_{in}^{pl}$ play the same role as the operational reference factors $\phi_{in}$ in Eq. (3). The number of coefficients in Eq. (A) can be significantly decreased by recalling that

[TABLE]

Thus, for the full-car model we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, rearranging Eq. (A) in a more suitable form we obtain

[TABLE]

where ${\bar{\varphi}}_{in}\left(\tau\right)$ is a lag-dependent operational reference factor, defined as

[TABLE]