*3.3. Wheel–Road Displacement Coupled Vibration*

It has been assumed that the vehicle is always in contact with the pavement during driving, and the vertical displacement of the pavement is divided into pavement flatness *zp*(*xi*) and pavement vibration vertical displacement *zr*(*xi*); then, the vertical displacement of the wheel can be expressed as:

$$z\_{li} = z\_{l'}(\mathbf{x}\_i) + z\_{l'}(\mathbf{x}\_i) \tag{18}$$

where *zp*(*xi*) is the road surface flatness value of tire *i* (*i* represents the wheel position) at the horizontal direction x of the road. Substituting Equation (14) into Equation (18), we get:

$$z\_{li} = \sum\_{n=1}^{N} A\_n \phi\_n(\mathbf{x}\_i) + z\_{l^\*}(\mathbf{x}\_i) \tag{19}$$

The contact between vehicle wheel and the road surface can be simplified according to Hertz's contact theory. Consider the coupling relationship between the wheel and the road surface to be a nonlinear contact force, as shown in Figure 7.

**Figure 7.** Wheel–road coupling vibration diagram.

According to the contact relation in Hertz's law, when further simplified the vertical contact force between the wheel and the road surface is expressed as:

$$F(t) = k\mu (Z\_w - Z\_r) \tag{20}$$

### *3.4. Establishment of Vehicle–Road Coupling Dynamic Analysis Model*

As mentioned above, the vehicle vibration balance equation and road vibration balance equation are combined to obtain the vibration coupling relationship model under the condition of displacement compatibility:


(21)

In the above equation, the total number of equations is N + 7, where N represents the modal equations of the road subsystem and generally takes the low-order vibration mode. The dynamic equilibrium equation will be solved by the Newmark-Beta method. The vibration equation of the standard vehicle–road system can be obtained by rearranging Equation (21): 

$$[M]\left\{\bar{Y}\right\} + [\mathcal{C}]\left\{\dot{Y}\right\} + [\mathcal{K}]\{Y\} = \{F\} \tag{22}$$

where [*M*] represents the total mass matrix, [*C*] represents the total damping matrix, [*K*] represents the total stiffness matrix, and [*F*] represents the total load matrix. The matrix expressions are as follows:

