Abstract

TECHNICAL PAPERS

Vehicle modeling by subsystems

Georg Rill

GeorgRill@aol.com FH Regensburg, University of Applied Sciences Galgenbergstr. 30 D-93053 Regensburg, Germany

ABSTRACT

Computer simulations have become very popular in the automotive industry. In order to achieve a good conformity with field test, sophisticated vehicle models are needed. A real vehicle incorporates many complex dynamic systems, such as the drive train, the steering system and the wheel/axle suspension. On closer inspection some force elements such as shock absorbers and hydro-mounts turn out to be dynamic systems too. Modern vehicle models consist of different subsystems. Then, each subsystem may be modeled differently and can be tested independently. If some subsystems are available as a set of nested models of different complexity it will be even possible to generate overall vehicle models which are well tailored to particular applications. But, the numerical solution of coupled subsystems is not straight forward. This paper shows that the overall vehicle model can be solved very effectively by suitable interfaces and an implicit integration algorithm. The presented concept is realized in the product ve-DYNA, applied worldwide by automotive companies and suppliers.

Keywords: Vehicle dynamics, vehicle model, axle modeling, drive train, multibody systems

Modeling Concept

For dynamic simulation the vehicles are usually modeled by multi body systems (MBS), van der Jagt (2000). Typically, the overall vehicle model is separated into different subsystems, Rauh (2003). Fig. 1 shows the components of a passenger car model which can be used to investigate handling and ride properties. The vehicle model consists of the vehicle framework and subsystems for the steering system and the drive train.

The vehicle framework represents the kernel of the model. It at least includes the module chassis and modules for the wheel/axle suspension systems. The vehicle framework is supplemented by modules for the load, an elastically suspended engine, and passenger/seat models. A simple load module just takes the mass and inertia properties of the load into account. To describe the sloshing effects of liquid loads dynamic load models are needed, Rill and Rauh (1992). The subsystems elastically suspended engine, passenger/seat, and in heavy truck models a suspended driver's cabin can all be handled by the presented generic free body model. For standard vehicle dynamics analysis the chassis can be modeled by one rigid body. For applications where the chassis flexibility has to be taken into account a suitable flexible frame model is presented. Most wheel/axle suspension systems can be described by typical multi body system elements such as rigid bodies, links, joints and force elements, Rill (1994). Using a modified implicit Euler algorithm for solving the dynamic equations, axle suspensions with compliancies and dry friction in the damper element can be handled without any problems, Rill (2004). Due to their robustness leaf springs are still a popular choice for solid axles. They combine guidance and suspension properties which causes many problems in modeling, Fickers and Richter (1994). A leaf spring model is presented in this paper which overcomes these problems.

The steering system at least consists of the steering wheel, a flexible steering shaft, and the steering box which may also be power-assisted. Neureder (2002) has developed a very sophisticated model of the steering system which includes compliancies, dry friction, and clearance.

Tire forces and torques have a dominant influence on vehicle dynamics. The semi-empirical tire model TMeasy has mainly been developed to meet both the requirements of user-friendliness and sufficient model accuracy, Hirschberg et. al. (2002). Complex tire models such as the FTire Model provided by Gipser (1998) can be used for special applications. The module tire also includes the wheel rotation which acts as input for the drive train model. The presented drive train model is generic. It takes lockable differentials into account, and it combines front wheel, rear wheel and all wheel drive. The drive train is supplemented by a module describing the engine torque. It may be modeled quite simply by a first order differential equation or by the enhanced engine torque module en-DYNA developed by TESIS.

Road irregularities and variations in the coefficient of friction present significant impacts on the vehicle. A road model generating a two-dimensional reproducible random profile was provided by Rill (1990).

This modeling concept is realized with a MATLAB/Simulink® interface in the product ve-DYNA which also includes suitable models for the driver, TESIS.

Module Flexible Frame

Multi Body Approach to First Eigenmodes

The chassis eigenmodes of most passenger cars start at f >20Hz. Hence, for standard vehicle dynamic analysis the chassis can be modeled as one rigid body. The lower chassis stiffness of trucks and pickups results in eigenmodes starting at ¦ » 10 Hz, Fig. 2.

The first eigenmodes consist of torsion and bending of the chassis. These modes can be approximated by a multi body chassis model where the chassis is divided into three parts, Fig. 3.

