Automotive dry clutch fully coupled transient tribodynamics

Determining the root causes of Noise, Vibration and Harshness (NVH) phenomena in modern automotive drivetrains is a task of critical importance. This research investigates the stability of dry clutch systems vibrational behaviour during engagement. A fully coupled dry clutch numerical model including the influence of friction is presented and validated using vehicle measurements. The clutch component frictional properties are measured using parts that exhibit aggressive NVH behaviour using representative tribometric experiments. The validated numerical tool highlights the occurrence of instabilities which are caused by modal couplings, particularly between the input shaft bending and clutch disc radial motions. Such a validated transient dynamics model of a dry clutch system has not hitherto been presented in the open literature.

Abstract Determining the root causes of Noise, Vibration and Harshness (NVH) phenomena in modern automotive drivetrains is a task of critical importance. This research investigates the stability of dry clutch systems vibrational behaviour during engagement. A fully coupled dry clutch numerical model including the influence of friction is presented and validated using vehicle measurements. The clutch component frictional properties are measured using parts that exhibit aggressive NVH behaviour using representative tribometric experiments. The validated numerical tool highlights the occurrence of instabilities which are caused by modal couplings, particularly between the input shaft bending and clutch disc radial motions. Such a validated transient dynamics model of a dry clutch system has not hitherto been presented in the open literature. Greek symbols a Mass proportional coefficient b Stiffness proportional coefficient f n Damping ratio € h Angular acceleration _ h Angular velocity h Angular displacement l k Kinetic coefficient of friction l st Static coefficient of friction x n Natural frequency
Many systems with frictional interfaces exhibit stick-slip behaviour, which often leads to self-excited vibrations [15]. The dynamics and oscillatory behaviour friction disc systems, such as brakes and clutches, can be examined by analysing their stability and investigating the oscillation amplitudes in the unstable regimes. Bostwick et al. [16] developed a torsional, two degree-of-freedom clutch model, to simulate the engagement process and investigate the self-excited oscillations generated as a result of friction coefficient variation. It is found that the amplitude of motion depends on the system parameters and not on the initial conditions or initial system energy. Moreover, self-excitation is shown to occur when the friction coefficient decreases with sliding speed (negative slope), something that Centea et al. [15] observed in their study also. Self-excited oscillations in brake systems have been analysed by Soobbarayen et al. [17] to characterize brake squeal. The results showed that squeal noise is generated due to the contact of the disc and the pad and the occurrence of self-excited vibrations. Nonlinearities involving friction can generate coalescence of modes that can lead to instability.
Systems exhibiting friction nonlinearities may lose their stability, accompanied by alterations in their qualitative structure at equilibrium positions (bifurcations) by variation of their nonlinear parameter values [18]. Kuznetsov et al. [19] showed that stability loss can occur under the two following mechanisms: (i) the fold bifurcation, when the real part of one eigenvalue becomes zero or positive and (ii) the Hopf bifurcation, when a complex eigenvalue pair becomes positive or zero. In a supercritical Hopf bifurcation a stable equilibrium switches to unstable and a stable limit cycle of finite amplitude appears. Elmer et al. [20] and Lim et al. [21] after investigating a sliding mass discussed the occurring bifurcations due to friction variation, which are linked with the change in the stability of the equilibrium positions and are meant to cause chaotic motion, stick-slip behaviour or oscillatory sliding of the mass. Finally, Putra and Nijmaijer [22] observed the limit cycling in a nonlinear system using numerical methods.
Mode coupling has been shown to be a mechanism leading to aggressive NVH during clutch actuation, particularly when there are high relative velocities between the mating surfaces of the annular discs [23,24]. Trinh et al. [25] and Wickramarachi et al. [26] developed dry clutch models to simulate the wobbling and axial motions of the pressure plate. A key assumption is that the pressure plate is set into wobbling motion due to initial conditions imposed. It was observed that this motion can lead to aggressive NVH behaviour during the engagement process. A six degree of freedom (dof) model considering the pressure plate motions was used to conduct a parametric study of the system stability. In order to simulate the pressure plate flexibility, the plate was divided into four segments connected by torsional springs. It was shown that mode coupling of the pressure plate tilting modes was responsible for the unstable behaviour. The behaviour could be moderated by modifying the pressure plate thickness, the coefficient of friction and the structural stiffness of the pressure plate. Moreover, Trinh et al. [27] studied the clutch system oscillatory motion using the limit cycles of the friction disc. The results showed that the selfexcited system under examination, tends to the limit cycles when static equilibrium is unstable.
Senatore et al. [28] created an eight dof analytical model to examine instabilities during the slipping phase of the clutch system. The model includes the key components of a clutch system (flywheel, clutch disc, pressure plate). A nonlinear profile of the cushion spring stiffness was introduced and a constant value for the coefficient of friction was assumed. The degrees of freedom introduced include the pressure plate tilting as well as the axial motion of the pressure plate, however the tilting degrees of freedom were uncoupled with the rest of the system. Initial conditions applied to impose a tilt on the pressure plate were introduced. Eigenvalue analysis was used to identify unstable modes (291.8-577.8 Hz between 0.3-0.6 mm cushion deflection). The associated mode responsible for this instability was identified as the pressure plate tilting motion.
Minas et al. [29] drafted a two dimensional, six dof analytical model to examine the oscillatory behaviour of the various clutch components during various pull away manoeuvres. The model used engagement data (friction torque, coefficient of friction) taken from a three-dof torsional model, simulating the engagement process. Different clutch pedal actuation times were selected for each engagement manoeuvre and the vibrations of the clutch subcomponents were compared. Coefficient of friction measurements from clutch systems that exhibit high frequency NVH issues behaviour were used as an input to the mathematical model. It was found that the radial motion of the clutch disc is responsible for instabilities produced with a natural frequency at 370 Hz. Moreover, during more aggressive pull away manoeuvres (reduced clutch pedal actuation time), greater clutch pedal travel is required to synchronize the angular velocity of the flywheel and the clutch disc. This leads to higher clamp load and higher friction forces. As a result, a higher level of instabilities appear which are associated with the clutch disc radial mode.
Aktir et al. [30] introduced a three-dimensional finite element model to predict high frequency NVH behaviour in dry clutches. The natural frequencies of the pressure plate and clutch disc were identified using impact hammer testing. The dynamic response of the model was validated against measurements of the whole clutch assembly excited by an electromagnetic shaker attached to the flywheel. After conducting stability analysis, it was found that high frequency NVH behaviour (up to 700 Hz) appears due to mode coupling of the input shaft bending motion and the clutch disc radial motion. A constant friction coefficient value was used during engagement, however sensitivity analysis showed that an increased coefficient of friction led to amplification of the clutch instabilities.
Fidlin et al. [31] presented a numerical study of friction disc wobbling dynamics, simulating the clutch disc and pressure plate assembly. A two-dof model was developed including the input shaft. It was found that the relationship between the friction circulatory terms (energy input due to the sliding contacts and damped energy) can determine the threshold of instability. Moreover, if there is no structural damping under consideration the system can be unstable if the slip between the two discs is small and the kinetic friction coefficient is large ([ 0.3). Limit cycles were presented to describe instabilities for different slip speeds. For higher slip speeds the limit cycle can present chaotic behaviour with increased vibration amplitudes compared to lower slip speeds. The length of the input shaft is also shown to affect the limit cycle (the size of the limit cycle increasing with increased length). The model proposed is shown to effectively predict NVH behaviour in clutches with natural frequencies between 250-400 Hz.
Herve et al. [32] developed a two-dof model, representing the clutch disc and the flywheel (including gyroscopic terms and friction), in order to explore mode coupling instabilities. The study included eigenvalue analysis to determine the system stability after considering gyroscopic and circulatory terms. The role of structural damping was also explored due to its influence on mode coupling instabilities [33]. It was found that altering the inertia and stiffness of the components that create their modes associated to themed coupling instabilities, affects the coalescence point of mode coupling instability, which can then be shifted to stabilize the equilibrium. Moreover, when proportional damping was included the gyroscopic action was shown to be negligible.
The published literature has shown the importance of the frictional interface on clutch stability, however commonly friction has been included as a constant coefficient. In this study a transient dynamics model considering the independent motions of the dry clutch sub-components has been developed and simulated during a vehicle pull away manoeuvre. The transient behaviour of the actual physical system is associated with the nature of the cushion spring force, diaphragm spring force and coefficient of friction. The former spring forces are nonlinear functions of the spring deflections, which change as the clutch engagement manoeuvre proceeds. The friction coefficient varies with the normal load and component relative speed, which again changes during the manoeuvre. Experimentally measured coefficient of friction is incorporated as a function of the slip speed and contact pressure. The validity of the clutch system model has been verified against measurements taken in-situ in a vehicle. Such a detailed analytical clutch model and validated complete simulation results have not been previously presented in the open literature.
This paper presents a validated dynamic dry clutch model for the investigation of unwanted NVH behaviours. Initially, the mathematical formulation of the transient dynamic behaviour of a dry clutch system is described in the Methodology section. In addition, the Methodology section describes the experimental approach used to determine key frictional characteristics required to model the clutch using components that previously presented aggressive NVH behaviour. Next the dynamic model is validated using on-vehicle experimental measurements. The Results section presents the dynamic instabilities that can be related to potential aggressive NVH behaviour. An examination of the transient behaviour of the clutch system during a typical engagement manoeuvre shows the instabilities present and mode coupling instabilities are examined. Finally, the key Conclusions and suggestions for future work are discussed.

Methodology
A typical dry clutch used in B-segment vehicles is studied in this work. A short description of the key components the system comprises, and their functions is provided herewith. The primary function of the clutch is the regulation of torque transmission from the engine to the drivetrain. The mechanism of torque transfer is through the frictional forces created at the pressure plate-clutch disc and flywheel-clutch disc interfaces, controlled through actuation of the clamp load. The flywheel is connected to the engine crankshaft and the clutch pressure plate, all of which rotate in unison, whereas the clutch disc is connected to the transmission input shaft, as shown in Fig. 1.
The methodology expounded in this research employs a transient dynamics model using lumped parameters to predict the oscillatory behaviour of the clutch components. The primary rotational and translational motions in the model are fully coupled through the frictional forces generated at the clutch disc interfaces. Experimental measurements were conducted to determine the frictional behaviour of clutch components that exhibit aggressive NVH behaviour, using a pin on disc test rig.

Transient dry clutch dynamics model
A thirteen dof model comprising torsional, axial, tilting/wobbling, radial and bending motions of the clutch components has been developed. A schematic representation of the model is shown in Fig. 1. The dofs of primary importance are the torsional dof of the clutch disc (h cd ) and flywheel ðh fw Þ, each of which has an associated effective mass moment of inertia I cd , and combined flywheel I fw À Á and engine I e ð Þ, respectively. The clutch disc has all six dofs of a rigid body to fully describe its behaviour. The dofs associated to the pressure plate are its tilting motion (rotation) about the X-and Y-axis and its translation in the Z-axis. The pressure plate inertia ðI pp Þ comprises the pressure plate cover (bolted on the flywheel and rotating with the same angular velocity, _ h fw ) and the diaphragm spring. Finally, the bending motion of the transmission input shaft, around X-and Y-axis, have been considered along with its torsional degree of freedom, h in .
The dofs utilised in the clutch model are depicted in the following vector: h pp;y : Tilting motion of the pressure plate around Y-axis. y cd;y : Radial motion of the clutch disc at Y axis. y cd;x : Radial motion of the clutch disc at X-axis. z cd : Axial motion of the clutch disc. h cd;x : Tilting motion of the clutch disc around X-axis. h cd;y : Tilting motion of the clutch disc around Y-axis. h insh;x : Bending motion of the input shaft around X-axis. h insh;y : Bending motion of the input shaft around Y-axis. Figure 2 depicts the employed degrees of freedom for each component: The dry clutch model includes nonlinear representations of the cushion spring (k cs Þ and diaphragm spring (k ds Þ stiffnesses [34]. The former corresponds to the clutch disc normal compression and reaction between the pressure plate and flywheel, and the latter describes the connection of the pressure plate with the pressure plate cover. Constant (linear) stiffness terms are considered for: pressure plate tilting (k tiltpp Þ, clutch disc tilting ðk tiltcd Þ, input shaft bending (k bend Þ and the clutch disc radial motion relatively to the input shaft ðk radial Þ. When the clutch disc is compressed in between the flywheel and pressure plate, friction is generated, as shown in Fig. 3. The friction force is a function of the kinetic friction coefficient calculated using the relative velocity of the contiguous surfaces and the normal force. Due to the clutch disc radial motion, eccentricity is generated between the clutch components (clutch disc and flywheel/pressure plate). This eccentricity and the relative radial velocities of the discs change the direction of the resultant friction force vector. However, the friction in radial direction is negligible compared to the friction generated due to the relative rotation of the components and consequently in this study only the latter is considered [29,35].
In order to simulate the tilting motion of the pressure plate, the radial motion of the clutch disc and the tilting behaviour of the clutch disc require an initial condition (minor variation in the application of the circumferential stiffness). To realise this initial condition a small deflection is applied on the input shaft and in addition the nonlinear stiffness profile of the lower and left hand stiffness is multiplied by a factor k = 0.9, whereas the upper and right domains are divided with the same factor [25][26][27].
The friction forces generated in each segment of the interfacial region can be calculated as: where l kin is the kinetic friction coefficient and the effective radius ðR ef Þ is calculated as [36,37]: The effective radius is a function of the inner and outer radii of the clutch disc (r in and r o , respectively). The equations of motion representing the clutch model can be written in matrix form: is the damping matrix and x f g is the vector described by Eq. (1). The inertial matrix is given below: The stiffness matrix is derived as follows: where K s is the symmetric structural stiffness matrix and K f is the asymmetric matrix, which originates from friction. Any gyroscopic terms are not considered in this study, as Senatore et al. [28] showed that the gyroscopic terms for the tilting modes of the pressure plate are uncoupled with the other dofs. The complete matrix [K] can be found in the ''Appendix''. The asymmetry in this matrix can lead to system instability for certain combinations of stiffnesses and friction coefficients [32,[38][39][40].
The external excitation vector is given by: where T in is the engine torque, T f is the friction torque generated due to the contact of the clutch components, T res is the resistive torque of the drivetrain and F ds is the preload from the diaphragm spring before the manoeuvre starts. The equations of motion were solved using a Runge-Kutta iterative method (ode45 function available in commercial software MATLAB). It is assumed that internal damping is low, thus proportional damping can be used as follows [7]: where [K] is the stiffness matrix and the coefficient b is calculated as follows: where x n are the natural frequencies calculated after solving the eigenproblem for the undamped system, using the (linearized) symmetric stiffness matrix [39] The term f n corresponds to the damping ratio for each particular mode shape of the system. In this study a damping ratio of 0.05% was selected at a natural frequency of 380 Hz [41]. The eigenproblem is solved as follows: where k ¼ x 2 n . The eigenvectors are given by: In order to conduct stability analysis, a vector solution of the following form is used: Substitution of the above function into the equations of motion give: The characteristic polynomial is then extracted as: The stability of the system can be examined after calculating the real parts of the complex eigenvalues. When the real part is positive, the system is unstable at this natural frequency. Conversely, a negative real part corresponds to a stable system mode [30,42]

Clamp Load and friction torque
When the surfaces of the clutch disc are compressed between the flywheel and the pressure plate the generated friction enables the transmission of the engine torque to the rest of the vehicle drivetrain. Assuming a uniform pressure distribution across the contacting surfaces, the friction torque T f À Á is given by [36,43]: where F cushion is the cushion spring reaction force (described further below). The dynamic friction coefficient l kin is determined from experimental measurements described in the next section. The dynamic coefficient of friction ðl kin Þ which couples the tilting, radial, bending and axial motions of the clutch components with the torsional degrees of freedom can be calculated as [9,10,29,44].
where m s is the gradient m s ¼ dl k d _ h [10] and l st is the static coefficient of friction, which can be predicted using the intercept of the linear fit of the kinetic During clutch engagement and disengagement, the total applied load on the clutch disc can be referred to as clamp load. This is the resultant of all forces acting simultaneously on the clutch disc surface. From the schematic representation of a typical dry clutch system shown in Fig. 4 it can be seen that the clamp load depends on the force exerted by the diaphragm spring, ðF diaphragm Þ, cushion spring, ðF cushion ) and straps (F strap Þ during the travel of the pressure plate.
The resultant applied force, F clamp , at the clutch disc surface is calculated as: The force from the straps is negligible compared to the force applied from the diaphragm and cushion springs and is therefore neglected in the current analysis [29]. The simulation of the clutch engagement manoeuvre uses an exponential profile for the cushion spring force, as well as for the diaphragm spring force. The force-deflection profiles (which are based on information received from industry) are shown in  [26]. The diaphragm spring force can be described as: where the coefficients describing the spring force are a ds =2905:56656, b ds ¼ 499:77521,c ds =3:04294 An example showing the nonlinearity of the diaphragm spring can be also found in [34]. The maximum compression of the cushion spring is considered to be at 0.6 mm, and at this point the clutch is described as closed.
The engine torque profile (T in Þ and the friction torque capacity of the clutch are shown in Fig. 6.
Finally, the resistive torque from the rest of the drivetrain is considered, assuming a flat road surface and negligible aerodynamic drag during vehicle pull away. The resistive torque T out ð Þcan be approximated by considering the rolling resistance from the vehicle tires F rol ð Þ: where i is the gear ratio (idle to 1st gear in this study), f is the coefficient of the rolling resistance, R wheel is the radius of the wheel and W is the vehicle weight [45].

Friction coefficient experimental measurements
A pin-on-disc tribometer was used to measure the coefficient of friction corresponding to various clutch engagement conditions. Representative contact pressures and slip velocities between the flywheel/pressure plate and clutch disc samples were used. Samples extracted directly from a clutch disc with dimensions 15 9 15 mm (height and width, respectively) were used for the experiments. The disc used for the current study is 100 mm in diameter and has a surface finish created to replicate that of a clutch pressure plate [41]. The dynamic friction coefficient was determined for the same sliding distance (150 m) for each sample (according to ASTM G115-10 standard). This distance was selected to limit any bulk temperature rise as a result of friction energy dissipation. Thus, the friction coefficient is a function of the component relative speed and normal load only. The results of the experimental measurements of the dynamic friction coefficient under representative contact pressures and relative velocity for the engagement manoeuvre considered are shown in Fig. 7. In order to replicate the kinematic behaviour at the interface the sliding velocity in the pin on disc rig is matched to the slip velocities that occur at the effective radius of the clutch disc during synchronisation.
The lowest friction coefficient values occur at high sliding velocities and lower clamp load. Conversely, the highest friction coefficient values are found at lower sliding velocities and higher clamp loads [46]. The slopes of the linear fits shown in Fig. 7 provide values for the friction coefficient gradient m s , which is used as an input to the simulations. The gradient is a function of the clamp load and thus, it changes accordingly during the engagement manoeuvre. The static coefficient of friction, l st , for the system under investigation is taken as 0.55, corresponding to linear projection at 0 rpm relative speed under 1200 N. The following Table 1 contains ranges of the input parameters used in the clutch mathematical model (exact values are not provided due to commercial sensitivity).
The validation of the dry clutch model was accomplished by comparing the numerical simulation results against experimental measurements obtained from a vehicle. The vehicle powertrain layout is shown in Fig. 8. The experiment starts with the vehicle traveling with an initial speed of 3.91 km/h and the clutch starts from a disengaged position with clutch pedal position at 60% compression. At time zero the pressure plate first touches the clutch disc and the engagement process starts. Gradually, from 0-4 s the driver releases the clutch pedal from 60% clutch pedal position to 51% and from 4-7.2 s the clutch pedal is held at this position. Figure 9 shows the clutch pedal position along with the corresponding cushion deflection. The torque profile during this time is shown in Fig. 10. The acceleration of the input shaft bending motion was measured using a triaxial piezoelectric accelerometer mounted in the vicinity of the input shaft bearing with a sampling frequency of 3 kHz. The distance of the bearing to the tip of the input shaft is L = 0.109 m (as shown in Fig. 8).
The engine torque profile shown in Fig. 10, is estimated using measured air mass, injected fuel mass data and exclusion of all auxiliary loads (alternator, A/C). The clutch pedal position was measured using an array of hall sensors and a magnetic target on the hydraulic master cylinder piston (which is proportional to the pedal pad travel). The sensor signal is processed by pulse code modulation method (PCM) and translated over the Controller Area Network (CAN).
In order to validate the dry clutch model, the input shaft bending motions ðh insh;x and h insh;y Þ were simulated using the engine torque (Fig. 10) and clutch pedal position as input data. The acceleration was calculated at the location of the input shaft bearing and the magnitude and frequency were compared with the experimental measurement. Due to the axisymmetric oscillatory behaviour only the bending motion in Y-axis will be presented below for the sake of brevity.
The magnitude of the input shaft bearing acceleration in the Y-axis is shown in Fig. 11a. During the   Pressure plate moment of inertia around X-, Y-axis I ppx ¼ I ppy 1 9 10 -2 -2 9 10 -2 kg m 2 Clutch disc mass m cd 0.7-1.8 kg Clutch disc moment of inertia around Z-axis I cd 4 9 10 -3 -6 9 10 -3 kg m 2 Clutch disc moment of inertia around X, Y axis I cdx ¼ I cdy   the acceleration magnitude is kept almost at the same level due to the increase in the clamp load. From 4 s onwards, the clutch pedal position is kept steady and decreasing amplitude is observed which leads to magnitudes close to 4 m/s 2 (peak to peak) at the end of the manoeuvre. Figure 10b shows the numerically predicted acceleration of the input shaft at the same location, along the same direction. The two sets of data match satisfactory, in terms of size and trend during the examined manoeuvre. The oscillations start from around 6 m/s 2 with an increasing rate until 0.5 s. From 0.5-4 s the magnitude of the acceleration is around 10 m/s 2 and when the clutch pedal is kept in the same position (4-7.2 s) the oscillations have a decreasing rate, leading to absolute magnitude of 4 m/s 2 .
In order to compare the corresponding frequency content, continuous wavelet transformations are shown for both signals in Fig. 12. The experimental measurements presented in the wavelet of Fig. 12a show clearly two natural frequencies during the 7.2 s measurement period and one natural frequency which is at the same level with the lowest one. The natural frequency near 400 Hz corresponds to the clutch tilting motion along Y-axis. A secondary frequency is shown at around 270 Hz corresponding to input shaft bending dof, this natural frequency slowly increases during the 5 s period. Figure 12b shows a wavelet of the same motion as captured by the numerical model, with similar frequency content. In both wavelets the intensity of the natural frequency starts decaying after 5 s which is the time at which the clutch pedal position is fixed.
Due to the axisymmetric design of the clutch system, the acceleration of the input shaft in X direction has similar magnitude and frequency content with that in Y direction. Thus, vibration data for only one of these directions is presented hereafter. In order to examine the system stability, the eigenvalue problem is solved for the complete clutch model. Figure 13 depicts the natural frequencies of the key modes that contribute to the system instability during the engagement manoeuvre presented above and the clutch pedal travel variation from Fig. 9. The cushion deflection step for the stability analysis was 0.01 mm. The instability region starts after 0.13 mm cushion deflection and lasts until 0.2 mm cushion deflection. It is noted that the wavelet frequency content for both the experiment and numerical model contains the predicted unstable mode. Moreover, the eigenvalues corresponding to the clutch disc (radial motion) and the input shaft bending motion have symmetric real parts (Fig. 13). This is an indication of mode coupling occurrence. A detailed analysis of a complete pull away manoeuvre is presented in the following section.

Results and discussion
In this section the oscillations of every clutch component will be investigated to examine both the frequency content and stability of the vibrations during the pull away manoeuvre. The relative angular velocities of the flywheel and input shaft are shown, where time zero corresponds to the first contact of the pressure plate and flywheel to the clutch disc, referred to as the 'kiss point' [43]. At this instant it is assumed there is a nominal cushion spring load of 70 N and as a result friction torque is generated between the clutch components. A typical initial flywheel velocity of 1570 rpm and the engine torque profile shown in  The friction coefficient variation during the engagement manoeuvre is shown in Fig. 15. After 1.51 s (synchronization time between the clutch disc and the flywheel), only the static friction is applied in the system. However, the pressure plate travels until 0.6 mm cushion deflection is achieved at 1.8 s. The rate that the cushion deflection is achieved is steady and is equal to t/3. Based on the experimental measurements of the friction coefficient, it can be noted that as the relative velocity between the flywheel/pressure plate and the clutch disc decreases the coefficient of friction increases. Between 0.4 and 0.47 mm of cushion spring deflection there is a slight drop in the friction coefficient. Within this deflection range, the spring load is between 1000-3000 N and the relative velocity is between 1200-400 rpm.
Simulation results for the pressure plate oscillations in the axial direction (z pp ) are shown in Fig. 16a. The amplitude of the oscillations decreases with time, eventually reaching the equilibrium position. After  1.8 s the clutch is fully closed and the pressure plate cannot move any further in the axial direction because the cushion spring is fully compressed [47]. The two frequencies contained in this motion are clearly shown in Fig. 16b. The primary frequency starting at 100 Hz corresponds to the axial motion of the pressure plate [48] and the secondary natural frequency which appears at 230 Hz is a contribution of the pressure plate tilting.
The time history and wavelet for the axial motion of the clutch disc are shown in Fig. 17a, b, respectively. The frequency at 120 Hz corresponds to the axial motion of the clutch disc and the secondary frequency   The pressure plate tilting motion about the Y-axis is shown in Fig. 18a. The oscillations occur about an equilibrium position of 3 9 10 -3 rad due to the circumferential variation prescribed to the cushion and diaphragm spring forcing. The tilting motion around Y has similar behaviour. The corresponding frequency content presented in the wavelet of Fig. 18b shows contributions from natural frequencies between 200-300 Hz, which corresponds to the tilting of the pressure plate and a second frequency centred at 400 Hz, which is due to the tilting of the clutch disc. Senatore et al. [28] reported increasing natural frequency for the tilting of the pressure plate, starting at 270 Hz.
The clutch disc oscillatory motion in radial direction (Y-axis) is shown in Fig. 19a. The corresponding wavelet is in Fig. 19b, where three frequencies are shown. The lowest starts at around 200 Hz and contributes until 1.2 s whereas the other two start at 230 Hz (pressure plate tilting mode) and 260 Hz, the first contributing until 0.3 s and the second until 0.6 s. The one starting at 200 Hz corresponds to the clutch disc radial motion. The highest, at 260 Hz, is a contribution from the input shaft bending mode. It will be shown further after conducting stability analysis that these two motions create a mode coupling instability.
The time history of the clutch disc tilting motion about the Y-axis is shown in Fig. 20a. The pressure plate equilibrium position is at 3 9 10 -3 rad due to the variation in circumferential cushion and diaphragm spring forces. After 1.2 s the vibrations have decreased significantly and only small oscillations around the equilibrium position remain. The corresponding frequency content is shown in Fig. 20b. Three main natural frequencies can be observed: the first starting at 230 Hz corresponds to the input shaft bending motion and the other two are contributions from the clutch disc tilting motion (around 400 Hz) and the clutch disc axial motion (around 100 Hz).
The bending motion of the input shaft around the Y-axis is shown in Fig. 21a. The frequency content of   Fig. 21b. Between 0-0.2 s a contribution from the clutch disc tilting motion (natural frequency around 380 Hz) can be observed. The other frequency contributions between 0-1.5 s correspond to the input shaft bending motion (starting at 260 Hz), and the clutch disc radial motion (200 Hz). It is shown later in the stability analysis between 0.4-0.6 s that a mode coupling instability at around 300 Hz occurs.
The last task is to conduct a stability analysis of the clutch system oscillatory behaviour during engagement. The real parts of the eigenvalues as well as the natural frequencies of the unstable modes are shown in Fig. 22a, b, respectively. The analysis is conducted for   Input shaft bending motion around Y axis and clutch disc radial motion along Y axis during the manoeuvre a Natural frequencies and b Real parts of the corresponding eigenvalues 1.8 s, which reflects the total duration of the manoeuvre until the clutch locks up. Between 0-0.4 s the frequencies corresponding to the bending motion of the input shaft around the Y axis (h insh;y ; mode1Þ and the clutch disc radial motion along the Y axis (Y cd;y ; mode2Þ are distinct. The real parts of these two modes shown in Fig. 22b are negative (stable modes) until 0.4 s. When the natural frequencies coalesce the two modes become coupled, which results in system instability. They reach a Hopf bifurcation point at 0.4 s [40,49]. The natural frequency of the input shaft bending motion increases at a higher rate than that of mode2, leading to modal decoupling at 0.6 s.
In order to further explore the behaviour during the unstable region identified in Fig. 22, limit cycles for the stable region 0.9-1.1 s and the unstable region 0.4-0.6 s are presented in Fig. 23. In this work the nonlinearities of the friction force, cushion spring stiffness and diaphragm spring stiffness lead to the Hopf bifurcation of the system between 0.4-0.6 s of the engagement manoeuvre under examination. The limit cycle of the unstable region shows greater amplitudes and more erratic behaviour compared to the limit cycle of the stable region (Fig. 23). In the limit cycle of the unstable region it can be noted that the motion is more erratic compared to the limit cycle of the stable region. The size of the former is larger, and the cycle is asymmetric in contrast to the stable cycle. Between 0.4-0.6 s the slipping speed between the pressure plate/flywheel and the clutch disc is higher and similar limit cycle behaviour is observed by Fidlin et al. [31].
A second stability analysis with different system parameters was conducted to explore other mode coupling occurrences that may lead to instability due to design variants of the system. The pressure plate mass and clutch disc mass were reduced by 54% and 36%, respectively. This reduction in the mass of the forementioned components also affects the inertia and the thickness of pressure plate and clutch disc [50]. Those limits were selected in order to be in agreement with experiential design rules associated with pressure plate and clutch disc inertia and thickness. The real parts of the eigenvalues and corresponding natural frequencies are presented in Figs. 24 and 25. Three areas of interest relating to the coupling of vibration modes are highlighted.
The first unstable region occurs between 0.3-0.4 s during which two modes are coupled at 349 Hz. Mode 1 starting at 349 Hz corresponds to the input shaft bending motion, whereas mode 2 starting at around 300 Hz corresponds to the clutch disc radial motion. After 0.4 s the two natural frequencies are again uncoupled (with the real parts of the eigenvalues becoming negative).
The second unstable region occurs at 0.8 s as a result of coupling between modes 3 and 4 (which correspond to the tilting motions of the clutch disc around Y and X axis). These modes coalesce at 0.8 s before uncoupling at 1.1 s leaving the system in a stable state.
Finally, the natural frequency of mode 4 increases and at 1.2 s it couples with the natural frequency of mode 5, which corresponds to the tilting motion of the pressure plate. Thus, the system becomes unstable until 1.5 s, when the clutch components are synchronized (locked clutch).

Conclusions
A novel fully coupled thirteen-degree of freedom dry clutch transient dynamic model has been presented to simulate the oscillations generated during the engagement phase. Coupling between the torsional, radial, bending and axial motions of the clutch components is established and studied. The inclusion of these motions has not been hitherto reported in clutch dynamics open literature other than involving complex finite element models of the clutch components/assembly. The natural frequencies of the system are determined, as well as the oscillatory behaviour of the components during the engagement manoeuvre. The clutch model results successfully identify unstable regions for different input clutch parameter combinations. Thus, it is a useful tool to examine NVH associated issues and lead towards clutch design with smooth NVH performance.
The analysis of a typical engagement manoeuvre has revealed mode coupling between the input shaft and the clutch disc radial motions, leading to system instability near 300 Hz. The self-excited oscillations in the unstable regimes are discussed after presenting the corresponding limit cycles. This is the result of the nonlinear forces applied to the clutch disc and the generated friction force among the contact surfaces. However, the amplitude of the oscillations did not increase significantly due to the fast transition from instability to stability, as well as the action of damping. The axial degrees of freedom present natural frequencies at around 100 Hz, which is also in qualitative agreement with the open literature. Altering the geometrical parameters of the clutch components, has led to previously unidentified occurrences of system instability due to mode coupling. Pressure plate tilting modes coupled with clutch disc's tilting modes can introduce instabilities to the dry clutch system. Although such mode couplings have been reported in the literature, no study has considered a combination of all the possible clutch component motions.
The future work will consider factorial analysis of each clutch component parameter to stipulate the effect that potential changes have on NVH intensity.
In which the terms of the matrices above are described by the following equations: