Consideration of the problem about vibration of automobile in one fourth model with taking road deformation and the loss of contact into account

2.1. Contact states, geometry of contact area and relations of quantities

Fig. 2 shows three possible relative positions between a rounded wheel and a road with planar surface. In Fig. 2(a), the wheel lies in contact state with the road surface and it is deformed under the contact force $Q$. The shape of the contact area is accepted to be rectangular [6]. For simplicity, we add here two assumptions as follows: a) the off-contact-area part of wheel profile is still exactly rounded with the radius unchanged; b) the width of the wheel is not changed under loading.

With the assumptions mentioned above, the contact area in this relative position has the dimensions as $b_L\times d_C$ with the length $d_C$ determined as:

where $\mathrm{\Delta }z$ is vertical deformation of the wheel (or tyre).

Once the linearity in behavior of the wheel is accepted, between the contact force $Q$ and vertical deformation $\mathrm{\Delta }z$ exists the following relation:

where $k$ and $c$ are stiffness and damping coefficient of the spring-damper couple which represents the elastic and damping properties of the wheel, respectively.

Eqs. (1) and (2) describe the relations of wheel vertical deformation with the contact area length and the contact force.

The second relative position of the wheel-road couple is shown in Fig. 2(b). In this situation, the wheel still lies in the contact state with the road surface but has no deformation. The contact area now reduces to a straight line segment with the length $b_L$ and the contact force is equal to zero.

In Fig. 2(c), the wheel completely loses contact with the road surface and it is not deformed. The contact area is absent and the contact force is also equal to zero.

Because the value of wheel vertical deformation and that of the contact force are simultaneously differ from zero in the first case and simultaneously equal to zero in the two others, so we can unify the equations which express the relation of contact force $Q$ and wheel vertical deformation $∆z$ for three cases in Fig. 2 as follows:

Parameter $s$ in Eq. (3) is taken the value of 0 for cases where $∆z=$0, $Q=$0, and taken the value of 1 in cases $∆z\ne$0, $Q\ne$0. Because $s$ reflects the state of contact or separation of the wheel-road couple, it is called as the contact state parameter. The introduction of the contact state parameter facilitates making the differential equations of motion of the considered mechanical system as presented in the next section.

In the process of vibration of automobiles, the values of $Q$ and $∆z$ change with respect to time. They can be determined only by solving the differential equations of motion of the mechanical system under consideration.

Fig. 2Three possible relative positions between a rounded wheel and a planar road surface: a) contact with wheel deformed; b) contact with wheel undeformed; c) no contact

2.2. Distribution of pressure in the contact area

The distribution of pressure in the contact area between the wheel and the road surface depends on many factors such as the geometrical and physical properties of the tire, those of the road and the speed of movement, etc. In fact, the distribution of pressure is very complicated [6]. However, for simplicity in considering vertical vibrations of automobiles, we assume that the distribution of pressure is symmetrical about the coordinate plane $Oyz$, a vertical plane passing through the axle of the wheel as shown in Fig. 3 (where the axis $Oy$ is not described), and moreover, the pressure distribution function does change only in the $x$-direction, not change in the $y$-direction [6, 8].

Fig. 3Pressure distribution laws assumed: 1 – even or rectangle, 2 – parabola, 3 – cosine, 4 –squared cosine)

Accordingly, the pressure distribution function in the contact area can be written in the following form:

where $P\left(t\right)$ is a time function that reflects the change of pressure with respect to time, $U\left(x\right)$ is an $x$ – variable function which represents the changing of pressure in the $x$ direction, $U\left(x\right)$ is a dimensionless function that depends on the pressure distribution law chosen. The four typical forms of function $U\left(x\right)$ proposed in the reference [8] are plotted in Fig. 3 and consist of the following laws:

– Even distribution: $U\left(x\right)=1\left(\forall x\right)$;

– Parabola distribution: $U\left(x\right)=1-{\left(\frac{2x}{d_c}\right)}^2$.

– Cosine distribution: $U\left(x\right)=\mathrm{c}\mathrm{o}\mathrm{s}\frac{\pi x}{d_c}$.

– Squared cosine distribution: $U\left(x\right)=\mathrm{c}\mathrm{o}{\mathrm{s}}^2\frac{\pi x}{d_c}=\frac{1}{2}\left(1+\cos \frac{2\pi x}{d_c}\right)$.

From the force equilibrium in vertical direction, we have:

where $R_C$ is the reaction force from the road upwards to the wheel and:

For rectangle, parabola, cosine and squared cosine distributions, the values of $I_0$ can be obtained from Eq. (6) is $d_C$, $2d_C/3$, $2d_C/\pi$ and $d_C/2$, respectively.

