Dynamic Instability Investigation of the Automotive Driveshaft’s Forced Torsional Vibration Using the Asymptotic Method

Featured Application

The first-order asymptotic method (FOAM) is used to analyze the dynamic instability in the area of principal parametric resonance (APPR) for the automotive driveshaft’s (AD) forced torsional vibration (ADFTV).

Abstract

This paper aims to investigate using FOAM to analyze the dynamic instability in the APPR for ADFTV based on a dynamic model (DMADFTV). The DMADFTV considers the following aspects: AD kinematic nonuniformity (ADKN), AD geometric nonuniformity (ADGN) of inertial characteristics for the spinning movements (ICSM) of the AD elements (ADE), and the excitations induced by the gearbox–internal combustion engine modulations. The DMADFTV is considered the already-designed dynamic model developed by the first author of the ADFTV in a previous publication. This DMADFTV was used to compute the stationary frontiers of instability and the nonstationary spectral velocity amplitude (NSVA) versus nonstationary spectral amplitude (NSA) in the configuration space in transition through APPR, using the FOAM. The use of FOAM is much more versatile, from the analytical point of view, than the method of multiple scales and allows the computation of the NSA and the NSVA in the APPR. In contrast, these computations cannot be performed using the harmonic balance method. MATLAB Software R 2017 was developed based on DMADFTV and used the FOAM to compute the stationary frontiers of instability and the NSVA versus the NSA in transition through APPR for the ADFTV. The numerical results were compared with the experimental and numerical data published in the literature, finding agreements. The computation of the NSVA versus NSA in the configuration space using FOAM represents a method of detection of the chaotic manifestation of ADFTV.

1. Introduction

The present work investigates the ADFTV, taking into consideration the following physical characteristics and phenomena: the fact that the links of the ADE do not manifest the property of isometry due to the ADKN and ADGN [1,2], and the ADE’s links are even considered CVJs (constant velocity joints); the effect of the torque loads, such as the entry gearbox torques [3] (p. 360); and the effects of ADGN of the ICSM, and the effects of the elasticity and damping in torsion for the ADE’s link-joints, which are functions of the ADE twisting angle. The ADKN and ADGN that generate the isometry nonuniformity properties of AD (INP-AD) are mentioned in [2,4] and experimentally studied by Steinwede in his Ph.D. thesis [5] (pp. 77–94). This phenomenon of INP-AD induces a nonlinear behavior of driveshafts generating forced parametric vibrations in the frequency range of 0.1–12 kHz, as shown by published experiments and data in [5]. Mazzei and Scott were the first to study the dynamic behavior of driveshafts, and they discovered through experiments the nonlinear parametric vibrations (NPV) of a universal AD joint [6]. Their research revealed and illustrated experiments concerning the INP-AD of quasi-CVJ driveshaft transmission [7] by Browne and Palazzolo. Xu et al. studied the amplitude variations as functions of frequency for the ADFTV considering the complete powertrain [8], developing a dynamic model with 29 DOF (degrees of freedom). Idehara et al. developed a model for the nonlinear ADFTV 3 DOF [9]. The rigidity and the damping of ADFTV were determined for calibration. Unfortunately, the DMADFTV developed in [9] does not have a nonlinear characteristic involving ICSM, ADKN, and ADGN for the ADE, respectively. Feng et al. studied and analyzed in [10] the dynamics of a driveshaft system with an interval of uncertainty, while Tiberiu-Petrescu and Petrescu realized a survey regarding the theoretical modelization of the geometrical structure and its influence on the kinematics of a universal joint [11]. In [12], Ertürka et al. developed an analysis of the NVH for a truck driveshaft using an SSA (Six Sigma Approach), while in his master thesis, Kamalakkannan modeled and simulated the kinematics and dynamics of the drive train [13]. The design and analysis of AD were investigated by Kishore et al. [14], and Shao et al. developed a model control for vehicle rollover prevention with time-varying speed [15]. Previously, Deng, Zhao, et al. [16] studied and implemented a hierarchical synchronization control for the active rear-axle steering system. The investigation of high-frequency torsional vibration of automotive drive lines (ADL) was done by Farashindianfar et al. using genetic algorithms [17]. At the same time, Komorska et al. researched diagnosing defects for ADL by analyzing induced vibration signals generated by such defect phenomena [18]. In [19], Alugongo studied the dynamic parametric behavior of the Cardan shaft, also performing a sensitivity analysis.

This study aimed to investigate the dynamic instability of the ADFTV for the ADE by developing a detection method for chaotic manifestation. The detection method represents the construction of a phase portrait of the NSVA as a function of the NSA in transition through APPR for the ADFTV. Such a picture’s phase portrait analysis shall be conducted to confirm the chaotic manifestation of the ADFTV, such as Jacobi stability analysis or Lyapunov stability analysis.

Many researchers have used the Jacobi stability analysis of dynamical systems in the last two decades based on the Kosambi–Cartan–Chern (KCC) theory. Böhmer et al. studied the Jacobi stability analysis in gravitation and cosmology using the KCC formalism in Riemannian and Finslerian spaces [20], while Harko et al. using the same KCC formalism determined the geometric invariants and the Berwald’s connection for the Lorenz system to make the Jacobi stability analysis [21]. In recent years, research concerning Jacobi stability analysis was carried out by Wang et al. concerning the onset of chaos in a two-degree-of-freedom mechanical system [22] and by Blaga et al. concerning the classical restricted three-body problem [23]. Mallikarjun and Narasimhamurthy investigated in [24] the metrics of Finsler Space for nonholonomic Finsler frame, such as Matsumoto metric and infinite series of $\left(\alpha ,\mathsf{\beta }\right)$ metric with applications to field theory. Roy et al. studied the chaos of hidden strange attractors in nonlinear autonomous systems based on KCC by determining the geometric invariants [25].

The Lyapunov stability analysis for mechanical systems was studied and developed in [26] (pp. 98–127, pp. 298–308) based on the method of Lyapunov exponents and the Kaplan–Yorke conjecture stating “that a system’s complete Lyapunov spectrum may be used to determine the fractal dimension” [26] (p. 301).

An AD is a homokinetic mechanical system designed to transmit a spinning moment from the engine crankshaft to the automotive wheel, as shown in Figure 1. The present paper considers the morphology of an AD in use on the models of TOYOTA’s heavy-duty off-road SUV (Sport Utility Vehicle); such morphology is illustrated in Figure 2 and Figure 3.

Looking inside the ADE, this kind of transmission, illustrated in Figure 2, has as its basic components (a) the AD’s global bowl, (b) the AD’s axis, and (c) the AD’s tulip. The TTJ’s (tulip–tripod joint-link) components are presented in Figure 3A, having the following parts:

• the AD–MA (automotive’s driveshaft-middle ax),
• the tripod linked with tulip and assembled using splines on the middle ax,
• the tulip’s bell (TB),
• the tulip ax (TA), assembled with splines in the geared box.

Figure 3B details the BIRJ’s (bowl–inner race joint-link) components. It shows the following components:

• the AD–MA,
• the inner race (IR) assembled with splines on the AD–MA,
• the BB (bowl’s bell) linked with IR through balls,
• the BA (bowl’s ax) linked to the SW (steering wheel) with splines.

2. The Analysis of DMADFTV

Figure 4 shows a simplified picture of the AD, where each AD’s element has attached a Cartesian reference system such as

−

$X_1{Y}_1{Z}_1$ attached to the GT (global tulip: TB and TA), having a rigid rotation with ${\mathsf{\phi }}_1$ around axis $X_1$,

−

