Investigation of complex nonlinear dynamic behaviors observed in a simplified driveline system with multistage clutch dampers

The results of the harmonic balance method (HBM) for a nonlinear system generally show nonlinear response curves with primary, super-, and sub-harmonic resonances. In addition, the stability conditions can be examined by employing Hill’s method. However, it is difficult to understand the practical dynamic behaviors with only their stability conditions, especially with respect to unstable regimes. Thus, the main goal of this study is to suggest mathematical and numerical approaches to determine the complex dynamic behaviors regarding the bifurcation characteristics. To analyze the bifurcation phenomena, the HBM is first implemented utilizing Hill’s method where various local unstable areas are found. Second, the bifurcation points are determined by tracking the stability variational locations on the arc-length continuation scheme. Then, their points are defined for various bifurcation types. Finally, the real parts of the eigenvalues are analyzed to examine the practical dynamic behaviors, specifically in the unstable regimes, which reflect the relevance of various bifurcation types on the real part of eigenvalues. The methods employed in this study successfully explain the basic ways to examine the bifurcation phenomena when the HBM is implemented. Thus, this study suggests fundamental method to understand the bifurcation phenomena using only the HBM with Hill’s method.

System modeling with piecewise type nonlinearities. The principle of piecewise-type nonlinearities, marked with a red dotted line in Fig. 1b, are depicted in Fig. 2. The basic profile of piecewise-type nonlinearities induced by the multi-staged clutch dampers is shown in Fig. 2a, where k Ci (i = 1, 2, 3, 4) are the stages of stiffness. In general, the clutch force is affected by various factors, and the input torque is transferred from the engine to the rest of the driveline system 6 . First, the clutch torque T S induced by multi-staged linear springs is considered, as shown in Fig. 2b, where φ p1 and − φ n1 are the transition angles on the positive and negative sides, respectively. Second, another clutch torque T H is induced by the dry friction between the clutch plate and friction materials, as depicted in Fig. 2c. Third, the preload, T SPr , must also be considered based on the design concepts along with various practical conditions. As Fig. 2b-d show only the two stages of relevant components to clutch forces, more stages of clutch force components will be included for the practical system. The other parameters will be explained while deriving the mathematical description below.
For the employed piecewise-type nonlinearities, the nonlinear force f n θ f ,θ f (or T C ) can be developed by mathematical formulation, as depicted in Fig. 2a 6,25 . First, the clutch torque T S θ f is derived from the stiffness with a smoothing factor σ C employed as 1 × 10 3 , as follows: (1) T sn(i) = θ f + φ n(i) tanh σ C θ f + φ n(i) − 1 . www.nature.com/scientificreports/ Here, 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 by a smoothing factor σ H of 0.1.
Here, 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 the torque calculated using Eqs. (1)-(6), the preload T Pr must be considered as a function of θ 1pr .
Here, T SPr is the total clutch torque induced by the pre-load; T Pr1 (or T Pr2 ) is the positive (or negative) torque induced by the pre-load; and φ Pr is the angle located at the pre-load. Overall, the total clutch torque is estimated by the summation of T S θ f , T H θ f ,θ f , and T SPr θ 1pr from Eqs. (1)-(7), as follows: The properties of the clutch torque f n θ f ,θ f (or T C ) employed are listed in Table 1. Finally, the practical clutch torque profile can be calculated using Eqs. (1)- (8), as illustrated in Fig. 2a.
Basic equation and development of HBM. The basic equation of motion for the torsional system, which includes the nonlinear function f n θ f ,θ f shown in Fig. 1b, is derived as follows: (8) f n θ f ,θ f = T C θ 1pr ,θ 1pr = T S θ 1pr + T H θ 1pr ,θ 1pr + T SPr θ 1pr .  www.nature.com/scientificreports/ Here, T E (t) and T D are the sinusoidal input and drag torques, respectively. In this study, the input torque is given by Fourier coefficients based on 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 with the harmonic index of the HBM. The input torque profiles employed are listed in Table 2.
In this study, the drag torque is assumed as T D = T m under steady-state conditions. To find out the nonlinear dynamic behaviors, this study will employ the harmonic balance method (HBM) since the HBM generally give the frequency and time domain solutions simultaneously under the steady-state and frequency sweeping conditions. Also, while the simulation is conducted, stability conditions can be efficiently determined by employing the Hill's method 5,6,27 . The Galerkin scheme of Eq. (10) is expressed as follows: Then, its corresponding formulae are defined as follows: sin(kψ 0 ) · · · 1 · · · cos(kψ 1 ) sin(kψ 1 ) · · · . . . . . .
1 · · · cos(kψ N−2 ) sin(kψ N−2 ) · · · 1 · · · cos(kψ N−1 ) sin(kψ N−1 ) · · · www.nature.com/scientificreports/ Likewise, its nonlinear and input functions are defined as follows: The relevant variables used are ωt = ψ and ω = ω ω n , where ω n is the non-dimensionalized time scale and normalized frequency value with the natural frequency; T = ητ is the concerned time period with respect to 0 ≤ t < T → 0 ≤ ψ < 2π ω n ; η is the sub-harmonic index; τ is the fundamental excitation frequency; and k is the incremental index as k = ω n , 2ω n , 3ω n . . . . By employing the relationship between θ (t) = dθ dt = ω dθ dψ = ωθ ′ and θ (t) = ω 2 θ ′′ , the overall Galerkin scheme of the basic equation of Eq. (11) is expressed as follows: To determine the solutions for θ c and ̟ for each step, the Newton-Raphson method is implemented under the condition → 0 , where is considered as a function of θ c and ̟ such that θ c , ̟ . A detailed explanation about these parameters and their use can be found in previous studies 6 . Initial results with frequency upsweeping condition Figures 3 and 4 show the initial results obtained by the HBM with η = 2 and N max = 24 under the frequency upsweeping condition. Figure 3 compares the results based on two different methods such as the HBM and the numerical simulation (NS), which reflects good correlations except for the sub-harmonic resonance area. Moreover, the Hill's method is employed to determine the stability condition as shown in Fig. 4 5,6,17 . For this study, the valid components of the input torque vector with η = 2 and N max = 24 are given as F m = 168.9 , F a(2) = −87.97 , and F b(2) = 235.65 in Eq. (18). Here, the employed values, η =2 and N max = 24 were enough to project the nonlinear effect when the simulations from the HBM are compared with the numerical simulation (NS) 5,6,27 . In general, the employment of values for η and N max differs depending on the various types of nonlinearities. Also, the other issues such as a convergence and number of index numbers for η and N max can be referred 5,6,27 .
To examine the sub-harmonic responses, the components of the input torque pertaining to the sub-harmonic locations, such as F a(1) and F b(1) , are F a(1) = εF a(2) and F b(1) = εF b (2) . Here, ε is taken as equal to 1 × 10 −5 . To obtain the sub-harmonic responses, a small range of excitation input values are artificially employed 18 . Figure 4a shows the simulated results, including both super-and sub-harmonic resonances, marked with red dotted circles. In addition, the resonant areas of (b) and (c) in the red dotted circles are more concentrated in Fig. 4b,c, respectively. When the super-and sub-harmonic resonances are examined, the stability variational conditions, such as stable to unstable (STU) and unstable to stable (UTS) conditions, are clearly observed. In addition, the STU and UTS conditions are simultaneously considered with the direction of the arc-length solution, as indicated by (1) and (2) in Fig. 4b,c. For example, the regimes indicated as (1) in Fig. 4b,c show that the directions of solutions are changed from facing the same direction as frequency upsweeping to facing the opposite direction under STU or UTS conditions.
T . www.nature.com/scientificreports/ However, the regimes indicated as (2) always follow the same direction as frequency upsweeping for both STU and UTS conditions. Based on prior studies, the regimes marked as (1) are expected to clearly show that the jumping phenomenon occurs; the phenomenon is normally known as the saddle-node points that are observed at slightly different locations based on the frequency upsweeping or down-sweeping conditions. In general, a representative saddle-node bifurcation is clearly observed in the primary resonant area, as shown in Fig. 4a. However, regarding the regimes (2) with unstable conditions, as shown in Fig. 4b,c, it is difficult to anticipate specifically whether the simulated results under STU occur in the practical system or they just project theories from the HBM using Hill's method. Thus, more techniques must be considered to understand the nonlinear dynamic behaviors with respect to the bifurcation characteristics.

