An Analytical Approach to Gear Mesh Dynamics for the Sustainable Design of Agricultural Machinery Drive Systems

1. Introduction

Gear meshing is a complex process in which precision of calculations, design, machining, and assembly stands extremely important. The proper mating between gears not only assures lesser chance of failures of other elements of the machine, such as bearings [1]. It also extends the working life of the gears themselves and improves safety of operation [2]. Thus, throughout the years a visible effort is put into analysing meshing performance in different types of gear arrangements [3,4,5]. Optimal distribution of forces and lubrication decreases overall levels of different loads emerging on the surfaces and edges of the gears which effect in wear and pitting [6]. In some specific applications, new types of teeth are introduced to improve general efficiency, such as concave-convex ball teeth for hydraulic pumps [7], as well as novel drives [8,9] or arrangements [10,11]. Multiple methods are present in the field of gear meshing analysis, such as tooth contact analysis or finite element analysis. Interesting modifications to them are presented continuously [12], and new methods are introduced simultaneously [13,14] resulting in better performance achievable.

In comparison with passenger cars, the users of agricultural equipment are subjected to completely different spectrum of vibroacoustic (VA) signals. Two main reasons for that might be highlighted, namely the working conditions and design requirements, which coincide strictly with each other. To be able to undertake various tasks such as hauling, plowing or harvesting in various weather and on adverse surfaces, a proper engine torque and chassis/instrumentation durability is required, yet unnecessary complexity should be avoided. This results in high rigidity of underbody, most often with complete lack or very simplistic suspension, combined with big-diameter tires made of hard rubber with aggressive type of tread. To conclude, big unsprung mass, high chassis stiffness, big-displacement compression ignition engines, multitude of geared mechanisms and types of ground excitations create extremely hazardous work environment, as basically only the seat and it’s suspension isolates the operator from harmful vibroacoustic signals [15,16,17,18,19]. But the non-direct influence on the environment and surroundings should not be omitted [20,21]. Formerly rural areas close to the big cities are being urbanized with growing pace, hence neighbourhood of big estates of apartments or single-family houses with farmlands and farms is not unusual sighting anymore. As the farm works follow specific schemes which cannot be disturbed (in terms of type of the equipment, season and time of its use) it is crucial to reduce all types of pollution, including vibroacoustic, to ensure some kind of balance between a job that has to be done, and living comfort of the residents. Lately, reduction of engine noises whether by dampening internal combustion engines or switching to electric motors introduced new challenges, as the noise generated by transmissions and gears became more evident for operators and surroundings [22,23,24,25].

As it was clearly stated above, vibrations are an inherent part of agricultural machinery operation which affect the operators and environment [15]. To complete the picture, an influence of the vibrations on the machine itself should also be highlighted [17,26]. The presence of vibrations can indicate both normal operation of the drive system [27,28] and the beginning of a failure [29,30]. Therefore, it is crucial to be able to recognize and interpret the characteristics of generated vibrations. One of the main sources of vibrations in agricultural machinery are gearboxes [24]. Their operation is associated with the occurrence of parametric vibrations [28,31,32,33,34,35]. These are vibrations whose frequency changes over time. In addition to gearboxes, forced vibrations can be generated by other elements of the drive system, such as: driveshafts (mass imbalance, bending, or damage to the shafts lead to the generation of centrifugal forces, which cause vibrations) [36,37,38], bearings (bearing damage, such as wear of the raceway or rolling elements, causes increased noise and vibrations) [1,37]. Analysis of vibrations generated by the drive system of an agricultural machine allows for early detection of potential failures and prevention of more serious damage [39,40,41]. Thanks to appropriate sensors and analysis (supported by analytical solutions), it is possible to track changes in vibration characteristics and, based on this, make decisions regarding the need for inspection or repair [27,42]. Vibrations in the drive systems of agricultural machinery are a complex phenomenon resulting from many factors. Understanding their causes and characteristics is crucial for ensuring the reliability and durability of machines [43]. Gearboxes are the main source of parametric vibrations, while other elements of the drive system generate mainly forced vibrations [31,44].

Regardless of their type and source, noise and vibrations generated by gearboxes pose a significant challenge to sustainable development. Their negative impact on the environment and human health is undeniable. High noise levels contribute to stress, sleep disturbances, and prolonged exposure can lead to hearing loss [18,21,45,46,47]. Reducing noise and vibrations generated by gearboxes is of great importance for the quality of people’s lives. Silence and peace are essential for the proper functioning of the human body. Reducing noise in the workplace improves working conditions, increases efficiency, and reduces the risk of occupational diseases [48,49,50]. In addition, noise reduction in the residential environment contributes to improving people’s mental and physical health. Vibrations, in turn, can damage machinery and structures, as well as affect work comfort [36]. There are many ways to reduce noise and vibrations generated by gearboxes. It’s possible to apply damping materials [46,51], vibration isolation [2,52,53], structural optimization [54,55,56] or lubrication [55,57,58]. However, there are many scenarios when it’s not possible to elimite VA energy emission. Vibration energy can be treated as an aspect of sustainable development. Gearbox vibrations, although undesirable due to noise, can be a valuable source of energy. Thanks to technological development, it is possible to harvest vibration energy and convert it into electrical energy. This opens up new perspectives for sustainable development, enabling a reduction in the consumption of conventional energy and limiting the emission of harmful substances [24,59]. The problem of noise and vibrations generated by gearboxes requires a comprehensive approach. On the one hand, efforts should be made to reduce them in order to improve the quality of the environment and human health [18,20,45]. On the other hand, vibrations can be a valuable source of energy, which opens up new opportunities for sustainable development. Therefore, research on ways to reduce noise and vibrations and to effectively utilize vibration energy is extremely important [21,60].

Therefore, the problem of vibration reduction (VR) and simultaneous energy harvesting (EH) is becoming increasingly urgent in the context of technological development. Many stuctures are exposed to vibrations that can lead to failures and energy losses. Given that, the search for effective methods of vibration reduction and efficient energy harvesting from this phenomenon is the subject of intensive research [52,61]. There are many vibration reduction methods, such as passive (vibration dampers, vibration isolation, high damping materials) [62], active (control systems that generate forces opposite to the forces exciting vibrations) [41], and hybrid (a combination of passive and active methods) [63]. The choice of a suitable method depends on many factors, such as: vibration frequency, vibration amplitude, operating environment, and cost and mass requirements. Energy harvesting from vibrations involves converting the mechanical energy of vibrations into electrical energy. Few types of phenomena can be useful in order to reach this goal: piezoelectricity (the use of materials that generate an electric voltage when deformed), electromagnetism (the use of the electromagnetic induction phenomenon), or electrostatics (the use of electric charges). Analytical solutions can significantly facilitate the process of selecting vibration reduction and energy harvesting methods [61]. Thanks to them, it is possible to quickly assess the effectiveness of different solutions, optimize system parameters or minimize the number of costly experiments [57]. It allows to state that vibration reduction and energy harvesting is very complex and requires interdisciplinary knowledge. Thanks to the development of new materials, computational methods, and sensor technologies, it is possible to create increasingly advanced solutions [2]. Analytical solutions play a key role in this process, enabling the rapid and efficient design of vibration reduction and energy harvesting systems.

Agricultural machines, including tractors and harvesters, produce complex vibration profiles. These vibrations originate from both internal machine dynamics, such as the drive train, and external interactions, including terrain variability [26,64]. The primary sources of these vibrations are internal combustion engines, gearboxes, and other drive train components, which generate low-frequency oscillations that propagate through various mechanical components. The attachement of supplementary devices further amplifies these vibrations, particularly through three-point linkage systems, exacerbating operator exposure during tasks such as tillage [56,65,66]. The operational activities such as tilling, harvesting, and soil preparation amplify these vibrations through dynamic interactions with the terrain, and towed or mounted implements [47,50,53,67].

Vibration propagation in agricultural machinery represents a complex phenomenon rooted in the interactions among various mechanical and structural components [64]. Vibrations generated in core components like the engine and gearbox propagate through the drive shafts, framework, and connected parts. It can be termed as transmission through structural elements (machinery components). The gearboxes, in particular, are critical nodes of vibration due to their dynamic mechanical interactions. Studies indicate that gear meshing imperfections and load variations in gearboxes can significantly amplify vibration transmission.

Many research text delves into the propagation mechanisms, impacts, and mitigation strategies, focusing on the intricate role of drive systems, especially gearboxes, in these processes. The transmission of vibrations through the drive train and gearboxes hastens wear on bearings, seals, and interconnected components. Research highlights that repeated vibratory stresses lead to microstructural fatigue, precipitating premature failures (accelerated component wear). Prolonged vibration exposure can induce structural degradations (cracks) and deformations in structural elements. This mechanical degradation is particularly pronounced in the connections between drive systems and the chassis, where stress concentrations are the highest.

Designing agricultural machinery focusing on vibration reduction differs significantly from designing cars or other vehicles. Unlike conventional vehicles, agricultural machines are exposed to a diverse array of operational conditions, demanding robust structural frameworks. These frames, while providing essential strength and support, often contribute to the efficient transmission of vibrations from the source, such as the drive train, to the entire machinery body.

Vibration in agricultural machinery propagates through the entire structure. This process is a result of complex interactions between various machine components [64,68]. Understanding the nature and transmission of these vibrations is crucial for mitigating their impact on operators’ health and ensuring optimal performance [46,66,69]. The operational environment of agricultural machinery subjects operators to significant whole-body vibration (WBV). Vibrations predominantly occur in the frequency range of 1–10 Hz, aligning with the resonance frequencies of human organs, which can lead to amplified mechanical responses and potential health issues [49,59]. In particular, the vertical axis (z-axis) vibrations tend to be most pronounced and pose significant ergonomic challenges with significant health implications when combined [47,53].