Free Body Motions

The position and orientation of the reference frame x_c, y_c, z_c which is fixed to the center body with respect to the inertial frame x₀, y₀, z₀ is given by the rotation matrix

and the position vector

where the comma separated subscript 0 indicates that the coordinates of the vector from to are expressed in the inertial frame. The generalized coordinates roll, pitch, and yaw angle a_0C, b_0C, g_0C as well as the coordinates x_0C, y_0C, z_0C of the vector r_0C,0 describe the free body motion of the vehicle.

Modal Coordinates

The motions of the front and rear body relative to the center body are small compared to the free body motions of the center part. Hence, the linearized rotation matrices

and the position vectors

are used to describe the orientation and position of the front and rear body relative to the center part. The vectors and denote the initial position of the front and rear body.

The generalized coordinates

describe the motions of the front and rear body relative to the center body. These motions are approximated by n_M eigenmodes e₁, e₂, ... en_M, now

where m₁ , m₂ , ..., m_nM are modal coordinates, and E_F and E_R represent 6´n_M matrices containing the eigenmodes.

Generalized Coordinates

The flexible chassis is modeled by 3 rigid bodies here. The orientation and the position of the bodies are described by free body motions and modal coordinates

where the 6 free body motions and the n_M modal coordinates are collected in the vector y_C. The dimension of depends on the number n_M of modal coordinates, n_y = 6 + n_M.

Equations of Motion

To generate the equations of motion Jordain's Principle with generalized speeds is used. For a multi body system consisting of m rigid bodies it results in a set of two first order differential equations

The kinematic matrix K follows from the definition of the generalized speeds. The elements of the mass matrix M are given by

where m_k is the mass and Q_k the inertia tensor of body k. Finally, the components of the generalized forces and torques are defined by

where F_Ak, T_Ak denote the forces and torques applied to body k and , are remaining parts of the accelerations which do not depend on the derivatives of the generalized speeds.

Applied Forces and Torques

The forces and torques applied to the bodies can be written as

and

where the superscripts ext and cmp denote external and compliance forces and torques.

Applying Jordain's Principle, one part within the equations of motion describes the whole chassis motion. The compliance forces and torques are internal forces for the whole chassis and, therefore, do not show up in the corresponding parts of the generalized force vector.

If we assume that the compliance forces and torques are proportional to the motions of the front and rear body then, we will get

where c_CF and c_CR mark 6´6 stiffness matrices. The modal coordinate approximation Eq. 6 results in

Within Jordain's Principle the compliance forces and torques are reduced to generalized forces which are calculated by

where

are n_M´n_M stiffness matrices. which are defined by the modal stiffnesses c_F1, c_F2, ... cFn_M and c_R1, c_R2, ... cRn_M. Thus, to describe the motions of a flexible chassis only some eigenmodes and modal stiffnesses have to be provided.

Results

Depending on the vehicle layout, a flexible frame has a significant influence on the driving behavior, Fig 4. The rear axle of the considered bus is guided by four links. Here, the arrangement of the links generates a steering effect which depends on the roll angle of the rear part of the chassis and, therefore, also on the torsional stiffness of the chassis.

Module Leaf Spring

Modeling Aspects

Poor leaf spring models approximate guidance and suspension properties of the leaf spring by rigid links and separate force elements, Matschinsky (1998). The deformation of the leaf springs must be taken into account for realistic ride and handling simulations.

Within ADAMS leaf springs can be modeled with sophisticated beam-element models, ADAMS/Chassis 12.0. But, according to Fickers (1994) it is not easy to take the spring pretension into account. To model the effects of a beam, ADAMS/Solver uses a linear 6 -dimensional action-reaction force (3 translational and 3 rotational) between two markers. In order to provide adequate representation for the nonlinear cross section, usually 20 elements are used to model one leaf spring. A subsystem consisting of a solid axle and two beam-element leaf spring models would have ¦ = 6+2 * (20*6) = 246 degrees of freedom. In addition, the beam-element leaf spring model results in extremely stiff differential equations. These and the large number of degrees of freedom slow down the computing time significantly.

For real time applications the leaf springs must be modeled by a simple, but still accurate model. Fig. 5 shows a model of a solid axle with leaf spring suspension, which is typical for light truck rear axle suspension systems. There are no additional links. Hence, only the forces and torques generated by leaf spring deflections guide and suspend the axle.

