Stability and bifurcation analysis of super- and sub-harmonic responses in a torsional system with piecewise-type nonlinearities

The nonlinear dynamic behaviors induced by piecewise-type nonlinearities generally reflect super- and sub-harmonic responses. Various inferences can be drawn from the stability conditions observed in nonlinear dynamic behaviors, especially when they are projected in physical motions. This study aimed to investigate nonlinear dynamic characteristics with respect to variational stability conditions. To this end, the harmonic balance method was first implemented by employing Hill’s method, and the time histories under stable and unstable conditions were examined using a numerical simulation. Second, the super- and sub-harmonic responses were investigated according to frequency upsweeping based on the arc-length continuation method. While the stability conditions vary along the arc length, the bifurcation phenomena also show various characteristics depending on their stable or unstable status. Thus, the study findings indicate that, to determine the various stability conditions along the direction of the arc length, it is fairly reasonable to determine nonlinear dynamic behaviors such as period-doubling, period-doubling cascade, and quasi-periodic (or chaotic) responses. Overall, this study suggests analytical and numerical methods to understand the super- and sub-harmonic responses by comparing the arc length of solutions with the variational stability conditions.

In the investigation of nonlinear dynamic responses, the stability conditions can contain important information, which may induce dynamic variations ranging from simple periodic responses to complex dynamic behaviors, including period-doubling, quasi-periodic, and chaotic motions [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] . To examine the complex dynamic behaviors, this study employs the harmonic balance method (HBM) using Hill's method, which reveals the stability conditions of the system 16,17 . In addition, a numerical simulation (NS) was conducted to examine the system responses in the phase plane diagrams and Poincare maps. The comparisons of nonlinear motions based on the HBM with phase diagrams and Poincare maps focused on the variational stability regimes allow the physical dynamic behaviors to be more comprehensible.
With regard to the methods employed in this study, various works have been conducted previously. For example, Peng et al. suggested nonlinear output frequency response functions (NOFRFs) as an HBM using the Duffing oscillator 18 . To implement NOFRFs in strong nonlinear equations, they employed the Volterra series to extend the classic FRF to the nonlinear case. Chen et al. used the incremental harmonic (IHB) method to investigate the limit cycle oscillation of a two-dimensional airfoil with parameter variability in an incompressible flow 19 . Here, the strong nonlinear cubic stiffness subject to either non-probability but bounded uncertainty or bounded stochastic parameters was estimated using the IHB method. Duan et al. suggested an excitation perturbation method to investigate sub-harmonic resonance 20 . To capture the sub-harmonic effects, the authors modified the input conditions. For example, the relevant sub-harmonic input terms were artificially included, which triggered the corresponding sub-harmonic responses. Additionally, various approaches have been developed by employing a multiterm HBM to examine the nonlinear frequency responses in a torsional system with clearance-type nonlinearity [21][22][23][24][25][26][27] .
To advance the existing pool of knowledge based on the HBM and its relevant techniques, this study suggests a method to investigate the nonlinear dynamic characteristics that occur in the super-and sub-harmonic regimes, which mostly focus on the analysis of the stability conditions along the arc-length direction under the frequency upsweeping condition. Thus, the specific objectives of this study are as follows: (1) to investigate the dynamic characteristics in the super-and sub-harmonic areas with variational stability conditions using the HBM. The stability variations along the direction of the arc length can present the close relationships of the system response to various nonlinear motions such as period-doubling, period-doubling cascade, and chaotic which was not examined in the prior studies 17,28 ; (2) to analyze the nonlinear dynamic responses by comparing the bifurcations with the stability conditions. To understand the complex dynamic conditions, FFT results, phase diagrams, and Poincare maps are also compared with the bifurcations. To resolve these issues, this study focuses on a specific multi-staged clutch damper model in a torsional system with 1DOF 17,28 , which will lead to understanding the vibration issues caused in the vehicle driveline system and give the better guideline to reduce the annoying noise and vibrational problems.

