Stability and local bifurcation analyses of two-wheeled trailers considering the nonlinear coupling between lateral and vertical motions

The nonlinear dynamics of two-wheeled trailers is investigated using a spatial 4-DoF mechanical model. The non-smooth characteristics of the tire forces caused by the detachment of the tires from the ground and other geometrical nonlinearities are taken into account. Beyond the linear stability analysis, the nonlinear vibrations are analyzed with special attention to the nonlinear coupling between the vertical and lateral motions of the trailer. The center manifold reduction is performed leading to a normal form up to third degree terms. The nature of the emerging periodic solutions, and, thus, the sense of the Hopf bifurcations are verified semi-analytically and numerically. Simplified models of the trailer are also used in order to point out the practical relevance of the study. It is shown that the linearly independent pitch motion affects the sense of the Hopf bifurcations at the linear stability boundary. Namely, the constructed spatial trailer model provides subcritical bifurcations for higher center of gravity positions, while the commonly used simplified mechanical models explore the less dangerous supercritical bifurcations only. Domains of loss of contact of tires are also detected and shown in the stability charts highlighting the presence of unsafe zones. Experiments are carried out on a small-scale trailer to validate the theoretical results. A good agreement can be observed between the measured and numerically determined critical speeds and vibration amplitudes.


Introduction
Vehicle handling and stability are critical factors in road transportation; hence, they became relevant research topics a long time ago. First, analyses of linear vehicle models were in focus (see, for example, [22,36]), which provided the basic rules of vehicle design. However, some of the phenomena of the field remained unexplored since they are related to the nonlinear dynamics of the vehicles. Thanks to the developed mathematical tools and computational methods, the nonlinear analysis could come to the front in the academic area and started to be utilized in the industry too. Stability charts, phase portraits and nonlinear bifurcation diagrams were composed for the most intricate nonlinear problems of vehicle dynamics, like the shimmy motions of steered wheels (see [2][3][4]19,21,25]) and aircraft landing nose gears (see [33,34]), the wobble motion of motorcycles and bicycles (see [9,17,18,20,23,26,27]). Vehicle handling in close to critical maneuvers is also in the focus of recent studies (see [12,24]), which allow the development of the new generation of vehicle motion control that utilizes the results of nonlinear analyses.
In this paper, the very interesting and complex stability problem of trailers is under investigation. Trailers are very commonly used on the roads in very different combinations with respect to the mass ratio of the towing vehicle and the trailer. Several road accidents happen due to the not appropriately chosen amount of payload and/or payload position/distribution, which easily leads to the so-called snaking and/or rocking motions of the trailer. These unwanted vibrations are also often associated with high towing speeds and are often generated by some initial impulses/perturbations caused by the driver or other circumstances (such as wind gusts or uneven roads), see [1,5,7,16,28,35]. In the past, theoretical and numerical investigations have been made on the stability properties of car-trailer combinations considering different level of complexity of the investigated mechanical models, see the 3-DoF model of [5,16], the 4-or 6-DoF models of [1] or the 32-DoF model of [28]. A systematic experimental investigation by Darling et al. [7] was also carried out using a car-adjustable trailer combination and drew the conclusions that the key factors that affect the stability properties are the yaw inertia, the mass distribution and the position of the trailer axle.
In theoretical stability analyses, the bounce, roll and pitch motions of the trailer are often separated from the lateral dynamics, which is a commonly used and accepted approach in vehicle dynamics. Namely, the use of in-plane models is a conventional solution in the field. Sometimes, the single track model, the socalled bicycle model, is also applied to simplify the analysis, see, e.g., [2][3][4]19]. As a result, the governing equations have a relatively simple form and analytical calculations can be performed regarding the nonlinear properties of car-trailer systems. These models partially explain unsafe parameter zones and related nonlinear phenomena. For example, it is shown in [2,4] that the Hopf bifurcations emerging at the critical towing speed are subcritical. Namely, the global stability of the rectilinear motion is not ensured in the linearly stable towing speed range, and large enough perturbations may lead to unwanted large amplitude vibrations of the trailer.
In this study, a spatial mechanical model of a twowheeled trailer with non-rigid wheel suspensions is used, where all the yaw, pitch and roll motions are considered, see Model A in Fig. 1a. In order to have a relatively low degree-of-freedom model, the towing car is imitated by a lateral spring and damping at the king pin. This simplification may seem inappropriate, but in this study, we focus on the effect of the pitch motion on the nonlinear vibrations. The main contribution of the paper is the analysis of the nonlinear coupling between the vertical and lateral motions of the trailer. Our hypothesis is that the linearly independent pitch motion affects the local bifurcations, namely it changes the sense of the Hopf bifurcations at the linear stability boundaries. In order to verify this assumption, three models are compared to each other in our analysis, see Fig. 1. These models represent different levels of complexity, beginning with the most complex, spatial model (Model A) with 4 DoF, see panel (a). Model B in panel (b) is a reduced model with blocked bounce/pitch motion, but it considers the effect of the roll motion on the stability. Model C is the simplest, 2-DoF inplane model with blocked pitch and roll motions, see panel (c).
The rest of the paper is organized as follows. The mechanical model of towed, two-wheeled trailers and the main sources of the nonlinearities are presented in Sect. 2. In Sect. 3, the governing equations are linearized around the rectilinear motion and linear stability charts are shown for the mechanical models of Fig. 1. The effects of the wheel suspension parameters and the payload positions on the linear stability are analyzed. The nonlinear vibrations are investigated in Sect. 4. Namely, we perform local bifurcation analysis with the center manifold reduction and numerical investigations using continuation. The limit cycles of the three different models are compared to each other both analytically and numerically. Unsafe zones and regions where the loss of contact of tires happens are also shown for the spatial model. Based on the mechanical model, a small-scale experimental rig is designed on which some measurements are performed. The experimental setup and experimental results are shown in Sect. 5. Finally, our theoretical, numerical and experimental results are concluded in Sect. 6.