$X_2{Y}_2{Z}_2$ attached to the MA, having a rigid rotation with ${\mathsf{\phi }}_2$ around axis $X_2$,

−

$X_3{Y}_3{Z}_3$ attached to the GB (global bowl: BB and BA), having a rigid rotation with ${\mathsf{\phi }}_3$ around axis $X_3$, while ${\mathsf{\beta }}_1$ is the angle between the longitudinal directions of the GT and the MA regarding the axis $Z_1$, ${\mathsf{\beta }}_1\in \left[0°,15°\right]$, and ${\mathsf{\beta }}_2$ is the angle between the longitudinal directions of the GB and the MA regarding the axis $Z_3$, ${\mathsf{\beta }}_2\in \left[0°,47°\right]$. The angle ${\mathsf{\beta }}_1$ between the longitudinal directions of the GT and the MA regarding the axis $Z_1$, is given by the internal length of the trucks inside the TB, the external diameter of the MA, and the internal diameter of the TB (see Figure 5a), conforming to the designer calculus [27] (pp. 214–230). The angle ${\mathsf{\beta }}_2$ between the longitudinal directions of the GB and the MA regarding the axis $Z_3$ is given by the internal length of the trucks inside the IR, the external diameter of the MA, and the internal diameter of the BB (see Figure 5b), conforming to the designer calculus [27] (pp. 153–208).

An AD is a multibody driveline transmission from the automotive geared box to the wheel with the geometrical details presented in Figure 5.

To determine the differential equations of the ADFTV based on the variational approach of Hamilton’s principle (VAHP), it is mandatory to realize the reduction of the AIMM (axial inertial mass moment) and the AIGM (axial inertial geometric moment) for the GT and the GB to the longitudinal axis of the MA, respectively $X_2$, as illustrated in Figure 5. The AIMM and the AIGM of the GT and the GB consider the angles ${\mathsf{\beta }}_1$ and ${\mathsf{\beta }}_2$ as well as the specific geometry dimensions of each GT’s components: TB, TA reduced in the tripod’s CM (center mass); and GB’s components: BB, BA reduced in the inner race’s CM, as expressed by the Equations (A1)–(A4) in Appendix A. In Appendix A, Equation (A5) expressed the AIMM and the AIGM for the AD’s MA, while Equation (A6) expressed the stiffnesses and the dampings for the GT and the GB.

Figure 6 illustrates the DMADFTV (dynamic model of the ADFTV).

The DMADFTV comprises the following elements: the GT, the MA, and the GB, the GT being linked in torsion with the MA through TTJ–CVJ while the GB is linked in torsion with the MA through BIRJ–CVJ. The geared box injects in the AD a moment produced by the internal combustion engine, yielding [3] (p. 361)

$$M_{Gb}={\overline{M}}_e\left[1+{\chi }_e\mathit{cos}\left(n{\Omega }_{e}t\right)\right],n\in N,{\Omega }_e={\displaystyle \frac{\pi n_e}{30}},$$

(1)

where ${\chi }_e$ is the engine moment nonuniformity, having values 0.10–0.12 [3] (p. 363), ${\overline{M}}_e$ is the moment’s amplitude, and $n_e$ is the engine crankshaft rpm. The moment induced at the wheel is a moderately impulsive kind and has the expression

$$M_w=M_H\left[1+q_3{t}^{q_1}{e}^{-q_{2}t}\right],$$

(2)

where $M_H$ is the moment of adhesion [27] (p. 130), where $q_i,i=\overline{1,3},q_1\gg q_2,q_3\ge 1.1$ are experimental constants depending on the shock applied to the wheel by the excitation [28].

3. The Equations of ADFTV and the Mathematical Solutions of FOAM

For the DMADFTV illustrated in Figure 6, based on the VAHP [29] (pp. 272–295) results,

$$\mathsf{\delta }{\displaystyle \underset{{\mathrm{t}}_1}{\overset{{\mathrm{t}}_2}{\mathsf{\int }}}\left(\mathsf{\Pi }-\mathrm{T}-\mathrm{W}\right)\mathrm{d}\mathrm{t}=0,}$$

(3)

where $\mathsf{\Pi }$ is strain energy of the DMADFTV that includes the springs and the dampers, as presented in Appendix A [29] (p. 274), [30] (pp. 610–613); $T$ is the kinetic energy of the DMADFTV [29] (p. 274), [30] (p. 719); and W is the work performed by the external torques. In Appendix B, the significance of the terms is shown in [1]. It becomes evident that the MA is a fixed–fixed uniform beam linked in spinning to the TTJ at $x=0$ and at the BIRJ at $x=L_{Ms}$ and, therefore, ${\mathsf{\phi }}_2\left(\mathrm{x},\mathrm{t}\right)$ has the expression presented in Appendix B (see ([30] (p. 720)), where ${\mathsf{\phi }}_{2n}$ is the amplitude, ${\mathsf{\Theta }}_{1n}$ and ${\mathsf{\Theta }}_{2n}$ are the phase angles for the general solution ${\mathsf{\phi }}_2(x,t)$ and ${\omega }_n$ is the free natural frequency of torsion vibrations of the midshaft. After mathematical calculus yielding the ADFTV’s equations,

$$\begin{array}{l}{\ddot{\mathsf{\phi }}}_1+2\mathsf{\zeta }{\overline{\mathsf{\Omega }}}_1\sqrt{{\displaystyle \frac{1+{\mathrm{a}}_1\cos 2{\mathsf{\phi }}_1}{\left(1+{\mathrm{a}}_2\cos 2{\mathsf{\phi }}_1\right)\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}\cos 2{\mathsf{\phi }}_1\right)}}}{\dot{\mathsf{\phi }}}_1+{\overline{\mathsf{\Omega }}}_1^2{\displaystyle \frac{1+{\mathrm{a}}_1\cos 2{\mathsf{\phi }}_1}{\left(1+{\mathrm{a}}_2\cos 2{\mathsf{\phi }}_1\right)\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}\cos 2{\mathsf{\phi }}_1\right)}}{\mathsf{\phi }}_1=\\ ={\displaystyle \frac{{\overline{\mathrm{M}}}_{\mathrm{e}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}}{\displaystyle \frac{1+{\mathrm{\chi }}_{\mathrm{e}}\cos (\mathrm{n}{\mathsf{\Omega }}_{\mathrm{e}}\mathrm{t})}{\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}\cos 2{\mathsf{\phi }}_1\right)}}\left(1-3{\mathrm{A}}_{\mathrm{T}\mathrm{T}\mathrm{r}}\sin 3{\mathsf{\phi }}_1\right)+\\ +{\displaystyle \sum _{\mathrm{n}}{\mathrm{\Phi }}_{2\mathrm{n}}\left[\frac{{\mathrm{I}}_{01\mathrm{M}\mathrm{s}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}{\mathrm{\omega }}_{\mathrm{n}}^2\cos {\mathsf{\Theta }}_{1\mathrm{n}}+\frac{\mathrm{G}{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{M}\mathrm{s}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}\frac{{\mathrm{\omega }}_{\mathrm{n}}}{\mathrm{c}}\sin {\mathsf{\Theta }}_{1\mathrm{n}}\right]\frac{\cos \left({\mathrm{\omega }}_{\mathrm{n}}\mathrm{t}-{\mathsf{\Theta }}_{2\mathrm{n}}\right)}{\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}\cos 2{\mathsf{\phi }}_1\right)}},\mathrm{n}=1,2,3,\dots ,\end{array}$$

(4)