The European Directive 2002/44/EC sets exposure action and limit values for WBV, indicating the importance of adhering to safety thresholds to prevent health risks [45]. Research has shown that many current agricultural practices exceed these thresholds, highlighting the need for enhanced compliance measures and operator training. Prolonged exposure to WBV has well-documented adverse effects, including musculoskeletal disorders, spinal degeneration, and lower back pain [46,48]. Moreover reduces performance and chronic health issues such as cardiovascular strain [67,69]. Studies indicate that the resonance frequencies between 4–7 Hz, which coincide with human body responses, exacerbate these conditions. Research on seat and cabin suspension systems has shown that without appropriate damping measures, the transmission of vibrations from the drive train and other machine components can significantly compromise operator health [26,45,46,59,70].

The ISO 2631 standard specifies vibration thresholds to protect against these effects. However, research indicates that even with standard-compliant equipment, drivers can experience discomfort, particularly at frequencies resonant with the human body’s natural response (4–7 Hz), leading to fatigue and drowsiness, which could impair driver performance linked to cognitive effects [26,51,56,69].

Vibrations are not confined to solid structures. They can also generate acoustic waves, which propagate through the air, manifesting as noise (acoustic wave propagation). This acoustic transmission is particularly notable in older machinery models where cabin insulation is inadequate. The repercussions of vibration propagation extend beyond machine components, affecting the machinery’s operational integrity and the well-being of operators [48,50]. Vibrations convert into noise, creating an auditory disturbance both for the operator and the surrounding environment [45,49]. Elevated noise levels, often exceeding permissible exposure limits, can impair hearing and lead to long-term occupational health risks. That is the key impacts include.

The insights mentioned above allow us to state that understanding the processes of vibration propagation in agricultural machinery is crucial for designing safe and comfortable machines. The use of appropriate materials, optimization of the design, and the application of vibration damping systems are just some of the methods that can significantly reduce vibration levels and improve operator working conditions.

Understanding vibration propagation in agricultural machinery, particularly through drive trains and gearboxes, is essential to designing equipment that ensures both machine longevity and operator safety. The integration of advanced damping systems and innovative gearbox designs holds promise for reducing vibration levels. In this respect, it is particularly important to reduce the impact of gearbox vibration on the wear and tear of bearings, seals and connected parts. It must be borne in mind that gearboxes are critical vibration nodes due to their dynamic internal interactions resulting from the nature of their operation. Gears generate parametric vibrations as a result of time-varying gear meshing. It is therefore extremely important to carry out research into the effects of gear meshing on load variations in gearboxes, which contribute to significantly amplifying the propagation of vibrations to the entire machine structure. Continuous research emphasizing vibration reduction and effective propagation control is part of the policy of sustainable design of drive systems for agricultural machinery. This has a direct bearing on increasing the durability of these machines and operator comfort, fostering a safer and more efficient farming practices, and more efficient agricultural work environment. These are undoubtedly key aspects of sustainable design of agricultural machinery.

Modelling of gear dynamics is important field of theoretical engineering (designing) [5,71]. Numerous models allow to simulate the operation of gearboxes under various conditions, analyze their strength, and predict potential failures [63]. Important case is parametric vibration, which is characterized by time-varying system parameters such as stiffness and damping [28,31,34,35].

Parametric vibrations differ significantly from other types of vibrations because their characteristics depend not only on the properties of the system but also on its operating conditions. Understanding the mechanisms of the generation and propagation of these vibrations is crucial for ensuring the reliability and durability of gearboxes. Parametric vibrations can lead to resonances, which in turn can cause damage to structural components [38,41].

Applications of dynamic models are as follows:

1. Design: dynamic models allow for the optimization of gearbox designs to minimize vibration and noise [72]. They can be used to select appropriate materials, element geometries, and bearing parameters to ensure a long service life of the device [71,73];

2. Diagnostics: analysis of vibration signals obtained from an operating gearbox allows for the early detection of faults such as tooth wear, shaft clearances, or bearing damage. This makes it possible to plan preventive maintenance and avoid costly downtime [27,30,37,43,72];

3. Reducing downtime: early detection and elimination of the causes of parametric vibrations significantly reduces the number of failures and increases machine availability [39].

Dynamic modeling of gearboxes, with a particular focus on parametric vibrations, is an invaluable tool for engineers. It enables to better understand the complex phenomena occurring in these devices and, consequently, to design more reliable and durable gearboxes. In addition, modeling enables early fault detection, which significantly reduces downtime and optimizes maintenance costs.

Although numerical methods provide very broad set of approaches to analyze the problem in question and are widely used, analytical solutions where applicable are advantageous in many aspects [74] [75]. Firstly, properly defined and assessed formulas ensure repeatability of obtainable results. Secondly, even though analytical solution primarily might be more time-consuming than its numerical counterparts, it shortens the overall analysis time as a simple substitution of different parameter values is sufficient enough to obtain new results. On contrary each change of parameters with numerical approach requires running the simulations once again. Given so, a demand arises for precise computational methods based on the determination of analytical solutions also in the aspect of gear meshing analysis [27,41,76]. Unfortunately, obtaining such a solution for a complex process of gear meshing is difficult, especially in the case of multiple degree of freedom transmissions with high-speed gears, which are subjected to broad changes of rotational speeds, heavy loads, and numerous excitation sources, as every pair of working gears have time varying mesh stiffness (TVMS). Given that, numerical methods are far more popular, but still approaches to obtain hybrid analytical-numerical [77] and fully analytical solutions are regarded. One of the recent works introduced an improved analytical model for calculating mesh stiffness, based on contact mechanics and gear engagement principle, which correctly calculates varying stiffness under non-stationary conditions [52]. The model is based on evaluating the potential energy of the non-uniform cantilever beam consisted of four parts: local contact energy (Hertz contact energy), bending energy, axial compressive energy and shear energy with additional consideration of addendum modification, and was verified experimentally. Similar topic is analyzed in [78], where the empirical ISO (ISOSM) standard model and deformation energy based model (DEM) are compared with finite element method results. It was proved, that due to results being highly approximated in case of ISOSM, it should be applied rather to gears with standard profile, whereas DEM shows much higher accuracy thanks to possibility of including exact geometry of the analysed gearing. Another examples introduce the kinematics and the dynamical effects of centrifugal forces for high-speed transmissions [79], time-varying backlash [80], models for meshing stiffness, load sharing and transmission errors [27,62,81,82], load sharing and friction torque of involute spur and helical gears considering a non-uniform line stiffness and line load [83]. TVMS is also a great tool for analysing and predicting defects in gears, such as spalling, pitting, and cracking, and analytical approaches are also utilized in this aspect [60,84,85]. The variable stiffness raises an interest not only in the case of analysing the dynamics of mating gears [86,87], but also when focusing on the gear shaping process. In [88] the authors introduced a novel TVMS method, which improves contact stiffness model considering processing textures. It might be noticed, that as mostly the solutions introduced throughout the years derive from basic cantilever beam approach, still major modifications of often different sources are needed to properly represent specific behaviour. In the present study, a formal analysis of the possibility to utilize TVMS-based parametric vibration solutions of non-linear equations of motion with Multiple Scale Method (MSM) is conducted, and how this approach could possibly simplify the process of identifing and including different phenomena and geometries in analytical analysis.

The primary objective of this research is to advance the understanding and application of dynamic models for gears within agricultural machinery drive trains by developing analytical solutions. Despite significant advancements in vibration analysis, there is a notable scarcity of comprehensive research that addresses the analytical modeling of gear dynamics, particularly using advanced mathematical techniques such as the Multiple Scale Method. This research aims to bridge this gap by starting with the simplest gear models and progressively scaling up to more complex systems.

The focus of this study is on employing the MSM to derive analytical solutions that enable the determination of vibration frequency ranges and the identification of emitted components. This mathematical approach offers a generalized framework for understanding vibration propagation and can significantly enhance the precision of vibration management strategies. By leveraging these insights, the research aims to contribute to the improved design of agricultural machinery drive trains, emphasizing vibration reduction and efficient propagation control. Ultimately, this work seeks to enhance these machinery’s durability and operator comfort (with particular concern for his welfare), fostering safer and more effective agricultural practices.

2. Dynamic Model of Gear Meshing and Intertooth Force

Gears in normal mesh (where engagement factor ${\epsilon }_{\alpha }$ falls between 1 and 2) experience both single-tooth and double-tooth contact along their engagement path. To understand the overall stiffness of normal mesh, it is first necessary to analyse the stiffness of individual teeth as they come into contact. This can be done analytically for a single pair of teeth (denoted as $k_p\left(t\right)$). Method B from ISO 6336 standard [73] can be used to find the maximum stiffness ($k_{th}^{\prime }$) for this single pair. Equation (1) is then applied for the complete depiction, as described in detail in reference [89].

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k_p={\displaystyle \frac{20k_{th}^{\prime }}{17{\epsilon }_{\alpha }}}\left(0.55+{\displaystyle \frac{1.8t}{{\epsilon }_{\alpha }{t}_z}}-{\displaystyle \frac{1.8t^2}{{\epsilon }_{\alpha }^2{t}_z^2}}\right)\end{array}\end{array}$$

(1)

where:

$k_p$

- stiffness of single pair of teeth;