*C*<sup>1</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2*ξ*1*ω*1 + 4 ∑ *j*=1 Φ1 # *xj* \$ *Ctj*<sup>Φ</sup>1 # *xj* \$ 4 ∑ *j*=1 Φ1 # *xj* \$ *Ctj*<sup>Φ</sup>2 # *xj* \$ ··· 4 ∑ *j*=1 Φ1 # *xj* \$ *Ctj*<sup>Φ</sup>*n* # *xj* \$ 4 ∑ *j*=1 Φ2 # *xj* \$ *Ctj*<sup>Φ</sup>1 # *xj* \$ 2*ξ*2*ω*2 + 4 ∑ *j*=1 Φ2 # *xj* \$ *Ctj*<sup>Φ</sup>2 # *xj* \$ ··· 4 ∑ *j*=1 Φ2 # *xj* \$ *Ctj*<sup>Φ</sup>*n* # *xj* \$ . . . . . . ... . . . 4 ∑ *j*=1 Φ*n* # *xj* \$ *Ctj*<sup>Φ</sup>1 # *xj* \$ 4 ∑ *j*=1 Φ*n* # *xj* \$ *Ctj*<sup>Φ</sup>2 # *xj* \$ ··· 2*ξnωn* + 4 ∑ *j*=1 Φ*n* # *xj* \$ *Ctj*<sup>Φ</sup>*n* # *xj* \$ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *C*<sup>2</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ 000 − *Ct*1Φ1(*<sup>x</sup>*1) − *Ct*2Φ1(*<sup>x</sup>*2) − *Ct*3Φ1(*<sup>x</sup>*3) − *Ct*4Φ1(*<sup>x</sup>*4) 000 − *Ct*1Φ2(*<sup>x</sup>*1) − *Ct*2Φ2(*<sup>x</sup>*2) − *Ct*3Φ2(*<sup>x</sup>*3) − *Ct*4Φ2(*<sup>x</sup>*4) . . . . . . . . . . . . . . . . . . . . . 000 − *Ct*1Φ*n*(*<sup>x</sup>*1) − *Ct*2Φ*n*(*<sup>x</sup>*2) − *Ct*3Φ*n*(*<sup>x</sup>*3) − *Ct*4Φ*n*(*<sup>x</sup>*4) ⎤ ⎥ ⎥ ⎥ ⎦ , *C*<sup>3</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 000 0 000 0 000 − *Ct*1Φ1(*<sup>x</sup>*1) − *Ct*1Φ2(*<sup>x</sup>*1) ··· − *Ct*1Φ*n*(*<sup>x</sup>*1) − *Ct*2Φ1(*<sup>x</sup>*2) − *Ct*2Φ2(*<sup>x</sup>*2) ··· − *Ct*2Φ*n*(*<sup>x</sup>*2) − *Ct*3Φ1(*<sup>x</sup>*3) − *Ct*3Φ2(*<sup>x</sup>*3) ··· − *Ct*3Φ*n*(*<sup>x</sup>*3) − *Ct*4Φ1(*<sup>x</sup>*4) − *Ct*4Φ2(*<sup>x</sup>*4) ··· − *Ct*4Φ*n*(*<sup>x</sup>*4) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *C*<sup>4</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *c*11 *c*12 *c*13 −*cs*1 −*cs*2 −*cs*3 −*cs*4 *c*21 *c*22 *c*23 *λ*1*cs*1 − *λ*2*cs*2 *λ*1*cs*3 − *λ*2*cs*4 *c*31 *c*32 *c*33 *λ*3*cs*1 *λ*3*cs*2 − *λ*4*cs*3 − *λ*4*cs*4 −*cs*1 *λ*1*cs*1 *λ*3*cs*1 *cs*1 + *ct*1 000 −*cs*2 − *λ*2*cs*2 *λ*3*cs*2 0 *cs*2 + *ct*2 0 0 −*cs*3 *λ*1*cs*3 − *λ*4*cs*3 0 0 *cs*3 + *ct*3 0 −*cs*4 − *λ*2*cs*4 − *λ*4*cs*4 000 *cs*4 + *ct*4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *K*<sup>1</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *ω*<sup>2</sup> 1 + 4 ∑ *j*=1 Φ1 # *xj* \$ *Ktj*<sup>Φ</sup>1 # *xj* \$ 4 ∑ *j*=1 Φ1 # *xj* \$ *Ktj*<sup>Φ</sup>2 # *xj* \$ ··· 4 ∑ *j*=1 Φ1 # *xj* \$ *Ktj*<sup>Φ</sup>*n* # *xj* \$ 4 ∑ *j*=1 Φ2 # *xj* \$ *Ktj*<sup>Φ</sup>1 # *xj* \$ *ω*<sup>2</sup> 2 + 4 ∑ *j*=1 Φ2 # *xj* \$ *Ktj*<sup>Φ</sup>2 # *xj* \$ ··· 4 ∑ *j*=1 Φ2 # *xj* \$ *Ktj*<sup>Φ</sup>*n* # *xj* \$ . . . . . . ... . . . 4 ∑ *j*=1 Φ*n* # *xj* \$ *Ktj*<sup>Φ</sup>1 # *xj* \$ 4 ∑ *j*=1 Φ*n* # *xj* \$ *Ktj*<sup>Φ</sup>2 # *xj* \$ ··· *ω*<sup>2</sup> *n* + 4 ∑ *j*=1 Φ*n* # *xj* \$ *Ktj*<sup>Φ</sup>*n* # *xj* \$ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *K*<sup>2</sup> = ⎡ ⎢ ⎢ ⎢ ⎣ 000 − *Kt*1Φ1(*<sup>x</sup>*1) − *Kt*2Φ1(*<sup>x</sup>*2) − *Kt*3Φ1(*<sup>x</sup>*3) − *Kt*4Φ1(*<sup>x</sup>*4) 000 − *Kt*1Φ2(*<sup>x</sup>*1) − *Kt*2Φ2(*<sup>x</sup>*2) − *Kt*3Φ2(*<sup>x</sup>*3) − *Kt*4Φ2(*<sup>x</sup>*4) . . . . . . . . . . . . . . . . . . . . . 000 − *Kt*1Φ*n*(*<sup>x</sup>*1) − *Kt*2Φ*n*(*<sup>x</sup>*2) − *Kt*3Φ*n*(*<sup>x</sup>*3) − *Kt*4Φ*n*(*<sup>x</sup>*4) ⎤ ⎥ ⎥ ⎥ ⎦ ,

