Phase plane analysis applied to non-explicit multibody vehicle models

A new methodology for constructing stability maps (phase-plane analysis) is presented and validated for application to complex multibody vehicle models implemented in Multibody Dynamics simulation software (Adams®). Traditional methodologies are developed to be applied to explicit mathematical models. Given the complexity of some special multibody systems, particularly in vehicle dynamics, simplifications are needed to apply this stability analysis technique. The main limitation when using simplified models is the need to neglect components which could have a significant influence on the dynamic behavior of the system and therefore on its stability. In the proposed methodology it is not necessary to have access to explicit mathematical models of multibody systems. Thus, the stability map of a vehicle model can be constructed by considering highly nonlinear dynamic elements, such as tires and silent-blocks components, modeled using the nonlinear finite element technique.


β V
Vehicle sideslip angle measured in the vehicle-fixed reference frame ω V x , ω V y , ω V z Components of vehicle angular velocity (vehicle-fixed reference frame) x V , y V , z V Axis of the vehicle-fixed reference frame Components of vehicle velocity in the vehicle-fixed reference frame F V , T V Applied force and torque (vehicle-fixed reference frame) Resultant torque of tire forces and moments in the z axis of the vehicle body reference frame

Introduction
Following a multibody dynamics approach, the equations of a mechanical system are a set of nonlinear differential-algebraic equations. This form of the equations allows the accurate systematic analysis of complex mechanical systems in different scenarios. However, for the stability analyses of the system it is necessary to have access to the linearized system dynamics [1].
The linearization of the dynamics equations can be performed in different ways, depending on the strategy selected to formulate the Differential Algebraic Equations (coordinate selection and kinematic constraints formulation) [1]. Several methods for the linearization of multibody dynamics have been published in the literature [2][3][4][5][6].
The stability analysis of a dynamic system about an equilibrium point involves the eigenvalue analysis of the corresponding linearized equations of motion. In the literature, several methods, such as the Lyapunov second method, are applied for analyzing the stability of dynamic systems [7], including both linear and nonlinear systems. The main difficulty in applying these methods is that the dynamic systems may need to be explicitly expressed [8].
Escalona [9] proposed a stability analysis method based on two coordinate projections and a special eigenvalue analysis of differential algebraic equations. It is applied to the stability analysis of steady motions of a vehicle in order to improve the computational cost, because applying Floquet's theory to large multibody systems involves very high computational costs.
For some applications, it is fundamental to have suitable stability maps available, which are typically constituted by curves in a phase plane (phase portrait). These methodologies allow comprehensive characterization of the dynamic response of a nonlinear system, which often has multiple steady-state solutions, providing an easy way to analyze the steady-state solution that a set of initial conditions will converge to.
For example, the stability analysis of vehicles, using stability maps (phase-plane analysis), is applied to improve the designs or to the development of control systems strategies [10][11][12].
Due to the necessity that the dynamic systems are explicitly expressed to construct these maps, several authors proposed methods based on extremely simplified vehicle models [10][11][12][13][14][15][16][17][18][19]. Therefore, instead of a complex multibody vehicle model with nonlinear components (bushing, dampers, structural behavior, etc.), the most frequently used simplified vehicle model is a two-state vehicle dynamics model. This model is known as the bicycle model or the two-wheel model (single-track model). To incorporate the nonlinearities of the tires, these models use the brush tire model [10][11][12]17] or the Magic Formula (Pacejka) tire model [13-16, 18, 19].
Flavio Farroni [20] proposed an improved model, a quadricycle planar model, characterized by two states referring to in-plane vehicle body motions (lateral and yaw motions). The model implements a nonlinear tire model for lateral behavior and two different analytical methods, the phase plane and the handling diagram, are applied to better analyze the in-curve vehicle lateral stability. With the same vehicle and tire model, Wencai Sun [21] analyzed the stability region of the vehicle, using the phase plane, employing the GA (Genetic Algorithm) method to solve for the bifurcation points.
Craig [22] extended the phase portrait to three states in order to include the longitudinal vehicle dynamics (drive/braking control inputs). Craig used a four-wheel vehicle model with a Brush tire model.
One recent application of phase plane analysis, using a single-track bicycle model, is to develop the path-tracking controller for autonomous vehicles [23,24].
The main problem with using simplified models for the stability analysis of complex multibody models is the fact that they neglect components which could have a significant influence on the dynamic behavior of the system and therefore on its stability. A full vehicle model contains an important number of bodies connected not only by kinetic joints, but also by nonlinear springs, shock absorbers, silent-blocks (bushings), elastic bodies (anti-roll bars, etc.) and highly nonlinear tires [25]. In addition, tire behavior is affected by variables such as temperature, internal pressure and wear, with a significant impact on the vehicle performance [26,27].
As the extensive bibliography reveals [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], the simplified models neglect these elements or the highly nonlinear dependence of the tires with combined slip (longitudinal/lateral), temperature, or internal pressure. A local linear model can be accurate enough to analyze the behavior of the vehicle in the vicinity of the equilibrium point and at low lateral acceleration, but it is not suitable if the distance between the different equilibrium points is large or the lateral accelerations are medium or high [21]. It means that the accuracy of the stability analysis of these models is clearly compromised.
For future improvements to vehicle control systems, it will be necessary to know and monitor the nonlinear dynamics of the vehicle, taking into account simultaneous inputs [22].
In order to improve the understanding of the nonlinear dynamics of the vehicle, in this paper a new methodology is presented to build stability maps (phase-plane) of multibody vehicle models implemented in a commercial software. The main advantage of the proposed methodology is that it is not necessary to have access to the explicit mathematical models of the multibody systems (the mathematical model does exist inside the simulation software, but the equations are not visible) and that it can be applied to any type of vehicle modeled in Multibody Dynamics (MBD) simulation software.
The existing MBD simulation software (e.g., the one used in this work, Adams ® ) is a simulation, analysis, and optimization tool, with high accuracy in modeling complex systems. This type of software allows the user to implement complex multibody models, avoiding mathematical formulation but including all dynamic properties [25].
Therefore, the proposed methodology allows taking advantage of the capabilities of this type of software, such as the possibility of including flexible bodies (linear or nonlinear), nonlinear elements (bushings, dampers, etc.), complex tire models, or to perform co-simulation with other software (FEM, real-time data, etc.).
Furthermore, the methodology easily allows the implementation of simultaneous driver inputs (steering, drive/brake forces, etc.) or control inputs (individual control of wheel torques, etc.). This paper does not aim to improve the traditional theories on the dynamics of nonlinear systems (Poincaré-Bendixson, etc.), but an advance in the construction of stability maps (phase-plane analysis) of complex multibody vehicle models.
The methodology was applied and validated with a complete three-dimensional model of a vehicle implemented in Adams ® (Fig. 2). This multibody system is easily implemented and simulated in a Multibody System Dynamics simulation software. In addition, it can be easily modified to analyze its sensitivity to design variables.

