Stability issues of two-wheeled trailers caused by the nonlinear coupling between the 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 and other geometrical nonlinearities are taken into account. Beyond the linear stability analysis, the nonlinear vibrations are analyzed with special respect 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, thus, the sense of the Hopf bifurcations are veriﬁed semi-analytically and numerically. Simpliﬁed models of the trailer are also used in order to point out the practical relavance of the study. Experiments are carried out on a small-scale trailer to validate the theoretical results.

plored 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 utilizied in the industry too. Stability charts, phase portrait and nonlinear bifurcation diagrams were composed for the most intricate nonlinear problems of vehicle dynamics, like the the shimmy motions of steered wheels (see [2][3][4]19,22,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 maneuovres are also in the focus of recent studies (see [12,24]), which allow the developement 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 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 so-called 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 cartrailer 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 is used, where all the yaw, pitch and roll motions are considered, see Model A in Fig. 1(a). In addition, the wheel suspensions are also modeled in our analysis. In order to have a relatively low degree-offreedom 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 in-plane 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 Sec. 2. In Sec. 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  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 happen 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 Sec. 5. Finally, our theoretical, numerical and exerimental results are concluded in Sec. 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 Fig. 2. 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 l and by the lateral damping c l . 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 corresponding mathematical formulas on a manageable level.
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. 2(a). In order to avoid the presense 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 [22]. 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 rigour. 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. 2(b)) 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 l = −k l u − c lu acting at point A due to the deformation of the lateral spring and damper.
The lateral tire forces are calculated based on Pacejka's Magic Formula [22] 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 µ(α) is where B is the stiffness factor, C is the shape factor, D is the peak factor and E is the curvature factor. Although the variation of the camber angle, which obviously occurs due to the rolling motion of the trailer, can be taken into account via the variation of these factors, we neglect this effect in our investigation.
On the contrary, we pay special attention to the rocking motion of the trailer, namely, we take into account that the tires can detach 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. 2(b). To handle the non-smooth characteristic of the suspension force, we use the so-called Heaviside step-function H(s) [8], often defined as H(s) = s −∞ δ(s)ds, where δ(s) is the Diracdelta 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 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 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 vehile 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 subsytem, one can determine the characteristic equation of the system in the form of a 6th 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 towing speed, which is an essential parameter in stability analyses of vehicles. On the other hand, one can eveluate 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 l and damping c l in our model with respect to the critical towing speed detected in [7]. Note that these parameters are used to imitate the interaction between the trailer and the towing vehicle, and their values could be identified more precisely by means of specific experiments. We collect all the parameters in Tab. 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.
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 inertia of the trailer may depend on these paremeters, and to consider these effects we use the approximatations 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. 5(a), 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. 5(b). 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 mechan-ical models of the car-caravan combination describe this phenomenon better and may give somewhat different results with respect to the effect of e, see, e.g. [3,7].

Bifurcation analysis
In practice, the linear stability analysis given in Sec. 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 power-series 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. (19). 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 eightdimensional 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 the linear part. Let the eigenvalues and the eigenvectors of A be denoted by λ j and s j , respectively. Here, we investigate the case when six eigenvalues have non-zero imaginary parts and two real eigenvalues exist, which case holds for the parameter setup analyzed in this section. At the linear stability boundary, a complex conjugate pair of eigenvalues exists with zero real part (see Fig. 5): 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 , (25) which has the following structure due to the linear independence of the lateral and vertical motions: In the matrix, symbols refer to non-zero 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 aṡ where T −1 AT is the so-called Jordan matrix, which enables the separation of (28) aṡ 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. (29) and (30), 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 non-zero 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. (33): where D xc refers to the partial derivative with respect to x c . We substitute the expression of the quadratic surface of the center manifold into Eqs. (29) and (30), thus, we obtain the following set of equations: Let us substitute Eq. (35) into (34): We substract Eq. (36) from Eq. (37): are substituted into the mixed second degree terms (e.g. x 1 x 7 ) mentioned above.
In our system, after substituting the calculated coefficients into Eq. (29), 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é-Ljapunov 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é-Ljapunov 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é-Ljapunov 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 Tab. 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 Sec. 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 presense 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 Sec. 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. 6(c), 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 Sec. 2.2.
We approximate the Heaviside function by where ε is the so-called smoothing parameter. The smaller ε is, the more sudden the switching is. Using this function, the ReLU(s) in Eq. (4) can be also formulated as ReLU(s) = H(s) · s. Thus, the originally non-smooth governing equations are approximated by their smoothed version during the continuation. In order to tend to the real scenario, the limit case ε → 0 is investigated by choosing smaller and smaller ε values as long as the numerical stability of the computation is ensured. Accordingly, parameter ε = 10 −3 is 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 presense 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 8(a) 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 nonsmooth nature of the suspension force.
In Fig. 8(b), we present an example for the case where supercritical Hopf bifurcation and unsafe zone coexist. In Fig. 8(c), a simpler structure of the limit cycles can be observed, namely, subcritical Hopf bifurcation leads to the unsafe zone. Figure 8(d) 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é-Ljapunov 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   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 wheel suspension and the loss of contact of tires occur for all investigated parameters. We will show more details about these scenarios in Sec. 5 with respect to our experiments.

Experiments
Based on the mechanical model described in Sec. 2, a small-scale experimental rig was designed and manufactured, see Fig. 9. The trailer is placed on a conveyor 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. 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. The parameters B, C, D and E of the Magic Formula (3) were identified via the fitting of the formula to measurement data, see Fig. 10.
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  Fig. 10: The fitted 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.
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 Tab. 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 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. 11(c), 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 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. 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. 12(ab)) are also plotted by dashed black curves. A good stable (stable rectilinear motion) Fig. 11: Experimental results. The detected amplitude u max (a) and freqency f (b) of the vibrations versus the towing speed v. Black crosses refer to measurement data, blue curves represent numerically determined stable solutions. Thick segments refer to motions with loss of contact of tires. Stability chart (c) for the damping c of the suspension; the dashed horizontal line corresponds to the experimental setup. Numerically determined min and max values of the suspension displacements d (d) and the suspension force F s (e). 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 re- sults, 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. 12(c). 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 occurance 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 were 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 limitcycles were compared to the measured vibration amplitudes and frequencies., and a good qualitative agreement was established.
By considering our realistic trailer parameters given in Tab. 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 Sec. 3 therefore, the coefficients of the third degree terms of the reduced system are The value of the Poincaré-Ljapunov parameter is ∆ = 6.475 · 10 −4 > 0 and the Hopf bifurcation at the linear stability boundary is subcritical.