The position of the axle center A and the orientation of an axle fixed reference frame x_A, y_A, z_A are described relative to a chassis fixed frame x_B, y_B, z_B by the displacements x, h, z and the rotation angles a, b, g which are collected in the 6´1 axle position vector

Similar to Fickers (1994) each leaf spring is modeled by five rigid bodies which are connected to each other by spherical joints, Fig. 5.

Each leaf spring is connected to the frame via the front leaf eye X. Furthermore each leaf spring is attached to the shackle at Y, and again to the frame at Z. In C the center part of each leaf spring is rigidly connected to the axle. The front eye bushings are modeled by spring/damper elements in x-, y-, and z-direction. The shackles are modeled by radial and a lateral spring/damper elements. Within each leaf spring the angles j₁, y1, and j₂, y2, describe the motions of part P-Q and part R-S relative to the center part. The outer parts Q-X and P-Y perform their rotations j₃, y3, and j₄, y4, relative to part P-Q and part R-S. As each leaf spring element is considered as a rigid rod, the roll motions can be neglected. The angles are collected in 4´1 position vectors

where y_1F, y_2F and y_1R, y_2R describe the momentary shape of the front and the rear part of the left (1) and the right (2) leaf spring.

A fully dynamic description of a solid axle with two five link leaf spring models would result in ¦ = 6+2*8=22 degrees of freedom. Compared to the beam-element model this is a really significant reduction.

But a dynamic description of the five link leaf spring model still includes some high frequent modes which will cause problems in the numerical solution of the equations of motion. As mass and inertia properties of the leaf spring model parts are small compared to the solid axle, a quasi static solution of the internal leaf spring deflection should be accurate enough within the overall vehicle model.

A quasi static solution provides the position vectors of the leaf spring parts as functions of the axle position vector, y_1F = y_1F (y_A), y_1R = y_1R (y_A) y_2F = y_2F (y_A), y_2R = y_2R (y_A). Hence, the subsystem solid axle with two leaf springs has only f=6 degrees of freedom.

Initial Shape and Pretension

At first it is assumed that the leaf spring is located in the xz-plane of the leaf spring fixed frame x_L, y_L, z_L and its shape in the design position can be approximated by a circle which is fixed by the points X, C and Y. By dividing the arc X -Y into 5 parts of equal length the position of the links P, R, S, Q and the initial values of the angles j₀₁, y₀₁, j₀₂, y₀₂, j₀₃, y₀₃, j₀₄, y₀₄ can be calculated very easily.

In design position each leaf spring is only preloaded by a vertical load which results in zero pretension forces in the y_L-direction, , and zero pretension torques around the z_l-axis, , , , . In addition the torques around the x_l-axis vanish, , , , .

To transfer the vertical preload F₀ to the front eye bushing and the shackle, the joints P, Q, R, S must provide torques around the y_L-axis, Fig. 6. The pretension forces in the front eye bushing , and in the shackle F_0S, can easily be calculated from the equilibrium conditions of the five link leaf spring model,

where u_YZ refers to the unit vector in the direction of the shackle, and , are the x and z components of the vector from pointing from X to Y. The pretension torques in the leaf spring joints around the y_L-axis, , , , follow from

where r_ij , i = P,Q,R,S , j = X,Y are vectors pointing from i to j.

Compliance

The leaf spring compliance is defined in the design position by the vertical and the lateral stiffness, C_V and C_L. In Fig. 7a the leaf spring is approximated by a beam which is supported on both ends and is loaded in the center by the force F.

The deflection and the force F are related to each other via the stiffness c

If we transfer the beam model to the five link leaf spring model and look at the front half, Fig. 4, one will get

where a defines the length of one link, and small deflections in the x_L, z_L plane were assumed. The torques around the y_L-axis in the joints P and Q are proportional to the deflection angles j₁ and j₃

The equilibrium condition results in

The leaf spring bending mode due to a single force can be approximated very well by a circular arc. Hence, the relative angle between connected links is equal, j₁=j₃=j and Eq. 23 can be simplified to w=3aj or or . From Eqs. 24 and 25 follows

Using Eq. 22, one finally obtains