$$\begin{array}{l}{\ddot{\mathsf{\phi }}}_3+2\mathsf{\zeta }{\overline{\mathsf{\Omega }}}_3\sqrt{{\displaystyle \frac{1+{\mathrm{a}}_3\cos 2{\mathsf{\phi }}_3}{\left(1+{\mathrm{a}}_4\cos 2{\mathsf{\phi }}_3\right)\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}\cos 2{\mathsf{\phi }}_3\right)}}}{\dot{\mathsf{\phi }}}_3+{\overline{\mathsf{\Omega }}}_3^2{\displaystyle \frac{1+{\mathrm{a}}_3\cos 2{\mathsf{\phi }}_3}{\left(1+{\mathrm{a}}_4\cos 2{\mathsf{\phi }}_3\right)\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}\cos 2{\mathsf{\phi }}_3\right)}}{\mathsf{\phi }}_3=\\ =-{\displaystyle \frac{{\mathrm{M}}_{\mathrm{H}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}}{\displaystyle \frac{\left[1+{\mathrm{q}}_3{\mathrm{t}}^{{\mathrm{q}}_1}{\mathrm{e}}^{-{\mathrm{q}}_2\mathrm{t}}\right]}{\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}\cos 2{\mathsf{\phi }}_3\right)}}\left(1-3{\mathrm{A}}_{\mathrm{B}\mathrm{I}\mathrm{r}}\sin 3{\mathsf{\phi }}_3\right)+\\ +{\displaystyle \sum _{\mathrm{n}}{\left(-1\right)}^{\mathrm{n}}{\mathrm{\Phi }}_{2\mathrm{n}}\left[\frac{{\mathrm{I}}_{02\mathrm{M}\mathrm{s}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}{\mathrm{\omega }}_{\mathrm{n}}^2\cos {\mathsf{\Theta }}_{1\mathrm{n}}+\frac{\mathrm{G}{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{M}\mathrm{s}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}\frac{{\mathrm{\omega }}_{\mathrm{n}}}{\mathrm{c}}\sin {\mathsf{\Theta }}_{1\mathrm{n}}\right]\frac{\cos \left({\mathrm{\omega }}_{\mathrm{n}}\mathrm{t}-{\mathsf{\Theta }}_{2\mathrm{n}}\right)}{\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}\cos 2{\mathsf{\phi }}_3\right)}},\mathrm{n}=1,2,3,\dots ,\end{array}$$

(5)

where ${\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}$ is the GT’s AIMM concerning the axis X₂ for ${\mathrm{\beta }}_1=0°$, ${\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}$ is the GT nonuniformity (as in Appendix B), ${\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}$ is the GB’s AIMM concerning the axis X₂ for ${\mathrm{\beta }}_2=0°$, ${\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}$ is the GB nonuniformity (as in Appendix B), and ${\overline{\mathsf{\Omega }}}_1$ and ${\overline{\mathsf{\Omega }}}_3$ are the circular eigenfrequencies in torsion of the GT and GB as functions of ${\mathsf{\beta }}_1,{\mathsf{\beta }}_2$, having the mathematical expressions

