Investigation on dynamics of a three-directional coupled vehicle-road system

1. Introduction

With the increase in road traffic and vehicle loads, premature damage of the asphalt pavement on expressways has become more prevalent, which has been greatly reducing the pavement’s effective lifetime. In China, some road damage occurs within the first 2 or 3 years after the road use [1]. Road damage is seen to cause an increase in dynamic tire forces, which may worsen vehicle vibration and reduce the passenger’s ride comfort and safety. Accordingly, the increased vehicle vibration leads to an increased tire force, which may speed up road damage. Hence, the vehicle and road together form a coupled system and the dynamics of vehicle and road should be investigated simultaneously based on a vehicle-road coupled system.

However, vehicle dynamics and road dynamics are two separate subjects at present. In these two subjects, the vehicle-road interaction is investigated using different models and methods, and a lot of evaluation index of road damage and design rules of vehicle and road parameters have been proposed and concluded.

In vehicle dynamics, road surface roughness is generally regarded as excitation to vehicle, and the dynamic tire loads and the road-friendly characteristics of vehicles are investigated according to road damage evaluation criterions. The popular vehicle models used for researching road-friendly characteristics includes a quarter-vehicle model with 2-DOF, a half-vehicle model with 4- or 5-DOF, a full-vehicle model with many degrees of freedom (DOF) and a Functional Virtual Prototyping (FVP) model. The dynamic load coefficient (DLC) is one of the most widely adopted evaluation criterion for dynamic interaction [2, 3]. Cole and Cebon put forward a fourth power aggregate force to evaluate the road damage induced by tire dynamic forces, which is another popular criterion for road damage. They investigated the effect of heavy vehicle parameters on tire pressure and road damage and optimized the suspension [3, 4]. Markow and Myers investigated the influence of heavy vehicle characteristics on tire forces and pavement response. They found that the vehicle configuration and suspension system greatly influences tire forces and that the tire pressure has little effect on the dynamic load [5, 6]. Sun optimized a 3-DOF quarter-truck suspension system by minimizing the dynamic pavement load and found that a large tire pressure and a small suspension damping may increase dynamic tire forces [7]. Lu validated a FVP full-vehicle model with experimental data of a heavy vehicle and discussed the effects of vehicle speed, vehicle load and road surface roughness on tire forces and DLC [8]. Chen simulated the influences of semi-active MSD and PID control on a heavy truck’s road friendliness and ride comfort through co-simulations with Adams and MATLAB [9]. It should be noted that these investigations in the field of vehicle dynamics generally assumed that the pavement was rigid and stationary and did not take into account the effect of pavement vibration on vehicle response.

In road dynamics, vehicles are generally regarded as moving loads acting on the pavement, and roads are modeled as a beam, a plate, or a multi-layer system on an elastic or viscoelastic foundation. The road response such as stress, strain and displacements of the pavement and foundation are studied analytically or numerically. And the road fatigue lifetime is predicted using different evaluation criterions, including the tensile strain at the bottom of an asphalt surface and the compressive strain in soil foundation. The Strategic Highway Research Program (SHRP) in the U.S. performed comprehensive research on road dynamics and road damage potential of dynamic wheel loads [10-11]. Hardy used the mode-superposition theory and the integral transform to derive an analytical solution of the displacement of a beam on a viscoelastic foundation, and validated it by field experiments [12]. Deng investigated the dynamics of a beam on a Winkler or a Kelvin foundation under moving vehicle loads using an integral transform [13]. Using the mode superposition and the integral transform methods, Kim studied the response of a plate on a Winkler foundation under moving loads of varying amplitude. He analyzed the influence of wheel space and driving speed on the plate displacement and stress [14]. Jeongho modeled the base and subgrade layers as stress-dependent cross-anisotropic materials and assess the pavement response using finite-element (FE) analysis. He established equations correlating the critical strains to layer displacements, axle loading, offset distance, and layer moduli in order to evaluate the accelerated damage potential due to overweight truck loading [15]. Hryniewicz presented a new analytical solution for the dynamic response of an infinitely long Timoshenko beam resting on a nonlinear viscoelastic foundation using the wavelet expansion combined with Adomian’s decomposition [16]. Ding investigated the dynamic response of infinite Timoshenko beams supported by nonlinear viscoelastic foundations subjected to a moving concentrated force using the Adomian decomposition and a perturbation method in conjunction with complex Fourier transformation [17]. Fang studied the dynamic response of a thin-plate resting on a layered poroelastic half-space under a moving traffic load using Fourier transformation [18]. While the researches in road dynamics are useful for predicting road stresses and deformations, they only concentrate on modeling the road, but seldom model the vehicle suspension system. Consequently, the road dynamics due to vehicle vibration cannot be predicted accurately.

