Dynamic model of a logistic train with different steering systems and tire models

Abstract

Logistic train is one of the recent solutions applied in the intralogistics. It bases on the concept of the “milk-run” delivery method. The main problem in the implementation of intralogistics trains is the possibility to determine whether the train is able to pass the given path without a colliding with surrounding objects. Dynamic model presented in this paper was performed for the three most common steering systems for the trolleys: double Ackermann steering system and two drawbar systems: conventional and virtual clutch. The dynamic model was developed following Lagrange’s theorem including the possibility of lateral slip. The system of differential equations was numerically solved. In order to calculate the lateral friction forces we used three different tire models: simplified, sigmoidal and Pacejka Magic formula. The results obtained for different tire models and steering systems are presented in the form of animations presenting train run in various conditions. The analyses have shown that the use of simplified tire models is justified under certain conditions.

Keywords:
logistic train; multibody; dynamic model; friction; vehicle; slip; steering system; double Ackermann; drawbar; trailer; trolley; tractor; tire model

Graphical Abstract

1 INTRODUCTION

The main goal of this paper is to develop a credible dynamic model of the multiple trailers on a tractor system for the production logistics application i.e. logistic train. This includes the most common steering system: virtual clutch and drawbar system, conventional clutch and drawbar system as well as double Ackermann steering. The objective of this paper is also to study the influence of different tire models and steering systems of trolleys on the slip effect. This model should contribute to a better understanding of an impact of geometrical relations and friction on the trajectory of such systems in motion. This can facilitate the logistic train design for improved maneuverability in narrow corridors through choice of applicable number of trailers and selections of the best geometric parameters of the trolley. It can also help to estimate collision risk with nearby static objects.

1.1 Steering systems for the logistic train

Logistic trains are applied means of transport across all sectors of production supply, especially in the intralogistics system (Lieb et al., 2019). A logistic train is usually called a milk-runner. The milk-run concept is derived from the method of delivering or receiving supplies in the dairy industry, where one tanker collects milk from many farmers. In the industry, one train supplies multiple workstations. The idea is to regularly visit assigned locations by one supplier in a single run (Kluska and Pawlewski, 2018). Milk run system minimizes the total distance travelled, as a single vehicle is use instead of many (Knez and Gajsek, 2015), especially in case of repeatable flow of materials between the same locations. The specificity of the train is that it transports more goods than other intralogistics means of transport per each run (Alnahhal et al., 2014).

The discussed train consists of a tractor and a certain number of towed trailers also called logistics trolleys. The tractor is usually equipped with Ackermann steering system. The most common steering systems for trolleys include: a conventional or virtual clutch as a drawbar system and double Ackermann steering system (Paszkowiak et al., 2019). The main differences in these systems are presented in Figure 1.

Figure 1
A comparison of three kinematics of steering systems for logistic trolleys.

The classic solution for connecting the trolleys is a long drawbar at the front the unit. This conventional connection could be potentially ineffective as it was showed that another two types of steering systems can perform better when passing in narrow corridors. The virtual clutch has a longer drawbar at the rear of the trolley. Rear wheels are fixed, whereas front ones are castor wheels. A more complex system is a double Ackermann steering, in which a number of degrees of freedom is higher than for others, assuming the same number of trolleys. This is due to the presence of a rotary drawbar, which turn determine a turn front wheels of the trolley. Rear wheels rotate in the opposite direction than front, which direction of turning is consistent with drawbar. A graphical comparison of kinematics of these steering system was presented in Figure 2 (Paszkowiak et al., 2019).

Figure 2
A comparison of the kinematics of the three steering systems.

1.2. Vehicle dynamics modelling

Vehicle dynamics are designed based on the assumption that a rigid vehicle behaves similar to a flat object moving on a plane. The modeled vehicle moves in a planar motion with three degrees of freedom (DOFs) (Jazar, 2009). Dynamic models of vehicles can be expressed based on Lagrange’s theorem or Newton’s second law (Tzafestas, 2014; Ahmad Abu Hatab, 2013). These models are applicable in numerous fields including: mobile robots (Keymasi Khalaji et al., 2014), cars (Altche et al., 2017), trucks (Adamiec-Wójcik et al., 2019) agricultural vehicles (Kayacan et al., 2013), airplanes (Gao and Jia, 2017), accident reconstruction (Brach and Brach, 2009), autonomous vehicles (Kong et al., 2015) and racing games (Hulme et al., 2009). The most common models involve nonholonomic constraints (Bloch et al., 2015). The idea in those approaches is based on Chaplygin sleigh model, where lateral friction forces and lateral velocity are not considered for wheels (Bravo-Doddoli and García-Naranjo, 2015). This assumption causes the trajectory of vehicle to be known and dependent on the kinematics, the drive power and the mass of the vehicle. In those models, it is assumed that the friction is not a relevant factor. Nonholonomic constraints appear with the assumption of a non-slip condition. Center points of the wheels cannot move along the direction of the wheel own axis of rotation (Lucet et al., 2008).

The solution that takes into account the slip effect involves the use of the holonomic system of a Lagrange formula or equations based on Newton’s second law. This is also possible when using the model for a non-holonomic system. In that case, forces acting on the system, as a result of lateral friction force, are expressed as Lagrangian multipliers. The torque is included in the input vector (Piyabongkarn et al., 2009). Many authors assume that the lateral friction force acting on a tire is proportional to the angle of slip. This assumption is valid for small values of side slip angles. They must be smaller than the values of slip angle corresponding to the peak lateral force (Kayacan et al., 2014). It allows determination of the slip angles without using trigonometric functions (Martinez et al., 2004). However, this assumption has some limitations. The model can represent a drive only in the smooth passing. It does not involve drifting or aggressive race (Li et al., 2014). The side-slip angles are determined by analyzing the constituents of the longitudinal and lateral velocity vectors of the wheel (Andrzejewski and Awrejcewicz, 2005).

In order to consider a non-linear slip, a tire model can be used. This model involves friction forces using a coefficient of friction and a normal force acting on a wheel. It also takes into account parameters of a tire, road and a condition of the slip between a tire and floor. Modeling of the tire can be analytical, numerical or parameterized. The analytical method is based on physical parameters and is simple to calculate. Nonetheless, such models do not accurately reflect the actual behavior of the tire (Mondek and Hromcik, 2017). Examples of these are: linear, Fiala and Brush models (Svendenius and Wittenmark, 2003). The simplest model of the tire exists in two versions. One of them is expressed as a sine function. In that case, the maximum of lateral friction force is for side slip angle equal to 90 degrees (Paszkowiak and Bartkowiak, 2019). However, for most tires such relationship does not exist in practice. Thus, a linear model in which a maximum of a lateral friction force occurring for smaller slip angle, is used more often (Mondek and Hromcik, 2017). This happens for side slip angle reaching around 4 degrees. For greater value of this angle the friction force remains constant (Heisler, 2002). Numerical methods use look up tables and there is no explicit mathematical form of equations. In these models, geometry and material of tire are considered. Numerical models derived from the detailed solution of finite difference and finite element multidimensional models of a tire and are represented as a set of tables. These models are difficult to be applied for model-based control development (Baranowski et al., 2012). Parametrized models are the most commonly used for vehicle dynamics simulation studies and the control design. They require experiments for a specific tire in order to estimate their parameters via statistical fitting. A notable example is Pacejka Magic formula, developed by Hans Bastian Pacejka (Pacejka et al., 2012). This tire model requires many experimental parameters and coefficients, which are complicated to obtain (Ortiz et al., 2006; Cabrera et al., 2018).

There are a multitude of tire models, each of them having advantages and disadvantages, depending on the specific application. Gang Liu et al. (2013) compared the simplified piecewise linear tire model with the Magic formula on the example of a bicycle model representing a car model. This comparison, in which vehicle trajectories and speed were compared, proved that the simplified model could be as accurate as the bicycle model with the MF. Comparisons of tire models involve fitting the model to the actual curve, represented by the Magic formula. Similar approach was shown by Lu et al. (2019). They compared a standard and modified LuGre tire models and Pacejka Magic formula. LuGre is one of the most popular analytical tire model (Liang et al., 2008). It can accurately describe the stick-slip motion and frictional hysteresis. They also showed that the modified LuGre dynamic tire model could also be applicable for analysis of vehicle system (Lu et al., 2019). All these aforementioned works show that new tire models, with better applicability, are still being sought.

Due to the fact that the logistic train is a new solution, its dynamics were not considered taking into account an accurate model of the tire. Paszkowiak and Bartkowiak (2019), presented a brief comparison of dynamics of a logistic train with different steering systems. They analyzed slip effect using only visualization of trajectory. Input parameters and conditions assumed heavily loaded trailers. Therefore, a very small friction coefficient was considered in order to trigger slip effect (Paszkowiak and Bartkowiak, 2019).

