Hindawi Publishing Corporation Differential Equations and Nonlinear Mechanics Volume 2007, Article ID 49157, 13 pages doi:10.1155/2007/49157

Research Article

Numerical Effectiveness of Models and Methods of Integration of the Equations of Motion of a Car

Marek Szczotka, Szymon Tengler, and Stanislaw Wojciech

Received 3 January 2006; Revised 19 October 2006; Accepted 2 November 2006 Recommended by Roger Grimshaw

The paper presents models of car dynamics with varying complexity. Joint coordinates and homogenous transformations are used to model the motion of a car. Having formulated the models of the car, we discuss the influence of the complexity of the model on numerical efficiency of integrating the equations describing car dynamics. Methods with both constant and adaptive step size have been applied. The results of numerical calculations are presented and conclusions are formulated.

Copyright © 2007 Marek Szczotka et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

From the mechanical point of view, vehicles are complex multibody systems with the treelike structures. Their modelling has to take into account not only the complex structure of kinematic chains but also flexibility of some elements and joints, friction in joints, and sophisticated relations that describe the forces acting on the wheels from the road surface. The complexity of car models depends on applications and is strongly influenced by availability of parameters. In the paper, we discuss three kinds of models with varying complexity.

The most complex models of a car may consist of very many subsystems, the motion of which can be described by a large number of degrees of freedom. These models allow structural models of suspensions and drive systems of a car to be taken into consideration. Usually the models are formulated using absolute coordinates and special professional software with a large scale of generality (e.g., MSC Adams, DADS, SIMPACK) [1]. Such calculation models are very time-consuming, and require a large amount of specific data. On the other hand, they allow us to obtain results which are in good correspondence with experimental measurements [2]. Our model formulated using joint

coordinates and homogenous transformations, denoted in the paper by {C}, belongs to this group.

The second group consists of models which are fairly general but are characterized by a smaller number of degrees of freedom. Here, usually functional models of suspensions and drive systems are used [3, 4]. The number of necessary parameters is much smaller, but they require measurements in order to define proper functional characteristics of the subsystems. These models are widely applied in the initial phase of research, in design of active subsystems, and in control of vehicle motion. In our paper, the model denoted by {I} represents this group.

The simplest models are usually applied in real-time simulations and in control of vehicle motion. The simplifying assumptions can result in unsatisfactory correspondence of calculation and measurement results. This is the reason that such models are usually devoted to particular types of cars or their motion. The model denoted in the paper by {S} is from this group.

Apart from some models formulated for specific applications, the equations of motion of a car form a system of differential ordinary equations which is nonlinear and stiff. Using the models described, we will compare numerical methods of integrating the equations of motion and discuss whether they can be applied to real-time simulations of car motion.

2. Models of a car

The models presented in the paper are obtained using joint coordinates and homogenous transformations. Joint coordinates allow us to reduce the number of generalized coordinates of the system considered, and in contrast to absolute coordinates reduce the number of constraint equations. However, the mass matrix is not diagonal and its coefficients depend on generalized coordinates [5]. Homogenous transformations are often used in robotics [6, 7]. This method allows us to transform coordinates from local coordinate system to another coordinate system ^ (Figure 2.1) using only one mathematical operation:

r = Tr',

cfcO cfsOsp — sfcp cfsOcp + sfsp x0~

sfcO sfsOsp + cfcp sfsOcp — cfsp y0

—sO cOsp cOcp z0

0 0 0 1

R is the rotation matrix, r0 is the position vector, s means sin function, and c means cos function.

Figure 2.1. Transformation of coordinates.

In our paper, we use the Euler ZYX angles (f -yaw, O-pitch, p-roll) and thus the matrix R takes the form

Application of joint coordinates and homogenous transformations in formulating the equations of motion of multibody systems using Lagrange equations are described in detail in [5, 8], which also demonstrate how equations of motion change when a body is attached to an existing kinematic chain. Additional information which allows us to include the constraint equations in vehicle models is presented in [9].

Below, we present the models of car dynamics with varying levels of complexity.

Complex model {C}. The model is shown in Figure 2.2.

The model includes structural models of suspensions and drive system. Formulation of this model, the equations of motion, and the verification of results are described in detail in [9]. The car body is treated as a rigid body and its vector of generalized coordinates has the following form:

where Xn, yn, Zn are coordinates of the origin of local coordinate system {n} in global coordinate system {}, fn, 9n, pn are ZYX Euler angles.