Over the last few decades, some scholars began to research vehicle-road dynamic interaction simultaneously. They connected vehicle and road model by tire forces and proposed different vehicle-road coupled models. The responses of vehicle and road are calculated using Finite Element (FE), Galerin methods, mode superposition, integral transformation and numerical integral. Using two sets of models to simulate vehicle and road respectively, Markov summarized that the key factors influencing a stiff pavement were the configuration of the vehicle and axle, vehicle load, suspension properties, vehicle speed, road roughness and plate bending [19]. Collop proposed a whole-life pavement-performance model to predict pavement damage due to tire forces. A 2-DOF vehicle model was used to generate dynamic tire forces and two classes of pavement were researched [20]. Wu considered the dynamic vehicle-pavement-foundation interaction effects in the 3D finite-element algorithm. The moving vehicle loads are modeled as lumped masses each supported by a spring-dashpot suspension system and the effects of a finite element grid, pavement thickness and foundation stiffness on pavement displacement were evaluated [21]. Mamlouk pointed out that vehicle-road interaction can be applied to weigh-in-motion, pavement design, and vehicle regulation [22]. Considering the dynamic variation and complex distribution of the tire-pavement contact stress, together with vehicle speed and viscoelastic material properties, Siddharthan developed a 3D-MOVE procedure using the Fourier-Transformation method and simulated pavement responses to the moving tire loads [23]. Liu built a 2-DOF vehicle and a plate on a viscoelastic foundation, and derived the analytical solution of the pavement displacement under moving vehicle loads using an integral transform. He analyzed the effects of road roughness and foundation modulus on pavement displacement and strains, and found that a large foundation modulus may lead to a small strain [24]. Giuseppe studied the response of a beam on a viscoelastic foundation subjected to a single-DOF oscillator. He converted the coupled-system equations to dimensionless ones and calculated the response by the methods of mode superposition and numerical integration [25]. Metrikine studied the stability of a two-mass oscillator that moves along a beam on a viscoelastic half-space. Using the Laplace and the Fourier integral transforms, he derived expressions for the dynamic stiffness of the beam at the point of contact with the oscillator and concluded that a proper combination of negative damping mechanisms may effectively stabilize the system [26]. Shi formulated a three-dimensional vehicle-pavement coupled model to simulate the pavement dynamic loads induced by the vehicle-pavement interaction and analyzed the effects of road surface conditions, vehicle parameters, and driving speed on pavement dynamic loads [27]. Taheri proposed a mechanistic-empirical pavement damage model to predict changes in 3D road profiles due to dynamic axle loads [28]. Ding researched a coupled nonlinear vibration of vehicle-pavement system, which is composed of a Timoshenko beam resting on a six-parameter foundation and a spring-mass-damper oscillator [29]. In the previous work, the authors put forward a vertical coupled vehicle-road model and investigated the effects of vehicle-road interaction on vehicle and road responses [30-32]. However, the present studies on vehicle-road coupled systems only consider the vertical tire load but usually neglect the impact of lateral and longitudinal tire loads on vehicle-road interaction.

