Modeling and simulation of longitudinal dynamics coupled with clutch engagement dynamics for ground vehicles

During the engagement of the dry clutch in automotive transmissions, clutch judder may occur. Vehicle suspension and engine mounts couple the torsional and longitudinal models, leading to oscillations of the vehicle body that are perceived by the driver as poor driving quality. This paper presents an effective formulation for the modeling and simulation of longitudinal dynamics and powertrain torsional dynamics of the vehicle based on non-smooth dynamics of multibody systems. In doing so friction forces between wheels and the road surface are modeled along with friction torque in the clutch using Coulomb’s friction law. First, bilateral constraint equations of the system are derived in Cartesian coordinates and the dynamical equations of the system are developed using the Lagrange multiplier technique. Complementary formulations are proposed to determine the state transitions from stick to slip between wheels and road surface and from the clutch. An event-driven scheme is used to represent state transition problem, which is solved as a linear complementarity problem (LCP), with Baumgarte’s stabilization method applied to reduce constraint drift. Finally, the numerical results demonstrate that the modeling technique is effective in simulating the vehicle dynamics. Using this method stick-slip transitions between driving wheel and the road surface and from the clutch, as a form of clutch judder, are demonstrated to occur periodically for certain values of the parameters of input torque from engine, and static and dynamic friction characteristics of tire/ground contact patch and clutch discs.