The front suspension is a McPherson type. Its main elements (k = 1,2—left- or right-hand side ofthe vehicle) are the following:

(a) suspension column, whose relative rotations are described by the components of

cfcO cfsOsp - sfcp cfsOcp + sfsp R = sfcO sfsOsp + cfcp sfsOcp - cfsp . -sO cOsp cOcp

the vector:

(b) hub with two degrees of freedom in relative motion:

(Z,k) = JZ(Z,k) f (Z,k)j T,

Figure 2.2. Complex model of the car {C}: (a) general scheme; (b) front McPherson suspension; (c) rear suspension with trailing arms.

(c) wheel rotation angle (with respect to the hub):

JW ,k)

q—' = [e (W>k)J, (2.7)

(d) control arm angle of rotation with respect to the car body:

q(H k = y^H j, (2.8)

(e) spherical joints connecting the tie rod with the hub and the control arm with the hub (this means that kinematic closed loops occur in the model).

The rear suspension is composed of the following elements (k = 1,2):

(a) shock absorber, whose relative motion is described by angles:

q(T ,k) = (T k e (T k p(T k ]T, (2.9)

(b) piston rod with two degrees of freedom in relative motion:

q(F ,k) = ,k) y(F ,k) ]T, (2.10)

(c) trailing arm:

q(H ,k) = [p(H,k^, (2.11)

Figure 2.3. Intermediate model of a car {I}.

(d) wheel rotation angle:

,(W,k) =

(2.12)

(e) spherical joints that link trailing arms with absorbers. The vector of generalized coordinates of the whole system can be written as follows:

where qV is the vector of generalized coordinates of the vehicle, which includes elements of suspension (2.5)-(2.12) and the generalized coordinates of the car body (2.4), q^ is the vector of generalized coordinates of the steering system, and qD is the vector of generalized coordinates of the drive system.

The equations of motion and constraint equations form a system of 75 + nL + n^ + + n^) equations where n^L), n^R), n^D), nW) are the numbers of rigid finite elements used for discretization of the semiaxles ( L,R), drive shaft ( D), and the steering column

Intermediate model {I}. The model is described in detail in [4]. The vehicle is composed of the following subsystems ( Figure 2.3):

( i) car body, which has six degrees of freedom, (ii) suspensions, each of which has four degrees of freedom,

(iii) wheels with tyres,

(iv) steering mechanism (simplified model) with two degrees of freedom,

(v) drive system (optional).

The motion is described by 28 generalized coordinates, which are the components of the following vector:

(2.13)

6 Differential FnnatinTT; anrl Nnnlinpar TVTprhimir';

Figure 2.4. Simplified model of the car.

where qn is the vector of coordinates of the car body, as described in (2.4),

qZk = [x(z) yf 4Z) 4zf, k = 1,... ,4, (2.15)

x[z), y^, z[z) are displacements of suspensions, are steering angles, fk are wheel rotation angles, yL is the displacement of the rack in the steering gear, fW is steering wheel angle, and fD = [f1D) ■ ■ ■ fnD!]T is the vector of generalized coordinates of the drive system.

This model has been used to solve an optimization problem described in [4, 9]. The main purpose of this task was to find courses of drive torques of independent electric engines which ensure car stability in extreme conditions (such as a sudden change of lane).

Simplified model {S}. In this model, the car body is treated also as a rigid body. The vector of generalized coordinates has the form as in (2.4). The motion of wheels is described by generalized coordinates which are the components of the vector:

qw = [f 1 f2 f3 f4~\ , (2.16)

where fi is the rotational angle of wheel i.

The vector of generalized coordinates of the vehicle has the following components:

qU = [xn yn zn dn fn f1 f 2 f3 . (2.17)

In this model, vertical flexibility of the suspension is taken into account by modification of stiffness coefficients of the tyres (Figure 2.4).

The equations of motion of all models presented can be written in the following general form of differential-algebraic equations:

M(q)q + D(q)R = F(f, q, q), (2.18a)

O(q) = 0, (2.18b)

where R is the vector of unknown reaction forces.

The constraint equations (2.18b) can be written in a differential form, and then the standard numerical procedures for integrating a set of differential ordinary equations can be applied.

3. Results of calculations

The Dugoff-Uffelmann model of tyres [10] has been applied in all models presented. The data used are appropriate for a small car.

All models presented were verified experimentally. In Figure 3.1, we compare the results of measurements with those obtained by calculations using the complex model {C} [2] (B-measurements, W-calculations).

In Figure 3.2, the test vehicle with measurement instrumentation is presented. The measured accelerations presented in Figure 3.1, were obtained from sensors placed on car body and suspension, as presented on photos in Figures 3.2(c) and 3.2(d). Motions (on straight and circular path) over one, two, or three identical obstacles (200 mm width and 25 mm height) with sharp edges have been analyzed.

The main vehicle parameters (used in all models analysed in the paper) are shown in Table 3.1.

The calculation results obtained for models with varying complexity are presented in Figure 3.3. The motion on a ^-split surface is considered. It can be observed that all models give similar results.

The Runge-Kutta method of the fourth order with constant step size (h = 0.001 second) has been applied in all cases.

In the paper, we present the calculation results obtained for the following methods of integration of differential equations [11] describing the car motion.

(A) Methods with constant step size:

(i) Runge-Kutta method of the fourth order (RK),

(ii) implicit Euler method (IE).

(B) Methods with adaptive step size:

(i) Runge-Kutta-Fehlberg method (RKF),

(ii) Bulirsch-Stoer method (BS).

(C) Methods for stiff differential equations with adaptive step size:

(i) Resenbrock method (RS),

(ii) Bulirsh-Stoer-Daufhardt method (BSD).

The calculations have been carried out for three different cases which are connected with a sudden change in vehicle motion. For methods with constant step size, we have considered the following values of the step size: h = 0.001 for RK (Runge-Kutta of the fourth order) and h = 0.01 for IE (implicit Euler) when the models {I} or {S} have been used. For the IE method and the model {C}, a special time step has to be selected. The solution is stable only if a small time step is used. For the IE method, we have to calculate the inverse of the mass matrix, so the time of simulation rises significantly due to the increase in the number of the equations in the system. The calculations allow us to conclude that the IE method is not effective when it is used to solve the equations of motion of the complex model {C}. Below, we present the calculation results for three different manoeuvres of a car.

200 150 ~ 100 50 0

-50 -100 -150„

3 obstacles on straight way, 60km/h

|t-| —r —-f —-

____1 _ j 1 ~ lu l\ ' 1

/vjl f y III ^

- _ J^—S (1 - h /fj li Vr

-----W" V ], \ L' * \r\h r r

'll / Mi1 /»» 1 1

In ^ TP t" fl [ ( I \

8 0.9 1 1.1 1.2 1

Time (s)

15 10 5 0 5 10 15

1.4 1.6 Time (s)

- W-calculations -W-calculations

- B-measurements ---B-measurements

(a) (b)

Figure 3.1. Comparison of measurement and calculation results: (a) vertical acceleration of the front right hub; (b) vertical acceleration of the car body.

Figure 3.2. Experimental verification: (a) measurement devices; (b) Daytona© steering wheel (angle and torque measurement); (c) car body sensor; (d) front suspension sensor.

3.1. Lane change manoeuvre. It has been assumed that steering angles of wheels of the vehicle change as shown in Figure 3.4(a).In Figure 3.4(b), we present the lateral displacement of the car body (its centre of gravity-CG) for all models. One can see that all models give similar results.

The comparison of numerical effectiveness of the methods used to solve this task is presented in Table 3.2.

Table 3.2 presents the number of calls (of the function that generates set of (2.18)) and calculation time (on PC with Intel Pentium IV 2.4 GHz processor) for each model presented in Section 2, and for each method of integration.

Table 3.1. The main vehicle parameters.

Parameter

Car body mass

Car body moments of inertia I(XX), I(YY) , I(ZZ) CG height

Unsprung mass front/rear Wheel\track base Unloaded tire radius Front\rear susp. coeff of stiffness Steering rack ratio Scrub radius

650 (kg) 300, 800, 850 (kg m2) 530 (mm) 20.5/18 (kg) 1320\2200 (mm) 300 (mm) 30\40 (kN/m) 75 (mm/rad) ~ 4 (mm)

p = 0.9

p = 0.35

--- Model C

-Model I

-■-■ Model S

(a) (b)

Figure 3.3. The influence of model complexity on calculation results: (a) ,w-split surface; (b) yaw velocity of the car body.

3.2. Travel over an obstacle. The second case analysed concerns a vehicle moving over an obstacle with the left-hand side wheels. The dimensions of the obstacle are presented in Figure 3.5(a). The obstacle causes a lateral displacement of the vehicle (it starts turning) if the driver does not take any action to correct the direction. This displacement is shown in Figure 3.5(b).

In this case, simulation time has been limited to 5 seconds. The results are presented in Table 3.3.

For intermediate model {I}, the step size h = 0.005 has to be applied to the IE method. With this step size, the IE method requires much more time for computing than the standard RK routine. It is important to note that the Runge-Kutta method requires the shortest calculation time for both {C} and {I} models.

Steered wheels input

123456789 10 Time (s)

30 25 ^ 20 £ 15 ^ 10 5 0

1 2 3 4 5 6 7 8 Time (s)

-- - ModelS -■-■■ Modell -Model C

Figure 3.4. (a) Assumed steering wheel angles; (b) lateral displacement of the car body CG.

Table 3.2. Comparison of methods for lane change manoeuvre.

Model c i s

Method Number of Calculation Number of Calculation Number of Calculation

function calls time min/s function calls time min/s function calls time min/s

RK4 50 000 5/36 50 000 4/00 50 000 0/30

IE — — 61000 4/01(* ) 21021 0/11(*)

BSD 256029 54/18 43 713 4/47 881 0/01

ROS 1 581255 151/32 442 255 63/52 31742 0/18

BS 147257 34/42 63 445 8/11 21 105 0/012

RKF 103615 23/07 118 250 9/28 11696 0/07

(*)-step size h = 0.01.

3.3. Braking on a ^-split surface. The last example analysed is related to braking on a ^-split surface. The situation (surfaces with different coefficients of adhesion on the road) is presented in Figure 3.3(a). The courses of braking torques (identical for each wheel) are presented in Figure 3.6.

As in the two previous cases, we have performed simulations for this task using all models and methods considered, and the results are presented in Table 3.4.

The time step h = 0.001 has to be selected for the IE method in ensure stability. This proves that the IE method is fast, but only for very simple models.

Below, we present results of calculations obtained when commercially available code has been applied (MSC.ADAMS). The case of motion analyzed was identical to our B example (moving over obstacle), and complex model {C} of vehicle (strictly equivalent model developed in MSC.ADAMS/Car) has been passed to the solver.

As can be seen in Table 3.5, GSTIFF solver is very effective implementation. CON-STANT_BDF solver, also designated to stiff systems, even if not as fast as GSTIFF, is mentioned to be more stable [12]. Nonstiff solver in ADAMS such as RKF45 uses appropriate

Car body CG lateral

Time (s)

......RK 4 .....ROS

-BSD ----RKF

(a) (b)

Figure 3.5. Travel over an obstacle: (a) obstacle dimensions; (b) lateral displacement of car body CG caused by obstacle.

Time (s)

Figure 3.6. The course of braking torque applied.

Table 3.3. Effectiveness of the methods for travel over an obstacle.

Model c i s

Number of Calculation Number of Calculation Number of Calculation

Method function calls time min/s function calls time min/s function calls time min/s

RK4 25 000 5/36 25 000 1/54 25 000 0/15

IE — — 61000 4/13(**) 10 521 0/6(*)

BSD 256029 24/12 32 025 2/20 1026 0/06

ROS 1 581255 123/43 232 800 14/42 24 984 0/17

BS 147257 17/02 59 458 3/42 14 776 0/10

RKF 103615 12/31 57717 3/35 8311 0/09

(*) time step h = 0.01 (s) (**) time step h = 0.005 (s).

Table 3.4. Comparison of methods for braking on a ^-split surface.

Model c i s

Method Number of Calculation Number of Calculation Number of Calculation

function calls time min/s function calls time min/s function calls time min/s

RK4 15000 3/24 15000 1/26 15000 0/09

IE — — 196685 5/40(**) 6300 0/03(*)

BSD 60874 8/16 18717 1/22 1708 0/02

ROS 580412 61/22 365847 42/32 25342 0/24

BS 35831 3/54 34556 2/20 6777 0/10

RKF 66838 8/09 26695 1/45 5566 0/08

(*)-step size h = 0.01 (**)-step size h = 0.001.

Table 3.5. Comparison of simulation time in MSC.ADAMS\Car (model {C}, obstacle passing).

Solver type min/s

GSTIFF 2\45

CONSTANT_BDF 11\49

RKF45 unstable

conversion of DAE into ODEs system of equation. This solver is not recommended: in our case, the simulation takes very long calculation time (even a couple of hours) or becomes unstable, depending on parameters. However, taking into account differences in modelling methods used (general coordinates in ADAMS, joint coordinates in our models), presented results can be treated only as a "loose comparison." At this point, we are not able to make any strict comparison MSC.ADAMS versus our models, since both codes are completely different in modelling approach, numerical methods applied, and so forth.

4. Conclusions

It can be concluded that there is not one specific method that allows us to quickly obtain a stable and accurate solution for all models. The complexity of the vehicle model strongly influences numerical effectiveness of methods considered.

For the complex model {C}, the best method is the Runge-Kutta method with constant step size. However, if the manoeuvres analysed are relatively moderate, the BSD method for stiff systems is much more effective. We have analysed an example in which BSD takes only 14 seconds of calculation time while the RK requires 5 minutes and 36 seconds. In the case of shock phenomena (such as travel over an obstacle), the best results were obtained when a method with a constant step size was applied. In other cases, we suggest to use BSD for stiff systems as a good and stable method.

For the intermediate model {I}, all of the methods have a similar effectiveness, except the ROS method. In this case, as for the complex model {C}, the effectiveness of the methods selected strongly depends on the manoeuvre analysed. If simulation does not contain any abrupt changes, the BSD method is the fastest.

The conclusions for model {S} are different. Here the IE method gives the most rapid solution. The BSD method also gives similar results, but in some applications a constant step size method is preferred. The simplified model {S} analysed in the paper can be applied in a driving simulator, when calculations have to be performed in real time. For this application, the implicit Euler (IE) method should be selected. In driving simulator application time has to be predictable.

We hope that the examples presented address the problems of dynamic simulations of vehicles. For the cases of motion presented, standard procedures with an adaptive integration step require much longer time for solution than the simple Runge-Kutta method with the constant step size; this is because the Jacobin matrix required to perform integration in more advanced methods has to be calculated.

References

[1] E. J. Haug, Intermediate Dynamics, Prentice-Hall, Englewood Cliffs, NJ, USA, 1992.

[2] M. Szczotka, "Numerical analysis of the dynamics of vehicle: mathematical model and its verification," The Archive of Automotive Engineering, vol. 3, pp. 241-260, 2005.

[3] W. Grzegozek, The Modelling of Vehicle Dynamics with Stabilizing Control of Braking Forces, vol. 275 of Scientific Papers in Mechanics, Press of TU Cracow, Cracow, Poland, 2000.

[4] M. Szczotka and S. Wojciech, "Model for simulation of vehicle dynamics," The Archive ofMe-chanical Engineering, vol. 50, no. 4, pp. 347-362, 2003.

[5] I. Adamiec-Wojcik, Modelling Dynamics of Multibody Systems Using Homogenous Transformations, vol. 3, Scientific Papers of Bielsko-Biala University, Bielsko-Biala, Poland, 2003.

[6] J. J. Craig, Introduction to Robotics, Addison-Wesley, Reading, Mass, USA, 1988.

[7] R. P. Paul, Robot Manipulators: Mathematics, Programming, and Control, MIT Press, Cambridge, Mass, USA, 1981.

[8] E. Wittbrodt, I. Adamiec-Wojcik, and S. Wojciech, Dynamics of Flexible Multibody Systems. The Rigid Finite Element Method, Springer, New York, NY, USA, 2006.

[9] M. Szczotka, Modelling of vehicle motion taking into account various drive systems, Ph.D. thesis, University of Bielsko-Biala, Bielsko-Biala, Poland, 2004.

[10] H. Dugoff, P. S. Fancher, and L. Segel, "An analysis of tire traction properties and their influence on vehicle dynamics performance," in Proceedings ofFISITA International Automobile Safety Conference, vol. 79, Detroit, Mich, USA, May 1970, SAE 700377.

[11] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, Cambridge University Press, Cambridge, UK, 2nd edition, 1992.

[12] MSC.Software, 2003, MSC.ADAMS Documentation.

Marek Szczotka: Department of Mechanics and Informatics, University of Bielsko-Biala,

43-309 Bielsko-Biala, Poland

Email address: marek.szczotka@wp.pl

Szymon Tengler: Department of Mechanics and Informatics, University of Bielsko-Biala,

43-309 Bielsko-Biala, Poland

Email address: stengler@ath.bielsko.pl

Stanislaw Wojciech: BOSMAL Automotive Research and Development Center,

43-300 Bielsko-Biala, Poland

Email address: swojciech@ath.bielsko.pl

Copyright of Differential Equations & Nonlinear Mechanics is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.