$k_{th}^{\prime }$

- maximum stiffness;

${\epsilon }_{\alpha }$

- engagement factor;

$t_z$

- period of mesh stiffness variation;

  t

- time.

The equation of motion of gears system that was reduced to a single degree of freedom is given by the equation (2), where time-varying stiffness is given as a Fourier expansion of rectangular signal:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -F+c_g\dot{z}+m_{eq}\ddot{z}+k_g\left({\kappa }_{mesh}\epsilon +1\right)z=0\end{array}\end{array}$$

(2)

where:

  F

- excitation force;

$c_g$

- mesh damping;

$m_{eq}$

- equivalent mass;

$k_g$

- mesh stiffness;

${\kappa }_{mesh}$

- formula representing stiffness variation (Fourier expansion of rectangular signal);

$\epsilon$

- small parameter;

  z

- gear vibration displacement;

$\dot{z}$

- gear vibration velocity;

$\ddot{z}$

- gear vibration acceleration.

In order to obtain a beneficial data from the point of view of mathematical analysis of the mesh stiffness signal the anlytical method of Fourier series expansion is applied. Obtained expression is a result of numerical integration of the mentioned signal. Formulas regarding expansion coefficients are presented below. Additionally the Fast Fourier Transform (FFT) was involved as an efficient tool to determine the the influence of specific harmonics on the signal amplitude:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill k_p={\displaystyle \frac{a_0}{2}}+\sum _{n=1}^{\infty }\left(a_{n}cos\left({\displaystyle \frac{\pi n{k}_p}{T}}\right)+b_{n}sin\left({\displaystyle \frac{\pi n{k}_p}{T}}\right)\right)\end{array}\end{array}$$

(3)

where:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill a_0={\displaystyle \frac{\underset{-T}{\overset{T}{\int }}{k}_{p}dt}{T}}\end{array}\end{array}$$

(4)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill a_n={\displaystyle \frac{2\underset{-T}{\overset{T}{\int }}{k}_{p}cos\left(\frac{\pi nt}{T}\right)dt}{T}}\end{array}\end{array}$$

(5)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill b_n={\displaystyle \frac{2\underset{-T}{\overset{T}{\int }}{k}_{p}sin\left(\frac{\pi nt}{T}\right)dt}{T}}\end{array}\end{array}$$

(6)

From the properties of analyzed function it can be observed that only odd terms of the series will be non-zero. Additionaly, from the theorem about integration of odd functions in the symmetrical domain, which is the case here, it occurs that only cosine terms will be present. In such case, a large number of terms would be required for satisfactory approximation, so the FFT algorithm is used to analyze the singal in the frequency domain to obtain the amplitudes of particular harmonics.

The formula (1) includes the parameter $t_z$, which represents the period of mesh stiffness variation that corresponds to the time of meshing gear teeth. The pair of teeth is in contact at this time. This means that a pair of teeth is on a path of contact followed by the point of contact between two meshing gear teeth. The contact points of the teeth line up with the line-of-action (pressure line) during the rotation of the gear. In general, for involute gears the tooth-to-tooth force is always directed along the same line. That is, the line-of-action is constant – along which the force vector between two meshing gear teeth is directed. It is worth noting the crucially, that the comprehensive mesh stiffness of gear is the elastic response of the comprehensive meshing effect of each pair of teeth in the whole meshing area, which is directly related to the elastic deformation of single tooth, the comprehensive elastic deformation of single pair of teeth and the coincidence degree of gear system itself. For the spur gear system with coincidence degree greater than 1 and less than 2, the gear includes single-tooth meshing area and double-tooth meshing area in one meshing cycle. For such a case, the time course of mesh stiffness in the single-tooth meshing region (a single-tooth pair mesh stiffness) takes the shape of a parabolic function. However, in the double-tooth meshing area, the stiffness is calculated through superposition, treating both pairs as parallel-connected springs with varying characteristics. The method to calculate the combined mesh stiffness of two interacting pairs of teeth is to add the mesh stiffness of two pairs of teeth together. The comprehensive time-varying mesh stiffness causes changes with a constant period, as shown in Figure 1. A higher stiffness value corresponds to double-tooth meshing, while a lower stiffness value corresponds to single-tooth meshing. Here, four different signals are presented:

1. reference - typical parabolic shape;

2. wave-shaped - commonly used to represent the parabolic character of teeth meshing due to the sufficient representation of physical properties;

3. rectangular - rather less popular for gear meshing approximation, worse in accuracy in comparison with wave-shaped;

4. rectagular’s Fourier approximation - an approach proposed in the study.

The FFT spectra based on the time-course (Figure 1) are shown in Figure 2 and Figure 3. A detailed description of determining the overall time-varying mesh stiffness course has been provided in the following paper [90]. All the calculations were made for the data presented in the Table 1, based on research [71].

A series of numerical simulations of the equation (2) was conducted to evaluate the effectiveness of the adopted gear stiffness models (Figure 1). These simulations were made to initially check the correctness of the gear meshing displacements (Figure 4, Figure 5, and Figure 6) for each stiffness model.

Observations of graphs in Figure 4, Figure 5, and Figure 6 allow for comparison of the simulation results obtained based on solution of the gear system, where mesh stiffness is given consecutively as a wave signal (WS), rectangular signal (RS), and its Fourier approximation (AP). This comparison aims to conclude the differences in obtained waveforms. All simulations were conducted for the same period to avoid doubts about the frequencies of obtained signals and allow for pure comparison of waveform shape and amplitude. In general, the amplitudes corresponding to single-tooth and two-teeth contact reach similar values, but bigger discrepancies are clearly visible when referring to the shape of the plots. Qualitatively, the RS and AP waveforms represent the physical behavior of the meshing gears with higher accuracy than WS. Moreover, AP waveform includes signal modulation, enhancing the overall information carried.

The apparent differences in the transient phases are striking when performing a comparative analysis of the RS (Figure 5) model and its Fourier approximation AP (Figure 6). This is due to the use of simplification in the case of approximated waveform (Figure 6). Namely, using Fourier series expansion results in smoother stiffness variation, while locally the susceptibility (mesh displacement) values can be artificially overstated (overloaded) or underestimated. In the case of the WS model (Figure 4), sudden changes result in longer transients, which are evident when the number of pairs of teeth in the contact line is changed. It is impossible to remove excess stiffness in the Fourier approximated model due to the Gibbs phenomenon, which is an immanent feature of the Fourier series of a periodic function with low smoothness (C0 class).

The equivalent meshing force for the pair of gears reduced to a single degree of freedom is given by formula (7):

$$\begin{array}{c}\hfill \begin{array}{c}\hfill q_{eq}={\displaystyle \frac{F_c}{m_{eq}}}={\displaystyle \frac{c_g\dot{z}}{m_{eq}}}+{\displaystyle \frac{k_{g}z}{m_{eq}}}\end{array}\end{array}$$

(7)

where:

$q_{eq}$

- equivalent meshing force;

$F_c$

- meshing force.

Further transformation of the formula (7) allows to obtain the following form of the analyzed equation:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill q_{eq}={\omega }_0^{2}z+2h\dot{z}\end{array}\end{array}$$

(8)

where:

${\omega }_0^2$

- natural frequency;

  h

- damping coefficient.

The time-varying meshing force waveform based on the Fourier approximated model and including damping along the line-of-action is shown in Figure 7.

The apparent transients associated with the stabilisation of the inter-tooth force values are due to the damping occurring in the mesh system. The time course of meshing force changes with a constant period, which correlates with the mesh stiffness (Figure 6). The distribution of the meshing force is characterised by sudden increments of force when one pair of teeth starts to mesh. This corresponds to the temporal moment the meshing force takes its maximum value. A lower meshing force values occurs when two pairs of teeth are in contact along the line-of-action – specifically, from entering the line-of-action to exiting it and sudden increments value when one pair of teeth is in contact.

3. Analytical Solution

General formulation of the dynamic model provided in (2) can be simplified for analytical investigation. Equation (9) is given in more convenient way for further analysis, where natural frequency of the system is substituted with meshing frequency. To ensure that natural frequency will be much higher than meshing frequency a relationship between these two frequencies was selected as ${\omega }_0=16\omega$, hence the appearance of numerical values of 256 and 582 in consecutive parts of the equation. The small parameter in the considered system indicates the elements of the equation that influence the main part of the solution and the perturbative components. Using the Multiple Time Scales Method (MTSM), a perturbation system of the equations of motion was formulated. The resulting system of equations is recursive in nature. It is solved in a sequential way by looking for solutions to successive approximations.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 256{\omega }^{2}z+256\delta {\omega }^2{u}_s\epsilon +256\delta {\omega }^2\epsilon z+582.0{\omega }^2{u}_s\epsilon sin\left(\omega t\right)+582.0{\omega }^2\epsilon zsin\left(\omega t\right)+\ddot{z}=0\end{array}\end{array}$$

(9)

where:

$u_s$

- static deflection;

$\delta$

- detuning parameter;

$\omega$

- meshing frequency.

The solution is approximated by a Taylor polynomial. For the initial conditions, it is assumed that the solution can be approximated by a sum of functions, each of which depends on small parameter. This provides a recursive structure for solving the equations [91].

$$\begin{array}{c}\hfill \begin{array}{c}\hfill z=\epsilon z_1+{\epsilon }^2{z}_2+{\epsilon }^3{z}_3+z_0\end{array}\end{array}$$