Introduction
Clutch judder between friction pairs of sliding contact during the clutch engagement process is typically defined as a friction-induced vibration and even torsional self-excitation of the transmission system. These oscillations introduce undesired dynamic loads, increase clutch slip and wear of dry clutches and reduce driver comfort [1][2][3][4][5][6][7]. Multibody simulation tools are frequently used to assess the performance of a vehicle to offer the driver better driving comfort characterized by the longitudinal dynamics of vehicle. The acceleration of the vehicle is calculated by integrating the powertrain model with the equations governing the longitudinal behavior of the vehicle [4,6]. However, the simulation of a complete acceleration process during the clutch engagement is not straightforward, as the kinematic constraints imposed on the system change.
Since the 1970s, the clutch judder phenomenon has attracted the attention of many researchers [8]. Some conclusions have been made that clutch-induced drivetrain oscillations could be initiated by different mechanisms of excitation, for example, engine harmonics, difference between static and dynamic coefficient of clutch disc and normal force on clutch face [2,[14][15][16]. Mathematical modeling and numerical studies of powertrain systems considering the contact of clutch engagement dynamics have been presented by some authors. Multi-degree-of-freedom powertrain models are proposed assuming that the effect of the tire slip is neglected, the input torque from the engine is constant and the fluctuating friction torque acts as a source of excitation of the system's torsional vibration during clutch engagement. Numerical results obtained indicate that, for negative gradient of the clutch friction coefficient for a dry friction clutch, the system is subjected to the self-excitation and becomes very unstable for certain values of this gradient [9]. In fact, it is also possible to simplify the model of the powertrain transmission to a four-mass model without losing too much information with respect to clutch vibrations [10][11][12][13][14]. Therefore, to investigate the reduction of friction clutch vibrations, fewer degrees of freedom models are applied to represent the detailed powertrain system [10,11]. The simplified powertrain model enables identical system descriptions during slipping and sticking phase of the dry friction clutch and is used for studying the onset of stick-slip motions [12]. Furthermore, stability analysis of a linearized system during slipping phase has studied how friction characteristics affect system stability and therefore clutch judder or the likelihood of the occurrence of stick-slip. Torque transmitted by the friction clutch is typically numerically calculated using a trial-and-error method, suitable for the multibody systems with several frictionaffected contact points [2,10,[14][15][16]. In the aforementioned literature, most investigations on clutch dynamics only account for torsional degrees of freedom. In Ref. [17], the powertrain model established has been included in a simplified longitudinal vehicle model with only one degree of freedom. Penalty functions are applied to approximate the friction torque from the clutch discs, and they show some limitations by neglecting the stick-slip nature of interfacial friction of the clutch. In Ref. [18], research on the longitudinal dynamics has been undertaken using numerical methods, neglecting the influence of transmission oscillation on the longitudinal acceleration of an automobile. In other studies considering the rotational-translational coupling of tires and vehicle body, self-excited torsional oscillations of the driveline triggered by the friction characteristic of contact patch between tire and ground have been studied assuming that no relative slip between the driving and driven component of the clutch occurs during the process of the torsional oscillation of driveline [19].
In studies of longitudinal vehicle dynamics coupled with powertrain torsional vibration, piece-wise functions of the friction torque acting on the clutch discs, and the friction forces generated in the contact region between tires and road surface are integrated into the dynamic equations of the system. This results in a discontinuously varying structure of the dynamical equations [10,20]. The main difficulty for numerical calculations results from instantaneous changes in the contact forces at transitions from sliding to sticking or reversed sliding. Clocker and Pfeiffer [21][22][23] have presented the theory of non-smooth multibody dynamics with unilateral contacts defined in normal and tangential contact directions for applications similar to that described in this paper. The Linear Complementarity Problem (LCP) formulation embedded in event-driven schemes can be efficiently used to solve the contact-impact problem of multibody systems with frictional unilateral constraints, but the description of the rolling friction law has not been given in detail. In addition, two groups of numerical schemes have been proposed: namely event-driven and time-stepping schemes. Following the non-smooth dynamics approach mentioned above, event-driven schemes in combination with LCP formulations are applied to the study of longitudinal dynamics of the vehicle coupled with the torsional dynamics of the powertrain transmission. For example, the constraint-stabilized method for planar rigid multibody systems with friction-affected translational joints was proposed using non-smooth dynamical approach [24,25].
The Baumgarte method is widely applied to eliminate constraint violations that arise from direct integration of the equations of motion. A quantitative study of the influence of the Baumgarte parameters on the dynamic response of multibody systems is demonstrated in Ref. [32]. As practical examples, the contact forces to drive the motion of a dumbbell and corresponding to the problem of the Painlevé paradox are calculated via the LCP formulation [26][27][28]. Different continuous contact-impact force models for both spherical and cylindrical contact surface in multibody mechanical systems are reviewed in detail [29]. Furthermore, by using a continuous contact force model, the approach to modeling and analysis of flexible spatial multibody systems with clearance of cylindrical joints is presented [30,31]. Following the non-smooth dynamics approach, the event-driven scheme is devised for study of the longitudinal train dynamics, which involves alternative stick-slip phases; however, constraint drift is not considered [33]. LCP formulations of dry frictional contacts with constraint stabilization techniques are also effectively applied to calculate the frictional forces acting on the contact patch between wheels of tracked vehicle and rough road surface [34]. The idea of complementarity is used to model the torque and the relative velocity of one-way clutch but the calculation of the friction torque from the wet clutch, while transitions from stick to slip or slip to stick cannot be achieved, neglecting the reality that during the sticking phase the relative acceleration is not always equal to zero [35].
Based on the above philosophies and with consideration to the presented background, to analyze the dynamic behavior of powertrain transmissions and its impact on the longitudinal acceleration, a comprehensive vehicle model is presented in this paper. Herein, the torsional powertrain and longitudinal vehicle dynamics are combined through the rotationaltranslational coupling of tires and vehicle body. Based on the non-smooth dynamics approach, modeling and numerical algorithms have been proposed in which the friction force in the contact patch between wheels and road surface, as well as frictional torque between two clutch discs, are described as piece-wise functions and defined as generalized friction forces.
The rest of the paper is arranged as follows. In Sect. 2, the physical model of the longitudinal dynamics combined with the powertrain transmission dynamics is developed and described in detail. The governing dynamical equations of the model are derived using Lagrange's method with Cartesian coordinates. In Sect. 3, the corresponding relations are established using complementarities with respect to normal and tangential directions between wheels and road surface, as well as tangential direction between two clutch discs [10]. LCP formulations embedded in the event-driven algorithm are completed as referred to above, with the Baumgarte constraint stabilization technique employed to reduce the numerical drift effects. In Sect. 4, numerical results are obtained that indicate transitions from static to sliding friction and vice versa periodically occur in the contact patch between driving wheel and road surface and also between the clutch discs for some certain parameters of the system.

The model description
The elements of multibody systems for vehicle modeling to study longitudinal dynamics are illustrated in Fig. 1. These include rigid body (the vehicle body), coupling suspension springs and dampers, as well as kinematical connecting in translation to the wheels. The tire deformation in the contact patch is considered [36]. The vehicle body is linked to the front or rear wheel by the spring-damper suspension system with ideal translational joints constraints, therefore the constraints resulting from the connecting elements are holonomic. Generally, the axis of translational joints is not vertical to the vehicle body. The coupling and connection elements, respectively, generate internal forces and torques between the bodies of the system or external forces with respect to the environment. Both are considered as massless elements. Consequently, for the vehicle dynamics herein, it is assumed that the tire shows elastic behavior only in the contact zone, friction force and rolling resistance that act on and the rest of the body are considered to be rigid [37].
It is well established that friction clutch vibrations take place in a low-frequency range, around 5 to 15 Hz. Therefore, modeling can be achieved by considering a limited number of degrees of freedom and only one contact [10]. The powertrain system is therefore modeled as a four-degree-freedom (4-DOF) lumped parameter inertia-damper-spring system as schematically detailed; see Fig. 2 [15]. On the assumption that the sample vehicle is a frontwheel-driven, the driveline is composed of the elastic half-shaft represented by an equivalent stiffness coefficient k 2 and the damping factor c 2 . The parameters k 1 and c 1 are the stiffness coefficient and damping factor of the input shaft, respectively. The engine inertia J e is driven by the engine torque T e . J f is the flywheel inertia. T f is the torque transmitted by the clutch and defined as friction torque exerting on the clutch discs. J c is the equivalent inertia of the clutch disc and the transmission. J c1 is the inertia of the front wheel. By definition M f 1 is the rolling resistance acting on the contact patch between the front wheel and road surface. In addition, a flat plane approximates the road surface near the contact patch [18].

Equations of motion
A set of Cartesian coordinates q = [x 1 y 1 θ 1 x c1 y c1 θ c1 x c2 y c2 θ c2 θ c θ f θ e ] T is used for the mathematical description of longitudinal dynamics of vehicle regarded as the bilater-  (x c1 , y c1 , θ c1 ) and (x c2 , y c2 , θ c2 ) are the coordinates of the mass center of the front and rear wheel. θ c , θ f , θ e are the representations of angular displacements of the clutch discs, flywheel and engine, respectively.
In the following the subscript i ∈ 1, 2 denotes contact point of the constraints from the translational joints between the vehicle body and wheels. The kinematic joints can be described by a set of linear holonomic algebraic equations as The Lagrange multiplier technique is used to derive the equations of motion at the acceleration level in the following form: where M denotes the symmetric and positive definite mass matrix, Q denotes the vector containing all smooth external forces, Φ q is the Jacobi matrix of the constraint of (1). λ = [λ 1 λ 2 ] T is the Lagrange multiplier vector and Q f is the force/torque vector which contains friction force and rolling resistance generated at the ground/tire interface, as well as friction torque from the clutch discs. Wheel friction forces are modeled with the classical Coulomb's friction law [38]. The friction law of Coulomb states that the sliding friction force is proportional to the normal force of a contact. The amount of the static friction force is less than or equal to the maximum static friction force which is also proportional to the normal force F Ni . The sliding friction force has the opposite direction of the relative velocity of the frictional contact. Consequently, it can be of the form where μ 0i is the coefficients of static friction, μ i is the coefficients of sliding friction, v ri = x ci − Rθ ci . The pressure distribution in the contact patch between wheel and road surface implies that the normal force F Ni acts through the center of pressure a distance forward of the wheel center. For equilibrium, a couple exists that must oppose the tire load and its reaction acting down through the wheel center. The couple reacting against the wheel is the rolling resistance represented with a rolling resistance coefficient. The rolling resistance M f i below is defined according to [18]: where δ i is the representation of the rolling resistance coefficient, F Ni is represented as the vertical force acting on the tire. The calculation of the normal force F Ni acting in the tire contact patch has a contribution due to stiffness F Nik and a contribution due to damping F Nic . These forces acting in the normal directions are expressed as [18] F Ni = K y ci δ y ci + C y ciδ y ci δ y ci ≥ 0 0 δ y ci < 0 where m t = mass of tire, K y ci = radial tire stiffness, ζ = radial damping ratio, δ y ci = R − y ci = tire penetration,δ y ci = −ẏ ci = rate of change of tire penetration, and C y ci = 2.0ζ m t K y ci .
Using a Coulomb friction model, the clutch friction torque is [16] where T c is the slipping friction torque, T s is the maximum static friction torque,θ r = θ f −θ c . The generalized friction force vector Q f can be finally expressed in matrix form where As a result, (2) can be rewritten as As a consequence, (8) describes the system under the influence of the tangential contact forces proportional to the multipliers λ T , which has to be determined by contact laws in the following section.

Numerical algorithm
The main difficulty to solve (8) results from instantaneous changes in the contact forces between the wheels and road surface and from the clutch at transitions from sliding to sticking or transitions from sticking to sliding. Moreover, if the contacts are kinematical coupled, the contact forces influence each other. Thus, each new contact generally affects all of the other contact constraints and produces jumps in contact forces. Due to these step changes, induced transitions in the states of contacts may occur, for which either are transitions to sliding or to separation. Formulations for the complementarity conditions of transitions from sticking to sliding in the multibody vehicle system are formulated. This evaluates the transition problem and avoids the combinatorial problem associated with the testing of all possible contact state combinations for the solution without contradiction to the contact laws.

Complementarity conditions of generalized friction forces
Firstly, a complementary formulation for Coulomb's friction law is formulated by introducing the definition of the friction saturation in positive and negative direction [23]: The positive and negative part of the tangential acceleration is expressed aṡ Next, in the same manner, a complementary formulation for rolling resistance friction's law has to be expressed. Similarly, by the definition of the rolling resistance saturation in positive and negative direction, it can be described in the following form: The positive and negative parts of the angular acceleration of the front and rear wheel are defined asθ Equations (9) and (11) can be finally written as (10) and (12) can be finally given in the unified form where:g H 1 =v r1 ,g H 2 =θ c1 ,g H 3 =v r2 ,g H 4 =θ c2 . Therefore, the corresponding relations have to be completed by the complementarities with respect to tangential direction between the wheels and road surface. It is given in the formg Equation (15) is As mentioned above, the torque transmitted through the clutch (both in slipping and stick state) is indicated by T f . The friction torque saturation in positive and negative direction is firstly stated as The positive and negative parts of the relative angular acceleration between two clutch discs are defined asθ Equation (16) also can be expressed as where λ H i = T f . Equation (17) also can be expressed as As a consequence, this transmission behavior is described by the complementary transmission law:g which is restricted toġ T i = 0 (ġ T i =θ r ) (i = 5).

LCP formulations
Combining (3), (4) and (6), the detailed expression is found for the generalized friction forces of (7): where As a result In that way, another form can be given where W H = W T B and B is represented as the vector which is composed of the column vectors of vector S which are not equal to a zero vector. Consequently, the following equation is obtained To keep the constraint violations under control, the Baumgarte stabilization method is used here. This method allows constraints to be slightly violated before corrective actions can take place in order to force the violation to vanish. The goal of Baumgarte method is to replace the differential equation (1) by the following equation: Equation (22) can also be expressed as Substituting the above equation into (23) results in where Therefore, (8) can be written in another form: By usingg T to indicate a transition from sliding (ġ T = 0) to stick or rolling (ġ T = 0), the tangential accelerations in the contact patch between wheels and road surface and also between two clutch discs can be written in the following form: (27) Ifġ T i = 0, (27) can be given byg (14), (19) and (26) into (28) results in (13) and (18), it is found that Therefore, whileg T i = 0, the following method to solve the planar frictional contact problem between the wheels and road surface and also between two clutch discs is proposed by transforming it into a linear complementarity problem. Equations (29) and (30) are integrated in matrix form: Equation (31) is a LCP in standard form and can be solved directly by a pivoting algorithm like Lemke's method or iteratively, with a Block-Gauss-Seidel relaxation scheme [39]. The event-driven schemes detect changes of the constraints (events), such as stick-slip transitions between driving wheel and road surface and that between two clutch discs, and they resolve the exact transition times. Between these events the motion of the system is smooth and can be computed by a standard ODE/DAE integrator with root-finding. If an event occurs, the integration stops and the computation of the contact forces is performed by solving a LCP.

Numerical results
The physical model is illustrated in Figs. 1 and 2. It is assumed that the axis of the translational joints is vertical to the vehicle body shown as Fig. 1. The values of the basic parameters for the investigated vehicle model are accepted as follows [36]: Baumgarte parameters α and β are chosen as positive constants [32]: The Stribeck model is described as the following equation [43]: T S = 100 N·m, T C = 80 N·m,θ s = 10 rad/s For simulation, the engine torque is modeled [44]: T e = T 0 1 + 1.8 · sin(2θ e t) + 0.7 · sin(4θ e t)

Case 1
The following initial conditions are used for the numerical integration: θ e =θ f =θ c =θ c1 = 20 rad/s,ẋ 1 =ẋ c1 =ẋ c2 = Rθ c1 = Rθ c2 = 5 m/s The vehicle is assumed to move over a flat level road surface with dry asphalt and concrete along a straight line. Simulations are performed using the friction coefficient determined via (3) with a μ 0 value of 0.8 and a μ value of 0.75 [37]. Typical value of rolling resistance coefficient δ is set to 0.005 m. The torque from engine T e is assumed as a constant value of 200 N·m. The maximum static friction torque T s of the clutch is large enough to ensure that the slipping motion between two discs will never occur once the clutch is locked. Figures 3, 4 and 5, respectively, show the pitch angular rotation of the vehicle body θ 1 , the relative velocity V r1 and frictional force F f 1 of contact patch between the driven wheel and road surface, and the relative velocity V r2 and frictional force F f 2 of contact patch between the rear wheel and road surface. Based on the numerical analysis, it has been found that the wheels are purely rolling and the torsional vibration from the powertrain transmission system leads to the fluctuation of frictional force between wheels and road surface during the launch phase (see Figs. 4 and 5). Figure 3 illustrates that the effect of vibration absorber results in the decreased amplitudes of pitch angular θ 1 which   Fig. 6, the normal force exerted on the front wheel gradually decreases, and, on the contrary, the force acted on the rear wheel increases. In the final phase, the observation F N1 < F N2 is qualitatively identical to the simulation results. It is evident from Fig. 7 that the clutch plates with zero relative velocity stick with each other, maintained by static friction torque.

Case 2
The vehicle is still assumed to move over a horizontal road surface with dry asphalt and concrete. The numerical solution is achieved using the friction coefficient determined via (3) with a μ 0 value of 0.8 and a μ value of 0.75.
Typical value of rolling resistance coefficient δ is set to 0.005 m. The output torque from engine T e is set to constant value 323.47 N·m, so that the impact of engine forced vibration on clutch dynamics is excluded [40]. This value for the engine torque is identified through a number of iterations as the point at which the output torque is high enough to initiate stick-slip in the clutch.
For the assessment of driving comfort influenced by the variation of the friction coefficient against the relative angular velocity in the frictional clutch contact, the parameters T S and T C determined via (5) are with a value of 400 N·m and 320 N·m, respectively. As can be seen in Fig. 8, during steady state phase, friction-induced stick-slip vibrations (clutch judder) occur periodically in the powertrain transmission system, which yield the periodic oscillations of pitch angular acceleration that can characterize driving comfort. This behavior results from the transient response of the multiple flexible-body system interacting with the sticking and slipping friction torques, resulting in multiple transitions into and out of sticking in the clutch. Furthermore, the torsional vibration behavior of the clutch discs is interpreted as a measure of the clutch judder. The torque variations in the driveline will vary the drive forces at the tire-ground contact and thus may act directly to generate longitudinal acceleration oscillations in the vehicle body (see Figs. 9 and 10). Figures 11 and 12 indicate that the front and rear wheels are purely rolling on the road surface (ġ T 1 = 0 andġ T 3 = 0).
The output torque from engine T e is set to constant value 100 N·m. The Coulomb friction model and the Stribeck model are used to describe the friction torque acting on the clutch Fig. 11 Steady state response of V r1 and F f 1 : zoom on V r1 and F f 1 Fig. 12 Steady state response of V r2 and F f 2 : zoom on V r2 and F f 2 discs, respectively. As can be seen in Fig. 13, the Stribeck effect increases the amplitude of relative angular velocity of clutch discs. Compared with the Coulomb friction model, during the slipping phase of clutch, the negative gradient of dynamic friction torque with relative angular velocity will worsen the clutch oscillation.

Case 3
The vehicle is assumed to move over a horizontal road surface with ice along a straight line, as such the friction coefficient for the road will be reduced. Simulations are carried out using the friction coefficient determined via (3) with a μ 0 value of 0.1 and a μ value of 0.07. Typical value of rolling resistance coefficient δ is still set to 0.005 m. The output torque from engine T e is set to −168.48 N·m. This value is chosen as it represents the smallest value at which periodic stick-slip can be identified in the simulation results. The negative torque implies that the vehicle is a rear-drive-wheel vehicle moving in the opposite direction to cases 1 and 2. The parameters determined via (5) with a T S value of 400 N·m and a T C value of 320 N·m. In Fig. 14, on steady state conditions, the simulation results for the relative angular velocity and clutch torque are shown, demonstrating that the clutch discs remain stick with each other all the time. As depicted in Figs. 1 and 2, during driving the transmission   front wheel and road surface is influenced by the difference of static and dynamic friction coefficients of the road surface. This interacts with torsional vibrations of the powertrain and the variation of pitch angular acceleration to produce a time varying wheel force, resulting in stick-slip. Thus, the above numerical results obtained can be used for interpretation of power hop phenomena with pitch oscillation superimposed on the forward motion of vehicles [41].
As the effect of engine harmonics on oscillations of driveline is taken into consideration, for simulations the parameter T 0 is set to 60 N·m. As shown in Fig. 18, the fluctuating engine torque acts as a source of the excitation of torsional vibration in the driveline system. Figure 19 shows the drift of constraint equations achieved by the Baumgarte stabilization technique can be limited to a small scale ( Φ < 10 −4 m).

Conclusions
This paper presents an effective model and formulation of the longitudinal dynamics of the vehicle affected by the torsional dynamics of powertrain transmission. The system is modeled as a multi-rigid-flexible-body system with bilateral constraints. The proposed model is composed of vehicle body, wheels and suspension system as well as a simplified powertrain system with clutch. Based on the non-smooth dynamics approach, the dynamical equations of the system with normal and tangential contact forces are derived at the acceleration-force level using the Lagrange multiplier technique, and geometric constraints of the translational joints are considered as bilateral constraints. Friction forces between wheels and road surface, as well as friction torque from the clutch, are included as the generalized frictional forces, which lead to the discontinuity of the dynamical equations.
Firstly, the Coulomb friction law is applied to the frictional contacts between wheels and road surface. Coulomb's law and Stribeck's friction model are used to describe the friction torque acting on the clutch, respectively. Complementary relations between friction saturation and tangential accelerations are formulated. Additionally, the complementary formulations of friction between wheels and road surface are established to determine state transitions from pure rolling to slipping. Next, by introducing the Baumgarte stabilization method into the non-smooth equations of motion, the problem of the stick-slip transitions between wheels and road surface and in the clutch is solved as a LCP. Based on eventdriven schemes, the numerical algorithm efficiently detects the stick-slip transitions, which periodically occur between wheels and road surface and in the clutch for specified input torque and friction coefficients. By analyzing the results numerically obtained under certain conditions, it demonstrates that the stick-slip motions (clutch judder) between two clutch discs are triggered and significantly influenced by the clutch friction torque. Compared with the Coulomb friction model, the Stribeck effect increases the amplitude of relative angular velocity of clutch discs. Throughout the three cases studied, it is demonstrated that the occurrence of periodical stick-slip motions between wheels and road surface is closely dependent on torsional dynamics of powertrain transmission and the friction characteristics of road surface. The influence of engine harmonic on the powertrain dynamics has been numerically investigated.
The approach proposed in this paper is also suitable for the study of the longitudinal dynamics of the rear-drive-wheel or all-drive-wheel vehicles. It is used to investigate the friction-induced stick-slip phenomena in the contact patch between wheels and road surface, and also on the contact interface between two clutch discs. Furthermore, regarding the friction models mentioned above, Coulomb model shows acceptable results, but the friction behavior in the facing materials is highly dependent on relative velocity and temperature of facing materials. Therefore, one alternative for future studies is the use of different friction models in the dry clutch and tire modeling [42,43].