$$\begin{gathered} {\overline{\mathsf{\Omega }}}_1={\displaystyle \frac{\mathrm{c}}{{\mathrm{L}}_{\mathrm{T}}}}{\left[{\displaystyle \frac{1-\frac{{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{A}\mathrm{T}}\mathrm{\rho }{\mathrm{L}}_{\mathrm{A}\mathrm{T}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}}{1+{\left(\frac{{\mathrm{L}}_{\mathrm{A}\mathrm{T}}}{{\mathrm{L}}_{\mathrm{T}}}\right)}^2\left(\frac{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}{{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{A}\mathrm{T}}\mathrm{\rho }{\mathrm{L}}_{\mathrm{A}\mathrm{T}}}-1\right)}}\right]}^{1/2}, \\ {\overline{\mathsf{\Omega }}}_3={\displaystyle \frac{\mathrm{c}}{{\mathrm{L}}_{\mathrm{B}}}}{\left[{\displaystyle \frac{1-\frac{{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{A}\mathrm{B}}\mathrm{\rho }{\mathrm{L}}_{\mathrm{A}\mathrm{B}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}}{1+{\left(\frac{{\mathrm{L}}_{\mathrm{A}\mathrm{B}}}{{\mathrm{L}}_{\mathrm{B}}}\right)}^2\left(\frac{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}{{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{A}\mathrm{B}}\mathrm{\rho }{\mathrm{L}}_{\mathrm{A}\mathrm{B}}}-1\right)}}\right]}^{1/2}. \end{gathered}$$

(6)

For the dynamic instability analysis for ADFTV, different approaches can be used, such as the harmonic balance method (HBM), asymptotic method (AM) [31] (pp. 299–393), and the method of multiple scales (MMS) [32] (pp. 424–427). The HBM is very efficient but allows only the computation of the stationary amplitude of vibrations; the AM permits the calculation of the amplitude and phase angle for both stationary and nonstationary oscillations, while MMS raises computational difficulties due to zeroing the secular term and supplementary analysis for the detuning parameter [31] (pp. 424–427). For the above reasons, the present work uses the asymptotic method approach (AMA) for the solutions of Equations (4) and (5). The use of first-order approximation of AM makes the investigation of dynamic instability of ADFTV in both stationary and nonstationary torsional vibrations in the APPR possible. A small positive parameter must be introduced for both stationary and nonstationary torsional vibrations; the slowing time $\mathsf{\tau }=\mathsf{\epsilon }\mathrm{t}$, ε being a positive small number, as mentioned in [30] (p. 299). These equations must be modified to apply AM before raising $\tau$ in Equations (4) and (5). For the AM’s first-order approximation (FOAM), the last terms in Equations (4) and (5) are neglected because the excitation frequency in APPR for GT and the GB is not in the range of free natural frequencies for torsional vibrations of the midshaft. The development gives the approximations used in the mathematical procedure of FOAM in Taylor’s series of the following functions:

$$\begin{array}{l}{\displaystyle \frac{1+{\mathrm{a}}_1\cos 2{\mathsf{\phi }}_1}{\left(1+{\mathrm{a}}_2\cos 2{\mathsf{\phi }}_1\right)\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}\cos 2{\mathsf{\phi }}_1\right)}}\simeq 1+2{\mathrm{\mu }}_1\cos 2{\mathsf{\phi }}_1\Rightarrow {\mathrm{\mu }}_1\simeq \left(1+{\mathrm{\alpha }}_1\right){\mathrm{\alpha }}_2,\\ {\displaystyle \frac{{\overline{\mathrm{M}}}_{\mathrm{e}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}}{\displaystyle \frac{1+{\mathrm{\chi }}_{\mathrm{e}}\cos (\mathrm{n}{\mathsf{\Omega }}_{\mathrm{e}}\mathrm{t})}{\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}\cos 2{\mathsf{\phi }}_1\right)}}\left(1-3{\mathrm{A}}_{\mathrm{T}\mathrm{T}\mathrm{r}}\sin 3{\mathsf{\phi }}_1\right)\simeq {\mathrm{\Gamma }}_1{\mathsf{\phi }}_1^3\Rightarrow {\mathrm{\Gamma }}_1=9.1125{\displaystyle \frac{{\mathrm{A}}_{\mathrm{T}\mathrm{T}\mathrm{r}}}{{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}}}{\displaystyle \frac{{\overline{\mathrm{M}}}_{\mathrm{e}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}},\end{array}$$

(7)

$$\begin{array}{l}{\displaystyle \frac{1+{\mathrm{a}}_3\cos 2{\mathsf{\phi }}_3}{\left(1+{\mathrm{a}}_4\cos 2{\mathsf{\phi }}_3\right)\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}\cos 2{\mathsf{\phi }}_3\right)}}\simeq 1+2{\mathrm{\mu }}_3\cos 2{\mathsf{\phi }}_3\Rightarrow {\mathrm{\mu }}_3\simeq \left(1+{\mathrm{\alpha }}_3\right){\mathrm{\alpha }}_4,\\ -{\displaystyle \frac{{\mathrm{M}}_{\mathrm{H}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}}{\displaystyle \frac{\left[1+{\mathrm{q}}_3{\mathrm{t}}^{{\mathrm{q}}_1}{\mathrm{e}}^{-{\mathrm{q}}_2\mathrm{t}}\right]}{\left(1+{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}\cos 2{\mathsf{\phi }}_3\right)}}\left(1-3{\mathrm{A}}_{\mathrm{B}\mathrm{I}\mathrm{r}}\sin 3{\mathsf{\phi }}_3\right)={\mathrm{\Gamma }}_3{\mathsf{\phi }}_3^3+{\mathrm{\Gamma }}_5{\mathsf{\phi }}_3^5,\\ \Rightarrow {\mathrm{\Gamma }}_3=-9.1125{\displaystyle \frac{{\mathrm{A}}_{\mathrm{B}\mathrm{I}\mathrm{r}}}{{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}}}{\displaystyle \frac{{\mathrm{M}}_{\mathrm{H}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}},{\mathrm{\Gamma }}_5=-18.225{\displaystyle \frac{{\mathrm{A}}_{\mathrm{B}\mathrm{I}\mathrm{r}}}{{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}}}{\displaystyle \frac{{\mathrm{M}}_{\mathrm{H}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}},\end{array}$$

(8)

whose coefficients ${\mathrm{\alpha }}_{\mathrm{i}},\mathrm{i}=1,2,3,4$ are presented in Appendix B. In Equations (7) and (8), the terms ${\mathrm{\mu }}_1={\mathrm{f}}_1\left({\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}},{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{A}\mathrm{T}},\mathrm{\rho },{\mathrm{L}}_{\mathrm{A}\mathrm{T}},{\left({\displaystyle \frac{{\mathrm{L}}_{\mathrm{A}\mathrm{T}}}{{\mathrm{L}}_{\mathrm{T}}}}\right)}^2,{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}\right)$, ${\mathrm{\mu }}_3={\mathrm{f}}_3\left({\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}},{\mathrm{J}}_{{\mathrm{X}}_2\mathrm{A}\mathrm{B}},\mathrm{\rho },{\mathrm{L}}_{\mathrm{A}\mathrm{B}},{\left({\displaystyle \frac{{\mathrm{L}}_{\mathrm{A}\mathrm{B}}}{{\mathrm{L}}_{\mathrm{B}}}}\right)}^2,{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}\right)$ express the internal parametric excitations of GT and GB (see Equations (A13)–(A16) in Appendix B), while the terms ${\mathrm{\Gamma }}_1={\mathrm{g}}_1\left({\displaystyle \frac{{\mathrm{A}}_{\mathrm{T}\mathrm{T}\mathrm{r}}}{{\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}}},{\displaystyle \frac{{\overline{\mathrm{M}}}_{\mathrm{e}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{T}}}}\right)$, ${\mathrm{\Gamma }}_3={\mathrm{g}}_3\left({\displaystyle \frac{{\mathrm{A}}_{\mathrm{B}\mathrm{I}\mathrm{r}}}{{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}}},{\displaystyle \frac{{\mathrm{M}}_{\mathrm{H}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}}\right)$, ${\mathrm{\Gamma }}_5={\mathrm{g}}_5\left({\displaystyle \frac{{\mathrm{A}}_{\mathrm{B}\mathrm{I}\mathrm{r}}}{{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}}},{\displaystyle \frac{{\mathrm{M}}_{\mathrm{H}}}{{\overline{\mathrm{I}}}_{{\mathrm{X}}_2\mathrm{G}\mathrm{B}}}}\right)$ express the external cubic and quintic external excitations for GT and GB (see Equations (2) and (A10) in Appendix B). The functions ${\mathrm{f}}_{\mathrm{i}},\mathrm{i}=1,2$ and ${\mathrm{g}}_{\mathrm{i}},\mathrm{i}=1,2,3$ have the variables mentioned above represent the physical sources of these internal and external excitations. Injecting the relations (7) in Equation (4) and repeating the same manipulation by injecting the relations (8) in Equation (5) yields

$${\ddot{\mathsf{\phi }}}_1++{\overline{\mathsf{\Omega }}}_1^2{\mathsf{\phi }}_1=-2\mathsf{\zeta }{\overline{\mathsf{\Omega }}}_1(1+{\mathrm{\alpha }}_1+{\mathrm{\alpha }}_2\cos 2{\mathsf{\phi }}_1){\dot{\mathsf{\phi }}}_1-{\overline{\mathsf{\Omega }}}_1^2(2{\mathrm{\alpha }}_1+{\mathrm{\alpha }}_1^2+2{\mathrm{\mu }}_1\cos 2{\mathsf{\phi }}_1){\mathsf{\phi }}_1+{\mathrm{\Gamma }}_1{\mathsf{\phi }}_1^3,$$

(9)

$${\ddot{\mathsf{\phi }}}_3++{\overline{\mathsf{\Omega }}}_3^2{\mathsf{\phi }}_3=-2\mathsf{\zeta }{\overline{\mathsf{\Omega }}}_3(1+{\mathrm{\alpha }}_3+{\mathrm{\alpha }}_4\cos 2{\mathsf{\phi }}_1){\dot{\mathsf{\phi }}}_3-{\overline{\mathsf{\Omega }}}_3^2(2{\mathrm{\alpha }}_3+{\mathrm{\alpha }}_3^2+2{\mathrm{\mu }}_3\cos 2{\mathsf{\phi }}_3){\mathsf{\phi }}_3+{\mathrm{\Gamma }}_3{\mathsf{\phi }}_3^3+{\mathrm{\Gamma }}_5{\mathsf{\phi }}_3^5.$$

(10)

Equations (9) and (10) represent the FTV of the GT and GB in APPR that include the following effects: the kinematic and geometric nonuniformity of quasi-homokinetic tripod’s joints, the nonuniformity of AIGM and AIMM for the elements of AD and the nonuniformity of external harmonic excitations of the engine, and the reaction torque induced by the wheel. The excitation frequencies ${\mathsf{\eta }}_1,{\mathsf{\eta }}_3$ of the GT and the GB of AD in the APPR must be defined by the expressions

$${\mathsf{\eta }}_3={\displaystyle \frac{\mathrm{d}{\mathsf{\Theta }}_3}{\mathrm{d}\mathrm{t}}}={\displaystyle \frac{2\mathrm{d}{\mathsf{\phi }}_3}{\mathrm{d}\mathrm{t}}}\simeq 2{\overline{\mathsf{\Omega }}}_3,{\mathsf{\eta }}_1={\displaystyle \frac{\mathrm{d}{\mathsf{\Theta }}_1}{\mathrm{d}\mathrm{t}}}={\displaystyle \frac{2\mathrm{d}{\mathsf{\phi }}_1}{\mathrm{d}\mathrm{t}}}\simeq 2{\overline{\mathsf{\Omega }}}_1,$$

(11)

and introducing the slowing time by imposing the condition that all the coefficients in Equations (9) and (10) are small yields

$${\displaystyle \frac{\mathrm{d}{\mathsf{\Theta }}_1}{\mathrm{d}\mathrm{t}}}={\mathsf{\eta }}_1(\mathrm{\tau }),{\mathsf{\Theta }}_1=2{\mathsf{\phi }}_1,\mathrm{\xi }=\mathrm{\epsilon }\mathrm{\xi },{\mathrm{\mu }}_1=\mathrm{\epsilon }{\mathrm{\mu }}_1,{\mathrm{\alpha }}_1=\mathrm{\epsilon }{\mathrm{\alpha }}_1,{\mathrm{\alpha }}_2=\mathrm{\epsilon }{\mathrm{\alpha }}_2,{\mathrm{\Gamma }}_1=\mathrm{\epsilon }{\mathrm{\Gamma }}_1,$$

(12)

$${\displaystyle \frac{\mathrm{d}{\mathsf{\Theta }}_3}{\mathrm{d}\mathrm{t}}}={\mathsf{\eta }}_3(\mathrm{\tau }),{\mathsf{\Theta }}_3=2{\mathsf{\phi }}_3,\mathrm{\xi }=\mathrm{\epsilon }\mathrm{\xi },{\mathrm{\mu }}_3=\mathrm{\epsilon }{\mathrm{\mu }}_3,{\mathrm{\alpha }}_3=\mathrm{\epsilon }{\mathrm{\alpha }}_3,{\mathrm{\alpha }}_4=\mathrm{\epsilon }{\mathrm{\alpha }}_4,{\mathrm{\Gamma }}_3=\mathrm{\epsilon }{\mathrm{\Gamma }}_3,{\mathrm{\Gamma }}_5=\mathrm{\epsilon }{\mathrm{\Gamma }}_5.$$

(13)

Neglecting the terms that contain ${\mathsf{\epsilon }}^2$ after injecting relations (12), (13) in Equations (9) and (10) results in

$${\ddot{\mathsf{\phi }}}_1++{\overline{\mathsf{\Omega }}}_1^2{\mathsf{\phi }}_1=\mathrm{\epsilon }\left[-2\mathsf{\zeta }{\overline{\mathsf{\Omega }}}_1{\dot{\mathsf{\phi }}}_1-{\overline{\mathsf{\Omega }}}_1^2(2{\mathrm{\alpha }}_1+2{\mathrm{\mu }}_1\cos 2{\mathsf{\phi }}_1){\mathsf{\phi }}_1+{\mathrm{\Gamma }}_1{\mathsf{\phi }}_1^3\right],$$

(14)

$${\ddot{\mathsf{\phi }}}_3++{\overline{\mathsf{\Omega }}}_3^2{\mathsf{\phi }}_3=\mathrm{\epsilon }\left[-2\mathsf{\zeta }{\overline{\mathsf{\Omega }}}_3{\dot{\mathsf{\phi }}}_3-{\overline{\mathsf{\Omega }}}_3^2(2{\mathrm{\alpha }}_3+2{\mathrm{\mu }}_3\cos 2{\mathsf{\phi }}_3){\mathsf{\phi }}_3+{\mathrm{\Gamma }}_3{\mathsf{\phi }}_3^3+{\mathrm{\Gamma }}_5{\mathsf{\phi }}_3^5\right].$$

(15)

Analyzing Equations (14) and (15), it can be remarked that the global terms on the right-hand side (RHS) are excitations, while in contrast, the left-hand side (LHS) global parts are linear oscillators. For FOAM, there were meant to be analytical solutions for Equations (14) and (15), having the expressions

$${\mathsf{\phi }}_1={\mathrm{\Phi }}_1\cos \left({\displaystyle \frac{1}{2}}{\mathsf{\Theta }}_1+{\mathrm{\Psi }}_1\right),{\mathsf{\phi }}_3={\mathrm{\Phi }}_3\cos \left({\displaystyle \frac{1}{2}}{\mathsf{\Theta }}_3+{\mathrm{\Psi }}_3\right),$$

(16)

where ${\mathrm{\Phi }}_1,{\mathrm{\Phi }}_3$ are the amplitudes and ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_3$ are the phase angles for global tulip and global bowl, given by the Equations

$$\begin{gathered} \left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_1}{\mathrm{d}\mathrm{t}}}=\mathrm{\epsilon }{\mathrm{D}}_1\left(\mathrm{\tau },{\mathrm{\Phi }}_1,{\mathrm{\Psi }}_1\right)\\ {\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_1}{\mathrm{d}\mathrm{t}}}={\overline{\mathsf{\Omega }}}_1-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_1\left(\mathrm{\tau }\right)+\mathrm{\epsilon }{\mathrm{E}}_1\left(\mathrm{\tau },{\mathrm{\Phi }}_1,{\mathrm{\Psi }}_1\right)\end{array}\right., \\ \left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_3}{\mathrm{d}\mathrm{t}}}=\mathrm{\epsilon }{\mathrm{D}}_3\left(\mathrm{\tau },{\mathrm{\Phi }}_3,{\mathrm{\Psi }}_3\right)\\ {\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_3}{\mathrm{d}\mathrm{t}}}={\overline{\mathsf{\Omega }}}_3-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_3\left(\mathrm{\tau }\right)+\mathrm{\epsilon }{\mathrm{E}}_3\left(\mathrm{\tau },{\mathrm{\Phi }}_3,{\mathrm{\Psi }}_3\right)\end{array}\right.. \end{gathered}$$