In this paper, the authors would like to present a dynamic model of a multibody system, i.e. a tractor towing multiple trailers, which can allow comparison of different steering systems with different tire models. The aim is to conduct the slip analysis using visualization and output parameters, which can estimate the intensity of the slip. This can facilitate the indication the steering system which is the most and least susceptible to slipping and also determine whether simplified tire models provide similar results to more complex approaches. The aim is also to show an alternative way to determine the side-slip angle without analyzing the velocity vectors. This is possible by calculating the time derivative of the wheel position, what simplifies the determination of the slip angle. At an initial stage of this research, the results of the new and classical methods were calculated. No significant differences were observed in those pretests. This testifies that use of the proposed approach is justifiable.

2 DYNAMIC MODEL OF A LOGISTIC TRAIN

2.1 Equation of motion

For each of three steering systems the main equation of motion remains the same. Due to the identical kinematic model of conventional and virtual clutch they also have the same dynamic model. The main equation of motion for holonomic system by the Lagrange formula can be given as

$M\left(q\right)\ddot{q}+C\left(q,\dot{q}\right)=Q\left(q\right)$ (1)

The generalized coordinates vector $q$ is a column matrix and it contains n number of elements in accordance to the number of DOFs. This vector can be presented as

$q={\left[\begin{array}{cccc}{q}_1& q_2& \cdots & q_n\end{array}\right]}^T$ (2)

Left side of main equation (1) is expressed as column matrix

$\zeta \left(q,\dot{q}\right)={\left[\begin{array}{cccc}{\zeta }_1& {\zeta }_2& \cdots & {\zeta }_n\end{array}\right]}^T$ (3)

According to Lagrange formula it can be presented as

$\zeta \left(q,\dot{q}\right)=\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{q}}\right)-\frac{\partial \mathcal{L}}{\partial q}$ (4)

where $\mathcal{L}$ is a Lagrangian of the system. For our case, the potential energy is null. It causes Lagrangian $\mathcal{L}$to be equal to the kinetic energy of the system $T$. System inertia matrix can be expressed as

$M\left(q\right)=\left[\begin{array}{cccc}\frac{\partial {\zeta }_1}{\partial {\ddot{q}}_1}& \frac{\partial {\zeta }_1}{\partial {\ddot{q}}_2}& \cdots & \frac{\partial {\zeta }_1}{\partial {\ddot{q}}_n}\\ \frac{\partial {\zeta }_2}{\partial {\ddot{q}}_1}& \frac{\partial {\zeta }_2}{\partial {\ddot{q}}_2}& \cdots & \frac{\partial {\zeta }_2}{\partial {\ddot{q}}_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial {\zeta }_n}{\partial {\ddot{q}}_1}& \frac{\partial {\zeta }_n}{\partial {\ddot{q}}_2}& \cdots & \frac{\partial {\zeta }_n}{\partial {\ddot{q}}_n}\end{array}\right]$ (5)

The centripetal and Coriolis matrix can be denoted as

$C\left(q,\dot{q}\right)=\zeta \left(q,\dot{q}\right)-M\left(q\right)\ddot{q}$ (6)

2.2 Friction forces

The tire model in the dynamic model is included as a function of lateral friction. The necessary parameters for determining the forces include load force, friction coefficient and sideslip angle. This is performed for each wheel separately. The idealized model of the wheel is shown in Figure 3. The lateral friction force is calculated based on the load force, the slip angle and the coefficient of friction between floor and wheel. The slip angle depends on the direction of the linear velocity of the wheel center $V_s$. For an empirical model, such as Pacejka Magic formula, there are more factors, but in the presented case they have a constant value.

Figure 3
An idealized model of a wheel.

2.2.1 Simplified model

Lateral friction forces acting on wheels cannot be treated as a product of a lateral friction coefficient and a load force. The reason for this are the specific properties of the tire, mainly shape and material. As a consequence, lateral stiffness changes as a function of the slip angle. Nevertheless, the simplified tire model presented in this section contains such a relationship between lateral friction forces and side-slip angle. Here, the tire is treated as solid. Lateral friction force acting on the wheel can be described as

$F_y\left(F_n,\mu ,\alpha \right)=F_n\mu \sin \left(\alpha \right)$ (7)

where $F_n$, $\mu$, $\alpha$ refer to the individual wheel and they are respectively: lateral friction force, friction coefficient and side slip angle at the moment. In this model, considering material and geometric characteristics of a specific tire might be limited as only manipulations of friction coefficient are possible. This makes this model potentially accurate in a narrow range of a side slip angle. However, because of its undeniable simplicity, it was decided to have it included in this study.

2.2.2 Sigmoid function

Experimental studies presented by Cabrera et al. (2018) and Baffet et al. (2006) showed that the relationship between friction force and slip angle was non-linear. The function consists of three distinct phases of vehicle handling behaviors, such as: linear, transient and saturation. In the linear phase, the lateral friction forces are linear relative to the slip angle. This occurs for small values of the sideslip angle, usually less than 4 degrees (Heisler, 2002). In the next phase, the transient function is non-linear and the force reaches its maximum. The last phase is a saturation, in which friction force decreases. The second and third phases occur when the vehicle approaches the physical limit of adhesion (Baffet et al., 2006). In our sigmoid tire model, we adopted similar stages with one exception. When the force reaches its peak, it does not change as the angle of slip increases (Figure 4). The sigmoid function is used instead of the piecewise linear function. The piecewise function can cause problems when solving numerically. A sudden change in the partial function can complicate the solution of a differential equation for some special conditions. Therefore, the continuous function is used as a sigmoidal function. In this case, the function increases linearly in the first stage. When the slip angle exceeds 8$°$, the friction force begins to smoothly reach its maximum value. For greater values of the slip angle, the friction does not change. The friction force as a sigmoid function can be expressed as

$F_y\left(F_n,\mu ,\alpha \right)=\frac{2F_n\mu }{1+e^{-g_s\sin \left(\alpha \right)}}-F_n\mu$ (8)

where $g_s$ is parameter of sigmoid function and can be obtained from below equation by numerical solving

$\frac{2}{1+e^{-g_s\sin \left({\alpha }_{lin}\right)}}-1=\frac{1}{1.1}$ (9)

where ${\alpha }_{lin}$ is an indicative value of angle, where the linear relationship ends. In presented model, ${\alpha }_{lin}$ is equal to $\frac{\pi }{15}$.

2.2.3 Pacejka Magic formula

The one of the most popular way to determine friction forces is Pacejka Magic formula, also known as MF. Lateral friction force is a function of normal force, friction coefficient and slip-side angle for an analyzed wheel. That function consist of numerous parameters, which make the model more realistic. These parameters describe the tire and the conditions in which it operates. This relation can be described as a function of load, friction coefficient and slip angle in the form

$F_y\left(F_n,\mu ,\alpha \right)=p_D\sin \left[p_{C}arctg\left\{p_B\left(\alpha \frac{180}{\pi }+p_H\right)-p_E\left(p_B\left(\alpha \frac{180}{\pi }+p_H\right)-arctg\left(p_B\left(\alpha \frac{180}{\pi }+p_H\right)\right)\right)\right\}\right]+p_V$ (10)

where the coefficient and tire/road parameters are described in ^{Appendix A}.

2.2.4 Comparison of tire models

The presented tire models are compared graphically in Figure 4. To visualize the fundamental differences between those approach, an arbitrary load of 250 N was assumed which acted vertically on a tire. The lateral friction force for MF increases steeply with angle for $\mu =0.1;0.5;1$ (Figures 4a, 4b, 4c). For the highest coefficient of friction, it grows more gently than the sigmoid function (Figure 4d). Only the Pacejka Magic formula function changes the shape of the curve depending on the coefficient of friction. The force increases slower with the angle in a simplified tire model than the other two. Saturation effect is included in MF only. The geometric and material parameters in Pacejka model were assumed based on Pacejka et al. (2012) and they do not reflect any specific tire. This paper does not focus on the characteristics of tires but rather on the dynamic model of logistic train and effects of different friction model of tires. We did not want to include any specific tire example in order to make this case more generalized.

Figure 4
Comparison tire model curves for parameters from ^{Appendix A} and also $F_n=250N$ and lateral friction coefficient equal to respectively: a) $\mu =0.1$, b) $\mu =0.5$, c) $\mu =1$, d) $\mu =10$. Please note that the vertical scale is different in each graph.

2.3 Dynamic model of a tractor

In our model, the tractor is treated as a tricycle, where a single front wheel is powered (Yavin, 2007). The input parameters for this system are: drive torque $\tau$ and steering angle of the front wheel ${\psi }_{f0}$. The notation of parameters is given in ^{Appendix B}. Basic geometrical parameters are shown in Figure 5a.

For the tractor, generalized coordinate vector contains three constituents: coordinates of $B_0$ point and the angle between main axis of tractor relative to the horizontal axis. The steering angle is an input parameter. Therefore, it is not included in this vector. Vector can be described as: $q={\left[\begin{array}{ccc}x& y& {\beta }_0\end{array}\right]}^T$.

Figure 5
a) geometrical parameters, b) friction forces and slip-side angles for a tractor.

The configuration of the tractor can be expressed with generalized coordinate vector. In order determine location of relevant points tractor dimensions are needed. Table 1 presents position equations for the most important points of tractor. $P_0$ is a center of mass.