This work tries to connect the vehicle and road by vertical, lateral and longitudinal tire forces, and put forward the methodology of three-directional (3D) coupled vehicle-road system modeling and simulating. Based on a 3D vehicle-road coupled model, the differences of 3D coupled, vertical coupled and uncoupled model are compared for pavement displacements and vehicle responses. The effects of vehicle maneuver parameters and road surface properties on 3D vehicle-road interaction are also analyzed. The proposed 3D-coupled vehicle-road model and the simulation method are able to research the vehicle-road interaction and will be beneficial to vehicle and road dynamics study.

2. Modeling the three-directional coupled vehicle-road system

A three-directional (3D) coupled vehicle-road system is proposed in this work as shown in Fig. 1. A 23-DOF vehicle and a double-layer rectangular thin plate on viscoelastic foundation with four simply supported boundaries are employed to model vehicle and pavement. The upper layer of the plate models asphalt topping, the lower layer models base course, and the viscoelastic foundation stands for subgrade of the road. The material of topping and base course is assumed to be isotropy and elastic. The vehicle is able to moving on the road in different maneuver conditions, such as straight-line running at a constant speed, braking and turning. The origin of vehicle coordinate system ($x$, $y$, $z$) is located in the intersection between the roll axle and the vertical line passing vehicle center of gravity. $\psi$ is the vehicle’s yaw angle which describe the orientation of the vehicle. The static ground coordinate system ($X$, $Y$, $Z$) is also shown in Fig. 1, taking the stress neutral layer of pavement as coordinate axial $X$, the lateral direction as coordinate axial $Y$ and the vertical direction as coordinate axial $Z$.

Fig. 1The 3D coupled vehicle-road model

2.1. Equations of motion for the vehicle

A nonlinear full-body model with 23-DOF for a three-axial heavy-duty truck is presented, as shown in Fig. 2. The vehicle is front wheel steered and rear wheel driven. $x$, $y$, $\psi$ represent the longitudinal and lateral displacements and the yaw angle of the full vehicle respectively. $z_c$, ${\theta }_c$, ${\varphi }_c$, $z_b$, $\theta$, $\varphi$ stand for the vertical, pitch and roll displacements of driver cab and vehicle body. $z_{ui}$, ${\varphi }_{ui}$ ($i=$1,…, 3) denote the vertical and roll displacements of three wheel axles. ${\theta }_{p1}$, ${\theta }_{p2}$ are the pitch angles of the left and right balancing pole on rear suspension. $M$, $M_b$, $M_c$ denote the masses of full vehicle, vehicle body and driver cab respectively. $h_s$, $l_s$ are the vertical and longitudinal distance from the sprung mass center of gravity to the coordinate origin. $d_{t1}$, $d_{t2}$, $d_{t3}$ are the front and rear wheel track width. $d_{s1}$, $d_{s2}$, $d_{s3}$ are the lateral distance between left and right spring on front or rear suspension.

The movements of the heavy-duty vehicle are coupled with each other greatly. Eqs. (1)-(3) give the longitudinal, lateral and yaw dynamics of the full vehicle. Eqs. (4)-(6) give the vertical, roll and pitch dynamics of the sprung mass:

1

$M(\ddot{x}-\dot{y}\dot{\psi })+M_b{\dot{z}}_b\dot{\theta }+M_b{h}_s\ddot{\theta }-M_b{l}_s({\dot{\psi }}^2+{\dot{\theta }}^2)=(F_{tx11}+F_{tx12})\mathrm{c}\mathrm{o}\mathrm{s}\delta$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}-\left(F_{ty11}+F_{ty12}\right)\mathrm{s}\mathrm{i}\mathrm{n}\delta +\sum _{i=2}^3\sum _{j=1}^2{F}_{txij},$

2