where the beam stiffness C was replaced by the vertical leaf spring stiffness c_V. Assuming symmetry, the stiffnesses in the rear joints are given by c_j2 = c_j1 and c_j4 = c_j3. The stiffnesses around the vertical axis c_y1, c_y2, c_y3 and c_y4 can be calculated in a similar way. The torsional stiffness of the leaf spring is neglected in this approach.

Actual Shape

The energy of a flexible system achieves a minimum value, E ® Min, in an equilibrium position. The energy of the five link leaf spring model is given by

where w_X is the 3´1 displacement vector and c_B is the 3´3 stiffness matrix of the front eye bushing, w_SR, w_SL are the radial and lateral shackle displacements, and c_SR, c_SLdenote the corresponding stiffnesses.

According to Eqs. 18 and 19, the actual shape of the leaf spring is determined by the position vectors y₁ = [j1,y1,j3,y3]^T and y₂ = [j2,y2,j4,y4]^T. If the leaf spring energy becomes a minimum, the following equations will hold

As the shackle displacements w_SR, w_SL do not depend on y₁ and the front bushing displacement vector w_X does not depend on y₂ the conditions in Eq. 29 form two independent sets of nonlinear equations ¦₁(y₁,y_A) = 0 and ¦₂(y₂,y_A) = 0, where y_A denotes the dependency of the actual position and orientation of the solid axle. These equations are solved iteratively by the Newton-Algorithm. Starting with initial guesses , one gets an improvement by solving the linear equations

Here, the Jacobians can be calculated analytically.

Leaf Spring Reaction Forces

The actual forces in the front leaf eye bushing are given by

where F_0B denotes the pretension force and C_B, d_B are 3´3 matrices, characterizing the stiffness and damping properties of the front leaf eye bushing. The displacement vector w_x in the front leaf eye bushing depend on the generalized coordinates y₁ and y_A which describe the actual shape of the front leaf spring part and the actual position and orientation of the solid axle. By solving Eq. 30, y₁ is given as a function of y_A. Hence, w_x only depends on y_A and its derivative can be determined by

where

The radial and lateral components of the shackle forces can be calculated from

where F_0S represents the pretension force, u_SR, u_SL are unit vectors in the radial and lateral shackle direction, and c_SR, c_SL, d_SR, d_SL are constants characterizing the stiffness and damping properties of the shackle. The shackle displacements w_SR and w_SL depend on the generalized coordinates y_2R and y_A which describe the actual shape of the rear leaf spring part and the actual position and orientation of the solid axle. Similarly to Eq. 32 the displacement velocities are given by

Finally, the shackle force reads as

Forces Applied to the Axle

The leaf springs act like generalized force elements in this approach, Fig. 8. Guidance and suspension of the solid axle is done by the resulting force

and the resulting torque

where r_AB1 = r_AB1(y_A) , ... r_AS2(y_A) describe the momentary position of the front eye bushings and the shackles relative to the axle center.

As the forces in the front eye bushings F_B1, F_B2 and the shackle forces F_S1, F_S2 depend on the axle state y_A, _A only, the resulting force and the resulting torque are also functions of the axle state only. Since, hereby each leaf spring acts as a generalized force element, it can easily be integrated into the vehicle framework. By suppressing high frequent leaf spring eigenmodes, it is perfectly adopted to real-time application.

Bending Modes

The quasi-static approach reproduces all significant bending modes of the leaf spring, Fig.9. A leaf spring is stiffer in the lateral direction than in the vertical direction. Hence, a displacement in the front eye bushing is noticeable only on lateral leaf spring deflections.

Model Performance

The five link leaf spring model was integrated into a ve-DYNA Ford Transit vehicle model.

Using this model at the rear axle instead of a poor kinematic approach means only 85% more computer run time. Hence, real time applications are still possible. The simulation results are in good conformity to measurements, Fig. 10. The nonlinearity in the spring characteristics is caused by an additional bump stop and by the change of the shackle position during jounce and rebound. Obviously, the five link model is accurate enough.

Free Body Module

Position and Orientation

To describe the momentary state of the body E the frame x_E y_E, z_E located in the center of gravity is used. In addition, sensor points S monitor position, velocity, and acceleration at specific body points, Fig. 11.