Table 1
Position equation of point of the tractor.

Point	Coordinate	Expression
$B_0$	$x_{B0}$	$x$
	$y_{B0}$	$y$
$P_0$	$x_{P0}$	$x+c_0\cos \left({\beta }_0\right)$
	$y_{P0}$	$y-c_0\sin \left({\beta }_0\right)$
$A_0$	$x_{A0}$	$x+h_0\cos \left({\beta }_0\right)$
	$y_{A0}$	$y-h_0\sin \left({\beta }_0\right)$
$D_0$	$x_{D0}$	$x-d_0\cos \left({\beta }_0\right)$
	$y_{D0}$	$y+d_0\sin \left({\beta }_0\right)$
$E_0$	$x_{E0}$	$x+\frac{1}{2}{v}_0\sin \left({\beta }_0\right)$
	$y_{E0}$	$y+\frac{1}{2}{v}_0\cos \left({\beta }_0\right)$
$F_0$	$x_{F0}$	$x-\frac{1}{2}{v}_0\sin \left({\beta }_0\right)$
	$y_{F0}$	$y-\frac{1}{2}{v}_0\cos \left({\beta }_0\right)$

The kinetic energy of the tractor $T_b$ can be denoted as

$T_b=\frac{m_b}{2}\left({\dot{x}}_{P0}^2+{\dot{y}}_{P0}^2\right)+\frac{1}{2}{I}_b{\dot{{\beta }_0}}^2+m_{w0}\left({\dot{x}}_{}^2+{\dot{y}}_{}^2\right)+\left(I_{w0}+m_{w0}{\left(\frac{v_0}{2}\right)}^2\right){\dot{{\beta }_0}}^2+\frac{m_{w0}}{2}\left({\dot{x}}_{A0}^2+{\dot{y}}_{A0}^2\right)+I_{w0}{\dot{{\beta }_0}}^2$ (11)

The side-slip angle is obtained for each wheel separately. It is the ratio of longitudinal and lateral velocity of the wheel. The idea of marking of the side slip angles is presented in Figure 5b. The side-slip angles for a rear wheels ${\alpha }_{rL0}$, ${\alpha }_{rP0}$ and the front wheel ${\alpha }_{f0}$ are expressed in the arctg function as

${\alpha }_{rL0}=-arctg\left(\frac{{\dot{y}}_{E0}}{{\dot{x}}_{E0}}\right)-{\beta }_0$ (12)

${\alpha }_{rR0}=-arctg\left(\frac{{\dot{y}}_{F0}}{{\dot{x}}_{F0}}\right)-{\beta }_0$ (13)

${\alpha }_{f0}=-arctg\left(\frac{{\dot{y}}_{A0}}{{\dot{x}}_{A0}}\right)-{\beta }_0-{\psi }_{f0}$ (14)

Forces acting on the tractor are depicted in Figure 5b. These forces can be expressed in a form of a vector in a local coordinate system of the tractor:

$F_v=\left(\begin{array}{c}{F}_{\tau }\cos \left({\psi }_{f0}\right)+F_{yf0}\sin \left({\psi }_{f0}\right)+F_{xD0}\\ F_{yrR0}+F_{yrL0}-F_{\tau }\sin \left({\psi }_{f0}\right)+F_{yf0}\cos \left({\psi }_{f0}\right)+F_{yD0}\end{array}\right)$ (15)

where: $F_{\tau }=\tau /r_{0f}$ is a driven force, $F_{xD0}$ and $F_{yD0}$ are constituents of $F_{vD0}$. For the tractor without trailers $F_{vD0}$ is a null vector, $F_{yf0},$ $F_{yrR0},$ $F_{yrL0}$ are lateral friction forces acting at the contact point of the wheel with the ground and can be calculated according to functions

$F_{yf0}=F_y\left(F_{nf0},{\mu }_{f0},{\alpha }_{f0}\right)$ (16)

$F_{yrR0}=F_y\left(F_{nr0},{\mu }_{r0},{\alpha }_{rR0}\right)$ (17)

$F_{yrL0}=F_y\left(F_{nr0},{\mu }_{r0},{\alpha }_{rL0}\right)$ (18)

The normal forces acting on the wheels can be denoted as

$F_{nf0}=g\left(\frac{c_0{m}_b}{h_0}+m_{w0}\right)$ (19)

$F_{nr0}=g\left(\frac{\left(h_0-c_0\right)m_b}{2h_0}+m_{w0}\right)$ (20)

where g is a gravity acceleration. In order to transform these forces to the global coordinate system, $F_{vG}$ can be calculated using rotation matrix $R\left(-{\beta }_0\right)$: $F_{vG}=R\left(-{\beta }_0\right)\cdot F_v$. For the tractor, a vector of generalized force can be denoted as $Q\left(q\right)={\left[\begin{array}{ccc}{F}_x& F_y& F_{\beta 0}\end{array}\right]}^T,$ where $F_x$ is a first and $F_y$ is a second constituents of $F_{vG}$. $F_{\beta 0}$ is a sum of momenta for generalized coordinate ${\beta }_0$:

$F_{\beta 0}=F_{yrR0}{c}_0+F_{yrL0}{c}_0+F_{\tau }\sin \left({\psi }_{f0}\right)\left(h_0-c_0\right)-F_{yf0}\cos \left({\psi }_{f0}\right)\left(h_0-c_0\right)+F_{yD0}\left(d_0+c_0\right)$ (21)

2.4 Dynamic model of conventional and virtual clutch with drawbar system

As kinematic models for conventional and virtual clutch are identical, one dynamic model was used for both of them as well. For this systems, we assume that front wheels have no influence on trajectory. In order to simplify general descriptions, we assume that subscript $i$ represents i-th number of trailer for the whole train. Basic geometrical parameters for analyzed cases are presented in Figure 6. For the system which consists of $n_t$ number of trailer, the configuration can be expressed with generalized coordinate vector as

$q={\left[\begin{array}{ccccccc}x& y& {\beta }_0& {\beta }_1& {\beta }_2& \cdots & {\beta }_{n_t}\end{array}\right]}^T$. The position equation of the most important points of a trailer are shown in Table 2.

Table 2
Position equation of point of the trailer with drawbar system.

Point	Coordinate	Expression
$B_i$	$x_{Bi}$	$x-{\displaystyle \sum }_{j=0}^{i-1}{d}_j\cos \left({\beta }_j\right)-{\displaystyle \sum }_{j=1}^i(h_j+d_{fj})\cos \left({\beta }_j\right)$
	$y_{Bi}$	$y+{\displaystyle \sum }_{j=0}^{i-1}{d}_j\sin \left({\beta }_j\right)+{\displaystyle \sum }_{j=1}^i(h_j+d_{fj})\sin \left({\beta }_j\right)$
$P_i$	$x_{Pi}$	$x_{Bi}+c_i\cos \left({\beta }_i\right)$
	$y_{Pi}$	$y_{Bi}-c_i\sin \left({\beta }_i\right)$
$A_i$	$x_{Ai}$	$x_{Bi}+h_i\cos \left({\beta }_i\right)$
	$y_{Ai}$	$y_{Bi}-h_i\sin \left({\beta }_i\right)$
$D_i$	$x_{Di}$	$x_{Bi}-d_i\cos \left({\beta }_i\right)$
	$y_{Di}$	$y_{Bi}+d_i\sin \left({\beta }_i\right)$
$E_i$	$x_{Ei}$	$x_{Bi}+\frac{1}{2}{v}_i\sin \left({\beta }_i\right)$
	$y_{Ei}$	$x_{Bi}+\frac{1}{2}{v}_i\cos \left({\beta }_i\right)$
$F_i$	$x_{Fi}$	$x_{Bi}-\frac{1}{2}{v}_i\sin \left({\beta }_i\right)$
	$y_{Fi}$	$x_{Bi}-\frac{1}{2}{v}_i\cos \left({\beta }_i\right)$

The kinetic energy of the i-th trailer $T_{ti}$ can be denoted as

$T_{ti}=\frac{1}{2}{m}_{ti}\left({\dot{x}}_{Pi}^2+{\dot{y}}_{Pi}^2\right)+\frac{1}{2}{I}_{ti}{\dot{\beta }}_i^2+m_{wi}\left(\left({\dot{x}}_{Bi}^2+{\dot{y}}_{Bi}^2\right)+\left({\dot{x}}_{Ai}^2+{\dot{y}}_{Ai}^2\right)\right)+2\left(I_{wi}+m_{wi}{\left(\frac{1}{2}{v}_i\right)}^2\right){\dot{\beta }}_i^2$ (22)

Total energy of the system can be denoted as the sum of the energy of individual units: $T=T_b+{\displaystyle \sum }_1^{n_t}{T}_{ti}$.

Figure 6
Drawbar system: a) basic geometrical parameters b) characteristic points.

Side-slip angles and friction forces acting for this system are presented in Figure 7. The side-slip angles for rear wheels of the i-th trailer can be expressed in the arctg function as