Examination of primary, super-, and sub-harmonic resonances
As a first step to examine the nonlinear dynamic behaviors, different resonant areas such as primary, super-, and sub-harmonic resonances are determined, as shown in Fig. 5. To define the various resonant regimes, basic linear concepts for the half-power method can be employed 26 . For example, the relationships between the Q factor location and the half power location lead to the following basic equations:  www.nature.com/scientificreports/ Here, ω n , ω 1 , and ω 2 are the natural frequency and the lower and higher values of half power locations for the linear system, respectively. In addition, ω i (with i = 1, 2) is the normalized frequency value of ω n . Thus, this study employs the same concepts to define different resonant areas based on Eqs. (22) and (23). Figure 5a illustrates the comparisons of the HBM results with various damping ratios using backbone curves. Here, B ω is the basic backbone curve, which is estimated under no input or damping conditions 27 . In addition to the basic B ω curve, two more backbone curves transferred from B ω along the ω axis can be calculated as follows: Figure 5a shows three cases of HBM results for ζ = 0.03, 0.05, and 0.1, via a comparison of the backbone curve B ω and the transferred backbone curves. Here, assuming that the concepts in Eq. (20)-(23) for the linear system can be approximately employed for the nonlinear system, the cross points of B ω1 and B ω2 with respect to the HBM result with ζ = 0.05 are considered as ω 1 and ω 2 as indicated in Fig. 5a,b. Based on Eq. (24), the local resonances are obtained, which are normally defined as primary, super-, and sub-harmonic resonances, along with different frequency response regimes. Thus, when the bifurcation point (BP) is estimated by focusing on the frequency upsweeping condition in this study, the simulation can simultaneously recognize their relevant regimes for the primary, super-, and sub-harmonic responses by employing the values ω 1 and ω 2 calculated from Eq. (24). Figure 5b shows the different resonant areas referenced by ω 1 and ω 2 . For example, the areas marked with (1), (2), and (3) indicate the primary, super, and sub-harmonic resonances, respectively. In addition, the primary resonance is normally affected by the saddle-node BP, as shown in Fig. 5b.