The frame B is fixed to the vehicle. The suspension of body E on the vehicle, frame B may consist of force elements and/or rubber mounts. The road-fixed frame 0 is considered as an inertial frame. The position of frame B with respect to the road-fixed inertial frame 0 is given by the position vector

The orientation of the frame axes is described by a rotation matrix. Three elementary rotations are put together. The sequence

results in

Hence, position and orientation of the vehicle-fixed reference frame are described by 6 generalized coordinates x_B,y_B,z_B and a_B, b_B, g_B.

The position and orientation of the elastically suspended body with respect to the reference frame B is defined by

and

Generalized Speeds

The velocity of the reference frame B with respect to the inertial frame 0 is given by

The velocity denoted in the inertial frame can be transformed to the reference frame

By this, the orthogonality of the rotation matrix

was already taken into consideration.

The angular velocity of the reference frame B with respect to the inertial frame may be expressed directly in reference frame B

The 6 components of n_0B,B and w_0B,B will be chosen as generalized speeds now. First order kinematical differential equations connect this speeds with the derivatives of the

and

where the solution of Eq. 48 is given by

The momentary state of the reference frame B is fully characterized by 6 generalized coordinates x_B, y_B, z_B, a_B, b_B, g_B and 6 generalized speeds n_0Bx, n_0By, n_0Bz, w_0Bx, w_0By, w_0Bz.

The velocity and the angular velocity of the elastically suspended body with respect to the inertial frame 0 is given by

where the derivative of the position vector and the angular velocity of the elastically suspended body follow from the Eqs. 41 and 42. They read as

and

By using the components of the velocity

and the angular velocity

as generalized speeds, Eq. 50 can be written as a set of kinematical differential equations

and

Where as the 6 generalized coordinates x_E, y_E, z_E, a_E, b_E, g_E describe the position and orientation of frame E relative to frame B, the 6 generalized speeds n_0Ex, n_0Ey, n_0Ez, w_0Ex, w_0Ey, w_0Ez are the components of the absolute velocity and angular velocity of body E.

Accelerations

The accelerations of body E with respect to the inertia frame 0 can be expressed in reference frame B. They read as

where

and

follow from the Eqs. 53 and 54.

Force Elements

If a force element is attached to the chassis at point i and to the body at point j, the momentary position of force element ij will be defined by

where

r_Bi,K , r_Ej,K are given by data, and r_BE,B follows from Eq. 41. The actual length can be calculated from

and the unit vector

describes the momentary direction of the force element. If denotes the initial length of the force element, the displacement of the force element will be formed by

The displacement velocity follows from

Using Eq. 60, Eq. 61, and

where

and

where ¦ describes an arbitrary spring/damper characteristic.

Equations of Motion

Applying liner and angular momentum to the elastically suspended body, one obtains

and

where m_E, Q_E,B denote mass and inertia tensor of the free body, F_E,B, T_E,B are the resulting forces and torques applied to the free body, and g_,B is the vector of gravity expressed in the body fixed reference frame. These equations are coupled with the chassis equations of motion only by the applied forces and torques. Due to the particular choice of generalized speeds, no mass or inertia coupling terms appear.

By using this modeling technique, Seibert and Rill (1998) showed that the comfort of a passenger car is significantly influenced by the engine suspension system. The free body model can also be used to model an elastically suspended driver's cab, Rill (1993).

Subsystem Drive Train

Generic Model Structure

The subsystem drive train, Fig. 12, interacts on one hand with the engine and on the other hand with the wheels. Hence, the angular velocities of the wheels w1, , w4 and the engine or the gear output angular velocity w₀ respectively are input quantities. For this reason, the calculation of the engine torque and the dynamics of the wheel rotation are performed in other subsystems. Via the tire forces and torques, the drive train is coupled with the steering system and the vehicle frame work.

The drive train model includes three lockable differentials. The angular velocities of the drive shafts w_S1: front left, w_S2: front right, w_SF: front, w_SR: rear, w_S3: rear left, w_S4: rear right are used as generalized coordinates.

The torque distribution of the front and rear differential is 1:1. If r_F and r_R are the ratios of the front and rear differential, one will get

The torque distribution of the center differential is formed by

where t_F, t_R denote the torques transmitted to the front and rear drive shaft, and µ is a dimensionless drive train parameter. A value of µ = 1 means front wheel drive, 0< µ <1 stands for all wheel drive, and µ = 0 is rear wheel drive. If the ratio of the center differential is given by r_C then