${\alpha }_{rLi}=-arctg\left(\frac{{\dot{y}}_{Ei}}{{\dot{x}}_{Ei}}\right)-{\beta }_i$ (23)

${\alpha }_{rRi}=-arctg\left(\frac{{\dot{y}}_{Fi}}{{\dot{x}}_{Fi}}\right)-{\beta }_i$ (24)

Figure 7
Friction forces and slip side angles.

$V{s}_{Fi}$ and $V{s}_{Ei}$ are the directions of velocity of the wheel for i-th trailer. The forces can be expressed as vectors in a local coordinate system of i-th unit: $F_{vi}={\left[\begin{array}{cc}{F}_{xi}& F_{yi}\end{array}\right]}^T$, where $F_{xi}$, $F_{yi}$ is a sum of forces acting longitudinally (subscript x) and laterally (subscript y) relative to the unit. These sums can be noted as $F_{xi}=F_{xDi}$, and $F_{yi}=F_{yDi}+F_{yrRi}+F_{yrLi}$. For a last trailer $F_{xDi}$ and $F_{yDi}$ do not exist, so they are equal to zero. Forces $F_{xi}$, $F_{yi}$ also act on front trailer (i-1), but they must be transformed to the local coordinate system of the antecedent trailer by rotation matrix $R\left({\beta }_{i-1}-{\beta }_i\right)$

$F_{vD\left(i-1\right)}=R\left({\beta }_{i-1}-{\beta }_i\right)\cdot F_{vi}$ (25)

This can be presented in another form

$F_{vD\left(i-1\right)}={\left[\begin{array}{cc}{F}_{xD\left(i-1\right)}& F_{yD\left(i-1\right)}\end{array}\right]}^T$ (26)

The resulting vector involves forces acting on point $D_{i-1}$. The lateral friction force acting on a wheel can be calculated according to the following functions

$F_{yrRi}=F_y\left(F_{nri},{\mu }_{ri},{\alpha }_{rRi}\right)$ (27)

$F_{yrLi}=F_y\left(F_{nri},{\mu }_{ri},{\alpha }_{rLi}\right)$ (28)

The normal force acting on the wheel can be denoted as

$F_{nri}=g\left(\frac{\left(h_i-c_i\right)m_{ti}}{2h_i}+m_{wi}\right)$ (29)

For this system, the generalized force vector can be denoted as $Q\left(q\right)={\left[\begin{array}{ccccccc}{F}_x& F_y& F_{\beta 0}& F_{\beta 1}& F_{\beta 2}& \cdots & F_{\beta n_t}\end{array}\right]}^T$.

$F_{\beta i}$ for trailers can be expressed as a sum of moment for generalized coordinate ${\beta }_i$:

$F_{\beta i}=(F_{yrRi}+F_{yrLi})\left(h_i+d_{fi}\right)+F_{yDi}\left(d_i+h_i+d_{fi}\right)$ (30)

where $F_{yDi}$ does not occur for last trailer.

2.5 Dynamic model of double Ackermann steering system

A higher number of articulations (DOFs) causes a greater number of generalized coordinates for the double Ackermann steering system. Similarly to the previous ones, the position equations of the characteristic points of a trailer are shown in Table 3.

Table 3
Position equation of point of the trailer with double Ackermann steering system.

Point	Coordinate	Position equation
$B_i$	$x_{Bi}$	$x-d_0\cos \left({\beta }_0\right)-{\displaystyle \sum }_{j=1}^{i-1}{d}_j\cos \left({\beta }_{j2}\right)-{\displaystyle \sum }_{j=1}^i{d}_{pj}\cos \left({\beta }_{j1}\right)-{\displaystyle \sum }_{j=1}^i{h}_j\cos \left({\beta }_{j2}\right)$
	$y_{Bi}$	$y+d_0\sin \left({\beta }_0\right)+{\displaystyle \sum }_{j=1}^{i-1}{d}_j\sin \left({\beta }_{j2}\right)+{\displaystyle \sum }_{j=1}^i{d}_{pj}\sin \left({\beta }_{j1}\right)+{\displaystyle \sum }_{j=1}^i{h}_j\sin \left({\beta }_{j2}\right)$
$P_{i1}$	$x_{Pi1}$	$x_{Ai}+c_{i1}\cos \left({\beta }_{i1}\right)$
	$y_{Pi2}$	$y_{Ai}-c_{i1}\sin \left({\beta }_{i1}\right)$
$P_{i2}$	$x_{Pi1}$ $x_{Bi}+c_{i2}\cos \left({\beta }_{i2}\right)$	$y_{Pi2}$
	$y_{Bi}-c_{i2}\sin \left({\beta }_{i2}\right)$	$A_i$
$x_{Ai}$	$x_{Bi}+h_i\cos \left({\beta }_{i2}\right)$	$y_{Ai}$
	$y_{Bi}-h_i\sin \left({\beta }_{i2}\right)$	$D_i$
$x_{Di}$	$x_{Bi}-d_i\cos \left({\beta }_{i2}\right)$	$y_{Di}$
	$y_{Bi}+d_i\sin \left({\beta }_{i2}\right)$	$E_i$
$x_{Ei}$	$x_{Bi}+\frac{1}{2}{v}_i\sin \left({\beta }_{i2}\right)$	$y_{Ei}$
	$y_{Bi}+\frac{1}{2}{v}_i\cos \left({\beta }_{i2}\right)$	$F_i$
$x_{Fi}$	$x_{Bi}-\frac{1}{2}{v}_i\sin \left({\beta }_{i2}\right)$	$y_{Fi}$
	$y_{Bi}-\frac{1}{2}{v}_i\cos \left({\beta }_{i2}\right)$	$G_i$
$x_{Gi}$	$x_{Ai}-\frac{1}{2}{v}_i\sin \left({\beta }_{i2}\right)$	$y_{Gi}$
	$y_{Ai}-\frac{1}{2}{v}_i\cos \left({\beta }_{i2}\right)$	$H_i$
$x_{Hi}$	$x_{Ai}+\frac{1}{2}{v}_i\sin \left({\beta }_{i2}\right)$	$y_{Hi}$
	$y_{Ai}+\frac{1}{2}{v}_i\cos \left({\beta }_{i2}\right)$	${\beta }_{i2}$

Comparing with the previous systems, a vector of generalized coordinates has been extended by angles defining the position of the drawbar ${\beta }_{i2}$. In this system, angles defining the position of the trolley are marked as ${\beta }_{i1}$. Using the above assumptions, a vector of generalized coordinates can be obtained:

$q={\left[\begin{array}{cccccccccc}x& y& {\beta }_0& {\beta }_{11}& {\beta }_{12}& {\beta }_{21}& {\beta }_{22}& \cdots & {\beta }_{n_{t}1}& {\beta }_{n_{t}2}\end{array}\right]}^T$. Due to the presence of four steering wheels, it is necessary to determine the steering angle of wheel, individually. For the double Ackermann steering system this can be denoted as

${\psi }_{fLi}=arctg\left(\frac{h_i}{v_i+h_{i}ctg\left({\beta }_{i1}-{\beta }_{i2}\right)}\right)and{\psi }_{fRi}=arctg\left(\frac{h_i}{-v_i+h_{i}ctg\left({\beta }_{i1}-{\beta }_{i2}\right)}\right)$ (31)

For rear wheels value of the steering angle is the same as in the front, but direction is opposite, so it can be expressed as ${\psi }_{rLi}={\psi }_{fLi}$, ${\psi }_{rRi}={\psi }_{fRi}$.

The kinetic energy of the i-th trailer $T_{ti}$ can be expressed as

$T_{ti}=\frac{1}{2}{m}_{dpi}\left({\dot{x}}_{Pi1}^2+{\dot{y}}_{Pi1}^2\right)+I_{dpi}{\dot{\beta }}_{i1}^2+\frac{1}{2}{m}_{ti}\left({\dot{x}}_{Pi2}^2+{\dot{y}}_{Pi2}^2\right)+\frac{1}{2}{I}_{ti}{\dot{\beta }}_{i2}^2+m_{wi}\left(\left({\dot{x}}_{Bi}^2+{\dot{y}}_{Bi}^2\right)+\left({\dot{x}}_{Ai}^2+{\dot{y}}_{Ai}^2\right)\right)+2\left(I_{wi}+m_{wi}{\left(\frac{1}{2}{v}_i\right)}^2\right){\dot{\beta }}_{i2}^2$ (32)

Total energy of the system can be denoted as the sum of the energies of individual units: $T=T_b+{\displaystyle \sum }_1^{n_t}{T}_{ti}$. Geometrical parameters for these systems are presented in Figure 8.

Figure 8
Double Ackermann steering system a) Basic geometrical parameters b) characteristic points.

Side-slip angles and friction forces acting on the systems are presented in Figure 9. The side slip angles for this system of the i-th trailer can be expressed as

${\alpha }_{fLi}=-arctg\left({\dot{y}}_{Hi}/{\dot{x}}_{Hi}\right)-{\beta }_{i2}-{\psi }_{fLi}$ (33)

