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

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.


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.

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: where R is the rotation matrix, r 0 is the position vector, s means sin function, and c means cos function.
Marek Szczotka et al. 3 In our paper, we use the Euler ZYX angles (ψ-yaw, θ-pitch, ϕ-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: (a) suspension column, whose relative rotations are described by the components of the vector: (b) hub with two degrees of freedom in relative motion: (c) wheel rotation angle (with respect to the hub): (d) control arm angle of rotation with respect to the car body: (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: (b) piston rod with two degrees of freedom in relative motion: (c) trailing arm: (d) wheel rotation angle: (e) spherical joints that link trailing arms with absorbers.The vector of generalized coordinates of the whole system can be written as follows: where q V 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 K is the vector of generalized coordinates of the steering system, and q D is the vector of generalized coordinates of the drive system.The equations of motion and constraint equations form a system of 75 + n ϕ are the numbers of rigid finite elements used for discretization of the semiaxles (L,R), drive shaft (D), and the steering column [9].{I}.The model is described in detail in [4].The vehicle is composed of the following subsystems (Figure 2.3):

Intermediate model
(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: where q n is the vector of coordinates of the car body, as described in (2.4), k are displacements of suspensions, δ (z) k are steering angles, ϕ k are wheel rotation angles, y L is the displacement of the rack in the steering gear, ϕ W is steering wheel angle, and nD ] 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).{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:

Simplified model
where ϕ i is the rotational angle of wheel i.
The vector of generalized coordinates of the vehicle has the following components: (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(t,q, q), (2.18a) 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.

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.001second) 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.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.Marek Szczotka et al. 9

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.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   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.

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.

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

3. 1 .
Lane change manoeuvre.It has been assumed that steering angles of wheels of the vehicle change as shown in Figure3.4(a).In Figure3.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 Table3.2.

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

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 Figure3.3(a).The courses of braking torques (identical for each wheel) are presented in Figure3.6.
4))where x n , y n , z n are coordinates of the origin of local coordinate system {n} in global coordinate system {}, ψ n , θ n , ϕ n are ZYX Euler angles.

Table 3 .
2. Comparison of methods for lane change manoeuvre.

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

Table 3 .
4.Comparison of methods for braking on a μ-split surface.