Problem formulation and its basic formulation
System modeling and piecewise-type nonlinearities. Figure 1a depicts a nonlinear torsional system with piecewise-type nonlinearities. Here, the nonlinear system with 1DOF is a part of a driveline focused on the flywheel and multi-staged clutch dampers based on prior studies 17,28 . The effective nonlinearities are based on asymmetric torsional springs with dry friction materials. To analyze the nonlinear force f n θ f ,θ f (or T C ) induced by the piecewise-type nonlinearities, mathematical formulations can be derived by considering the hysteresis effects as illustrated with the red dotted arrows in Fig. 1b 17,28 .
Here,θ f and θ f represent the angular displacement and velocity of the flywheel (subscript f ), respectively, as indicated in Fig. 1a. First, the clutch torque T S θ f induced by the stiffness with a smoothening factor σ C of 1 × 10 3 is derived as follows: 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) ) denotes 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) ) represents the ith transition angle of the positive (or negative) side. Furthermore, the clutch torque T H induced by hysteresis is considered with a smoothening factor σ H of 0.1. Here, T SPr is the total clutch torque induced by the preload, T Pr1 (or T Pr2 ) denotes the positive (or negative) torque induced by the preload, and φ Pr represents the angle located at the preload. 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)-(3) as follows: In addition, this study employed relatively higher value for σ C and lower one for σ H based on the prior studies 17,28 since the linear torsional springs changes suddenly and the hysteresis effects are reflected smoothly 15,17,28 . The profiles of the clutch torque f n θ f ,θ f (or T C ) are listed in Table 1 17,28 . Figure 1b illustrates the 1st to 4th stages of the stiffness and hysteresis areas. To investigate the dynamic characteristics, the parameters employed for the torsional system shown in Fig. 1a are as follows 7 : inertia of flywheel, I f = 1.38 × 10 −1 kg·m 2 ; viscous damping, c f = 1.59 N m s/rad. Based on the single-degree-of-freedom system shown in Fig. 1a, the equation of motion with the nonlinear function f n θ f ,θ f is derived as follows: Here, T E (t) and T D are the sinusoidal input and drag torques, respectively. In general, the input torque is given by the 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 denote the excitation frequency and phase angle, respectively; and N max represents 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.
T pi cos iω p t + ϕ pi . Its corresponding formulae are defined below.
Likewise, its nonlinear and input functions are defined as follows.
To calculate the solutions for θ c in Eq. (10) with their relevant parameter ̟ , the Newton-Raphson method is used based on the condition → 0 , where is considered as a function of θ c and ̟ , i.e., � θ c , ̟ . Using the Newton-Raphson method, the solutions for θ c and ̟ for each step were determined, as presented in previous publications 16,17 . Examination of initial results according to frequency upsweeping. Figure 2 shows the initial results obtained by the HBM with η = 2 and N max = 24 under the frequency upsweeping condition. In order to capture the nonlinear responses induced by the asymmetrical piecewise type nonlinearities, sufficient number such as N max = 24 is employed based on the prior studies 16,17,28 . Also, this study focuses on frequency upsweeping conditions to examine the nonlinear dynamic characteristics. To examine the stability conditions, Hill's method was employed in this study 16,17,27 . Here, the red dotted lines marked with (b) and (c) in Fig. 2a indicate the super-and sub-harmonic response regimes. To obtain the sub-harmonic responses, a small range of excitation input values is employed artificially 10 . For this study, the valid components of the input torque vector with η = 2 and N max = 24 are expressed as F m = 168.9 , F a(2) = −87.97 and F b(2) = 235.65 in Eq. (9). To trigger 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 expressed as F a(1) = εF a(2) and F b(1) = εF b(2) , respectively. Here, the value employed for ε is 1 × 10 −5 . Figure 2b,c) show the system responses in the super-and sub-harmonic regimes, respectively, where the red circles and blue cross lines indicate stable and unstable conditions, respectively. In these regimes, the stability conditions clearly show variation, such as from stable to unstable (STU) and from unstable to stable (UTS) along the arc length under the frequency upsweeping condition. For example, the locations marked with (1) in Fig. 2b,c illustrate that the direction of the solutions changes abruptly from upsweeping to reverse (or from reverse to upsweeping) with the STU (or the UTS). However, the location marked with (2) indicates that the direction of the solutions follows the same trend as the frequency upsweeping with the STU (or the UTS). Thus, as suggested in the specific objectives with regard to the direction of arc length and their relevant stability sin(kψ 0 ) · · · 1 · · · cos(kψ 1 ) sin(kψ 1 ) · · · . . . . . . Examination of super-harmonic responses.   Figure 4a clearly indicates the superharmonic response, which is also observed in the FFT results, as shown in Fig. 5a, because the nonlinear response is located in the stable regimes, as shown in Fig. 3. As the frequency range changes to that of regimes (B), (C), and (D), their dynamic behaviors become complex, as evident from the time histories shown in Fig. 4. When the responses of regimes (B), (C), and (D) are examined carefully, it is observed that the dynamic characteristics are closely related to the stability conditions. For example, when the time histories of regimes (C) and (B) (or (D)) are compared with each other, the dynamic behavior in regime (C) clearly includes the number of harmonic components dependent on the fundamental frequency rather than the ones in regime (B) (or (D)). When the dynamic response of regime (C) is considered with respect to the stability conditions, it reflects the stable condition. However, the dynamic responses in regimes (B) and (D) were estimated as unstable conditions.
In addition, the FFT results in Fig. 5 indicate the relevance of the dynamic characteristics corresponding to the stability conditions. For instance, Fig. 5c shows that the system responses consist of harmonic terms relevant to their fundamental harmonic components such as ω = 5.2 Hz (or ̟ = 0.285). However, the system responses in the unstable regimes (B) and (D) reflect various numbers of harmonic components, as shown in Fig. 5b,d. The dynamic behavior in regime (D), shown in Fig. 5d, exhibits the presence of the most complex components.
The dynamic behaviors under various stability conditions with respect to the STU and UTS can be thoroughly examined when their responses are projected in phase diagrams and Poincare maps, as shown in Figs. 6 and 7. To determine the phase diagrams and Poincare maps using the NS, this study employed the modified Runge-Kutta method 29 . The system responses were estimated for 500 cycles until the transient response effects were completely removed, and the last 200 cycles were considered to obtain the results. The system response in regime (A), shown show more complex circles. When the dynamic behavior in the stable regime (C) is examined, it is evident that the phase diagram is also composed of a finite number of circles, and the number of dots in the Poincare map is relatively low, e.g., two dots, as shown in Fig. 7c. However, the dynamic responses in the unstable regimes (B) and (D) show a greater number of circles, even though the overall shape of the circle is evident from Fig. 6b, d. These dynamic characteristics are projected well in the Poincare maps, as shown in Fig. 7b, d. In particular, the dynamic behavior in regime (D) shows a chaotic phenomenon, which is reflected by various numbers of dotted lines, as shown in Fig. 7d. Based on the examination of the dynamic responses along different frequency regimes by considering the stability conditions, their dynamic behaviors can be understood more clearly when they are projected in bifurcation diagrams. To create the bifurcation diagrams, 200 cycles of RMS values were obtained from each cycle.
Overall, when the dynamic behaviors in regimes (A), (B), (C), and (D), shown in Figs. 4, 5, 6 and 7, are compared with the bifurcation diagram shown in Fig. 8, the response in each regime can be determined in a clear way. For example, regime (A) does not show any bifurcation, which has already been confirmed as a pure super-harmonic response. As the response transitions into regimes (B), (C), and (D), the relevant bifurcation of these regions is clearly observed with complex characteristics. When the stability conditions are compared with each bifurcation status, regime (B) corresponding to the unstable condition shows a period-doubling cascade immediately after the period-doubling located between regimes (A) and (B). Regime (C) returns to the period-doubling status with the stable condition. When the regime is changed to (D) under the STU condition, the bifurcation becomes more complicated and reflects the chaotic response. Table 3 summarizes the dynamic characteristics under various conditions based on the stability and arc-length direction. In the table, the symbols "+" and "−" indicate the same and reverse directions of frequency upsweeping.  A) and (B), where the direction of arc length is reversed. In Fig. 9, regime (A) shows the fundamental harmonic response area immediately before the sub-harmonic responses occur.