${\alpha }_{fRi}=-arctg\left({\dot{y}}_{Gi}/{\dot{x}}_{Gi}\right)-{\beta }_{i2}-{\psi }_{fRi}$ (34)

${\alpha }_{rLi}=-arctg\left({\dot{y}}_{Ei}/{\dot{x}}_{Ei}\right)-{\beta }_{i2}+{\psi }_{rLi}$ (35)

${\alpha }_{rRi}=-arctg\left({\dot{y}}_{Fi}/{\dot{x}}_{Fi}\right)-{\beta }_{i2}+{\psi }_{rRi}$ (36)

Figure 9
Forces and slip angles.

$V{s}_{Fi}$, $V{s}_{Ei}$, $V{s}_{Gi}$, $V{s}_{Hi}$ are the directions of velocity of the wheel for i-th trailer. The forces can be expressed as vector in a local coordinate system of i-th unit: $F_{vi2}={\left[\begin{array}{cc}{F}_{xi2}& F_{yi2}\end{array}\right]}^T$, where $F_{xi2}$, $F_{yi2}$ is a sum of forces acting longitudinally (subscript x) and laterally (subscript y) relative to the unit. These sums can be expressed as

$F_{xi2}=F_{xDi}-F_{yrRi}\sin \left({\psi }_{rRi}\right)-F_{yrLi}\sin \left({\psi }_{rLi}\right)+F_{yfRi}\sin \left({\psi }_{fRi}\right)+F_{yfLi}\sin \left({\psi }_{fLi}\right)$ (37)

$F_{yi2}=F_{yDi}+F_{yrRi}\cos \left({\psi }_{rRi}\right)+F_{yrLi}\cos \left({\psi }_{rLi}\right)+F_{yfRi}\cos \left({\psi }_{fRi}\right)+F_{yfLi}\cos \left({\psi }_{fLi}\right)$ (38)

For a last trailer $F_{xDi}$ and $F_{yDi}$ are not involved, therefore they are equal to zero. Forces $F_{xi2}$, $F_{yi2}$ also are acting on front drawbar of the trailer, but they must be transformed to its local coordinate system by rotation matrix $R\left({\beta }_{i1}-{\beta }_{i2}\right)$:

$F_{vAi}=R\left({\beta }_{i1}-{\beta }_{i2}\right)\cdot F_{vi2}$ (39)

This can be presented in another form as

$F_{vAi}={\left[\begin{array}{cc}{F}_{xAi}& F_{yAi}\end{array}\right]}^T$ (40)

$F_{vi1}=F_{vAi}$, where $F_{vi1}={\left[\begin{array}{cc}{F}_{xi1}& F_{yi1}\end{array}\right]}^T$. Forces $F_{xi1}$, $F_{yi1}$ also act on front trailer (i-1), but they must be transformed to the local coordinate system of antecedent trailer by rotation matrix $R\left({\beta }_{\left(i-1\right)2}-{\beta }_{i1}\right)$:

$F_{vD\left(i-1\right)}=R\left({\beta }_{\left(i-1\right)2}-{\beta }_{i1}\right)\cdot F_{vi1}$ (41)

This can be presented in another form as

$F_{vD\left(i-1\right)}={\left[\begin{array}{cc}{F}_{xD\left(i-1\right)}& F_{yD\left(i-1\right)}\end{array}\right]}^T$ (42)

The resulting vector involves forces acting on point $D_{i-1}$. The lateral friction force acting on wheel can be calculated according to functions

$F_{yfRi}=F_y\left(F_{nfRi},{\mu }_{fi},{\alpha }_{fRi}\right)$ (43)

$F_{yfLi}=F_y\left(F_{nfLi},{\mu }_{fi},{\alpha }_{fLi}\right)$ (44)

$F_{yrRi}=F_y\left(F_{nrRi},{\mu }_{ri},{\alpha }_{rRi}\right)$ (45)

$F_{yrLi}=F_y\left(F_{nrLi},{\mu }_{ri},{\alpha }_{rLi}\right)$ (46)

The normal force acting on the wheel was determined, considering also the mass of the drawbar. We assume that the entire mass of the drawbar is associated with a trailer, for which the drawbar is located at the front. This case is a statically indeterminate system. In order to express the normal force, we consider three static equations and one equation based on the principle of least work (Garcia et al., 2015):

$\{\begin{array}{c}\left(F_{nrLi}+F_{nrRi}-2g\cdot m_{wi}\right)h_i-g(m_{ti}\left(h_i-c_{i2}\right)+m_{dpi}\left(c_{i1}\cos \left({\beta }_{i1}-{\beta }_{i2}\right)\right))=0\\ \left(F_{nfRi}+F_{nrRi}-2g\cdot m_{wi}\right)v_i-g\left(\frac{m_{ti}{v}_i}{2}-m_{dpi}\left(c_{i1}\sin \left({\beta }_{i1}-{\beta }_{i2}\right)+\frac{v_i}{2}\right)\right)=0\\ \left(F_{nfLi}+F_{nrLi}-2g\cdot m_{wi}\right)v_i-g\left(\frac{m_{ti}{v}_i}{2}+m_{dpi}\left(c_{i1}\sin \left({\beta }_{i1}-{\beta }_{i2}\right)-\frac{v_i}{2}\right)\right)=0\\ F_{nrLi}+F_{nfRi}=F_{nrRi}+F_{nfLi}\end{array}$ (47)

The load forces can be obtained by solving the system of equations (47) and expressed as

$F_{nfRi}=\frac{g}{2}\left(\frac{m_{ti}{c}_{i2}}{h_i}+m_{dpi}\left(\frac{c_{i1}\cos \left({\beta }_{i1}-{\beta }_{i2}\right)}{h_i}+\frac{c_{i1}\sin \left({\beta }_{i1}-{\beta }_{i2}\right)}{v_i}+1\right)+2m_{wi}\right)$ (48)

$F_{nfRi}=\frac{g}{2}\left(\frac{m_{ti}{c}_{i2}}{h_i}+m_{dpi}\left(\frac{c_{i1}\cos \left({\beta }_{i1}-{\beta }_{i2}\right)}{h_i}-\frac{c_{i1}\sin \left({\beta }_{i1}-{\beta }_{i2}\right)}{v_i}+1\right)+2m_{wi}\right)$ (49)

$F_{nrRi}=\frac{g}{2}\left(m_{ti}\left(1-\frac{c_{i2}}{h_i}\right)+m_{dpi}\left(-\frac{c_{i1}\cos \left({\beta }_{i1}-{\beta }_{i2}\right)}{h_i}+\frac{c_{i1}\sin \left({\beta }_{i1}-{\beta }_{i2}\right)}{v_i}\right)+2m_{wi}\right)$ (50)

$F_{nrLi}=\frac{g}{2}\left(m_{ti}\left(1-\frac{c_{i2}}{h_i}\right)+m_{dpi}\left(-\frac{c_{i1}\cos \left({\beta }_{i1}-{\beta }_{i2}\right)}{h_i}-\frac{c_{i1}\sin \left({\beta }_{i1}-{\beta }_{i2}\right)}{v_i}\right)+2m_{wi}\right)$ (51)

For this system generalized force vector can be denoted as

$Q\left(q\right)={\left[\begin{array}{cccccccccc}{F}_x& F_y& F_{\beta 0}& F_{\beta 11}& F_{\beta 12}& F_{\beta 21}& F_{\beta 22}& \cdots & F_{\beta n_{t}1}& F_{\beta n_{t}2}\end{array}\right]}^T$. $F_{{\beta }_{i1}}$ for the drawbar can be denoted as:

$F_{\beta i1}=F_{yAi}{d}_{pi}$ (52)

$F_{\beta i2}$ for trailers can be expressed as a sum of moment for generalized coordinate ${\beta }_i$:

$F_{\beta i2}=\left(F_{yrRi}\cos \left({\psi }_{rRi}\right)+F_{yrLi}\cos \left({\psi }_{rLi}\right)\right)h_i+\left(F_{yrRi}\sin \left({\psi }_{rRi}\right)-F_{yrLi}\sin \left({\psi }_{rLi}\right)-F_{yfRi}\sin \left({\psi }_{fRi}\right)+F_{yfLi}\sin \left({\psi }_{fLi}\right)\right)\frac{1}{2}{v}_i+F_{yDi}\left(d_i+h_i\right)$ (53)

where $F_{yDi}$ is not present for the last trailer.

3 SIMULATIONS RESULTS

In this section, the results of numerical simulations are presented. For each system, simulations were conducted for three different tire models as discussed in section 2. For simulation, it was assumed that the logistic train consists of a tractor and four trailers. The geometrical parameters and mass properties of the system are given in Table 4. In order to cause slippage, we assumed that trolleys were lightly loaded. This resulted in the greater impact of the trailer own weight, especially the drawbar, on the system dynamics. This is due to the different location of the center of gravity of the trailer for conventional and virtual clutches. The assumption was that the trolley center of mass, without taking into account the drawbar, was ideally in the center of the trailer. Furthermore, the drawbar mass caused a change in the center of gravity of the location trolley. Therefore, for double Ackermann steering system, the trolley center of gravity is in the middle of the trailer. The drawbar for this system was treated as a separate link. In order to ensure no-slip motion for the tractor, under any circumstances, we assumed a very high value of friction coefficient and simplified tire model for each numerical experiment ($\mu =1000$). The friction coefficient for the trailer wheels for each system was 10; 1; 0.5; 0.1. In addition, one more simulation for each system was performed. It was an ideal ride for the whole train. For this variant, $\mu =1000$ and the simplified tire model were assumed for each wheel. The initial conditions for all scenarios are the same for each simulation and are $q=q_0$ and $\dot{q}={\dot{q}}_0$, where $q_0,{\dot{q}}_0$ are null vectors. The input values are the steering angle of the front wheel and the driving torque applied to the front wheel. These input values are given as a step signal. The steering angle of the front wheel of the tractor ${\psi }_{f0}$ is constant and equal $\pi /18$. The tractor drive torque $\tau$ is a function of time, and its value has been chosen so that the visualization could clearly show the differences between rides and causes slippage. For $t<3$ $\tau =200Nm$ and otherwise $\tau =40Nm$. The logistic train running time in the simulation was 15 seconds. Numerical simulations were performed using Wolfram Mathematica. The results are obtained by solving the system of differential equations using the residual method (Chehab and Laminie, 2005). The length of the train was chosen so that trolley wheel axes were spaced by the same distance from the tractor in the initial position for each system (Figure 10). The visualization did not include the front axle with wheels for conventional and virtual clutch (colored with grey). The results of the experiments are presented in Figures 11, 12 and 13.

Table 4
Geometric parameters and mass properties.

PARAMETERS		Conventional Clutch	Virtual Clutch	Double Ackermann
Mass [kg]	$m_b$	1100	1100	1100
	$m_{w0}$	2	2	2
	$m_{dpi}$	-	-	5
	$m_{ti}$	50	50	45
	$m_{wi}$	2	2	2
Moment of inertia $\left[kg{m}^2\right]$	$I_b$	30	30	30
	$I_{w0},I_{wi}$	0.02	0.02	0.02
	$I_{dpi}$	-	-	0.1
	$I_{ti}$	1	1	1
Dimensions [m]	$h_0$	1	1	1
	$v_0$	0.5	0.5	0.5
	$c_0$	0.3	0.3	0.3
	$c_i,c_{i1}$	0.4	0.3	0.35
	$c_{i2}$	-	-	0.15
	$r_{0f}$	0.2	0.2	0.2
	$r_{0r}$	0.1	0.1	0.1
	$d_0,d_i$	0.2	0.6	0.2
	$h_i$	0.7	0.7	0.7
	$d_{fi},d_{pi}$	0.6	0.2	0.6
	$v_i$	0.3	0.3	0.3
	$r_{ir}$	0.08	0.08	0.08

Figure 10
Visualization of a train in an initial position with: a) conventional clutch, b) virtual clutch, c) double Ackermann steering system.

