Real-Time Multibody Simulation of Vehicle Wheel Suspensions of Different Topologies with Elastokinematic Properties

Driving simulators are becoming increasingly important in vehicle development, as they allow test drivers to evaluate the driving characteristics of a vehicle at an early stage of the development process. Vehicle dynamics models used for such real-time applications are usually simpliﬁed to achieve sufﬁciently low computation times, especially with respect to the elastokinematic properties of the wheel suspension. In this contribution, a real-time capable multibody simulation model for wheel suspensions taking into account their elastokinematics is presented. Linear-implicit integration methods enable real-time capable integration of the underlying numerically stiff equation of motion without major model simpliﬁcations. It is shown how the computational effort of the linear-implicit integration can be reduced by explicitly taking the suspension topology into account. A versatile process for the linearization of the equation of motion and a corresponding data model are presented. They allow to simulate different suspension topologies without requiring high implementation effort since only the model parametrization has to be adapted. With the presented methods, both a double wishbone and a multi-link suspension are simulated with high accuracy in real-time on a standard PC.


Introduction
The driving characteristics of passenger vehicles are significantly influenced by their wheel suspension.The wheel suspension guides the wheels relative to the vehicle body and thus defines the wheel position.The wheel position in turn affects the tire forces acting from the road onto the vehicle and thereby the vehicle's motion.The wheel position is not only defined by the kinematics but also by the elastokinematics of the suspension.The bodies of the wheel suspension are often connected by elastic bushings [1].Due to their elastic properties, the wheel suspension features some compliance under the influence of external forces.This results in additional wheel position changes that are superimposed on the purely kinematic guidance of the wheel.In automotive engineering, a variety of wheel suspension concepts with different topologies has emerged [1].An important part of their development is the design of their elastokinematic properties.It aims either to compensate the resulting wheel position changes or to convert them into desired wheel position changes to improve the vehicle driving characteristics [1].In addition to the objective evaluation of the driving characteristics, the subjective perception of test drivers plays a crucial role in vehicle dynamics development.Therefore, testing on driving simulators is an important tool in the development process.Vehicle dynamics models for such real-time applications are usually simplified to meet the required computation times.Simplifications often include the elastokinematic properties of the wheel suspension.In order to extend the validity of driving simulators and their applicability in the development process, vehicle dynamics models without such simplifications are desirable.However, real-time multibody simulation of vehicle wheel suspensions with elastic bushings is challenging.Elastic bushings in wheel suspensions have high stiffness and damping properties.Together with the low inertia of the suspension bodies, high-frequency and/or strongly damped motions occur in the suspension system.In contrast, the dynamics of, e.g., the wheel travel motion is comparatively slow.This means that the system features dynamics on different time scales.From a numerical point of view, the system's equation of motion has a high numerical stiffness, making stable and real-time capable numerical integration challenging.
In literature, one can find many approaches for real-time multibody simulation of wheel suspension with special focus on their elastokinematics.This includes, for example, modeling the wheel suspension as a macro joint [2], quasi-static approximation of the fast motions of the suspension bodies [3], or reducing the stiffness of bushings through additional stiffnesses and dampings [4].These approaches have in common that the dynamics and/or the elastokinematics of the suspension are modified or simplified.Alternatively, linear-implicit integration methods promise real-time numerical integration of numerically stiff equations of motion without the need for major model simplifications.Typical examples are the linear-implicit Euler method, Rosenbrock-Wanner or W-methods and the non-iterative HHT-α method [5,6,7].
In this contribution, the authors present a multibody model for real-time simulation of the dynamics and elastokinematics of wheel suspensions.Here, elastic bushings in the wheel suspension are modeled as stiff, (non)linear force elements.The underlying numerically stiff equation of motion is integrated with linear-implicit integration methods.Through that, real-time simulation without simplifications regarding the elastic bushings is achieved.An approach for further reduction of computation times based on the suspension topology is presented.The model is implemented in Matlab/Simulink.Its data model and implementation allow to simulate wheel suspensions with different topologies.For this, only the parametrization of the model has to be adapted.The paper is structured as follows.In section 2, the multibody systems of two representative wheel suspensions and their numerically stiff equation of motion are introduced.In section 3, challenges in real-time capable integration are addressed and the linear-implicit integration method LSRT2 is described.For linear-implicit integration, the system's equation of motion has to be linearized.The data model and calculation process for this linearization are explained in section 4. Section 5 presents the results in terms of computation times and model accuracy for a dynamic load case.The paper ends with the summary in section 6.

