**1. Introduction**

The present work presents a perturbation approach, specifically the asymptotic method approach (AMA) [1], to investigate the principal parametric resonance (PPR) for the forced bending vibrating (FBV) movements of the elements of an automotive driveshaft (AD), with the AMA being a powerful tool for the investigation of vibrations induced by shocks, as mentioned by Webber in the literature [2]. To analyze the PPR region, a model for FBV movements of the AD elements was designed, with a specific multibody dynamic structure considering the following aspects: the geometric nonuniformity of the AD elements and

**Citation:** Bugaru, M.; Vasile, O. Modeling and Analysis of FBV Movements for Automotive Driveshafts in the PPR Region. *Appl. Sci.* **2022**, *12*, 3237. https://doi.org/ 10.3390/app12073237

Academic Editors: Jin-Gyun Kim, Jae Hyuk Lim and Peter Persson

Received: 9 March 2022 Accepted: 21 March 2022 Published: 22 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the shock excitation induced by the road geometry. The excitations of the FBV movements of the AD elements are due to the shock excitations acting on the automotive wheel, being an excitation induced by the road geometry [3]. The geometric nonuniformities of the AD elements are undoubtedly the leading cause of nonlinear parametric vibrations of the AD in the range of 0.1–15 kHz, as established by experimental data in the literature [4] (pp. 98–123). The first researchers who considered the specific geometry of the AD elements were Mazzei and Scott. They first analyzed the nonlinear dynamic behavior of the AD elements in the PPR region [5]. The lateral vibration of the AD was analyzed by Browne and Palazzolo for superharmonic resonance in [6]. In [7], Xia et al., studied the bending vibrations of a shaft for 4WD (four-wheel drive) drivetrains to design the control methods at low frequencies. Analysis of the parametric bending vibrations for a cardan shaft was performed by Alugongo, and the data are presented in [8]. The NVH phenomena related to the vehicle driveline architecture were investigated by Wellmann et al. in [9] within different frequency ranges. At the same time, in [10], the same authors realized the driveline integration process.

In [11], Yang et al., performed dynamic analysis and vibration tests on the carbon-fiberreinforced plastics drive-line transmission for machine tools. Yao studied the vibration control of driveshafts considering the Sommerfeld effect on multiple linked shafts, and the data are presented in [12]. Jadhav et al. [13] conducted an extended analysis of the vibrations for a driveshaft with a crack using experimental modal analysis and finite element analysis (FEA) for isotropic materials such as structural steel. The increased interest in achieving a high standard of comfort in the automotive industry led to developing research on dynamic vibration absorbers for the bending vibration of the vehicle propeller shaft [14]. In the same area, Wu et al. investigated the use of photonic crystals in the vibration dampers of the AD [15]. In recent years, composite materials gained advancements in industrial utilization. Therefore, analysis of their dynamic behavior was needed, so Prakash and Sinha estimated the deflection and the natural frequencies of an AD made of such materials using the FEA [16]. The authors of [17] analyzed an AD carrying torque from the principal car differential to the rear wheel differential using the FEA for lateral and torsional vibration, considering the conventional and composite materials. Alam et al. in [18] performed a detailed investigation of the vibration characteristics of the composite AD using the finite element method (FEM). This paper aims to investigate the dynamic behavior of the AD elements for the FBV movements in the PPR region. This investigation supposes the design of a model for FBV movements for an AD's elements, the computation of the stationary and nonstationary amplitude spectrum in the PPR region, and the analysis of the dynamic instability frontiers in the same area. The novelty for this investigation is the design of the model of FBV movements of the AD elements and the use of a perturbation approach (AMA) to compute the amplitude, the phase angle, and the dynamic instability frontiers for both stationary and nonstationary motions. The designed model for FBV movements for an AD's elements contains the phenomena of geometric inertia nonuniformity of the AD and the shock excitation of the AD due to the road geometry. An AD is a multibody mechanical structure mechanism that allows transmitting a moment from the gearbox to the wheel, as shown in Figure 1. The principal components of such an AD, designed for heavy-duty SUVs, are presented in Figure 2 and consist of the bowl (a) fixed with the steering wheel and the midshaft (b) that links both the bowl and the tulip (c).

**Figure 1.** The ADs are mounted in the gearbox.

**Figure 2.** The details of an AD.

The elements of the tulip–tripod joint is presented in Figure 3, where the elements include the midshaft (a), on which the tripod (b) is mounted through the splines and linked with the tulip's bell (c), and the tulip ax (d), which is fixed in the gearbox through the splines.

**Figure 3.** Details of tulip-tripod joint (tripod fixed by splines on the midshaft).

The elements of the bowl–inner race joint can be seen in detail in Figure 4, showing the midshaft (a) on which is mounted the inner race (b) is fixed through splines, linked with the bowl's bell (c), and the bowl ax (d) that is fixed to the steering wheel through the splines.

**Figure 4.** Details of the bowl-inner race (inner race fixed by splines on the midshaft) joint.

#### **2. Computation of the Inertial Characteristics of the AD Elements for the FBV Movements**

To use the AMA, it is necessary to compute the equations of the FBV movements for the AD elements, so it must be established from the beginning the movements of each AD element. Therefore, it was considered that the tulip and the bowl had deflections and rotations specific to the rigid body movements. In contrast, the midshaft had continuous deflections and rotations typical to a continuous elastic beam. These assumptions agreed with the technical reality because the rigidity of the tulip and the bowl are more significant compared with the midshaft's rigidity. Each AD element has its referential system, presented in Figure 5, namely the cartesian systems X1Y1Z1, X2Y2Z2, X3Y3Z3.

**Figure 5.** Schematic representation of an automotive driveshaft.

The planes (X1Y1),(X2Y2),(X3Y3) are coplanar planes, so the axes Z1,Z2,Z3 are parallel, and the movements considered for each AD element for the FBV movements are as follows (see Figure 5):


To compute the equations of the FBV movements of the AD elements using Hamilton's principle, it was necessary to reduce the mass and geometric moments of inertia of the tulip and bowl to the cartesian system of reference of the midshaft X2Y2Z2 using the variation of the geometric moments of inertia concerning the concurrent axis and parallel axis (Steiner's theorem). The geometric tulip's design (see Figure 3) is composed of the tulip's bell (TB) and tulip ax (TA), while the geometric bowl's design (see Figure 4) is composed of the bowl's bell (BB) and bowl ax (BA). The geometric and mass moments of inertia for the tulip JY2T,JZ2T, IY2T,IZ2T reduced to the Y2,Z2 axes of the midshaft are given by the following equations:

$$\mathbf{J}\_{\mathbf{Y}\_2\mathbf{T}} = \mathbf{J}\_{\mathbf{Y}\_2\mathbf{T}\mathbf{B}} + \mathbf{J}\_{\mathbf{Y}\_2\mathbf{T}\mathbf{A}\mathbf{A}}.\tag{1}$$

$$\mathbf{J}\_{\mathbf{Z}\_2\mathbf{T}} = \mathbf{J}\_{\mathbf{Z}\_2\mathbf{T}\mathbf{B}} + \mathbf{J}\_{\mathbf{Z}\_2\mathbf{T}\mathbf{A}}.\tag{2}$$

$$\mathbf{I}\_{\rm Y\_2T} = \mathbf{J}\_{\rm Y\_2TB} \mathbf{\rho} \mathbf{L\_{TB}} + \mathbf{J}\_{\rm Y\_2TA} \mathbf{\rho} \mathbf{L\_{TA}} \tag{3}$$

$$\mathbf{J}\_{\mathbf{Z}\_2\mathbf{T}} = \mathbf{J}\_{\mathbf{Z}\_2\mathbf{T}\mathbf{B}}\mathbf{\rho}\mathbf{L}\_{\mathbf{TB}} + \mathbf{J}\_{\mathbf{Z}\_2\mathbf{TA}}\mathbf{\rho}\mathbf{L}\_{\mathbf{TA}\mathbf{A}} \tag{4}$$

$$\mathbf{J}\_{\rm Y\_2TB} = 0.5(\mathbf{J}\_{\rm IT} + \mathbf{J}\_{\rm ZT}) \left[ 1 + \sin^2 \beta\_1 + \chi\_{\rm nT} \cos(2\varphi\_1) \cos^2 \beta\_1 \right] + \frac{\mathbf{S}\_{\rm T} \mathbf{L}\_{\rm TB}^2}{12} \cos^2 \beta\_1 + \mathbf{S}\_{\rm T} (\mathbf{d}\_{\rm CT})^2,\tag{5}$$

$$\chi\_{\rm nT} = \frac{\mathbf{J}\_{\rm 1T} - \mathbf{J}\_{\rm 2T}}{\mathbf{J}\_{\rm 1T} + \mathbf{J}\_{\rm 2T}},\tag{6}$$

$$\mathbf{J}\_{\rm Y\_2TA} = \frac{\pi \mathbf{d}\_{\rm TA}^4}{64} \left( 1 + \sin^2 \beta\_1 \right) + \frac{\pi \mathbf{d}\_{\rm TA}^2}{4} \frac{\mathbf{L}\_{\rm TA}^2}{12} \cos^2 \beta\_1 + \frac{\pi \mathbf{d}\_{\rm TA}^2}{4} (\mathbf{L}\_{\rm TB} + 0.5 \mathbf{L}\_{\rm TA})^2,\tag{7}$$

$$\mathbf{J}\_{\rm Z\_2TB} = 0.5(\mathbf{J}\_{\rm I\Gamma} + \mathbf{J}\_{\rm Z\Gamma})[1 - \chi\_{\rm n\Gamma}\cos(2\varphi\_1)] + \frac{\mathbf{S}\_{\rm T}\mathbf{L}\_{\rm TB}^2}{12} + \mathbf{S}\_{\rm T}(\mathbf{d}\_{\rm CT})^2,\tag{8}$$

$$\mathbf{J}\_{\rm Z\_2TA} = \frac{\pi \mathbf{d}\_{\rm TA}^4}{64} + \frac{\pi \mathbf{d}\_{\rm TA}^2}{4} \frac{\mathbf{L}\_{\rm TA}^2}{12} + \frac{\pi \mathbf{d}\_{\rm TA}^2}{4} (\mathbf{L}\_{\rm TB} + 0.5 \mathbf{L}\_{\rm TA})^2,\tag{9}$$

where J1T, J2T are the principal geometric moments of inertia for the tulip's bell, JY2TB, JY2TA, JZ2TB,JZ2TA are the geometric moments of inertia for the tulip's elements reduced to the axes Y2, Z2, ρ is the mass density of the AD elements, dCT is the distance between the center mass of the tulip and the tripod's center mass, ST is the surface of the cross-section for the tulip's bell, *χ*nT is the nonuniformity of the geometric moments of inertia for the tulip (see Figure 3), LTB is the length of the tulip's bell, LTA is the length of the tulip ax, and dTA is the diameter of the tulip ax. The geometric and mass moments of inertia for the bowl JY2B, JZ2B, IY2B, IZ2B reduced to the Y2, Z2 axes of the midshaft are given by the following equations:

$$\mathbf{J}\_{\mathbf{Y}\_2\mathbf{B}} = \mathbf{J}\_{\mathbf{Y}\_2\mathbf{B}\mathbf{B}} + \mathbf{J}\_{\mathbf{Y}\_2\mathbf{B}\mathbf{A}}.\tag{10}$$

$$\mathbf{J}\_{\mathbf{Z}\mathbf{2}\mathbf{B}} = \mathbf{J}\_{\mathbf{Z}\mathbf{2}\mathbf{B}\mathbf{B}} + \mathbf{J}\_{\mathbf{Z}\mathbf{2}\mathbf{B}\mathbf{A}\mathbf{A}}.\tag{11}$$

$$\mathbf{I}\_{\mathbf{Y}\_2\mathbf{B}} = \mathbf{J}\_{\mathbf{Y}\_2\mathbf{B}\mathbf{B}} \boldsymbol{\uprho} \mathbf{L}\_{\mathbf{B}\mathbf{B}} + \mathbf{J}\_{\mathbf{Y}\_2\mathbf{B}\mathbf{A}} \boldsymbol{\uprho} \mathbf{L}\_{\mathbf{B}\mathbf{A}\mathbf{A}} \tag{12}$$

$$\mathbf{I}\_{\mathbf{Z}\_2\mathbf{B}} = \mathbf{J}\_{\mathbf{Z}\_2\mathbf{B}\mathbf{B}} \boldsymbol{\uprho} \mathbf{I}\_{\mathbf{B}\mathbf{B}} + \mathbf{J}\_{\mathbf{Z}\_2\mathbf{B}\mathbf{A}} \boldsymbol{\uprho} \mathbf{I}\_{\mathbf{B}\mathbf{A}\mathbf{A}}.\tag{13}$$

$$\mathbf{J}\_{\rm Y\_2BB} = 0.5(\mathbf{J}\_{\rm IB} + \mathbf{J}\_{\rm BB}) \left[ 1 + \sin^2 \beta\_2 + \chi\_{\rm nB} \cos(2\varphi\_3) \cos^2 \beta\_2 \right] + \frac{\mathbf{S}\_{\rm B} \mathbf{L}\_{\rm BB}^2}{12} \cos^2 \beta\_2 + \mathbf{S}\_{\rm B} (\mathbf{d}\_{\rm CB})^2,\tag{14}$$

$$\chi\_{\rm nB} = \frac{\mathbf{J}\_{\rm 1B} - \mathbf{J}\_{\rm 2B}}{\mathbf{J}\_{\rm 1B} + \mathbf{J}\_{\rm 2B}},\tag{15}$$