Mechanical model and generalized coordinates
The spatial mechanical model (represented as Model A in Fig. 1) of two-wheeled trailers can be seen in  The trailer is moving in the (X, Y, Z ) ground-fixed coordinate system. The towing speed v is kept constant in the X direction, while the movement of the king pin (point A) is not constrained in the Y direction. The trailer is considered as a rigid body with mass m and mass moment of inertia tensor J C at the center of gravity C with respect to the moving reference frame (x, y, z). The longitudinal and vertical positions of the center of gravity are given relative to the axle of the wheels by the parameters e and h, respectively. We consider that the trailer is symmetric with respect to the (x, z) plane, i.e., point C is located in this center plane of the trailer. The constant vertical distance between the king pin and the ground is denoted by h 0 , the track width is 2b and the caster length, i.e., the longitudinal distance between the axle of the massless wheels and the king pin, is l. In the mechanical model, the pitch dynamics of the trailer is considered; accordingly, the overall stiffness and the damping of the tires and their suspensions are denoted by the parameters k and c, respectively. The interaction between the towing vehicle and the trailer is modeled here by the lateral stiffness k lat and by the lateral damping c lat . However, this assumption seems to be an oversimplification, but it is often used in the literature (see [31,32]) in order to keep the complexity of vehicle models and the cor- Fig. 2 The spatial, 4-DoF model of towed two-wheeled trailers. a Axonometric view with the generalized coordinates and geometrical parameters. b The interpretations of the distances d R and d L and the lateral tire forces responding mathematical formulas on a manageable level. Accordingly, the number of parameters is moderate in our model. However, in order to help the understanding of the following sections, all the parameters with their notations and names are listed in Table 1, close to the stability charts where their numerical values are first referred.
The spatial motion of the trailer would have six degrees of freedom without any constraints. However, the constant towing speed assumption manifests in the rheonomic geometric constraint X A = vt. In addition, the vertical position of the towing point is fixed Z A ≡ h 0 . As a result, the system has n = 6 − 2 = 4 degrees of freedom, and the vector of generalized coordinates can be composed as: where ψ is the yaw angle, ϑ is the pitch angle, ϕ is the roll angle of the trailer, and u is the lateral displacement of the king pin, see Fig. 2a.
In order to avoid the presence of kinematic constraints, which could be originated in the rolling of the wheels, we rather use an elastic tire model that is based on the creep-force idea considering the side slip of the wheels [21]. Thus, no kinematic constraint is taken into account; correspondingly, the system is holonomic and the equations of motion can be derived with the Lagrange equation of the second kind [11].