Multibody systems of suspension systems
In the following, we consider two different suspension systems, a double wishbone suspension and a multi-link suspension.Fig. 1 shows the simplified topology diagrams of these systems.In the double wishbone suspension, the wheel mount is connected to the vehicle body via the upper wishbone, the track rod, the trailing arm and the lower wishbone.The spring-damper unit consists of the damper tube, the piston rod and the spring which is divided into an upper and a lower part.The unit is connected to the vehicle body at the upper end and connected to the lower wishbone at the lower end.An anti-roll bar connects the left and right side (not shown) of the suspension.It is divided into two rigid bodies with a rotational spring in between.On its end, the anti-roll bar is connected to the damper tube via the anti-roll bar link.In total, the model contains 24 rigid bodies.All bodies are solely connected by force elements.The connections are represented by the black lines in Fig. 1, where one line stands for at least one force element.Elastic bushings are modeled as force elements with linear or nonlinear stiffness and damping characteristic curves in all translational and rotational directions.Bushing force elements with appropriate stiffness and damping values are also used to replace the few kinematic joints present in the model, e.g. the translational joint between the damper tube and the piston rod.The force characteristics of the spring-damper unit is modeled with point-to-point force elements.Furthermore, bump and rebound stops with highly progressive stiffness limit the wheel travel.In total, the model contains 55 force elements.The topology of the multi-link suspension differs from the double wishbone suspension, cf.Fig. 1.Here, the wheel mount is connected to the subframe via several control arms.Furthermore, the connections of the damper tube and the anti-roll bar are different.The multi-link suspension consists of 27 rigid bodies and 63 force elements.
As already mentioned, all bodies are connected solely by force elements without any kinematic constraints.According to that, the equation of motion can directly be formulated as a nonlinear ordinary differential equation in the form with the initial conditions r 0 and v 0 at time t = t 0 .Here, the vectors r and v denote the vectors of generalized coordinates and velocities, respectively.They are composed of the corresponding generalized vectors of the individual bodies as where m stands for the total numer of bodies.The generalized vectors of an indivial body l are defined as 0 x l , 0 y l , 0 z l and 0 v x,l , 0 v y,l , 0 v z,l are translational displacements and velocities of the body in the global coordinate system.Rotational kinematics is described through the Euler angles α z,l , α y ,l and α x ,l and the angular velocities B ω x,l , B ω y,l , B ω z,l in the body coordinate system.Corresponding to the number of bodies in the suspension system, the double wishbone suspension has 144 degrees of freedom and the multi-link suspension 162.The kinematics matrix K(r) describes the linear transformation of the generalized velocities to the time derivative of the generalized coordinates.The vector of generalized accelerations q(r, v, u) is a function of the generalized coordinates and velocities as well as further excitations, which are denoted by u.In the following, we consider excitations in the form of wheel forces.
Because of the high stiffness and damping properties of the bushings, the equation of motion features a high numerical stiffness.Numerical stiffness can be quantified through the ratio max i,j of the eigen values λ of the linearized equation of motion.A system is considered to be numerically stiff if the real parts of all eigen values are smaller than zero and the ratio max i,j (|λ i |/|λ j |) is much higher than one.This ratio is higher than 10 5 for the considered suspension systems, indicating a high numerical stiffness.
3 Linear-implicit integration of the equation of motion Numerically stiff multibody systems are challenging in terms of real-time capable numerical time integration.In real-time applications such as driving simulators, model outputs must be provided at a fixed, predefined output rate.Their calculation, i.e. the numerical time integration in this case, must be performed in each time step within the specified step size.On the one hand, implicit integration methods can enable stable integration of numerically stiff differential equations but may require the iterative solution of nonlinear equation systems.This involves a high computational effort that can even change during simulation depending on state and excitation of the system.Thus, numerical integration cannot be guaranteed to be completed within the step size specified by the real-time application.On the other hand, explicit integration methods usually have poor stability properties so that stable integration of numerically stiff differential equations is not possible.In contrast, linear-implicit integration methods are a suitable approach for stable and real-time capable integration of numerically stiff multibody systems [5].They are based on implicit integration methods but the equation of motion is integrated in a linearized form.Through that, the nonlinear equation systems to be solved for implicit integration are transformed to linear equation systems.
These can be solved with less effort which is even the same in each time step.Linear-implicit integration methods are therefore well suited for real-time applications.Many linear-implicit integration methods can be found in literature.Some of these were compared by the authors for simulation of the double wishbone suspension [8].Here, the integration method LSRT2 [9] was found to be a good compromise between real-time capability and accuracy.It belongs to the class of Rosenbrock-Wanner methods, which are the linear-implicit version of Runge-Kutta methods [7].LSRT2 is a single-step integration method with s = 2 stages and high numerical damping [9].The method is based on the linearization of the equation of motion according to with q lin = q(r lin , v lin , u lin ) and the Jacobians in the linearization state r lin , v lin and u lin .Based on the known generalized coordinates r k and velocities v k in time step k, the values in the next time step k + 1 are calculated according to Here, ∆t int is the constant integration step size and K lin is the kinematics matrix in the linearization state.For calculating ∆v k,i , the linear equation system has to be solved.The matrix I denotes the identity matrix.Afterwards, ∆r k,i is calculated according to The stage vectors rk,i , ṽk,i and ũk,i are defined as The parameters of the method s, b i , c i , a i j and γ i j are given in Table 1.
Table 1.Parameters of the method LSRT2 [9] s b 2 ) The choice of the linearization state r lin , v lin , u lin determines the computational effort and the accuracy of the method.The linearization could be performed only once at the beginning of a simulation with a suitably chosen linearization state.However, as the multibody systems considered here are strongly nonlinear, this approach is not appropriate.Indeed, the linearization has to be updated during the simulation.One obvious possibility is to update the linearization in each integration step.Alternatively, several integration steps can be performed between an update of the linearization.In this way, the accuracy of the integration can be increased without increasing the computational effort too much.However, in our previous research it was found that integration with step size ∆t int = 1.0 ms and linearization in each integration step is a good choice for simulation of the double wishbone suspension with the LSRT2 method [8].
In each stage i = 1...s of an integration step, linear equation systems with the coefficient matrix must be solved (cf.(7)).The solution is done by means of LU decomposition and forward/backward substitution.Because of γ 11 = γ 22 , the coefficient matrix is the same in both stages of an integration step.Furthermore, the coefficient matrix only changes when the linearization (4) and thus the Jacobians J r,lin and J v,lin as well as the kinematics matrix K lin are updated.Only then, the coefficient matrix and its LU decomposition must be recalculated.
From (11) it can be observed that the sparsity pattern of the coefficient matrix results from the sparsity pattern of the matrices J r,lin , J v,lin and K lin .Thus, it is strongly connected to the topology of the wheel suspension.Fig. 2 shows the sparsity patterns for the double wishbone and the multilink suspension (cf.topology diagrams in Fig. 1).The sparsity patterns are shown for the complete suspensions while the topology diagrams only show the left side of the suspensions.From Fig. 2 it can be seen that the coefficient matrix is a block matrix with the colored blocks.The selected colors of the main-diagonal blocks refer to the colors in the topology diagrams in Fig. 1.The bodies that connect the left and right side of the suspension are colored blue.These are the antiroll bar and, in the case of the multi-link suspension, also the subframe.The subsystem wheel, wheel mount and control arms is colored in red and the spring-damper subsystem is colored in yellow.The non-zero elements outside the matrix diagonal results from the connections between these subsystems.In the double wishbone suspension, the anti-roll bar is connected to the springdamper unit on both sides.The subsystem wheel, wheel mount and control arms is also connected to the spring-damper unit, but not to the anti-roll bar.Accordingly, there are only non-zero blocks on the main and on the lower and upper diagonal of the block matrix.The same applies for the multi-link suspension.This sparsity pattern can be utilized to solve the linear equation system (10) more efficiently.For that, the linear equation system is written with the coefficient matrix as a block matrix: This linear equation system is then divided into several smaller equation systems that are solved subsequently.Firstly, the linear equation system with the coefficient matrices and the right-hand side vectors are solved.Afterwards, the equation systems are solved.This procedure requires the matrices A ii and Ãii to be invertible.Savings of computation time result from the fact that several small linear equation systems are solved instead of one large equation system.This requires less computational effort for the LU decomposition and forward/backward substitution.Further savings are achieved if only the non-zero elements of the matrices A i j (i = j) are considered in the matrix multiplications in equations ( 14)-( 16).It should be noted that the partitioning of the coefficient matrix according to Fig. 2 may not be the best in terms of computational efficiency.However, it is closely related to the suspension topology and therefore straightforward to formulate.4 Data model and linearization of the equation of motion For the linear-implicit integration method introduced in the previous section, the generalized accelerations q lin , the Jacobians J r,lin , J v,lin , J u,lin and the kinematics matrix K lin in the linearization state r lin , v lin , u lin have to be determined.In this section, the underlying data model and the calculation process are described.Fig. 3 shows the data model for parametrizing a wheel suspension system.In the data model, all model elements are defined with their numerical parameters.In addition, each model element is assigned a unique index for identification.With these indices, the suspension topology is represented in the data model, as described in the following.Each marker is a local coordinate system on a body, cf.Fig. 3 (right), and thus uniquely allocated to that body.According to that, the index of that body is included in the marker properties (cf.Fig. 3 (left)).The force elements, e.g. the bushing in Fig. 3 (right), are attached to the markers.Exactly two markers are assigned to one force element.In turn, each marker may be assigned to one force element only.
The corresponding marker indices are therefore included in the force element properties.In this way, information about which bodies are connected by which force elements, i.e. the suspension topology, is included in the data model.Different suspension systems can be parametrized by adding further bodies and connecting them with each other according to the suspension topology.For calculating and linearizing the equation of motion, the parameters of all bodies, markers, bushing and P2P force elements are listed in tables.The calculation process is shown in Fig. 4. The process starts with the calculation of the body kinematics.All bodies are solely connected through force elements but not kinematically.Thus, their kinematics is calculated separately based on their generalized coordinates and velocities as well as their parameters listed in the body parameters table.The calculated kinematic quantities are stored in another table, in which each entry is assigned to one body together with its index.In the second step, the kinematics of the markers is calculated one after the other.Required quantities of the body kinematics are taken from the body kinematics table according to the body index stored in the marker properties.Again, the results of the marker kinematics are listed in a table, in which each entry is assigned to one marker together with its index.The same logic is used to calculate the force elements, surface and volume forces.With these results, the Newton-Euler equations are set up for each body separately to obtain its generalized accelerations.The vectors of generalized accelerations of all bodies are then written into the global vector of generalized accelerations q lin .The global Jacobians J r,lin , J v,lin and J u,lin are obtained analytically.In all calculation steps in Fig. 4, not only the kinematic and kinetic quantities but also their partial derivatives are calculated and output.They are composed in the last calculation step and written into the global Jacobians.The global kinematics matrix is also set up in the last calculation step.This calculation process and its implementation can be used for arbitrary wheel suspension topologies.Also, the calculation steps for the linear-implicit integration of the equation of motion (cf.section 3) are independent of the topology, as long as the sparsity pattern of the coefficient matrix according to Fig. 2 is not explicitly considered when solving the linear equation systems.This means that only the parametrization according to Fig. 3 must be adapted to simulate wheel suspension systems with different topologies.