The Eqs. (3) and (5) allow to deduce a relation between function $P\left(t\right)$ and vertical deformation $∆z$ of the wheel:

It is noted that when the wheel is not in contact with the road surface, both $s$ and $I_0$ are simultaneously equal to zero. At this time, although the Eq. (7) of $P\left(t\right)$ has an infinite form 0/0, but the value of contact force $Q$ is equal to zero.

3.1. One-fourth vibration model of automobiles

An one-fourth vibration model of automobiles with road deformation taken into account is shown in Fig. 4(a). The automobile is modeled as a vertical vibration system of two masses and the deformable road is model as an elastic beam which is simply supported at the two ends and lies on the Kelvin’s visco-elastic ground.

Vertical vibration of the automobile takes place at the middle of the beam. The effect of automobile movement on vibration is represented by time changing in depth or height of the center point of the contact area due to road roughness. Deformation of the road with respect to time is represented by vibration of the beam.

The beam with the length $L$ has uniform rectangular cross-section whose width and height are denoted as $b$ and $h$, respectively. Material of the beam is assumed to be homogeneous, isotropic, and has mass density $\rho$ and Young’s modulus $E$.

The position presented in Fig. 4(a) of the vibration system is a so-called natural position where the wheel is in contact with the road surface and all the springs in the model are completely free. In this figure, $m_1$ and $m_2$ are the suspended and unsuspended masses of the vehicle in one-fourth model; $\left\{k_1,c_1\right\}$ is the spring-damper couple which represents the vehicle suspension; $\left\{k_2,c_2\right\}$ is the spring-damper couple representing the elastic and damping properties of the wheel; $k_S$ and $c_S$ are the ground stiffness and damping coefficient per unit area. For simplicity, $\left\{k_1,c_1\right\}$ and $\left\{k_2,c_2\right\}$ are also denoted for the stiffness and damping coefficient of the corresponding spring-damper couples. Vibrations of the two masses and the elastic beam are respectively expressed by the displacement functions $u_1=u_1\left(t\right)$, $u_2=u_2\left(t\right)$ and $w=w\left(x,t\right)$ those are measured from the natural position.

Fig. 4One-fourth vibration model of automobile with road deformation taken into account: a) vibration model and characteristic parameters; b) force diagram of masses

3.2. Differential equations of motion of the vibration system

3.2.2. Differential equation of motion of the road

The differential equation of motion of the road can be made by considering the equilibrium of a typical beam element which has the length of $dx$ and locates at the position determined by $x$ coordinate as shown in Fig. 5.

Fig. 5Force diagram of a typical beam element in vibration model

Forces acting on the element include:

– Axial forces $\left\{N,N+\left(\frac{\partial N}{\partial x}\right)dx\right\}$, shear forces $\left\{Q,Q+\left(\frac{\partial Q}{\partial x}\right)dx\right\}$ and bending moments around the y-axis $\left\{M,M+\left(\frac{\partial M}{\partial x}\right)dx\right\}$ at the left and the right sections.

– The force of gravity $dG$, $dG=rgbhdx$ ($g$ – gravitational acceleration).

– Contact force $dQ$ which may exists in the contact area, $dQ=p(x,\mathrm{}t)b_{L}dx$. If the element is not in contact with the wheel then $p\left(x,t\right)$ and $dQ$ are equal to zero.

– Inertial force $dF$, $dF=\rho bh\frac{{\partial }^{2}w}{\partial t^2}dx$.

– Reaction force from the ground $d{F}_S$, $d{F}_{\mathrm{S}}=k_{\mathrm{S}}w(x,\mathrm{}t)bdx+c_{\mathrm{S}}\frac{\partial w(x,\mathrm{}t)}{\partial t}bdx$.

The force equilibrium in horizontal and vertical directions, and the moment equilibrium about the center point of the left cross-section of the element combining with neglecting the infinitesimal quantities lead to these equations:

14

$\frac{\partial N}{\partial x}=0,$

15

$\rho h\frac{{\partial }^{2}w}{\partial t^2}+c_S\frac{\partial w}{\partial t}+k_{S}w\left(x,\mathrm{}t\right)-\frac{1}{b}\frac{\partial Q}{\partial x}-p\left(x,\mathrm{}t\right)=\rho gh,$

16

$Q=\frac{\partial M}{\partial x}.$

Eq. (14) gives $N=$0 for all $x$, $t$ since no external force acts in the $x$-direction.

By substituting Eq. (16) into (15) and noting that $M=-EI\left(\frac{{\partial }^{2}w}{\partial x^2}\frac{{\partial }^{2}w}{\partial x^2}\right)$ [6], one obtains the differential equation of motion of the beam:

17

$\rho h\frac{{\partial }^{2}w}{\partial t^2}+c_S\frac{\partial w}{\partial t}+k_{S}w\left(x,\mathrm{}t\right)+\frac{1}{b}EI\frac{{\partial }^{4}w}{\partial x^4}-p\left(x,\mathrm{}t\right)=\rho gh,$

where $E$ – the Young’s modulus of beam material, $I$ – the inertial moment of bending around the $y$-axis of beam cross-section, $I=b{h}^3/12$.

The differential equations of motion of the vehicle-road system are combination of ordinary differential Eqs. (12), (13) and partial differential Eq. (17).

If the loss of contact is not taken into account, the differential equations of motion of the vehicle-road coupled system can be obtained from the equations mentioned above by fixing $s=$1. In case road deformation is neglected, $w\left(x,t\right)=$0 for all $x$ and $t$, Eq. (17) is obviously satisfied and the differential equations of motion of the mechanical system reduces to the differential equations of motion Eqs. (12), (13) of the automobile only.

3.3. Transforming the differential equations of motion to a system of all ordinary differential equations (ODE)

In order to obtain the response of the automobile under consideration according to some initial conditions, it is necessary to solve the differential equations of motion mentioned above. Solution given here to reaching that purpose is applying the Bubnov – Galerkin’s method to transform the partial differential Eq. (17) to a system of all ordinary differential equations.

Firstly, the function $w\left(x,t\right)$ is approximated by the series of $N$ terms as:

where $T_l\left(t\right)$ ($l=$ 1, 2, ..., $N$) are time functions to be found, and $X_l\left(x\right)=\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2l-1)\pi x}{L}$ are functions those satisfy the boundary conditions of the beam, ${w\left(x,t\right)‖}_{x=0}={w\left(x,t\right)‖}_{x=L}=$ 0, and match the symmetry of $w\left(x,t\right)$ about the middle cross-section of the beam.

The functions $X_l\left(x\right)$ in Eq. (18) are also linearly independent and have the property of orthogonality as follows:

Substituting the Eq. (18) of $w\left(x,t\right)$ and the Eq. (4) of $p\left(x,t\right)$ into Eq. (17) leads to the following equation:

In succession, multiplying both two sides of Eq. (20) by the expression: $X_k\left(x\right)=\mathrm{s}\mathrm{i}\mathrm{n}\frac{(2k-1)\pi x}{L}$ ($k=$ 1, 2, ..., $N$).

Then integrating the two sides of the obtained equation over the entire length of the beam, paying attention to the orthogonality Eq. (19), we obtain a system of $N$ ordinary differential equations of the form:

where $k=$1, 2, ..., $N$ and:

Because the functions to be determined are $u_1\left(t\right)$, $u_2\left(t\right)$ and $T_l\left(t\right)$ ($l=$1, 2, ..., $N$), so it is needed to express $P\left(t\right)$ in Eq. (21) in terms of those unknown functions. The expression of $P\left(t\right)$ can be obtained by using the Eq. (9) of $\mathrm{\Delta }z$ in Eq. (7):

The expressions of $w_C$ and ${\dot{w}}_C$ in Eq. (23) can also be written in terms of the unknown functions by using the Eq. (18) of $w_C$ and Eq. (10):

Substituting Eqs. (24) and (25) into Eq. (23) and putting the result into Eq. (21), then taking some arrangements, we obtain the following equations:

where:

With the Eqs. (24), (25) of $w_C$ and ${\dot{w}}_C$, Eq. (13) is also written in terms of the unknown functions ($u_1$, $u_2$ and $T_l\left(t\right)$) as follows:

The original differential Eqs. (12), (13), (17) has now been transformed into a system of all ordinary differential Eq. (12), (28) and (26) those can be called as the transformed differential equations and written in matrix form:

where $\overrightarrow{q}$ is the vector of generalized coordinates, $\overrightarrow{F}$ is the vector of excitation force and [$M$], [$C$], [$K$] are the mass, damping and stiffness matrices, respectively. The dimension of the two vectors is ($N+2$)×1 and that of the three matrices is ($N+2$)×($N+2$). Those vectors and matrices can be written in detail as follows:

3.4. Determination of static displacements of the masses and elastic beam

Determination of static displacements of the components is needed in case the vehicle is moving on a completely smooth road when entering a rough road and static displacements are chosen as the corresponding initial conditions to consider vertical vibration of the mechanical system.