$M(\ddot{y}+\dot{x}\dot{\psi })-M_b{\dot{z}}_b\dot{\phi }-M_b{h}_s\ddot{\phi }=(F_{tx11}+F_{tx12})\mathrm{s}\mathrm{i}\mathrm{n}\delta$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\left(F_{ty11}+F_{ty12}\right)\mathrm{c}\mathrm{o}\mathrm{s}\delta +\sum _{i=2}^3\sum _{j=1}^2{F}_{tyij},$

3

$I_z\ddot{\psi }+2I_{zb}\dot{\phi }\dot{\theta }-(I_{bxz}+M_b{l}_s{h}_s)\ddot{\phi }-M_b{l}_s(\ddot{y}+\dot{x}\dot{\psi }-{\dot{z}}_b\dot{\phi })=l_1\left[\right(F_{tx11}+F_{tx12})\mathrm{s}\mathrm{i}\mathrm{n}\delta$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+(F_{ty11}+F_{ty12})\mathrm{c}\mathrm{o}\mathrm{s}\delta ]-\left(l_2-\frac{l_3}{2}\right)\mathrm{}\sum _{j=1}^2{F}_{ty2j}-\left(l_2+\frac{l_3}{2}\right)\mathrm{}\sum _{j=1}^2{F}_{ty3j}$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\frac{d_{t1}}{2}\left[\right(-F_{tx11}+F_{tx12})\mathrm{c}\mathrm{o}\mathrm{s}\delta +(F_{ty11}-F_{ty12}\left)\mathrm{s}\mathrm{i}\mathrm{n}\delta \right]$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\frac{d_{t2}}{2}\left(-F_{ty21}+F_{ty22}-F_{ty21}+F_{ty22}\right)+\sum _{i=1}^3\sum _{j=1}^2{T}_{zij},$

4

$M_b({\ddot{z}}_b-\dot{x}\dot{\theta }+\dot{y}\dot{\phi })-M_b{h}_s({\dot{\theta }}^2+{\dot{\phi }}^2)-M_b{l}_s\ddot{\theta }-(F_{c1}+F_{c2}+F_{c3}+F_{c4})$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\left(F_{s11}+F_{s12}+F_{s21}+F_{s22}+F_{s31}+F_{s32}\right)=-M_{b}g,$

5

$(M_b{h_s}^2+I_{bx})\ddot{\phi }-(I_{bxz}+M_b{l}_s{h}_s)(\ddot{\psi }+2\dot{\phi }\dot{\theta })-M_b{h}_s(\ddot{y}+\dot{x}\dot{\psi }-{\dot{z}}_b\dot{\phi })$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+(F_{c2}+F_{c4}-F_{c1}-F_{c3})\frac{d_c}{2}+(F_{s11}-F_{s12})\frac{d_{s1}}{2}+(F_{s21}-F_{s22})\frac{d_{s2}}{2}$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+(F_{s31}-F_{s32})\frac{d_{s3}}{2}=(z_b-z_{t11}+R_{11})(F_{tx11}\mathrm{s}\mathrm{i}\mathrm{n}\delta +F_{ty11}\mathrm{c}\mathrm{o}\mathrm{s}\delta )$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+(z_b-z_{t12}+R_{12})(F_{tx12}\mathrm{s}\mathrm{i}\mathrm{n}\delta +F_{ty12}\mathrm{c}\mathrm{o}\mathrm{s}\delta )+\sum _{i=2}^3\sum _{j=1}^2{F}_{tyij}\left(z_b-z_{tij}+R_{ij}\right),$

6