Results of sub-harmonic responses
These are clearly observed in Figs. 10 and 11, respectively. Among the sub-harmonic response areas, i.e., regimes (B), (C), and (D), regimes (C) and (D) pertain to the unstable condition. The arc length in the unstable regimes (C) and (D) follows the same direction as the frequency upsweeping. As the frequency regime transitions from (A) to (B), the dynamic behaviors clearly reflect the sub-harmonic effects in regime (B). When the time histories and FFT results shown in Figs. 10a,b and 11a,b are compared, the variation in the dynamic characteristics In particular, the dynamic response in regime (D) includes various harmonic components, as shown in Fig. 11d. The dynamic behaviors in the sub-harmonic area with respect to STU and UTS were also analyzed using phase diagrams and Poincare maps, as shown in Figs. 12 and 13. From the phase diagrams, the number of circles increases as the frequency range transitions from regime (A) to regime (D). In particular, regime (D) shows a large number of periodic circles, where chaotic motions are expected. These variational dynamic conditions were examined based on the Poincare maps, as shown in Fig. 13. As the regimes transition from (A) to (D), the number of periodic points increases. Figure 13d shows the chaotic motions clearly as variational periodic points are passed through.
Overall, when the dynamic behaviors in regimes (A), (B), (C), and (D) in Figs. 10, 11, 12 and 13 are compared in a bifurcation diagram, as shown in Fig. 14, the strong relationship between the dynamic behaviors and stability conditions is evident, which was examined previously for the super-harmonic area. For example, regime (A) only shows the fundamental harmonic response. However, bifurcation is revealed when the frequency  Table 3.

Conclusion
This study analyzed the nonlinear dynamic characteristics of super-and sub-harmonic response areas. To understand the nonlinear motions related to the unstable regimes, this study suggests the use of stability and bifurcation analysis. First, the dynamic characteristics in the super-and sub-harmonic areas according to the variations in stability conditions were investigated using the HBM. Here, the stability variations along the direction of the arc length were proved to have a close relationship with nonlinear dynamic behaviors such as period-doubling, period-doubling cascade, and chaotic motions. Second, this study analyzed the nonlinear dynamic responses by comparing the bifurcations with the stability conditions. To understand the variation of dynamic behaviors, FFT results, phase diagrams, and Poincare maps were also compared. Furthermore, to investigate and understand the nonlinear dynamic behaviors under the STU and UTS conditions by projecting them in bifurcation diagrams, this study determined the relevant arc-length direction empirically based on the NS results. However, the method to verify the dynamic characteristics with respect to the arclength direction and stability conditions can be defined analytically, which will be the future scope of this work.