Figure 11
Train trajectories for virtual clutch steering system, where Red, Blue, Green, Orange, Pink and Black are for different time values in seconds: 0, 3, 6, 9, 12, 15 respectively.

Figure 12
Train trajectories for conventional clutch steering system, where Red, Blue, Green, Orange, Pink and Black are for different time values in seconds: 0, 3, 6, 9, 12, 15 respectively.

Figure 13
Train trajectories for double Ackermann steering system, where Red, Blue, Green, Orange, Pink and Black are for different time values in seconds: 0, 3, 6, 9, 12, 15 respectively.

The obtained trajectories differ significantly in individual simulations. The high friction coefficient (µ=10) ensures no slip. A low coefficient of friction (0.1 and 0.5) causes the slippage of each system with simplified tire model. Although this is least noticeable with the double Ackermann steering system. In contrast, this effect is the most evident in the case of conventional and virtual clutch. Based on the visualization, it can be concluded that a conventional clutch is more susceptible to slipping under these conditions. The slipping in all systems causes the trolleys to move away from the center of the circular trajectory due to centrifugal force. The double Ackermann steering system is less prone to the slip than other systems. The reason for this is the influence of the front wheel and the possibility of turning. In this case, the presence of four active wheels leads to the situation in which the entire weight of the trolley can affect friction. For trailers with double Ackermann steering friction forces act on the trolley both in longitudinal and lateral directions. In the other two systems, only lateral friction forces act on the trailer. In those cases, the castor wheels located at the front cause the situation in which a large part of the trolley weight does not affect lateral friction (only rear wheels are involved). For the smallest considered coefficient of friction, it can be seen that Magic formula tire model is the most slip-resistant. Although in the case of double Ackermann steering system, it does not provide an unambiguous assessment. In that case, the train traveled a greater distance, but in the last stage of motion the trolleys are positioned more transversely to the track when compared with other tire models. For higher coefficient of friction, there is no difference between Magic formula and the sigmoid model. The simplified model seems to be the least stable.

To analyze the influence of different tire models on slip behavior for a higher friction coefficient in a quantitative way, the length of the path without slip was examined. For this purpose, ideal trajectories without slipping were used as references. It was assumed that the slip effect was significant when the positions of the last trolley differed by 200 mm. This value was chosen to allow a clear comparison of simulations. The comparison of the percentage of the path without significant slip impact is presented in Figure 14. This Figure shows for which system the slip begins earlier. For given conditions, a train with conventional clutch slips first, then the virtual clutch and the last one is the double Ackermann steering system. For $\mu =10$, the sigmoid tire model is the most slip-resistant. Magic formula is slightly more prone to slip. For other lower coefficients of friction, the Magic formula is definitely the most slip-resistant, while the simplified model is the least stable. This confirms the dependence of the curves as shown in Figure 4.

The slip start time cannot be the only parameter used to assess the slip. Therefore, the maximum values of the difference in the position of the last trolley were compared with the ideal trajectories (Figure 15).

Figure 14
Percent of path without a significant influence of slip, where V - virtual clutch, C - conventional clutch, A - double Ackermann steering system.

Figure 15
Maximum values of the difference in the position of the last trolley reference to the ideal path, where V - virtual clutch, C - conventional clutch, A - double Ackermann steering system.

For each case, as shown in Figure 15, the trend is similar to the one in Figure 14, but with some exceptions. In the case of the lowest coefficient of friction used, the virtual clutch slipped more heavily for the sigmoid tire model than a conventional clutch. This is indicated by a slightly higher value of position difference. Comparison of tire models for the smallest coefficient of friction is not similar as in the previous analysis. In that case, the conventional clutch slipped less with the sigmoid tire model than with the Magic formula tire model. However, the virtual clutch slipped more with the sigmoid tire model than with the simplified one. The moment at which the position of a trailer deviates mostly from its ideal, non-slip, path is important to determine the intensity of the slip. Therefore, the waveforms of the deviations from the ideal path covered by the last trailer are shown in Figure 16. Figure 16a shows that waveforms for $\mu =0.1$ are fluctuating. It is also noticeable for the simplified model in Figures 16b and 16c. For $\mu =0.5$ and 1, functions are similar for the sigmoidal and MF models. It is not noticeable for $\mu =10$ (Figure 16d), where the waveforms for the MF model are more similar to the simplified model than to the sigmoidal one. These fluctuations arise during a slip. The more susceptible the system is to slippage, the greater there fluctuations.

Figure 16
The waveforms of the difference in the position of the last trolley reference to the ideal (without slip) path, where V - virtual clutch, C - conventional clutch, A - double Ackermann steering system.

Another important parameter for assessing slip is the value of the slip angle. Similar to previous analyses, the maximum slip angle for the last trailer was recorded during the entire ride. Therefore, the absolute values of maximum angle values can be compared (Figure 17). If the value is 90$°$ or close, a full side slip occurs. Figure 17a shows that for $\mu =0.1$ full slip does not occur only in the case of the double Ackermann steering with a sigmoidal tire model. For other variants, with the same values of coefficient of friction, the slip angle is close to 90$°$. In Figures 17b, c and d, trends are similar to previous analyses, but here for a high value of the friction coefficient the differences are more discernible.

The velocity of the train has a significant influence on the intensity of the slip. Due to the occurrence of lateral friction forces that decelerate the train, the maximum speed of the tractor differs for each variant. Figure 18 shows the maximum speed of point $A_0$. In any case, the tractor speed was less than 12 km/h. For the sigmoid and Magic formula tire models, the speed was less than 11 km/h only for $\mu =0.1$. For $\mu =0.5$ and 1, the speed of a tractor with a simplified model was definitely smaller than for other tire models. Assuming that intense slip causes deceleration, the values in Figure 18 for $\mu =0.1$ have the same trend as the values in Figure 15, where the greatest slip occurred for double Ackerman steering system with a simplified tire model. Low tractor speed can cause a significant difference in the position of the last trailer relative to the ideal path shown in Figure 15. This can only be seen for these two analyses. This cannot be seen in any way considering only Figures 14 and 17.