Nonlinearities
Details of the derivation of the governing equations can be found in [14,15]. Here, we present the essential assumptions of the mechanical model, pointing out the different sources of the relevant nonlinearities and nonsmoothness. On the one hand, geometric nonlinearities are present in the model, but we do not emphasize them here, since they are straightforwardly provided by the derivation of the equations of motion with mathematical rigor. On the other hand, nonlinearities originated in the tire forces and in the suspension system have crucial roles in the nonlinear behavior of the trailer. Thus, they are in focus in this subsection.
The active forces acting on the trailer are the gravitational force G; the lateral tire forces F tire R and F tire L (see Fig. 2b) acting at the (possible) contact points R and L of the right and the left tires, respectively; the wheel suspension forces F R and F L ; the lateral force F lat = −k lat u − c latu acting at point A due to the deformation of the lateral spring and damper. The effect of the tires is taken into account by means of the lateral tire forces only, which are calculated based on Pacejka's Magic Formula [21] via the side slip angles α R and α L of the right and left wheels, respectively. Thus, where refers to R or L accordingly to the right or left wheel. With this formula, we assume that the lateral tire force depends linearly on the vertical load N = F cos ϑ cos ϕ, which is calculated as the vertical Zcomponent of the suspension force F . Namely, we neglect the effect of the unsprung mass by considering zero masses for the wheels. The characteristics μ(α) are where B is the stiffness factor, C is the shape factor, D is the peak factor, and E is the curvature factor. In practice, these parameters are not constant, and they may depend on the vertical load, the static and dynamic coefficients of friction, the temperature, the camber angle, etc., see [21]. For the sake of simplicity, we neglect these dependencies and consider fixed values in our analysis. We also notice that we neglect the self-aligning moments of the tires, which have negligible effect on the final results.
On the contrary, we pay special attention to the rocking motion of the trailer, namely we take into account that the tires can detach from the ground and the vertical load N becomes zero. As mentioned above, this corresponds to zero suspension forces in our model. Hence, we introduce the piece-wise smooth formula: for both of the right F R = F s (d R ) and the left F L = F s (d L ) suspension forces, where the distances d R and d L are measured between the points R and R and L and L , respectively, as shown in Fig. 2b. To handle the non-smooth characteristic of the suspension force, we use the so-called Heaviside step-function H (s) [8], where δ(s) is the Dirac-delta function. Thus, in Eq. (4), L min corresponds to the minimal distance, where the spring of the suspension is compressed to its minimal length. Namely, for d < L min , we consider higher stiffness and damping for the suspension. The distance L max relates to the maximal distance, where the suspension of the wheel is totally expanded. Accordingly, and where L 0 is the free length of the spring. To imitate the collision in the suspension, one should choose k ∞ k and c ∞ c. Note that we did not take into account that the suspension dampers can be highly nonlinear. We used this simplification because most caravans only have leaf springs without any dampers.
Note that the ReLU(s) = max{0, s} function in Eq. (4) guarantees that the suspension force cannot be negative even when the suspension expands rapidly (i.e.,ḋ 0). For the steady-state conditionḋ = 0, the schematic characteristics of the suspension forces can be seen in Fig. 3.
Considering all the above-mentioned nonlinearities, the equations of motion can be written in the general form: where M(q) is the mass matrix and C(q,q) is a vector containing the remaining parts of the equation including all the nonlinear and non-smooth characteristics.