Results
In this section, the results in terms of computation times and accuracy for a dynamic load case are presented.

Computation times
The real-time factor is used to evaluate the computation times of the model.It is the ratio between the required computation time and the specified duration of a simulation, i.e. a value below one means real-time capability.To determine the real-time factor, the Simulink model was compiled with the Simulink Coder and C code was generated [10].Intel oneAPI Library was incorporated in the compilation process.It contains implementations of Basic Linear Algebra Subprograms (BLAS) and Linear Algebra Package (LAPACK) [11].In the compilation process, Simulink Coder produces BLAS and LAPACK calls to speed up linear algebra operations.The produced code was executed on a standard desktop computer with an Intel Core i7-3770 CPU at 3.4 GHz.Numerical integration was performed with the LSRT2 method with step size ∆t int = 1.0 ms and linearization in each integration step.Table 2 shows the obtained real-time factors.Results with and without utilizing the sparsity pattern of the coefficient matrix (11) (also cf.Fig. 2) are given.The real-time factors of the multi-link suspension are higher than the factors of the double wishbone suspension because of its higher number of model elements and degrees of freedom.However, all real-time factors are smaller than one, meaning simulation of both suspension systems is real-time capable.Utilizing the sparsity pattern reduces the required computation time by up to 20%.This results in real-time of 0.59 and 0.73 for the double wishbone and the multi-link suspension, respectively, which are well below one.