In the study by Paszkowiak and Bartkowiak (2019), where other condition were used, they showed that a conventional clutch is more prone to slipping than a virtual clutch. In that paper, the trolleys were more heavily loaded, the drawbar was longer and center of gravity was the same for both systems. If the trolley is fully loaded, the weight of the drawbar has little effect on the center of gravity. In this case, it can be assumed that the center of gravity is in the geometric center of the trolley. Another set of simulations was performed to check whether different locations of the centers of gravity made the conventional clutch more slippery than the virtual clutch. The results for the new conditions are shown in the Figures 19 and 20. To ensure a similar tractor trajectory with the same input parameters, the train weight remained the same. Only a few parameters were changed: $m_{dpi}=0$, $m_{ti}=50kg$, $c_i$ and $c_{i1}$ is equal to 350 mm. These parameters have been changed for all systems. The comparison of the percentage of the path without significant influence of slip and the maximum values of the difference in the position of the last trolley relative to the ideal path were presented.

Figure 17
Maximum values of the slip angle of the last trolley, where V - virtual clutch, C - conventional clutch, A - double Ackermann steering system, a) $\mu =0.1$, b) $\mu =0.5$, c) $\mu =1$, d) $\mu =10$.

Figure 18
Maximum speed of $A_0$ point of a tractor.

Figure 19
Percent of path without a significant influence of slip neglecting weight of the drawbar, where V - virtual clutch, C - conventional clutch, A - double Ackermann steering system.

Neglecting the weight of the drawbar cause the conventional clutch to be more slip-resistant, except for $\mu =0.1$. The differences between the virtual and conventional clutch are not as significant as in the previous conditions. Therefore, if the front drawbar for conventional and rear drawbar for virtual is longer, the differences are greater. This make the virtual clutch less stable and the conventional clutch more stable.

Figure 20
Maximum values of the difference in the position of the last trolley reference to the ideal path neglecting weight of the drawbar, where V - virtual clutch, C - conventional clutch, A - double Ackermann steering system.

Magic formula takes into account more parameters. This make this model more accurate. Other models used in study are less complex. However, the use of these models may justify application. One of the arguments is the calculation time. The comparison of the average time (with ± 1 × standard deviation) needed to solve the system of differential equations is presented in Figure 21.

Figure 21
Average calculation time - needed to solve system of differential equations a) conventional and virtual clutch, b) double Ackermann steering system.

Due to the same dynamic model, the calculation time for virtual and conventional clutch is presented in the single chart (Figure 21a). For the double Ackermann steering system, the model is more complex than for the drawbar system. In particular the number of differential equations is almost twice as greater. The dynamic model with the Magic formula tire model definitely needs the most time to solve the equations. Other models need less time to solve the equations, and a system with a simplified model needs the least. However, the differences between simplified and sigmoid are not as great as when comparing the Magic formula with other simple models. This is caused by the fact that Pacejka model includes more parameters. This results in a more complicated expression of the lateral friction force.

4 CONCLUSION

In this paper, we showed dynamic models of the three most popular steering systems for logistic trains with three different tire models. The results allow a comparison of their behavior in the presence of a slip effect. The slip effect affects the analyzed systems differently. The conducted study shows that several parameters are needed for the analysis of the slip. The comparison of drawbar systems helped to notice that the length of the longer drawbar and the center of gravity have a significant impact on slip. In the presented conditions, Magic formula tire model proved to be the least prone to the slip except the highest used coefficient of friction. Sigmoid model is the second in this hierarchy. For the smallest coefficient of friction, comparisons do not converge in every case. In such situations, braking is intense and depends on the moment when the slip starts and its speed, what causes discrepancies in results even with minimal changes in initial conditions and input parameters. The comparison showed that the double Ackermann, out of the three compared steering systems, is the least prone to slippage. Only one variant ($\mu =0.1$ with simplified tire model) of this system slipped more than other systems. This might be due to the presence of more DOFs than in drawbar systems. It also shows that comparing the effect on the slip of different systems is not conclusive for low coefficient of friction. The comparison of trajectories where there is no full slip is clearer. Taking into account mass of the drawbar and unloaded trolleys, the virtual clutch is more stable than a conventional clutch. Due to the change in the location of gravity center, the weight of the trailer caused changes in the system dynamics. The greater distance from the center of gravity to the rear axle of the trolley, the smaller the proportion of mass in friction. For fully loaded trolleys, the weight of the drawbar can be neglected, and then the conventional clutch is more stable. The number of active wheels is a factor favoring the double Ackermann steering system. That means that in this system the entire weight of the trolley causes friction forces acting on all four wheels. Due to the front castor wheels in other two systems, only the weight acting on the rear axle of the trailer participated in the friction forces, because the front castor wheels adjust their movement to the actual direction of motion. Another factor favoring double Ackermann may be due to the fact that this system has active swivel wheels, which create longitudinal forces acting on the trolley.

Analysis of calculation time showed that if accuracy was not as important as time, it is worth using a simplified method as a simplified or sigmoid model. The presented dynamic model can be used in the development of models of automatically or autonomously controlled logistic train, if the accuracy related to slip effect is not needed. The computational time of the dynamic model is a crucial parameter, especially in the application for real-time controllers with hardware-in-the-loop solutions.

The dynamic model of the train can be used for other analyses. These include kinematic analysis, which can facilitate to choose the best system for non-slip rides or determining geometrical parameters of the trailers which are optimal for a given trajectory. It is also possible to determine whether there is a potential collision with the surroundings before the application of the model in a real logistic train.

Acknowledgements

We would like to thank Dr. Marcin Pelic for his thoughtful suggestions that helped to improve our research. This research was funded by the Polish Ministry of Science and Higher Education as a part of research subsidy 0614/SBAD/1529.

References