Conventionally, the subscript or superscript “0” is used to refer to the values of the corresponding quantities in static state. Because the phenomenon of losing contact and vertical vibrations of the whole system do not occur in static state so $s=$1 and ${\ddot{\overrightarrow{q}}}_0=\overrightarrow{0}$, ${\dot{\overrightarrow{q}}}_0=\overrightarrow{0}$. At this time, the spring $k_2$ is subjected to the force equal to the total weight of the two masses $m_1$ and $m_2$, so its static deformation is $\mathrm{\Delta }{z}_0=\frac{\left(m_1+m_2\right)g}{k_2}$. As the result, the corresponding length of the contact area can be calculated as:

The values $I_k^0$ of the quantities $I_k$ in Eq. (22) can be determined according to the pre-chosen law of pressure distribution. Then, we can calculate the value ${\left[K\right]}_0$ of the matrix [$K$] by using Eqs. (27) and (34), then deduce static displacements ${\overrightarrow{q}}_0$ based on Eq. (29) as follows:

where ${\overrightarrow{F}}_0$ is the value of $\overrightarrow{F}$ in static state. It is determinned from Eq. (31) by setting $r_C={\dot{r}}_c$= 0.

The displacement expression of the beam is deduced from Eq. (18):

where the values of quantities $T_l^0$ are last $N$ components of vector ${\overrightarrow{q}}_0$ which has been calculated as in Eq. (36).

3.5. Method for numerically solving the differential equations of motion

Differential Eq. (29) are ODEs with the matrices [$C$] and [$K$] depending on time. Let’s consider the case where the vehicle is moving on a smooth road when entering the rough road and the time point $t=$0 is chosen as the time point of transition between the two road profiles. In this situation, the relevant initial values are taken as those of static state.

The procedure for numerically solving the differential equations of motion Eq. (29) can be described as follows:

1) Assigning values to the quantities $g$, $m_1$, $m_2$, $k_1$, $k_2$, $k_S$, $c_1$, $c_2$, $E$, $b$, $h$, $I$, $L$, $r_0$, $b_L$ and the vehicle speed of movement.

2) Describing the profile of road surface as the time functions $r_C=r_C\left(t\right)$, ${\dot{r}}_c={\dot{r}}_c\left(t\right)$.

3) Choosing the pressure distribution function $U\left(x\right)$ and the value of $N$, the number of terms of the series which is used to approximate function $w\left(x,t\right)$.

4) Selecting the computing time interval [0, $t_{max}$] and the computing time step $\mathrm{\Delta }t$.

5) Assigning $i=$ 0, $t_i=$ 0, $s=$1, ${\overrightarrow{q}}_i={\overrightarrow{q}}_0$, ${\dot{\overrightarrow{q}}}_i=\overrightarrow{0}$ for starting step of computing.

6) Calculating the values of $d_C^0$, I₀ and $I_k^0$ by using Eqs. (35), (6) and (22).

7) Calculating the values of matrices ${\left[C\right]}_i$, ${\left[K\right]}_i$ according to Eqs. (33), (34) and the exciting vector ${\overrightarrow{F}}_i$ according to Eq. (31).

8) Calculating vectors ${\overrightarrow{q}}_{i+1}$, ${\dot{\overrightarrow{q}}}_{i+1}$ and ${\ddot{\overrightarrow{q}}}_{i+1}$ using a known method such as Runge-Kutta’s method.

9) Determining displacement function $w_{i+1}\left(x,t\right)$of the beam according to Eq. (18) and deducing the values of $w_C^{i+1}$, ${\dot{w}}_C^{i+1}$.

10) Calculating $Q^{i+1}=k_2[u_2^{i+1}-(w_C^{i+1}+r_C^{i+1}\left)\right]+c_2[{\dot{u}}_2^{i+1}-({\dot{w}}_C^{i+1}+{\dot{r}}_C^{i+1}\left)\right]$ and following one from the two possibilities:

– If $Q^{i+1}>$0 then the loss of contact does not occur, so $s=$1. At this point of time, it is needed to calculate the values of $d_C^{i+1}$, $I_0^{i+1}$, $I_k^{i+1}$.

– If $Q^{i+1}\le$0 then the loss of contact occurs or starts to occur and $s=$0. Now $p\left(x,t\right)=$0 and the terms concerning $p\left(x,t\right)$ in Eq. (26) are absent and ${\mu }_k^{i+1}=$0.

11) Assigning $i=i+1$, $t_i=t_i+∆t$ and repeating the process of calculation from step 7. The process of calculation ends when reaching the condition $t_i>t_{max}$.

At the end of the process of calculation, we obtain the following functions $u_1\left(t\right)$, $u_2\left(t\right)$, $T_1\left(t\right)$, $T_2\left(t\right)$, …, $T_N\left(t\right)$ together with their first and second derivatives. The displacement function $w\left(x,t\right)$ which describes the road deformation at each time point is determined by applying Eq. (16).