(17)

Differentiating concerning the time, the system of Equation (17) yields

$$\begin{gathered} \left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Phi }}_1}{\mathrm{d}{\mathrm{t}}^2}}=\mathrm{\epsilon }{\displaystyle \frac{\partial {\mathrm{D}}_1}{\partial {\mathrm{\Psi }}_1}}\left({\overline{\mathsf{\Omega }}}_1-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_1\right)+{\mathrm{\epsilon }}^2\dots \dots \dots \dots \dots \dots \dots \dots \\ {\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Psi }}_1}{\mathrm{d}{\mathrm{t}}^2}}=\mathrm{\epsilon }\left[-{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial {\mathsf{\eta }}_1}{\partial \mathrm{\tau }}}+{\displaystyle \frac{\partial {\mathrm{E}}_1}{\partial {\mathrm{\Psi }}_1}}\left({\overline{\mathsf{\Omega }}}_1-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_1\right)\right]+{\mathrm{\epsilon }}^2\dots \end{array}\right., \\ \left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Phi }}_3}{\mathrm{d}{\mathrm{t}}^2}}=\mathrm{\epsilon }{\displaystyle \frac{\partial {\mathrm{D}}_3}{\partial {\mathrm{\Psi }}_3}}\left({\overline{\mathsf{\Omega }}}_3-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_3\right)+{\mathrm{\epsilon }}^2\dots \dots \dots \dots \dots \dots \dots \dots \\ {\displaystyle \frac{{\mathrm{d}}^2{\mathrm{\Psi }}_3}{\mathrm{d}{\mathrm{t}}^2}}=\mathrm{\epsilon }\left[-{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial {\mathsf{\eta }}_3}{\partial \mathrm{\tau }}}+{\displaystyle \frac{\partial {\mathrm{E}}_3}{\partial {\mathrm{\Psi }}_3}}\left({\overline{\mathsf{\Omega }}}_3-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_3\right)\right]+{\mathrm{\epsilon }}^2\dots \end{array}\right.. \end{gathered}$$

(18)