holds.

Equation of Motion

The equation of motion for the drive train is derived from Jordain's Principle, which reads as

where Q_i is the inertia of body i , _i denotes the time derivative of the angular velocity, t_i is the torque applied to each body, and dw_idescribe the variation of the angular velocity. Applying Eq. 75 for the different parts of the drive train model results in

Collecting all terms with dw_S1, dw_S2, dw_SF, dw_SR, dw_S3, dw_S4 and using the abbreviations n=1-µ, , and finally leads to three blocks of differential equations

which describe the dynamics of the drive train. Due to its simple structure, an extension to a 6x6 or 8x8 drive train will be straight forward.

Drive Shaft Torques

The torques in the drive shafts are given by

and c_S0, c_S1, c_S2, c_S3, c_S4, c_SF, c_SR denote the stiffnesses of the drive shafts. The first order differential equations can be arranged in matrix form

where

represents the vector of the angular velocities,

contains the torsional angles in the drive shafts,

is the excitation vector, and

forms a 7x6 distribution matrix.

Locking Torques

The differential locking torques are created by an enhanced dry friction model consisting of a static and a dynamic part

The dynamic parts are modeled by a torque proportional to the differential output angular velocities

where d_LF, d_LC, d_LR are damping parameters which have to be chosen appropriately. In steady state operating conditions, the static parts , , will provide torques, even if the differential output angular velocities are equal. From the Eqs. 80, 81, and 82, one gets

By this locking torque model, the effect of dry friction inside the differentials can also be taken into account.

Numerical Solution

The equations of motion 80, 81, and 82 can be combined in a matrix differential equation

where w, Dj are given by the Eqs. 85, 86, and the mass matrix M is built by three 2x2 submatrices

where the elements of M_F, M_C, and M_R follow from the Eqs. 80, 81, and 82. The vector of the generalized torques is written as

Because the model also includes the high frequent drive shaft vibrations the differential equations for the drive train are stiff. Hence, implicit integration algorithm should be used for the numerical solution. Vehicle dynamic equations can be solved very effectively by a modified implicit Euler algorithm, Rill (2004).

The implicit Euler-formalism for Eq. 92 and Eq. 84 results in

where h is the integration step size, and the superscripts k and k+1 indicate the states at t and t+1. Applying the Taylor-expansion to q at and w^k, one gets

By using the Eqs. 84 and 96 the second term on the right side can be written as

Now, the implicit algorithm in Eq. 95 can be approximated by

which finally results in

where the partial derivatives ¶q/¶Dq and ¶q/¶q can be calculated quite easily.

Partial Derivatives

Only the dynamic locking torques , and depend on the angular velocities. Hence, one obtains

The change of q with respect to Dj leads to a 6x7 matrix

The term which is finally needed in Eq. 100 is symmetric and reads as

where the abbreviations

were used.

System Performance

Locking the differential improves the traction of a vehicle. In Fig. 13 the simulation results of a vehicle with rear wheel drive starting on a µ-split surface are shown.

At first all differentials are unlocked. The left rear wheel which is running on a low µ-plate immediately starts spinning. At t = 2.5 s the rear differential is locked. Now, the locking torque, which is generated by the drive train model forces both wheels to run with the same angular velocity.

Conclusion

Vehicle modeling by subsystems make a large variety of applications possible. The combination of simple subsystems and modules results in a vehicle model with a minimum number of data and a very good run time performance. Such "light models" can be used to develop enhanced control strategies for electronic safety devices. Depending on the focus of interest, more and more subsystems and modules may be replaced by enhanced ones. Then, sophisticated design studies or a comfort analysis are possible. If the modified implicit Euler algorithm is also applied to the critical subsystems drive train and steering system the numerical solution of the overall vehicle model still will not be time consuming.

References

[1] ADAMS/Chassis 12.0 Reference Guide.

[15] tesis: www.tesis.de.

Paper accepted: June, 2005.

Presented at XI DINAME – International Symposium on Dynamic Problems of Mechanics, February 28th - March 4th, 2005, Ouro Preto. MG. Brazil.

Technical Editors: J.R.F. Arruda and D.A. Rade.

Publication Dates

History