Investigation of saddle-node bifurcation point
As a second step to determine the BPs, the locations where the stability conditions vary can be determined. In general, bifurcation occurs when the real parts of the eigenvalues are zero 28 . In addition, if the variational condition of arc-length directions at a specific point is examined, then the bifurcation characteristics in the next steps of the solutions can be determined in terms of saddle-node or other types of bifurcations. To determine the saddle-node BP, the direction of the arc-length vector can be estimated. Here, by assuming the physical dynamic behaviors at the saddle-node BP are directly correlated with the jumping phenomena, as clearly observed in Figs. 4 and 5, the angle at this point is deduced to be 90° under STU or UTS conditions. To calculate and determine the saddle-node point, the angle between the arc-length and unit vectors can be estimated. Figure 6 depicts the various cases of arc-length vectors and their angles with unit vectors.
Here, the unit vectors with the 2ηN max + 1 number of components along with the size of HBM matrices and vectors are defined as follows: Thus, the unit vectors are always in the horizontal direction, which are referenced to calculate the angles of the arc-length vector s i (with i = a1, b1, and a) from u ω or −u ω at each step, as indicated in Fig. 6. Here, θ a1 , θ b1 , θ a , and θ b are certain points of solution that determine the BPs. To calculate the angles constructed by each arc-length with the unit vectors, the basic concept of inner product can be used, as follows: While ϕ i is estimated, u ω or −u ω must be chosen based on the frequency up-or down-sweeping conditions, respectively. For the locations of θ a1 and θ b1 , the composed angles ϕ a1 (or ϕ b1 ) between s a1 (or s b1 ) and u ω (23)  www.nature.com/scientificreports/ (or −u ω ) are much less than 90°. However, other points, such as θ a and θ b show much stiffer variational changes in angles, then their angles approach 90°. For example, when the arc-length vector around θ a changes its direction from the frequency upsweeping to the opposite direction, s k a at the prior step finally reaches s k+1 a at the current step, where the composed angle between s k+1 a and u ω becomes 90°. Again, when the arc-length vector around θ b changes its direction from the opposite to the frequency upsweeping direction, s k b goes to s k+1 b , which is at 90° with u ω . Based on the given variational conditions, the angle index can be determined as follows: Because ϕ i is calculated numerically while the HBM is conducted, the angles obtained are approximately 90°. Thus, the absolute value |ϕ i | is used to determine the difference from 90°. Here, |ϕ i | is considered as approximately 90° if �ϕ < ε ϕ , where the index number ε ϕ is π/180 rad. Thus, if the condition �ϕ < ε ϕ is satisfied, this point is determined as a saddle-node point. Overall, the examination of �ϕ with respect to the stability conditionsnamely, STU and UTS-leads to two different bifurcation conditions, such as saddle-node and other types of bifurcations, as summarized in Table 3. Figure 7a clearly shows the saddle-node and other bifurcation points determined by Eqs. (25)- (27). In addition, the distribution of �ϕ under the frequency upsweeping condition is shown in Fig. 7b, where the saddle-node points indicated by red rectangles are located near 90°. Except for the saddle-node points, other BPs illustrated with green circles are located below or away from 90°; their relationships are described in Eq. (27). Figure 8a-c show the BPs in more detail by focusing on the primary, super-, and sub-harmonic resonant areas, as indicated by (1), (2), and (3) in Fig. 7a, respectively. From the results illustrated in Fig. 8, it can be observed that the saddle-node BPs are located at the jump phenomena points, whereas the other BPs follow the frequency upsweeping direction. In addition, the multiple saddle-node points found at the primary resonance are due to numerical problems during the HBM process, which will be further investigated in a future work.