$(M_b{h}_s^2+I_{by})\ddot{\theta }-M_b{h}_s^{}{l}_s({\dot{\psi }}^2+{\dot{\theta }}^2)+M_b{h}_s(\ddot{x}-\dot{y}\dot{\psi }+{\dot{z}}_b\dot{\theta })-M_b{l}_s({\ddot{z}}_b+g-\dot{x}\dot{\theta }+\dot{y}\dot{\phi })$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+M_b{l}_{{\mathrm{}}_s}^2\ddot{\theta }+(F_{c1}+F_{c2})(l_4+l_5)+(F_{c3}+F_{c4})(l_4-l_6)-(F_{s11}+F_{s12})l_1$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+(F_{s21}+F_{s22}+F_{s31}+F_{s32})l_2=(z_b-z_{t11}+R_{11})(F_{tx11}\mathrm{c}\mathrm{o}\mathrm{s}\delta -F_{ty11}\mathrm{s}\mathrm{i}\mathrm{n}\delta )$
$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\left(z_b-z_{t12}+R_{12}\right)\left(F_{tx12}\mathrm{c}\mathrm{o}\mathrm{s}\delta -F_{ty12}\mathrm{s}\mathrm{i}\mathrm{n}\delta \right)+\sum _{i=2}^3\sum _{j=1}^2{F}_{txij}\left(z_b-z_{tij}+R_{ij}\right),$

where $\delta$, $T_{zij}$ and $R_{ij}$ are the steering angle of front wheel, the tire aligning torque, and the tire effective radius respectively. $F_{cs}$ ($s=$1,…, 4) denotes the suspension forces between driver cab and vehicle body and can be expressed as:

7