• Adamiec-Wójcik, I., Drąg, Ł., Grzegożek, W. and Wojciech, S. (2019). Calibration of an Articulated Vehicle Model and Analysis of Friction Model in the Connection Between Two Vehicle Units. Journal of Computational and Nonlinear Dynamics, 14(5).
• Ahmad Abu Hatab, R. (2013). Dynamic Modelling of Differential-Drive Mobile Robots using Lagrange and Newton-Euler Methodologies: A Unified Framework. Advances in Robotics & Automation, 02(02).
• Alnahhal, M., Ridwan, A. and Noche, B. (2014). In-plant milk run decision problems. 2014 International Conference on Logistics Operations Management.
• Altche, F., Polack, P. and de La Fortelle, A. (2017). A simple dynamic model for aggressive, near-limits trajectory planning. 2017 IEEE Intelligent Vehicles Symposium (IV).
• Andrzejewski, R. and Awrejcewicz, J., 2005. Nonlinear Dynamics Of A Wheeled Vehicle. Boston, MA: Springer US.
• Baffet, G., Charara, A. and Stephant, J. (2006). Sideslip angle, lateral tire force and road friction estimation in simulations and experiments. 2006 IEEE Conference on Computer Aided Control System Design, 2006 IEEE International Conference on Control Applications, 2006 IEEE International Symposium on Intelligent Control.
• Baranowski, P., Bogusz, P., Gotowicki, P. and Małachowski, J. (2012). ASSESSMENT OF MECHANICAL PROPERTIES OF OFFROAD VEHICLE TIRE: COUPONS TESTING AND FE MODEL DEVELOPMENT. Acta Mechanica et Automatica, 6(2), pp.17-22.
• Bloch, A., Krishnaprasad, P. and Murray, R. (2015). Nonholonomic Mechanics and Control. New York, NY: Springer New York, pp.26-30.
• Brach, R. and Brach, M. (2009). Tire Models for Vehicle Dynamic Simulation and Accident Reconstruction. SAE Technical Paper Series.
• Bravo-Doddoli, A. and García-Naranjo, L. (2015). The dynamics of an articulated n-trailer vehicle. Regular and Chaotic Dynamics, 20(5), pp.497-517.
• Cabrera, J., Castillo, J., Pérez, J., Velasco, J., Guerra, A. and Hernández, P. (2018). A Procedure for Determining Tire-Road Friction Characteristics Using a Modification of the Magic Formula Based on Experimental Results. Sensors, 18(3), p.896.
• Chehab, J. and Laminie, J. (2005). Differential equations and solution of linear systems. Numerical Algorithms, 40(2), pp.103-124.
• Gang Liu, Hongbin Ren, Sizhong Chen and Wenzhu Wang (2013). The 3-DoF bicycle model with the simplified piecewise linear tire model. Proceedings 2013 International Conference on Mechatronic Sciences, Electric Engineering and Computer (MEC).
• Gao, J. and Jia, H. (2017). Control Research for a Small Fixed-Wing UAV During Ground Taxiing. Journal of Harbin Institute of Technology (New Series), 24(2), pp.51-57.
• Garcia, J., Wood, G. and Barrera-Mora, F. (2015). Reactions on Rigid Legs of Rectangular Tables. Applied Mathematics, 06(03), pp.599-616.
• Heisler, H. (2002). Advanced vehicle technology. Burlington: Elsevier, pp.270-310.
• Hulme, K., Kasprzak, E., English, K., Moore-Russo, D. and Lewis, K. (2009). Experiential Learning in Vehicle Dynamics Education via Motion Simulation and Interactive Gaming. International Journal of Computer Games Technology, 2009, pp.1-15.
• Jazar, R., 2009. Vehicle Dynamics. New York, NY: Springer.
• Kayacan, E., Kayacan, E., Ramon, H. and Saeys, W. (2013). Modeling and identification of the yaw dynamics of an autonomous tractor. 2013 9th Asian Control Conference (ASCC).
• Kayacan, E., Kayacan, E., Ramon, H. and Saeys, W. (2014). Nonlinear modeling and identification of an autonomous tractor-trailer system. Computers and Electronics in Agriculture, 106, pp.1-10.
• Keymasi Khalaji, A., Rahimi Bidgoli, M. and Moosavian, S. (2014). Non-model-based control for a wheeled mobile robot towing two trailers. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, 229(1), pp.97-108.
• Kluska, K. and Pawlewski, P. (2018). The use of simulation in the design of Milk-Run intralogistics systems. IFAC-PapersOnLine, 51(11), pp.1428-1433.
• Knez, M. and Gajsek, B. (2015) Implementation of in-plant milkrun system for material supply in lean automotive parts manufacturing. The International Conference on Logistics & Sustainable Transport 2015.
• Kong, J., Pfeiffer, M., Schildbach, G. and Borrelli, F. (2015). Kinematic and dynamic vehicle models for autonomous driving control design. 2015 IEEE Intelligent Vehicles Symposium (IV).
• Li, B., Yang, X. and Yang, J. (2014). Tire Model Application and Parameter Identification-A Literature Review. SAE International Journal of Passenger Cars - Mechanical Systems, 7(1), pp.231-243.
• Liang, W., Medanic, J. and Ruhl, R., 2008. Analytical dynamic tire model. Vehicle System Dynamics, 46(3), pp.197-227.
• Lieb, C., Prinz, T., Gunthner, W. and Fottner, J. (2019). Modeling and Simulation of Generic Handling Operations in In-Plant Milk-Run Systems. 2019 8th International Conference on Industrial Technology and Management (ICITM).
• Lu, Y., Zhang, J., Yang, S. and Li, Z. (2019). Study on improvement of LuGre dynamical model and its application in vehicle handling dynamics. Journal of Mechanical Science and Technology, 33(2), pp.545-558.
• Lucet, E., Grand, C., Salle, D. and Bidaud, P. (2008). Stabilization Algorithm for a High Speed Car-Like Robot Achieving Steering Maneuver. 2008 IEEE International Conference on Robotics and Automation.
• Martinez, J., Avila, J. and Canudas de Wit, C. (2004). A New Bicycle Vehicle Model with Dynamic Contact Friction. IFAC Proceedings Volumes, 37(22), pp.625-630.
• Mondek, M. and Hromcik, M. (2017). Linear analysis of lateral vehicle dynamics. 2017 21st International Conference on Process Control (PC).
• Ortiz, A., Cabrera, J., Guerra, A. and Simon, A. (2006). An easy procedure to determine Magic Formula parameters: a comparative study between the starting value optimization technique and the IMMa optimization algorithm. Vehicle System Dynamics, 44(9), pp.689-718.
• Pacejka, H., Besselink, I. and Pacejka, H. (2012). Tire and vehicle dynamics. Amsterdam: Butterworth-Heinemann.
• Paszkowiak, W., Bartkowiak, T. (2019). Dynamics of logistic train. 2019 15^th International Conference Dynamical Systems - Theory and Applications.
• Paszkowiak, W., Bartkowiak, T. and Pelic, M. (2019). Kinematic model of multiple trailers on a tractor system for production logistics applications. Archives of Mechanical Technology and Materials, 39(1), pp.16-25.
• Piyabongkarn, D., Rajamani, R., Grogg, J. and Lew, J. (2009). Development and Experimental Evaluation of a Slip Angle Estimator for Vehicle Stability Control. IEEE Transactions on Control Systems Technology, 17(1), pp.78-88.
• Svendenius, J. and Wittenmark, B. (2003). Brush tire model with increased flexibility. 2003 European Control Conference (ECC).
• Tzafestas, S., 2014. Introduction To Mobile Robot Control. Amsterdam: Elsevier.
• Yavin, Y. (2007). Modelling the motion of a trolley-like car: The case where the front wheel is constraint to roll on a given plane curve. Mathematical and Computer Modelling, 46(3-4), pp.307-315.

Appendix A: Coefficients and parameters for Pacejka Magic formula

Table A1
Coefficients for Pacejka Magic formula.

Coefficient	Name	Parameters	Formula
$p_{dfn}$	Normalized change in vertical load	-	$\frac{F_n-F_{nN}}{F_{nN}}$
$p_C$	Shape factor	${\rho }_0$	$p_C={\rho }_{CY1}$
$p_D$	Peak factor	${\rho }_1,{\rho }_{15}$	$p_D=\frac{p_{\mu }{F}_n}{1000}$
$p_{\mu }$	Peak friction coefficient		$\left({\rho }_{DY1}+{\rho }_{DY2}{p}_{dfn}\right)\left(1-P_{DY3}{\gamma }^2\right)$
$p_{BCD}$	Stiffness	${\rho }_3,{\rho }_4,{\rho }_5$	$\frac{F_{nN}{\rho }_{KY1}}{1000}\sin \left(2arctg\left(\frac{F_n}{{\rho }_{KY2}{F}_{nN}}\right)\right)\left(1-{\rho }_{KY3}\left|\gamma \frac{180}{\pi }\right|\right)$
$p_B$	Stiffness factor	$p_{BCD},p_C,p_D$	$\frac{p_{BCD}}{p_C{p}_D}$
$p_E$	Curvature factor	${\rho }_6,{\rho }_7,{\rho }_{16},{\rho }_{17}$	$\left({\rho }_{EY1}+{\rho }_{EY2}{p}_{dfn}\right)\left(1-\left({\rho }_{EY3}+{\rho }_{EY4}\gamma \frac{180}{\pi }\right)sign\left(\alpha \frac{180}{\pi }+p_H\right)\right)$
$p_H$	Horizontal shift	${\rho }_8,{\rho }_9,{\rho }_{10}$	${\rho }_{HY1}+{\rho }_{HY2}{p}_{dfn}+{\rho }_{HY3}\gamma \frac{180}{\pi }$
$p_V$	Vertical shift	${\rho }_{11},{\rho }_{12},{\rho }_{13},{\rho }_{14}$	$\frac{F_n}{1000}({\rho }_{VY1}+{\rho }_{VY2}{p}_{dfn})$

Table A2
Parameters for Pacejka Magic formula.

Parameter	Role	Value	Parameter	Role	Value
${\rho }_{CY1}$	Shape factor	1.3	${\rho }_{KV1}$	Change of stiffness with slip	1100
${\rho }_{DY1}$	Lateral friction	$\mu$	${\rho }_{KY2}$	Change of progressivity of stiffness / load	10
${\rho }_{DY2}$	Load influence on lateral friction coefficient	0	${\rho }_{KY3}$	Camber influence on stiffness	0
${\rho }_{DY3}$	Camber influence on lateral friction coefficient	0	${\rho }_{HY1}$	Horizontal shift at load equal to 0 and camber equal to 0	0
${\rho }_{EY1}$	Curvature factor	-2	${\rho }_{HY2}$	Load influence on horizontal shift	0
${\rho }_{EY2}$	Curvature change with load	0	${\rho }_{HY3}$	Camber influence on horizontal shift	0
${\rho }_{EY3}$	Curvature shift	0	${\rho }_{VY1}$	Vertical shift	0
${\rho }_{EY4}$	Curvature change with camber	0	${\rho }_{VY2}$	Vertical shift at load equal to 0	0
$F_{nN}$	Nominal load	4000	$\gamma$	Camber angle	0

Appendix B: Notation

Notation

Parameter	Description
$M\left(q\right)$	System inertia matrix
$C\left(q,\dot{q}\right)$	Centripetal and Coriolis matrix
$Q\left(q\right)$	Generalized forces vector
$\mathcal{L}$	Lagrangian
$q$	Generalized coordinate vector
$T$	Kinetic Energy
$\psi$	Steering angle
$\alpha$	Side-slip angle
$\beta$	Yaw angle
$F$	Force
$m$	Mass
$I$	Moment of inertia around the vertical axis
$\mu$	Lateral friction coefficient
$\tau$	Torque applied to the front wheel of the tractor
$t$	Time
$Subscript$
$i$	Unit number
$f$	Front
$r$	Rear
$R$	Right
$L$	Left
$w$	Wheel
$n$	Normal (force)
$b$	Tractor
$t$	Trailer
$\tau$	Driven

Edited by

Editor:

Jan Awrejcewicz.

Publication Dates

• Publication in this collection
  10 Mar 2021
• Date of issue
  2021

History

• Received
  05 June 2020
• Reviewed
  22 Sept 2020
• Accepted
  23 Nov 2020
• Published
  08 Jan 2021