Methodology
The phase-plane plot is applied to two-state variable systems in order to analyze their stability. A complex multibody model has more than two state variables but, for many systems, its stability can be analyzed using only two state variables of one body. This is the case of ground vehicles where the phase portrait describes vehicle sideslip angle (β V ) and vehicle yaw rate (ω V z ) or vehicle body side slip angle (β V ) and its changing rate (β V ) in a planar diagram [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].
The sideslip angle is defined as follows: with v V x and v V y being the longitudinal and lateral velocities ( Fig. 1) in a noninertial vehiclefixed reference (x V , y V ).
A phase-plane plot for these two-state variables consists of curves of one state variable versus the other (ω V z (t) vs. β V (t)), where each curve is based on a different initial condition (β V (0), ω V z (0)). At each point in this phase-plane, a corresponding vector (β V ,ω V z ) is associated. By repeating this for sufficiently many points (grid) in the state space, a vector field plot is obtained. From the field plot, a phase trajectory can be sketched out for any initial condition.
The points on the phase-plane plot where β V = 0 ∧ω V z = 0 the system is considered to be in equilibrium. For a vehicle traveling along a curve at constant speed, this is equivalent to steady state motion. An analysis of the phase plane around that point would allow an assessment of whether the equilibrium is stable, unstable, or indifferent.
For the other points in the phase-plane plot, the system is not in equilibrium, which implies that the vehicle is moving in a transient mode. This movement is characterized both by its kinematic condition (β V , ω V z ), and by its movement tendency, which is analyzed as a function of the values ofβ V andω V z . Both obtaining and characterization of the equilibrium points, as well as the study of the motion tendencies at any other point (β V , ω V z ), have been widely analyzed in the state-of-the-art investigations, as discussed above, but always requiring the explicit formulation of the equations of the system. In general, if the equations of the system , the values β V ,ω V z can be obtained. When a dynamic system (e.g., a road vehicle) is implemented in MBD simulation software, a "black box" (this term is used under the consideration that the equations of the system are not explicitly known) is available, which can perfectly characterize the dynamic response of the system to defined inputs. If the inputs are time invariant, the system will reach, after a transient process, a steady state regime characterized by a kinematic condition (β V , ω V z ) in which the energy of the system is minimal. Taking all this into account, obtaining the equilibrium points (steady state circulation) of a vehicle is quite straightforward using an MBD simulation software by setting initial (and time invariant) conditions for the speed and the steering wheel angle. When the simulation is performed, after a transient period, a steady state regime characterized by its kinematics β V , ω V z is reached. This implies that the system is free, i.e., in "unforced" conditions. That is, there are no additional external forces or torques (F V , T V ). In case of applying external actions (forces or torques), the kinematic condition (β V , ω V z ) obtained in the steady state will be different and will be a function (f ) of the input (F V , T V ), Taking all of the above into account, what is proposed in the methodology developed in this work is a modification of the original system (vehicle driving at constant speed with a fixed steering wheel angle) by including external forces and torques. These forces and torques will be called "virtual" (F V v , T V v ), since they do not exist in the real system (vehicle), but they allow (force) the system to be taken to a certain kinematic condition (β V , ω V z ). In the real system, this condition represents a point where the vehicle is not in equilibrium. However, in the modified system (virtual system), this kinematic condition is reached in a steady state, forced by the added external forces and torques, after a transient period.
As will be demonstrated in Sect. 2.2, these virtual forces and torques, which force the virtual system to reach a certain kinematic condition (β V , ω V z ), are related to the response β V ,ω V z of the real system if it were in this initial condition. So, the presented method proposes that the vector field be obtained by the introduction of a virtual force and torque (in the multibody model of the vehicle) in order to impose each kinematic (initial) condition (β V , ω V z ). For a multibody vehicle model traveling at a constant longitudinal velocity and keeping the tires in contact with the road, the application of an external lateral force at ground level and an external yaw moment allow reaching a condition that could occur while driving.
The longitudinal movement of the vehicle is controlled by the traction forces (longitudinal tire forces) that allow maintaining a constant longitudinal velocity. In the virtual system, the vertical forces on the vehicle body and the pitch and roll torques are balanced by the ground reactions (tire forces, linked to the vehicle body through the steering and suspension mechanisms). Any situation other than equilibrium is impossible or would mean an accident (e.g., rollover) and is not the subject of this methodology. Therefore, it is not necessary to introduce a force to reach equilibrium in this direction.
Therefore, the vehicle sideslip angle and yaw rate (kinematic condition) depend on the virtual force and moment and the lateral forces and self-aligning torques of the tires. Thus, the stability of a vehicle traveling at constant longitudinal velocity can be analyzed using these two state variables (ω V z vs.β V ). As detailed in the following sections, this virtual force and torque allows the calculation of the corresponding vector (β V ,ω V z ) for each (β V , ω V z ) point.

