Examination of sub-harmonic responses along with various initial conditions induced by multi-staged clutch damper system

Using the harmonic balance method to investigate the nonlinear dynamic behaviors pertaining to sub-harmonic responses is difficult compared with that of super-harmonic cases because of the limitations of the HBM. Since sub-harmonic motions differ under various initial conditions, difficulties can arise when this method is used to calculate all possible solutions within sub-harmonic resonances. To explore complex dynamic behaviors in sub-harmonic resonant areas, this study suggests mathematical and numerical techniques to estimate sub-harmonic responses depending on various initial conditions. First, sub-harmonic responses are calculated under various excitation conditions relevant to the sub-harmonic input locations of the HBM formula. Second, the HBM results are verified by comparing them with the results of the numerical simulation (NS) under various initial conditions with respect to different frequency up-sweeping paths. Finally, the positive real part of the eigenvalues is examined to anticipate bifurcation characteristics, which reflect the relevance of the complex dynamic behaviors in the eigenvalues’ unstable solutions. Overall, this study successfully proves that the techniques and methods described are suitable for examining complex sub-harmonic responses, and suggests basic ideas for analyzing nonlinear dynamic behaviors in sub-harmonic resonances using the HBM.

Nonlinear dynamic responses observed in a practical system, as illustrated in Fig. 1, reflect various types of dynamic phenomena, generally called super-and sub-harmonic, periodic, quasi-periodic, and chaotic. These responses are normally signified by various types of bifurcations in super-and sub-harmonic resonant areas. Complex nonlinear dynamic behaviors in sub-harmonic resonances are particularly difficult to find, relative to other resonances such as primary and super-harmonic areas, when the harmonic balance method (HBM) is employed.
The methods for solving complex nonlinear problems using the HBM have been reported for decades  . For example, Peng et al. implemented nonlinear output frequency response functions (NOFRFs) in strong nonlinear equations by applying the Volterra series to extend the classic frequency response function (FRF) to a nonlinear case 1 . Al-shyyab and Kahraman investigated sub-harmonic and chaotic motions in a multi-mesh gear train using a nonlinear time-varying dynamic model 2 . Here, nonlinear dynamic motions were simulated with the multi-term harmonic balance method and correlated with the direct numerical integration results. Chen et al. used the incremental harmonic balance (IHB) method to investigate the limit cycle oscillation of a two-dimensional airfoil with parameter variability in an incompressible flow 3 , and further utilized the IHB to estimate the strong nonlinear cubic stiffness, subject to either (1) non-probability but bounded uncertainty, or (2) bounded stochastic parameters. Kim et al. proposed and verified a multi-term harmonic balance method by including adaptive arc-length continuation and stability calculation capabilities 4 , which aided them in developing nonlinear frequency response calculations of a torsional system with clearance-type nonlinearity. Miguel et al. suggested the closed-form solutions for the Bouc-Wen and LuGre models by developing a smoothing procedure with the harmonic balance method 5 . The bifurcation analysis techniques with the employment of the harmonic balance method were suggested by Detroux et al. and Xie et al. 6,7 . The basic mathematical models of the harmonic balance method by using Galerkin forms were investigated by the prior researches [8][9][10][11] . Masiani et al. suggested masing model to analyze the hysteretic behavior of the elements with the multi-component harmonic balance method 12 . Raghothama 16 . To determine the stability of system responses, the Hill's method was suggested and employed 17,18 . Comparin and Singh suggested the driveline model embedded by the clearance type nonlinearities and investigated its complex system responses 19 . Sundararajan and Noah examined the nonlinear responses in rotor systems by employing shooting/arc-length continuation method 20,21 . However, some difficulties remain in capturing the complex motions of sub-harmonic responses, especially when their dynamic behaviors differ depending on various initial conditions. Thus, to investigate sub-harmonic responses, the sub-harmonic index must be used, and basic matrix formulae are constructed 8,18 . To capture sub-harmonic responses, further techniques must be utilized, such as numerical modification of the input conditions 18 . This study expands upon the theoretical and numerical methodology of prior studies which could not clearly determined the sub-harmonic responses, by investigating the nonlinear dynamic behaviors in subharmonic resonant areas sensitive to initial conditions using the HBM 8 . The specific objectives are as follows: first, to investigate nonlinear dynamic motions under various input excitation conditions when the sub-harmonic input components of the HBM are artificially defined; second, to verify sub-harmonic responses of the HBM by comparing the numerical simulation (NS) under various initial conditions and different frequency up-sweeping paths; and finally, to examine the eigenvalues' unstable solutions for bifurcation characteristics, which will provide a guide for anticipating bifurcation characteristics while the stability conditions are determined simultaneously based on the Hill's method. Therefore, understanding the dynamic behaviors at sub-harmonic resonance areas will give the reasonable insights to resolve the severe vibrational problems such as gear rattle under the vehicle coast conditions 8,22 . To analyze these main issues, this study focuses on a torsional system with one degree-offreedom (DOF) affected by piecewise-type nonlinearities, as presented in Fig. 2. This system has been simplified from the physical system shown in Fig. 1.  www.nature.com/scientificreports/