Linear stability analysis
For the linear stability analysis of the rectilinear motion of the trailer, the equations of motion are linearized around q(t) ≡ q 0 . For the sake of simplicity, we consider the case q 0 = 0, namely Due to the symmetry of our model, ϕ 0 = 0 is satisfied for any parameter setup. Moreover, we tune the free length of the spring in the suspensions to maintain zero pitch angle ϑ 0 = 0 in steady state, i.e., We assume that none of the switches in Eq. (4) occur in case of small amplitude vibrations around the rectilinear motion. This can be guaranteed by the proper choices of L min and L max . Using the perturbation y(t) = q(t) − q 0 around the rectilinear motion, the linearized equations of motion can be written as: where M lin is the mass matrix, C lin is the damping matrix and K lin is the stiffness matrix of the linearized system: where we use the mass moments of inertia about the (x, y, z) axes with respect to the king pin: In Eqs. (12) and (13), we also introduce the nominal vertical load and the tire related damping in order to shorten the formulas. Note that the stiffness matrix K lin is asymmetric due to the terms related to the lateral tire forces. See also [3] for similar asymmetry properties introduced by tire forces.
As it can be seen in Eqs. (11), (12) and (13), the linearized system can be separated into two subsystems: One differential equation can be disjointed, namely the pitch motion can be analyzed alone (1-DoF subsystem with generalized coordinate ϑ), while the remaining equations are coupled (3-DoF subsystem with generalized coordinates ψ, ϕ and u). Hence, the pitch motion does not affect the linear stability of the rectilinear motion of the trailer. In vehicle dynamics, this type of separation of the lateral (handling) and the vertical dynamics (bouncing and pitching) is commonly used for linear analysis. For the dynamics of the trailer, this approach is also confirmed by our mechanical model.
Using the exponential trial solution of the 3-DoF subsystem, one can determine the characteristic equation of the system in the form of a sixth-order polynomial with respect to the characteristic root. The stability of the rectilinear motion can be investigated by means of the Routh-Hurwitz criteria [11] based on the coefficients of the characteristic equation. On the one hand, no closed-form formula can be derived for the critical  Table 1 providing the critical towing speed 29.9 m/s marked by red dots towing speed, which is an essential parameter in stability analyses of vehicles. On the other hand, one can evaluate the Routh-Hurwitz criteria numerically for a specific set of parameters. Here, we utilize the parameter values of the experimental study of Darling et al. [7], which were also used in [3] to validate an in-plane car-trailer combination model. Some of the parameters, which are also necessary for our spatial trailer model, are not involved in the analysis [7]. Hence, we tried to identify and/or estimate these parameters based on the available information and our physical sense. Moreover, we tuned the lateral stiffness k lat and damping c lat in our model with respect to the critical towing speed and vibration frequencies detected in [7]. Note that these parameters are used to imitate the interaction between the trailer and the towing vehicle, and one can tune them based on the natural frequency of the vibration mode that is related to the stability loss of the car-trailer system. We collect all the parameters in Table 1 called as realistic trailer setup.
For this setup, our mechanical model provides a Hopf bifurcation (dynamic/oscillatory loss of stability) at the critical towing speed 29.9 m/s, which agrees with the theoretical calculations of [3]. Namely, the real part of the rightmost complex conjugate pair of characteristic roots is negative for v < 29.9 m/s, i.e., the rectilinear motion is linearly stable and the oscillations decay in time. For v > 29.9 m/s, these roots have positive real parts; accordingly, the rectilinear motion is linearly unstable and the oscillations of the trailer increase in time.  Table 1 In addition, one can also investigate the effect of different parameters on the linear stability. In Fig. 4, we construct linear stability charts in wide speed range for the stiffness k and damping c of the wheel suspension. We use these parameters because they can be varied independently from the other parameters of the trailer, and their effects are not considered in the in-plane models (like Model C). In the figure, the original setup is illustrated by dashed horizontal lines. Shaded areas correspond to linearly unstable rectilinear motion, where snaking, rocking, roll-over of the trailer may appear in practice; the white areas correspond to linearly stable rectilinear motion. The theoretical vibration frequency related to the stability boundary is not plotted in the figure, but it is in the range of 0.8 − 1.3 Hz for our system parameters. These values have good agreement with the experimental results in [7].
As it can be observed in Fig. 4, the investigated parameters influence the value of the critical speed strongly. The stiffness k of the wheel suspension has a relevant role since small variation of the original value can significantly increase the critical speed. Thus, these stability charts may suggest suitable parameter values for the trailers and can support the engineering design process.
Road accidents involving trailers are often related to inappropriate values of the mass moment of inertia and/or the vertical position of the payload. Hence, in Fig. 5, stability charts are also constructed with respect to the vertical h and longitudinal e positions of the center of gravity C. Note that the mass moments of iner-tia of the trailer may depend on these parameters, and to consider these effects we use the approximations J Cx = m(4b 2 + 4h 2 )/6, J Cy = m(l 2 + 4h 2 )/6 and J Cz = m(l 2 + 4b 2 )/6. These formulas are based on the mass and mass moment of inertia data given by Darling et al. [7], namely, J Cz provides the same value as their experimental setup and J Cx and J Cy are formulated using the same approach. As shown in Fig. 5a, the critical towing speed can be relevantly increased by increasing the height h. However, the critical speed drops again for high payload positions, what has good agreement with our physical sense. The effect of the longitudinal position e is illustrated in Fig. 5b. The plotted stability boundary suggests that the larger e is, the higher the critical speed is. Note that our model does not take into account that the so-called nose weight (the load on the towing hook) increases together with e, which may cause other stability issues. More detailed mechanical models of the car-caravan combination describe this phenomenon better and may give somewhat different results with respect to the effect of the longitudinal position e, see, e.g., [3,7].

Bifurcation analysis
In practice, the linear stability analysis given in Sect. 3 does not provide enough information about the behavior of the trailer since the nonlinear vibrations often have a key role in the dynamics. In several cases, the presence of unstable limit cycles in linearly stable speed ranges may cause safety critical large amplitude vibrations when large perturbations occur. These limit cycles can relate to subcritical Hopf bifurcations emerging at the linear stability boundaries. In this section, first, the local bifurcation of the system is investigated by means of the center manifold reduction [6,29], which is accomplished semi-analytically only due to the complexity of the system. Then, we perform numerical bifurcation analysis using continuation techniques in order to obtain information about the global dynamics of the trailer.

Local bifurcation analysis
Although the theory of the center manifold reduction [6,29] is well known, it is rarely applied to higher dimensional systems. Due to the complexity of the method, even the sense of the Hopf bifurcation cannot be determined analytically for such systems, because at a certain point of the algorithm, one should evaluate the formulas numerically. The scenario is the same for our mechanical model, too. However, as we will show, essential information is provided by the center manifold reduction with respect to the nonlinear coupling between the vertical and the lateral motions of the trailer. Such information may never be obtained if numerical continuation is applied only.
Although the algorithm of the center manifold reduction is very conventional, we provide all the details of the calculation in order to emphasize the key steps where the above-mentioned nonlinear coupling emerges and can be identified.
For the local bifurcation analysis, we rewrite the equations of motion into a system of first order differential equations. To do this, let us introduce the vector of state variables where y is the perturbation around the rectilinear motion (cf. Eq. (10)). Based on this, the governing equations (7) can be rewritten and expanded in powerseries form up to the third degree terms: where the coefficient matrix A reads The vectors F 2 and F 3 contain the second and third degree terms, respectively. Since we are interested in the local bifurcation, the non-smooth characteristics of the suspension forces (4) do not appear in Eq. (20). But it is worth to mention that our system is asymmetric even in the linear domain due to the characteristics of the lateral tire forces (see Eq. (13)). The second degree terms (F 2 = 0) are also related to the tire forces, and they make coupling between the linearly independent lateral (described by ψ, ϕ and u) and the vertical (described by ϑ) motions of the trailer. For example, the first component of the vector F 2 contains a term related to ϑψ. By means of the center manifold reduction, our eight-dimensional system is reduced to a two-dimensional one [6]. First, the normal form of the governing equations is determined by transforming the equations of motion to the basis of the eigenvectors of while the other eigenvalues are with σ j < 0 and ω j > 0. Thus, the eigenvectors corresponding to the complex eigenvalues are: where the overline refers to the complex conjugate.
Since the lateral and the vertical (pitch) motions are not coupled linearly through A (cf. Eqs. (11)(12)(13)), the eigenvectors corresponding to these linearly independent motions are perpendicular to each other. Let us denote λ 7,8 and s 7 = s 8 to be the eigenvalues and the eigenvectors of the 1-DoF subsystem of the pitch motion, respectively. Namely, s 1 · s 7 = s 3 · s 7 = s 4 · s 7 = s 5 · s 7 = 0. The columns of the transformation matrix can be composed by means of the eigenvectors (for more details, see [13]): T = Re s 1 Im s 1 s 3 s 4 Re s 5 Im s 5 Re s 7 Im s 7 , (26) which has the following structure due to the linear independence of the lateral and vertical motions: In the matrix, symbols refer to nonzero elements. By means of this transformation matrix, we introduce the new variable x = [x 1 x 2 . . . x 8 ] T as Note that after the transformation, the new variables x 7 and x 8 still describe the pitch motion of the trailer. The system in the new basis can be formulated as: where T −1 AT is the so-called Jordan matrix, which enables the separation of (29) as: where the subscript c refers to the so-called center manifold [6,29] with the center variables, while s denotes the stable variables, which are in our case: In Eqs. (30) and (31), matrices B and C are respectively. Vectors f and g contain second and third degree terms. Again, we notice that f have mixed second degree terms related to x 7 , x 8 and the center variables x 1 , x 2 . Thus, terms like x 1 x 7 , x 1 x 8 , x 2 x 7 and x 2 x 8 with nonzero coefficients occur in f. Later, we will show how the pitch motion affects the sense of the Hopf bifurcation through these mixed terms.
In order to eliminate the stable variable x s and to describe the dynamics of the system on the center manifold by means of the center variables x c , the center manifold can be approximated as a purely quadratic surface: where h k = c k20 x 2 1 + c k11 x 1 x 2 + c k02 x 2 2 for k = 3, . . . , 8 contains 18 unknown constant coefficients. We take the time derivative of Eq. (34): where D x c refers to the partial derivative with respect to x c . We substitute the expression of the quadratic surface of the center manifold in Eqs. (30) and (31); thus, we obtain the following set of equations: Let us substitute Eq. (36) in (35): We subtract Eq. (37) from Eq. (38): 2 ) are substituted in the mixed second degree terms (e.g., x 1 x 7 ) mentioned above.
In our system, after substituting the calculated coefficients into Eq. (30), namely using the second degree approximation of the center manifold, the reduced order system is given by Thus, there are no second degree terms in the nonlinear part anymore, and the Poincaré-Lyapunov parameter can be calculated via the third degree terms only, see [30]: The sense of the bifurcation can be identified by investigating the sign of the Poincaré-Lyapunov parameter Δ.
If Δ is positive, then the corresponding Hopf bifurcation is subcritical and the emerging periodic solutions are unstable. If Δ is negative, then the corresponding Hopf bifurcation is supercritical and the emerging periodic solutions are stable [30]. As mentioned before, due to the complexity of the governing equations and the large number of parameters of our mechanical model, no closed-form formula can be derived for the Poincaré-Lyapunov parameter. However, for a certain parameter setup, the sign of Δ, thus the sense of the Hopf bifurcation can be obtained semi-analytically. Here we consider our realistic trailer parameters given in Table 1 with the height h = 1 m of the center of gravity. We apply the same approximation formulas for the mass moments of inertia like in Sect. 3 to assume the dependence of these parameters on the value of h. For this setup, the above detailed center manifold reduction is applied (see Appendix A for details) providing Δ = 6.475 · 10 −4 > 0, namely the Hopf bifurcation at the linear stability boundary is subcritical. This is the more dangerous case from engineering point of view, since an unstable periodic solution coexists with the stable rectilinear motion.
If the nonlinear effect of the pitch motion is neglected, i.e., the nonlinear coupling between the vertical (pitch) and the lateral motions is omitted by considering h 7 = h 8 = 0 (see Appendix B), then we obtain Δ = −8.031 · 10 −4 < 0. Therefore, the Hopf bifurcation is supercritical, which suggests that using models neglecting pitch motion may not identify the presence of the more dangerous subcritical Hopf bifurcation case.
In Fig. 6, stability charts are constructed in the plane of the towing speed v and the height h of the center of gravity. The subcritical and supercritical Hopf bifurcations are illustrated by dashed red and solid blue curves, respectively. The sense of the Hopf bifurcation is determined by means of numerical continuation, detailed later in the next subsection. Panel (a) corresponds to our spatial trailer model (Model A), while panel (b) is related to the simplified model (Model B) in which the pitch motion is blocked. As shown in Sect. 3, the linear stability boundaries coincide for both models. However, the sense of the Hopf bifurcation disagrees for certain parameter setups. In the figures, we marked the original h = 0.21 m and the modified h = 1 m trailer setups with dashed horizontal lines. For h = 1 m, subcritical and supercritical bifurcations take place for Model A and Model B, respectively. Thus, the numerical continuation confirms the results of the semi-analytical bifurcation analysis.
In Fig. 6c, the stability chart of the in-plane model (Model C) is constructed, which does not describe the dependence of the critical towing speed on the vertical position of the center of gravity. By the way, this model provides supercritical Hopf bifurcation at a smaller critical speed, namely, v = 23.8 m/s.