Simulation of a dynamic load case
To show the accuracy of the integration method, a dynamic load case was performed with both suspension models.We considered a force step input of the longitudinal wheel force according to where the positive direction points in driving direction.The upward vertical wheel force was F z (t) = 5000 N.All wheel forces were applied on both wheels at the respective wheel center.Integration of the equation of motion was done with the LSRT2 method with step size ∆t int = 1.0 ms and linearization in each integration step.For comparison, reference results were obtained with the Matlab/Simulink solver ode23tb [12].Simulations with this solver require computation times far from real-time, but serve as a reference for evaluating the accuracy of the LSRT2 method.Fig. 5 shows the simulation results for both suspension systems.After the wheel force step, weakly damped oscillations of the longitudinal wheel displacement and the bushing force occur.According to the elastokinematic properties of the suspensions, the wheel displacement tends towards a static equilibrium with a displacement of a few millimeters.The deviations of the results obtained with LSRT2 from the reference results are within the line width.Table 3 gives the normalized root mean square errors of the results compared to the reference results.The values refer to the time range shown in Fig. 5 and are normalized by the mean value of the reference results.The error values are in the order of 1% for both suspension systems.From that we can conclude that the simulation results of the suspension systems obtained with the LSRT2 method have a good accuracy.

Conclusion
This contribution presents the real-time capable simulation of numerically stiff multibody systems on the example of vehicle wheel suspensions.Here, high numerical stiffness results from elastic bushings connecting the suspension bodies.Real-time capable simulation without model simplifications is achieved using linear-implicit integration methods, in this case the method LSRT2.These methods transform the nonlinear equation systems that must be solved for implicit integration into linear equation systems.It is shown how the suspension topology can be taken into account in the solution of these linear equation systems for further reduction of the computational effort.With the underlying data model, wheel suspension systems with arbitrary topologies can be parametrized.Furthermore, the linearization process allows to linearize the equation of motion of arbitrary suspension systems.The applicability of the approach is demonstrated using two suspensions systems with different topologies, a double wishbone and a multi-link suspension.With Table 3. Normalized root mean square errors of the results shown in Fig. 5 Double Wishbone Suspension Multi-Link Suspension Wheel displacement 0.00243 0.01286 Force 0.00751 0.01545 real-time factors of 0.59 and 0.73, both suspension systems can be simulated in real-time.The simulation results for the considered dynamic load case have a high accuracy with an error in the order of 1%.

Figure 1 .
Figure 1.Simplified topology diagrams of the double wishbone (left) and the multi-link suspension (right).Only the left side of the suspension systems is shown.

Figure 2 .
Figure 2. Sparsity patterns of the coefficient matrix (11) for the double wishbone (left) and multi-link suspension (right)

Figure 3 .
Figure 3. Simplified data model for parametrizing a suspension system (left) and connection of two bodies with a force element (right)

Figure 4 .
Figure 4. Process for calculating and linearizing the equation of motion

Figure 5 .
Figure 5. Longitudinal wheel displacement and force in a bushing of the double wishbone and multi-link suspension obtained with LSRT2 compared to the reference with ode23tb

Table 2 .
Obtained real-time factors