Physical system and its problem formulation
Physical system and its modeling. Figure 1 illustrates a physical driveline system based on the frontengine and front-wheel (FF) layout 22 . The entire system can be considered to have three main parts: the engine, which is the power generating system; the multi-staged clutch dampers, which operate as the vibration isolator; and the transmission and driveline subsystems, which transfer the torque from the engine into the wheel.
Assuming the inertia value of the transmission and driveline systems' lumped mass is relatively higher than the engine inertia, the physical system can be simplified to a single DOF system, as illustrated in Fig. 2a 8,22 . The multi-staged clutch dampers have the properties listed in Table 1 8,22 , and are the main nonlinearities for the system, whereas the transmission and driveline systems are assumed as the ground. Using the values provided in Table 1 which were measured from the clutch damper experimental setup, the clutch torque profile is drawn, as shown in Fig. 2b. Therefore, based on prior studies, the nonlinear one DOF torsional system can be considered a part of the driveline by focusing on the rotational motion observed in the crankshaft and flywheel's lumped mass with multi-staged clutch dampers 8,22,23 . In this study, the inertia of the flywheel and crankshaft's lumped mass, I f = 1.38 × 10 −1 kg m 2 , and viscous damping, c f = 1.59 N m s/rad 8 were the parameters implemented for the torsional system shown in Fig. 2a, and the input torque, T E , drag torque, T D , and clutch torque, f n θ f ,θ f define the additional parameters. Here, θ f and θ f are the angular displacement and velocity of the flywheel (subscript f), respectively, as indicated in Fig. 2a. Here, the employed system parameters are measured and given based on the manual transmission type of driveline system 23 .
Development of a mathematical model with piecewise-type nonlinearities. The equation of motion for the one DOF system shown in Fig. 2a is derived as follows: Here, f n θ f ,θ f is the nonlinear function, which will be explained in terms of piecewise-type nonlinearities. T E (t) and T D are the sinusoidal input and drag torque, respectively. The excitation of the system in terms of its input torque is given by Fourier coefficients from the measured data as follows: Here, T m and T pi are the mean and alternating parts of the input torque, respectively; ω p and ϕ pi are the excitation frequency and phase angle, respectively; and N max is the maximum number of harmonics correlated to the harmonic index of the HBM. The input torques are described using the properties listed in Table 2 and, assuming the system is under steady-state conditions, the drag torque can be expressed as T D = T m , also the employed number of alternating part of toque is 1 with T pi (i = 1). Here, the torque profile employed as shown in Table 2 is measured from the engine dynamometer setup under the WTO (wide open throttle) conditions based on the 4-cylinder engine 23 .
In addition, the clutch torque is dependent on various factors, while the input torque is transferred from the engine into the rest of the driveline system 8 . The basic profile of piecewise-type nonlinearities induced by multistaged clutch dampers is shown in Fig. 2b, where k Ci (i = 1,2,3,4) indicates the stage of stiffness. The main factors to consider for the clutch torque are the clutch torque T S from multi-staged linear springs, a second clutch torque T H from the dry friction between the clutch plate and friction materials, and the total preload effect T SPr due to T pi cos iω p t + ϕ pi Table 1. Properties for the piecewise type nonlinearities based on the practical system 8,23 .