$\left\{\begin{array}{l}{F}_{c1}=K_{c1}\left(z_c-{\theta }_c{l}_5-z_b+\left(\theta -{\theta }_0\right)\left(l_4+l_5\right)-\frac{\left(\phi -{\phi }_c\right)d_c}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{c1}\left({\dot{z}}_c-{\dot{\theta }}_c{l}_5-{\dot{z}}_b+\dot{\theta }\left(l_4+l_5\right)-\frac{\left(\dot{\phi }-{\dot{\phi }}_c\right)d_c}{2}\right),\\ F_{c2}=K_{c2}\left(z_c-{\theta }_c{l}_5-z_b+\left(\theta -{\theta }_0\right)\left(l_4+l_5\right)+\frac{\left(\phi -{\phi }_c\right)d_c}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{c2}\left({\dot{z}}_c-{\dot{\theta }}_c{l}_5-{\dot{z}}_b+\dot{\theta }\left(l_4+l_5\right)+\frac{\left(\dot{\phi }-{\dot{\phi }}_c\right)d_c}{2}\right),\\ F_{c3}=K_{c3}\left(z_c+{\theta }_c{l}_6-z_b+\left(\theta -{\theta }_0\right)\left(l_4-l_6\right)-\frac{\left(\phi -{\phi }_c\right)d_c}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{c3}\left({\dot{z}}_c+{\dot{\theta }}_c{l}_6-{\dot{z}}_b+\dot{\theta }\left(l_4-l_6\right)-\frac{\left(\dot{\phi }-{\dot{\phi }}_c\right)d_c}{2}\right),\\ F_{c4}=K_{c4}\left(z_c+{\theta }_c{l}_6-z_b+\left(\theta -{\theta }_0\right)\left(l_4-l_6\right)+\frac{\left(\phi -{\phi }_c\right)d_c}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{c4}\left({\dot{z}}_c+{\dot{\theta }}_c{l}_6-{\dot{z}}_b+\dot{\theta }\left(l_4-l_6\right)+\frac{\left(\dot{\phi }-{\dot{\phi }}_c\right)d_c}{2}\right).\end{array}\right.$

For this heavy-duty vehicle, two hydraulic dampers are fixed on the left and right front suspension and the tandem balanced suspension doesn’t have any damper. The damper in vehicle suspension is one of the major factors influencing the vehicle handling and ride comfort and shows apparent nonlinearity, asymmetry and hysteresis. Many dynamic models for vehicle suspension damper have been proposed, including the parametric model [33], the equivalent parametric model [34, 35] and the fitted model [36, 37]. The fitted model is quite suitable to modeling the ascertained damper and the function of restoring force relate to velocity are fitted by the test data. In this work, the fitted model is used. The damping force of the damper on front axle can be described by a nonlinear segmented model:

8

$F_d=\left\{\begin{array}{ll}{F}_{01}\mathrm{s}\mathrm{g}\mathrm{n}\left(v_d\right),& v_d>v_{\mathrm{l}\mathrm{i}\mathrm{m}1},\\ C\left(1+\beta \mathrm{s}\mathrm{g}\mathrm{n}\left(v_d\right)\right)v_d{\left|v_d\right|}^n,& v_{\mathrm{l}\mathrm{i}\mathrm{m}2}\le v_d\le v_{\mathrm{l}\mathrm{i}\mathrm{m}1},\\ F_{02}\mathrm{s}\mathrm{g}\mathrm{n}\left(v_d\right),& v_d<v_{\mathrm{l}\mathrm{i}\mathrm{m}2},\end{array}\right.$

where $v_d$, $C$, $\beta$ and $n$ are the relative velocity of cylinder and plunger, the damping coefficient, the asymmetry ratio and the exponent respectively. $F_{01}$, $F_{02}$, $v_{\mathrm{l}\mathrm{i}\mathrm{m}1}$ and $v_{\mathrm{l}\mathrm{i}\mathrm{m}2}$ are the damping force and relative velocity when the damper reaching saturation in tension or compression process. Eq. (8) is composed of three equations and modeled from the test data [38]. The first and third equations are dry friction damping models, which define the damping forces when the damper plunger arrives at the above or below limit position. The second equation describes the hysteresis viscous damping force when the plunger moves between the two limit positions. Here, $C=$30893, $\beta =\mathrm{}$0.56, $n=$0.16, $F_{01}=$4119 N, $F_{02}=$726 N, $v_{\mathrm{l}\mathrm{i}\mathrm{m}1}=$0.12 m/s, $v_{\mathrm{l}\mathrm{i}\mathrm{m}2}=$0.08 m/s. These parameters are obtained by fitting the measured data of the damper’s dynamic test. The test’s description and the related results were given by the author’s another work [38].

Fig. 2The full-body heavy vehicle model with 23-DOF

Using model Eq. (8) to calculate the damping force, the front suspension forces are expressed by:

9

$F_{s11}=K_{s11}\left(z_b-\left(\theta -{\theta }_0\right)l_1-z_{u1}+\frac{\left(\phi -{\phi }_{u1}\right)d_{s1}}{2}\right)+F_{d11},$
$F_{s12}=K_{s12}\left(z_b-\left(\theta -{\theta }_0\right)l_1-z_{u1}-\frac{\left(\phi -{\phi }_{u1}\right)d_{s1}}{2}\right)+F_{d12},$

where $F_{d11}$, $F_{d12}$ are the left and right damping force of front suspension. The relative velocities of left and right damper are obtained by $v_{d1}=({\dot{z}}_b-\dot{\theta }{l}_1-{\dot{z}}_{u1}+(\dot{\phi }-{\dot{\phi }}_{u1})d_{s1}/2)$ and $v_{d2}=({\dot{z}}_b-\dot{\theta }{l}_1-{\dot{z}}_{u1}-(\dot{\phi }-{\dot{\phi }}_{u1})d_{s1}/2)$.

In order to represent the frictional property of leaf spring, the damping forces of tandem balanced suspension are modeled linearly. The suspension forces between middle or rear axle and vehicle body are expressed by:

10

$\left\{\begin{array}{l}{F}_{s21}=K_{s21}\left(z_b+\left(\theta -{\theta }_0\right)l_2-\frac{{\theta }_{p1}{l}_3}{2}-z_{u2}+\frac{\left(\phi -{\phi }_{u2}\right)d_{s2}}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{s21}\left({\dot{z}}_b+\dot{\theta }{l}_2-\frac{{\dot{\theta }}_{p1}{l}_3}{2}-{\dot{z}}_{u2}+\frac{\left(\dot{\phi }-{\dot{\phi }}_{u2}\right)d_{s2}}{2}\right),\\ F_{s22}=K_{s22}\left(z_b+\left(\theta -{\theta }_0\right)l_2-\frac{{\theta }_{p2}{l}_3}{2}-z_{u2}-\frac{\left(\phi -{\phi }_{u2}\right)d_{s2}}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{s22}\left({\dot{z}}_b+\dot{\theta }{l}_2-\frac{{\dot{\theta }}_{p2}{l}_3}{2}-{\dot{z}}_{u2}-\frac{\left(\dot{\phi }-{\dot{\phi }}_{u2}\right)d_{s2}}{2}\right),\\ F_{s31}=K_{s31}\left(z_b+\left(\theta -{\theta }_0\right)l_2+\frac{{\theta }_{p1}{l}_3}{2}-z_{u3}+\frac{\left(\phi -{\phi }_{u3}\right)d_{s3}}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{s31}\left({\dot{z}}_b+\dot{\theta }{l}_2+\frac{{\dot{\theta }}_{p1}{l}_3}{2}-{\dot{z}}_{u3}+\frac{\left(\dot{\phi }-{\dot{\phi }}_{u3}\right)d_{s3}}{2}\right),\\ F_{s32}=K_{s32}\left(z_b+\left(\theta -{\theta }_0\right)l_2+\frac{{\theta }_{p2}{l}_3}{2}-z_{u3}-\frac{\left(\phi -{\phi }_{u3}\right)d_{s3}}{2}\right)\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+C_{s32}\left({\dot{z}}_b+\dot{\theta }{l}_2+\frac{{\dot{\theta }}_{p2}{l}_3}{2}-{\dot{z}}_{u3}-\frac{\left(\dot{\phi }-{\dot{\phi }}_{u3}\right)d_{s3}}{2}\right).\end{array}\right.$

Eq. (11) gives the pitch motions of the balancing poles:

11

$I_{pi}{\ddot{\theta }}_{pi}+(F_{s3i}-F_{s2i})\frac{l_3}{2}=0.$

Here, the subscript $i=$1, 2 stands for the left and right balancing pole of rear suspension.

The vertical, roll and pitch motions of the cab are formulated by:

12

$\left\{\begin{array}{l}{M}_c{\ddot{z}}_c+\left(F_{c1}+F_{c2}+F_{c3}+F_{c4}\right)=-M_{b}g,\\ I_{cx}{\ddot{\varphi }}_c+(F_{c1}+F_{c3}-F_{c2}-F_{c4})\frac{d_c}{2}=0,\\ I_{cy}{\ddot{\theta }}_c-(F_{c1}+F_{c2})l_5+(F_{c3}+F_{c4})l_6=0.\end{array}\right.$

The vertical and roll motions of three axles and the wheel rotating rates are given by:

13

$\left\{\begin{array}{l}{M}_{ui}{\ddot{z}}_{ui}-F_{si1}-F_{si2}=F_{tzi1}+F_{tzi2}-M_{ui}g,\\ I_{ui}{\ddot{\mathrm{\varphi }}}_{ui}+\left(F_{si2}-F_{si1}\right)\frac{d_{si}}{2}=\left(F_{tzi1}-F_{tzi2}\right)\frac{d_{ti}}{2}+\left(F_{tyi1}+F_{tyi2}\right)R,\end{array}\right.$

14

$I_{ij}{\dot{\omega }}_{ij}=T_{sij}-T_{bij}-R_{ij}\cdot F_{txij},$

where, $T_{sij}$ and $T_{bij}$ ($i=$1,…, 3, $j=$1, 2) are the driving torque and braking torque of six wheels. The subscript $i$ stands for the front, middle or rear axle, and $j$ stands for the left or right wheel.