Abstract
This article considers lateral vibration of an automobile in a socalled quasiplanar model where both the loss of contact and road deformation are taken into account. The automobile with dependent suspension is modeled as a vibration system which has two masses and four degrees of freedom. The deformed road is modeled as an elastic beam which has uniform rectangular crosssection, is simply supported at the two ends and lies on the Kelvin's viscoelastic ground. The loss of contact and the change in dimensions of contact areas are considered. The differential equations of motion of the vehicleroad coupled system which contains a partial differential equation are transformed into a set of all ordinary differential equations by applying the BubnovGalerkin’s method. A procedure for numerically solving the transformed differential equations of motion is proposed. Some illustrating results coming from numerical consideration are also presented in the paper.
Highlights
 We proposed a new of vehicle vibration model which called quasiplanar.
 The change in dimension of the contact area between wheel with road surface is considered.
 By applying the BubnovGalerkin's method, we get a set of all ordinary differential equations.
 Some illustrating numerical example are presented to confirm the dynamic response of vehicle with taking and not taking wheel separation.
1. Introduction
Vibration of automobiles while moving on rough roads appears naturally and frequently. Among components of vibration, the ones which cause vertical vibration of points on the automobile are more considerable. When the level of vibration exceeds a definite threshold, one or more wheels may separate from the road surface and the state of contact between the wheels and the road surface is broken. This phenomenon is called as the loss of contact, or wheel separation. The loss of contact reduces controllability of the automobile in both velocity and direction, therefore, the safety of movement. For this reason, it is needed to take the loss of contact (wheel separation) into account while considering vibration of automobiles.
Without wheel separation, the vehicle dynamics attract a lot of attention from the researchers. In [1, 2], the road damage caused by heavy commercial vehicles is investigated using a quartervehicle model with two degree of freedom (DOF) or a half vehicle with four DOF. Using a four DOF model, the dynamics of vehicles with nonlinear spring and damping element is considered by Zhu and Ishitobi in [3]. The researches on vehicletrack and vehiclebridge interaction have developed rapidly, in which the vehicle and the road are linked by the tire forces. In [4], a vibration seven DOF model of heavy vehicleroad coupled system with independent suspensions is presented, where the pavementfoundation is modeled as a doublelayer plate on a linear viscoelastic foundation. Using the Galerkin’s method and quick direct integral method, the dynamical behaviour of the vehiclepavementfoundation coupled system is numerically investigated. Li et al. [5] analyzed nonlinear responses of a seven DOF vehicle moving along a simply supported doublelayer rectangular plate on a nonlinear viscoelastic foundation. The effects of system parameters on vertical acceleration of vehicle body and pavement displacements are also studied.
However, wheel separation has still appeared in very few researches on dynamics of vehicleroad coupled system. Cheng et la [6] investigated the vehiclebridge interaction taking into account the effects of separation and impact on reestablishment of contact between the moving vehicle and bridge. The dynamic response of a stationary structure made of several Timoshenko beams under moving oscillator with the separation and reattachment phenomena is studied by Baeza and Ouyang in [7]. Stăncioiu et al. [8] inspected the vibration of a continuous beam resting on elastic supports excited by a moving twoaxle system with separation. The dynamic response of the system can be very different from the case of not taking separation into account. Base on linear complementary relation linking relative displacement and interaction force between the wheel and the bridge at contact points, Zhu et al. [9] established a model to consider separation in a vehiclebridge system. The influence of velocity, the vehicle to bridge mass ratio and the road roughness on separation is studied. A two DOF quartercar model and a simply supported Euler beam representing the bridge to analyze dynamic interaction of vehiclebridge system with separation are used in [10]. In addition, a nonlinear multiplespring tire is modeled to describe the contact surface of tire and road instead of a point and the form of governing equations does not change whatever the tire separates from the bridge or not.
Vibration of automobiles with wheel separation taken into account has been also considered in some different models by the own authors of this paper. The quartercar models without and with road deformation have been considered in [11, 12]. Vehicle vibrations in halfcar models can be considered in two types, longitudinal and lateral vibrations. In [13], the authors investigated the longitudinal vibration of twoaxle automobiles with wheel separation, but road deformation was ignored. Vibration of twoaxle automobiles in spatial (fullcar) model without taking road deformation has been considered in [14].
The lateral vibration is a common component of vehicle vibration in general and has a significant effect on the passenger’s ride comfort and safety of vehicles. In fact, this component of vibration even occurs more frequently than the longitudinal one. As a continuation of our researches, this paper deals with the lateral vibration of automobiles in halfcar model with road deformation and wheel separation considered. The automobile with dependent suspension is modeled as a vibration system which has two masses and four degrees of freedom, while the deformed road is modeled as an elastic beam which has uniform rectangular crosssection and is simply supported at the two ends and lies on the Kelvin's viscoelastic ground. Although movements of the masses in the physical model take place in a vertical plane, the model is still called as a quasiplanar model because the change in lengths of rectangular contact areas takes place in the direction of movement that is perpendicular to the moving plane, and so does the change in pressure distribution law. The differential equations of motion of the vehicleroad coupled system including a partial differential equation are transformed into a set of all ordinary differential equations by applying the BubnovGalerkin’s method. Some illustrating results obtained from a numerical example are also presented to confirm the occurrence of separation and the difference of dynamic response of vehicle with taking and not taking wheel separation.
2. Vibration model and differential equations of motion
2.1. Assumptions
The lateral vibration model of automobile with dependent suspensions is created with the following assumptions:
– Body and axle of the automobile are assumed to be absolute rigids.
– Each wheel can be represented by a spring  damper couple and behavior of all spring – damper couples in the model is linear.
– Deformed road can be modeled as an elastic beam on Kelvin’s viscoelastic ground and the middle crosssection of the beam passes through the two centers of gravity.
– Lines of action of the elastic and damping forces in the same spring – damper couple coincide and perpendicular to the nominal road surface.
– Each wheel has the width unchanged while loaded and the wheel circumference part off the contact area is still round in the process of vibration.
– The contact area between each wheel and road surface (if exists) is a rectangle whose width is unchanged and equal to the original width of the wheel.
– The road profiles at two wheel tracks and the form of pressure distribution in the contact areas can be described mathematically.
– Kinematic excitation from the road upward to each wheel is calculated with using vertical fluctuation of the center point of the contact area only (not using other points).
– Automobile moves in line with a constant velocity.
2.2. Vibration model of the mechanical system
Fig. 1 shows lateral vibration model of an automobile where road deformation is taken into account. The automobile with dependent suspensions has two masses (body 1 and axle 2) whose centers of gravity are ${S}_{b}$, ${S}_{c}$ and the inertial characteristics are (${M}_{b}$, ${J}_{b}$), (${M}_{c}$, ${J}_{c}$), respectively. Deformed road is modeled as an elastic beam which lies on the Kelvin's viscoelastic ground and is simply supported at the two ends.
Fig. 1Lateral vibration model of automobile taking account of road deformation (viewed from the front)
The beam with uniform rectangular crosssection has the dimensions ${L}_{Bn}$, ${b}_{B}^{x}$ and ${h}_{B}$ for the length, the width and the height, respectively. Material of the beam is assumed to be identical, isotropic and has the mass density $\rho $ and Young’s module $E$. The stiffness and damping coefficient per an area unit of the ground are ${k}_{S}$ and ${c}_{S}$, respectively.
Two suspension springdamper couples are denoted by the notations for their stiffness and damping coefficients as (${k}_{T}$, ${c}_{T}$). Similarly, two wheel springdamper couples are denoted as (${k}_{L}$, ${c}_{L}$).
The points numbered as 1, 2, 1', 2', 1", 2" are the mounting points of the spring damper couples on the body and the axle while points ${D}_{1}$ and ${D}_{2}$ are the centers of the contact areas in cases the loss of contact does not appear.
The effect of movement on vibration of the system is expressed by the changes in height and depth of the two centers of the contact areas (${r}_{D1}$, ${r}_{D2}$) and their derivatives with respect to time (${\dot{r}}_{D1}$, ${\dot{r}}_{D2}$). Deformation of the road is represented by the displacement function of the beam.
Vibration of the automobile and the elastic beam are expressed in terms of the following functions of time as ${u}_{b}={u}_{b}\left(t\right)$, ${\psi}_{b}={\psi}_{b}\left(t\right)$, ${u}_{c}={u}_{c}\left(t\right)$, ${\psi}_{c}={\psi}_{c}\left(t\right)$ and $w=w\left(y,t\right)$ those correspond to the displacements which are measured from the socalled natural position where both wheels still contact with the nominal road surface while all springdamper couples in the model have no deformation.
Figure 1 also presents two typical geometrical parameters of the automobile as $b$ and $c$, and the ycoordinates of points ${D}_{1}$ and ${D}_{2}$.
2.3. Forces acting on the vibration system
Vertical displacements of mounting points of the springdamper couples in the system (points numbered as 1, 2, 1', 2', 1", 2" in Fig. 1) can be determined as follows:
Fig. 2Force diagrams of body and axle of the automobile
After freeing the body and the axle of automobile from the springdamper constraints, we get their force diagrams as shown in Fig. 2.
In the diagrams, $\left\{{G}_{b},{G}_{c}\right\}$ – forces of gravity, $\left\{{F}_{T1},{F}_{T2},{F}_{L1},{F}_{L2}\right\}$ – resultant forces in four springdamper couples, subscripts "1" and "2" imply the forces corresponding to the wheels on the left and the right side, respectively.
Resultant forces in two suspension springdamper couples can be determined as:
Vertical deformations of the wheels in case the loss of contact does not appear can be calculated as:
where ${w}_{D1}$, ${w}_{D2}$ – vertical displacements of the representing beam at the center points (${D}_{1}$, ${D}_{2}$) of the contact areas; ${r}_{D1}$, ${r}_{D2}$ – the heights or depths (with respect to nominal surface of the road) of ${D}_{1}$ and ${D}_{2}$. Basing on Fig. 1, we can determine the ycoordinates of points ${D}_{1}$, ${D}_{2}$ and quantities ${w}_{D1}$, ${w}_{D2}$ as follows:
The resultant forces in the wheel springdamper couples (equal to the contact forces) can be calculated according to the following formulas [1114]:
where ${s}_{1}$, ${s}_{2}$ are the contact state parameters whose values are determined by the sign of the socalled verifying values of contact forces at the two wheels:
Particularly, the $j$th wheel contacts the road, ${\stackrel{}{Q}}_{j}\ge $ 0 and ${s}_{j}=$1; it separates from the road if ${\stackrel{}{Q}}_{j}<$ 0 and ${s}_{j}=$0 ($j=$1, 2).
2.4. Differential equations of motion of the system
2.4.1. Differential equations of motion of automobile
The equations of motion of automobile can be obtained by applying the Newton’s second law to the body and the axle of automobile using the force diagrams in Fig. 2 and relations Eqs. (2), (3), (5):
2.4.2. Differential equation of motion of deformed road
The mentioned equation can be generated by considering the equilibrium of a typical element of the beam. The force diagram of the beam element which lies at coordinate $y$ and has the axle length $dy$ is shown in Fig. 3.
Forces acting on the beam element consist of inertial force ($d{F}_{qt}$), resultant force of Kelvin’s viscoelastic ground ($d{F}_{S}$), force of gravity ($dG$), contact force from the wheels ($d{Q}_{w}$), the shear forces and the bending moments in the two end crosssection $\left\{Q,Q+\left(\partial Q/\partial y\right)dy,M,M+\left(\partial M/\partial y\right)dy\right\}$.
Fig. 3a) Force diagram of a typical beam element, b) the dimension of beam crosssection
a)
b)
The expressions of four first forces can be written as:
where $w=w\left(y,t\right)$ – displacement function of the beam, ${d}_{cj}$ – the length of the $j$th contact area ($j=$1, 2) in the direction of movement, $p\left(x,y,t\right)$ – function of pressure distribution on the upper surface of the beam due to the load coming from the automobile, $g$ – gravitational acceleration.
Fig. 4Description of function px,y,t on beam a) and four types of pressure distribution in xdirection (1 – constant, 2 – parabolic, 3 – cosine, 4 – cosine squared)
The function of pressure distribution $p\left(x,y,t\right)$ is defined over the whole top surface of the beam but really exists in the contact areas of the wheels (Fig. 4). Moreover, this function can be taken only based on observing the fact and using assumptions. This paper still applies the expression of $p\left(x,y,t\right)$ which is proposed in [12] and has the form:
In Eq. (12), the expression of $p\left(x,y,t\right)$ is assumed unchanged in ydirection in the domain of each contact area, ${P}_{j}\left(t\right)$ ($j=$1, 2) are unknown functions of time to be found and ${U}_{j}\left(x\right)$ are prechosen expressions to describe the change of $p\left(x,y,t\right)$ in $x$direction over the length ${d}_{cj}$ of each contact area. Hence, the expressions ${U}_{j}\left(x\right)$ really represent the laws of pressure distribution in the contact areas.
Four typical types of pressure distribution which can be applied in considering vibration of automobiles are proposed in [12] and graphically described in Fig. 4(b). They consist of constant (rectangular), parabolic, cosine and cosinesquared distributions. In the coordinate system $Oxz$ with the origin coinciding with the center of contact area and the $x$axis being the direction of movement, the corresponding expression of ${U}_{j}\left(x\right)$ is 1, 14${\left(x/{d}_{cj}\right)}^{2}$, $\mathrm{c}\mathrm{o}\mathrm{s}(\pi x/{d}_{cj})$, ${\mathrm{c}\mathrm{o}\mathrm{s}}^{2}(\pi x/{d}_{cj})$.
With the expression of function $p\left(x,y,t\right)$ as taken in Eq. (12), the formula of $d{Q}_{w}$ in Eq. (11) becomes:
where:
For the elements which do not lie under the contact areas, we have $d{Q}_{w}=$0 because $p\left(x,y,t\right)=$0.
In case wheel separation does not appear, the values of ${I}_{0}^{\left(j\right)}$ depend on ${d}_{cj}$ and can be determined by using the used type of pressure distribution. For the four types of pressure distribution mentioned above, the values of ${I}_{0}^{\left(j\right)}$ are ${d}_{cj}$, 2${d}_{cj}$/3, 2${d}_{cj}$/$\pi $ and ${d}_{cj}$/2, respectively.
The equilibrium of forces in $z$direction and the moment equilibrium about the center of left crosssection of the element lead to these equations:
where $M=EI\frac{{\partial}^{2}w}{\partial {y}^{2}}$ for EulerBernoulli’s beam theory.
Putting the expression of $Q$ from Eq. (16) into Eq. (15) and using the expressions Eq. (11) and Eq. (13) of forces then taking some needed arrangements, we obtain the differential equation of motion of the beam which represents deformed road:
The differential equations of motion of the mechanical system are the combination of ordinary differential Eqs. (710) and partial differential Eq. (17).
2.5. Simpler cases of differential equations of motion of the system
a) If the loss of contact is ignored, we can fix ${s}_{1}={s}_{2}=$1. Eqs. (7), (8) and (17) are unchanged while Eq. (9), (10) become:
b) If the deformation of road is not taken into account, we have $w\left(y,t\right)=$ 0 for every $y$ and $t$. In this situation, it is allowed to set ${k}_{S}=\infty $, ${c}_{S}=\infty $ and Eq. (17) becomes an identity, Eqs. (7) and (8) are unchanged while Eqs. (9) and (10) have the new forms as:
$={M}_{c}g+{s}_{1}\left({k}_{L}{r}_{D1}+{c}_{L}{\dot{r}}_{D1}\right)+{s}_{2}\left({k}_{L}{r}_{D2}+{c}_{L}{\dot{r}}_{D2}\right),$
c) If road deformation and wheel separation are ignored, the differential equations of motion of the vehicleroad coupled system are reduced to those of only automobile with the addition of setting ${s}_{1}={s}_{2}=$1. In this situation, Eqs. (7) and (8) are still unchanged while Eqs. (9) and (10) respectively becomes:
3. Solving the differential equations of motion
3.1. Transforming the differential equations of motion into ordinary differential equations
The presence of a partial differential equation does not allow to solve the original differential equations of motion of the vehicleroad coupled system until now. In order to overcome this difficulty, here the BubnovGalerkin’s method is applied to transform those differential equations of motion into a system of all ordinary differential equations (ODEs).
According to the BubnovGalerkin’s method, we need to fulfil a succession of operations as follows:
1) Approximating the displacement function of the beam $w\left(y,t\right)$ with a series of $N$ terms as:
where ${T}_{l}\left(t\right)$ are unknown functions of time to be found, ${Y}_{l}\left(y\right)=\mathrm{s}\mathrm{i}\mathrm{n}\left[\right(2l1)\pi y/{L}_{Bn}]$ – shape functions satisfying boundary conditions of the beam, $w(y,t){}_{y=0}=w(y,t\left)\right{}_{y={L}_{Bn}}=0$.
It is important to note that the functions ${Y}_{l}\left(y\right)$ are linearly independent and have the property of orthogonality as:
The approximation of function $w\left(y,t\right)$ as in Eq. (24) allows to define the vector of generalized coordinates:
2) Putting the Eq. (24) of $w\left(y,t\right)$ into Eq. (17) leads to the following equation:
3) Assigning to variable $k$ consecutively the values of the set {1, 2, ..., $N$}. With each value, multiplying Eq. (27) by $\mathrm{s}\mathrm{i}\mathrm{n}\left[\right(2k1)\pi y/{L}_{Bn}]$ then taking the integration of the two sides with respect to variable $y$ from 0 to ${L}_{Bn}$ (over the length of the beam), noting the orthogonality Eqs. (25), we obtain a system of $N$ ODEs as follows:
In order to calculate the integration in Eq. (28), it is important to see the form Eq. (12) of function $p\left(x,y,t\right)$ and its description in Fig. 4(a). The values of ${y}_{1}$, ${y}_{2}$, ${y}_{3}$, ${y}_{4}$ can be expressed in terms of the $y$coordinates of points ${D}_{1}$, ${D}_{2}$ (centers of contact areas) and the width of the wheels as follows: ${y}_{1}={y}_{D1}\mathrm{0,5}{b}_{L}$, ${y}_{2}={y}_{D1}+\mathrm{0,5}{b}_{L}$, ${y}_{3}={y}_{D2}\mathrm{0,5}{b}_{L}$, ${y}_{4}={y}_{D2}+\mathrm{0,5}{b}_{L}$
The relations in Eq. (29) allow to calculate the integration in Eq. (28) as follow:
$+{\int}_{{y}_{3}}^{{y}_{4}}{P}_{2}\left(t\right){I}_{0}^{\left(2\right)}\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2k1)\pi y}{{L}_{Bn}}dy=\frac{2{L}_{Bn}}{(2k1)\pi}\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2k1)\pi {b}_{L}}{2{L}_{Bn}}$
$\times \left[{P}_{1}\left(t\right){I}_{0}^{\left(1\right)}\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2k1)\pi {y}_{D1}}{{L}_{Bn}}+{P}_{2}\left(t\right){I}_{0}^{\left(2\right)}\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2k1)\pi {y}_{D2}}{{L}_{Bn}}\right]$
$=\frac{1}{2}{L}_{Bn}{b}_{B}^{x}\left[{\beta}_{k}^{\left(1\right)}{P}_{1}\left(t\right)+{\beta}_{k}^{\left(2\right)}{P}_{2}\left(t\right)\right],$
where:
Substituting the results in Eq. (30) into Eq. (28) leads to the following equation:
4) Expressing functions ${P}_{1}\left(t\right)$, ${P}_{2}\left(t\right)$ and Eq. (32) in terms of generalized coordinates.
This operation is needed because functions ${P}_{1}\left(t\right)$ and ${P}_{2}\left(t\right)$ are generated in the process of transforming the partial differential Eq. (17) into ODEs. To reach the purpose, we consider the equilibrium of the forces acting on the $j$th wheel ($j=$ 1, 2) in vertical direction (see Fig. 4(b)):
where ${R}_{j}$ – reaction force from road and ${F}_{Lj}$ – resultant force in the $j$th wheel springdamper couple.
The reaction force ${R}_{j}$ can be calculated based on the prechosen pressure distribution function as:
or:
In order to find the expressions of forces ${F}_{Lj}$, we firstly use Eq. (24) to calculate ${w}_{D1}$, ${w}_{D2}$, ${\dot{w}}_{D1}$, ${\dot{w}}_{D2}$ as:
${w}_{D2}={\sum}_{l=1}^{N}{T}_{l}\left(t\right)\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2l1)\pi {y}_{D2}}{{L}_{Bn}}={\sum}_{l=1}^{N}{\chi}_{l}^{\left(2\right)}{T}_{l}\left(t\right),$
${\dot{w}}_{D1}={\sum}_{l=1}^{N}{\dot{T}}_{l}\left(t\right)\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2l1)\pi {y}_{D1}}{{L}_{Bn}}={\sum}_{l=1}^{N}{\chi}_{l}^{\left(1\right)}{\dot{T}}_{l}\left(t\right),$
${\dot{w}}_{D2}={\sum}_{l=1}^{N}{\dot{T}}_{l}\left(t\right)\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2l1)\pi {y}_{D2}}{{L}_{Bn}}={\sum}_{l=1}^{N}{\chi}_{l}^{\left(2\right)}{\dot{T}}_{l}\left(t\right),$
where:
Eqs. (5) and (35) allow to get the expressions of ${F}_{Lj}$:
Using Eqs. (33), (34), (37) we obtain:
$+\left.{c}_{L}\left({\sum}_{l=1}^{N}{\chi}_{l}^{\left(2\right)}{\dot{T}}_{l}\left(t\right)+{\dot{r}}_{D2}{\dot{u}}_{c}b{\dot{\psi}}_{c}\right)\right].$
Putting the expressions of ${P}_{1}\left(t\right)$ and ${P}_{2}\left(t\right)$ from Eq. (38) into Eq. (32) and taking some needed arrangements lead to the following equation:
$+{\sum}_{l=1}^{N}\left[{\delta}_{kl}{c}_{S}+\left({\mu}_{k}^{\left(1\right)}{\chi}_{l}^{\left(1\right)}+{\mu}_{k}^{\left(2\right)}{\chi}_{l}^{\left(2\right)}\right){c}_{L}\right]{\dot{T}}_{l}\left(t\right)({\mu}_{k}^{\left(1\right)}+{\mu}_{k}^{\left(2\right)}){k}_{L}{u}_{c}$
$+({\mu}_{k}^{\left(1\right)}{\mu}_{k}^{\left(2\right)}){k}_{L}b{\psi}_{c}+{\sum}_{l=1}^{N}\left[{\delta}_{kl}{H}_{k}+\left({\mu}_{k}^{\left(1\right)}{\chi}_{l}^{\left(1\right)}+{\mu}_{k}^{\left(2\right)}{\chi}_{l}^{\left(2\right)}\right){k}_{L}\right]{T}_{l}\left(t\right)$
$=\frac{4\rho g{h}_{B}}{\left(2k1\right)\pi}{\mu}_{k}^{\left(1\right)}\left({k}_{L}{r}_{D1}+{c}_{L}{\dot{r}}_{D1}\right){\mu}_{k}^{\left(2\right)}\left({k}_{L}{r}_{D2}+{c}_{L}{\dot{r}}_{D2}\right),k\mathrm{}=\mathrm{}1,\mathrm{}2,\mathrm{}...,\mathrm{}N,$
where:
and ${\delta}_{kl}$ – the Kronecker’s operator defined as ${\delta}_{kl}=$ 1 if $k=l$, and ${\delta}_{kl}=$ 0 if $k\ne l$.
5) Expressing Eqs. (9) and (10) in terms of the generalized coordinates. This purpose can be reached by putting the expressions of functions ${w}_{D1}$, ${w}_{D2}$, ${\dot{w}}_{D1}$, ${\dot{w}}_{D2}$ from Eq. (35) into Eqs. (9), (10) as follows:
$+{\sum}_{l=1}^{N}({s}_{1}{\chi}_{l}^{\left(1\right)}{s}_{2}{\chi}_{l}^{\left(2\right)}){c}_{L}b{\dot{T}}_{l}\left(t\right)2{k}_{T}{c}^{2}{\psi}_{b}({s}_{1}{s}_{2}){k}_{L}b{u}_{c}$
$+[2{k}_{T}{c}^{2}+({s}_{1}+{s}_{2}\left){k}_{L}{b}^{2}\right]{\psi}_{c}+{\sum}_{l=1}^{N}({s}_{1}{\chi}_{l}^{\left(1\right)}{s}_{2}{\chi}_{l}^{\left(2\right)}){k}_{L}b{T}_{l}\left(t\right)$
$={s}_{1}b\left({k}_{L}{r}_{D1}+{c}_{L}{\dot{r}}_{D1}\right)+{s}_{2}b\left({k}_{L}{r}_{D2}+{c}_{L}{\dot{r}}_{D2}\right).$
Now the original differential equations of motion of the mechanical system have been transformed into a system of all ordinary differential Eqs. (7), (8), (41), (42) and (39). These equations can be called as the transformed differential equations of motion (TDEM). The TDEM can be solved numerically because each from them has been expressed in terms of the generalized coordinates.
3.2. Matrix form of TDEM
The TDEM can be written in matrix form as:
where $\overrightarrow{q}$ – vector of generalized coordinates of the form Eq. (26), $\overrightarrow{F}$ – excitation vector; [$M$], [$C$], [$K$] – mass, damping and stiffness matrices, respectively. The two vectors have the size as $\left(4+N\right)\times 1$ while the matrices are all square of order $\left(4+N\right)$.
The excitation vector has the form as:
where:
The mass matrix [$M$] has the diagonal form as:
The stiffness matrix [$K$] contains ($4+N$) rows, each of which has ($4+N$) elements and can be successively written as follows:
$\left\{{K}_{3}\right\}=\{2{k}_{T},\mathrm{0,2}{k}_{T}+({s}_{1}+{s}_{2}){k}_{L},({s}_{1}{s}_{2}){k}_{L}b,[\left({s}_{1}{\chi}_{1}^{\left(1\right)}+{s}_{2}{\chi}_{1}^{\left(2\right)}\right){k}_{L},$
$({s}_{1}{\chi}_{2}^{\left(1\right)}+{s}_{2}{\chi}_{2}^{\left(2\right)}){k}_{L},\dots ,({s}_{1}{\chi}_{N}^{\left(1\right)}+{s}_{2}{\chi}_{N}^{\left(2\right)}){k}_{L}\left]\right\},$
$\left\{{K}_{4}\right\}=\{0,2{k}_{T}{c}^{2},({s}_{1}{s}_{2}){k}_{L}b,2{k}_{T}{c}^{2}+({s}_{1}+{s}_{2}){k}_{L}{b}^{2},\left({s}_{1}{\chi}_{1}^{\left(1\right)}{s}_{2}{\chi}_{1}^{\left(2\right)}\right){k}_{L}b,$
$({s}_{1}{\chi}_{2}^{\left(1\right)}{s}_{2}{\chi}_{2}^{\left(2\right)}){k}_{L}b,\dots ,({s}_{1}{\chi}_{N}^{\left(1\right)}{s}_{2}{\chi}_{N}^{\left(2\right)}){k}_{L}b\left]\right\},$
$\left\{{K}_{4+k}\right\}=\{\mathrm{0,0},\left({\mu}_{k}^{\left(1\right)}+{\mu}_{k}^{\left(2\right)}\right){k}_{L},\left({\mu}_{k}^{\left(1\right)}{\mu}_{k}^{\left(2\right)}\right){k}_{L}b,[{\delta}_{k1}{H}_{k}+\left({\mu}_{k}^{\left(1\right)}{\chi}_{1}^{\left(1\right)}+{\mu}_{k}^{\left(2\right)}{\chi}_{1}^{\left(2\right)}\right){k}_{L},$
${\delta}_{k2}{H}_{k}+({\mu}_{k}^{\left(1\right)}{\chi}_{2}^{\left(1\right)}+{\mu}_{k}^{\left(2\right)}{\chi}_{2}^{\left(2\right)}){k}_{L},\dots ,{\delta}_{kN}{H}_{k}+({\mu}_{k}^{\left(1\right)}{\chi}_{N}^{\left(1\right)}+{\mu}_{k}^{\left(2\right)}{\chi}_{N}^{\left(2\right)}){k}_{L}\left]\right\},$
$k=\mathrm{1,2},\dots ,N.$
The damping matrix [$C$] has the form similar to that of the stiffness matrix [$K$] so we can obtain matrix [$C$] from matrix [$K$] by substituting the notations $\left\{{H}_{k},{k}_{T},{k}_{L}\right\}$ with these ones $\left\{{c}_{S},{c}_{T},{c}_{L}\right\}$, respectively.
3.3. Initial conditions
In almost problems on vibration of automobiles moving on horizontal roads, the initial conditions commonly relate to the fact that the vehicle is running on a definitely even road surface then entering an uneven one. In this situation, theoretically, there is no vibration occurring until the initial point of time, so the values of generalized coordinates, velocities, accelerations and all quantities concerned can be determined from the static state where the automobile lies on a horizontal road with even surface. Therefore, the initial conditions for this situation can be written as follows:
The value of ${\overrightarrow{q}}_{0}$ can directly be deduced from Eq. (43) using the data in Eq. (48):
where $[K{]}_{0}$, ${\overrightarrow{F}}_{0}$ – the values of the stiffness matrix and excitation vector at the time point $t=$ 0.
At the initial point of time, because the vehicle is running on the even road surface so we have:
Putting data Eq. (50) into Eq. (45) we obtain the values of excitation vector at the initial time point:
In order to determine $[K{]}_{0}$, we firstly generate the static equilibrium equations of the automobile and deduce:
Now we can determine the values of quantities concerned with the elements of matrix $[K{]}_{0}$, such as static deformations of the wheels ($\mathrm{\Delta}{z}_{L1}^{0}$, $\mathrm{\Delta}{z}_{L2}^{0}$), the lengths of two contact areas (${d}_{c1}^{0}$, ${d}_{c2}^{0}$), etc.
Once the value of ${\overrightarrow{q}}_{0}$ has been found, we can calculate displacement function of the beam in static state by using Eq. (24) as:
where ${T}_{l}^{0}$ – elements from 5th to (4+$N$)th of vector ${\overrightarrow{q}}_{0}$.
3.4. Procedure for numerically solving the TDEM
Matrix form Eq. (43) of TDEM contains matrices [$K$] and [$C$] which continuously change with respect to time. Hence, at every step in the process of numerically solving this equation, we need to recalculate the elements of these matrices. A procedure for numerically solving Eq. (43) can be proposed as follows:
1) Assigning the values of quantities concerned with automobile (${M}_{b}$, ${J}_{b}$, ${M}_{c}$, ${J}_{c}$, ${k}_{T}$, ${k}_{L}$, ${c}_{T}$, ${c}_{L}$, $b$, $c$, ${r}_{w}$, ${b}_{L}$), viscoelastic ground (${k}_{S}$, ${c}_{S}$), elastic beam ($\rho $, $E$, ${L}_{Bn}$, ${b}_{B}^{x}$, ${h}_{B}$, $I$), velocity of movement ($V$), number $N$ of terms in the series Eq. (24) and acceleration of gravity ($g$).
2) Describing road profile as functions of time ${r}_{D1}\left(t\right)$, ${r}_{D2}\left(t\right)$, ${\dot{r}}_{D1}\left(t\right)$, ${\dot{r}}_{D2}\left(t\right)$ and choosing the type of pressure distribution function $U\left(x\right)$ for consideration.
3) Calculating the values of ${y}_{D1}$, ${y}_{D2}$, ${H}_{k}$, ${\chi}_{k}^{\left(1\right)}$,${\chi}_{k}^{\left(2\right)}$.
4) Setting the initial conditions as in Section 3.3.
5) Specifying the time interval for consideration [0, ${t}_{max}$], the time step of computation $\mathrm{\Delta}t$ and the point of time ${t}_{0}$ at which the uneven road surface really starts.
6) Assigning ${s}_{1}={s}_{2}=$1, ${\ddot{\overrightarrow{q}}}_{0}=0$, ${\dot{\overrightarrow{q}}}_{0}=0$ for the initial point of time and calculating the initial value ${\overrightarrow{q}}_{0}$ of vector $\overrightarrow{q}$.
7) Assigning $i:=$0, ${t}_{i}:=$0, ${\ddot{\overrightarrow{q}}}_{i}:={\ddot{\overrightarrow{q}}}_{0}$, ${\dot{\overrightarrow{q}}}_{i}:={\dot{\overrightarrow{q}}}_{0}$, ${\overrightarrow{q}}_{i}:={\overrightarrow{q}}_{0}$ for the initial point of computation.
8) Determining the values ${\ddot{\overrightarrow{q}}}_{i+1}$, ${\dot{\overrightarrow{q}}}_{i+1}$, ${\overrightarrow{q}}_{i+1}$ at the time point ${t}_{i+1}=\left({t}_{i}+\mathrm{\Delta}t\right)$ of the ($i+1$)th point of computation using a suitable method such as Newmark’s, RungeKutta’s, etc.
9) Calculating the values of quantities ${w}_{Dj}$, ${\dot{w}}_{Dj}$, ${r}_{Dj}$, ${\dot{r}}_{Dj}$, ${u}_{{j}^{\text{'}\text{'}}}$, ${\dot{u}}_{{j}^{\text{'}\text{'}}}$ ($j=$ 1, 2) at the ($i+1$)th point of computation using Eq. (27), functions ${r}_{D1}\left(t\right)$, ${r}_{D2}\left(t\right)$, ${\dot{r}}_{D1}\left(t\right)$, ${\dot{r}}_{D2}\left(t\right)$ and the results obtained in step 8.
10) Calculating verifying values of contact forces ${\stackrel{}{Q}}_{1}$, ${\stackrel{}{Q}}_{2}$ according to Eq. (6) and deducing the values of quantities ${s}_{1}$, ${s}_{2}$, ${\mu}_{k}^{\left(1\right)}$, ${\mu}_{k}^{\left(2\right)}$.
11) Determining the values ${\left[K\right]}_{i+1}$, ${\left[C\right]}_{i+1}$, ${\overrightarrow{F}}_{i+1}$ of matrices [$K$], [$C$] and vector $\overrightarrow{F}$ at the ($i+1$)th point of computation.
12) Assigning $i:=i+1$, ${t}_{i}:={t}_{i}+\mathrm{\Delta}t$ and repeating the process of computation, starting from step 8. The process of computation stops when condition ${t}_{i}>{t}_{max}$ is reached.
The results directly obtained are functions of time such as generalized coordinates, generalized velocities, generalized accelerations, contact forces, the total time of losing contact, etc. which reflect vibration response of the automobile in particular conditions which concern with of the road profile and the velocity of movement.
4. Example for illustration
This section presents some typical results obtained from numerical consideration on vibration of the vehicleroad coupled system where the initial conditions mentioned in Section 3.3 are applied. The results coming from two cases of taking and not taking wheel separation into account are compared to show the differences.
The situation under consideration is that the automobile successively passes a single bump on the right track and another one on the left track after having passed a distance ${x}_{0}$ measured from the initial time point ($t=$0) as seen in Fig. 5. Profiles of the two bumps are chosen as parabolas whose lengths and heights are (${l}_{1}$, ${h}_{1}$) and (${l}_{2}$, ${h}_{2}$), respectively. The distance between the two bumps in the direction of movement is denoted as $d$.
Fig. 5 is a general description of some cases. If setting ${h}_{1}=$0 or ${h}_{2}=$0, the excitation has the type of a single pulse on the left or the right track only. If taking $d<$0, the automobile reaches the bump on the left track before reaching the bump on the right one. If assigning $d=$0, ${l}_{1}={l}_{2}$ and ${h}_{1}={h}_{2}$, we have the situation of a speed bump which crosses over the whole width of road. If choosing ${x}_{0}=$0, the initial point of time ($t=$0) is also the time when the vehicle starts passing the left bump.
Fig. 5Geometrical description of road profile
The input data for numerical computation are taken as:
– Parameters concerned with the automobile [15]: $b=$ 0.90 m, $c=$ 0.60 m, ${r}_{w}=$ 0.45 m, ${b}_{L}=$ 0.25 m, ${M}_{b}=$ 2150 kg, ${J}_{b}=$ 650 kg.m^{2}, ${M}_{c}=$ 660 kg, ${J}_{c}=$ 720 kg.m^{2}, ${k}_{T}=$ 250000 N/m, ${k}_{L}=$ 800000 N/m, ${c}_{T}=$ 1500 N.s/m, ${c}_{L}=$ 62000 N.s/m.
– Quantities belonging to the ground and the beam [16]: ${k}_{S}=$ 48000000 N/m^{3}, ${c}_{S}=$ 30000 N.s/m^{3}, $E=$ 1.6×10^{9} N/m^{2}, $\rho =$ 2500 kg/m^{3}, ${L}_{Bn}=$ 15 m, ${b}_{B}^{x}=$ 0.45 m, ${h}_{B}=$ 0.50 m.
– Excitation data (geometries of bumps): ${h}_{1}=$ 0.10 m, ${l}_{1}=$ 0.75 m, ${h}_{2}=$ 0.12 m, ${l}_{2}=$ 0.80 m, $d=$1.0 m.
– Velocity of movement of automobile: $V=$10 km/h.
– Parameters involved computation: ${t}_{max}=$ 5 s, ${t}_{0}=$ 0.5 s, $\mathrm{\Delta}t=$0.0001 s, $N=$ 5 (${x}_{0}=Vt$).
– Type of pressure distribution in contact areas: cosine.
The obtained results are graphically showed in Fig. 69 where the changes in some vibration characteristics with respect to time are plotted. Some comparisons between the two cases of taking and not taking account of wheel separation are given in Figs. 6, 7 and 8.
Fig. 6Linear vibration of automobile body
Fig. 7Angular vibration of automobile body
Fig. 8Vertical acceleration of automobile body
Fig. 9The changes in contact forces at two wheels
It can be seen from the figures that:
a) Differences in behaviour of the mechanical system between the two cases of taking and not taking account of wheel separation are clearly (see Figs. 6, 7 and 8).
b) In the situation of consideration, the amplitude of vibration in the case of not taking wheel separation is greater than the amplitude in the other case. The reason may concern with the ways of transmitting potential energy which has been accumulated by the springs in compressing periods to the body in stretching ones.
c) Wheel separation really occurs in the situation of consideration. This can be identified by flat pieces (on the abscissa) of the curves in Fig. 9. The calculated total times of losing contact are 0.0416 s and 0.0458 s for the right and the left wheels, respectively.
d) After each period of wheel separation, the contact forces rapidly increase (Fig. 9) due to converting kinematic energy of falling into potential energy of the springs.
e) Wheel separation is quite easy to occur. This can be recognized by looking at the values of the parameters which describe the road profile and especially the small value of velocity of movement.
5. Conclusions
The article has created a physical model for considering the quasiplanar lateral vibration of the automobiles with dependent suspensions where all three wheel separation, road deformation and the change in dimension of contact areas have been taken into account. The original differential equations of motion of the vehicleroad coupled system with the presence of a partial differential equation have been transformed into a set of all ordinary differential equations by applying the BubnovGalerkin’s method. The authors have also proposed a procedure for numerically solving the transformed differential equations of motion. Vibration of automobile, according to the model introduced above, has been considered in the case in which an automobile moves through two bumps (each bump on one track) of parabolic profiles. The results obtained from numerical computation have proved that the difference between taking and not taking wheel separation into account is significant and the loss of contact is easy to appear while automobiles run on rough roads. Hence, one can affirm the need of taking wheel separation into account in considering vertical vibration of automobiles.
References

