Application of Multiple-Scales Method for the Dynamic Modelling of a Gear Coupling

Abstract

Thin-walled gears, designed for aeronautical applications, have shown very rich dynamics that must be investigated in advance of the design phase. One of the signatures of their dynamics is coupling due to the meshing teeth which stand-alone gear models cannot capture. This paper aims to investigate the dynamics of thin-walled gears considering time-varying coupling due to the gear meshing. Each gear is modelled with lumped parameters according to a local rotating reference system and the coupling is modelled by a traveling meshing stiffness. The set of equations of motion is solved by the non-linear Method of Multiple-Time-Scales (MMTS). MMTS is a very powerful technique that is widely used to solve perturbation problems in many fields of mathematic and physics. In the analyzed numerical test case, the relevance of gear coupling is demonstrated as well as the capability of the MMTS to capture the fundamental features of the system dynamics. In this study the analytical methodology, which uses MMTS, allows for the calculation of the forced response of the system made of two meshing gears despite the presence of a parametric quantity, i.e., the mesh stiffness. The calculation is performed in the frequency domain using modal coordinates, which ensures a fast computation. The result is compared with time domain analysis for validation purposes.

1. Introduction

The need to identify a non-linear methodology for a dynamic study of two meshing gears moves from the evidence of some critical resonances occurring during operations, which cannot be investigated by analyzing a single gear considered as a stand-alone component, but it requires the analysis of the overall system which can be made of two or more than two meshing gears (planetary system), where time-varying parameters and non-linearities appear in the equations of motion. In practice, it is experimentally verified that one gear can influence the dynamics of the other meshing gears, under certain conditions, causing unexpected resonances, which are dangerous for the overall system. Then, a dynamic coupling is established between the meshing gears. This phenomenon is mainly due to the fluctuation of the mesh stiffness at the meshing teeth, which varies because of the different mesh conditions and contact points during meshing. Fluctuations of the mesh stiffness can induce severe instability conditions and affect also the resonances of the system. The phenomenon of dynamic coupling can be experimentally verified in industrial applications, in particular for aeronautical applications where the gears, having specific mechanical characteristics and working at critical speed regimes, show mutual interactions, which largely affect the forced response of the system. In more detail, the dynamic coupling causes critical resonances on a gear, which are induced by the excitation of the mode shapes of the meshing gear. The presence of these mutual interactions among the components leads to a different study of the system, which must include all the gears involved in the interactions. The dynamic coupling is a direct consequence of the variations of tooth flexibility (mesh stiffness) because of the different contact conditions during rotation, but also because of the variations of the contact ratio. Thus, time-varying mesh stiffness causes the system to be non-linear. It is worth remembering that here the term “non-linear” is adopted to highlight the fact that the system is not the usual Linear Time-Invariant System (LTIS), but a Linear Time-Variant System (LTVS). Being the system of the type LTVS, it is not possible to compute the forced response by inverting the dynamic stiffness matrix as for a LTIS, since the latter is a time-variant parameter inside the equation of motion. Nevertheless, the dependent variable of the equation of motion does not show elevation to powers or other non-linearities (e.g., Duffing equation). Then, the superimposition effect principle is valid for such a system and it will be used in the following discussion. As a consequence, the adoption of a “non-linear” method is needed to compute the forced response of the system. The dynamics of gear systems have been extensively studied by researchers for decades and still represent an important matter of interest for the understanding of phenomena affecting the dynamics of such systems. The evaluation of the mesh stiffness variations and the related non-linear aspects have a primary importance for researchers who provide several modelling solutions of the phenomenon [1] for mathematical definition. According to the different types of modelling of the mesh stiffness, many works provided different methodologies for the computation of the response of the system, according to different levels of complexity. Most of the works focused on the combined effect of mesh stiffness variation and backlash between the meshing teeth, which affects largely the response of the meshing gears, developing non-linear methodologies for the iterative and numerical computation of the response [2,3,4,5,6,7,8,9]. Other studies focused on the analysis of the instability conditions, which can be caused by the fluctuations of the mesh stiffness involving sometimes wide operational speed ranges of the gears and can lead to failures. Works on instability provide an analytical solution, using perturbation methods (e.g., method of multiple time-scales, MMTS) to establish relations between the analyzed instability conditions and the entity of the mesh stiffness fluctuations [10]. Recently, instability analyses and forced response studies were extended to more complex systems like planetary gear systems [11,12,13,14]. Most of the cited works analyze the dynamics of a gear system by considering the gears as rigid bodies connected by the mesh stiffness and introducing the transmission error between two meshing gears, which consider the fluctuations of an equivalent tooth compliance that excites the system.

In this paper, the aim is to consider the gears as compliant bodies and compute analytically the forced response of the system excited both parametrically and externally. The backlash phenomenon is not considered at this stage in order to focus the attention on the phenomenon of the dynamic coupling and on the method to be developed to study the phenomenon without the nonlinearity introduced by intermitting contacts during meshing. Here, transmission error cannot be used anymore since the gear bodies are considered as compliant. The gears, which constitute the overall system, are linked together by means of a time-variant mesh stiffness, which acts on the nodes of the teeth, where the contact takes place. In other words, the system sees both a parametric excitation and an external force exciting the system. The methodology developed here applies the Method of Multiple-Time-Scale (MMTS) to compute the frequency response of a single mesh gear pair, modelled with lumped parameters, and investigate the dynamic coupling, which is established between the gears, verifying the mutual interactions and resonances induced by the phenomenon. MMTS allows a good approximation to the solution of the problem by introducing “scales” variables, which will substitute the independent variable of the problem. The solution of the problem passes through the elimination of the so-called “secular terms”. This procedure represents a necessary solvability conditions for the solution of the problem. In this paper, numerical examples of forced response are reported, based on test cases. Upon these test-case analyses, the methodology is finally validated by means of direct time integration (DTI) of the non-linear equations of motion.

2. Model of the System

The system under analysis is made of two meshin_g gears (Gear-1 or G_1, and Gear-2 or G₂, Figure 1). For each gear, a local reference system rotating with the gear itself is defined. Each gear is divided into sectors (Z₁ for Gear-1, and Z₂ for Gear-2), one per each tooth (Figure 2). A gear ratio η can be defined for the system under analysis as the ratio between Z₁ and Z₂. Each sector is modelled as a lumped parameter model with two degrees of freedom (dof), or nodes, one for the tooth and one for the gear sector wheel (Figure 3). The rotation of each gear around its own axis is allowed (no radial or axial displacement are allowed, only tangential displacement is allowed). The latter assumption is reasonable for the case of thin walled spur gears where radial and axial displacements can be assumed as negligible. The sectors are then coupled together. The periodic coupling between the teeth of the two gears is modelled by a time-variant mesh stiffness K_M(t), described in more detail in Section 3.

As shown in Figure 3, the two gears are constrained to ground by means of the stiffness elements $k_{a,1}$ and $k_{a,2}$ respectively. The mechanical characteristics of mass $m_{b,1}$ and $m_{b,2},$ and stiffness $k_{b,1}$ and $k_{b,2}$ are associated to the teeth of the gears, whose displacement coordinates are $x_{G_1,t}$ and $x_{G_2,t}$ (respectively for the teeth of G₁ and the teeth of G₂). The mechanical characteristics of mass $m_{c,1}$ and $m_{c,2},$ and stiffness $k_{c,1}$ and $k_{c,2}$ are associated to the gear sector wheel, whose displacement coordinates are $x_{G_1,c}$ and $x_{G_2,c}$ (respectively for the sector wheel of G₁ and the sector wheel of G₂). In Figure 3, a force acting on the meshing teeth dof (equal but with opposite direction for the tooth of G₁ and the tooth of G₂) is shown, which represents the force establishing between them when the meshing couple is in contact. Obviously, when the couple is not meshing, no forces are exchanged. Then, the force travels along the circumference of each gear as the system rotates. The excitation force will be discussed more in detail in Section 5.