Property Stage Value
Torsional stiffness, k Ci (Linearized in a piecewise manner) (Nm/rad) www.nature.com/scientificreports/ design concepts under various practical conditions. The mathematical model of nonlinear torque f n θ f ,θ f (or T C ) can be derived from prior studies 8,22,23 . First, the clutch torque T S θ f from the stiffness with a smoothing factor σ C of 1 × 10 3 is described as follows: where k C(N) (or k C(i) ) is the Nth (or ith) stage of the clutch stiffness (with subscript N or i), T sp(i) (or T sn(i) ) is the positive (or negative) side of the clutch torque induced by the stiffness at the ith stage (with subscript p or n), and φ p(i) (or − φ n(i) ) is the ith transition angle of the positive (or negative) side. Second, the clutch torque T H induced by dry friction is derived with a smoothing factor σ H of 0.1.
where H N (or H (i) ) is the Nth (or ith) stage of hysteresis (with subscript N or i), and T Hp(i) (or T Hn(i) ) is the positive (or negative) side of the clutch torque induced by hysteresis at the ith stage (with subscript p or n). In addition to calculating the torque using Eqs. (3) and (4), the preload T Pri (subscript i = 1(or 2) for the positive (or negative) value) must be considered as a function of θ 1pr .
Here, T SPr is the total clutch torque induced by the preload, T Pr1 (or T Pr2 ) is the positive (or negative) torque induced by the preload, and φ Pr is the angle located at the preload. Thus, the total clutch torque is estimated by the summation of T S θ f , T H θ f ,θ f , and T SPr θ 1pr from Eqs.  1) can be expressed as follows 8,22 : Its relevant terms are defined as follows: www.nature.com/scientificreports/ Further, the nonlinear and input functions are defined as follows: Within these equations, the relevant variables include ̟ t = ψ , the nondimensionalized time scale; ̟ = ω/ω n , the normalized frequency value with a natural frequency of ω n ; T = ητ , the concerned time period with respect to 0 ≤ t < T → 0 ≤ ψ < 2π/ω n ; η , a sub-harmonic index; τ , a fundamental excitation frequency; and k, the incremental index defined by k = ω n , 2ω n , 3ω n . . . . By employing the relationship between θ (t) = dθ dt = ̟ dθ dψ = ̟ θ ′ and θ (t) = ̟ 2 θ ′′ , the overall Galerkin scheme of Eq. (7) is expressed as follows: To determine the solutions of θ c and ̟ for each step, the Newton-Raphson method was implemented under the condition → 0 , where is considered as a function of θ c and ̟ such that � θ c , ̟ . Prior studies can be referred to for more derivation and descriptions on the HBM 8 . Figure 3 shows a comparison of the HBM and NS which was also studied in the prior study 8 . The HBM was conducted with η = 2 and N max = 12 under the frequency up-sweeping condition, Hill's method was utilized to determine the stability condition 4,8,17,24 , and the valid components of the input torque vector were expressed as F m = 168.9 , F a(2) = −87.97 , and F b(2) = 235.65 in Eq. (9d) 8,22,23 . Figure 3 also presents a correlation between the simulations resulting from the HBM and NS; however, the outcome of both at the sub-harmonic region marked with a red dotted line, as shown in Fig. 3b, are not matched with each other; the NS results simulate the sub-harmonic resonance well, while the HBM follows the normal harmonic response path with an indication of instability.