Constrained multibody systems
The differential-algebraic equations of motion for multibody systems, such as a vehicle, can be expressed as follows [28]: where q is the composite set of generalized coordinates for the system, Q is the vector of generalized forces (applied forces and torques), q is the Jacobian matrix of the vector of constraints, λ is the vector of Lagrange multipliers associated with the constraints present in the system, M is the system mass matrix, and γ is the vector obtained by the following acceleration equation: For k bodies in the system, q takes the following expression (where k is the number of bodies in the system): where q i is the generalized coordinate vector associated with the rigid body i in a system. For a 3D vehicle model, the differential-algebraic equations of motion can be expressed as follows [28]: where r is the composite set of generalized coordinates, defined by three Cartesian coordinates for each body of the vehicle (origin of the body-fixed coordinate system with respect to the global inertial frame XYZ), ω is the body angular velocity expressed in the body-fixed coordinate system (centroidal frame xyz), M is the mass matrix, J is the inertia tensor with respect to the centroidal frame; and F A T A are the applied forces and torques. Considering a vehicle model formed by k bodies, the components in the equation are: M ≡ diag (m 1 I 3 , m 2 I 3 , . . . , m k I 3 ) , The Lagrange multipliers, obtained along withr andω as the solution of system (6), are related to the reaction forces and torques produced by the constraints of the system.