Numerical bifurcation analysis
In order to validate the analytical calculations, numerical bifurcation analysis is performed with the help of DDE Biftool, see [10]. Namely, the Hopf bifurcations corresponding to the linear stability boundaries are located and the branches of the emerging periodic solutions are followed. As it was already shown in Fig. 6, the sense of the Hopf bifurcation is determined by this methodology.
Of course, numerical continuation is not limited to local bifurcation analysis, but it can also provide information about the large amplitude vibrations of the trailer. However, for this, one has to consider the nonsmooth nature of the governing equations detailed in Sect. 2.2.
We approximate the Heaviside function bỹ where ε is the so-called smoothing parameter. The smaller ε is, the more sudden the switching is. Using In order to tend to the real scenario, the limit case ε d → 0 and ε F → 0 is investigated by choosing smaller and smaller ε d and ε F values as long as the numerical stability of the computation is ensured. Accordingly, parameters ε d = 10 −3 m and ε F = 10 −3 N are used in our calculations, which leads to sudden changes within one time step of the detected periodic solutions related to the suspension forces (cf. Eq. (4)). We also take care how the number of mesh points and the degree of the interpolation function are selected in the collocation method. Namely, we increased the number of mesh points to 144, while linear interpolation was applied. This set of the parameters ensures an adequate convergence with acceptable computational time.
A nonlinear stability chart is constructed in Fig. 7, where the domains of different colors refer to different nonlinear properties of the trailer. In the white area, the rectilinear motion is globally stable. In the light and dark gray domains, the rectilinear motion is linearly unstable and a stable limit cycle takes place in the state space. We differentiate the limit cycles with respect to the minimum suspension forces (minimum vertical force on the tires) corresponding to them. In the dark gray zones, this minimum force is zero, namely the tires lose contact with the ground. The red unsafe zone refers to the parameter domain where the linearly stable rectilinear motion is not globally stable due to the presence of an unstable limit cycle. Namely, large enough perturbations may lead to large amplitude vibrations of the trailer. Such unsafe zones typically appear at subcritical Hopf bifurcations, but saddle-node bifurcations of periodic orbits can also provide such zones at supercritical Hopf bifurcations. This latter case also happens for our system parameters, see the figure.
To present more information about the nonlinear properties of the system, bifurcation diagrams related to dashed horizontal lines of Fig. 7 are plotted in Fig. 8. The amplitudes of the generalized coordinates ψ, ϑ, ϕ and u are plotted as functions of the towing speed v. Dashed red lines and solid blue lines refer to unstable and stable motions, respectively. Red dots denote the critical speeds where Hopf bifurcations take place. Black dots characterize the towing speed above which the loss of contact of the tires happen. This latter type of limit cycles is marked with thick blue lines. Based on the continuation, one can observe that symmetric solutions exist for the lateral motion (ψ, ϕ, u), while the pitch motion (ϑ) is asymmetric. Thus, both the min/max values of the periodic solutions are illustrated for ϑ. Figure 8a shows the results of the original parameters of [7], providing supercritical Hopf bifurcation without unsafe zone. More precisely, saddle-node bifurcations of periodic orbits are within the linearly unstable speed range. At the first sharp folding point on the periodic bifurcation branch, the full compression of the wheel suspension appears leading to a qualitative change of the motion of the trailer thanks to the non-smooth nature of the suspension force.
In Fig. 8b, we present an example for the case where supercritical Hopf bifurcation and unsafe zone coexist. In Fig. 8c, a simpler structure of the limit cycles can be  Figure 8d relates to the parameter setup for which we presented the results of the semi-analytical local bifurcation analysis. The sense of the Hopf bifurcation corresponding to the linear stability boundary is subcritical. As suggested by the small absolute value of the calculated Poincaré-Lyapunov parameter Δ = 6.475 · 10 −4 , the numerical continuation provides a steeply rising branch of unstable periodic solutions. The saddlenode bifurcation of periodic solutions emerges resulting in a narrow bistable speed region, which is the reason why no red unsafe zone can be observed in Fig. 7 for the larger values of h. Although our spatial trailer model (Model A) shows subcritical Hopf bifurcation instead of the supercritical one (related to Model B), this phenomenon may not have any practical relevance.
As shown by the bifurcation diagrams presented here, the non-smooth characteristic of the suspension forces influences relevantly the nonlinear behavior of the system. The full compression/expansion of the Fig. 9 The small-scale experimental setup wheel suspension and the loss of contact of tires occur for all investigated parameters. We will show more details about these scenarios in Sect. 5 with respect to our experiments. Fig. 10 The normalized tire force characteristic (black curve) μ(α) plotted for positive side slip angle α (see Eq. (3)). Blue dots refer to the mean value of measurement data. The standard deviation of the measurement points is marked by blue bars belt, whose speed can be tuned between 0 m/s and 5 m/s. The motion of the king pin is constrained in the towing direction by a roller bearing linear guide, while the king pin can move in the lateral direction against the spring and damper imitated by a rubber band in our experiment. In order to be able to compare the experimental results with the theoretical ones, the spring stiffness and the damping of the rubber band were measured. The position of the center of mass and the mass moment of inertia of the trailer can be tuned by the positioning of the dummy payloads. The caster length l can also be modified by shifting the wheel suspension along the longitudinal direction.
Measurements were carried out to determine the tire force characteristics of the experimental rig. A pair of wheels were towed on the conveyor belt by a rigid caster. The side slip angle α was measured, while a constant lateral force was applied to the wheels. The measurement points normalized by the vertical load are shown in Fig. 10. The parameters B, C, D and E of the Magic Formula (3) were identified via the fitting of the formula to the measurement data.
Since the mass moment of inertia depends on the position of the payload, we fixed the parameters e and h during the experiments, as well as the caster length l. Thus, we set different towing speed values for the conveyor belt and examined the linear stability and the nonlinear vibrations of the trailer. We measured the lateral displacement of the king pin with a laser sensor and collected and evaluated the data in Matlab.
Measurement data are compared to theoretical results in Fig. 11 for the parameter setup of the experimental rig (see Table 1). In panel (a) and (b), the measured amplitude u max of the lateral displacement and the frequencies f of detected vibrations are plotted versus the towing speed v, respectively. Measurement points denoted by black crosses are compared to the results of the numerical continuation shown by blue curves. Relatively good agreement can be observed for the amplitudes; however, some shift occurs with respect to the linear stability boundary. In addition, there is somewhat larger difference between the predicted and the measured frequencies.
In Fig. 11c, a stability chart is constructed by numerical continuation for the experimental setup in the plane of the towing speed v and the damping c of the wheel suspension. As shown, there is a domain where the  Fig. 11c loss of contact of the tires occurs. The experimental setup is marked with a horizontal dashed line, for which the minimum and maximum values of suspension distance d (see Fig. 3) are also plotted in panel (d), highlighting that neither full compression nor full expansion of the suspension happen for this setup. The minimum and maximum values of the suspension force F s are plotted in panel (e).
In Fig. 12, we present some numerical results for the parameter points M 1 , M 2 and P of Fig. 11. We plot the periodic solutions for all four generalized coordinates (ψ, ϑ, ϕ and u) by solid blue lines for one period T of the oscillation. The suspension displacements d and suspension forces F s are also illustrated with thin black and thick gray curves corresponding to the right and left wheels, respectively.
The measured time signals of the lateral displacement u at parameter points M 1 and M 2 (see Fig. 12ab) are also plotted by dashed black curves. In order to compare the shapes of the theoretical and the measured periodic motions meanwhile we neglect the dif-ference of the vibration frequencies (cf. Fig. 11b), we plot one period of the measured signal, too. A good agreement can be observed between the shapes of the measured and the theoretically calculated time histories. As mentioned before, a loss of contact happens for the parameter point M 2 , but the amplitude of the vibrations remains almost the same as for M 1 , no large amplitude, so-called rocking motion of the trailer occurs. This was also the case in the experiment, namely no gap between the wheels and the conveyor belt was visible. The wheel suspension was neither fully compressed nor fully expanded as suggested by the theoretical results, see the evolution of the suspension displacement in the figures.
Point P is in the parameter domain of Fig. 11 where full compression and full expansion of the wheel suspensions occur together with the loss of contact of the tires. The corresponding time signals are shown in Fig. 12c. In this case, the vibration amplitudes are significantly larger than in the former cases, namely a rocking motion may occur in practice.