Using the first system of equations expressed in Equations (17) and (18) and the supposed solution (16) in Equation (14), balancing the coefficients of the mathematical expressions $\mathrm{\epsilon }\cos \left({\displaystyle \frac{1}{2}}{\mathrm{\Phi }}_1+{\mathrm{\Psi }}_1\right),\mathrm{\epsilon }\sin \left({\displaystyle \frac{1}{2}}{\mathrm{\Phi }}_1+{\mathrm{\Psi }}_1\right)$ from the LHS with those from the RHS, and neglecting the terms containing $\sin ({\displaystyle \frac{{\mathrm{k}}_1}{2}}{\mathrm{\Phi }}_1+{\mathrm{\Psi }}_1),\cos ({\displaystyle \frac{{\mathrm{k}}_1}{2}}{\mathrm{\Phi }}_1+{\mathrm{\Psi }}_1),{\mathrm{k}}_1=2,3,4$ results in the differential system of equations (DSOE) for the GT:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial {\mathrm{D}}_1}{\partial {\mathrm{\Psi }}_1}}\left({\overline{\mathsf{\Omega }}}_1-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_1\right)-2{\mathrm{\Phi }}_1{\overline{\mathsf{\Omega }}}_1{\mathrm{E}}_1={\mathrm{\mu }}_1{\mathrm{\Phi }}_1{\overline{\mathsf{\Omega }}}_1^2\cos \left(2{\mathrm{\Psi }}_1\right)-{\displaystyle \frac{3}{4}}{\mathrm{\Gamma }}_1{\mathrm{\Phi }}_1^3\\ {\mathrm{\Phi }}_1{\displaystyle \frac{\partial {\mathrm{E}}_1}{\partial {\mathrm{\Psi }}_1}}\left({\overline{\mathsf{\Omega }}}_1-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_1\right)+2{\overline{\mathsf{\Omega }}}_1{\mathrm{D}}_1=-{\mathrm{\mu }}_1{\mathrm{\Phi }}_1{\overline{\mathsf{\Omega }}}_1^2\sin \left(2{\mathrm{\Psi }}_1\right)-2\mathrm{\xi }{\overline{\mathsf{\Omega }}}_1^2{\mathrm{\Phi }}_1\end{array}\right.,$$

(19)

and proceeding with the same mathematical manipulations in Equation (15) by using the second system of equations expressed in Equations (17) and (18), and the supposed solution (16), balancing the coefficients of the mathematical expressions $\mathrm{\epsilon }\cos \left({\displaystyle \frac{1}{2}}{\mathrm{\Phi }}_3+{\mathrm{\Psi }}_3\right),\mathrm{\epsilon }\sin \left({\displaystyle \frac{1}{2}}{\mathrm{\Phi }}_3+{\mathrm{\Psi }}_3\right)$ from the LHS with those from the RHS, and neglecting the terms containing $\sin ({\displaystyle \frac{{\mathrm{k}}_2}{2}}{\mathrm{\Phi }}_3+{\mathrm{\Psi }}_3),\cos ({\displaystyle \frac{{\mathrm{k}}_2}{2}}{\mathrm{\Phi }}_3+{\mathrm{\Psi }}_3),$${\mathrm{k}}_2=2,3,4,5$ results in the differential system of equations (DSOE) for the GB:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial {\mathrm{D}}_3}{\partial {\mathrm{\Psi }}_3}}\left({\overline{\mathsf{\Omega }}}_3-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_3\right)-2{\mathrm{\Phi }}_3{\overline{\mathsf{\Omega }}}_3{\mathrm{E}}_3={\mathrm{\mu }}_3{\mathrm{\Phi }}_3{\overline{\mathsf{\Omega }}}_3^2\cos \left(2{\mathrm{\Psi }}_3\right)-{\displaystyle \frac{3}{4}}{\mathrm{\Gamma }}_3{\mathrm{\Phi }}_1^3-{\displaystyle \frac{5}{8}}{\mathrm{\Gamma }}_5{\mathrm{\Phi }}_3^5\\ {\mathrm{\Phi }}_1{\displaystyle \frac{\partial {\mathrm{E}}_3}{\partial {\mathrm{\Psi }}_3}}\left({\overline{\mathsf{\Omega }}}_3-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_3\right)+2{\overline{\mathsf{\Omega }}}_3{\mathrm{D}}_3=-{\mathrm{\mu }}_3{\mathrm{\Phi }}_3{\overline{\mathsf{\Omega }}}_3^2\sin \left(2{\mathrm{\Psi }}_3\right)-2\mathrm{\xi }{\overline{\mathsf{\Omega }}}_3^2{\mathrm{\Phi }}_3\end{array}\right.,$$

(20)

which yields the solutions of DSOE (19) for the GT and (20) for the GB:

$$\left\{\begin{array}{c}{\mathrm{D}}_1=-{\displaystyle \frac{{\mathrm{\mu }}_1{\overline{\mathsf{\Omega }}}_1^2{\mathrm{\Phi }}_1}{{\mathsf{\eta }}_1}}\sin 2{\mathrm{\Psi }}_1-\mathrm{\xi }{\overline{\mathsf{\Omega }}}_1{\mathrm{\Phi }}_1\\ {\mathrm{E}}_1=-{\displaystyle \frac{{\mathrm{\mu }}_1{\overline{\mathsf{\Omega }}}_1^2}{{\mathsf{\eta }}_1}}\cos 2{\mathrm{\Psi }}_1+{\displaystyle \frac{3}{8}}{\displaystyle \frac{{\mathrm{\Gamma }}_1{\mathrm{\Phi }}_1^2}{{\mathsf{\Omega }}_1}}\end{array},\right.\left\{\begin{array}{c}{\mathrm{D}}_3=-{\displaystyle \frac{{\mathrm{\mu }}_3{\overline{\mathsf{\Omega }}}_3^2{\mathrm{\Phi }}_3}{{\mathsf{\eta }}_3}}\sin 2{\mathrm{\Psi }}_3-\mathrm{\xi }{\overline{\mathsf{\Omega }}}_3{\mathrm{\Phi }}_3\\ {\mathrm{E}}_3=-{\displaystyle \frac{{\mathrm{\mu }}_3{\overline{\mathsf{\Omega }}}_3^2}{{\mathsf{\eta }}_3}}\cos 2{\mathrm{\Psi }}_3+{\displaystyle \frac{3}{8}}{\displaystyle \frac{{\mathrm{\Gamma }}_3{\mathrm{\Phi }}_3^2}{{\mathsf{\Omega }}_3}}+{\displaystyle \frac{5}{16}}{\displaystyle \frac{{\mathrm{\Gamma }}_5{\mathrm{\Phi }}_3^4}{{\mathsf{\Omega }}_3}}\end{array}.\right.$$

(21)

Using the system of Equation (21) in the DSOE (17) and considering all the characteristics of the systems in real-time values yields the DSOEs of the analytical amplitudes and the phases angles of FTV for GT and GB in the APPR for each element:

$$\begin{gathered} \left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_1}{\mathrm{d}\mathrm{t}}}=-{\displaystyle \frac{{\mathrm{\mu }}_1{\overline{\mathsf{\Omega }}}_1^2{\mathrm{\Phi }}_1}{{\mathsf{\eta }}_1}}\sin 2{\mathrm{\Psi }}_1-\mathrm{\xi }{\overline{\mathsf{\Omega }}}_1{\mathrm{\Phi }}_1\\ {\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_1}{\mathrm{d}\mathrm{t}}}={\overline{\mathsf{\Omega }}}_1-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_1-{\displaystyle \frac{{\mathrm{\mu }}_1{\overline{\mathsf{\Omega }}}_1^2}{{\mathsf{\eta }}_1}}\cos 2{\mathrm{\Psi }}_1+{\displaystyle \frac{3}{8}}{\displaystyle \frac{{\mathrm{\Gamma }}_1{\mathrm{\Phi }}_1^2}{{\mathsf{\Omega }}_1}}\end{array},\right. \\ \left\{\begin{array}{c}{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_3}{\mathrm{d}\mathrm{t}}}=-{\displaystyle \frac{{\mathrm{\mu }}_3{\overline{\mathsf{\Omega }}}_3^2{\mathrm{\Phi }}_3}{{\mathsf{\eta }}_3}}\sin 2{\mathrm{\Psi }}_3-\mathrm{\xi }{\overline{\mathsf{\Omega }}}_3{\mathrm{\Phi }}_3\\ {\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_3}{\mathrm{d}\mathrm{t}}}={\overline{\mathsf{\Omega }}}_3-{\displaystyle \frac{1}{2}}{\mathsf{\eta }}_3-{\displaystyle \frac{{\mathrm{\mu }}_3{\overline{\mathsf{\Omega }}}_3^2}{{\mathsf{\eta }}_3}}\cos 2{\mathrm{\Psi }}_3+{\displaystyle \frac{3}{8}}{\displaystyle \frac{{\mathrm{\Gamma }}_3{\mathrm{\Phi }}_3^2}{{\mathsf{\Omega }}_3}}+{\displaystyle \frac{5}{16}}{\displaystyle \frac{{\mathrm{\Gamma }}_5{\mathrm{\Phi }}_3^4}{{\mathsf{\Omega }}_3}}\end{array}.\right. \end{gathered}$$