Cebon D. Interaction between heavy vehicles and roads. SAE Technical Paper 930001, 1993, https://doi.org/10.4271/930001.

Cole D. J. Truck suspension design to minimize road damage. Proceedings of the Institution of Mechanical Engineers, Vol. 210, Issue 22, 1996, p. 95107.

Zhu Q., Ishitobi M. Chaos and bifurcations in a nonlinear vehicle model. Journal of Sound and Vibration, Vol. 275, Issues 35, 2004, p. 11361146.

Yang S., Li S., Lu Y. Investigation on dynamical interaction between a heavy vehicle and road pavement. Vehicle System Dynamics, Vol. 48, Issue 8, 2010, p. 923944.

Li S., Lu Y., Li H. Effects of parameters on dynamics of a nonlinear vehicleroad coupled system. Journal of Computers, Vol. 6, No. 12, 2011, p. 26562661.

Cheng Y. S., Au E. T. K., Cheung Y. K., Zheng D. Y. On the separation between moving vehicles and bridge. Journal of Sound and Vibration, Vol. 222, Issue 5, 1999, p. 781801.

Baeza L., Ouyang H. Dynamics of a truss structure and its movingoscillator exciter with separation and impactreattachment. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, Vol. 464, Issue 2098, 2008, p. 25172533.