$$\mathbf{J}\_{\rm Y\_2BA} = \frac{\pi \mathbf{d}\_{\rm BA}^4}{64} \left( 1 + \sin^2 \beta\_2 \right) + \frac{\pi \mathbf{d}\_{\rm BA}^2}{4} \frac{\mathbf{L}\_{\rm BA}^2}{12} \cos^2 \beta\_2 + \frac{\pi \mathbf{d}\_{\rm BA}^2}{4} (\mathbf{L}\_{\rm BB} + 0.5 \mathbf{L}\_{\rm BA})^2,\tag{16}$$

$$\mathrm{J}\_{\mathrm{Z\_2BB}} = 0.5(\mathrm{J}\_{\mathrm{1B}} + \mathrm{J}\_{\mathrm{2B}})[1 - \chi\_{\mathrm{nB}}\cos(2\,\varphi\_{\mathrm{3}})] + \frac{\mathrm{S\_{\mathrm{B}}}\mathrm{L}\_{\mathrm{BB}}^2}{12} + \mathrm{S\_{\mathrm{B}}}(\mathrm{d\_{CB}})^2,\tag{17}$$

$$\mathbf{J}\_{\rm Z\_2BA} = \frac{\pi \mathbf{d}\_{\rm BA}^4}{64} + \frac{\pi \mathbf{d}\_{\rm BA}^2}{4} \frac{\mathbf{L}\_{\rm BA}^2}{12} + \frac{\pi \mathbf{d}\_{\rm BA}^2}{4} (\mathbf{L}\_{\rm BB} + 0.5 \mathbf{L}\_{\rm BA})^2,\tag{18}$$

where J1B, J2B are the principal geometric moments of inertia for the bowl's bell, JY2BB, JY2BA, JZ2BB, JZ2BA are the geometric moments of inertia for the bowl's elements reduced to the axes Y2,Z2, dCB is the distance between the center mass of the bowl and the inner race's center mass, SB is the surface of the cross-section for the bowl's bell, *χ*nB is the nonuniformity of the geometric moments of inertia for the bowl (see Figure 4), LBB is the length of the bowl's bell, LBA is the length of the bowl ax, and dBA is the diameter of the bowl ax. By analyzing Equations (1)–(18), it can be remarked that the geometric and mass moments of inertia of the tulip and the bowl, reduced to the referential of the midshaft, are functions of the angles ϕ1, ϕ3, Φ1, Φ3, β1, β2, meaning that they depend on the dynamic behavior of the tulip and the bowl, with this aspect being a novelty in designing dynamic models for the AD elements' movements.

#### **3. The Dynamic Model for FBV Movements of an AD's Elements**

The dynamic model for the FBV (DMFFBV) movements of an AD's elements is presented in Figures 6 and 7.

**Figure 6.** Tulip-tripod joint part of DMFFBV: tulip, tulip-tripod joint, and midshaft.

**Figure 7.** Bowl-inner race joint part of DMFFBV: bowl, bowl-inner race joint, and midshaft.

By surveying the literature already presented in the introduction, it can be remarked that there has not been detailed research performed for the dynamic behavior of each AD element. The dynamic model enhanced by the present work allows the investigation of the dynamic behavior of each component of the AD, the internal mechanisms of excitation transmissions through the AD elements, as well as the causes of producing pitting and micro-cracks in the AD elements, as revealed by the experimental data in the literature [4] (pp. 88–94). The DMFFBV movements, presented in Figures 6 and 7, have three elements: a tulip–midshaft–bowl linked through two joints, the tulip-tripod joint (mounted on the midshaft (see Figure 3)), and the bowl–inner race joint (mounted at the other edge of the midshaft (see Figure 4)). These have the following dynamic characteristics:

1. The tulip has a stiffness k1, given by the serial springs k11, k12 (see Figure 6), and damping c1, given by the serial dampers c11, c12, for the bending vibration rigid movement of the tulip regarding the axis Z2 and a stiffness kt1, given by the serial springs kt11 , kt12 , as well as damping ct1, given by the serial dampers ct11 , ct12 , for the angular bending vibration rigid movement of the tulip regarding the axis Y2, given by the following relations:

$$\mathbf{k}\_{11} = \frac{3\mathbf{E}|\_{\mathbf{Z}\_2\mathbf{T}\mathbf{A}}}{\mathbf{L}\_{\mathbf{T}\mathbf{A}}^3}, \mathbf{k}\_{12} = \frac{3\mathbf{E}|\_{\mathbf{Z}\_2\mathbf{T}\mathbf{B}}}{\mathbf{L}\_{\mathbf{T}\mathbf{B}}^3}, \mathbf{k}\_1 = \frac{\mathbf{k}\_{11}\mathbf{k}\_{12}}{\mathbf{k}\_{11} + \mathbf{k}\_{12}}, \mathbf{c}\_1 = \frac{2\boldsymbol{\Delta}}{\sqrt{4\pi^2 + \boldsymbol{\Delta}^2}}\sqrt{\mathbf{k}\_1\mathbf{m}\_{\mathbf{T}\mathbf{A}}} \tag{19}$$

$$\mathbf{k}\_{\rm t1} = \frac{\mathbf{G} \mathbf{J}\_{\rm Y\_2 \rm TA}}{\mathbf{L}\_{\rm TA}}, \mathbf{k}\_{\rm t2} = \frac{\mathbf{G} \mathbf{J}\_{\rm Y\_2 \rm TB}}{\mathbf{L}\_{\rm TB}}, \mathbf{k}\_{\rm t1} = \frac{\mathbf{k}\_{\rm t1} \mathbf{k}\_{\rm t2}}{\mathbf{k}\_{\rm t1} + \mathbf{k}\_{\rm t2}}, \mathbf{c}\_{\rm t1} = \frac{2\boldsymbol{\Delta}}{\sqrt{4\pi^2 + \boldsymbol{\Delta}^2}} \sqrt{\mathbf{k}\_{\rm t1} \mathbf{l}\_{\rm Y\_2} \mathbf{T}}, \tag{20}$$

where E is Young's modulus, G is the shearing modulus, mT is the tulip's mass, and Δ is the logarithmic decrement of the free bending vibrations of the tulip (Δ = 0.001–0.01) [19,20].

2. The uniform midshaft (see Figures 2–7) in continuous FBV movement is assimilated with a uniform Timoshenko beam simply supported at both ends by elastic supports (the tulip-tripod and inner race-bowl joints are elastic supports for the midshaft), having at x = 0 a tripod (see Figure 3) fixed on the midshaft through splines and elastically linked in the tulip-tripod joint with the tulip (see Figure 6) and on the left-hand side at x = LMs an inner race (see Figure 4) fixed on the midshaft through splines and elastically linked in the bowl-inner race joint with the bowl (see Figure 7), with the inertial characteristics given by the relations below:

$$\mathbf{J}\_{\rm Y\_2Tr} = 0.5(\mathbf{J}\_{\rm ITr} + \mathbf{J}\_{\rm ZTr})[1 + \chi\_{\rm nTr}\cos(2\,\varphi\_2)],\\\mathbf{J}\_{\rm Z\_2Tr} = 0.5(\mathbf{J}\_{\rm ITr} + \mathbf{J}\_{\rm ZTr})[1 - \chi\_{\rm nTr}\cos(2\,\varphi\_2)],\tag{21}$$

$$\mathbf{J}\_{\mathbf{Y}\_{2}\mathbf{I}\mathbf{r}} = 0.5(\mathbf{J}\_{\text{1Ir}} + \mathbf{J}\_{\text{2Ir}})[1 + \chi\_{\text{nIr}}\cos(2\varphi\_{2})],\\\mathbf{J}\_{\mathbf{Z}\_{2}\mathbf{I}\mathbf{r}} = 0.5(\mathbf{J}\_{\text{1Ir}} + \mathbf{J}\_{\text{2Ir}})[1 - \chi\_{\text{nIr}}\cos(2\varphi\_{2})].\tag{22}$$

where J1Tr, J2Tr are the principal geometric moments of inertia for the tripod, J1Ir, J2Ir are the principal geometric moments of inertia for the inner race, *χ*nTr and *χ*nIr are the geometric nonuniformities of the tripod and inner race, and JY2Tr,JZ2Tr,JY2Ir,JZ2Ir are the geometric moments of inertia of the tripod and inner race concerning the axes Y2, Z2;

3. The bowl has a stiffness k2, given by the serial springs k21, k22 (see Figure 7), and damping c2, given by the serial dampers c21, c22, for the bending vibration rigid movement of the bowl regarding the axis Z2 and a stiffness kt2, given by the serial springs kt21 , kt22 , as well as damping ct2, given by the serial dampers ct21 , ct22 , for the angular bending vibration rigid movement of the bowl regarding the axis Y2, given by the following relations:

$$\mathbf{k}\_{21} = \frac{\mathbf{3E}\mathbf{J}\_{\mathbf{Z}\cdot\mathbf{B}\mathbf{A}}}{\mathbf{L}\_{\mathbf{B}\mathbf{A}}^{3}}, \mathbf{k}\_{22} = \frac{\mathbf{3E}\mathbf{J}\_{\mathbf{Z}\cdot\mathbf{B}\mathbf{B}}}{\mathbf{L}\_{\mathbf{B}\mathbf{B}}^{3}}, \mathbf{k}\_{2} = \frac{\mathbf{k}\_{21}\mathbf{k}\_{22}}{\mathbf{k}\_{21} + \mathbf{k}\_{22}}, \mathbf{c}\_{2} = \frac{2\boldsymbol{\Delta}}{\sqrt{4\pi^{2} + \boldsymbol{\Delta}^{2}}}\sqrt{\mathbf{k}\_{2}\mathbf{m}\_{\mathbf{B}\prime}}\tag{23}$$

$$\mathbf{k}\_{\mathbf{l}\_{\rm 21}} = \frac{\mathbf{G} \mathbf{J}\_{\rm Y\_2 \rm BA}}{\mathbf{L}\_{\rm BA}}, \mathbf{k}\_{\rm t\_{\rm 22}} = \frac{\mathbf{G} \mathbf{J}\_{\rm Y\_2 \rm BB}}{\mathbf{L}\_{\rm BB}}, \mathbf{k}\_{\rm t2} = \frac{\mathbf{k}\_{\rm t\_{\rm 21}} \mathbf{k}\_{\rm t\_{\rm 22}}}{\mathbf{k}\_{\rm t\_{\rm 21}} + \mathbf{k}\_{\rm t\_{\rm 22}}}, \mathbf{c}\_{\rm 2} = \frac{2\boldsymbol{\Delta}}{\sqrt{4\pi^2 + \boldsymbol{\Delta}^2}} \sqrt{\mathbf{k}\_{\rm 2} \mathbf{I}\_{\rm Y\_2} \mathbf{B}}, \tag{24}$$

where mB is the bowl's mass.

4. The tulip–tripod joint in FBV movement (see Figures 3 and 6) realizes the elastic link between the tulip and the midshaft through the stiffness kTTr and the damping cTTr for the bending vibrating movements concerning the Z2 axis and the angular stiffness ktTTr and angular damping ctTTr for the bending vibrating movements concerning Y2 axis.

5. The bowl–inner race joint in FBV movement (see Figures 4 and 7) realizes the link between the bowl and the midshaft through the stiffness kIrB and the damping cIrB for the bending vibrating movements concerning the Z2 axis and the angular stiffness ktIrB and the angular damping ctIrB for the bending vibrating movements concerning the Y2 axis. The wheel induces excitations as a moderate impulsive shock force Fs acting in the Z2 direction, and the excitation load can be expressed as

$$\mathbf{F\_S = \overline{F}\_S \left[ 1 + \mathbf{q\_3} t^{\mathbf{q\_1}} e^{-\mathbf{q\_2}t} \right]},\tag{25}$$

where FS is the amplitude of the shock on the bowl's longitudinal axis X3 transmitted from the wheel axis and qi , i = 1, 3, are experimental constants, depending on the type of shock applied at the wheel by the road excitation [4] (pp. 142–172). All the computations were performed considering that the tulip was a cantilever beam fixed in the gearbox with simple elastic supports in the tulip–tripod joint. The bowl was a cantilever beam fixed in the steering wheel by simple elastic supports in the bowl-inner race joint. As mentioned in the literature [4] (pp. 142–172), the shock excitation loads produce huge automotive stress solicitations in the car suspension system and the rim-tire system. These two systems can absorb 90% of the shock energy. Therefore, only 10% of the shock acts on the AD elements as a variation of the quantity of movement during a very short time, estimated at 0.001 s. Because of this phenomenon, the shock is distributed to each AD element, inducing for the bowl stress in the region of the splines, for the midshaft the initial velocity of the "ends", and for the tulip stress in the region of the splines (see Figure 2). Due to this mechanism of excitation, it will be considered that the stress in the tulip's ax and bowl's ax imposes the geometric capability dimensions of the geometric elements dTA, dBA. At the same time, the midshaft dynamic behavior will be an excitation element, like the internal tuned damper (ITD) mentioned in the literature [6], for the tulip and the bowl. The novelty of this DMFFBV of an AD's elements, as can be noted from Equations (19)–(25), consists of linking the rigid movements of the bowl and the tulip with the continuous movement of the midshaft. In addition, another novelty is that the dynamic characteristics

concerning stiffness and damping are expressed for the bending displacements and the bending rotation displacements for each AD element, and the serial stiffness and damping of the bowl and tulip are considered, keeping in mind that each of these elements has two distinctive parts: the ax and "bell". In the meantime, the stiffness and damping for the tuliptripod and bowl-inner race joints were considered for the relative bending displacements and the relative bending rotation displacements of the movements between the bell of the tulip and the tripod, as well as for the inner race and the bell of the bowl.