(22)

The investigation of dynamic instability for the stationary FTV of AD consists of computing the dynamic instability boundary frontiers (DIBF) of GT and GB in the APPR. The stationary dynamic behavior of FTV for the GT and the GB is characterized by zeroing ${\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_1}{\mathrm{d}\mathrm{t}}},{\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_1}{\mathrm{d}\mathrm{t}}},{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_3}{\mathrm{d}\mathrm{t}}},{\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_3}{\mathrm{d}\mathrm{t}}}$ in DSOE (22) and expressing the trigonometric identity ${\sin }^2{\mathrm{\Psi }}_{\mathrm{i}}+{\cos }^2{\mathrm{\Psi }}_{\mathrm{i}}=1,\mathrm{i}=1,3$. This yields, for global tulip and global bowl,

$${\left({\displaystyle \frac{{\mathsf{\eta }}_1}{{\mathsf{\Omega }}_1}}\right)}^2\left[{\mathrm{\xi }}^2+{\left(1-{\displaystyle \frac{{\mathsf{\eta }}_1}{2{\mathsf{\Omega }}_1}}\right)}^2\right]-{\mathrm{\mu }}_1^2=0,{\left({\displaystyle \frac{{\mathsf{\eta }}_3}{{\mathsf{\Omega }}_3}}\right)}^2\left[{\mathrm{\xi }}^2+{\left(1-{\displaystyle \frac{{\mathsf{\eta }}_3}{2{\mathsf{\Omega }}_3}}\right)}^2\right]-{\mathrm{\mu }}_3^2=0,$$

(23)

in the spaces $\left({\mathsf{\eta }}_1,\mathrm{\xi },{\mathrm{\mu }}_1\right)$,$\left({\mathsf{\eta }}_3,\mathrm{\xi },{\mathrm{\mu }}_3\right)$, that is, the space (excitation frequency, damping ratio, excitation parameter).

As can be remarked for Equation (23), the instability frontiers for stationary FTV (SFTV) of GT and GB do not contain the amplitudes ${\mathrm{\Phi }}_1,{\mathrm{\Phi }}_3$ because the frontiers of instability are defined as the limits for which the amplitudes may initiate, and the amplitudes must be put to zero in APPR. For the dynamic instability analysis of nonstationary FTV (NFTV), the DSOE (22) was used to compute the spectral amplitudes (SA) and the spectral velocity amplitudes (SVA) of GT and GB in transition through APPR. For dynamic instability analysis in transition through APPR, the sweep of rate (SOR) of the excitation frequencies for global tulip and global bowl were considered to have logarithmic variations, yielding thus the equations

$${\mathsf{\eta }}_1={\overline{\mathsf{\eta }}}_1{\mathrm{\kappa }}_1^{\mathrm{t}(1+\mathrm{n})},{\mathrm{\kappa }}_1\ge 1,{\mathsf{\eta }}_3={\overline{\mathsf{\eta }}}_3{\mathrm{\kappa }}_3^{\mathrm{t}(1+\mathrm{n})},{\mathrm{\kappa }}_3\ge 1.$$

(24)

where ${\overline{\mathsf{\eta }}}_1,{\overline{\mathsf{\eta }}}_3$ are the starting frequencies of the logarithmic frequency balayage (LFB) for global tulip and global bowl, and n is the frequency’s SOR. To integrate the DSOE (22) numerically, MATLAB software (version R 2017) was developed based on function ode 45 (Runge–Kutta integration method of fifth order). Thus, the SA and SVA were computed, making the determination of the graph’s velocity amplitude versus amplitude for the GT and GB possible. The graph’s velocity amplitude versus amplitude represents the phase portrait of NFTV for GT and GB, which is the detection method of possible chaotic dynamic behavior.

4. Results and Discussion

Figure 7 illustrates the variation of the AD’s GT and GB natural frequencies based on Equation (6).

Comparing the obtained results with the published ones in [5] (pp. 119–139), the conclusion is clear that the DMADFTV-computed circular eigenfrequencies agree with the data in the literature. The potential sources of discrepancies between the theoretical predictions and experimental measurements mentioned before might be manufacturing tolerances, material inconsistencies, or assembly imperfections. These practical technical aspects could explain the differences between the simulated results and the experimental data published in [5] (pp. 119–139), which are less than 3% for the eigenfrequency of the GT and less than 4% for the GB. The free natural frequencies of the GT and the GB that are computed based on Equation (6) are not sensitive to the logarithmic decrement $\mathsf{\xi }$ that participates in Equations (14) and (15) in the RHS as an attenuation term of the excitations. The nonuniformities of the AIMM and AIGM for the GB and the GT are expressed in the terms ${\mathrm{\chi }}_{\mathrm{n}\mathrm{B}}$, ${\mathrm{\chi }}_{\mathrm{n}\mathrm{T}}$ (see Equations (A1)–(A3) in Appendix A) and ${\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}$, ${\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}$ (see Equations (A13) and (A14) in Appendix B), and the angles ${\beta }_1,{\beta }_2$ influence the free eigenfrequencies computed with Equation (6).

Based on Equation (23), MATLAB software was developed to calculate the dynamic instability boundary frontiers (DIBF) for the stationary GT’s FTV and GB’s FTV; the graphs are presented in Figure 8 and Figure 9.