Phase plot of the multibody vehicle model
Considering as state variables the vector (ṙ m i , ω n i ), corresponding to the velocity variables of the body i coordinates m and n, of a body-fixed reference frame, the phase portrait can be generated using the method of building the vector field diagram. As mentioned above, in the case of analyzing the stability of a road vehicle, the variables to be considered are the vehicle body sideslip angle (β V ) and vehicle yaw rate (ω V z ). In the traditional approach, simplified models are used to provide an explicit formulation of the equations of the system. This allows the algebraic solution of the system of equations where the values of β V and ω V z are delimited in a predefined interval with maximum and minimum values, which define the state space under study.
Solving system (8) for enough points in the state space, a vector field diagram can be obtained. Thus, setting a pair of values (β V , ω V z ) that identify the kinematics of the vehicle body, a pair of values β V ,ω V z are obtained that allow assessing the dynamic response and the trend of the handling response of the vehicle.
In the case of a detailed approximation to the dynamics of the system, by means of a mathematical model based on MBD software, an explicit formulation of the equations is not available. To solve this issue, a methodology based on the introduction of a virtual force and torque (applied force and torque external to the system) that provide the system with the desired kinematic response (β V , ω V z ) is proposed. This virtual force and torque are introduced into the complete highly nonlinear model, without simplifications, to obtain the characterization of its dynamic response β V ,ω V z by solving with an MBD software. So, the methodology proposed to obtain the vector field is based on the introduction of a virtual force and torque vectors (F v , T v ) in the system (6) in order to impose each kinematic condition β V , ω V z , with the virtual system in dynamic equilibrium (steady state response). Assuming the road is flat, for a local-level reference frame with the origin at ground level (Fig. 2), the virtual force F ll v and torque T ll v are: Solving the system of Eqs. (6) at a constant longitudinal velocity v V x , when the vehicle reaches dynamic equilibrium, the kinematic condition (β V , ω V z ) is obtained. In this condition, the tire forces and moments (longitudinal, vertical, lateral, etc.), other external forces (e.g., aerodynamic) or the vehicle body angles (pitch and roll) are the same as if the condition is kinematically imposed.
The numerical integration of the differential-algebraic equations of motion is carried out in an MBD software. For each point β V , ω V z of the state space, the vector field associated is a function of the virtual force/torque imposed to reach that point, so This vector β V ,ω V z , as a function of the force and virtual torque, must be obtained. For this purpose, the whole vehicle is considered as an unconstrained system. The Newton-Euler equations of motion in the vehicle body frame are: The force equation in lateral direction is After the introduction of a virtual force and torque, the system reaches the dynamic equilibrium for each initial condition, sov V y = ω V x = v V z = 0. Equation (12) can be simplified as follows: where F V vy is the component of the virtual force in the y axis of the vehicle body frame. With the system in dynamic equilibrium, if the virtual force disappears, the equation in a lateral direction will be Therefore, the virtual force (F V vy ) is related to the inertia effect andβ V can be obtained as follows:v The virtual torque is applied in the vertical axis. It is assumed that the axes of the vehicle's reference system are the principal axes of inertia, so the equation in this direction is The system of Eq. (6) was solved to reach the dynamic equilibrium for each initial condition, soω V z = ω V x = ω V y = 0. Equation (16) can be simplified as follows: where T V vy is the component of the virtual torque in the z axis of the vehicle body frame. The virtual torque (T V vy ) is related to the inertia effect, so if this torque disappears, the equation will be Therefore,ω V z is evaluated asω Once the values of β V ,ω V z have been obtained, by means of Eq. (15) and (19), for each of the points defined in the state space the phase plot can be drawn.
The application of a virtual force and torque to achieve a dynamic equilibrium is valid for conditions far from the equilibrium points. It is a methodology applied to multibody vehicle models, so these situations must be possible to achieve during real driving, and the multibody model is able to simulate them.
The main limitation of the methodology is in those situations where the vehicle reaches the rollover situation. In such cases the virtual forces and torques do not allow an equilibrium condition to be reached. These situations are not the subject of the methodology. Its development is intended to be applied to the development of stability controls in a sliding situation and not in a rollover situation.