Examination of eigenvalues for the bifurcation characteristics
The saddle-node BP is clearly revealed along with the estimation of the arc-length direction based on frequency upsweeping in this study, as shown in Figs. 5, 7, and 8. However, other BPs relevant to the period-doubling, quasiperiodic, and chaotic cascades are not clearly observed from the previous techniques. To analyze other types of bifurcations, the real parts of eigenvalues Re(λ) will be examined. To determine the relationships between Re(λ) and the bifurcation characteristics, the distribution of positive Re(λ) with respect to the overall number of Re(λ) is first examined. The positive Re(λ) is correlated with the unstable condition of the system responses. Second, a reference value of Re(λ), such as the maximum positive Re(λ), is investigated over the range of frequencies for the super-and sub-harmonic resonances. Figure 9 shows the two properties previously described. Figure 9a shows the distribution of positive Re(λ) values over all numbers of Re(λ), calculated as follows:         Fig. 13. However, E d and E max in regimes (C) and (D), where more complex dynamic behaviors are observed, show higher values, as shown in Fig. 12b,c. For instance, E d values in regimes (C) and (D) vary between 4.5 and 5. In addition, while the nonlinear dynamic behaviors are more complex, the level of E d reaches approximately 0.5. Moreover, E max increases when the dynamic behaviors change from period-doubling cascade to chaotic response, as shown in Fig. 13.

Conclusion
This study investigated nonlinear dynamic responses by examining bifurcation characteristics and their relevant properties. In particular, to anticipate complex dynamic behaviors, arc-length vectors and properties of eigenvalues were estimated by employing the HBM. The contributions of this study are as follows: First, this study  www.nature.com/scientificreports/ suggested the numerical techniques to determine BPs on the HBM solutions by estimating the direction of the arc-length vector. From the given results, the BPs were defined as either saddle-node BPs or other types of BPs. Second, the relationships between the bifurcation characteristics and the properties of the eigenvalue real parts, such as E d and E max , were investigated. Based on these results, while the bifurcation became more complex, E d and E max became higher than those in other areas.
There are several possibilities for future research based on the unfinished tasks in this study. For example, the primary resonance has multiple saddle-node points, as shown in Figs. 7a and 8a. In addition, the bifurcation phenomena above ̟ = 1.06 are still to be analyzed using the results shown in Fig. 4a, which were not fully examined here. These issues will be investigated as the next stage of our research.