The physical displacement vector (with the corresponding size) is:

$$\left\{x\right\}={\left\{{\left\{x_{G_1}\right\}}^T,{\left\{x_{G_2}\right\}}^T \right\}}^T{}_{1xN},$$

(1)

where,

$${\left\{x_{G_1}\right\}}^T={\left\{{\left\{x_{G_1,c}\right\}}^T, {\left\{x_{G_1,t}\right\}}^T\right\}}^T{}_{1 x 2Z_1},$$

(2)

$${\left\{x_{G_2}\right\}}^T={\left\{{\left\{x_{G_2,c}\right\}}^T, {\left\{x_{G_2,t}\right\}}^T\right\}}^T{}_{1 x 2Z_2},$$

(3)

with $N=2 Z_1+2 Z_2$. $\left\{x\right\}$ including Gear-1 displacement coordinates, $x_{G_1}$, subdivided into $x_{G_1,c}$ which indicates the displacement of the nodes of the gear wheel, and $x_{G_1,t}$ indicating the displacement of the teeth. The same holds for $x_{G_2}$ of Gear-2. The equation of motion, in matrix form, can be written in general as:

$$M \ddot{x}+\widehat{C} \dot{x}+\widehat{K}\left(t\right) x=\widehat{F}\left(t\right),$$

(4)

where $M$ is the mass matrix; $\widehat{C}$ is the damping matrix; $\widehat{K}\left(t\right)$ is the stiffness matrix which includes time-variant parameters, corresponding to the mesh stiffness K_M(t) used to couple the two gears; $\widehat{F}\left(t\right)$ is the force vector containing non-zero values for the teeth dof. As anticipated before, $\widehat{F}\left(t\right)$ is a time-variant vector since the mesh force passes from one tooth to another as the system rotates. Then, each tooth is periodically subjected to a force excitation due to the meshing, where the period is equal to the rotation period of the gear. In the next section, the mesh stiffness will be discussed, then the assembly of the matrices will be presented in Section 4.

3. Definition of Mesh Stiffness

During meshing, many factors can induce fluctuations of the stiffness characteristics of the teeth. As explained before, the fluctuations of the mesh stiffness can be due to different contact conditions given by different contact ratios and contact positions along the tooth face. Hertzian contact phenomena can also influence the stiffness of the teeth. The combined effect of all these fluctuation sources produces a time history of the mesh stiffness acting on a single tooth. In this paper the time history of the mesh stiffness is not investigated in detail and it is approximated to a rectangular waveform traveling from one tooth to another one. In more detail, the mesh stiffness, which couples the n^th tooth pair, is assumed to have a constant value k_t when the n^th tooth pair is in contact and a null value when the contact is missing. Within the meshing time interval, the constant value, k_t, assumed by the mesh stiffness can represent an equivalent mean value of a real trend during meshing. Since rotating reference systems are used, in each gear the mesh stiffness rotates with the same speed as the gear but in the opposite direction. In Figure 4 the time history of a generic mesh stiffness K_M(t) is shown with equivalent value k_t and unitary contact ratio, acting on a single tooth of a gear, with Z sectors (teeth) rotating at certain speed with revolution period T. In this qualitative example of mesh stiffness, one can distinguish between the meshing time interval, when the tooth is in contact, from the rest of the time history when the tooth is not in contact and the mesh stiffness assumes a null value. Once a full revolution is performed, so after a period T, the tooth experiences again the mesh stiffness.

The rectangular waveform can be translated into a sum of harmonics by developing the Fourier Series of time trend:

$$K_M\left(t\right)=K_c+{{\displaystyle \sum }}_{s=1}^{\infty }[K_{Va}^s\cos \left(s\Omega t\right)+K_{Vb}^s\sin \left(s\Omega t\right)],$$

(5)

where $K_c$ is the mean value of the function, $s$ is the harmonic index, $\Omega$ is the speed of the gear ($T=2\pi /\Omega$), while $K_{Va}^s$ and $K_{Vb}^s$ are the coefficients of the “cosine” harmonics and the “sine” harmonics, respectively. The expression of the Fourier series can be further manipulated by means of the Euler formula, to redefine Equation (5) as the real-valued form of the complex notation of the Fourier Series which will be used in the following discussion (Equation (6)).

$$K_M\left(t\right)=K_c+K_V\left(t\right)=K_c+{{\displaystyle \sum }}_{s=1}^{\infty }\left[K_V^s{e}^{i s\Omega t}+\overline{K_V^s} e^{-i s\Omega t}\right],$$

(6)

where:

$$K_V^s=\frac{1}{2}\left(K_{Va}^s-i{K}_{Vb}^s\right); \overline{K_V^s},\mathrm{complex}\mathrm{conjugate}\mathrm{of}{K}_V^s.$$

(7)

In Figure 5a,b a numerical example is shown (with Z₁ = 10, Z₂ = 20 and k_t = 10⁶ N/m ) where T₁ and T₂ are the revolution periods of Gear-1 and Gear-2 respectively, and they are different from one another since Z₁ ≠ Z₂.

By looking at the previous example, it is easy to note that the periodic time history of the mesh stiffness of the two gears is different. As a matter of fact, the rectangular waveform is the same, but its period differs in the two gears, being T₁ for Gear-1 and T₂ for Gear-2. Here, this concept is clarified with an example. Let us consider the previous system with Z₁ = 10 teeth and Z₂ = 20 teeth. In Figure 6a the couple 1-1 (tooth-1 of G₁—tooth-1 of G₂) starts meshing for an angle θ₁ = 0° of G₁. After a full revolution of G₁ (Figure 6e) tooth-1 of G₁ meshes again but now with tooth-11 of G₂. The same is for the other couples 2-2 and 2-12 (Figure 6b–f) and so on. Thus, a certain couple (i-j) meshes with a base period that is twice the base period T₁. Of course, the latter relation changes for systems with different number of teeth.

4. Equations of Motion and Construction of the Matrices

Having adopted a lumped parameters model, it is convenient to write, for clarity’s sake, an equation of motion (considering at this stage the un-damped unforced system) of one dof connected to the i_th tooth of Gear-1. Then, let us consider the previous example of a system constituted by one gear (Gear-1) having 10 teeth (Z₁) and a second gear (Gear-2) having 20 teeth (Z₂). Let us write the equation of motion of tooth-1 of Gear-1. By looking at Figure 7 below, neglecting at this stage the presence of damping and excitation force, the resulting equation of motion is:

$$m_b{\ddot{x}}_{t_1,1}+k_b\left(x_{t_1,1}-x_{c_1,1}\right)+K_{M_{1,1}}\left(t\right) \left(x_{t_1,1}-x_{t_2,1}\right)+K_{M_{1,11}}\left(t\right) \left(x_{t_1,1}-x_{t_2,11}\right)=0,$$

(8)

where $m_b$ refers to mass of the teeth; $k_b$, nominal stiffness of the teeth; $K_{M_{1,1}}\left(t\right)$ and $K_{M_{1,11}}\left(t\right)$ mesh stiffness coupling the 1-1 teeth pair and 1-11 teeth pair respectively.

It is possible to write the latter equation in a clearer form as follows:

$$m_b{\ddot{x}}_{t_1,1}+k_b{x}_{t_1,1}+\left(K_{M_{1,1}}\left(t\right)+ K_{M_{1,11}}\left(t\right)\right)x_{t_1,1}-K_{M_{1,1}}\left(t\right) x_{t_2,1}-K_{M_{1,11}}\left(t\right) x_{t_2,11}-k_b{x}_{c_1,1}=0.$$

(9)

In Equation (9), two types of time-variant stiffness can be highlighted. The first type includes $K_{M_{1,1}}\left(t\right)$ and $K_{M_{1,11}}\left(t\right)$. The two terms refer to a specific teeth pair that meshes during operation. They are characterized by a rectangular waveform having a period $T_{pair}$ that is strictly dependent on the number of teeth of the two gears. It is easy to derive $T_{pair}$ and to relate it to the revolution periods of the two gears, $T_1$ and $T_2$ respectively. In general, $T_{pair}$ is a multiple of both the periods $T_1$ and $T_2$. It is convenient to write the latter relation in a mathematical form:

$$T_{pair}=n_{T_1}{T}_1=n_{T_2}{T}_2$$

(10)

where the coefficients $n_{T_1}$ and $n_{T_2}$ define the multiple of the respective revolution periods. In the example analyzed here it is easy to note that $n_{T_1}=2$ and $n_{T_2}=1$. In other words, a specific pair meshes after two revolution of the Gear-1 as well as after one revolution of Gear-2. The other type of variable stiffness term which appears into Equation (9) is the term $K_{M_{1,1}}\left(t\right)+K_{M_{1,11}}\left(t\right)$. This sum creates a rectangular waveform different from the previous one. As a matter of fact, $K_{M_{1,1}}\left(t\right)+K_{M_{1,11}}\left(t\right)$ is a waveform with the same constant value k_t but with a period exactly equal to the revolution period $T_1$, as shown in Figure 8.

Then, for all the equations of motion of the dof belonging to Gear-1, two different sets of harmonics appear into the equation. The fundamental frequencies, which describe the two sets, are respectively:

$${\mathsf{\Omega }}_1=\frac{2\pi }{T_1},$$

(11)

$${\mathsf{\Omega }}_3=\frac{2\pi }{T_{pair}}$$

(12)