Analyzing the graphs presented in Figure 8 and Figure 9, it can be remarked that the folded surface of the DIBF manifests symmetry concerning the plane (excitation frequency, damping ratio), having the same manifestation as for the global tulip and the global bowl. Because the damping ratio $\mathsf{\xi }$ is an isotropic physical property of the material, analyzing the needed limit values to obtain DIBF for the global tulip and the global bowl, the condition is that the steel composition must satisfy the condition $\mathsf{\xi }\ge 0.0044$ if it is used for the fabrication of the driveline of the same steel. It can also be remarked that the folded surfaces illustrated in Figure 8 and Figure 9 close, between their two branches, the instability volume in the space $\left({\mathsf{\eta }}_1,\mathsf{\xi },{\mathrm{\mu }}_1\right)$ for the global tulip and the space $\left({\mathsf{\eta }}_3,\mathsf{\xi },{\mathrm{\mu }}_3\right)$ for the global bowl. Figure 10, Figure 11 and Figure 12 illustrate the graphs for investigating the dynamic instability of NFTV for the global tulip in transition through APPR ${\mathsf{\eta }}_1\simeq 6131.6 \mathrm{H}\mathrm{z}$, considering the worst case for the global tulip ${\mathrm{\beta }}_1=15°$, ${\mathrm{\mu }}_1=0.1001$, ${\overline{\mathsf{\Omega }}}_1=19263 \mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$, ${\mathrm{\chi }}_{\mathrm{G}\mathrm{T}\mathrm{n}}=0.0073$, ${\mathrm{\Gamma }}_1=1.9298\cdot {10}^4 {\mathrm{s}}^{-2}\mathrm{r}\mathrm{a}{\mathrm{d}}^{-2}$, $\mathrm{n}=0.000008$ [octaves/min], and the damping ratio $\mathsf{\xi }=\mathrm{0.0016\dots 0.0043}$. MATLAB software was developed based on DSOE (22) for NFTV of the global tulip to calculate ${\mathrm{\Phi }}_1\left({\mathsf{\eta }}_1\right),{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_1}{\mathrm{d}\mathrm{t}}}\left({\mathsf{\eta }}_1\right)$.

For the graphs illustrated in Figure 10, Figure 11 and Figure 12, the excitation frequency has a logarithmic balayage for NFTV of the GT in transition through APPR. In contrast, the characteristics of internal and external excitations ${\mathsf{\mu }}_1,{\chi }_{\mathrm{G}\mathrm{T}\mathrm{n}},{\mathsf{\Gamma }}_1$ are constants. In Figure 10 and Figure 11, it is remarked that the nonstationary SA and SVA ${\mathrm{\Phi }}_1\left({\mathsf{\eta }}_1\right),$ ${\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_1}{\mathrm{d}\mathrm{t}}}\left({\mathsf{\eta }}_1\right)$ for NFTV of the GT manifest the “beatings” and a rapid decrease of the extreme values with the damping’s increases in the range $(\mathrm{1.6\dots \dots 4.3})\times {10}^{-3}$. The greatest value of the nonstationary SA and SVA for NFTV of the GT is reached when ${\mathsf{\eta }}_1=6136 \mathrm{H}\mathrm{z}$.

Figure 12 shows that the “look” is similar to the strange attractor’s representations. A chaos detection method was developed for AD by employing the FOAM for the GT’s NFTV; a chaos detection method was developed for AD’s NFTV.

For the graphs illustrated in Figure 13, Figure 14 and Figure 15, the excitation frequency has a logarithmic balayage for the GB’s NFTV in transition through APPR ${\eta }_3\simeq 5180 \mathrm{H}\mathrm{z}$, considering the worst case for the GB ${\mathrm{\beta }}_2=47°,$ ${\mathsf{\mu }}_3=0.0065,$ ${\overline{\mathsf{\Omega }}}_3=16273.5 \mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s},$ ${\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}}=0.0492,$ ${\mathrm{\Gamma }}_3=-9.9638\times {10}^4 {\mathrm{s}}^{-2}\mathrm{r}\mathrm{a}{\mathrm{d}}^{-2},$ ${\mathrm{\Gamma }}_5=-1.9928\times {10}^5 {\mathrm{s}}^{-2}\mathrm{r}\mathrm{a}{\mathrm{d}}^{-4},$ $\mathrm{n}=0.00001$ [octaves/min], and the damping ratio $\mathsf{\xi }=\mathrm{0.0016\dots 0.006}$. MATLAB software was developed based on DSOE (22) for nonstationary NFTV of the GB to calculate ${\mathrm{\Phi }}_3\left({\mathsf{\eta }}_3\right),{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_3}{\mathrm{d}\mathrm{t}}}\left({\mathsf{\eta }}_3\right)$.

For the graphs illustrated in Figure 13, Figure 14 and Figure 15, the excitation frequency has a logarithmic balayage for NFTV of the GB in transition through APRR. In contrast, the characteristics of internal and external excitations ${\mathsf{\mu }}_3,{\mathrm{\chi }}_{\mathrm{G}\mathrm{B}\mathrm{n}},{\mathrm{\Gamma }}_3,{\mathrm{\Gamma }}_5$ are constants. In Figure 13 and Figure 14, it is remarked that the nonstationary SA and SVA ${\mathrm{\Phi }}_3\left({\mathsf{\eta }}_3\right),{\displaystyle \frac{\mathrm{d}{\mathrm{\Phi }}_3}{\mathrm{d}\mathrm{t}}}\left({\mathsf{\eta }}_3\right)$ for NFTV of GB manifest the “beatings” and a rapid decrease of the extreme values with the damping’s increases in the range $(\mathrm{1.6\dots \dots 6.0})\times {10}^{-3}$. The greatest value of the nonstationary SA and SVA for NFTV of the GB is reached when ${\mathsf{\eta }}_3=5180 \mathrm{H}\mathrm{z}$.

It is remarked that the “pictures” in Figure 15 are like those specific to a strange attractor. In this way, the aim of using FOAM for NFTV of the GB was reached, FOAM being thus a detection method for chaos manifestation of NFTV for the elements of AD. Comparing Figure 10, Figure 11 and Figure 12 with Figure 13, Figure 14 and Figure 15 shows the same trend for the global tulip and the global bowl regarding the NFTV of AD’s elements in transition through APPR. A possible chaotic dynamic behavior for these elements was detected using the FOAM to investigate NFTV for AD’s elements.

5. Conclusions

This paper deals with the dynamic instability analysis for the AD in transition through APPR. To achieve its goals, first, the study modified an existent model for ADFTV by designing a DMADFTV of an AD that includes the nonuniformity of AIGM, the nonuniformity of AIMM, the quasi-homokinetic CVJ joints, and the external harmonic excitations induced by the gearbox and the reactive torque generated by the wheel. The specific assumptions made when reducing the AIMM and the AIGM to the longitudinal axis of the GB and GT do not impact the generality of the DMADFTV. Based on this DMADFTV, the stationary dynamic instability of the ADE (GT, GB) was investigated by computing the instability frontiers using FOAM. Using FOAM with the second-order differential Equation (9), expressing the GT’s FTV, and Equation (10), expressing the GB’s FTV, results in the nonstationary amplitudes and nonstationary phase angles in the APPR given by the DSOE (22). Based on the LFB with a specific SOR, the phase portraits NSVA–NSA for the GT and GB were computed, and this mathematical approach represented the detection method of chaotic AD’s FTV. Steinwede demonstrated in [4], using sophisticated experiments, that the chaotic consequences are the cracks, the micro-cracks, and the pitting that were detected on the tripod’s TTJ, the balls’s BIRJ, and the IR’s BIRJ; this is mandatory for the chaos manifestation control in the design of AD’s components. The FOAM used in this study represents only a method of detecting the possible chaotic manifestation of ADFTV. The confirmation of the chaotic manifestation of ADFTV is given only by a study of the chaotic manifestation based on an investigation of Lyapunov’s or Jacobi’s stability.

In future research, the authors hope to develop the use of Jacobi’s stability analysis and the KCC theory, as mentioned in [20,21,22,23], to control the chaotic manifestation of GT’s FTV and GB’s FTV in transition through other possible resonances such as internal, superharmonic, combination, and subharmonic, as mentioned by Nayfeh in [32].