Results
In order to apply the proposed methodology, a multibody vehicle model is implemented and its simulation is performed in Adams ® commercial software developed by MSC Software Corporation, as shown in Fig. 2.
The vehicle model implemented in Adams ® has 65 moving parts, 3 flexible bodies, 26 nonlinear bushings, and 59 joints, which gives a model with 144 independent degrees of freedom. A detailed model of the front and rear suspension system has been implemented. The front suspension is McPherson type with a lower control arm and an antiroll bar, modeled as a flexible body (Fig. 3).
The rear suspension is a compound crank axle suspension mechanism, with a twistable cross-member (V shape) and an antiroll bar, both modeled as flexible bodies (Fig. 4).
For tires, the Pacejka PAC2002 tire model is used, with cornering and longitudinal slip considered simultaneously with combined slip and transient behavior.
Inputs required for simulating the model are the steering-wheel angle and drive torque. The first is constant and the drive torque is applied to maintain the longitudinal velocity constant. To properly define it, a PID (Proportional Integral Derivative) controller is used.
Several simulations were carried out at constant speeds and for different steering angles. The methodology described in point 2 has been applied in order to obtain the phase-plane. Figure 5 shows the phase-plane at 30 km/h for four steering angles (mean front wheels steering angles 0, 0.1, 0.2, and 0.3 rad).
As can be seen, there is one stable equilibrium point and a saddle equilibrium point for each steering angle. In the case of 0.3 rad of steering, the stable equilibrium point and the saddle equilibrium point are close together.
In order to verify the identification of the phase-plane stable equilibrium points, an openloop step steer maneuver of the vehicle model has been carried out (Fig. 6). At a constant velocity (30 km/h), the input steer angle is as shown in Fig. 6. As a result, the equilibrium points for each steering angle are shown in the following Table 1: The equilibrium points are the same as those identified in the phase-plane (Fig. 5).
The main advantage of the new methodology is that it allows the phase-plane to be defined whilst taking into consideration all the nonlinearities of the model. This is especially important when the vehicle is close to its limits, i.e., when one or more tires are saturated or close to saturation. This is the test case with a 0.3 rad steering angle. As can be seen in Fig. 7, the response of the vehicle is not the same as in the previous steps. The vehicle slip angle passes through a maximum value, stabilizing at a lower value.
To analyze why this vehicle slip reduction occurs, the lateral forces of each tire are checked. As can be seen in Figs. 7 and 8, the inner rear tire of the curve (rear left) is saturated. It reaches a maximum value of lateral force because it exceeds its slip limit (Fig. 9), as a consequence of a reduction in the vertical load.
This phenomenon is only identifiable if the nonlinearity of its rear axle is considered, in addition to the tires themselves. The twistable cross-member affects the kinematic and dynamic behavior of the compound crank axle suspension [29]. Due to this and other nonlinear elements (bushings, etc.), it shows a nonlinear relationship between the wheel angles (toein and camber) and the wheel displacement [30]. This nonlinear relationship significantly influences tire response.  Simplified models, such as those currently used, neglect these nonlinearities and do not consider the influence of each tire individually. Therefore, they would not be able to take these effects into consideration when determining the phase plot.
Therefore, this limitation of traditional methods conditions the validity of the phase-plot in the vicinity of the vehicle limits. Furthermore, this questions the validity of the phase-   Obviously these traditional methods have a lower computational cost. In any case, the application for which the methodology is developed does not require the simulation to run in real time.

Conclusions
The stability maps are a tool that allows the stability of complex systems to be analyzed. The traditional methods for their construction, curves in a phase plane (phase portrait), require that the dynamic systems be explicitly expressed. Therefore, in the traditional approach, simplified models are used that ignore components which could have a significant influence on the dynamic behavior of the system and, therefore, on its stability.
Road vehicles are among the complex systems to which stability maps are applied, for example, for the development of control systems. Simplified vehicle models are used that do not take nonlinear components into consideration. Examples are nonlinear bushing, antiroll bars, twistable cross members, tires, etc., which are linearized or neglected. This paper proposes a methodology for the construction of stability maps, that does not need explicit equations and that considers all the nonlinearities of the model. The methodology is based on the introduction of virtual forces, in the multibody model, that allow a dynamic equilibrium of the system to be reached under different conditions. The methodology is applied to a passenger car. A complex virtual model was created using Adams ® commercial software. The model includes different nonlinear components: antiroll bars, coil springs, dampers, tires, and a twistable cross-member in the rear axle.
Several simulations were carried out at constant speeds and for different steering angles. The phase-planes for each steering angle have been obtained. To validate the identification of the phase-plane stable equilibrium points, an open-loop step steer maneuver of the vehicle model has been carried out. The correlation has been proven.
The open-loop maneuver shows the influence of the nonlinearities at high lateral acceleration, when one or more tires are saturated or close to saturation. With a simplified model, the behavior of the system would be different, and it would not allow the real behavior to be reproduced.
Therefore, this limitation of traditional methods conditions the validity of the phase-plot in the vicinity of the vehicle limits. That is, when the point of stability is very close to nonequilibrium points. Furthermore, this questions the validity of the phase-plot at points that are far from the equilibrium point, even when this equilibrium point is far from the vehicle limits.
Three-dimensional multibody vehicle models are usually completely modeled in Multibody Dynamics simulation software, such as Adams ® from MSC Software, for simulation, analysis, and optimization, including highly nonlinear dynamic elements, such as tires and silent-blocks. Therefore, applying the phase-plane to these non-explicit models is needed so as to be able to continue improving the designs and control systems on these models. Funding Note Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Conflicts of interest/Competing interests The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.