Investigation of sub-harmonic responses under the frequency up-sweeping condition.
Based on a previous study 18 , numerical techniques must exist to determine sub-harmonic resonant behaviors; hence, to reveal these sub-harmonic responses, the input torque values pertaining to the sub-harmonic locations, such as F a(1) and, F b(1) were given by F a(1) = εF a(2) and F b(1) = εF b (2) . Here, 1 × 10 −5 is artificially employed for ε to avoid triggering the system responses corresponding to the normal input torque. Figure 4a shows that the simulated results successfully included sub-harmonic resonances marked with blue dotted circles. In general, www.nature.com/scientificreports/ the unstable conditions of the HBM are closely related to complex dynamic behaviors, such as super-or sub-harmonic resonances and bifurcations. However, the stable conditions always indicate responses from the practical system that follow their relevant excitation harmonics such as in Fig. 4b, where a specific region under the stable condition is observed around the concentrated sub-harmonic resonance area at ̟ = 1.22 . To examine its nonlinear dynamic behaviors, a NS was conducted by artificially applying the initial conditions, θ f (0) = −1.89(rad) and θ f (0) = −1.04(rad/s) , with the assumption that these conditions can be abruptly changed physically. The blue rectangles and the green diamonds indicated in Fig. 4b show the frequency down-and up-sweeping results, respectively, from the stable region in the sub-harmonic resonance at ̟ = 1.22 . Furthermore, the areas marked with (1) and (2)  To capture all the relevant sub-harmonic resonances in this regime, the appropriate input torque values for the sub-harmonic components in Eq. (9c) can be artificially obtained using various values of ε . To trigger their subharmonic resonances, ε must be increased within a certain range of values assumed to be nearly zero numerically. Figure 5 compares the three results of the HBM using 1 × 10 −5 , 1 × 10 −3 and 2.9 × 10 −3 as values of ε . Assuming the baseline is the HBM simulation with ε = 1 × 10 −5 , the results with ε = 1 × 10 −3 and ε = 2.9 × 10 −3 follow the nonlinear system response of the baseline well. However, as displayed in Fig. 5, the simulation with ε = 2.9 × 10 −3 was the only one to successfully capture all sub-harmonic resonances. Figure 6a shows the stability conditions with the arc-length continuation schemes, and Fig. 6b compares the two results from the HBM and NS. The latter are calculated using the initial conditions at ̟ = 1.22 of the upper sub-harmonic resonance regime and three different solution paths are presented, with paths (A), (B), and (C) defined in this study as the upper, lower, and normal paths. Complex dynamic behaviors in sub-harmonic resonant areas can be analyzed using information shown in Fig. 6b. First, as indicated in Fig. 3, the NS with frequency up-sweeping from ̟ = 0.93 and frequency down-sweeping from ̟ = 1.7 , follows path (C). Second, the nonlinear solutions on paths (A) and (B) can be obtained by utilizing each stable regime's particular initial conditions, as shown in Fig. 6a. Third, the dynamic behaviors between ̟ = 1 and ̟ = 1.7 are sensitive to various initial conditions; hence, the system responses can be abruptly changed among paths (A), (B), and (C). The initial conditions of a physical system can be determined depending on various external factors, such as sudden changes in drag torques and velocities due to fluctuations in road conditions. For example, vibro-impacts such as gear rattles occur with different driving conditions 23 . While gear rattles generally arise around acceleration ranges of 1800 rpm, the same vibro-impacts  www.nature.com/scientificreports/ are also found at higher velocity ranges of 3200 rpm when the system is under coast conditions. This, as shown in the nonlinear resonant issues in Fig. 6, demonstrates that the same dynamic behaviors can occur under different initial conditions. Examination of bifurcation characteristics with the real part of eigenvalues. As shown in Figs. 3, 4, and 6, the stability conditions, which were determined during the HBM, are examined, and the real parts of the eigenvalues are considered. Applying the basic procedures of Hill's method, the system responds as follows.
Here, θ f (t) and ξ f (t) are the periodic (particular) and perturbation (homogeneous) parts of the solutions, respectively 4,8,22 , and ξ f (t) consists of p(t) and e t , the periodic and decay terms, respectively. Based on the calculated eigenvalues i with i = 1, 2, . . . , 2(2ηN max + 1) , at least one positive value among Re( i ) makes the system unstable. In addition,Re( i ) = 0 satisfies the condition that the system responses fall into bifurcation, as indicated in Fig. 8a 7,24 . From this point onward, this study assumes that the properties of the positive Re( i ) signify the bifurcation characteristics, and by focusing on the frequency up-sweeping direction, this bifurcation is assumed to occur in the regimes determined as unstable. Thus, only these unstable regimes will be examined, with the exception of the arc-length solutions' reverse direction. Though various conditions of Re( i ) are investigated, they are generally considered in three forms: Re( i ) < 0 , Re( i ) > 0 , and Re( i ) ≫ 0 . Under the condition Re( i ) < 0 , its relevant transient responses ξ f (t) reach zero, which stabilizes the system response and sustains the particular solution, θ f (t) = θ fp (t) as depicted in Fig. 8b. However, Re( i ) > 0 and Re( i ) ≫ 0 cause the system to fall into diverging conditions by preventing the system responses from maintaining that particular solution, as illustrated in Fig. 8c and d. Re( i ) > 0 , rather than Re( i ) ≫ 0 , is assumed to gradually diverge the system response from the periodic solution of θ f (t) = θ fp (t) ; hence, it is assumed that Re( i ) > 0 is relevant for period-doubling or period-halving cascades. On the other hand, it is assumed that Re( i ) ≫ 0 allows the system responses to change rapidly into more complicated phenomena, such as quasi-periodic and chaotic.
The relationships between each unstable solution of the HBM and the properties of Re( i ) are shown in Fig. 9. To implement the previously suggested ideas, two properties, E d the distribution of positive Re( i ) and    www.nature.com/scientificreports/ is examined, the value of E d is higher than 0.45 and the value of E max is higher than 5, and the system responses are much more complex with quasi-periodic or chaotic motions. Moreover, the system responses are smoothly connected between regimes (a) and (c), even though it follows the arc-length path (A), as indicated in Fig. 6b.
Since the two stable regions in (a) and (c) are located near the unstable responses of (b), instead of obeying the arc-length path (A), the dynamic behaviors occur by following the order of: (1) a stable response in regime (a), (2) an unstable response in regime (b), and (3) a stable response in regime (c). Due to a series of smooth changes in the stability conditions from regimes (a) to (c), regime (c) is expected to be under pure harmonic responses, and thus, the calculated results of E d and E max can only be ignored for regime (c). Additionally, the characteristics of the bifurcation diagrams observed in Figs. 10, 11, and 12, also correlate well when focused on regimes (a)-(c) in particular. The green dotted circles and ellipses in regime (d) are the reverse path results of those for the lower branch of sub-harmonic resonance, which is not the main concern of this study. An examination of the calculated values E d and E max must be carefully carried out when investigating the nonlinear dynamic behaviors for regime (d), since these values are taken from two different paths (A) and (B), as shown in Fig. 6b. When the bifurcation characteristics shown in Figs. 11 and 12 are compared, clear differences in  www.nature.com/scientificreports/ the two system responses are observed depending on the initial conditions. For example, Fig. 11 shows the results when the initial conditions are set at ̟ = 1.22 on path (A), as shown in Fig. 4. From this point, the solutions follow the frequency down-sweeping direction, and first show the period-halving and period-halving cascades which fall into chaotic in regime (d). When the solutions follow the frequency up-sweeping direction, they enter the period-doubling cascade and finally reach the pure harmonic response regime. When these complex bifurcation phenomena are compared with Fig. 9b Fig. 7b; thus, the bifurcations only occur for the arc-length paths (B) and (C). The two bifurcation diagrams in Figs. 11 and 12, with the exception of the responses around the area at ̟ = 1.22 are correlated with each other and through careful examination of Fig. 9, 10, 11 and 12 the results of regimes (c) and (d) provide the basic ideas required to map the complex nonlinear dynamic behaviors. First, the system responses can be abruptly changed under various initial conditions, and should follow different solution paths, despite the potential for a smooth change in the initial conditions along the normal solutions path (C), as shown in Fig. 6. Second, the bifurcation characteristics can also be defined under various initial conditions and can be anticipated along with the relevant eigenvalues E d and E max in unstable conditions, as shown in Fig. 9.

Conclusion
This study investigated the nonlinear dynamic behaviors in sub-harmonic resonant areas under various initial conditions. Depending on these different initial conditions, the system responses followed different paths of the arc-length solutions from the HBM. Regarding the complex nonlinear phenomena, the contributions of this research are as follows. First, this study suggested a method to determine all possible sub-harmonic resonances by employing various excitation conditions, successfully simulating two branches of sub-harmonic resonances using the HBM with a variety of initial conditions. Second, this study investigated various paths of nonlinear responses under different initial conditions by employing the HBM in comparison with the NS. These results show the multitude of system responses, while the dynamic behaviors follow different arc-length solution paths. Finally, the bifurcation characteristics that occur in the unstable solutions were examined based on the eigenvalues in terms of E d and E max . This provided basic ideas for predicting the bifurcation characteristics numerically.
The current study has dealt with a simple torsional system with 1DOF which is a part of whole driveline system 23 . As a further study, the conceptual approach employed for this study can be extended to investigate the severe vibro-impact problems which occur in the driveline system concerned with the gear rattle or the gear whine.

Data availability
The datasets generated and/or analysed during the current study are not publicly available due to its relatively large size but are available from the corresponding author on reasonable request.