Stăncioiu D., Ouyang H., Mottershead J. E. Vibration of a continuous beam with multiple elastic supports excited by a moving twoaxle system with separation. Meccanica, Vol. 44, Issue 3, 2008, p. 293303.

Zhu D. Y., Zhang Y. H., Ouyang H. A linear complementarity method for dynamic analysis of bridges under moving vehicles considering separation and surface roughness. Computers and Structures, Vol. 154, 2015, p. 135144.

Zhang Y., Zhao H., Lie S. T. A nonlinear multispring tire model for dynamic analysis of vehiclebridge interaction system considering separation and road roughness. Journal of Sound and Vibration, Vol. 436, 2018, p. 112137.

Cuong P. M., Dung T. Q., Nam L. H. Consideration on vibration of automobiles in quartercar model with the loss of contact taken into account. Journal of Sciences and Technology, Vol. 183, 2017, p. 6473.

Ham V. C., Cuong P. M., Dung T. Q. Consideration of the problem about vibration of automobiles in onefourth model with taking road deformation and the loss of contact into account. Journal of Vibroengineering, Vol. 22, Issue 4, 2020, p. 945958.

Ham V. C., Cuong P. M., Dung T. Q. Consideration on vibration of automobiles in planar model taking account of the loss of contact. Proceedings of the National Conference on Mechanics and Transportation Engineering, Ho Chi Minh City, Vietnam, Vol. 2, 2017, p. 326332, (in Vietnamese).

Ham V. C., Cuong P. M. Consideration on vibration of automobiles in spatial model with the loss of contact taken into account. International Journal of Applied Engineering Research, Vol. 15, Issue 6, 2020, p. 594599.

Lap V. D. Handbook for looking up technical specifications of automobiles. Le Quy Don Technical University, 2004, (in Vietnamese).

Shaopu Y., Liqun C., Shaohua L. Dynamics of VehicleRoad Coupled Systems. First Edition, SpringerVerlag Berlin Heidelberg, 2015.