*K*<sup>3</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 000 0 000 0 000 − *Kt*1Φ1(*<sup>x</sup>*1) − *Kt*1Φ2(*<sup>x</sup>*1) ··· − *Kt*1Φ*n*(*<sup>x</sup>*1) − *Kt*2Φ1(*<sup>x</sup>*2) − *Kt*2Φ2(*<sup>x</sup>*2) ··· − *Kt*2Φ*n*(*<sup>x</sup>*2) − *Kt*3Φ1(*<sup>x</sup>*3) − *Kt*3Φ2(*<sup>x</sup>*3) ··· − *Kt*3Φ*n*(*<sup>x</sup>*3) − *Kt*4Φ1(*<sup>x</sup>*4) − *Kt*4Φ2(*<sup>x</sup>*4) ··· − *Kt*4Φ*n*(*<sup>x</sup>*4) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *K*<sup>4</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *k*11 *k*12 *k*13 −*ks*<sup>1</sup> −*ks*<sup>2</sup> −*ks*<sup>3</sup> −*ks*<sup>4</sup> *k*21 *k*22 *k*23 *λ*1*ks*1 − *λ*2*ks*2 *λ*1*ks*3 − *λ*2*ks*4 *k*31 *k*32 *k*33 *λ*3*ks*1 *λ*3*ks*2 − *λ*4*ks*3 − *λ*4*ks*4 −*ks*<sup>1</sup> *λ*1*ks*1 *λ*3*ks*1 *ks*1 + *kt*1 000 −*ks*<sup>2</sup> − *λ*2*ks*2 *λ*3*ks*2 0 *ks*2 + *kt*2 0 0 −*ks*<sup>3</sup> *λ*1*ks*3 − *λ*4*ks*3 0 0 *ks*3 + *kt*3 0 −*ks*<sup>4</sup> − *λ*2*ks*4 − *λ*4*ks*4 000 *ks*4 + *kt*4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , {*F*} = ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ 4 ∑ *j*=1 Φ1 # *xij* \$% *λj λ ms* + *mtj* + *mlj g* − *Ktjzp* # *xj* \$ − *Ctj* . *zp* # *xj* \$ − *mlj* .. *Zp* # *xj* \$& 4 ∑ *j*=1 Φ2 # *xij* \$% *λj λ ms* + *mtj* + *mlj g* − *Ktjzp* # *xj* \$ − *Ctj* . *zp* # *xj* \$ − *mlj* .. *Zp* # *xj* \$& . . . 4 ∑ *j*=1 Φ*n* # *xij* \$% *λj λ ms* + *mtj* + *mlj g* − *Ktjzp* # *xj* \$ − *Ctj* . *zp* # *xj* \$ − *mlj* .. *Zp* # *xj* \$& 0 0 0 *Kt*1*zp*(*<sup>x</sup>*1) + *Ct*1 . *zp*(*<sup>x</sup>*1) *Kt*2*zp*(*<sup>x</sup>*1) + *Ct*2 . *zp*(*<sup>x</sup>*1) *Kt*3*zp*(*<sup>x</sup>*1) + *Ct*3 . *zp*(*<sup>x</sup>*1) *Kt*4*zp*(*<sup>x</sup>*1) + *Ct*4 . *zp*(*<sup>x</sup>*1) ⎫ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ {*Y*} = ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *A*1 *A*2 . . . *An Zs ψ θ Zt*1 *Zt*2 *Zt*3 *Zt*4 ⎫ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ , ! . *Y* " = ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ . *A*1 . *A*2 . . . . *An*. *Zs* . *ψ* . *θ* . *Zt*1 . *Zt*2 . *Zt*3 . *Zt*4 ⎫ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ , ! .. *Y* " = ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ .. *A*1 .. *A*2 . . . .. *An* .. *Zs* .. *ψ* .. *θ* .. *Zt*1 .. *Zt*2 .. *Zt*3 .. *Zt*4 ⎫ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭ [*C*] = *C*<sup>1</sup> *C*<sup>2</sup> *C*<sup>3</sup> *C*<sup>4</sup> , [*K*] = *K*<sup>1</sup> *K*<sup>2</sup> *K*<sup>3</sup> *K*<sup>4</sup> 