Following the same approach, it is possible to write the equations of motion for the dof belonging to Gear-2. For those equations, still two types of variable stiffness can be distinguished into the resulting equation (here the equation of motion is not reported because is similar to Equation (9). The first type is again the rectangular waveform of the single pair of teeth described by the period $T_{pair}$, so by the fundamental frequency ${\mathsf{\Omega }}_3$. The second type of time-variant stiffness term is the sum of all the other pair waveforms involved into the equation, but now the resulting rectangular waveform has a base period equal to the revolution period of Gear-2, i.e., $T_2$. Therefore, another set of harmonics appears in the equations of motion of the system and it is described by the fundamental frequency ${\mathsf{\Omega }}_2$, which is the speed of Gear-2

$${\mathsf{\Omega }}_2=\frac{2\pi }{T_2}$$

(13)

Finally, the construction of the stiffness matrix $\widehat{K}\left(t\right)$ allows to write the full stiffness matrix in the next form (Equation (14)), by properly distinguishing between the terms referred to the three different sets of harmonics, having fundamental frequencies ${\mathsf{\Omega }}_1$, ${\mathsf{\Omega }}_2$ and ${\mathsf{\Omega }}_3$:

$$\widehat{K}\left(t\right)=\left(K_c+\widehat{K_{V1}}\left(t\right)+\widehat{K_{V2}}\left(t\right)+\widehat{K_{V3}}\left(t\right)\right),$$

(14)

$$\widehat{K_{V1}}\left(t\right)={{\displaystyle \sum }}_{s_1=1}^{\infty }[K_{V1}^{s_1}{e}^{i{s}_1{\Omega }_{1}t}+\overline{K_{V1}^{s_1}}{e}^{-i{s}_1{\Omega }_{1}t}],$$

(15)

$$\widehat{K_{V2}}\left(t\right)={{\displaystyle \sum }}_{s_2=1}^{\infty }[K_{V2}^{s_2}{e}^{i{s}_2{\Omega }_{2}t}+\overline{K_{V2}^{s_2}}{e}^{-i{s}_2{\Omega }_{2}t}],$$

(16)

$$\widehat{K_{V3}}\left(t\right)={{\displaystyle \sum }}_{s_3=1}^{\infty }[K_{V3}^{s_3}{e}^{i{s}_3{\Omega }_{3}t}+\overline{K_{V3}^{s_3}}{e}^{-i{s}_3{\Omega }_{3}t}].$$

(17)

The mass matrix $M$ and the stiffness matrix $\widehat{K}\left(t\right)$ are obtained by writing the equation of motion for each dof of the system as shown in Equation (9) according to the matrix notation. The damping matrix $\widehat{C}$ will be defined later in Section 6.1 assuming a modal damping ratio $\varsigma$ to the mode shapes of the system.

5. Excitation Force

A mesh force applied to two meshing teeth excites the system. The mesh force F(t) of Figure 3 corresponds to a torque (for example applied on G₁) and a resistant torque (applied on G₂) acting on the system. Since the mesh force rotates along the gears as the system rotates, each tooth of a gear undergoes periodically the same mesh force, with a period that is equal to the rotation period of the gear itself ($T_1$ for teeth of G₁ and $T_2$ for the teeth of G₂) each tooth is subjected to the same force with a time delay. During meshing, the value of the mesh force is assumed to be constant. Then, considering the generic i_th tooth of a gear, the mesh force $F_G{}_i\left(t\right)$ acting on it will have the same trend of the mesh stiffness in the time domain (rectangular waveform, Figure 4). As a consequence, the force vector $\widehat{F}\left(t\right)$ in the equation of motion Equation (4) contains, in correspondence to the teeth dof, the mesh force trends $F_G{}_i\left(t\right)$, properly phased in the time domain according to the position of the tooth which is considered (e.g., if tooth-1 of G₁ undergoes $F_{1_1}\left(t\right)$, tooth-2 of G₁ undergoes $F_{1_2}\left(t\right)=F_{1_1}\left(t+∆T\right)$, being $∆T$ the meshing time interval). Finally, the force vector $\widehat{F}\left(t\right)$ can be represented as

$$\widehat{F}\left(t\right)={\left\{\begin{array}{c}\begin{array}{c}{\left\{0\right\}}_{Z_1\times 1}\\ {\left\{\begin{array}{c}{F}_{1_1}\\ ⋮\\ F_{1_{2Z_1}}\end{array}\right\}}_{Z_1\times 1}\end{array}\\ \begin{array}{c}{\left\{0\right\}}_{Z_2\times 1}\\ {\left\{\begin{array}{c}{F}_{2_1}\\ ⋮\\ F_{2_{2Z_2}}\end{array}\right\}}_{Z_2\times 1}\end{array}\end{array}\right\}}_{N\times 1}$$

(18)

As for the mesh stiffness, it is convenient to express the mesh force as a Fourier series, since it is periodic in the time domain. This will be advantageous for the computation of the forced response (Section 6), since each harmonic will be considered separately, computing the overall response as a superimposition of the effects due to each harmonic. As a matter of fact, the system is a linear time-variant system whereby the superimposition principle is valid. Then, let us write the expression of the force function, in the Fourier series. It is convenient to distinguish between the force acting on the teeth nodes of Gear-1 from the force acting on the teeth nodes of Gear-2. The two force sets, expressed in Fourier Series, have fundamental frequencies that are the speed of the Gear-1 (${\Omega }_1$) for the excitation force terms acting on the teeth nodes of Gear-1 and the speed of Gear-2 (${\Omega }_2$) for the excitation force terms acting on the teeth nodes of Gear-2. As a consequence, let us write the general expression of the force acting on the i_th tooth of G₁ and the force acting on the j_th tooth of G₂ respectively in Equations (19) and (20):

$$F_{1_i}\left(t\right)=F_{1_i,0}+{{\displaystyle \sum }}_{k_1=1}^{\infty }[F_{1_i}{}_a^{k_1}\cos \left(k_1{\Omega }_{1}t\right)+F_{1_i}{}_b^{k_1}\sin \left(k_1{\Omega }_{1}t\right)];$$

(19)

$$F_{2_j}\left(t\right)=F_{2_i,0}+{{\displaystyle \sum }}_{k_2=1}^{\infty }[F_{2_j}{}_a^{k_2}\cos \left(k_2 {\Omega }_{2}t\right)+F_{2_j}{}_b^{k_2}\sin \left(k_2{\Omega }_{2}t\right)].$$

(20)

As for the mathematical expression of the mesh stiffness (Section 3), The Fourier Series of the mesh force can be written by means of the exponential notation. Equations (19) and (20) become:

$$F_{1_i}\left(t\right)=F_{1_i,0}+{{\displaystyle \sum }}_{k_1=1}^{\infty }[ F_{1_i}^{k_1} e^{i{k}_1{\Omega }_{1}t}+\overline{F_{1_i}^{k_1}}{e}^{-i k_1{\Omega }_{1}t}],$$

(21)

$$F_{2_i}\left(t\right)=F_{2_i,0}+{{\displaystyle \sum }}_{k_1=1}^{\infty }[ F_{2_i}^{k_2} e^{i{k}_1{\Omega }_{1}t}+\overline{F_{2_i}^{k_2}}{e}^{-i k_2{\Omega }_{2}t}],$$

(22)

With

$$F_{G_i}^{k_G}=\frac{1}{2}\left(F_{G_i}{}_a^{k_G}-i F_{G_i}{}_b^{k_G}\right),$$

(23)

$$\overline{F_{G_i}^{k_G}}=\frac{1}{2}\left(F_{G_j}{}_a^{k_G}+i F_{G_j}{}_b^{k_G}\right).$$

(24)

Being the forcing functions $F_{1_i}\left(t\right)$ and $F_{2_j}\left(t\right)$ expressed as a series of harmonics with frequencies $k_1{\Omega }_1$ and $k_2{\Omega }_2$ respectively, it is convenient to separate the force vector $F\left(t\right)$ into a sum of two vectors whose components can be expressed as harmonics having frequencies $k_1{\Omega }_1$ and $k_2{\Omega }_2$, respectively ${\widehat{F}}_1\left(t\right)$ (Equation (25)) and ${\widehat{F}}_2\left(t\right)$ (Equation (26)). Thus, let us represent the two vectors as follows:

$${\widehat{F}}_1\left(t\right)={\left\{\begin{array}{c}\begin{array}{c}{\left\{0\right\}}_{Z_1\times 1}\\ {\left\{\begin{array}{c}{F}_{1_1}\\ ⋮\\ F_{1_{2Z_1}}\end{array}\right\}}_{Z_1\times 1}\end{array}\\ {\left\{0\right\}}_{2Z_2\times 1}\end{array}\right\}}_{N\times 1},$$

(25)

$${\widehat{F}}_2\left(t\right)={\left\{\begin{array}{c}{\left\{0\right\}}_{2Z_1\times 1}\\ \begin{array}{c}{\left\{0\right\}}_{Z_2\times 1}\\ {\left\{\begin{array}{c}{F}_{2_1}\\ ⋮\\ F_{2_{2Z_2}}\end{array}\right\}}_{Z_2\times 1}\end{array}\end{array}\right\}}_{N\times 1}.$$

(26)

The overall forcing function vector $F\left(t\right)$ is

$$\widehat{F}\left(t\right)={\widehat{F}}_1\left(t\right)+{\widehat{F}}_2\left(t\right).$$

(27)

6. Forced Response Computation with MMTS

Section 6 deals with the development of a general analytical solution of the forced response of the system under exam using MMTS. MMTS is a very used technique able to obtain approximations of solutions to non-linear problems. It works by substituting different “scales” variables (according to the level of approximation the user desires) to the independent variable of the equation, treating them as independent variables. Being a frequency-based method, it allows the study of the response of a complex system in a very flexible way, reducing considerably the computational time, with respect to other time-based methods which provide a numerical and iterative solution to the problem. Before going through the discussion of MMTS, it is advantageous to rewrite the equation of motion Equation (4) in terms of modal coordinates. The use of modal coordinates leads to a great simplification of the problem by decoupling the equation of motions. Moreover, it allows the direct evaluation of the effect of a harmonic component of the excitation force on the respective mode shape that is excited by that component. Thus, the following Paragraph 6.1 is dedicated to the transformation of the equation of motion in modal coordinates while the mathematical development of MMTS is discussed in Paragraph 6.2.

6.1. Modal Analysis and Transformation of the Equation of Motion in Modal Coordinates

Since the stiffness matrix $\widehat{K}\left(t\right)$ contains time-variant elements, the computation of modal analysis cannot be performed in general. For this reason, the computation of the modal analysis is performed on the mean part of the stiffness matrix $K_C$ (Equation (14)) ($K_C$ is a symmetric matrix). Thus, let us consider the time-invariant (unforced and undamped) part of the equation of motion Equation (4) (Equation (28)) and perform the modal analysis (Equation (29)).

$$M\ddot{x}+K_{C}x=0;$$

(28)

$$\det \left(K_C-{\omega }^{2}M\right)=0, eigenvalue problem \to {\omega }_n^2, {\psi }_n, n=1÷N;$$

(29)

$$\begin{array}{c}\Psi ={\left[{\psi }_1, \dots ,{\psi }_n,\dots ,{\psi }_N\right]}_{N\times N}, modal matrix\\ {\omega }_n, natural frequency\end{array}$$

From modal analysis, the natural frequencies ${\omega }_n$ and mode shapes ${\psi }_n$ of the (mean) overall system are obtained. Then, let us apply Direct Modal Transformation (DMT, Equation (30)) to the equation of motion using the modal matrix $\Psi$ and multiply the latter by the transpose of the modal matrix ${\Psi }^T$. The resulting equation written in modal coordinates ($u$, vector of modal coordinates) is reported in Equation (31).

$$x=\Psi u.$$

(30)

$$M_{mod} \ddot{u}+{\widehat{C}}_{mod} \dot{u}+K_{c,mod} u+\widehat{D}\left(t\right) u=\widehat{P}\left(t\right),$$

(31)

with: $M_{mod}$, modal mass matrix; ${\widehat{C}}_{mod}$, modal damping matrix; $K_{c,mod}$, modal mean stiffness matrix;

$$\widehat{D}\left(t\right)={\Psi }^T\widehat{K_V}\left(t\right)\Psi ;$$

(32)

$$\widehat{P}\left(t\right)={\Psi }^T\widehat{F}\left(t\right)$$

(33)

If the mode shapes are normalized with respect to the unitary modal masses, $M_{mod}$ is equal to the identity matrix and $K_{c,mod}$ is a diagonal matrix having on its main diagonal the eigenvalues of the system ${\omega }_n^2$. As anticipated in Section 2, damping is not modelled physically into the model. For the sake of simplicity, damping is introduced by means of the modal damping ratio ${\varsigma }_n$ to be associated to each n_th mode shape. It is possible to obtain the modal damping ${\widehat{c}}_n$ as:

$${\widehat{c}}_n=2\frac{{\varsigma }_n}{\sqrt{k_{mod}{}_n \times m_{mod}{}_n}}=2\frac{{\varsigma }_n}{\sqrt{{\omega }_n^2}},=1÷N;$$

(34)

with $m_{mod}{}_n$ modal mass of the n_th mode shape ($m_{mod}{}_n=1$, if the mode shapes are normalized with respect to the unitary modal masses) and $k_{mod}{}_n$ modal stiffness of the n_th mode shape ($k_{mod}{}_n={\omega }_n^2$, if the mode shapes are normalized with respect to the unitary modal masses). Then, the modal damping matrix ${\widehat{C}}_{mod}$ is a diagonal matrix having on its main diagonal the modal damping ${\widehat{c}}_n$, satisfying the following relation that allows computation of the damping matrix in physical coordinates.

$$C={\Psi }^T{}^{-1} {\widehat{C}}_{mod} {\Psi }^{-1}$$

(35)

The diagonalization of most of the matrices inside the equation of motion allows us to write singularly the equations of motion in modal coordinates (Equation (36)). The only term that is not diagonalized is $\widehat{D}\left(t\right)$ being a non-symmetric matrix that was not involved in the modal analysis. As a consequence, in the n_th equation of motion in modal coordinates $\widehat{D}\left(t\right)$ must be expressed as sum of the products of the elements ${\widehat{D}}_{nr}\left(t\right)$ (elements of the matrix $\widehat{D}\left(t\right)$ at n_th row and r_th column) times the r_th modal displacement $u_r$.

$$\ddot{u_n}+{\widehat{c}}_n\dot{u_n}+{\omega }_n^2 u_n+{\displaystyle {\displaystyle \sum }_{r=1}^N}\left\{{\widehat{D}}_{nr}\left(t\right) u_r\right\}={\widehat{P}}_n\left(t\right), n=1÷N$$

(36)

The Equation (36) represents the useful equation for the development of MMTS by means of which it is possible to compute the modal response $u_n$, in the frequency domain, of the generic n_th mode shape subjected to a modal force ${\widehat{P}}_n\left(t\right)$. The development of MMTS will be presented in the next Paragraph 6.2 starting from the single n_th equation of motion in modal coordinates (Equation (35)).

6.2. Forced Response Computation Using MMTS

MMTS operates by substituting the independent time variable $t$ with the time scales $t_0$, $t_1$, $t_2$, …, where $t_0={\epsilon }^{0}t$, $t_1={\epsilon }^{1}t$ and $t_2={\epsilon }^{2}t$ and ε is the scale factor which describes the time scale. The relation between the original time variable and the new time scales is expressed in Equation (37). Here, the series approximating the old variable $t$ is truncated at the second order (power ${\epsilon }^2$):

$$t=t_0+t_1+o\left({\epsilon }^2\right)=t+\epsilon t+o\left({\epsilon }^2\right)=t+\tau +o\left({\epsilon }^2\right).$$

(37)

As a consequence, the dependent variable $u_n\left(t\right)$ as well as the derivative operators need to be written (Equations (38)–(40)), in the new time scale variables:

$$u_n=u_{n0}\left(t,\tau \right)+\epsilon u_{n1}\left(t,\tau \right)+o\left({\epsilon }^2\right),n=1÷N;$$

(38)

$$\frac{d}{dt}\Rightarrow \frac{\partial }{\partial t}+\epsilon \frac{\partial }{\partial \tau }+o\left({\epsilon }^2\right),$$

(39)

$$\frac{d^2}{d{t}^2}\Rightarrow \frac{{\partial }^2}{\partial t^2}+2\epsilon \frac{{\partial }^2}{\partial t\partial \tau }+o\left({\epsilon }^2\right).$$

(40)

The solution to the problem requires the proper manipulation of the equations of motion. Recalling Equation (36), the manipulation consists of associating some terms of the equation to the coefficient $\epsilon$ (Equations (41)–(43)). Let us associate the time-variant part of the stiffness matrix, $\widehat{D}\left(t\right)$, the damping matrix ${\widehat{c}}_n$ and the modal force ${\widehat{P}}_n\left(t\right)$ to the scale factor $\epsilon$:

$$\widehat{D}\left(t\right)=\epsilon D\left(t\right),$$

(41)

$${\widehat{c}}_n=\epsilon c_n;$$

(42)

$${\widehat{P}}_n\left(t\right)=\epsilon P_n\left(t\right).$$

(43)

Substituting Equations (41)–(43) into Equation (36), the latter becomes:

$$\ddot{u_n}+\epsilon c_n\dot{u_n}+{\omega }_n^2 u_n+{{\displaystyle \sum }}_{r=1}^N\left\{\epsilon D_{nr}\left(t\right) u_r\right\}=\epsilon P_n\left(t\right), n=1:N.$$

(44)

Once the new equation of motion Equation (44) is defined, the MMTS operates by substituting the new expression of the modal response $u_n\left(t\right)$ (Equation (38)) and the new derivative operators (Equations (39) and (40)) into Equation (44). The resulting equation of motion will be an equation where all the terms are characterized by a multiplying coefficient that is in general a power of the scale factor $\epsilon$ (${\epsilon }^n$). Then, a separation of the terms according to the power of $\epsilon$ is performed, by creating $n$ different equations. The new equations are reported below (Equations (45) and (46)).

• Equation corresponding to the power ${\epsilon }^0$:

  $$\ddot{u_{n0}}+{\omega }_n^2{u}_{n0}=0.$$

  (45)
• Equation corresponding to the power ${\epsilon }^1$:

  $$\ddot{u_{n1}}+{\omega }_n^2{u}_{n1}=-2\frac{{\partial }^2{u}_{n0}}{\partial t\partial \tau }-{{\displaystyle \sum }}_{r=1}^N\left\{D_{nr}\left(t\right) u_{r0}\right\}-c_n\frac{\partial u_{n0}}{\partial t}+P_n\left(t\right).$$

  (46)

The solution of Equation (45) can be written in the general form

$$u_{n0}=A_n\left(\tau \right)e^{i{\omega }_{n}t}+\overline{A_n}\left(\tau \right)e^{-i{\omega }_{n}t}=A_n\left(\tau \right)e^{i{\omega }_{n}t}+CC,$$

(47)

where $A_n\left(\tau \right)$ is a function of the time variable $\tau$. Substituting Equation (47) into Equation (46) and adopting the notation in Equation (48), one obtains Equation (49):

$$\dot{(\dots )}=\frac{\partial (\dots )}{\partial t};{(\dots )}^{\prime }=\frac{\partial (\dots )}{\partial \tau };$$

(48)

$$\begin{array}{c}\ddot{u_{n1}}+{\omega }_n^2{u}_{n1}=-2\left[i{\omega }_n{A}_n^{\prime }{e}^{i{\omega }_{n}t}+CC\right]-{{\displaystyle \sum }}_{r=1}^N\left\{D_{nr}\left(t\right) \left[A_r{e}^{i{\omega }_{r}t}+CC\right]\right\}-\\ c_n\left[i{\omega }_n{A}_n{e}^{i{\omega }_{n}t}+CC\right]+P_n\left(t\right),\end{array}$$

(49)

where the quantity $CC$ represents the complex conjugate of the previous term.

In order to solve the N equations of motion (Equation (49)) in the modal coordinate $u$ it is necessary to develop the terms $D_{nr}\left(t\right)$ and $P_n\left(t\right)$ in their harmonic series:

$$D\left(t\right)=D_1\left(t\right)+D_2\left(t\right)+D_3\left(t\right),$$

(50)

$$\begin{array}{c}D\left(t\right)={\displaystyle {\displaystyle \sum }_{s_1=1}^{\infty }} [D_1^{s_1} e^{i s_1{\Omega }_{1}t}+\overline{D_1^{s_1}}{e}^{-i{s}_1{\Omega }_{1}t}]+{\displaystyle {\displaystyle \sum }_{s_2=1}^{\infty }} [D_2^{s_2} e^{i s_2{\Omega }_2 t}+\overline{D_2^{s_2}}{e}^{-i s_2{\Omega }_{2}t}]\\ +{\displaystyle {\displaystyle \sum }_{s_3=1}^{\infty }} [D_3^{s_3} e^{i s_3{\Omega }_{3}t}+\overline{D_3^{s_3}}{e}^{-i s_3{\Omega }_{3}t}]\end{array}$$

(51)

$$P\left(t\right)=P_1\left(t\right)+P_2\left(t\right),$$

(52)

$$\begin{array}{c}P\left(t\right)=P_{1,0}+P_{2,0}+{{\displaystyle \sum }}_{k_1=1}^{\infty } [P_1^{k_1} e^{i{k}_1{\Omega }_{1}t}+\overline{P_1^{k_1}}{e}^{-i k_1{\Omega }_{1}t}]\\ +{{\displaystyle \sum }}_{k_2=1}^{\infty } [P_2^{k_2} e^{i{k}_2{\Omega }_{2}t}+\overline{P_2^{k_2}}{e}^{-i k_2{\Omega }_{2}t}].\end{array}$$

(53)

It is worth to remind that the fundamental frequency ${\Omega }_1$ is the speed of G₁, the fundamental frequency ${\Omega }_2$ is the speed of G₂, linked to ${\Omega }_1$ through the gear ratio $\mathsf{\eta }=\frac{Z1}{Z2}$ and the fundamental frequency ${\Omega }_3$ is the frequency whereby a generic pair of teeth meshes (see Equations (11)–(13)). Substituting the expressions Equation (51) and Equation (53) into Equation (49), one obtains the extended p_th equation of motion in modal coordinates with all the time-variant parameters developed inside (see Appendix A, Equation (A2)). This equation contains all the harmonics of the excitation force, but they can be treated singularly by computing the forced response to each harmonic component of the force. Finally, a sum of all the response contributions can be performed, thanks to the superimposition effect principle, so to compute the overall multi-harmonic response. Thus, the following discussion analyzes the forced response to the generic k_th harmonic component of $P_1\left(t\right)$ acting on the teeth of G₁. The same approach is used to compute the forced response to the second set of harmonics $P_2\left(t\right)$ acting on the nodes of G₂, but here they are not treated for sake of clarity. Let us consider the following equation (Equation (54)) of the p_th modal equation of the system, where only the generic k_th harmonic component of the force function $P_1\left(t\right)$ is considered:

$$\begin{array}{c}\ddot{u_{p1}}+{\omega }_p^2{u}_{p1}=-2i{\omega }_p{A}_p^{\prime }{e}^{i{\omega }_{p}t}\\ -{{\displaystyle \sum }}_{r=1}^N{{\displaystyle \sum }}_{s_1=1}^{\infty }[ D_1{}_{pr}^{s_1} A_r e^{i\left(s_1{\Omega }_1 + {\omega }_r\right)t}+ D_1{}_{pr}^{s_1} \overline{A_r} e^{i\left(s_1{\Omega }_1 - {\omega }_r\right)t} ]\\ -{{\displaystyle \sum }}_{r=1}^N{{\displaystyle \sum }}_{s_2=1}^{\infty }[ D_2{}_{pr}^{s_2} A_r e^{i\left(s_2\mathsf{\eta }{\Omega }_1 + {\omega }_r\right)t}+ D_2{}_{pr}^{s_2} \overline{A_r} e^{i\left(s_2\mathsf{\eta }{\Omega }_1 - {\omega }_r\right)t} ]\\ -{{\displaystyle \sum }}_{r=1}^N{{\displaystyle \sum }}_{s_3=1}^{\infty }[ D_3{}_{pr}^{s_3} A_r e^{i\left(s_3\frac{{\Omega }_1}{n_{T_1}}+{\omega }_r\right)t}+ D_3{}_{pr}^{s_3} \overline{A_r} e^{i\left(s_3\frac{{\Omega }_1}{n_{T_1}}-{\omega }_r\right)t} ]\\ -i{c}_p{\omega }_p{A}_p{e}^{i{\omega }_{p}t}+P_1{}_p^{k_1} e^{i k_1{\Omega }_{1}t}+CC.\end{array}$$

(54)

It is now possible to investigate and remove the unwanted secular terms, inside the latter equation. The elimination of secular terms represents a solvability condition for the solution of the problem, because of the additional freedom introduced with the new independent variables. In order to eliminate secular terms, the resonant terms of each equation need to be forced to zero. The discussion upon secular terms research and elimination is faced more in detail in the Appendix B. Two types of resonant terms can be distinguished inside the equation Equation (54) which can give secular terms: the first type gives “exact” secular terms and they are reported in Equations (55) and (56); the second type can gives nearly secular terms when the excitation frequency $k_1{\Omega }_1$ approaches to ${\omega }_p$, and they are reported in Equations (57)–(59).

$$-2i{\omega }_p{A}_p^{\prime }{e}^{i{\omega }_{p}t},$$

(55)

$$-i{c}_p{\omega }_p{A}_p{e}^{i{\omega }_{p}t},$$

(56)

$$D_1{}_{pr}^{s_1} \overline{A_r} e^{i\left(s_1{\Omega }_1-{\omega }_r\right)t},$$

(57)

$$D_2{}_{pr}^{s_2} \overline{A_r} e^{i\left(s_2 \mathsf{\eta }{\Omega }_1-{\omega }_r\right)t},$$

(58)

$$D_3{}_{pr}^{s_3} \overline{A_r} e^{i\left(s_3 {\frac{{\Omega }_1}{n_T}}_1-{\omega }_r\right)t}.$$

(59)

Since we are interested in the computation of the p_th modal response, it is convenient to introduce an auxiliary frequency variable $\sigma$ to express the neighborhood of the excitation frequency $k_1{\Omega }_1$ to the p_th natural frequency ${\omega }_p$:

$$k_1{\Omega }_1={\omega }_p+\epsilon \sigma ;$$

(60)

$${\Omega }_1=\frac{{\omega }_p}{k_1}+\epsilon \frac{\sigma }{k_1}.$$

(61)

By properly substituting the frequency variable $\sigma$ and performing secular terms elimination by equating to zero the sum of all the possible secular terms, in their critical conditions (see Appendix B for a detailed development), one obtains the following equation in the unknown $A_p$, which is the amplitude of the p_th modal response $u_{p0}$:

$$\begin{array}{c}-2i{\omega }_p{A}_p^{\prime }-D_1{}_{pp}^{\left(2k_1\right)} \overline{A_p} e^{i 2 \sigma \tau } -D_2{}_{pp}^{\left( 2\mathsf{\eta }{k}_{1 }\right)} \overline{A_p} e^{i 2\sigma \tau }-D_3{}_{pp}^{\left( 2 n_{T_1}{k}_1\right)} \overline{A_p} e^{i 2\sigma \tau }\\ -i{c}_p{\omega }_p{A}_p+P_1{}_p^{k_1} e^{i \sigma \tau }=0.\end{array}$$

(62)

Now, let

$$A_p=a_p e^{i\sigma \tau }.$$

(63)

Substituting Equation (63) into Equation (62), it follows

$$-2i{\omega }_p\left(a_p^{\prime }+i\sigma a_p\right)e^{i\sigma \tau }-\mathcal{D} \overline{a_p} e^{i\sigma \tau }-i{c}_p{\omega }_p{a}_p e^{i\sigma \tau }+P_1{}_p^{k_1} e^{i \sigma \tau }=0,$$

(64)

where

$$\mathcal{D}=D_1{}_{pp}^{\left(2 k_1\right)}+D_2{}_{pp}^{\left( 2 \mathsf{\eta }{k}_{1 }\right)}+D_3{}_{pp}^{\left( 2 n_{T_1}{k}_1\right)}.$$

(65)

In order to have a steady-state solution $a_p^{\prime }=\frac{\partial a_p}{\partial \tau }$ has to be null. By eliminating the common term $e^{i\sigma \tau }$, the equation Equation (64) becomes:

$$2\sigma {\omega }_p{a}_p-\mathcal{D} \overline{a_p} -i{c}_p{\omega }_p{a}_p +P_1{}_p^{k_1}=0.$$

(66)

From Equation (66) it is possible to derive analytically an expression of $a_p$ as a function of the frequency variable $\sigma$. As a consequence, the analytical solution of the modal response $u_{p0}$ (Equation (67)) due to the k_th harmonic of the excitation force $P_1\left(t\right)$ is derived.

$$u_{p0}=A_p\left(\tau \right)e^{i{\omega }_{p}t}+CC=a_p\left(\sigma \right) e^{i \left({\omega }_p+\sigma \epsilon \right) t }+CC.$$

(67)

Since $a_p$ is a complex quantity, it is convenient to express $a_p$ according to its real and imaginary parts (Equations (68)–(70)). The analytical expression of the real and imaginary parts is derived as a function of $\sigma$:

$$a_p=a_R+i a_I;$$

(68)

$$a_R=\frac{{\mathcal{P}}_I\left({\mathcal{D}}_I-c_p{\omega }_p\right)+{\mathcal{P}}_R\left(2\sigma {\omega }_p+{\mathcal{D}}_R\right)}{{\omega }_p^2\left(c_p^2+4{\sigma }^2\right)-{\mathcal{D}}_R^2-{\mathcal{D}}_I^2},$$

(69)

$$a_I=\frac{a_R\left({\mathcal{D}}_I+c_p{\omega }_p\right)-{\mathcal{P}}_I}{2\sigma {\omega }_p+{\mathcal{D}}_R},$$

(70)

where:

$$P_1{}_p^{k_1}={\mathcal{P}}_R+i {\mathcal{P}}_I,$$

(71)

$$\mathcal{D}={\mathcal{D}}_R+i {\mathcal{D}}_I,$$

(72)

The latter expressions (Equations (68)–(72)) allow the computation of the modal response of the p_th mode shape in the frequency domain. Each mode shape, connected to a specific nodal diameter of the system, is excited in resonance by some EO of the mesh force, according to the law reported in the following equation Equation (73), from gear dynamics theory [9]:

$$EO=mZ\pm ND , m=1,2,\dots .$$

(73)

Thus, the construction of the multi-harmonic forced response is computed by considering, one by one, each EO of the mesh force, associating it to the mode shapes, which are described by the relative nodal diameter ND, excited by the selected EO and computing the modal responses. Once the modal responses in the frequency domain are computed for all the EO, they are transformed through DMT (Equation (30)), passing from modal coordinates to physical coordinates, and the forced response of a given node is developed in the time domain (in our case the nodes of the teeth of the gears) as the sum of all the mode shapes contributions (i.e., DMT). The validity of the superimposition effect principle (as explained in the Introduction) allows the summation of all the responses due to the different EO in the time domain. The result will be a multi-harmonic response, developed in the time domain. Since the response is not described by a single harmonic it is not possible to acquire the amplitude of the time domain response. Thus, the latter will be expressed by acquiring the peak-to-peak measure of the time domain trend as function of a reference frequency (which can be the speed of G₁ or the speed of G₂ as well) defining uniquely the excitations (both parametric and external) of the overall system. As a matter of fact, setting a certain speed for G₁ means also setting the speed of G₂ since the speed of the gears are linked by the gear ratio. As consequence, the mesh stiffness and mesh force, representing the parametric and the force excitations respectively, directly depend on the speed of the two gears. Finally, the reference frequency describes uniquely the operational conditions of the overall system. In the next section an example of forced response is computed on a dedicated test case and a comparison with Direct Time Integration (DTI) method is made to validate the MMTS methodology, developed into this paper.

7. Forced Response Computed on Test-Cases

In this Section, a study of the forced response of a test case is presented. The aim is to show when MMTS is applied for such applications and why it is convenient to use it. Such a problem can be studied, on the other hand, by Direct Time Integration (DTI) of the equations of motion but this is a very time-consuming method which makes difficult a detailed study of the forced response over a wide range of operational frequency. Anyway, here DTI is used to validate the methodology that is developed in this paper. More in detail, given a certain system model, characterized by given mass, stiffness and damping matrices, two parallel studies are developed on the system, applying MMTS and DTI respectively. The validation of the MMTS methodology is made by comparing the peak-to-peak (P2P) measures of the multi-harmonic response developed in the time domain, calculated with the two methods. Here, the P2P result is plotted against a reference frequency, which is decided at the beginning of the calculation and defines uniquely all the excitations of the system (both parametric and external excitation). The reference frequency chosen for the P2P plots is the speed of Gear-1, ${\Omega }_1$. As a matter of fact, the speed of Gear-2, ${\Omega }_2$, is directly connected to ${\Omega }_1$ through the gear ratio η. Then, all the parametric and external excitations are directly defined. Through the test-case analysis, the dynamic coupling phenomenon is investigated, by remarking its causes and consequences. It is demonstrated that the dynamic coupling, caused by the presence of a time-variant mesh stiffness, leads to a nodal coupling of certain nodal diameters of the meshing gears. To clearly note the phenomenon, a specific test case is built. In more detail, the example which has been used several times in this paper is considered. That is the case of a couple of gears (G₁ and G₂) having respectively Z₁ = 10 teeth and Z₂ = 20 teeth. The system model of each gear, as it was introduced at the beginning of the paper (Section 2) is constituted by two nodes per sector of the gear (the number of sectors is equal to the number of teeth). As consequence, each gear has number of dof equal to twice the number of its teeth. Being the dimension of the gear model twice the number of the sectors, two modal families of natural frequencies can be derived by performing modal analysis of both the gears considering them as stand-alone components. In Figure 9a,b the frequency vs. nodal diameter diagram of the two gears (considered as stand-alone components) is reported, given the mechanical characteristics of the gears shown in table

Table 1. Mechanical characteristics of G1 model.

GEAR-1
Z1	10	Teeth
ND	[0,...,5]
Mechanical characteristics
mb	0.2	(kg)
mc	1	(kg)
ka	107	(N/m)
kb	107	(N/m)
kc	108	(N/m)
damping ratio	0.005	(-)

and

Table 2. Mechanical characteristics of G2 model.

GEAR-2
Z2	20	Teeth
ND	[0,...,10]
Mechanical characteristics
mb	0.2	(kg)
mc	2	(kg)
ka	107	(N/m)
kb	107	(N/m)
kc	108	(N/m)
damping ratio	0.005	(-)

.

The gears are coupled by means of a mesh stiffness of the same type described in Section 3. So, it assumes for the n_th pair of meshing teeth a constant value equal to k_t when the pair is in contact, it assumes a null value when it is not. In this test-case k_t is equal to 10⁶ N/m. The mesh stiffness causes a modal coupling between some modes characterized by specific nodal diameters of the two gears respectively. It is good to remember that the mesh stiffness does not have a remarkable influence on the natural frequencies, which remains practically equal to those of the two gears considered as stand-alone components, even though it may have an important influence on the dynamic response of the overall system. In such a case the numbers of teeth Z₁ and Z₂ has a strong relation with each other. This condition emphasizes the nodal coupling between specific nodal diameters of the gears. It is worth to remind that such a system represents an unusual system because, in practice, gear systems are never designed with such numbers of teeth to avoid the same couple of teeth meshes too often. Nevertheless, this choice, which does not affect the validity of the methodology, aims to boost the effect of the analyzed phenomenon of the dynamic coupling so to better understand it. In Figure 10 an example of nodal coupling is reported for the system under analysis. That is the case of a nodal coupling between the nodal diameter ND-5 of G₁ and the nodal diameter ND-10 of G₂. By recalling the vector of the physical coordinates of the system defined in Section 2 (Equation (1)), the mode shape shown in Figure 10 contains both the nodal diameters. It means that an excitation of ND-5 of G₁ affects the vibration of G₂ which will vibrate at the same natural frequency with a ND-10 shape. It is important to remark that the ND-5 of G₁ (considered as stand-alone component) has a natural frequency ${\omega }_5=1118 \mathrm{Hz}$ (see Figure 9a). When the G₁ is coupled to G₂ by means of the mesh stiffness, the ND-5 of G₁ still has natural frequency ${\omega }_5$, but the mode shape of the coupled system associated to that natural frequency shows an ND-5 mode shape coupled to an ND-10 of G₂. In other words, it is numerically demonstrated that a mesh stiffness with such a value of k_t causes the coupling between the nodal diameters of the gears without changing remarkably the natural frequencies with respect to those of the gears (considered as stand-alone components). The choice to keep a value of k_t that does not cause a remarkable change in the natural frequencies is a reasonable assumption that verifies what is experimentally found in the industrial applications. As a matter of fact, real test cases are characterized by mode shapes showing modal couplings between nodal diameters at given natural frequencies, which are practically the same of those of the stand-alone gears. Thus, the test case under exam aims at simulating a real coupled system where the dynamic characteristics of the mode shapes and natural frequencies remain practically unchanged.

As it was described in Section 5, the external force acting on the system is a mesh force which travels from one tooth to another one with the mesh stiffness. In other words, a force of value F_m acts on the teeth nodes (with same value but with opposite direction) of a specific n_th teeth pair when it is in contact. So, as for the mesh stiffness, the Fourier series of the mesh force is studied, as described in Section 5. In Figure 11a,b the harmonic content (or Engine Orders, EO) of the forces, acting on the two gears respectively, is shown in terms of amplitudes of the various EO.

As anticipated in Section 6.2 (Equation (73)), a harmonic index EO of the travelling force excites mode shapes characterized by a specific ND. Since the ND under analysis for G₁ is 5, the EO excitation that have been selected are: 5, 15, 25, 35. MMTS allows the computation of the modal response of the mode shape of interest due to the selected EO. The forced response in the physical coordinates is easily derived through Direct Modal Transformation (DMT). Here, the forced response (expressed as the P2P measure of the multi-harmonic response developed in the time domain) of the teeth of the two gears is computed, in a given operational speed range (the reference speed is the speed of G₁) where the excitation of the mode shape in Figure 8 occurs. What is expected is to see a resonance of the G₁ due to the excitation of the ND-5 by some EO of the mesh force and an “induced” resonance of the G₂ due to the action of the latter EO exciting the G₁. The fact that the second resonance is induced by the first one is demonstrated by the fact that no excitation of the ND of G₂ should be present for that operational frequency conditions. As a matter of fact, by looking at the Campbell diagrams of the gears, you can note that for G₁ (Figure 12a) the involved EO crosses the natural frequency line in that operational speed range, while for G₂ (Figure 12b) no crossing of the natural frequency lines occurs by the involved EO. Thus, the conclusion is that the second resonance on G₂ is caused by the first one on G₁.

The P2P measure of the multi-harmonic response of the two gears computed using MMTS is reported in Figure 13. Here, also the P2P measure computed by means of DTI is shown in order to make a comparison between the two results. It is worthy to note that there is a big difference in terms of computational time for the construction of the forced response using the two methodologies (MMTS and DTI). In more detail, a DTI study can take some hours, and the computational time can increase considerably as the number of dof of the system increases as well as the resolution of the operational speed range and the integration time interval increase. As consequence, the computational efforts can be practically unsustainable for systems with a large number of dof and very low damping ratios whereby DTI must integrate for larger time interval in order to reach a steady-state response, where the transient part is completely extinguished (i.e., a real test-case of gear system where damping ratio is lower than 0.1%). On the other hand, MMTS, operating in the frequency domain, results in a very fast calculations of the forced response which can take few minutes and then it allows to select the EO of interest which have an influence in a given operational speed range, neglecting the minor effects of other EO and so reducing considerably the computational efforts. In both the analysis, two resonance peaks are visible in a given operational speed range (speed of G₁ from 210 Hz to 240 Hz). The blue curves are the resonances of G₁ computed respectively with MMTS and DTI. The same for the red curves. The resonance of G₂ is induced by the resonance of G₁, as anticipated before. By looking at the figures Figure 14a,b showing respectively the FFT of the time domain responses, computed through DTI, of the gears (respectively shown in the figures Figure 15a,b), it is easy to note that the main harmonic component of the multi-harmonic response in both cases is exactly the natural frequency of the mode shape shown in Figure 10 which couples ND5 of G₁ to ND10 of G₂. This represents an additional proof that the resonance of G₂ is directly induced by the excitation of that single mode shape by the mesh force on the G₁. This is a clear example of dynamic coupling between two meshing gears and MMTS allowed to forecast the resonance of G₂ induced by the excitation of the ND of G₁ and this could not be possible if you have considered the gears as stand-alone components. In that case, no interaction between the ND of the gears can be studied.

8. Conclusions

The objective of this paper is to investigate the mutual interactions (dynamic coupling) which affect the response of a couple of meshing gears, by developing a methodology able to compute easily, with limited computational efforts, the forced response of the gears, without loosing the generality and complexity of the system. As a matter of fact, here the gears are considered as compliant bodies. As opposed to other methodologies which were developed in the past, whereby the assumption of the gears as rigid bodies needed to be supported by the introduction of the transmission error to simulate the compliance of the gears. Here, the challenge is to couple two compliant gears (whose dynamic characteristics are automatically included) and to investigate how the dynamics of a gear interacts with the other one when phenomena of mesh stiffness fluctuations occur. In addition to that, the methodology provides guidelines for an analytical solution to the problem, allowing the researcher to compute the forced response of a complex system undergoing both a parametric and external excitation. The choice of MMTS for the mathematical solution of the linear time-variant problem is addressed to its capability to provide an analytical solution to the problem, by strongly simplifying it. This is a great advantage for applications like gear coupling where the simplicity of the methodology (MMTS) compensates for the complexity of the system and allows for the analysis of the behavior of the gears in a considerably wide range of operation.