(10)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill z=\epsilon z_1\left(t_0,t_1,t_2\right)+{\epsilon }^2{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^3{z}_3\left(t_0,t_1,t_2\right)+z_0\left(t_0,t_1,t_2\right)\end{array}\end{array}$$

(11)

where:

$z_0,z_1,z_2,z_3$

- sequence of solution approximations.

Formulas that represent arguments of approximation (time scales) are as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill t_0=t\end{array}\end{array}$$

(12)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill t_1=t\epsilon \end{array}\end{array}$$

(13)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill t_2=t{\epsilon }^2\end{array}\end{array}$$

(14)

Consecutive derivatives of the functions are calculated and substituted into the equation (9). By grouping the elements with respect to small parameter a set of equations can be formulated. Their solution estimates the approximation functions. At this stage, the solution was approximated for different orders of $\epsilon$. These approximations provide increasing accuracy with respect to the original nonlinear equation.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \dot{z}=\epsilon {\dot{z}}_1+{\epsilon }^2{\dot{z}}_2+{\dot{z}}_0\end{array}\end{array}$$

(15)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \ddot{z}=\epsilon {\ddot{z}}_1+{\epsilon }^2{\ddot{z}}_2+{\ddot{z}}_0\end{array}\end{array}$$

(16)

In the MTSM, the function is approximated in regard to different time scales (denoted as $t_0$, $t_1$, $t_2$, ..., $t_n$). Each time scale represents a different rate of change in the system’s dynamics. The resulting partial functions are substituted back into original function and derivatives are calculated which later can be substituted into equation of motion.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill z=\epsilon z_1\left(t_0,t_1,t_2\right)+{\epsilon }^2{z}_2\left(t_0,t_1,t_2\right)+z_0\left(t_0,t_1,t_2\right)\end{array}\end{array}$$

(17)

$$\begin{gathered} \begin{array}{c}\hfill \begin{array}{c}\hfill \dot{z}={\dot{t}}_0{\displaystyle \frac{d}{d{t}_0}}{z}_0\left(t_0,t_1,t_2\right)+{\dot{t}}_1{\displaystyle \frac{d}{d{t}_1}}{z}_0\left(t_0,t_1,t_2\right)+{\dot{t}}_2{\displaystyle \frac{d}{d{t}_2}}{z}_0\left(t_0,t_1,t_2\right)+\epsilon {\dot{t}}_0{\displaystyle \frac{d}{d{t}_0}}{z}_1\left(t_0,t_1,t_2\right)+\epsilon {\dot{t}}_1{\displaystyle \frac{d}{d{t}_1}}{z}_1\left(t_0,t_1,t_2\right) \\ +\epsilon {\dot{t}}_2{\displaystyle \frac{d}{d{t}_2}}{z}_1\left(t_0,t_1,t_2\right)+{\epsilon }^2{\dot{t}}_0{\displaystyle \frac{d}{d{t}_0}}{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^2{\dot{t}}_1{\displaystyle \frac{d}{d{t}_1}}{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^2{\dot{t}}_2{\displaystyle \frac{d}{d{t}_2}}{z}_2\left(t_0,t_1,t_2\right)\end{array}\end{array} \end{gathered}$$

(18)

$$\begin{array}{cc}\hfill \ddot{z}& ={\dot{t}}_0^2{\displaystyle \frac{d^2}{d{t}_0^2}}{z}_0\left(t_0,t_1,t_2\right)+{\dot{t}}_1^2{\displaystyle \frac{d^2}{d{t}_1^2}}{z}_0\left(t_0,t_1,t_2\right)+{\dot{t}}_2^2{\displaystyle \frac{d^2}{d{t}_2^2}}{z}_0\left(t_0,t_1,t_2\right)+{\ddot{t}}_0{\displaystyle \frac{d}{d{t}_0}}{z}_0\left(t_0,t_1,t_2\right)+{\ddot{t}}_1{\displaystyle \frac{d}{d{t}_1}}{z}_0\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +{\ddot{t}}_2{\displaystyle \frac{d}{d{t}_2}}{z}_0\left(t_0,t_1,t_2\right)+\epsilon {\dot{t}}_0^2{\displaystyle \frac{d^2}{d{t}_0^2}}{z}_1\left(t_0,t_1,t_2\right)+\epsilon {\dot{t}}_1^2{\displaystyle \frac{d^2}{d{t}_1^2}}{z}_1\left(t_0,t_1,t_2\right)+\epsilon {\dot{t}}_2^2{\displaystyle \frac{d^2}{d{t}_2^2}}{z}_1\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +\epsilon {\ddot{t}}_0{\displaystyle \frac{d}{d{t}_0}}{z}_1\left(t_0,t_1,t_2\right)+\epsilon {\ddot{t}}_1{\displaystyle \frac{d}{d{t}_1}}{z}_1\left(t_0,t_1,t_2\right)+\epsilon {\ddot{t}}_2{\displaystyle \frac{d}{d{t}_2}}{z}_1\left(t_0,t_1,t_2\right)+{\epsilon }^2{\dot{t}}_0^2{\displaystyle \frac{d^2}{d{t}_0^2}}{z}_2\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +{\epsilon }^2{\dot{t}}_1^2{\displaystyle \frac{d^2}{d{t}_1^2}}{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^2{\dot{t}}_2^2{\displaystyle \frac{d^2}{d{t}_2^2}}{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^2{\ddot{t}}_0{\displaystyle \frac{d}{d{t}_0}}{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^2{\ddot{t}}_1{\displaystyle \frac{d}{d{t}_1}}{z}_2\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +{\epsilon }^2{\ddot{t}}_2{\displaystyle \frac{d}{d{t}_2}}{z}_2\left(t_0,t_1,t_2\right)+2{\dot{t}}_0{\dot{t}}_1{\displaystyle \frac{d^2}{d{t}_{1}d{t}_0}}{z}_0\left(t_0,t_1,t_2\right)+2{\dot{t}}_0{\dot{t}}_2{\displaystyle \frac{d^2}{d{t}_{2}d{t}_0}}{z}_0\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +2{\dot{t}}_1{\dot{t}}_2{\displaystyle \frac{d^2}{d{t}_{2}d{t}_1}}{z}_0\left(t_0,t_1,t_2\right)+2\epsilon {\dot{t}}_0{\dot{t}}_1{\displaystyle \frac{d^2}{d{t}_{1}d{t}_0}}{z}_1\left(t_0,t_1,t_2\right)+2\epsilon {\dot{t}}_0{\dot{t}}_2{\displaystyle \frac{d^2}{d{t}_{2}d{t}_0}}{z}_1\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +2\epsilon {\dot{t}}_1{\dot{t}}_2{\displaystyle \frac{d^2}{d{t}_{2}d{t}_1}}{z}_1\left(t_0,t_1,t_2\right)+2{\epsilon }^2{\dot{t}}_0{\dot{t}}_1{\displaystyle \frac{d^2}{d{t}_{1}d{t}_0}}{z}_2\left(t_0,t_1,t_2\right)+2{\epsilon }^2{\dot{t}}_0{\dot{t}}_2{\displaystyle \frac{d^2}{d{t}_{2}d{t}_0}}{z}_2\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +2{\epsilon }^2{\dot{t}}_1{\dot{t}}_2{\displaystyle \frac{d^2}{d{t}_{2}d{t}_1}}{z}_2\left(t_0,t_1,t_2\right)\hfill \end{array}$$

(19)

where:

$\dot{t_0},\dot{t_1},\dot{t_2}$

- first time derivatives of time scales;

$\ddot{t_0},\ddot{t_1},\ddot{t_2}$

- second time derivatives of time scales.

The resulting approximated equation of motion is as follows:

$$\begin{array}{cc}\hfill \epsilon {\displaystyle \frac{d^2}{d{t}_0^2}}{z}_1\left(t_0,t_1,t_2\right)& +{\epsilon }^2{\displaystyle \frac{d^2}{d{t}_1^2}}{z}_0\left(t_0,t_1,t_2\right)+{\epsilon }^2{\displaystyle \frac{d^2}{d{t}_0^2}}{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^3{\displaystyle \frac{d^2}{d{t}_1^2}}{z}_1\left(t_0,t_1,t_2\right)+{\epsilon }^4{\displaystyle \frac{d^2}{d{t}_2^2}}{z}_0\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +{\epsilon }^4{\displaystyle \frac{d^2}{d{t}_1^2}}{z}_2\left(t_0,t_1,t_2\right)+{\epsilon }^5{\displaystyle \frac{d^2}{d{t}_2^2}}{z}_1\left(t_0,t_1,t_2\right)+{\epsilon }^6{\displaystyle \frac{d^2}{d{t}_2^2}}{z}_2\left(t_0,t_1,t_2\right)+2\epsilon {\displaystyle \frac{d^2}{d{t}_{1}d{t}_0}}{z}_0\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +2{\epsilon }^2{\displaystyle \frac{d^2}{d{t}_{2}d{t}_0}}{z}_0\left(t_0,t_1,t_2\right)+2{\epsilon }^2{\displaystyle \frac{d^2}{d{t}_{1}d{t}_0}}{z}_1\left(t_0,t_1,t_2\right)+2{\epsilon }^3{\displaystyle \frac{d^2}{d{t}_{2}d{t}_1}}{z}_0\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +2{\epsilon }^3{\displaystyle \frac{d^2}{d{t}_{2}d{t}_0}}{z}_1\left(t_0,t_1,t_2\right)+2{\epsilon }^3{\displaystyle \frac{d^2}{d{t}_{1}d{t}_0}}{z}_2\left(t_0,t_1,t_2\right)+2{\epsilon }^4{\displaystyle \frac{d^2}{d{t}_{2}d{t}_1}}{z}_1\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +2{\epsilon }^4{\displaystyle \frac{d^2}{d{t}_{2}d{t}_0}}{z}_2\left(t_0,t_1,t_2\right)+2{\epsilon }^5{\displaystyle \frac{d^2}{d{t}_{2}d{t}_1}}{z}_2\left(t_0,t_1,t_2\right)+256{\omega }^2{z}_0\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +256{\omega }^2\epsilon z_1\left(t_0,t_1,t_2\right)+256{\omega }^2{\epsilon }^2{z}_2\left(t_0,t_1,t_2\right)+256\delta {\omega }^2{u}_s\epsilon +256\delta {\omega }^2\epsilon z_0\left(t_0,t_1,t_2\right)\hfill \\ \hfill & +256\delta {\omega }^2{\epsilon }^2{z}_1\left(t_0,t_1,t_2\right)+256\delta {\omega }^2{\epsilon }^3{z}_2\left(t_0,t_1,t_2\right)+582.0{\omega }^2{u}_s\epsilon sin\left(\omega t_0\right)\hfill \\ \hfill & +582.0{\omega }^2\epsilon z_0\left(t_0,t_1,t_2\right)sin\left(\omega t_0\right)+582.0{\omega }^2{\epsilon }^2{z}_1\left(t_0,t_1,t_2\right)sin\left(\omega t_0\right)\hfill \\ \hfill & +582.0{\omega }^2{\epsilon }^3{z}_2\left(t_0,t_1,t_2\right)sin\left(\omega t_0\right)+{\displaystyle \frac{d^2}{d{t}_0^2}}{z}_0\left(t_0,t_1,t_2\right)=0\hfill \end{array}$$

(20)

The ordering and separation of equation (20) in terms of the power of a small parameter $\epsilon$ leads to obtaining a recursive sequence of linear equations of motion leading to the solution of the nonlinear equation. The zeroth-order approximate linear equation is given in (21):

$$\begin{array}{c}\hfill \begin{array}{c}\hfill 256{\omega }^2{z}_0\left(t_0,t_1,t_2\right)+{\displaystyle \frac{d^2}{d{t}_0^2}}{z}_0\left(t_0,t_1,t_2\right)=0\end{array}\end{array}$$

(21)

The next component of a recursive sequence of linear equations of motion leading to the solution of the considered nonlinear one is given in ():

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -{\displaystyle \frac{d^2}{d{t}_0^2}}{z}_1\left(t_0,t_1,t_2\right)-2{\displaystyle \frac{d^2}{d{t}_{1}d{t}_0}}{z}_0\left(t_0,t_1,t_2\right)-256{\omega }^2{z}_1\left(t_0,t_1,t_2\right)-256\delta {\omega }^2{u}_s-256\delta {\omega }^2{z}_0\left(t_0,t_1,t_2\right)-582.0{\omega }^2{u}_{s}sin\left(\omega t_0\right)-582.0{\omega }^2{z}_0\left(t_0,t_1,t_2\right)sin\left(\omega t_0\right)=0\end{array}\end{array}$$

(22)

The system was solved and an important issue was to remove the secular terms. When the secular terms were removed, the following solution was obtained in the form of a power series of a small parameter:

$$\begin{array}{cc}\hfill z=& 2.59u_s{\epsilon }^2+C_{1}sin\left(8\delta \omega t\epsilon \right)cos\left(16\omega t\right)+C_{1}sin\left(16\omega t\right)cos\left(8\delta \omega t\epsilon \right)+C_{2}cos\left(16\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -C_{2}sin\left(16\omega t\right)sin\left(8\delta \omega t\epsilon \right)-\delta u_s\epsilon +1.0{\delta }^2{u}_s{\epsilon }^2-2.28u_s\epsilon sin\left(\omega t\right)-2.64u_s{\epsilon }^{2}cos\left(2\omega t\right)\hfill \\ \hfill & +9.39C_1\epsilon sin\left(15\omega t\right)sin\left(8\delta \omega t\epsilon \right)+9.39C_2\epsilon sin\left(8\delta \omega t\epsilon \right)cos\left(15\omega t\right)\hfill \\ \hfill & +9.39C_2\epsilon sin\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)+37.7C_2{\epsilon }^{2}sin\left(18\omega t\right)sin\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & +8.82C_1\epsilon sin\left(17\omega t\right)sin\left(8\delta \omega t\epsilon \right)+8.82C_2\epsilon sin\left(8\delta \omega t\epsilon \right)cos\left(17\omega t\right)\hfill \\ \hfill & +8.82C_2\epsilon sin\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)+4.57\delta u_s{\epsilon }^{2}sin\left(\omega t\right)+45.5C_2{\epsilon }^{2}sin\left(14\omega t\right)sin\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -9.39C_1\epsilon cos\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)-37.7C_1{\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(18\omega t\right)\hfill \\ \hfill & -37.7C_1{\epsilon }^{2}sin\left(18\omega t\right)cos\left(8\delta \omega t\epsilon \right)-37.7C_2{\epsilon }^{2}cos\left(18\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -8.82C_1\epsilon cos\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)-45.5C_1{\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(14\omega t\right)\hfill \\ \hfill & -45.5C_1{\epsilon }^{2}sin\left(14\omega t\right)cos\left(8\delta \omega t\epsilon \right)-45.5C_2{\epsilon }^{2}cos\left(14\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & +4.28C_1\delta {\epsilon }^{2}cos\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)+4.84C_1\delta {\epsilon }^{2}cos\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -4.28C_1\delta {\epsilon }^{2}sin\left(17\omega t\right)sin\left(8\delta \omega t\epsilon \right)-4.28C_2\delta {\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(17\omega t\right)\hfill \\ \hfill & -4.28C_2\delta {\epsilon }^{2}sin\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)-4.84C_1\delta {\epsilon }^{2}sin\left(15\omega t\right)sin\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -4.84C_2\delta {\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(15\omega t\right)-4.84C_2\delta {\epsilon }^{2}sin\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \end{array}$$

(23)

where:

$C_1,C_2$

- integration constants.

It is necessary to determine the relative velocity from the first derivative for further calculations:

$$\begin{array}{cc}\hfill \dot{z}=& 16.0C_1\omega cos\left(16\omega t\right)cos\left(8\delta \omega t\epsilon \right)+5.27\omega u_s{\epsilon }^{2}sin\left(2\omega t\right)-16.0C_1\omega sin\left(16\omega t\right)sin\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -16.0C_2\omega sin\left(8\delta \omega t\epsilon \right)cos\left(16\omega t\right)-16.0C_2\omega sin\left(16\omega t\right)cos\left(8\delta \omega t\epsilon \right)-2.28\omega u_s\epsilon cos\left(\omega t\right)\hfill \\ \hfill & +141.0C_1\omega \epsilon sin\left(8\delta \omega t\epsilon \right)cos\left(15\omega t\right)+141.0C_1\omega \epsilon sin\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & +141.0C_2\omega \epsilon cos\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)+679.0C_1\omega {\epsilon }^{2}sin\left(18\omega t\right)sin\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & +679.0C_2\omega {\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(18\omega t\right)+679.0C_2\omega {\epsilon }^{2}sin\left(18\omega t\right)cos\left(8\delta \omega t\epsilon \right)+4.57\delta \omega u_s{\epsilon }^{2}cos\left(\omega t\right)\hfill \\ \hfill & +150.0C_1\omega \epsilon sin\left(8\delta \omega t\epsilon \right)cos\left(17\omega t\right)+150.0C_1\omega \epsilon sin\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & +150.0C_2\omega \epsilon cos\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)+637.0C_1\omega {\epsilon }^{2}sin\left(14\omega t\right)sin\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & +637.0C_2\omega {\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(14\omega t\right)+637.0C_2\omega {\epsilon }^{2}sin\left(14\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -141.0C_2\omega \epsilon sin\left(15\omega t\right)sin\left(8\delta \omega t\epsilon \right)-679.0C_1\omega {\epsilon }^{2}cos\left(18\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -150.0C_2\omega \epsilon sin\left(17\omega t\right)sin\left(8\delta \omega t\epsilon \right)-637.0C_1\omega {\epsilon }^{2}cos\left(14\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & +72.8C_2\delta \omega {\epsilon }^{2}sin\left(15\omega t\right)sin\left(8\delta \omega t\epsilon \right)+72.8C_2\delta \omega {\epsilon }^{2}sin\left(17\omega t\right)sin\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -72.8C_1\delta \omega {\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(15\omega t\right)-72.8C_1\delta \omega {\epsilon }^{2}sin\left(8\delta \omega t\epsilon \right)cos\left(17\omega t\right)\hfill \\ \hfill & -72.8C_1\delta \omega {\epsilon }^{2}sin\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)-72.8C_1\delta \omega {\epsilon }^{2}sin\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \\ \hfill & -72.8C_2\delta \omega {\epsilon }^{2}cos\left(15\omega t\right)cos\left(8\delta \omega t\epsilon \right)-72.8C_2\delta \omega {\epsilon }^{2}cos\left(17\omega t\right)cos\left(8\delta \omega t\epsilon \right)\hfill \end{array}$$

(24)

The MTSM made it possible to obtain approximate relationships describing the speed and relative position of the gear (measured along the line-of-action). The obtained dependencies have the character of a series of small parameter. This makes it possible to estimate the size of the error occurring in the approximate solution. In the case under consideration, it is at the level of ${\epsilon }^2$.

The obtained formulas allow analytical study of the dynamics of the considered system despite its approximate nature. In particular, to analyze the forces and slip rates occurring in the system.

4. Sensitivity Analysis

Several analyses were run to compare the signals of displacements z, velocities $\dot{z}$, and accelerations $\ddot{z}$. The objective of this comparison was to identify (demonstrate) which of the time-varying quantities better represent the dynamic behaviour of the system, i.e. vibration amplitudes or interim states (transients). It was assumed that the analyses would be performed for selected values of mesh stiffness, mesh damping and gear inertia. The corresponding results were compared in the time-domain (Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14) and frequency-domain (Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19). Reference values of the system’s parameters are presented in the Table 2.

The numerical investigation of these approximated analytical solutions began with a stiffness analysis. Figure 8 shows the time-varying displacement for eight mesh stiffness values.

Analysis of the vibration displacement of the gear system leads to the following conclusion. Observations of the system vibrations displacement based on the time-domain signals in Figure 8 and Figure 11) indicate the amplitude dependence on mesh stiffness. Lower stiffness results in greater vibration amplitude. The displacement amplitudes decrease consistently as gear mesh stiffness increases. A periodic pattern is evident in displacement. This has to do with the parametric nature of the gear system as evidenced by the meshing ratio, which indicates the time-varying number of teeth pairs in contact along the line-of-action. There is a clear periodic change in the displacement value when the second pair of teeth comes into contact, and this occurs when transitioning from single-tooth to double-tooth meshing. The time course of the displacement (tooth deflections) changes with a constant period, which correlates with the mesh stiffness. This pattern is due to the time when the second pair of teeth are in contact – specifically, from entering the line-of-action to exiting it. A lower displacement occurs when two pairs of teeth are in contact (within the line-of-action), resulting in higher mesh stiffness. Conversely, a higher displacement occurs when only one pair of teeth is in contact (along the line-of-action), leading to lower mesh stiffness. The occurrence of periodicity confirms the inherent time-varying nature of gear meshing and the associated occurrence of parametric vibration of the gear system

The mesh stiffness has a direct effect on the displacement of the gear system. Lower stiffness leads to higher amplitudes, indicative of weaker resistance to deformation. On the other hand, higher stiffness corresponds to more rapid stabilization and smaller oscillatory amplitudes. The trend suggests higher stiffness contributes to reduced system displacement vibrations. The signals show periodic vibrations with consistent periods across different stiffness values, indicating that gear system frequency remains largely unaffected by stiffness within the examined range.

Observations from the frequency-domain signals in Figure 15 allow to conclude that peaks in the displacement frequency spectrum appear at specific harmonic frequencies, corresponding to the gear mesh frequency ($f_z$ = 632.7 Hz) and its multiples. Gear systems concentrate vibratory energy at fundamental harmonics in the frequency domain. Peaks in the frequency spectrum diminish as mesh stiffness values increase. Notably, lower stiffness values show pronounced displacement energy concentrations, whereas higher stiffness disperses this energy. At higher values, the displacement energy reflects greater damping effects from increased stiffness. Increased mesh stiffness reduces displacement amplitude, suggesting a beneficial effect in terms of mechanical stability.

Observations of the system vibrations velocity based on the time-domain signals in Figure 9 and Figure 12 indicate the velocity signals are highly dynamic and exhibit a more nuanced representation of transient behavior than displacement. Velocity amplitudes increase in transients when there is a step-change in the mesh stiffness resulting from the time-varying number of teeth pairs along the line of action. However, velocity oscillations exhibit gradually reducing peak amplitudes as the mesh stiffness value increases. The frequency of these oscillations increases slightly for multiples of the mesh stiffness value. Higher values lead to quicker damping of velocity peaks in transients, reflecting enhanced energy dissipation dynamics. The zero-crossing frequency remains stable, consistent with the displacement analysis. It can be clearly seen the vibration velocity pattern corresponds to the time-varying number of teeth pairs in contact along the line-of-action. In steady states, the velocity oscillations amplitudes are constant, but slightly greater in the case of single-tooth meshing than in the case of double-tooth meshing.

The frequency-domain signals in Figure 16 and Figure 17 allow to conclude the velocity spectrum reveals sharp peaks at higher harmonic frequencies. These peaks become more distinct as stiffness increases. Similar to displacement, velocity spectrum peaks diminish with increasing stiffness. Increased gear mesh stiffness significantly dampens velocity oscillations and redistributes energy across a wider frequency range indicating a shift in system dynamic behavior towards higher-frequency components. Velocity analysis provides insights into the system’s dynamic transitions. Increased stiffness reduces peak velocity amplitudes, showcasing better control over vibrational dynamics. The dominance of higher harmonics in the frequency domain highlights the coupling effects of damping and inertia within the gear system. Velocity is a critical indicator of dynamic stability, showing smoother and milder transitions at higher stiffness. On the other hand, the Figure 16 shows that velocity oscillations are characterised by rapid (sharp) amplitude peaks in transients.

Observations of the system vibrations acceleration as the time-domain signals (Figure 10, Figure 13, Figure 14) show that the acceleration amplitudes exhibit a constant value regardless of the increase in mesh stiffness. The acceleration amplitude modulates, i.e. increases in transients when there is a step-change in the mesh stiffness resulting from a time-varying number of teeth pairs in contact along the line-of-action. The sharp increase in the amplitude of the mesh vibrations acceleration is particularly noticeable when there is a transition from single-tooth to double-tooth meshing (while the second pair of teeth additionally comes into contact along the line-of-action) – and somewhat less when there is a transition from double-tooth to single-tooth meshing (while the second pair of teeth leaves the line-of-action). The behavior remains consistent as can be seen in Figure 13 and Figure 14. To summarize, the periodic gear system excitation in transients is due to the peculiarities of parametric meshing during gear operation. Nevertheless, for different values of mesh stiffness, the vibration acceleration in the transient states is characterised by a constant peak amplitude, although the frequency of vibration increases.

Figure 18 and Figure 19 show the system’s vibration acceleration waveforms in the frequency domain. Peak acceleration for $k_g=3.0e+8$ $\mathrm{N}$/$\mathrm{m}$ was estimated for $664.741$ $\mathrm{m}$/${\mathrm{s}}^2$ at $5699.43$ $\mathrm{Hz}$. For analysis based on $k_g$ curve for value $4.0E+8$ $\mathrm{N}$/$\mathrm{m}$ it was observed, that acceleration was no higher than $615.442$ $\mathrm{m}$/${\mathrm{s}}^2$ and appeared at $5699.43$ $\mathrm{Hz}$. Simulations carried out for ${\mathrm{ks}}_g=5.0e+8$ $\mathrm{N}$/$\mathrm{m}$ showed, that the system responded in frequnecy domain with peak acceleration of $568.196$ $\mathrm{m}$/${\mathrm{s}}^2$ at $6999.3$ $\mathrm{Hz}$. Simulations for $k_g=6.0e+8$ $\mathrm{N}$/$\mathrm{m}$ resulted in $574.233$ $\mathrm{m}$/${\mathrm{s}}^2$ peak acceleration value appearing at $9499.05$ $\mathrm{Hz}$. Basing on numerical research carried out for $k_g=7.0e+8$ $\mathrm{N}$/$\mathrm{m}$ peak acceleration took value of $600.714$ $\mathrm{m}$/${\mathrm{s}}^2$ at $9499.05$ $\mathrm{Hz}$. Simulations carried out for $k_g=8.0e+8$ $\mathrm{N}$/$\mathrm{m}$ showed, that the system responded in frequnecy domain with peak acceleration of $598.76$ $\mathrm{m}$/${\mathrm{s}}^2$ at $9499.05$ $\mathrm{Hz}$. Analysis of an acceleration curve representing system’s response in frequency domain for $k_g$ value $9.0E+8$ $\mathrm{N}$/$\mathrm{m}$ resulted in peak of $570.935$ $\mathrm{m}$/${\mathrm{s}}^2$ at $9499.05$ $\mathrm{Hz}$. Analysis of the curve for $k_g$ value $1.0E+9$ $\mathrm{N}$/$\mathrm{m}$ showed that peak of acceleration was not higher than $531.315$ $\mathrm{m}$/${\mathrm{s}}^2$ and was registered at $9499.05$ $\mathrm{Hz}$.

Observations from the frequency-domain signals in Figure 18 and Figure 19 allow to conclude the acceleration waveforms demonstrate significant peaks during transitions between single-tooth and double-tooth meshing (i.e., single-pair and double-pair teeth meshing). However, the transient behavior (transitional condition) remains consistent as can be seen in Figure 13, and also in Figure 11 and Figure 12. The frequency associated with the gear meshing ($f_z$ = 632.7 Hz) has a dominant peak value, which is clearly visible at lower values of mesh stiffness. This peak frequency value remains visible, regardless of the stiffness, indicating a fixed excitation, although acceleration spectra show concentrated energy at higher frequencies as mesh stiffness increases. While peak magnitudes decrease as mesh stiffness increases, the high-frequency concentration is more pronounced at higher stiffness levels. This is because an increase in mesh stiffness causes an increase in transient vibration frequency, i.e., during transitions from two-pair teeth meshing to one-pair teeth meshing. When a single pair of teeth is in contact at the line of action, the lower mesh stiffness for single-tooth meshing enhances frequency effects, increasing the risk of fatigue and damage. In general, energy distribution expands across the spectrum as stiffness increases, which amplifies damping-like effects, as previously observed. This reflects the enhanced system strength under high stiffness conditions. It is also important to remember that acceleration is a vital metric for assessing material degradation and fatigue.

As can be seen in Figure 5, and also in Figure 2 and Figure 3.

The general observations are as follows. This analysis demonstrates the critical role of mesh stiffness of gear vibrational behavior in the drive system. Increased mesh stiffness improves overall vibrational stability by reducing displacement, velocity, and acceleration amplitudes., indicating better vibrational control. Fundamental frequencies remain stable across stiffness levels, while energy redistribution occurs towards higher frequencies as stiffness increases. Higher mesh stiffness values lead to broader energy dispersion across frequency spectra, suggesting enhanced mechanical stability. The energy redistribution across higher frequencies implies reduced wear and tear from parametric resonance. The reduction of resonance effects, especially low-frequency resonance risks, promotes a longer machinery lifespan.

The important thing is that transitions between single-tooth and double-tooth meshing produce noticeable peaks in all vibratory metrics. These transient effects are more pronounced in systems with lower stiffness, highlighting the need for robust design against such dynamic behavior. The sharpness of the transitions in the approximated waveforms might be influenced by the Gibbs phenomenon, which is an inherent characteristic of Fourier series approximations of periodic functions with low smoothness.

Based on the analyses, it can be concluded that the displacements in the time domain give more intuitive results than the other dynamic quantities, i.e. velocity and acceleration. On the other hand, the analysis in the frequency domain is clearer for $\dot{z}$ and $\ddot{z}$, since the mean value that dominates the spectrum of the displacement z (frequency equal to zero) is eliminated. In summary, it was decided that further simulation studies would be carried out as follows: with displacements being analysed in the time domain and accelerations in the frequency domain.

Subsequently, a study was carried out on the system’s dynamic behavior for different values of the reduced mass of the gears. Overall, Figure 20 and Figure 21 demonstrate that larger mass values result in more pronounced dynamic responses. For an assumed constant mesh stiffness, an increase in mass results in less favourable system dynamics, which need more time to stabilise (vibration damping). As can easily be seen, the greater the mass for an assumed gear stiffness, the greater the vibration amplitude. Observations from time-domain signals (Figure 21 and Figure 22) indicate the displacement amplitudes increase consistently as the equivalent mass ($m_{eq}$) increases from 0.2 to 3.0. Lower masses result in higher displacement vibrations. The amplitude increase with increasing mass suggests a direct correlation between mass and the system’s ability to suppress vibrational responses. It is also easy to see the values of equivalent mass affect the instability of the oscillation period (sweeping). The mass affects both the amplitude and frequency of the oscillation. For smaller masses displacement signals show concentrated energy at lower amplitudes, reflecting a system prone to larger oscillatory behavior.

In Figure 23 and Figure 24 the system’s acceleration response in the frequency domain is shown. Basing on numerical research carried out for $m_{eq}=0.2$ $\mathrm{k}$$\mathrm{g}$ peak acceleration took value of $575.504$ $\mathrm{m}$/${\mathrm{s}}^2$ at $9499.05$ $\mathrm{Hz}$. It was observed that for $m_{eq}$ value $0.4$ $\mathrm{k}$$\mathrm{g}$ peak appeared at $8199.18$ $\mathrm{Hz}$ with a value of $513.426$ $\mathrm{m}$/${\mathrm{s}}^2$. Analysis of an acceleration curve representing system’s response in frequency domain for $m_{eq}$ value $0.6$ $\mathrm{k}$$\mathrm{g}$ resulted in peak of $501.152$ $\mathrm{m}$/${\mathrm{s}}^2$ at $6999.3$ $\mathrm{Hz}$. Simulations carried out for $m_{eq}=0.8$ $\mathrm{k}$$\mathrm{g}$ showed, that the system responded in frequnecy domain with peak acceleration of $526.302$ $\mathrm{m}$/${\mathrm{s}}^2$ at $5699.43$ $\mathrm{Hz}$. Analysis of the curve for $m_{eq}$ value 1 $\mathrm{k}$$\mathrm{g}$ showed that peak of acceleration was not higher than $470.504$ $\mathrm{m}$/${\mathrm{s}}^2$ and was registered at $5699.43$ $\mathrm{Hz}$. For analysis based on $m_{eq}$ curve for value $1.5$ $\mathrm{k}$$\mathrm{g}$ it was observed, that acceleration was no higher than $438.296$ $\mathrm{m}$/${\mathrm{s}}^2$ and appeared at $4399.56$ $\mathrm{Hz}$. Analyzed $m_{eq}$ spectra for value 2 $\mathrm{k}$$\mathrm{g}$ had peak of $297.313$ $\mathrm{m}$/${\mathrm{s}}^2$ at $3199.68$ $\mathrm{Hz}$. Peak value for $m_{eq}=2$ $\mathrm{k}$$\mathrm{g}$ was estimated at $297.313$ $\mathrm{m}$/${\mathrm{s}}^2$ for $f=$$3199.68$ $\mathrm{Hz}$. Basing on numerical research carried out for $m_{eq}=3$ $\mathrm{k}$$\mathrm{g}$ peak acceleration took value of $455.45$ $\mathrm{m}$/${\mathrm{s}}^2$ at $2499.75$ $\mathrm{Hz}$.

Observations from frequency-domain signals in Figure 23 and Figure 24 indicate the lower masses ($m_{eq}$ < 0.6 kg) exhibit significantly higher acceleration peaks (greater susceptibility to vibration). Additionally for smaller masses acceleration spectra exhibit a distribution of vibrational energy over a broader range of frequencies with a preponderance of higher frequency peaks. The increase in mass narrows the range of vibrational energy distribution toward lower frequencies. Higher masses ($m_{eq}$ = 1.5 kg, and more) show a marked reduction in peak acceleration.

Howaever it should be noted the observed trends are specific to the analyzed gear system and operating conditions. The sensitivity analysis of the gear system to an increase in mass was carried out for specific parameters (Table 1) and reference values of simulation parameters (Table 2). Normally, it is the case that an increase in may is correlated with an increase in the geometric dimensions of the gearbox. This is accompanied by an increase in stiffness and damping. Thus, increased equivalent mass significantly enhances vibrational stability, reducing displacement and acceleration magnitudes under dynamic conditions act as natural dampers.

In these simulation studies, the increase in mass was implemented with an assumed value of stiffness and damping, which determines the dynamic behavior of the transmission system. Above a certain value of equivalent mass, the stability of the solution of the system mathematical model is lost. Figure 21 shows that in the case of a mass of $m_{eq}$ = 3 $\mathrm{k}$$\mathrm{g}$ the displacement periodically each obtain minimum negative values, what should be interpreted as a loss of contact between teeth (lack of meshing). Meanwhile, the model does not consider backlash. Therefore, the mass of 3kg is the limit, the model is not applicable above this value. Sensitivity analysis reveals the identification of the system’s parameters is essential, and the data cannot be selected completely arbitrarily. The limit values of the parameters determine the range of model applicability in the physical sense.

Nevertheless, the analysis showed the system’s sensitivity to a change in mass and its importance in managing vibration dynamics. This should be considered in the sustainable design of drive system gears.

The graphs in Figure 25, Figure 26, and Figure 27 show the system’s acceleration response over time for different damping coefficients.

Displacement amplitudes decrease as the mesh damping coefficient $c_g$ increases from 3000 to 90000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$, indicating logarithmic decrement increases, hence stabilization occurs faster and oscillations decrease. Lower damping values allow for larger displacement oscillations. As can be seen from Figure 27 as the mesh damping coefficient ($c_g$) increases the subsequent amplitudes are much smaller, indicating that the system is being damped more rapidly. For lower damping ($c_g$ = 3000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ and $c_g$ = 6000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$), displacement amplitudes are significantly higher and exhibit slower decay. For higher damping ($c_g$ = 60000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$, and $c_g$ = 90000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$), displacement signals rapidly approach minimal amplitudes, reflecting efficient energy dissipation. In particular, a zoomed-in view in Figure 26 shows that higher damping coefficients suppress peaks more effectively, confirming the role of cg in moderating displacement magnitudes. Increasing the mesh damping coefficient reduces displacement amplitudes and accelerates vibrational energy dissipation. Higher cg values stabilize the system more effectively, with a clear threshold effect evident at $c_g$ > 12000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$.

The effect of varying mesh damping is particularly evident during a two-pair mesh period when both pairs of teeth are along the line-of-action. In this case, displacement values are lower, the vibration period decreases, and the number of oscillations increases. In contrast, higher displacement values occur when only one pair of teeth is in contact (along the line-of-action). In this case, the mesh damping is less effective due to the periodic excitations caused by the second pair of teeth entering and exiting the line-of-action.

In turn, Figure 28 and Figure 29 show the acceleration of the system vibrations in the frequnecy domain for different mesh damping coefficients. Simulations for $c_g=3000$ $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ resulted in $1291.753$ $\mathrm{m}$/${\mathrm{s}}^2$ peak acceleration value appearing at $9499.05$ $\mathrm{Hz}$. Peak acceleration for $c_g=6000$ $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ was estimated for $797.233$ $\mathrm{m}$/${\mathrm{s}}^2$ at $9499.05$ $\mathrm{Hz}$. Simulations carried out for $c_g=9000$ $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ showed, that the system responded in frequnecy domain with peak acceleration of $598.76$ $\mathrm{m}$/${\mathrm{s}}^2$ at $9499.05$ $\mathrm{Hz}$. Analyzed $c_g$ spectra for value 12000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ had peak of $478.432$ $\mathrm{m}$/${\mathrm{s}}^2$ at $9499.05$ $\mathrm{Hz}$. Peak value for $c_g=12000$ $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ was estimated at $478.432$ $\mathrm{m}$/${\mathrm{s}}^2$ for $f=$$9499.05$ $\mathrm{Hz}$. Analysis of the curve for $c_g$ value 18000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ showed that peak of acceleration was not higher than $337.441$ $\mathrm{m}$/${\mathrm{s}}^2$ and was registered at $9499.05$ $\mathrm{Hz}$. For analysis based on $c_g$ curve for value 30000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ it was observed, that acceleration was no higher than $209.347$ $\mathrm{m}$/${\mathrm{s}}^2$ and appeared at $9499.05$ $\mathrm{Hz}$. For $c_g$ value 60000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$, peak acceleration of $106.256$ $\mathrm{m}$/${\mathrm{s}}^2$ was observed at $9499.05$ $\mathrm{Hz}$. Analyzed $c_g$ spectra for value 90000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ had peak of $71.217$ $\mathrm{m}$/${\mathrm{s}}^2$ at $5699.43$ $\mathrm{Hz}$. Peak value for $c_g=90000$ $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$ was estimated at $71.217$ $\mathrm{m}$/${\mathrm{s}}^2$ for $f=$$5699.43$ $\mathrm{Hz}$.

Similarily to the case of stiffness, the shape of the FFT plot remains almost unchaged. A predictible behavior of highest amplitude peaks is observed, where for increasing damping the accelertions decrease. The frequency at which highest peak appears remains constant except for the extremely high damping value $c_g>60000$ $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$.

Observations from frequency-domain signals in Figure 28 and Figure 29 show the acceleration peaks are highest for low damping coefficients ($c_g$ = 3000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$), indicating higher vibrational forces. As damping coefficient increases, peak acceleration values drop significantly, with the most substantial reductions observed for $c_g$ >= 18000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$. For lower damping, energy is concentrated at a fundamental frequency. Higher damping coefficients redistribute this energy across a broader range of frequencies, reducing resonance risks. At $c_g$ > 60000 $\mathrm{N}$$\mathrm{s}$/$\mathrm{m}$, peak acceleration values are minimal, with energy spread thinly across frequencies, confirming effective damping.

In general, Higher mesh damping coefficients enhance system stability by minimizing displacement and acceleration responses to dynamic loads. Acceleration amplitudes decrease sharply with higher damping coefficients. Fundamental frequencies remain stable across all damping coefficients, while higher $c_g$ redistributes vibrational energy across a broader spectrum. Energy redistribution across frequencies improves mechanical stability and mitigates resonance risks. This leads to increased resistance of the system to dynamic behavior.

Effective damping reduces vibrational energy more rapidly, promoting efficient system operation and reducing mechanical fatigue. Optimizing damping coefficients can mitigate vibrational forces, critical for improving the operational reliability and lifespan of machinery.

5. Application of Obtained Results for Energy Harvesting System

Piezoelectric energy harvesters are based on converting vibration displacements into electrical voltage. Therefore, they should be selected based on the analysis of the frequency of vibrations occurring in the system (related to vibration phenomena in the device).

The resonance frequency of the energy harvester system corresponds to the point at which the largest amount of electrical energy can be generated. Therefore, to achieve maximum efficiency of the energy harvester, a piezoelectric material with a resonance frequency close to the major vibration frequency in the system should be selected to ensure the resonance vibrations of the device, providing the largest amplitude at a given frequency.

To effectively select a piezoelectric material, it is necessary to thoroughly examine the vibration characteristics in the system. Both the highest vibration frequencies and locations should be identified. For this purpose, the following can be used: modal analysis, measurement methods, numerical simulations. In the case under consideration, a series of analyses were carried out, presented in the energy spectra Figure 30, Figure 31, Figure 32, Figure 33, Figure 34 and Figure 35 obtained with a formula:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill E_i=2{\pi }^2{A}_i^2{f}_i^2\end{array}\end{array}$$

(25)

where:

$E_i$

- energy of i-th frequecy stamp;

$A_i$

- amplitude of i-th frequency stamp;

$f_i$

- i-th frequency stamp.

Analysis of the plots (Figure 30, Figure 31, Figure 32, Figure 33, Figure 34 and Figure 35) allowed to draw conclusions about how the changes of values of subsequent parameters affect the energy generation. In the case of stiffness, the shape of the plots stays constant, with only difference in values of the energy - the higher the stiffness, the less energy is generated. Changes in the mass affect the energy waveforms both qualitatively and quantitatively. With the mass increasing, the overall values of energy also increase, but the span in which visible peaks appear shortens rapidly. When considering damping, with increase in it’s value, the generated energy decreases and bigger discrepancies in the plot shape are visible - the range of frequencies in which the highest peaks of energy appear gradually shortens and and shifts towards lower values.

In Figure 36, energy levels for reference data (Table 2) are shown. The highest energy peak appears at $5699.43$ $\mathrm{Hz}$, what constitutes the optimal working point of such pair of gears from the point of view of energy harvesting. Hence, as stated above the piezoelectric device should be chosen in such a manner, as to respond the best for said frequency value. If from any reason such application is not possible, further analysis of the plot shows at least five different points, at which the amount of generated energy would not be as high, but still reasonable in comparison with the highest peak. Application of piezoelectric device reposnding to frequencies higher than 1000 $\mathrm{Hz}$ would not provide satisfactory effect for such configuration of the system.

The harvester’s operating point is strongly related to the system parameters, so depending on their configuration, a different material and piezoelectric system solution should be selected. The analysis indicates that using analytical solutions allows the preliminary identification of the problem and estimation of the harvester’s performance. It should be emphasized that the presented analysis can be performed before conducting expensive experimental studies or at the design stage. That is why it is worth applying analytical solutions. Mathematical models allow for an initial estimation of the range of resonant frequencies and initial calculations of the profitability of using the energy harvester.

6. Conclusions

This study provides an analytical framework for understanding and optimizing gear mesh dynamics in machinery drive systems. The research bridges critical gaps in vibration analysis by employing the Multiple Scale Method (MSM) to model time-varying mesh stiffness and parametric vibrations. The main objective was to develop and validate analytical solutions and advance the understanding of gear’s dynamic models within drive trains. MSM is used to derive analytical solutions to determine the vibration frequency ranges and identify emitted components. The authors argue for the need to integrate such methods to optimize gearbox design, reduce vibrations, and improve operational reliability.

Using ISO 6336 standards and Fourier series, it formalizes stiffness variations during single-tooth and double-tooth meshing cycles. The model captures the periodic and parametric nature of gear vibrations, providing a solid mathematical foundation for analysis. Basd on this model the authors derived a series of approximations using MSM to solve nonlinear equations of motion for the gear system. These solutions reveal insights into displacement, velocity, and acceleration dynamics, emphasizing the role of damping and stiffness in modulating vibration amplitudes.

Based on the conducted simulation investigation of the gear system and comparison of the results for the analytical and numerical approaches, the following conclusions can be drawn:

1) The proposed model of varying stiffness based on the approximation of rectangular signal with Fourier series expansion allows for obtaining a more accurate representation of the meshing process.

2) The said approximated stiffness model stands as initially sufficient input for Hill’s equation, so the amplitudes and frequencies of the gear mechanism can be correctly identified.

3) The Hill’s equation approach itself enables vibration signals analysis of the working gear mechanism, despite the type of gear transmission system or gear meshing.

4) The analytical solution ensures the reliability and repeatability of results. A single run of the solving process allows results to be obtained for different parameters by substituting new values into the solution already given, which significantly reduces calculation time.

5) The proposed approach allows for quick and effective analysis of the vibracoustic properties of drive systems, including finding the operating point of machinery.

6) The possibility of energy-based determination of the machinery operating point is an important factor in the sustainable design of agricultural machinery. Using the analytical solution allows for a faster and more accurate selection of single gear (mass) and transmission (meshing stiffness, damping) physical parameters, so overall VA influence on operators and the environment can be diminished.

7) Energy analysis showed, that when considering different equivalent mass, mesh stiffness, and mesh damping values changes in these parameters influence the peak amplitude values and their position in the spectrum. Moreover, multiple operating points could be chosen, allowing for optimization regarding energy harvesting and other technical requirements for the system.

The analyses conducted showed that the operating point of the energy harvesting strongly depends on the system parameters. The presented results confirm the area application of the model for selecting energy harvesting parameters described in the introduction. Based on the simulations, it can be concluded that the model for selected gear parameters allows a preliminary determination of the amount of energy harvested and the operating point. Thus, it confirms the purpose of the article.

The general conclusion is that the Multiple Scale Method provides a robust analytical framework for predicting and managing parametric vibrations in gear systems. Increasing gear mesh stiffness effectively reduces vibrational amplitudes and shifts energy distribution to higher frequencies, enhancing stability. Optimised damping coefficients minimise transient oscillations, which contributes to improved machine durability and operator comfort.

The authors highlight the utility of their analytical solutions in enhancing gearbox configurations to minimize vibrations and improve durability and early detection of failures through vibration signal analysis. The application of MSM in gear mesh dynamics offers precision in predicting vibration behavior. The presented research results can be helpful in the sustainable design of gear meshing to optimise the drive systems dynamics, especially in agricultural machinery, where vibration control is important for the comfort and health of the operator. One can venture the far-reaching utilitarian conclusion that detailed mathematical modeling and a focus on computing performance can contribute to fostering safer and more efficient agricultural practices.