,

It is difficult to see from the above matrices that the total mass matrix [*M*], the total damping matrix [*C*], and the total stiffness matrix [*K*] change with the road position *x* (*<sup>x</sup>*1, *x*2, *x*3, and *x*4 represent the positions of the four wheels). Therefore, the vehicle road system coupling dynamic balance equations can be regarded as a highly second-order nonlinear differential equation. By solving the equation, the vertical vibration displacement *Zs*, vertical velocity . *Zs*, and vertical acceleration .. *Zs* of the vehicle body; the vertical vibration displacement *Zt*1, *Zt*2, *Zt*3, and *Zt*4 of each unsprung mass part of the vehicle body, the vertical velocity . *Zt*1, . *Zt*2, . *Zt*3, and . *Zt*4; and the vertical acceleration .. *Zt*1, .. *Zt*2, .. *Zt*3 and .. *Zt*4 can be obtained. At the same time, the nodal displacement *θ*, nodal velocity . *θ*, nodal acceleration .. *θ*, and deflection displacement *ψ* of the vehicle body can also be obtained using the deflection velocity . *ψ* and deflection acceleration .. *ψ*.

The vehicle is regarded as a multi-stiffness system. According to D'Alembert's principle, considering the dynamic balance conditions of each rigid body in the system, the road vibration balance equation and the vehicle dynamic balance equation are listed according to the coupling model under the contact force condition in Figure 8.

**Figure 8.** Coupling vibration relation model under contact force.

The dynamic balance equation of the vehicle–road system is obtained by combining Equation (20). The total number of equations is N + 7, where N represents the number of modal equations of the path sub-model.

$$[M]\left\{\bar{Y}\right\} + [\mathbb{C}]\left\{\dot{Y}\right\} + [K]\{Y\} = \{F\} \tag{23}$$

The Newmark-β and Park integral methods are used to solve the above equations, and the responses of the vehicle and road in the time domain, such as the force, displacement, and vibration velocity, can be obtained [27].

### **4. Theoretical Analysis of Vehicle–Road System**

According to the vehicle–road coupling dynamic equation and the dynamic balance equation listed above, the following vehicle driving mass and road design parameters can be obtained using the Newmark-β method in a self-programmed program. It is assumed that the heavy vehicle has a mass of 20 t and runs at a speed of 80 km/h on a class B road. The road structure adopts planar 8-node PLANE82 unit, and the road model takes 100 m along the road longitudes. The calculated step time is 0.02. See Appendices A and B for detailed parameters.