Conclusions
In this study, we performed linear and nonlinear analyses of three possible mechanical models of trailers, see Fig. 1. The in-plane model, which is commonly used in vehicle dynamics, neglects the pitch and roll motions of the trailer, and it provides a smaller critical speed and supercritical Hopf bifurcation for the investigated parameter setup. Model B, the reduced model with blocked bounce/pitch motion, yields the same linear stability boundaries as the full spatial model (Model A). However, the linearly independent pitch motion affects the sense of the Hopf bifurcation at the linear stability boundary. The full spatial model provides subcritical Hopf bifurcation for higher center of gravity positions, but this does not lead to the occurrence of relevant unsafe parameter domains. Nevertheless, as the critical speed is crossed, large amplitude vibrations can suddenly appear, which means a relevant safety risk on roads that can be identified only by the use of the full spatial trailer model.
The different types of motions of the trailer were also explored by our study. Domains of loss of contact of tires were shown for the realistic trailer setup. An unsafe zone with practical relevance was detected at lower center of gravity positions, which significantly reduces the speed range, where the rectilinear motion is globally stable.
Laboratory experiments were used to validate the theoretical results. The numerically determined stable limit cycles were compared to the measured vibration amplitudes and frequencies, and a good qualitative agreement was established.
Funding Open access funding provided by Budapest University of Technology and Economics. This research was partly supported by the National Research, Development and Innovation Office under Grant No. NKFI-128422 and by the NRDI Fund (TKP2020 IES, Grant No. BME-IE-MIFM and TKP2020 NC, Grant No. BME-NC) based on the charter of bolster issued by the NRDI Office under the auspices of the Ministry for Innovation and Technology.

Declarations
Conflict of interest 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/.

Appendix A
By considering our realistic trailer parameters given in Table 1 with the height h = 1 m of the center of gravity and applying the same approximation formulas for the mass moments of inertia like in Sect. 3, the coefficient matrix A of Eq.
The value of the Poincaré-Lyapunov parameter is Δ = 6.475 · 10 −4 > 0, and the Hopf bifurcation at the linear stability boundary is subcritical.
The value of the Poincaré-Lyapunov parameter is Δ = −8.031 · 10 −4 < 0, and the Hopf bifurcation at the linear stability boundary is supercritical.