## **4. The Equations of FBV Movements of AD Elements Induced by Shock Road Excitations**

For the DMFFBV movements of AD elements presented in Figures 6 and 7, using the Hamilton's principle [21] (pp. 265–269) yields

$$\delta \int\_{\mathbf{P}\_1 \mathbf{P}\_2} \mathcal{L}(\mathbf{q}\_1, \dots, \mathbf{q}\_{\mathbf{n}'} \mathbf{q}\_{1'}, \dots, \mathbf{q}\_{\mathbf{n}'} \mathbf{t}) \mathbf{dt} = \mathbf{0},\tag{26}$$

where Lagrange's function <sup>L</sup>(q1, ...., qn, • q1, ....., • qn, t) depends on the generalized coordinates q1, ...., qn and the generalized velocities • q1, ....., • qn, while P1, P2 are two points in the spatial configurations <sup>−</sup>−−−−−−−→ • q1, ....., • qn <sup>=</sup> <sup>→</sup> Ξ(q1, ...., qn). Lagrange's function is given by the following equation:

$$\mathbf{L} = \mathbf{T} + \prod\_{\prime} \mathbf{I}\_{\prime} \tag{27}$$

where the potential energy ∏ for the DMFFBV movements of AD elements (see Figures 6 and 7) is given by the following generalized equation [22] (pp. 371–376) [23] (pp. 734–739):

∏ = <sup>1</sup> 2 LMs 0 EJY2Ms *<sup>∂</sup>*Φ<sup>2</sup> *∂*x 2 + kAG *<sup>∂</sup>*w2 *<sup>∂</sup>*<sup>x</sup> − Φ<sup>2</sup> 2 2 dx + <sup>1</sup> 2 c1 • w1 2 + c2 • w3 2 <sup>+</sup> ct1- • Φ<sup>1</sup> 2 <sup>+</sup> ct2- • Φ<sup>3</sup> 2 + 1 2 cTTr • w1 <sup>−</sup> *<sup>∂</sup>*w2 *<sup>∂</sup>*<sup>t</sup> (0, t) 2 <sup>+</sup> ctTTr- • <sup>Φ</sup><sup>1</sup> cos <sup>β</sup><sup>1</sup> <sup>−</sup> *<sup>∂</sup>*Φ<sup>2</sup> *<sup>∂</sup>*<sup>t</sup> (0, t) 2 <sup>+</sup> cIrB *<sup>∂</sup>*w2 *<sup>∂</sup>*<sup>t</sup> (LMs, t) <sup>−</sup> • w3 2 <sup>+</sup> ctIrB- *∂*Φ<sup>2</sup> *<sup>∂</sup>*<sup>t</sup> (LMs, t) <sup>−</sup> • Φ<sup>3</sup> cos β<sup>2</sup> 2 + + <sup>1</sup> 2 # k1w2 <sup>1</sup> + k2w2 <sup>3</sup> + kt1Φ<sup>2</sup> <sup>1</sup> + kt2Φ<sup>2</sup> 3 \$ + (28)

$$\frac{1}{2} + \frac{1}{2} \left[ \mathbf{k\_{T}Tr} \left( \mathbf{w}\_{1} - \mathbf{w}\_{2}(\mathbf{0}, \mathbf{t}) \right)^{2} + \mathbf{k\_{OTr}} \left( \boldsymbol{\Phi}\_{1} \cos \boldsymbol{\mathfrak{x}}\_{1} - \boldsymbol{\Phi}\_{2}(\mathbf{0}, \mathbf{t}) \right)^{2} + \mathbf{k\_{aR}} \left( \mathbf{w}\_{2}(\mathbf{L}\_{\text{Ma}}, \mathbf{t}) - \mathbf{w}\_{3} \right)^{2} \\ + \mathbf{k\_{aR}} \left( \boldsymbol{\Phi}\_{2} (\mathbf{L}\_{\text{Ma}}, \mathbf{t}) - \boldsymbol{\Phi}\_{2} \cos \boldsymbol{\mathfrak{x}}\_{2} \right)^{2} \right], \quad \mathbf{x}\_{i} \in \mathbb{R}^{3}$$

where A is the cross-section area of the midshaft, w2(x, t) is the bending deflection (including the shear deformation) of the midshaft concerning the Z2 axis, Φ<sup>2</sup> (x, t) is the rotation of the cross-section of the midshaft, and concerning the Y2 axis, due only to the pure bending deflection, k is the shear correction factor, which in the literature [24] is in the range of 0.64–0.846, LMs is the length of the midshaft, and JY2Ms is the geometric moment of inertia of the midshaft concerning the Y2 direction given by the following equations:

$$\mathbf{J}\_{\mathbf{Y}\_2\mathbf{M}\mathbf{s}} = \frac{\pi \mathbf{d}\_{\mathbf{M}\mathbf{s}}^4}{64},\tag{29}$$

$$\mathbf{J}\_{\rm Y\_2Ms} = \frac{\pi \left(\mathbf{d}\_{\rm cMs}^4 - \mathbf{d}\_{\rm iMs}^4\right)}{64},\tag{30}$$

where the midshaft is considered to have a circular or tubular uniform cross-section with a diameter dMs for the circular cross-section or the diameters deMs, diMs for the tubular cross-section. In Equation (28), the generalized Rayleigh's dissipation function [23] (p. 611) specific to the mathematical formulations of the Euler–Lagrange generalized approach was added due to the presence of dampers in the DMFFBV movements, given by the following equation:

$$\begin{split} \Lambda &= \frac{1}{2} \Big[ \mathbf{c}\_{1} \Big( \mathbf{W}\_{1} \Big)^{2} + \mathbf{c}\_{2} \Big( \mathbf{w}\_{3} \Big)^{2} + \mathbf{c}\_{3} \Big( \mathbf{\varPhi}\_{1} \Big)^{2} + \mathbf{c}\_{2} \Big( \mathbf{\varPhi}\_{3} \Big)^{2} + \mathbf{c}\_{\mathrm{Tlr}} \Big( \mathbf{\varPhi}\_{1} - \frac{\partial \mathbf{w}\_{3}}{\partial t}(0, \mathbf{t}) \Big)^{2} + \mathbf{c}\_{\mathrm{Tlr}} \Big( \mathbf{\varPhi}\_{1} \cos \beta\_{1} - \frac{\partial \mathbf{s}\_{2}}{\partial t}(0, \mathbf{t}) \Big)^{2} \Big] + \\ &+ \frac{1}{2} \Big[ \mathbf{c}\_{\mathrm{H}} \Big( \frac{\partial \mathbf{w}\_{2}}{\partial t}(\mathbf{L}\_{\mathrm{Mb}, \mathrm{t}}) - \mathbf{w}\_{3} \Big)^{2} + \mathbf{c}\_{\mathrm{H}\mathrm{B}} \Big( \frac{\partial \mathbf{s}\_{2}}{\partial t}(\mathbf{L}\_{\mathrm{Mb}, \mathrm{t}}) - \boldsymbol{\Phi}\_{3} \cos \beta\_{2} \Big)^{2} \Big]. \end{split} \tag{31}$$

The kinetic energy of the DMFFBV movements (see Figures 6 and 7) for an AD's elements is given by the following generalized equation [22] (p. 374) [23] (p. 721):

$$\begin{array}{c} \bullet \text{ T} = \frac{\bullet}{2\,\text{m}\,\text{r}\,\text{w}\_{1}^{2} + \frac{1}{2}\,\text{m}\,\text{r}\_{\text{f}}\left(\frac{\partial\text{w}\_{2}}{\partial\text{f}}(0,\text{t})\right)^{2} + \frac{\bullet}{2\,\text{m}\,\text{w}\_{3}^{2} + \frac{1}{2}\,\text{m}\,\text{r}\_{\text{f}}\left(\frac{\partial\text{w}\_{2}}{\partial\text{f}}(\text{L}\_{\text{Ms}},\text{t})\right)^{2} + \frac{1}{2}\,\text{I}\_{\text{Y}\_{2}}\bullet\text{s}\_{\text{f}}^{2} +\\ \bullet \text{ } + \frac{1}{2}\text{I}\_{\text{Y}\_{2}\text{Tr}}\left(\frac{\partial\text{\theta}\_{2}}{\partial\text{t}}(0,\text{t})\right)^{2} + \frac{1}{2}\text{I}\_{\text{Y}\_{2}\text{B}}\bullet\text{s}\_{3}^{2} + \frac{1}{2}\text{I}\_{\text{Y}\_{2}\text{Im}}\left(\frac{\partial\text{\theta}\_{2}}{\partial\text{t}}(\text{L}\_{\text{Ms}},\text{t})\right)^{2} + \int\_{0}^{\bullet} \frac{1}{2}\left[\rho\text{A}\left(\frac{\partial\text{w}\_{2}}{\partial\text{t}}\right)^{2} + \text{I}\_{\text{Y}\_{2}\text{Ms}}\left(\frac{\partial\text{\theta}\_{2}}{\partial\text{t}}\right)^{2}\right] d\text{x},\tag{32} \end{array} \tag{33}$$

where mTr is the mass of the tripod and mIr is the mass of the inner race (see Figures 6 and 7). Several mathematical manipulations that include integration by parts of the nonlinear system of equations with partial derivatives of the second degree in the unknowns w1(t), Φ1(t), w2(x, t), Φ2(x, t), and w3(t), Φ3(t), yield

$$\mathbf{r}\mathbf{r}\mathbf{r}\mathbf{w}\_1^\bullet + \mathbf{c}\_1\mathbf{w}\_1 + \mathbf{c}\mathbf{r}\mathbf{r}\left(\mathbf{w}\_1^\bullet - \frac{\partial \mathbf{w}\_2}{\partial \mathbf{t}}(\mathbf{0}, \mathbf{t})\right) + \mathbf{k}\_1\mathbf{w}\_1 + \mathbf{k}\_{\mathrm{T}\mathrm{Tr}}(\mathbf{w}\_1 - \mathbf{w}\_2(\mathbf{0}, \mathbf{t})) = \mathbf{0},\tag{33}$$

$$\mathbf{d}\_{\mathbf{Y}\_{2}\mathbf{T}}\overset{\bullet}{\Phi}\_{\mathbf{I}} + \mathbf{c}\_{\text{t}\mathbf{I}}\overset{\bullet}{\Phi}\_{\mathbf{I}} + \mathbf{c}\_{\text{t}\mathbf{T}\mathbf{T}\mathbf{r}}\left(\overset{\bullet}{\Phi}\_{\mathbf{I}}\cos\beta\_{\mathbf{I}} - \frac{\partial\Phi\_{\mathbf{2}}}{\partial\mathbf{t}}(\mathbf{0},\mathbf{t})\right) + \mathbf{k}\_{\text{t}\mathbf{I}}\Phi\_{\mathbf{I}} + \mathbf{k}\_{\text{t}\mathbf{T}\mathbf{T}\mathbf{r}}(\Phi\_{\mathbf{1}}\cos\beta\_{\mathbf{I}} - \Phi\_{\mathbf{2}}(\mathbf{0},\mathbf{t})) = \mathbf{0},\tag{34}$$

$$
\rho \mathbf{A} \frac{\partial^2 \mathbf{w}\_2}{\partial \mathbf{t}^2} - \mathbf{k} \mathbf{A} \mathbf{G} \left[ \frac{\partial^2 \mathbf{w}\_2}{\partial \mathbf{x}^2} - \frac{\partial \Phi\_2}{\partial \mathbf{x}} \right] = \mathbf{0},\tag{35}
$$

$$\mathrm{I}\_{\mathrm{Y\_2Ms}} \frac{\partial^2 \Phi\_2}{\partial \mathbf{t}^2} - \mathrm{EJ}\_{\mathrm{Y\_2Ms}} \frac{\partial^2 \Phi\_2}{\partial \mathbf{x}^2} - \mathrm{kAG} \left(\frac{\partial \mathbf{w\_2}}{\partial \mathbf{x}} - \Phi\_2\right) = 0,\tag{36}$$

$$\mathbf{m}\_{\rm B} \overset{\scriptstyle \mathbf{w}\_{\rm B}}{\mathbf{w}\_{\rm B}} + \mathbf{c}\_{2} \overset{\scriptstyle \mathbf{w}\_{\rm B}}{\mathbf{w}\_{\rm B}} + \mathbf{c}\_{\rm IrB} \left(\overset{\scriptstyle \mathbf{w}\_{\rm B}}{\mathbf{w}\_{\rm B}} - \frac{\partial \mathbf{w}\_{\rm B}}{\partial \mathbf{t}} (\mathbf{L}\_{\rm M5} \mathbf{t})\right) + \mathbf{k}\_{2} \mathbf{w}\_{\rm 3} + \mathbf{k}\_{\rm IrB} (\mathbf{w}\_{\rm 3} - \mathbf{w}\_{\rm 2} (\mathbf{L}\_{\rm M5} \mathbf{t})) = \mathbf{0},\tag{37}$$

$$\mathbf{L}\_{\mathbf{Y}\_{2}\mathbf{B}}\overset{\bullet}{\Phi}\_{3} + \mathbf{c}\_{\text{t}1}\overset{\bullet}{\Phi}\_{1} + \mathbf{c}\_{\text{dirB}}\left(\overset{\bullet}{\Phi}\_{3}\cos\beta\_{2} - \frac{\partial\Phi\_{2}}{\partial\mathbf{t}}(\mathbf{L}\_{\text{M5s}}\mathbf{t})\right) + \mathbf{k}\_{\text{t2}}\Phi\_{3} + \mathbf{k}\_{\text{dirB}}(\Phi\_{3}\cos\beta\_{2} - \Phi\_{2}(\mathbf{L}\_{\text{M5s}}\mathbf{t})) = \mathbf{0},\tag{38}$$

where the boundary conditions are

$$\text{w}\_2(0, \mathbf{t}) = 0, \frac{\partial \Phi\_2}{\partial \mathbf{x}}(0, \mathbf{t}) = 0, \text{ at } \mathbf{x} = \mathbf{0},\tag{39}$$

$$\text{sw2}(\text{L}\_{\text{M}\text{s}}, \mathbf{t}) = 0, \frac{\partial \Phi\_2}{\partial \mathbf{x}}(\text{L}\_{\text{M}\text{s}}, \mathbf{t}) = 0, \text{ at } \mathbf{x} = \text{L}\_{\text{M}\text{s}}.\tag{40}$$

$$\mathbf{k}\mathbf{r}\_{\rm TT}\mathbf{\overset{\bullet}{\mathbf{w}}}\_{1} + \mathbf{k}\_{\rm TT}\mathbf{w}\_{1} - \mathbf{k}\mathbf{A}\mathbf{G}\left(\frac{\partial\mathbf{w}\_{2}}{\partial\mathbf{x}}(\mathbf{0},\mathbf{t}) - \Phi\_{2}(\mathbf{0},\mathbf{t})\right) = \mathbf{0}, \text{ at } \mathbf{x} = \mathbf{0},\tag{41}$$

$$\mathbf{c}\_{\text{IrB}}\mathbf{\overset{\bullet}{\mathbf{w}}\_{3}} + \mathbf{k}\_{\text{IrB}}\mathbf{w}\_{3} - \mathbf{k}\mathbf{A}\mathbf{G} \left(\frac{\partial \mathbf{w}\_{2}}{\partial \mathbf{x}} (\mathbf{L}\_{\text{M\star}\prime}\mathbf{t}) - \Phi\_{2}(\mathbf{L}\_{\text{M\star}\prime}\mathbf{t})\right) = \mathbf{0}, \text{ at } \mathbf{x} = \mathbf{L}\_{\text{M\star}\prime} \tag{42}$$

$$\mathrm{I}\_{\mathrm{Y}\_{2}\mathrm{Tr}}\frac{\partial^{2}\Phi\_{2}}{\partial t^{2}}(0,\mathrm{t})-\mathrm{c}\_{\mathrm{t}\mathrm{T}\mathrm{Tr}}\Big(\stackrel{\bullet}{\Phi\_{1}\cos\beta\_{1}-\frac{\partial\Phi\_{2}}{\partial\mathbf{t}}}(0,\mathrm{t})\Big)-\mathrm{k}\_{\mathrm{t}\mathrm{T}\mathrm{Tr}}(\Phi\_{1}\cos\beta\_{1}-\Phi\_{2}(0,\mathrm{t}))=0,\ \mathrm{at}\,\mathrm{x}=0,\tag{43}$$

$$\mathrm{I}\_{\mathrm{Y}\_{2}\mathrm{Ir}}\frac{\partial^{2}\Phi\_{2}}{\partial\mathbf{t}^{2}}(\mathrm{L}\_{\mathrm{M\varsigma}},\mathbf{t})-\mathrm{c}\_{\mathrm{IIrB}}\Big(\stackrel{\bullet}{\Phi\_{3}\cos\beta\_{2}}-\frac{\partial\Phi\_{2}}{\partial\mathbf{t}}(\mathrm{L}\_{\mathrm{M\varsigma}},\mathbf{t})\Big)-\mathrm{k}\_{\mathrm{I}\mathrm{Tr}}(\Phi\_{3}\cos\beta\_{2}-\Phi\_{2}(\mathrm{L}\_{\mathrm{M\varsigma}},\mathbf{t}))=0,\ \mathrm{at}\,\mathrm{x}=\mathrm{L}\_{\mathrm{M\varsigma}}\tag{44}$$

Taking into account that the tulip was fixed in the gearbox (see Figure 6 on the left side) and had simple elastic supports in the tulip-tripod joint (see Figure 6), the midshaft was a uniform Timoshenko beam simply supported by elastic at both ends (see Figures 6 and 7) in both the tulip-tripod (left side) and inner race-bowl joints (right side), while the bowl was fixed (see Figure 7, right side) in the steering wheel and had simple

elastic support in the inner race-bowl joint (see Figure 7). The system is given by the Equations (33)–(38) together with the boundary conditions from Equations (39)–(44) to represent the nonlinear dynamic behavior of the AD elements in FBV induced by shock force through the wheel by road excitations. By analyzing Equations (35) and (36), it can be remarked that they represent a system of equations with partial derivatives for the bending-shearing vibrations of a uniform shaft that considers the effects of rotary inertia and shear deformation, with the midshaft being a Timoshenko beam. The boundary conditions are given by Equations (39)–(44) and link the bending-shearing vibrations of the midshaft with the tulip and the bowl through the tulip-tripod and bowl-inner race joints, inducing in the solutions of the system of Equations (33)–(38) the next phenomena: the joints of the driveshaft are quasi-isometric [25–27] (p. 78), with the effect of geometric nonuniformity of the inertia characteristics of the joints that vary with the rigid angle of rotation for each element of the driveshaft in the directions X1, X2, X3, the effects of the bending deflection and bending-twisting stiffness for the tulip and the bowl, the effects of the bending deflection and bending-twisting damping for each joint of the driveshaft, the rotary inertia effect in bending, and the shearing effect for the midshaft. The novelty of the Hamilton's Principle approach was that it was used for the first time to compute the equations of motion for each AD element and the determination of the boundary conditions, as is noted in Equations (26)–(44).

#### **5. The Analytical Solutions**

The starting point to solve the system differential equations of the FBV movements (SDEOFBVM) (Equations (33)–(38)) for the AD elements (ADEs) was to analyze the vibration mechanism of the midshaft as a Timoshenko beam simply supported at the ends (see Figures 8 and 9). For the midshaft element of the AD, it was considered that f(x, t) = 0.

**Figure 8.** The part of the DMFFBV for the midshaft.

**Figure 9.** Cross-section of the tulip-tripod joint.

In analyzing the tulip-tripod-midshaft and bowl-inner race-midshaft joints, it became evident that the midshaft was simply supported at both ends where the concentrated mass of the tripod and the inner race mTr, mIr were placed, being a uniform shaft linked in bending-shearing vibration with the tulip for x = 0 and with the bowl forx=LMs. Therefore, the general solutions of Equations (35) and (36) that satisfy the boundary conditions of Equations (39) and (40) expressed in normalized bending deflection are [23] (pp. 326–328)

$$\mathbf{w}\_{2}(\overline{\mathbf{x}},\mathbf{t}) = \sin(\mathbf{n}\pi\overline{\mathbf{x}}) \left[ \mathbf{C}\_{1\text{n}}\cosh(\omega\_{\text{n}1}\mathbf{t}) + \mathbf{C}\_{2\text{n}}\sinh(\omega\_{\text{n}1}\mathbf{t}) + \mathbf{C}\_{3\text{n}}\cosh(\omega\_{\text{n}2}\mathbf{t}) + \mathbf{C}\_{4\text{n}}\sinh(\omega\_{\text{n}2}\mathbf{t}) \right],\\\overline{\mathbf{x}} = \frac{\mathbf{x}}{\mathbf{L}\_{\text{M\\$}}}, \mathbf{n} = 1, 2, 3, \dots \text{ (45)}$$

$$\begin{aligned} \Phi\_{2}(\overline{\mathbf{r}},t) &= \text{tr}\pi\cos(n\pi\overline{\mathbf{x}}) \left[ \left(1 + \frac{\mu\omega\_{\text{th}}^{2}}{4\mathcal{E}\left(\frac{\mu\omega}{\hbar\omega}\right)^{2}}\right) \mathbb{C}\_{\text{in}}\cosh(\omega\_{\text{th}}t) + \mathbb{C}\_{\text{2b}}\sinh(\omega\_{\text{th}}t) \right] + \left(1 + \frac{\mu\omega\_{\text{th}}^{2}}{\hbar\mathcal{E}\left(\frac{\mu\omega}{\hbar\omega}\right)^{2}}\right) \left[ \mathbb{C}\_{\text{2b}}\cosh(\omega\_{\text{th}}t) + \mathbb{C}\_{\text{4b}}\sinh(\omega\_{\text{th}}t) \right], \\ \mathbb{Z} &= \frac{\mathbf{x}}{\mathbf{1}\_{\text{th}}}, \mathbf{n} = 1,2,3,\ldots \end{aligned} \tag{46}$$

where ωn1, ωn2 are the natural frequencies of the Timoshenko beam simply supported at both ends [23] (p. 378), the smaller value corresponds to the bending deformation mode, and the bigger value corresponds to the shear deformation mode, given by the following equation:

$$w\_{\rm n1,2} = \sqrt{\frac{1 + \frac{\mathbf{n}^2 \pi^2 \mathbf{I}\_{\rm Y2\rm Mz}}{\Lambda \mathbf{I}\_{\rm Ms}^2} [1 + \mathbb{E}/\kappa \mathbf{c}] \mp \sqrt{\left[1 + \frac{\mathbf{n}^2 \pi^2 \mathbf{I}\_{\rm Y2\rm Mz}^2}{\Lambda \mathbf{I}\_{\rm Ms}^2} [1 + \mathbb{E}/\kappa \mathbf{c}]\right]^2 - 4 \frac{\mathbf{E}\_{\rm Y2\rm Mz}^2 \mathbf{n}^4 \pi^4}{\mathbf{k} \mathbf{A}^2 \mathbf{G} \cdot \mathbf{I}\_{\rm Ms}^4}}{2 \frac{\mathbf{P}\_{\rm Ms}^2 \mathbf{M}}{\mathbf{k} \mathbf{A} \mathbf{G}}}}}} \}} \,\,\mathbf{n} = 1, 2, 3, \dots \tag{47}$$

Additionally, the constants C1n, C2n, C3n, C4n are determined by the initial conditions at t = 0. In our case, the constants are given by the equation [4] (pp. 142–172)

$$\mathbf{C\_{1n}} = \mathbf{0}, \mathbf{C\_{2n}} \simeq \mathbf{0}. 1 \frac{\mathbf{F\_{S}} \Delta \mathbf{t\_{s}}}{\omega\_{\text{n}1} \mathbf{L\_{Ms}M}}, \mathbf{C\_{3n}} = \mathbf{0}, \mathbf{C\_{4n}} = \mathbf{0}, \tag{48}$$

where Δts is the time interval for shock excitations due to the road geometry and M is the car's mass. Injecting the solutions given by Equations (45) and (46) in the boundary conditions of Equations (41) and (42) yields the following boundary conditions for the tulip and the bowl in bending:

$$\mathbf{c}\_{\rm TTr} \overset{\bullet}{\mathbf{w}}\_{\rm I} + \mathbf{k}\_{\rm TTr} \mathbf{w}\_{\rm I} = -\frac{\rho \mathbf{A} \mathbf{L}\_{\rm M\*}^2}{\mathbf{n} \pi} \boldsymbol{\omega}\_{\rm n1}^2 \mathbf{C}\_{\rm 2n} \sinh(\omega \mathbf{u}\_{\rm n1} \mathbf{t}), \mathbf{n} = 1, 2, 3, \dots \tag{49}$$

$$\mathbf{c\_{IrB}w\_3 + k\_{IrB}w\_3 = (-1)^{\text{n}} \frac{\rho A L\_{\text{Ms}}^2}{\mathbf{n}\pi} w\_{\text{n1}}^2 \mathbf{C\_{2n}} \text{sinch}(\omega\_{\text{n1}} \mathbf{t}), \mathbf{n} = 1, 2, 3, \dots \tag{50}$$

Equations (49) and (50) can now be injected in Equations (33) and (37) for the DMFFBV, thus yielding DMFFBV equations for the tulip and bowl:

$$\mathbf{m}\_{\rm T} \overset{\bullet \bullet}{\mathbf{w}\_{1}^{\bullet}} + \mathbf{c}\_{1} \overset{\bullet}{\mathbf{w}\_{1}^{\bullet}} + + \mathbf{k}\_{1} \mathbf{w}\_{1} = -\frac{\rho \mathbf{A} \mathbf{L}\_{\rm M\rm s}^{2}}{\mathbf{n} \pi} \boldsymbol{\omega}\_{\rm n1}^{2} \mathbf{C}\_{\rm 2n} \sinh(\omega\_{\rm n1} \mathbf{t}), \mathbf{n} = 1, 2, 3, \dots \tag{51}$$

$$\mathbf{m\_{B}}\overset{\bullet\bullet}{\mathbf{w}}\_{3} + \mathbf{c\_{2}}\overset{\bullet}{\mathbf{w}}\_{3} + \mathbf{+} \mathbf{k\_{2}}\mathbf{w}\_{3} = (-1)^{\mathbf{n}} \frac{\rho \mathbf{A} \mathbf{L}\_{\mathbf{Ms}}^{2}}{\mathbf{n} \boldsymbol{\pi}} \omega\_{\mathbf{n}1}^{2} \mathbf{C}\_{\mathbf{2n}} \sinh(\omega\_{\mathbf{n}1} \mathbf{t}), \mathbf{n} = 1, 2, 3, \dots \tag{52}$$

Taking into account that Browne and Palazzolo in [6] mentioned that excitation terms, as presented on the right side of Equations (51) and (52) could be expressed as cubic internal tuned dampers (ITDs), the forms of Equations (51) and (52) can be modified by dividing Equations (51) and (52) into the mass of the tulip mT and the mass of the bowl mB, respectively, expressing the stiffness and damping coefficients as functions of the geometrical characteristics given by Equations (19) and (23) coupled with the expressions of inertial characteristics provided by Equations (8), (9), (17), and (18). Equations (51) and (52) become the normalized differential equations in the time functions w1(t), w3(t):

$$\mathbf{w}\_1^\bullet + 2\boldsymbol{\xi}\boldsymbol{\Omega}\_1 \sqrt{\frac{1-\mathsf{C}\_1\cos(2\varphi\_1)}{1-\mathsf{C}\_2\cos(2\varphi\_1)}} \mathbf{w}\_1^\bullet + \boldsymbol{\Omega}\_1^2 \frac{1-\mathsf{C}\_1\cos(2\varphi\_1)}{1-\mathsf{C}\_2\cos(2\varphi\_1)} \mathbf{w}\_1 = -\boldsymbol{\Gamma}\_0 \mathbf{w}\_1 - \boldsymbol{\Gamma}\_1 \mathbf{w}\_1^\mathbf{\overline{}} - \boldsymbol{\Gamma}\_2 \mathbf{w}\_1^\mathbf{\overline{}} \tag{53}$$

$$\overset{\bullet \bullet}{\mathbf{w}\_3} + 2\underline{\epsilon}\Omega\_3\sqrt{\frac{1-\mathcal{C}\_3\cos(2\varphi\_3)}{1-\mathcal{C}\_4\cos(2\varphi\_3)}}\overset{\bullet}{\mathbf{w}\_3} + \Omega\_3^2\frac{1-\mathcal{C}\_3\cos(2\varphi\_3)}{1-\mathcal{C}\_4\cos(2\varphi\_3)}\mathbf{w}\_3 = -\Gamma\_5\mathbf{w}\_3 - \Gamma\_3\mathbf{w}\_3^3 - \Gamma\_4\mathbf{w}\_3^5 \tag{54}$$

where the coefficients C1, C2, C3, C4 are presented in Appendix A, and the damping ratio ξ is given by the following equation:

$$\mathcal{L} = \frac{\Delta}{\sqrt{4\pi^2 + \Delta^2}},\tag{55}$$

The natural frequency of the tulip in bending Ω<sup>1</sup> and the natural frequency of the bowl in bending Ω<sup>3</sup> are given by the following equations:

$$\Omega\_{1} = \sqrt{\frac{3\mathcal{E}}{\mathrm{m}\_{\mathrm{T}}\mathrm{L}\_{\mathrm{TA}}^{3}} \frac{\mathbf{b}\_{1}(1+\mathbf{a}\_{1})}{1+\frac{\mathbf{b}\_{1}}{\mathrm{L}\_{2}\mathrm{TA}}(1+\mathbf{a}\_{1})}}, \mathbf{a}\_{1} = \frac{0.5(\mathrm{J}\_{1\mathrm{T}}+\mathrm{J}\_{2\mathrm{T}})}{\mathrm{S}\_{\mathrm{T}}\left[\frac{\mathrm{L}\_{\mathrm{TA}}^{2}}{12}+\mathrm{d}\_{\mathrm{CT}}^{2}\right]}, \mathbf{b}\_{1} = \mathrm{S}\_{\mathrm{T}}\left[\frac{\mathrm{L}\_{\mathrm{TA}}^{2}}{12}+\mathrm{d}\_{\mathrm{CT}}^{2}\right], \tag{56}$$

$$\Omega\_3 = \sqrt{\frac{3\mathcal{E}}{\mathrm{mg}\,\mathrm{L}\_{\mathrm{BA}}^3} \frac{\mathrm{b}\_2(1+\mathrm{a}\_2)}{1+\frac{\mathrm{b}\_2}{\mathrm{L}\_{2\mathrm{BA}}}(1+\mathrm{a}\_2)}}, \mathrm{a}\_2 = \frac{0.5(\mathrm{J}\_{1\mathrm{B}}+\mathrm{J}\_{2\mathrm{B}})}{\mathrm{S}\_{\mathrm{B}}\left[\frac{\mathrm{L}\_{\mathrm{BA}}^2}{\mathrm{L}\_2}+\mathrm{d}\_{\mathrm{CB}}^2\right]}, \mathrm{b}\_2 = \mathrm{S}\_{\mathrm{B}}\left[\frac{\mathrm{L}\_{\mathrm{BA}}^2}{12}+\mathrm{d}\_{\mathrm{CB}}^2\right]. \tag{57}$$

The coefficients of the cubic and quintic terms expressed in Equations (53) and (54) can be mathematically explained by development in a Taylor infinite series induced by the excitation function sinh(x) (see Equations (53) and (54)) expressed approximately as

$$\sinh(\chi) \approx \chi + \frac{\chi^3}{3!} + \frac{\chi^5}{5!}. \tag{58}$$

In the PPR region, the excitation terms cos(2ϕ1(t)), sinh(ωn1t), and cos(2ϕ3(t)) must satisfy the Equations (59) and (60), as mentioned in [28] (p. 425), since this resonance is one of the most important, as mentioned in [5,6,28]:

$$
\sigma\_1 \simeq 2\Omega\_1 \simeq \frac{\text{d}\theta\_1}{\text{d}\mathfrak{t}} = \frac{2\,\text{d}\varphi\_1}{\text{d}\mathfrak{t}} \simeq \omega\_{\text{n1}}.\tag{59}
$$

$$
\eta\_3 \simeq 2\Omega\_3 \simeq \frac{\mathrm{d}\theta\_3}{\mathrm{d}\mathfrak{t}} = \frac{2\mathrm{d}\varphi\_3}{\mathrm{d}\mathfrak{t}} \simeq \omega\_{\mathrm{n1}}.\tag{60}
$$

where η<sup>1</sup> is the excitation frequency in the PPR region for the tulip, η<sup>3</sup> is the excitation frequency in the PPR region for the bowl, <sup>d</sup>ϕ<sup>1</sup> dt is the excitation induced in the FBV movement (Equation (53)) by the rigid twisting angle of the tulip concerning the X1 axis, and <sup>d</sup>ϕ<sup>3</sup> dt is the excitation induced in the FBV movement (Equation (54)) by the rigid twisting angle of the bowl concerning the X3 axis, as given by the following equation:

$$\log y(\mathbf{t}) = \varphi\_1(\mathbf{t}) + \frac{\mathsf{R}\_{\mathrm{Tlr}}}{2\mathrm{L}\_{\mathrm{Ma}}} \tan \beta\_1 \tan^2 \frac{\mathfrak{F}\_1}{2} \cos(3\varphi\_1) + \frac{\mathsf{R}\_{\mathrm{bft}}}{2\mathrm{L}\_{\mathrm{Mb}}} \tan \beta\_2 \tan^2 \frac{\mathfrak{F}\_2}{2} \cos \left[ 3 \left( \varphi\_1(\mathbf{t}) + \frac{\mathsf{R}\_{\mathrm{Tlr}}}{2\mathrm{L}\_{\mathrm{Ma}}} \tan \beta\_1 \tan^2 \frac{\mathfrak{F}\_1}{2} \cos(3\varphi\_1) \right) \right],\tag{61}$$

This considers the nonuniformity of the kinematic isometry of the AD [25,26]. In Equation (61), the terms RTTr,RIrB represent the radii of the tulip-tripod and inner racebowl joints. Considering the conditions imposed by Equations (59) and (60) together with Equation (58) and Equations (48)–(50), the coefficients of the cubic and quintic excitation terms on the right-hand side of Equations (53) and (54) are

$$
\Gamma\_1 = \frac{0.2\rho\mathcal{A}}{\pi\mathfrak{d}!} \frac{\overline{\mathcal{F}}\_{\mathsf{S}}\mathsf{A}\mathsf{t}\_{\mathsf{S}}}{\mathrm{M}\mathsf{m}\_{\mathsf{T}}} \Omega\_{1\prime}\Gamma\_2 = \frac{0.2\rho\mathcal{A}}{\pi\mathfrak{d}!} \frac{\overline{\mathcal{F}}\_{\mathsf{S}}\mathsf{A}\mathsf{t}\_{\mathsf{S}}}{\mathrm{M}\mathsf{m}\_{\mathsf{T}}} \Omega\_{1\prime} \frac{0.2\rho\mathcal{A}}{\pi} \frac{\overline{\mathcal{F}}\_{\mathsf{S}}\mathsf{A}\mathsf{t}\_{\mathsf{S}}}{\Omega\_{1}\mathrm{M}\mathsf{m}\_{\mathsf{T}}} \ll 1 \Rightarrow \Gamma\_0 \approx 0,\tag{62}
$$

$$
\Gamma\_3 = \frac{0.2\rho\text{A}}{\pi \text{3!}} \frac{\text{F}\_\text{S}\text{At}\_\text{s}}{\text{M}\text{m}\_\text{B}} \Omega\_\text{3},
\Gamma\_4 = \frac{0.2\rho\text{A}}{\pi \text{5!}} \frac{\text{F}\_\text{S}\text{At}\_\text{s}}{\text{M}\text{m}\_\text{B}} \Omega\_\text{3},
\frac{0.2\rho\text{A}}{\pi} \frac{\text{F}\_\text{S}\text{At}\_\text{s}}{\Omega\_\text{3}\text{M}\text{m}\_\text{B}} \ll 1 \Rightarrow \Gamma\_5 \approx 0.\tag{63}
$$

By analyzing Equations (53) and (54), they are a generalized form of Mathieu–Hill nonlinear equations coupled through Equation (61), describing the FBV movements of the tulip and the bowl induced by shock excitations due to the road geometry. These equations contain the following phenomena: the non-uniformity of geometric and kinematic isometry of the AD, the non-uniformity of the inertial characteristics for the tulip and the bowl and forced shock excitations due to road excitations. By analyzing the initial system of Equations (33)–(38) with the boundary conditions of Equations (39)–(44), it can be seen that Equations (34) and (38), together with the boundary conditions of Equations (43) and (44), were not used in the mathematical procedure solution because these equations represented the nonlinear dynamic behavior of the FBV angular shearing movements of the AD (the AD's elements of the tulip, midshaft, and bowl), and the present article is dealing with the nonlinear dynamic behavior of the FBV deflection movements of the ADE, namely the bowl, midshaft, and tulip. The FBV movements of the midshaft have the solution given by Equation (45) coupled with Equations (47) and (48). In contrast, the shearing FBV movements of the midshaft have the solution given by Equation (46) coupled with Equations (47) and (48), both representing the solution of Equations (35) and (36) from the SDEOFBVM (Equations (33)–(38)), with the boundary conditions of Equations (39) and (40) from the general system of Equations (39)–(44) that imposed the boundary conditions. Additionally, it is noted that for the nonlinear dynamic behavior of the AD, the dynamic behavior of the midshaft represents an excitation-like internal tuned damper (ITD) used by Browne and Palazzolo in the superharmonic nonlinear lateral vibrations (forced bending vibrations) of an AD [6]. This aspect represents a novelty in the investigation of the dynamic behavior of the ADE, because it allows for the computation of the equations of motion for the tulip and the bowl in a mathematical form adequate for performing a perturbation approach, as can be seen from Equations (53) and (54). Another aspect of this novelty is the capability of computing the natural free frequency in bending for the tulip and the bowl (see Equations (56) and (57)) and to express the excitation frequency in the PPR region, as is noted from Equations (59) and (60). In the meantime, a supplementary aspect of the novelty is the estimation of the shock excitation of the road geometry as cubic and quintic excitation terms in the equations of motion, as is shown by Equations (56) and (57) and the relations of Equations (62) and (63). For the first time, the effect of the nonuniformity of the kinematic isometry of the AD was considered, being expressed by Equation (61), which links the equations of motion for the FBVM regarding the tulip and the bowl.

#### **6. Analysis of the PPR Region for the Tulip and Bowl**

As already mentioned above, one of the most important resonant cases of an AD is the PPR [5,6], [28] (p. 425), and for this paper, the authors decided to investigate the dynamic behavior of FBV movement in the PPR region for the tulip and the bowl based on Equations (53) and (54). Similar experimental data for this case study are presented in the literature by Steinwede [4] (pp. 69–144). The study considered a tulip having geometric characteristics such as the geometric moment of inertia and nonuniformity of the geometric moments of inertia, as presented in Table 1 [29], for the tulip AD element. By comparing these presented geometric characteristics with those considered by Steinwede [4] (p. 111), it can be concluded that there is agreement. Using AUTOCAD software, J1T, J2T were computed based on the direct geometric characteristics LTB, LTA, LMs, dTA, dCT, ST and the general geometry of the tulip.

**Table 1.** Geometry characteristics of the tulip [29].


It can be remarked that the nonuniformity of the inertial geometric characteristic *χ*nT in Table 1 had another value than what was considered in literature because of the geometric characteristic in the cross-section of the tulip presented in Figure 9.

Table 2 presents the physical properties of the material of the tulip and the amplitude and duration of the shock, considering that the material was steel and iron cast. By comparing these presented material properties with those considered by Steinwede [4] (p. 112), it can be concluded that they were in very close agreement.

**Table 2.** Material properties of the tulip and shock amplitude.


To compute the amplitude of the FBV movements of the ADE tulip in the PPR region using a perturbation approach, several methods could be used: the method of harmonic balance [29] (p. 66), the asymptotic method [1] (pp. 299–393), or the method of multiple scales [28] (pp. 424–427). The method of harmonic balance is very efficient, as mentioned in the literature [30] (p. 63), [31], but it allows the computation of the amplitude versus the frequency excitation only for stationary (steady state) motion without the possibility of computing the phase angle [30] (pp. 63–68). Therefore, this method is limited to the steady state FBV. The multiple scales method is challenging to use due to the conditions of zeroing the secular term and the additional necessary study of the convergence of the detuning parameter. Considering these aspects, the authors chose to use the AMA in the first-order approximation. This would allow the investigation of the amplitude and phase angle versus the excitation frequency for the FBV movement of the tulip in the PPR region for both the stationary and nonstationary cases. To compute the solution of Equation (53), it was assumed that the slowing time was τ = εt, where ε is a small positive parameter [1] (p. 299). To introduce the slowing time, Equation (53) needed to be transformed to be used

in the AMA. The coefficients of the second and third term of Equation (53) on the left side can be expressed as

$$\begin{cases} \frac{1-\mathcal{C}\_{1}\cos(2\varphi\_{1})}{1-\mathcal{C}\_{2}\cos(2\varphi\_{1})} \approx 1-2\mu\cos(2\varphi\_{1}) = 1-2\mu\cos(\theta\_{1}),\\ \sqrt{\frac{1-\mathcal{C}\_{1}\cos(2\varphi\_{1})}{1-\mathcal{C}\_{2}\cos(2\varphi\_{1})}} \approx \left(1-2\mu\cos(2\varphi\_{1})\right)^{-\frac{1}{2}} \approx 1+\mu\cos(2\varphi\_{1}) = 1+\mu\cos(\theta\_{1}),\\ \mu = \frac{1}{2}(\mathcal{C}\_{2}-\mathcal{C}\_{1}),\end{cases} \tag{64}$$

where in Equation (64), the condition given by Equation (59) was considered. With the relations in Equation (64), Equation (53) becomes

$$\mathbf{w}\_1^\bullet \mathbf{w}\_1^\bullet + \Omega\_1^2 \mathbf{w}\_1 = -2\xi \Omega\_1 (1 + \mu \cos(\theta\_1)) \mathbf{w}\_1^\bullet + 2\mu \Omega\_1^2 \cos(\theta\_1) \mathbf{w}\_1 - \Gamma\_1 \mathbf{w}\_1^3 - \Gamma\_2 \mathbf{w}\_1^5. \tag{65}$$

The assumption that the damping ratio ξ, the excitation coefficient μ, and the coefficients of cubic and quintic nonlinearity Γ1, Γ<sup>2</sup> are small is incorporated in the analysis by representing these quantities in the following form:

$$
\xi = \varepsilon \xi, \mu = \varepsilon \mu, \Gamma\_1 = \varepsilon \Gamma\_1, \Gamma\_2 = \varepsilon \Gamma\_2. \tag{66}
$$

where ε is the same small positive parameter used to obtain the slowing time. It is also assumed that the excitation frequency η<sup>1</sup> and the excitation parameter μ vary slowly with time, such that

$$\frac{d\theta\_1}{dt} = \eta\_1(\tau), \mu = \mu(\tau). \tag{67}$$

Equation (65) becomes the following after neglecting the terms in ε2:

$$\mathbf{W}\_{1}^{\bullet \bullet} + \Omega\_{1}^{2} \mathbf{w}\_{1} = \varepsilon \left[ -2\xi\Omega\_{1}\mathbf{w}\_{1}^{\bullet} + 2\mu\Omega\_{1}^{2}\cos(\theta\_{1})\mathbf{w}\_{1} - \Gamma\_{1}\mathbf{w}\_{1}^{3} - \Gamma\_{2}\mathbf{w}\_{1}^{5} \right]. \tag{68}$$

Regarding Equation (68), the right-hand side represents a perturbation of the mathematical form H(w1, θ1), being a periodic function in θ<sup>1</sup> with period 2π, while the left-hand side of the equation is a linear oscillator. By considering all these physical considerations and confining our attention to the investigation of the PPR region, a solution for Equation (68) is sought after in the following form to the first-order approximation in ε:

$$\mathbf{w}\_1 = \mathbf{W}\cos\left(\frac{1}{2}\theta\_1 + \Psi\right),\tag{69}$$

where W and Ψ are functions of time defined by the system of differential equations:

$$\begin{cases} \begin{array}{c} \frac{\text{d}\text{W}}{\text{d}\text{t}} = \varepsilon \text{A}\_{1}(\text{\textpi}, \text{W}, \text{\textpi})\\ \frac{\text{d}\text{t}}{\text{d}\text{t}} = \Omega\_{1} - \frac{1}{2}\eta\_{1}(\text{\textpi}) + \varepsilon \text{B}\_{1}(\text{\textpi}, \text{W}, \text{\textpi}) \end{array} . \end{cases} \tag{70}$$

Differentiating the right-hand side of Equation (70) and expanding the results in powers of ε yields

$$\begin{aligned} \stackrel{\bullet \bullet}{W} &= \frac{d^2 \Psi}{dt^2} = \varepsilon \frac{\partial \Lambda\_1}{\partial \Psi} \left(\Omega\_1 - \frac{1}{2} \eta\_1\right) + \varepsilon^2 \dots \\\stackrel{\bullet \bullet}{\Psi} &= \frac{d^2 \Psi}{dt^2} = \varepsilon \left[ -\frac{1}{2} \frac{\partial \eta\_1}{\partial \tau} + \frac{\partial \mathcal{B}\_1}{\partial \Psi} \left(\Omega\_1 - \frac{1}{2} \eta\_1\right) \right] + \varepsilon^2 \dots \end{aligned} \tag{71}$$

When differentiating the right-hand side of Equation (69) while considering the systems of Equations (70) and (71) and retaining only the first-order terms, we obtained the following:

$$\mathbf{w}\_{1}^{\bullet} = -\mathbf{W}\Omega\_{1}\sin\left(\frac{1}{2}\theta\_{1} + \psi\right) + \varepsilon \left[\mathbf{A}\_{1}\cos\left(\frac{1}{2}\theta\_{1} + \psi\right) - \mathbf{W}\mathbf{B}\_{1}\sin\left(\frac{1}{2}\theta\_{1} + \psi\right)\right],\tag{72}$$

$$\mathbf{W}\_{1}^{\bullet \bullet} = -\mathsf{W}\boldsymbol{\Omega}\_{1} \cos\left(\frac{1}{2}\boldsymbol{\theta}\_{1} + \boldsymbol{\psi}\right) + \varepsilon \left[\frac{\partial\boldsymbol{\Lambda}\_{1}}{\partial\boldsymbol{\psi}}\left(\boldsymbol{\Omega}\_{1} - \frac{1}{2}\boldsymbol{\eta}\_{1}\right) - 2\mathsf{W}\boldsymbol{\Omega}\_{1}\boldsymbol{\Omega}\_{1}\right] \cos\left(\frac{1}{2}\boldsymbol{\theta}\_{1} + \boldsymbol{\psi}\right) - \varepsilon \left[\mathsf{W}\frac{\partial\mathsf{B}\_{1}}{\partial\boldsymbol{\psi}}\left(\boldsymbol{\Omega}\_{1} - \frac{1}{2}\boldsymbol{\eta}\_{1}\right) + 2\mathsf{D}\_{1}\mathbf{A}\_{1}\right] \sin\left(\frac{1}{2}\boldsymbol{\theta}\_{1} + \boldsymbol{\psi}\right). \tag{73}$$

Using Equations (72) and (73) in the basic form of Equation (68), equating the terms of the form ε cos<sup>1</sup> <sup>2</sup>θ<sup>1</sup> + ψ , ε sin<sup>1</sup> <sup>2</sup>θ<sup>1</sup> + ψ from the left-hand side of Equation (68) with the same terms from the right-hand side of the same equation and neglecting the overtones (the terms having cos <sup>n</sup> <sup>2</sup> θ<sup>1</sup> + ψ , sin <sup>n</sup> <sup>2</sup> θ<sup>1</sup> + ψ , n = 2, 3, 4, 5) yields the two following coupled first-order differential equations for the unknown functions A1 and B1:

$$\frac{\partial \mathcal{A}\_1}{\partial \boldsymbol{\Psi}} \left( \boldsymbol{\Omega}\_1 - \frac{1}{2} \boldsymbol{\eta}\_1 \right) - 2 \boldsymbol{\mathcal{W}} \boldsymbol{\Omega}\_1 \mathcal{B}\_1 = \mu \boldsymbol{\mathcal{W}} \boldsymbol{\Omega}\_1^2 \cos(2\boldsymbol{\psi}) - \frac{3}{4} \boldsymbol{\Gamma}\_1 \boldsymbol{\mathcal{W}}^3 - \frac{5}{8} \boldsymbol{\Gamma}\_2 \boldsymbol{\mathcal{W}}^5,\tag{74}$$

$$2\frac{\partial \mathcal{B}\_1}{\partial \psi} \left(\Omega\_1 - \frac{1}{2}\eta\_1\right) + 2\Omega\_1 \mathcal{A}\_1 = -\mu \mathcal{W} \Omega\_1^2 \sin(2\psi) - 2\xi \Omega\_1^2 \mathcal{W},\tag{75}$$

This system of Equations (74) and (75) has the following solutions:

$$\begin{cases} \mathbf{A}\_{1} = -\frac{\mu \Omega\_{1}^{2} \mathbf{W}}{\eta\_{1}} \sin 2\boldsymbol{\upmu} - \xi \Omega\_{1} \mathbf{W}\_{\prime} \\\ \mathbf{B}\_{1} = -\frac{\mu \Omega\_{1}^{2}}{\eta\_{1}} \cos 2\boldsymbol{\upmu} + \frac{3}{8} \frac{\Gamma\_{1} \mathbf{W}^{2}}{\Omega\_{1}} + \frac{5}{16} \frac{\Gamma\_{2} \mathbf{W}^{4}}{\Omega\_{1}} \end{cases} . \tag{76}$$

By introducing the solutions of Equation (76) into the system of Equation (70) and transforming all the system parameters back to their real-time values, the solution of Equation (65), representing the FBV movement of the tulip in the PPR (η<sup>1</sup> 2Ω1) is given by Equation (69) using the first-order approximation of the AMA, where the amplitude W and the phase angle ψ are given by integrating over time the first-order differential system:

$$\begin{cases} \begin{array}{c} \frac{\text{dW}}{\text{dt}} = -\frac{\mu \Omega\_1^2 \text{W}}{\eta\_1} \sin 2\psi - \xi \Omega\_1 \text{W} \\\ \frac{\text{d}\psi}{\text{dt}} = \Omega\_1 - \frac{1}{2}\eta\_1 - \frac{\mu \Omega\_1^2}{\eta\_1} \cos 2\psi + \frac{3}{8} \frac{\Gamma\_1 \text{W}^2}{\Gamma\_1} + \frac{5}{16} \frac{\Gamma\_2 \text{W}^4}{\Gamma\_1} \end{array} . \end{cases} \tag{77}$$

#### *6.1. Steady State Forced Bending Vibration of the Tulip*

For the stationary FBV movement for the tulip of the AD, it would be set to zero dW dt and <sup>d</sup><sup>ψ</sup> dt in Equation (77), and expressed from these equations, sin 2ψ and cos 2ψ yields the trigonometric equation sin2(2ψ) + cos2(2ψ) = 1, expressing in terms of the system in Equation (77) the bi-quartic equation in the unknown W the amplitude versus the excitation frequency:

$$
\lambda\_5 \mathcal{W}^8 + \lambda\_6 \mathcal{W}^6 + \lambda\_7 \mathcal{W}^4 + \lambda\_8 \mathcal{W}^2 + \lambda\_9 = 0,\tag{78}
$$

where the coefficients λ5–λ<sup>9</sup> are given in Appendix A. The solutions to Equation (78) represent the stationary amplitude of the forced bending vibration of the tulip in the PPR region. The stationary phase angle of the FBV movement for the tulip in the PPR region is given from the same system (Equation (77)) set to zero by expressing the tangent of the phase angle, yielding

$$\Psi = \frac{1}{2} \text{arctg} \left[ -\frac{\xi\_1}{1 - \frac{\Omega\_1}{2\Omega\_1} + \frac{3}{8}\frac{\Gamma\_1 W^2}{\Omega\_1^2} + \frac{5}{16}\frac{\Gamma\_2 W^4}{\Omega\_1^2}} \right]. \tag{79}$$

#### *6.2. Investigation of Dynamic Instability*

The investigation of the dynamic instability of the FBV movement for the tulip represents the computation of the boundaries of the principal parametric region of instability. The

base width of the stationary response is the only region in which vibrations may typically initiate. Setting to zero the amplitude in Equation (78) yields the following equation:

$$\left(1 - \frac{\eta\_1}{2\Omega\_1}\right)^2 + \xi^2 - \frac{1}{4}\mu^2 \left(\frac{2\Omega\_1}{\eta\_1}\right)^2 = 0,\tag{80}$$

Alternatively, it may yield

$$2\left(\frac{\eta\_1}{2\Omega\_1}\right)^4 - 2\left(\frac{\eta\_1}{2\Omega\_1}\right)^3 + \left(1 + \xi^2\right)\left(\frac{\eta\_1}{2\Omega\_1}\right)^2 - \frac{\mu^2}{4} = 0. \tag{81}$$

Equation (81) gives the boundaries of the principal parametric region of instability in the space <sup>η</sup><sup>1</sup> 2Ω<sup>1</sup> , ξ, μ . This is the investigation for the stationary FBV movement of the tulip. For the nonstationary FBV movement of the tulip, the investigation of the dynamic instability consisted of the analysis of the spectral graphs of the nonstationary amplitude and the first derivative of the nonstationary amplitude concerning time, given by the integration of the system of Equation (77). The graphed representation in the configuration space of the "speed" nonstationary amplitude versus the nonstationary amplitude dW dt , W is evidence of the transition of FBV movement for the tulip through the chaotic movement in the PPR region, and the results are presented in the next paragraph.

#### *6.3. Nonstationary FBVM of the Tulip and Bowl during Transition through the PPR*

The resonant regimes of the principal parametric vibration for the forced bending vibration of the tulip were determined. It may be of interest to examine the nature of the nonstationary vibrations performed by the system during a transition of the excitation frequency η<sup>1</sup> through the resonant regime. For this purpose, it was necessary to integrate the system of Equation (77) governing the nonstationary amplitude and the phase angle. This system of differential equations cannot be integrated into the closed form, and therefore, it was obvious to use numerical integration. It may be pointed that numerical integration of first-order equations governing the stationary amplitude and phase angle versus frequency excitation is a much simpler process than the integration of the original system of Equation (77), which is nonlinear since it must evaluate the envelope of an oscillatory function and not the function itself. In the present nonstationary analysis, the sweep of the excitation frequency is taken to be logarithmic and is expressed as

$$
\mathfrak{u}\_1(\mathfrak{t}) = \mathfrak{u}\_{10} \sigma^{\mathfrak{t}(1+\mathfrak{m})}, \sigma \ge 1,\tag{82}
$$

where η<sup>10</sup> is the initial frequency at t = 0 for the tulip and m is the sweep rate (in octaves/min). The differential system of Equation (77) is integrated numerically using a fourth-to-fifth order Runge–Kutta algorithm from the MATLAB software with a variable integration step of the prediction-correction type. The initial values of the variables W, ψ, and η<sup>1</sup> were chosen as a stationary state starting point. The investigation of the nonstationary forced bending vibration of the tulip was confined to the PPR region, and the results are presented in the next paragraph. For the bowl, the analysis of the PPR region was based on Equation (54) for FBV movement (FBVM), with the excitation frequency η<sup>3</sup> given by Equation (60), the natural frequency Ω<sup>3</sup> established by Equation (57), the cubic and quintic terms Γ3, Γ<sup>4</sup> given by Equation (63), and the excitation coefficient being μ = <sup>1</sup> <sup>2</sup> (C4 − C3), while the terms C3, C4 are presented in Appendix A. For the investigation of the nonstationary FBVM of the bowl, the system of differential equations is from Equation (77), using the excitation frequency η3, the natural frequency Ω3, and the cubic and quintic terms Γ3, Γ4, while the rate of sweep m is established as

$$
\mathfrak{u}^{\mathfrak{z}}(\mathfrak{t}) = \mathfrak{u}^{\mathfrak{z}0} \mathfrak{a}^{\mathfrak{t}(1+\mathfrak{m})}, \mathfrak{a} \ge 1,\tag{83}
$$

where η<sup>30</sup> is the initial frequency at t = 0 for the bowl. Table 3 presents the geometry and inertial characteristics of the bowl. Using AUTOCAD software, J1B,J2B were computed based on the direct geometric characteristics LBB, LBA, LMs, dBA, dCB, SB and the general geometry of the tulip. The material's physical properties are the same as those presented in Table 2.

**Table 3.** Geometry characteristics of the bowl.


**Figure 10.** Axonometric cross-section of the bowl-balls-inner race joint.

When comparing these presented geometry characteristics (see Table 3) with those considered by Steinwede [4] (p. 112), it can be concluded that they agree. The novelty for investigating the PPR region of an AD is that the AMA was used for such a multibody technical system for the first time. The AMA has no limitations on its use because it allows the development of analytical solutions for the amplitude, phase angle, and dynamic instability frontiers for stationary and nonstationary motion, as can be seen from Equations (65)–(83). Additionally, the AMA allows for significant precision for the computations because it can provide several levels of approximation. In our case, it was used for firstorder approximation. For the first time, the amplitude, phase angle, and dynamic instability frontiers for the stationary motion of the FBV of a tulip and bowl in the PPR region were computed, as can be seen by analyzing Equations (78), (79) and (81). In the meantime, the amplitude spectrum, "speed" amplitude spectrum, and representation in the phase space W, dW dt of the "speed" amplitude versus amplitude for the nonstationary motion in transition through the PPR region of the tulip and the bowl was determined for the first-time using Equations (77), (82) and (83).

## **7. Results and Discussion**

The first results obtained were represented by the natural free frequencies in bending and shearing of the midshaft, given by Equation (56):

$$\mathbf{v}\_{\mathbf{n}\_{1,2}} = \frac{\omega\_{\mathbf{n}1,2}}{2\pi}, \mathbf{n}\_{1,2} = 1, 2, 3, \dots \dots \tag{84}$$

Based on Equations (56) and (83), the natural free frequencies of the midshaft in bending-shearing is an infinite series covering a wide range of frequencies, and thus by comparing the numerical results obtained with the mentioned equations and the experimental data and numerical data (using FEA-ANSYS) presented in the literature [13] and analyzing the data presented in Table 4, agreement for the bending natural free frequencies of the midshaft of the AD was found.


**Table 4.** Comparison of the experimental natural frequencies (ENFs) of the driveshaft and numerical simulations using the finite element analysis (FEA) program Ansys of the natural frequencies (NSNFs) in [13] with those obtained for the midshaft based on Equations (57) and (83).

As can be seen from Table 4, the agreement of the results was for the free bending natural frequency, while the second free natural frequency represents the free shearing natural frequency that, as can be noticed, was like the principal parametric resonance compared with the bending. This confirms that the designed DMFFBVM of the AD was much more reliable than the FEA-ANSYS and modal analysis presented in the literature [13]. As mentioned in [13], the dynamic forced vibration of the AD was inducing cracks in the AD's structure, a fact also revealed by Steinwede [4]. Additionally, it was mentioned in [13] that experimental data were found for the midshaft, and this aspect revealed the fact that the midshaft had such a nonlinear dynamic behavior that represented the excitation element for the other components of the AD, namely the tulip and the bowl, playing the same role as the internal tuned damper (ITD) mentioned in the literature [8] for the superharmonic nonlinear lateral vibrations (bending vibrations) of the AD. The fundamental natural free bending vibrations of the tulip given by Equation (56) was ν<sup>1</sup> = 519.19 Hz, considering the geometrical characteristic and the physical properties presented in Tables 1 and 2 for the AD designed for a heavy duty SUV, which conformed with the experiments of Steinwede that imposed an investigation of the forced torsional or bending vibrations of a driveshaft in the range of 0.5–15 kHz [6] (p. 119). Unfortunately, no published experiments investigated the natural free frequency in bending only for the global tulip, midshaft, or global bowl. For the steady state FBV of the tulip and the bowl, specific MATLAB software was developed to compute the amplitude and phase angle in the region of principal parametric resonance based on Equations (78) and (79). Figures 11 and 12 present the variation of the normalized stationary amplitude for the FBV in the PPR region for the tulip of an AD, with this being around 1038.38 Hz.

When analyzing the normalized stationary parametric amplitude presented in the graphs in Figures 11 and 12 for ζ = 0.0016–0.0318 and *χ*nT = 0.25, it can be found that for the cases, we had a manifestation of a "soft spring" with one branch for η<sup>1</sup> ≤ 2Ω<sup>1</sup> that indicated the presence of interaction between the principal parametric resonance and the primary resonance, while for η<sup>1</sup> ≥ 2Ω1, two branches of a "hard spring" existed, which suggests the manifestation of interaction with the combination resonance or with the internal resonance for the tulip [28] (pp. 132–160). This aspect agrees with the experimental data in the literature [4] (pp. 130–144). It can also be remarked that the increase in the damping ratio in the range of 0.0016–0.0318 increased the minimum normalized stationary amplitude from 0.13 to 0.42 for one branch, and for the second branch, the increase in the minimum normalized stationary amplitude was from 0.32 to 0.425. This can be explained by the existence of dynamic instability in the region of principal parametric resonance, as mentioned in the literature [28] (pp. 132–160). Figures 13 and 14 present the variation of the stationary phase angle for the FBV in the PPR region for the tulip of an AD.

**Figure 11.** Normalized stationary parametric amplitude W of the tulip for *χ*nT = 0.25, where <sup>μ</sup> = 0.623 <sup>×</sup> <sup>10</sup><sup>−</sup>4. (**a**) <sup>ζ</sup> = 0.0016. (**b**) <sup>ζ</sup> = 0.0116.

**Figure 12.** Normalized stationary parametric amplitude W of the tulip for *χ*nT = 0.25, where <sup>μ</sup> = 0.623 <sup>×</sup> <sup>10</sup><sup>−</sup>4. (**a**) <sup>ζ</sup> = 0.0216. (**b**) <sup>ζ</sup> = 0.0318.

When analyzing the stationary parametric phase angle presented in the graphs in Figures 13 and 14 for ζ = 0.0016–0.0318 and *χ*nT = 0.25, it can be seen that for the cases, we had a manifestation of a "soft spring" with one branch for η<sup>1</sup> ≤ 2Ω1, while for η<sup>1</sup> ≥ 2Ω1, two branches of a "hard spring" exist, but with the increase in the damping ratio passing through the principal parametric resonance for η<sup>1</sup> > 2Ω1, and in the region of a "hard" spring, the two branches are approaching one against the other. In addition, the manifestation of increasing the value of the minimum of the second branch was sensitive to the increase in the value of the damping ratio in the range mentioned above (0.0016–0.0318). Through experiments, Steinwede demonstrated that the nonlinear parametric dynamic behavior of automotive driveshafts is like that of geared systems [4] (p. 117), [32] (pp. 132–160), and therefore, we observed similar pitting phenomena as well as micro-cracks inside the tulip and the bowl of the CVJ tulip-tripod and bowl-inner race joints [6] (pp. 88–94), which was also mentioned in the literature [13]. These aspects will "conduct" dynamic behavior through a chaotic dynamic that has, as a practical effect, an accelerating effect on pitting phenomena and cracking, as mentioned by Steinwede [4]

(pp. 88–94) or, in the worst case, results in the manifestation of cracks followed by failure (breaking) of the tulip [4] (p. 89). Figures 15 and 16 present the variation of the normalized stationary amplitude for the FBV in the PPR region for the bowl of an AD (around 6306.6 Hz). The fundamental natural frequency in bending for the bowl (3153.3 Hz) was computed with Equation (57). As can be seen, the natural frequency for the bowl was much bigger than the natural frequency of the tulip, and the PPR region for the bowl was in the PPR region for the forced torsional vibrations of the tulip, as mentioned in the literature [29].

**Figure 13.** Stationary parametric phase angle <sup>Ψ</sup> of the tulip for *<sup>χ</sup>*nT = 0.25, where <sup>μ</sup> = 0.623 <sup>×</sup> <sup>10</sup>−4. (**a**) ζ = 0.0016. (**b**) ζ = 0.0116.

**Figure 14.** Stationary parametric phase angle <sup>Ψ</sup> of the tulip for *<sup>χ</sup>*nT = 0.25, where <sup>μ</sup> = 0.623 <sup>×</sup> <sup>10</sup>−4. (**a**) ζ = 0.0216. (**b**) ζ = 0.0318.

**Figure 15.** Normalized stationary parametric amplitude W of the bowl for *χ*nB = 0.10, where <sup>μ</sup> = 0.754 <sup>×</sup> <sup>10</sup><sup>−</sup>4. (**a**) <sup>ζ</sup> = 0.00016 (**b**) <sup>ζ</sup> = 0.00016.

**Figure 16.** Normalized stationary parametric amplitude W of the bowl for *χ*nB = 0.10, where <sup>μ</sup> = 0.754 <sup>×</sup> <sup>10</sup><sup>−</sup>4. (**a**) <sup>ζ</sup> = 0.0216 (**b**) <sup>ζ</sup> = 0.0318.

When analyzing the normalized stationary parametric amplitude of the bowl presented in the graphs in Figures 15 and 16 for ζ = 0.0016–0.0318 and *χ*nB = 0.10, the same manifestation in the PPR region for the bowl as for the tulip can be seen but in a different range frequency. The manifestation of a "soft spring" with one branch for η<sup>3</sup> ≤ 2Ω<sup>3</sup> can be noted, indicating the presence of interaction between the principal parametric resonance and the primary resonance of the bowl, while for η<sup>3</sup> ≥ 2Ω3, two branches of a "hard spring" existed, which indicates the manifestation of interaction with the combination resonance or with the internal resonance for the bowl [28] (pp. 132–160). This aspect agrees

with the experimental data in the literature [4] (pp. 130–144). It can also be remarked that the increase in the damping ratio in the range of 0.0016–0.0318 increased the minimum normalized stationary amplitude of the bowl from 0.18 to 0.45 for one branch, and for the second branch, the increase in the minimum normalized stationary amplitude was from 0.29 to 0.49. This can be explained by the existence of dynamic instability in the region of principal parametric resonance, as mentioned in the literature [28] (pp. 132–160). Figures 17 and 18 present the variation of the stationary phase angle for the FBV in the PPR region for the AD bowl.

**Figure 17.** Stationary parametric phase angle <sup>Ψ</sup> of the bowl for *<sup>χ</sup>*nB = 0.10, where <sup>μ</sup> = 0.754 <sup>×</sup> <sup>10</sup>−4. (**a**) ζ = 0.00016 (**b**) ζ = 0.00016.

**Figure 18.** Stationary parametric phase angle <sup>Ψ</sup> of the bowl for *<sup>χ</sup>*nB = 0.10, where <sup>μ</sup> = 0.754 <sup>×</sup> <sup>10</sup>−4. (**a**) ζ = 0.0216 (**b**) ζ = 0.0318.

When analyzing the stationary parametric phase angle for the bowl presented in the graphs in Figures 17 and 18 for ζ = 0.0016–0.0318 and *χ*nB = 0.10, the same kind of manifestation as that for the tulip, but in a different range frequency, can be seen, and for the cases, we had a manifestation of a "soft spring" with one branch for η<sup>3</sup> ≤ 2Ω3, while for η<sup>3</sup> ≥ 2Ω3, two branches of a "hard spring" existed, but with the increase in the damping ratio, the passing through principal parametric resonance was found for η<sup>3</sup> > 2Ω3, and in the region of the "hard" spring, the two branches were approaching one against the other. For the bowl, this aspect of the two branches closing one against the other was more accentuated. Additionally, the manifestation of increasing the value of the minimum of the second branch was sensitive to the increase in the value of the damping ratio in the range mentioned above (0.0016–0.0318). Figure 19a,b illustrates the stationary dynamic instability region of the tulip in the space (η1, ξ, μ) using MATLAB software and Equation (81). When analyzing Figure 19a,b, it can be noticed that the two folded surfaces obtained were symmetrical concerning the plan given by (η1, ξ) for the excitation frequency and the damping ratio. In contrast, the excitation coefficient μ could be positive or negative. The folded surface of the dynamic instability "kept" inside the two branches the region where it manifested the stationary instability. It can also be seen that increasing the damping ratio had a stabilizing effect on the dynamics of the tulip, as expected. Figure 20a,b illustrates the stationary dynamic instability region of the bowl in the space (η3, ξ, μ) using MATLAB software and Equation (81).

**Figure 19.** Stationary dynamic instability surface frontier for the tulip. (**a**) μ > 0. (**b**) μ < 0.

By analyzing Figure 20a,b, it can be noticed that the two folded surfaces obtained were symmetrical concerning the plan given by (η3, ξ) for the excitation frequency and the damping ratio, while the excitation coefficient μ could be positive or negative. The folded surface of the dynamic instability "kept" inside the two branches the region where it manifested the stationary instability for the bowl. It can also be seen that increasing the damping ratio had a stabilizing effect on the dynamic of the bowl, as expected. The only difference for these phenomena for the bowl was that the manifestation was in another range frequency than that of the tulip: the range frequency given by the natural frequency of the bowl in bending.

Figure 21 illustrates the nonstationary amplitude spectrum W of the FBVM (forced bending vibrating movement) for the tulip during transition through the PPR region.

**Figure 20.** Stationary dynamic instability surface frontier for the bowl. (**a**) μ > 0. (**b**) μ < 0.

**Figure 21.** Nonstationary parametric amplitude spectrum W of the tulip for *χ*nT = 0.25, where <sup>μ</sup> = 0.623 <sup>×</sup> <sup>10</sup>−4. (**a**) <sup>ζ</sup> = 0.0016. (**b**) <sup>ζ</sup> = 0.0116. (**c**) <sup>ζ</sup> = 0.0216. (**d**) <sup>ζ</sup> = 0.0318. Here, m = 8 <sup>×</sup> <sup>10</sup>−<sup>5</sup> octaves/min.

It may be noted that while the excitation frequency was varied logarithmically through the resonant regime, the load parameter μ was kept constant. A vast number of numerical tests was performed to arrive at a reasonable value for the sweep rate of the excitation frequency m = 0.00008 octaves/min. When analyzing the nonstationary parametric amplitude spectrum presented in Figure <sup>20</sup> for <sup>ζ</sup> = (1.6–31.8) × <sup>10</sup>−<sup>3</sup> and *<sup>χ</sup>*nT = 0.25, the manifestation of beating effects as well as a significant decrease in the amplitude peaks with the increase in the damping ratio value in the range mentioned above can be seen. The maximum peak of the amplitude was obtained when the excitation frequency was in the close vicinity of PPR 1038 Hz. Each time the damping ratio increased at a rate of 0.01, the maximum peak amplitude decreased by more than five times. Figure 22 illustrates the nonstationary velocity amplitude spectrum dW dt of the forced bending vibrating movement (FBVM) for the tulip during transition through the PPR region.

**Figure 22.** Nonstationary parametric velocity amplitude spectrum dW dt of the tulip for *χ*nT**=** 0.25, where <sup>μ</sup> = 0.623 <sup>×</sup> <sup>10</sup>−4. (**a**) <sup>ζ</sup> = 0.0016. (**b**) <sup>ζ</sup> = 0.0116. (**c**) <sup>ζ</sup> = 0.0216. (**d**) <sup>ζ</sup> = 0.0318. Here, m=8 <sup>×</sup> <sup>10</sup>−<sup>5</sup> octaves/min.

When analyzing the graphs in Figure 22, the moment of passing in close vicinity of the PPR by a sudden increase in the peak velocity amplitude can be seen, while the manifestation of beatings was like the nonstationary amplitude spectrum (see Figure 20). In the range of 0.0016–0.0116, the increase in the damping ratio had the effect of increasing the maximum peak of the nonstationary velocity amplitude by approximately 15%, while in the range of 0.0116–0.0318, the increase in the damping ratio had the beneficial effect of decreasing the maximum peak by 72%. It can be concluded that the growth of the damping ratio had a beneficial impact through a reduction of the nonstationary velocity amplitude only in the range of 0.0116–0.0318. Figure 23 illustrates the graphs of the phase space W, dW dt for the nonstationary FBVM of the tulip in the PPR region.

**Figure 23.** Phase space W, dW dt of the tulip for *<sup>χ</sup>*nT **<sup>=</sup>** 0.25, where <sup>μ</sup> = 0.623 <sup>×</sup> <sup>10</sup>−4. (**a**) <sup>ζ</sup> = 0.0016. (**b**) <sup>ζ</sup> = 0.0116. (**c**) <sup>ζ</sup> = 0.0216. (**d**) <sup>ζ</sup> = 0.0318. Here,m=8 <sup>×</sup> <sup>10</sup>−<sup>5</sup> octaves/min.

By analyzing the graphs in Figure 23, a transition to chaotic behavior for the FBVM of the tulip due to the presence of limit cycles or even of strange attractors can be seen, but this last conclusion needs to be certified by detailed analysis using the methods of chaotic

movements. One aspect that is evident is that the transition through the PPR region for the tulip was an unstable one. Figure 24 illustrates the nonstationary amplitude spectrum W of the FBVM for the bowl during the transition through the PPR region.

**Figure 24.** Nonstationary parametric amplitude spectrum W of the bowl for *χ*nB = 0.10, where <sup>μ</sup> = 0.754 <sup>×</sup> <sup>10</sup>−4. (**a**) <sup>ζ</sup> = 0.0016. (**b**) <sup>ζ</sup> = 0.0036. (**c**) <sup>ζ</sup> = 0.0076. (**d**) <sup>ζ</sup> = 0.0096. Here, m = 1.4 <sup>×</sup> <sup>10</sup>−<sup>5</sup> octaves/min.

As for the tulip, a vast number of numerical tests was performed to arrive at a reasonable value for the sweep rate of the excitation frequency m = 1.4 × <sup>10</sup>−<sup>5</sup> octaves/min. By analyzing the nonstationary parametric amplitude spectrum presented in Figure 24 for <sup>ζ</sup> = (1.6–31.8) × <sup>10</sup>−<sup>3</sup> and *<sup>χ</sup>*nB= 0.10, the manifestation of beating effects as well as a significant decrease in the amplitude peaks with the increase in the damping ratio value

in the range of (1.6–11.6) × <sup>10</sup>−<sup>3</sup> can be found, while in the range of (11.6–31.8) × <sup>10</sup>−3, the decrease was more significant. The maximum peak of the amplitude was obtained when the excitation frequency was in close vicinity of the PPR for the bowl (6306 Hz). Each time the damping ratio increased at a rate of 0.01, the maximum peak amplitude decreased by more than 10 times (5 times for the tulip), because the bowl had a more considerable rigidity than the tulip. Figure 25 illustrates the nonstationary velocity amplitude spectrum dW dt of the FBVM for the bowl during the transition through the PPR region.

**Figure 25.** Nonstationary parametric velocity amplitude spectrum dW dt of the bowl for *χ*nT **=** 0.25, where <sup>μ</sup> = 0.754 <sup>×</sup> <sup>10</sup>−4. (**a**) <sup>ζ</sup> = 0.0016. (**b**) <sup>ζ</sup> = 0.0036. (**c**) <sup>ζ</sup> = 0.0076. (**d**) <sup>ζ</sup> = 0.0096. Here, m = 1.4 <sup>×</sup> <sup>10</sup>−<sup>5</sup> octaves/min.

When analyzing the graphs in Figure 25, the moment of passing in close vicinity of the PPR was marked by a sudden increase in the peak of velocity amplitude, while the manifestation of beatings was like that of the nonstationary velocity amplitude spectrum for the tulip (see Figure 22). The manifestation for the bowl was different than that for the tulip (see Figure 22) because from the beginning, in the whole range of the damping ratio (0.0016–0.0318), the increase in the damping ratio had the effect of decreasing the maximum peak of the nonstationary velocity amplitude by approximately 5 times in the range of 0.0016–0.0116, while in the range of 0.0216–0.0318, the increase in the damping ratio had the effect of decreasing the maximum peak by more than 8 times. It can be concluded that the growth of the damping ratio had the beneficial impact of reducing the nonstationary velocity amplitude in the whole range (0.0016–0.0318) of the damping ratio. Figure 26 illustrates the graphs of phase space W, dW dt for the nonstationary FBVM of the bowl during the transition through the PPR region.

**Figure 26.** Phase space W, dW dt of the bowl for *<sup>χ</sup>*nB **<sup>=</sup>** 0.10, where <sup>μ</sup> = 0.754 <sup>×</sup> <sup>10</sup>−4. (**a**) <sup>ζ</sup> = 0.0016. (**b**) <sup>ζ</sup> = 0.0036. (**c**) <sup>ζ</sup> = 0.0076. (**d**) <sup>ζ</sup> = 0.0096. Here, m = 1.4 <sup>×</sup> <sup>10</sup>−<sup>5</sup> octaves/min.

When analyzing the graphs in Figure 26, the same manifestation for the bowl as for the tulip (see Figure 23) can be seen, that being a transition to chaotic behavior for the FBVM of the bowl due to the presence of limit cycles or even strange attractors. However, this last conclusion needs to be certified by detailed analysis using the methods of chaotic movements. This manifestation is valid only in the damping ratio range of 0.0016–0.0096. One aspect that is evident is how the transition through the PPR region for the bowl was also an unstable one.

Unfortunately, there are no published studies which analyze in detail the dynamic behavior of each element of the AD for the FBVM (the tulip, bowl, and midshaft) apart from [4], as all the studies analyzed the global dynamic behavior of the automotive driveshaft. As can be seen from this paragraph, the use of AMA coupled with Hamilton's principle allowed the investigation of the stationary motion for the FBV movements of an AD's tulip and bowl in the PPR region by computing the amplitude (see Figures 10, 11, 14 and 15, the phase angle (see Figures 7, 12, 13 and 16), and the dynamic instability frontiers (see Figures 18 and 19) in the PPR region. In the meantime, the use of the AMA coupled with Hamilton's principle allowed the investigation of the nonstationary motion for the FBV movements of an AD's tulip and bowl in the transition through a PPR region by computing the amplitude spectrum (see Figures 21 and 24), the velocity amplitude spectrum (see Figures 22 and 25), and the velocity amplitude versus the amplitude in the phase space W, dW dt (see Figures 23 and 26). Figures 23 and 26 allowed the investigation of the dynamic instability in the transition through the PPR region. As is noted in Figures 21–26, the transition to the PPR region had an aspect of chaotic manifestation due to the manifestation of beating effects. To check whether this nonstationary dynamic behavior is deterministic chaos or a stochastic process, it would be necessary to use Lyapunov's exponents method coupled with the Poincare map method. All these aspects represent novelties in the dynamic instability investigation for the AD. Based on the phenomena presented above, the direction for future research involves the study of chaotic motion generated by the FBVM of the AD elements. Another future research direction is to investigate the dynamic instability in the proximity of each possible resonance type for the stationary and nonstationary FBVM of the AD elements, such as subharmonic resonance, superharmonic resonance, and combination and internal resonance.

#### **8. Conclusions**

The present work introduces a newly designed DMFFBVM for the AD, with the following phenomena being included for the first time: nonuniformity of the inertial characteristics of the AD's elements (see Equations (1)–(18)), serial stiffness and damping for the tulip and bowl (see Equations (19) and (20) as well as (23) and (24)), shock excitation due to the road geometry (see Equation (25)), and nonuniformity of the kinematic isometry (see Equation (61)). Based on this newly designed DMFFBVM, using Hamilton's principle coupled with the first-order approximation of the AMA, the stationary and nonstationary dynamic instability behavior of the AD elements were investigated in detail by computing the following:


an AD's durability. Moreover, the DMFFBVM must be added to the design algorithm for predicting the comfort elements of automobiles.

**Author Contributions:** Conceptualization, M.B.; methodology, M.B. and O.V.; software, M.B. and O.V.; validation, M.B.; formal analysis, M.B.; investigation, M.B.; resources, M.B. and O.V.; data curation M.B.; writing—original draft preparation, M.B.; writing—review and editing, M.B.; visualization, M.B.; supervision, M.B.; project administration, M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors of this article are thankful to the University POLITEHNICA of Bucharest for providing a serene environment and facilities for carrying out this research.

**Conflicts of Interest:** The authors declare no potential conflicts of interest with respect to the research, authorship, or publication of this article.
