Effects of Uncertain Parameters on the Dynamic Response of a Wind Turbine Gearbox

1. Introduction

The wind turbine is a structure of variable sizes, depending on the desired energy production capacity. This system is intended to produce the electric current obtained in power station or for self-consumption. Currently, a wind turbine is composed of three main parts: a horizontal axis rotor, a nacelle and a tower. Vibrations in one of these parts are propagated to the other through relative movements between mechanical parts such as shafts, gears, pinions, etc [1]. In fact, gear transmission mechanisms are widely used for the connection and transfer of power in wind power generation. According to the heavy loads and poor operating conditions, the gear transmission mechanisms are likely subjected to a strangle vibration. Therefore, along with the advancement of industrial technologies. The vibrational characteristics of gear drive mechanisms are subject to stringent performance requirements. In-depth investigations into their dynamic responses offer valuable theoretical foundations for failure detection, performance optimization, and dependability assessment of gear-based power transmission systems [2,3]. Due to production and machining errors, material imperfections, degradation over time, environmental conditions, and several additional influences, certain mechanical properties such as weight, rigidity, and rotational resistance tend to exhibit variability. To accurately assess the system’s dynamic performance, it is crucial to analyze the motion behavior of gear mechanisms under uncertain conditions [4,5]. Various methods exist in the literature which MCS simulation is widely used technique in this context [6] It can provide the total probability density function of all the system variables, but it is generally considered to be costly because a large number of samples are required to reach an acceptable accuracy. Parallel simulation [7] and appropriate orthogonal decomposition [8] represent some solutions that have been considered to overcome the drawbacks associated with the MCS calculation. These techniques consist of solving the deterministic problem repeatedly for random parameters’ values and extracting the closest solution to this problem. The MCS technique is widely used in this field, but remains generally excessively expensive, as a significant number of samples is necessary to obtain satisfactory results.

To achieve information about the uncertainty in terms of first and second statistical moments quickly enough, the perturbation method is effective [9]. In this case, the statistical moments are determined by expanding the stochastic quantities around their mean via a Taylor series expansion. However, its use is limited to small perturbations. Nevertheless, the PC method has already been shown to deal effectively for conflicts associated with uncertainty. It is a development in orthogonal polynomials based on random variables, aimed at approximating the overall distribution of output uncertainty. The main homogeneous polynomial chaos method is founded on Wiener’s homogeneous chaos theory [10,11]. It utilizes Hermite polynomials in conjunction with Gaussian random variables. The effectiveness of this approach has been demonstrated in various application fields, particularly in managing uncertainties in environmental and biological phenomena [12,13] and multi-body dynamical systems [14,15,16,17], the solution of ordinary and partial differential equations and parameters estimation.

Many works have been reached on the reliability of vibrating structures considering the uncertainties.Guerine and El Hami [1], examined the nonlinear dynamic behavior of a mechanical torque limiter designed to safeguard drivetrain components such as helical transmission gears—against overload. In a previous study, Guerine et al. (2017) [18] analyzed the dynamic response of wind turbine gear systems with uncertain yet bounded parameters, employing the interval analysis method. They also studied the stability of an eight-degree of freedom spur gear according to the coefficient of friction value. Wang et al. (2021) [19] studied an uncertain transient response of the dual disk rotor system using a non-probabilistic Chebyshev convex method (CCM). Zhu et al. (2018) [20] have studied a nonlinear coupled lateral-torsional dynamic model with 16 degrees of freedom for the gear-rotor-bearing system. They have considered nonlinear characteristics such as variable gear stiffness, backlash, transmission error, and friction force. Their research focused on the effects of the friction coefficient and mean load, revealing the complex influence of these factors on the behavior of the gear system. Giannella (2021) [21] presented a stochastic model, which considers various sources of uncertainty, to define the probability distribution of the residual fatigue life of a defective railway axle. Garoli and Castro (2019) [22] considered the consequences of temperature uncertainty on the dynamic coefficient of the bearing using the PC to solve the problem of the vortex instability of the oil film of the rotor-bearing system. To address the issue of reliability, specifically the estimation of very low exceedance probabilities, Wang and Zhang (2021) [23] analyzed the effect of stochastic interval variables such as meshing rigidity, driving moment, kinematic inaccuracy, and backlash on the oscillatory velocity spectrum and operational dependability of gear train mechanisms. In their research, the dynamic system incorporating interval uncertainties was resolved through the integration of the Chebyshev enclosure function technique and the Runge-Kutta numerical integration algorithm.

In addition, it should be noted that the literature has not thoroughly investigated the dynamic behavior of a two-stage gearbox used in wind turbines, particularly with the presence of uncertain parameters. Furthermore, most studies have only focused on analyzing a few uncertain parameters, such as the input torque or the mesh stiffness. However, to the best of our knowledge, there is a lack of research that considers the simultaneous effects of multiple uncertain parameters. In this study, we have considered seven uncertain parameters that have a significant impact on the dynamic response of the gearbox. Therefore, this research contributes to addressing this gap in the literature by providing a comprehensive analysis of the dynamic behavior of a two-stage gearbox under uncertainty, considering multiple uncertain parameters. In addition, the literature offers several approaches to addressing uncertainty, typically classified into four main categories. The first category is the deterministic MCS approach, which is commonly used as a reference due to its ability to produce reliable results that are close to reality. However, the MCS can be computationally intensive, limiting its applicability. To minimize the computational cost of simulations with uncertainties, other probabilistic approaches have been developed, such as the perturbation method, PC method, and the meta-model method. These methods provide efficient ways to generate accurate results. Finally, the literature also includes a possibilistic approach, such as the interval computing method. These methods provide an alternative framework for addressing uncertainty and offer valuable insight into various fields of study.

This study presents a proposed approach to model the dynamic behavior of a two-stage gearbox used in a horizontal wind turbine while accounting for uncertainties related to the mass, the moment of inertia, the bearing’s bending stiffness, the traction-compression stiffness, the damping coefficient and the wind speeds. The main innovation of this study lies in its ambitious approach to include as many parameters as possible that could potentially contribute to the random vibrations of the system. Setting it apart from the existing literature. By comprehensively analyzing this large number of variables and considering three differents standard deviations, our aim is to obtain a more realistic understanding of the behavior and operation of the studied system in the face of varying uncertainty conditions. Furthermore, by studying the effects of multiple sources of uncertainty on the system’s random vibrations with different standard deviations, we can define their influence on the dynamic response and identify the critical thresholds where the system’s behavior could be most sensitive to uncertainty variations. In this way, this approach provides a better understanding of the standard deviation thresholds that affect system robustness and stability under real operating conditions and directs our efforts towards improving the design to minimize the negative effects of undesirable standard deviations on system performance.

2. Overview of Mechanical Architecture

A wind energy conversion system is a rotary apparatus that extracts kinetic energy from a fluid medium (i.e., atmospheric airflow) via aerodynamic airfoils and transduces it into utilizable mechanical power [24]. Figure 1 illustrates the nonlinear dynamic representation of the analyzed horizontal-axis wind turbine (HAWT). The rotor assembly typically consists of a tri-blade configuration affixed to a central hub, which is mechanically coupled to a speed-increasing gearbox through a transmission shaft supported by an initial compliant bearing element.

The comprehensive dynamic formulation is composed of three subsystems $\left(\mathrm{j}=1,...,3\right)$, each supported by a compliant bearing assembly. For the first subsystem, the lateral (bending) stiffness of the bearing is denoted by${\mathrm{k}}_1^{\mathrm{x}}$, and the axial (traction-compression) stiffness is designated as ${\mathrm{k}}_1^{\mathrm{y}}$. Analogously, for the second subsystem, the stiffness parameters are${\mathrm{k}}_2^{\mathrm{x}}$and ${\mathrm{k}}_2^{\mathrm{y}}$.While for the third subsystem they are represented as ${\mathrm{k}}_3^{\mathrm{x}}$and ${\mathrm{k}}_3^{\mathrm{y}}$. The interconnecting shafts among these blocks exhibit torsional rigidity, quantified by ${\mathrm{k}}_1^{\mathsf{\theta }}$,${\mathrm{k}}_2^{\mathsf{\theta }}$ and ${\mathrm{k}}_3^{\mathsf{\theta }}$.It is assumed that the intermediate flexible couplings possess negligible inertial mass relative to that of the gearbox assembly. The angular displacements associated with each rotational component are denoted respectively as${\mathsf{\theta }}_{\left(11\right)}$,${\mathsf{\theta }}_{\left(1,2\right)}$,${\mathsf{\theta }}_{\left(2,1\right)}{\mathsf{\theta }}_{\left(2,2\right)}$,${\mathsf{\theta }}_{\left(3,1\right)}$and${\mathsf{\theta }}_{\left(3,2\right)}$. Moreover, the translational displacements of the bearing supports, in the plane orthogonal to the axis of rotation, are expressed as ${\mathrm{x}}_1$ and ${\mathrm{y}}_1$ for the first subsystem, ${\mathrm{x}}_2$ and ${\mathrm{y}}_2$ for the second, ${\mathrm{x}}_3$ and ${\mathrm{y}}_3$ for the third.

3. Equations of Motion

The dynamic behavior of the investigated system (Figure 1) is governed by the equations of motion formulated via the Lagrangian approach. These equations are represented in the canonical matrix format as follows:

$$\left[\mathrm{M}\right]\left\{\ddot{\mathrm{q}}\right\}+\left[\mathrm{C}\right]\left\{\dot{\mathrm{q}}\right\}+\left(\left[{\mathrm{K}}_{\mathrm{S}}\right]+\left[\mathrm{K}\left(\mathrm{t}\right)\right]\right)\left\{\mathrm{q}\right\}=\left\{{\mathrm{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\right\}$$

(1)

Where:

$$\left\{\mathrm{q}\left(\mathrm{t}\right)\right\}={\left\{{\mathrm{x}}_1,{\mathrm{y}}_1,{\mathrm{x}}_2,{\mathrm{y}}_2,{\mathrm{x}}_3,{\mathrm{y}}_3,{\mathsf{\theta }}_{11},{\mathsf{\theta }}_{12},{\mathsf{\theta }}_{21},{\mathsf{\theta }}_{22},{\mathsf{\theta }}_{31},{\mathsf{\theta }}_{32}\right\}}^{\mathrm{T}}$$

(2)

$⟨{\mathrm{L}}_{\mathsf{\delta }\mathrm{n}}⟩$ is a function of the n-level.

$$⟨{\mathrm{L}}_{\mathsf{\delta }1}⟩=\left[-\sin \left({\mathsf{\alpha }}_1\right),\cos \left({\mathsf{\alpha }}_1\right),\sin \left({\mathsf{\alpha }}_1\right),-\cos \left({\mathsf{\alpha }}_1\right),\mathrm{0,0},0,{\mathrm{r}}_{12}^{\mathrm{b}},{\mathrm{r}}_{21}^{\mathrm{b}},\mathrm{0,0},0\right]$$

(3)

$$⟨{\mathrm{L}}_{\mathsf{\delta }2}⟩=\left[-\sin \left({\mathsf{\alpha }}_2\right),\cos \left({\mathsf{\alpha }}_2\right),\sin \left({\mathsf{\alpha }}_2\right),-\cos \left({\mathsf{\alpha }}_2\right),\mathrm{0,0},0,{\mathrm{r}}_{22}^{\mathrm{b}},{\mathrm{r}}_{31}^{\mathrm{b}},\mathrm{0,0},0\right]$$

(4)

The parameters ${\mathrm{r}}_{12}^{\mathrm{b}}$ and ${\mathrm{r}}_{21}^{\mathrm{b}}$ denote the base radii of gears (12) and (21), respectively. Similarly, ${\mathrm{r}}_{22}^{\mathrm{b}}$and ${\mathrm{r}}_{31}^{\mathrm{b}}$ correspond to the base radii of gears (22) and (31), respectively. The pressure angles associated with the first and second gear stages are represented by${\mathsf{\alpha }}_1$ and ${\mathsf{\alpha }}_2$, respectively.

The global mass matrix of the mechanical system is expressed in the following block matrix form:

$$\left[{\mathrm{M}}_{\mathrm{G}}\right]=\left[\begin{array}{cc}{\mathrm{M}}_{\mathrm{L}}& 0\\ 0& {\mathrm{M}}_{\mathrm{A}}\end{array}\right]$$

(5)

Here, $\left[{\mathrm{M}}_{\mathrm{L}}\right]$denotes the mass matrix associated with the rotating components (wheels), and is defined as:

$$\left[{\mathrm{M}}_{\mathrm{L}}\right]=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathrm{m}}_1,{\mathrm{m}}_1,{\mathrm{m}}_{2,}{\mathrm{m}}_2,{\mathrm{m}}_3,{\mathrm{m}}_3\right)$$

(6)

${\mathrm{m}}_{\mathrm{j}}$is the mass of the bloc j.

$\left[{\mathrm{M}}_{\mathrm{A}}\right]$is the matrix of inertia and is written as follows:

$$\left[{\mathrm{M}}_{\mathrm{A}}\right]=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathrm{I}}_{11},{\mathrm{I}}_{12},{\mathrm{I}}_{21},{\mathrm{I}}_{22},{\mathrm{I}}_{31},{\mathrm{I}}_{32}\right)$$

(7)

${\mathrm{I}}_{\mathrm{ji}}$ are the polar inertia of the Wheel (ji) relative to the axis of rotation.

To attenuate the vibratory response of the system, a proportional damping matrix $\left[\mathrm{C}\right]$ was implemented, following the methodology outlined in [24,25].

$$\left[\mathrm{C}\right]=\mathsf{\eta }\left[\mathrm{M}\right]+\mathsf{\mu }\left(\left[{\mathrm{K}}_{\mathrm{S}}+\left[{\mathrm{K}}_{\mathrm{c}}\right]\right]\right)$$

(8)

Where $\left[{\mathrm{K}}_{\mathrm{c}}\right]$ is the mean component of mesh stiffness

With: $\mathsf{\eta }$ and $\mathsf{\mu }$ are the damping constants defined by: $\mathsf{\eta }=10^{-2};\mathsf{\mu }=10^{-5}$ [25].

The global stiffness matrix $\left[{\mathrm{K}}_{\mathrm{s}}\right]$is constructed from the superposition of the mean structural stiffness matrix, the torsional rigidity matrix $\left[{\mathrm{K}}_{\mathsf{\theta }}\right]$of the interconnecting shafts, and the elastic stiffness matrix of the bearing$\left[{\mathrm{K}}_{\mathrm{p}}\right]$ assemblies [25].

$$\left[{\mathrm{K}}_{\mathrm{s}}\right]=\left[\begin{array}{cc}{\mathrm{K}}_{\mathrm{p}}& 0\\ 0& {\mathrm{K}}_{\mathsf{\theta }}\end{array}\right]$$

(9)

The gear-mesh stiffness matrix $\left[\mathrm{K}\left(\mathrm{t}\right)\right]$ is expressed by:

$$\left[\mathrm{K}\left(\mathrm{t}\right)\right]=\left[\begin{array}{cc}\left[{\mathrm{K}}_1\left(\mathrm{t}\right)\right]& \left[{\mathrm{K}}_{12}\left(\mathrm{t}\right)\right]\\ \left[{\mathrm{K}}_{12}^{\mathrm{T}}\left(\mathrm{t}\right)\right]& \left[{\mathrm{K}}_2\left(\mathrm{t}\right)\right]\end{array}\right]$$

(10)

With:

$$\left[{\mathrm{K}}_1\left(\mathrm{t}\right)\right]=\left[\begin{array}{cccccc}{\mathrm{k}}_1{\mathrm{x}}_1^2& -{\mathrm{k}}_1\dot{\mathrm{x}}& -{\mathrm{k}}_1{\mathrm{x}}_1^2& {\mathrm{k}}_1\dot{\mathrm{x}}& 0& 0\\ {\mathrm{k}}_1\dot{\mathrm{x}}& {\mathrm{k}}_1{\dot{\mathrm{y}}}^2& {\mathrm{k}}_1\dot{\mathrm{x}}& -{\mathrm{k}}_1{\dot{\mathrm{x}}}^2& 0& 0\\ -{\mathrm{k}}_1{\mathrm{x}}_1^2& {\mathrm{k}}_1\dot{\mathrm{x}}& {\mathrm{k}}_1{\mathrm{x}}_1^2+{\mathrm{k}}_2{\mathrm{x}}_2^2& -{\mathrm{k}}_1\dot{\mathrm{x}}-{\mathrm{k}}_2\ddot{\mathrm{x}}& -{\mathrm{k}}_2{\mathrm{x}}_2^2& {\mathrm{k}}_2\ddot{\mathrm{x}}\\ {\mathrm{k}}_1\dot{\mathrm{x}}& {\mathrm{k}}_1{\dot{\mathrm{y}}}^2& -{\mathrm{k}}_1\dot{\mathrm{x}}-{\mathrm{k}}_2\ddot{\mathrm{x}}& {\mathrm{k}}_1{\dot{\mathrm{y}}}^2+{\mathrm{k}}_2{\dot{\mathrm{y}}}^2& {\mathrm{k}}_2\ddot{\mathrm{x}}& {\mathrm{k}}_2{\ddot{\mathrm{y}}}^2\\ 0& 0& -{\mathrm{k}}_2{\mathrm{x}}_2^2& {\mathrm{k}}_2\ddot{\mathrm{x}}& {\mathrm{k}}_2{\mathrm{x}}_2^2& -{\mathrm{k}}_2\ddot{\mathrm{x}}\\ 0& 0& {\mathrm{k}}_2\ddot{\mathrm{x}}& {\mathrm{k}}_2{\ddot{\mathrm{y}}}^2& {\mathrm{k}}_2\ddot{\mathrm{x}}& {\mathrm{k}}_2{\ddot{\mathrm{y}}}^2\end{array}\right]$$

(11)

$$\left[{\mathrm{K}}_{12}\left(\mathrm{t}\right)\right]=\left[\begin{array}{cccccc}0& {\mathrm{k}}_1{\mathrm{r}}_{12}^{\mathrm{b}}{\mathrm{x}}_1& {\mathrm{k}}_1{\mathrm{r}}_{21}^{\mathrm{b}}{\mathrm{x}}_1& 0& 0& 0\\ 0& {\mathrm{k}}_1{\mathrm{r}}_{12}^{\mathrm{b}}\dot{\mathrm{y}}& {\mathrm{k}}_1{\mathrm{r}}_{21}^{\mathrm{b}}\dot{\mathrm{y}}& 0& 0& 0\\ 0& {\mathrm{k}}_1{\mathrm{r}}_{21}^{\mathrm{b}}{\mathrm{x}}_1& {\mathrm{k}}_1{\mathrm{r}}_{12}^{\mathrm{b}}{\mathrm{x}}_1& {\mathrm{k}}_2{\mathrm{r}}_{22}^{\mathrm{b}}{\mathrm{x}}_2& {\mathrm{k}}_2{\mathrm{r}}_{31}^{\mathrm{b}}{\mathrm{x}}_2& 0\\ 0& -{\mathrm{k}}_1{\mathrm{r}}_{12}^{\mathrm{b}}\dot{\mathrm{y}}& -{\mathrm{k}}_1{\mathrm{r}}_{21}^{\mathrm{b}}\dot{\mathrm{y}}& -{\mathrm{k}}_2{\mathrm{r}}_{22}^{\mathrm{b}}\ddot{\mathrm{y}}& -{\mathrm{k}}_2{\mathrm{r}}_{31}^{\mathrm{b}}\ddot{\mathrm{y}}& 0\\ 0& 0& 0& -{\mathrm{k}}_2{\mathrm{r}}_{22}^{\mathrm{b}}{\mathrm{x}}_2& -{\mathrm{k}}_2{\mathrm{r}}_{31}^{\mathrm{b}}{\mathrm{x}}_2& 0\\ 0& 0& 0& {\mathrm{k}}_2{\mathrm{r}}_{22}^{\mathrm{b}}\ddot{\mathrm{y}}& {\mathrm{k}}_2{\mathrm{r}}_{31}^{\mathrm{b}}\ddot{\mathrm{y}}& 0\end{array}\right]$$

(12)

$$\left[{\mathrm{K}}_2\left(\mathrm{t}\right)\right]=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& {\mathrm{k}}_1{\mathrm{r}}_{12}^{\mathrm{b}}& {\mathrm{k}}_1{\mathrm{r}}_{12}^{\mathrm{b}}{\mathrm{r}}_{21}^{\mathrm{b}}& 0& 0& 0\\ 0& {\mathrm{k}}_1{\mathrm{r}}_{12}^{\mathrm{b}}{\mathrm{r}}_{21}^{\mathrm{b}}& {\mathrm{k}}_1{\mathrm{r}}_{21}^{\mathrm{b}}& 0& 0& 0\\ 0& 0& 0& {\mathrm{k}}_2{\mathrm{r}}_{22}^{\mathrm{b}}& {\mathrm{k}}_2{\mathrm{r}}_{22}^{\mathrm{b}}{\mathrm{r}}_{31}^{\mathrm{b}}& 0\\ 0& 0& 0& {\mathrm{k}}_2{\mathrm{r}}_{22}^{\mathrm{b}}{\mathrm{r}}_{31}^{\mathrm{b}}& {\mathrm{k}}_2{\mathrm{r}}_{31}^{\mathrm{b}}& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$$

(13)

Where:

$\left[{\mathrm{K}}_1\left(\mathrm{t}\right)\right]$ is the first-time varying stiffness.

$\left[{\mathrm{K}}_2\left(\mathrm{t}\right)\right]$is the second-time varying stiffness.

${\mathrm{x}}_1=\sin \left({\mathsf{\alpha }}_1\right),{\mathrm{x}}_2=\sin \left({\mathsf{\alpha }}_2\right),\dot{\mathrm{x}}=\left(\sin \left({\mathsf{\alpha }}_1\right)\cos \left({\mathsf{\alpha }}_1\right)\right),\ddot{\mathrm{x}}=\left(\sin \left({\mathsf{\alpha }}_2\right)\cos \left({\mathsf{\alpha }}_2\right)\right),\dot{\mathrm{y}}=\cos \left({\mathsf{\alpha }}_1\right),\ddot{\mathrm{y}}=\cos \left({\mathsf{\alpha }}_2\right)$The external vector force is defined by:

$$\left\{{\mathrm{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\right\}={\left\{\mathrm{0,0},\mathrm{0,0},\mathrm{0,0},{\mathrm{C}}_{\mathrm{a}\mathrm{e}\mathrm{r}}\left(\mathrm{t}\right),\mathrm{0,0},\mathrm{0,0},{\mathrm{C}}_{\mathrm{g}\mathrm{e}\mathrm{n}}\left(\mathrm{t}\right)\right\}}^{\mathrm{T}}$$

(14)

Where${\mathrm{C}}_{\mathrm{a}\mathrm{e}\mathrm{r}}\left(\mathrm{t}\right)$is the aerodynamic torque and${\mathrm{C}}_{\mathrm{g}\mathrm{e}\mathrm{n}}\left(\mathrm{t}\right)$is the generator torque.

Under ideal operational conditions, the system is subjected to both external and internal excitation sources. The primary contributors to these dynamic excitations include the stochastic variability of the wind energy input and the time-dependent fluctuation of the gear meshing stiffness [25]. In the current investigation, only aerodynamic excitation due to wind variability is taken into account [26]. Therefore, the wind turbine speed ${\mathsf{\Omega }}_{\mathrm{w}\left(\mathrm{t}\right)}$is considered as a fixed value for an average level of wind magnitude.

The wind velocity is expressed as follow [25].

$${\mathrm{V}}_{\mathrm{w}\mathrm{i}\mathrm{n}}\left(\mathrm{t}\right)={\mathrm{V}}_{\mathrm{m}\mathrm{o}\mathrm{y}}\left[\left.\right)\right]$$

(15)

Where $\left({\mathsf{\omega }}_1=\mathsf{\pi }\right)$ and $\left({\mathsf{\omega }}_2=2\mathsf{\pi }\right)$ represent the angular frequencies (rad/s).

Under these operating conditions, the instantaneous aerodynamic torque is defined as follows:

$${\mathrm{C}}_{\mathrm{a}\mathrm{e}\mathrm{r}}\left(\mathrm{t}\right)={\displaystyle \frac{1}{2{\mathsf{\Omega }}_{\mathrm{w}\left(\mathrm{t}\right)}}}\mathsf{\pi }{\mathsf{\rho }}_{\mathrm{a}\mathrm{i}\mathrm{r}}{\mathrm{R}}^2{\mathrm{V}}_{\mathrm{w}\mathrm{i}\mathrm{n}\left(\mathrm{t}\right)}^3{\mathrm{C}}_{\mathrm{p}}$$

(16)

Where:

${\mathrm{C}}_{\mathrm{p}}$,$R$and${\mathsf{\rho }}_{\mathrm{air}}$represent respectively the power coefficient, the rotor radius and the air density.

The generator torque ${\mathrm{C}}_{\mathrm{g}\mathrm{e}\mathrm{n}}\left(\mathrm{t}\right)$ is dependent on the aerodynamic torque ${\mathrm{C}}_{\mathrm{a}\mathrm{e}\mathrm{r}}\left(\mathrm{t}\right)$ and the gearbox ratio. It is defined by:

$${\mathrm{C}}_{\mathrm{g}\mathrm{e}\mathrm{n}}\left(\mathrm{t}\right)={\displaystyle \frac{-{\mathrm{C}}_{\mathrm{a}\mathrm{e}\mathrm{r}}\left(\mathrm{t}\right)}{{\mathsf{\eta }}_{\mathrm{g}\mathrm{e}\mathrm{a}\mathrm{r}}}}$$

(17)

where ${\mathsf{\eta }}_{\mathrm{g}\mathrm{e}\mathrm{a}\mathrm{r}}$ represent the gear ratio and is calculated as follows:

$${\mathsf{\eta }}_{\mathrm{g}\mathrm{e}\mathrm{a}\mathrm{r}}={\displaystyle \frac{{\mathrm{Z}}_{12}{\mathrm{Z}}_{21}}{{\mathrm{Z}}_{22}{\mathrm{Z}}_{31}}}$$

(18)

With:

${\mathrm{Z}}_{12}$,${\mathrm{Z}}_{21}$,${\mathrm{Z}}_{22}$and ${\mathrm{Z}}_{31}$are the teeth number of n gear stage.

The motion of the wind turbine blocks, which involve both linear and torsional displacements, can be mathematically described by a system of differential equations, as presented below.

$$|\begin{array}{l}m{\ddot{\mathrm{x}}}_1-{\mathrm{c}}^{\mathrm{x}}{\dot{\mathrm{x}}}_1+\sin \left({\mathsf{\alpha }}_1\right){\mathrm{c}}_1{\mathsf{\delta }}_1\left(\mathrm{t}\right)+{\mathrm{k}}^{\mathrm{x}}{\mathrm{x}}_1-{\mathrm{k}}_1{\mathsf{\delta }}_1\left(\mathrm{t}\right)\sin \left({\mathsf{\alpha }}_1\right)=0\\ m{\ddot{\mathrm{y}}}_1-{\mathrm{c}}^{\mathrm{y}}{\dot{\mathrm{y}}}_1+\cos ({\mathsf{\alpha }}_1\left){\mathrm{c}}_1\right(t\left){\mathsf{\delta }}_1\right(t)+{\mathrm{k}}^{\mathrm{y}}{\mathrm{y}}_1+{\mathrm{k}}_1(t\left){\mathsf{\delta }}_1\right(t\left)\cos ({\mathsf{\alpha }}_1\right)=0\\ m{\ddot{\mathrm{x}}}_2+{\mathrm{c}}^{\mathrm{x}}{\dot{\mathrm{x}}}_2-\sin ({\mathsf{\alpha }}_1\left){\mathrm{c}}_1\right(t\left){\mathsf{\delta }}_1\right(t)+\sin ({\mathsf{\alpha }}_2){\mathrm{c}}_2\left(t\right){\mathsf{\delta }}_2\left(t\right)+{\mathrm{k}}^{\mathrm{x}}{\mathrm{x}}_2-{\mathrm{k}}_1\left(t\right)\sin ({\mathsf{\alpha }}_1)+{\mathrm{k}}_2(t\left)\sin ({\mathsf{\alpha }}_2\right){\mathsf{\delta }}_2\left(t\right)=0\\ m{\ddot{\mathrm{y}}}_2+{\mathrm{c}}^{\mathrm{y}}{\dot{\mathrm{y}}}_2-{\mathrm{k}}_1\left(t\right){\mathsf{\delta }}_1\left(t\right)\cos ({\mathsf{\alpha }}_1)+\cos ({\mathsf{\alpha }}_2){\mathrm{c}}_2\left(t\right){\mathsf{\delta }}_2\left(t\right)+{\mathrm{k}}^{\mathrm{y}}{\mathrm{y}}_2-{\mathrm{k}}_1\left(t\right)\cos ({\mathsf{\alpha }}_1\left){\mathsf{\delta }}_1\right(t)+{\mathrm{k}}_2(t\left){\mathsf{\delta }}_2\right(t\left)\cos ({\mathsf{\alpha }}_2\right)=0\\ m{\ddot{\mathrm{x}}}_3-{\mathrm{c}}^{\mathrm{x}}{\dot{\mathrm{x}}}_3-\sin ({\mathsf{\alpha }}_2\left){\mathrm{c}}_2\right(t\left){\mathsf{\delta }}_2\right(t)+{\mathrm{k}}^{\mathrm{x}}{\mathrm{x}}_3-{\mathrm{k}}_2(t\left){\mathsf{\delta }}_2\sin ({\mathsf{\alpha }}_2\right)=0\\ m{\ddot{\mathrm{y}}}_3-{\mathrm{c}}^{\mathrm{y}}{\dot{\mathrm{y}}}_3-\cos ({\mathsf{\alpha }}_2\left){\mathrm{c}}_2\right(t\left){\mathsf{\delta }}_2\right(t)+{\mathrm{k}}^{\mathrm{y}}{\mathrm{y}}_3-{\mathrm{k}}_2(t\left){\mathsf{\delta }}_2\cos ({\mathsf{\alpha }}_2\right)=0\end{array}$$

(19)

$$|\begin{array}{l}{\mathrm{I}}_{\mathrm{b}}{\ddot{\mathsf{\theta }}}_{\left(\mathrm{1,1}\right)}+{\mathrm{k}}^{\mathsf{\theta }}({\mathsf{\theta }}_{\left(\mathrm{1,1}\right)}-{\mathsf{\theta }}_{\left(\mathrm{1,2}\right)})={\mathrm{C}}_{\mathrm{a}\mathrm{é}\mathrm{r}\mathrm{o}}\left(t\right)\\ {\mathrm{I}}_{\left(\mathrm{1,2}\right)}{\ddot{\mathsf{\theta }}}_{\left(\mathrm{1,2}\right)}+{\mathrm{r}}_{\left(\mathrm{1,2}\right)}^{\mathrm{b}}{\mathrm{c}}_1\left(t\right){\mathsf{\delta }}_1\left(t\right)-{\mathrm{k}}^{\mathsf{\theta }}({\mathsf{\theta }}_{\left(\mathrm{1,1}\right)}-{\mathsf{\theta }}_{\left(\mathrm{1,2}\right)})+{\mathrm{k}}_1\left(t\right){\mathrm{r}}_{\left(\mathrm{1,2}\right)}^{\mathrm{b}}{\mathsf{\delta }}_1\left(t\right)=0\\ {\mathrm{I}}_{\left(\mathrm{2,1}\right)}{\ddot{\mathsf{\theta }}}_{\left(\mathrm{2,1}\right)}+{\mathrm{r}}_{\left(\mathrm{2,1}\right)}^{\mathrm{b}}{\mathrm{c}}_1\left(t\right){\mathsf{\delta }}_1\left(t\right)-{\mathrm{k}}^{\mathsf{\theta }}({\mathsf{\theta }}_{\left(\mathrm{2,1}\right)}-{\mathsf{\theta }}_{\left(\mathrm{2,2}\right)})+{\mathrm{k}}_1\left(t\right){\mathrm{r}}_{\left(\mathrm{2,1}\right)}^{\mathrm{b}}{\mathsf{\delta }}_1\left(t\right)=0\\ {\mathrm{I}}_{\left(\mathrm{2,2}\right)}{\ddot{\mathsf{\theta }}}_{\left(\mathrm{2,2}\right)}+{\mathrm{r}}_{\left(\mathrm{2,2}\right)}^{\mathrm{b}}{\mathrm{c}}_2\left(t\right){\mathsf{\delta }}_2\left(t\right)-{\mathrm{k}}^{\mathsf{\theta }}({\mathsf{\theta }}_{\left(\mathrm{2,1}\right)}-{\mathsf{\theta }}_{\left(\mathrm{2,2}\right)})+{\mathrm{k}}_2\left(t\right){\mathrm{r}}_{\left(\mathrm{2,2}\right)}^{\mathrm{b}}{\mathsf{\delta }}_2\left(t\right)=0\\ {\mathrm{I}}_{\left(\mathrm{3,1}\right)}{\ddot{\mathsf{\theta }}}_{\left(\mathrm{3,1}\right)}+{\mathrm{r}}_{\left(\mathrm{3,1}\right)}^{\mathrm{b}}{\mathrm{c}}_2\left(t\right){\mathsf{\delta }}_2\left(t\right)-{\mathrm{k}}^{\mathsf{\theta }}({\mathsf{\theta }}_{\left(\mathrm{2,2}\right)}-{\mathsf{\theta }}_{\left(\mathrm{3,1}\right)})+{\mathrm{k}}_2\left(t\right){\mathrm{r}}_{\left(\mathrm{3,1}\right)}^{\mathrm{b}}{\mathsf{\delta }}_2\left(t\right)=0\\ {\mathrm{I}}_{\mathrm{b}}{\ddot{\mathsf{\theta }}}_{\left(\mathrm{3,2}\right)}+{\mathrm{k}}^{\mathsf{\theta }}({\mathsf{\theta }}_{\left(\mathrm{3,1}\right)}-{\mathsf{\theta }}_{\left(\mathrm{3,2}\right)})=-{\mathrm{C}}_{\mathrm{r}}\left(t\right)\end{array}$$

(20)

4. Uncertainty Methods

4.1. Monte Carlo Simulation

The MCS technique, formally was introduced in 1949 [27], is the most widely used technique for determining the progression of uncertainties in random parameters. Indeed, this method can be used for any model. Moreover, the statistical results of MCS are realistic, so it is always used as a reference method. The MCS is based on sampling N random samples (or realizations) from the probability density function of the random variables. Each realization constitutes of a deterministic challenge that must be solved through a deterministic method. To achieve good accuracy, the MCS technique requires a large number of samples. Therefore, if the processing time of a single sample is very high, the methodology will become very expensive in terms of computational time. Consequently, its use in complex problems is avoided [27].

In this work, simulations by MCS have been performed following the description in (Figure 2). It started by the identification of the uncertain parameters, then the random generation of this parameters in the model and the calculation of their responses. Finally, the calculation of their mean value and standard deviation. The execution of the simulation involving 100,000 Monte Carlo samples required approximately 48 hours of computational time. The resulting dataset is utilized as a benchmark solution for validating the outcomes obtained via the Polynomial Chaos Expansion (PCE) technique. It is important to note that the Monte Carlo Simulation (MCS) method is characterized by slow convergence behavior and typically necessitates a large number of realizations to achieve statistically significant results, with 100,000 samples being a common requirement.

4.2. Polynomial Chaos

This segment introduces a contemporary procedural framework grounded in the projection of Principal Components (PC). The technique involves projecting targeted stochastic resolutions utilizing orthogonal polynomial functions with orthonormal Gaussian variables [28]. The foundational polynomial’s characteristics are leveraged to derive a linear equation system via solution projection, resulting in a polynomial basis expansion for the resolution. This expansion facilitates the computation of the probabilistic solution’s statistical moments. By implementing this methodology, the dynamic response of the mechanical system is numerically evaluated.

The generalization of the PC method establishes can be resolved into a series of convergent polynomial functions (mean square terms) as follows:

$${\mathrm{V}}_{\mathrm{i}}\left(\mathrm{x},\mathsf{\xi }\right)={\sum }_{\mathrm{j}=0}^{\mathrm{\infty }}{\bar{\mathsf{\nu }}}_{\mathrm{ij}}{\mathsf{\varphi }}_{\mathrm{j}}\left(\mathsf{\xi }\right),\mathsf{\xi }=\left({\mathsf{\xi }}_1,{\mathsf{\xi }}_2,...\right)$$

(21)

Where x is the nodal displacement vector.

In which ${\mathsf{\varphi }}_{\mathrm{j}}\left(\mathsf{\xi }\right)$denote an orthogonal polynomial within independent random factors $\mathsf{\xi }$ with a given random density of probability$\mathrm{W}\left(\mathsf{\xi }\right)$.

$$⟨{\mathsf{\varphi }}_{\mathrm{i}},{\mathsf{\varphi }}_{\mathrm{j}}⟩=\int {\mathsf{\varphi }}_{\mathrm{i}}{\mathsf{\varphi }}_{\mathrm{j}}\mathrm{W}\left(\mathsf{\xi }\right)\mathrm{d}\left(\mathsf{\xi }\right)=\left\{\begin{array}{c}0\\ ⟨{\mathsf{\varphi }}_{\mathrm{i}},{\mathsf{\varphi }}_{\mathrm{j}}⟩\end{array}\begin{array}{c}\mathrm{i}\mathrm{f}\mathrm{i}\ne j\\ \mathrm{i}\mathrm{f}\mathrm{i}=\mathrm{j}\end{array}\right.$$

(22)

According to Askey’s scheme, Legendre polynomials are the most relevant for treating uncertainties. The uncertainty of the random eigenvalues is given by:

$${\mathrm{N}}_{\mathrm{p}}={\displaystyle \frac{\left(\mathrm{p}+\mathrm{r}\right)!}{\mathrm{P}!\mathrm{r}!}}-1$$

(23)

To determine the coefficients ${\bar{\mathrm{V}}}_{\mathrm{ij}}$, the non-intrusive approach was used. This approach is particularly efficient, since it only requires simulations performed on specific models of the random variables [29]. It consists in reducing to the minimum of the gap$\epsilon$produced through the solution of the random model and its approximation in the generalized polynomial chaos basis.

$$\mathsf{\epsilon }={\sum }_{\mathrm{k}=1}^{\mathrm{Q}}\left[{\mathrm{V}}_{\mathrm{i}}\left(\mathrm{x},{\mathsf{\xi }}^{\left(\mathrm{k}\right)}\right)-{\sum }_{\mathrm{j}=0}^{\mathrm{N}\mathrm{p}}{\bar{\mathrm{V}}}_{\mathrm{i}\mathrm{j}}\left(\mathrm{x}\right){\mathrm{f}}_{\mathrm{j}}\left({\mathsf{\xi }}^{\left(\mathrm{k}\right)}\right)\right]$$

(24)

In efforts to minimize vibrational amplitudes within the transmission mechanism, prior studies have demonstrated that the dynamic response of the system is influenced by a multitude of mechano-geometric parameters. These factors introduce variability that can stem from various sources, including inherent design and manufacturing tolerances, as well as external perturbations arising from environmental interactions. Moreover, during the system design phase, many parameters must be assigned uncertainty intervals to account for stochastic variations. This raises critical questions regarding the system’s dynamic displacement, dimensional stability, and operational longevity, as well as the validity of the defined uncertainty bounds. Specifically, it remains an open challenge whether a single uncertainty interval can sufficiently capture the dynamic behavior of a system governed by multiple interdependent uncertain parameters.

In the coming section, the physical parameters of mass, moment of inertia, the bearing’s bending stiffness, the traction-compression stiffness, the damping coefficient and the wind speeds for the two-stage gear box will be defined as a random variable using MC and PC techniques according to a procedure described in the flowchart of the numerical model in Figure 4. This study will allow us to examine the standard deviation of these parameters taking into account multivariate uncertainty.

5. Numerical Simulation

The numerical simulation process in MATLAB initiates with the assembly of displacement matrices corresponding to the degrees of freedom defined in Equation (21). To account for parametric uncertainties, variations in mass, moment of inertia, bearing bending stiffness, traction-compression stiffness, damping coefficients, and wind speeds are incorporated into the model. These stochastic parameters are then processed using the polynomial chaos expansion (PCE) methodology, which reconstructs the matrices by leveraging their mean values and prescribed standard deviations. This approach enables a probabilistic representation of the system’s response, ensuring robustness against the inherent variability of the input parameters.

The numerical simulation in this study is a combination of two approaches, which are presented in a coupled form as shown in (Figure 3). The simulation begins with identifying the uncertain parameters and their standard deviation. The MCS is then performed in three main steps: randomly generating uncertain parameters, propagating them through the model to calculate responses using a deterministic model, and finally calculating their mean value and standard deviation. Once the MCS is complete, the PC calculation begins with three main steps: generating the model, calculating the polynomial coefficients of the chaos, and calculating the probabilistic distribution of the mean values and standard deviations. The two approaches converge and the results are compared.

As part of this study, the different system parameters used in the simulation are presented in Table 1.

In the following simulations, uncertainty approaches are implemented to analyze the dynamic behavior of the system in the presence of uncertainties in the design parameters. A case study of a coupled form was analyzed to evaluate the effect of such uncertainty on the dynamic response of a two-stage gear transmission using the mentioned approaches (Figure 3). The dynamic response of the first bearing is studied in both directions, x and y. The simulation parameters used are shown in Table 1. The numerical results obtained through the PC methodology are evaluated and validated through a MCS comparison including up to 100 000 simulations.

The random variables considered below are respectively: the mass of the gear wheels ${\mathrm{m}}_{\mathrm{g}}$, the damping coefficients $\mathsf{\eta }=10^{-2}$,$\mathsf{\mu }=10^{-5}$, the inertia of the wheel ${\mathrm{I}}_{\mathrm{w}}$ , the inertia of the pinion ${\mathrm{I}}_{\mathrm{p}}$, and the driving inertia ${\mathrm{I}}_{\mathrm{m}}$.

$${\mathrm{m}}_{\mathrm{g}}={\mathrm{m}}_{\mathrm{g},0}+{\mathsf{\sigma }}_{{\mathrm{m}}_{\mathrm{g}}}\mathsf{\xi }$$

(25)

$$\mathsf{\eta }={\mathsf{\eta }}_0+{\mathsf{\sigma }}_{\mathsf{\eta }}\mathsf{\xi }$$

(26)

$$\mathsf{\mu }={\mathsf{\mu }}_0+{\mathsf{\sigma }}_{\mathsf{\mu }}\mathsf{\xi }$$

(27)

$${\mathrm{I}}_{\mathrm{w}}={\mathrm{I}}_{\mathrm{w},0}+{\mathsf{\sigma }}_{{\mathrm{I}}_{\mathrm{w}}}\mathsf{\xi }$$

(28)

$${\mathrm{I}}_{\mathrm{P}}={\mathrm{I}}_{\mathrm{P},0}+{\mathsf{\sigma }}_{{\mathrm{I}}_{\mathrm{P}}}\mathsf{\xi }$$

(29)

$${\mathrm{I}}_{\mathrm{d}}={\mathrm{I}}_{\mathrm{d}}+{\mathsf{\sigma }}_{{\mathrm{I}}_{\mathrm{d}}}\mathsf{\xi }$$

(30)

Where $\mathsf{\xi }$ represents the Gaussian random variable with zero mean; ${\mathrm{m}}_{\mathrm{g},0}$,${\mathrm{I}}_{\mathrm{w},0}$, ${\mathrm{I}}_{\mathrm{P},0}$, ${\mathsf{\eta }}_0$, ${\mathsf{\mu }}_0$, ${\mathrm{I}}_{\mathrm{d},0}$ denotes respectively the mean value.${\mathsf{\sigma }}_{{\mathrm{m}}_{\mathrm{g}}}$,${\mathsf{\sigma }}_{{\mathrm{I}}_{\mathrm{d}}}$, ${\mathsf{\sigma }}_{\mathsf{\eta }}$, ${\mathsf{\sigma }}_{\mathsf{\mu }}$, ${\mathsf{\sigma }}_{{\mathrm{I}}_{\mathrm{p}}}$, ${\mathsf{\sigma }}_{{\mathrm{I}}_{\mathrm{w}}}$ are the standard deviation of these parameters.

To evaluate the effect of the standard deviation and the dynamic response of the studied system, the standard deviations used in the simulations is$\mathsf{\sigma }=2\%$,

6. Results

6.1. Impact of Uncertain Gear Mass

Figure 4 and Figure 5 shows the mean values of the linear displacement along x-axis and along y-axis that are presented by three curves: the desicise graph and the response of the MCS simulation and the mean-value of the PC simulation with ${\mathsf{\sigma }}_{\mathrm{m}}=2\%$. By exploiting an uncertain masse in the simulation, it should be noted that the average values of ${\mathrm{x}}_1\left(\mathrm{t}\right)$ and ${\mathrm{y}}_1\left(\mathrm{t}\right)$ remain relatively stable, with the MCS and PC curves closely following the decisive graph, suggesting a limited influence of uncertainties on the means value.

Figure 4. Average value of the displacement ${\mathrm{x}}_1\left(\mathrm{t}\right)$ with ${\mathsf{\sigma }}_{\mathrm{m}}=2\mathrm{\%}.$

Figure 4. Average value of the displacement ${\mathrm{x}}_1\left(\mathrm{t}\right)$ with ${\mathsf{\sigma }}_{\mathrm{m}}=2\mathrm{\%}.$

According to the values given by the standard deviation, (Figure 6 and Figure 7), the influence of the uncertain variable is more pronounced on ${\mathrm{y}}_1\left(\mathrm{t}\right)$ than on ${\mathrm{x}}_1\left(\mathrm{t}\right)$.For ${\mathrm{x}}_1\left(\mathrm{t}\right)$, the standard deviation gradually increases to about $5\times {10}^{-5}$, indicating moderate variability due to uncertainties. In contrast, for ${\mathrm{y}}_1\left(\mathrm{t}\right)$, the standard deviation reaches approximately $\mathrm{1,4}\times {10}^{-5}$, showing much greater variability and a more significant influence of uncertainties.

6.2. Impact of Damping Coefficient

The comparison between the decisive graph and the two other simulations (Monte Carlo and Chaos Polynomial), (Figure 8 and Figure 9), shows that, despite the uncertainty in the damping coefficients, both methods effectively reproduce the overall trend of the decisive curve. The simulation curves exhibit a slight attenuation of the oscillations, indicating that the uncertainty in the damping coefficients primarily impacts local variations (peaks and troughs) without significantly affecting the overall mean value. In summary, while the uncertainty in the coefficients slightly reduces the precision of the oscillations, it does not alter the mean value of the system’s response.

According to the values given by the standard deviation (Figure 10 and Figure 11), the influence of the uncertain variable is more pronounced on $y_1\left(t\right)$ than on $x_1\left(t\right)$. For displacement along the x-axis, the standard deviation gradually increases up to approximately $5\times {10}^{-5}$, indicating moderate variability. In contrast, along the y-axis, the standard deviation reaches about$\mathrm{1,4}\times {10}^{-4}$, showing much greater variability and a significant influence of uncertainties.

6.3. Impact of the Pinion’s Moment of Inertia

Figure 12 and Figure 13 show the mean values of the linear displacement along x-axis and along y-axis that are presented by three curves: the desicise graph and the response of the MCS simulation and the mean-value of the PC simulation with ${\mathsf{\sigma }}_{\mathrm{m}}=2\%$. By exploiting an uncertain masse in the simulation, it should be noted that the average values of ${\mathrm{x}}_1\left(\mathrm{t}\right)$ and ${\mathrm{y}}_1\left(\mathrm{t}\right)$ remain relatively stable, with the MCS and PC curves closely following the decisive graph, suggesting a limited influence of uncertainties on the means value.

In the case where the pinion’s moment of inertia is considered uncertain, the results are presented in (Figure 14 and Figure 15). These results show equal values,$3\times 10^{-6}$, in terms of standard deviation levels. This would imply that the effect due to the uncertainty of the pinion’s moment of inertia is evenly distributed between the two axes. Such a situation would indicate that this uncertainty similarly affects displacements along both the x and y axes.

6.4. Impact of the Wheel’s Uncertain Inertia

By comparing the curves (Figure 16 and Figure 17), it is evident that the uncertainty in the wheel inertia has a more significant impact on the mean displacement values along the x-axis compared to the y-axis. This effect is visible through the amplitude of oscillation and the dispersion of results, particularly in the Monte Carlo simulation and the chaos polynomial simulation.

According to the values given by the standard deviation (Figure 18 and Figure 19), the influence of the uncertain variable is more pronounced on ${\mathrm{y}}_1\left(\mathrm{t}\right)$ than on ${\mathrm{x}}_1\left(\mathrm{t}\right)$. For displacement along the x-axis, the standard deviation gradually increases up to approximately $7\times 10^{-5}$, indicating moderate variability. In contrast, along the y-axis, the standard deviation reaches about$1.8\times 10^{-4}$, showing much greater variability and a significant influence of uncertainties.

6.5. Impact of Uncertain Driving Inertia

Figure 20 and Figure 21 show the mean values of the linear displacement along x-axis and along y-axis that are presented by three curves: the desicise graph and the response of the MCS simulation and the mean-value of the PC simulation with ${\mathsf{\sigma }}_{{\mathrm{I}}_{\mathrm{d}}}=2\%$. By exploiting an uncertain masse in the simulation, it should be noted that the average values of ${\mathrm{x}}_1\left(\mathrm{t}\right)$ and ${\mathrm{y}}_1\left(\mathrm{t}\right)$ remain relatively stable, with the MCS and PC curves closely following the decisive graph, suggesting a limited influence of uncertainties on the means value.

According to the values given by the standard deviation Figure 22 and Figure 23, the influence of the uncertain variable is more pronounced on ${\mathrm{y}}_1\left(\mathrm{t}\right)$ than on ${\mathrm{x}}_1\left(\mathrm{t}\right)$. For displacement along the x-axis, the standard deviation gradually increases up to approximately$2.5\times 10^{-5}$, indicating moderate variability. In contrast, along the y-axis, the standard deviation reaches about$7\times 10^{-5}$, showing much greater variability and a significant influence of uncertainties.

Table 2 evaluates the impact of uncertainties on the system, where the randomness of each parameter affects the system’s dynamic behavior along the (x) and (y) directions. The mass of the gearwheels has a significant effect on the (y) direction, likely influencing the alignment or vibration of the gears. The damping coefficients strongly impact the first two degrees of freedom ($x_1$ and $y_1$), which could relate to the primary stage of the gear system. The pinion inertia, representing the smaller gear, has a weak influence on both directions, while the wheel inertia, associated with the larger gear, shows a strong impact on both (x) and (y). The drive inertia, linked to the input shaft, has a weak effect on the (x) direction but a strong effect on the (y) direction, potentially affecting the system’s torsional stability. These insights are critical for optimizing the performance and reliability of the two-stage gear system, particularly in applications requiring precise motion control and minimal vibration.

7. Conclusion

This study deals with the dynamic behavior of a two-stage gearing system by integrating an uncertainty on the design parameters. This is performed through two approaches, the basic Monte Carlo tool which presents the most reliable results, and the probabilistic Polynomial Chaos methodology by applying three different standard deviation. A developed coupled model that takes into account the interaction between five parameters, namely gear mass, damping coefficient, pinion’s moment of inertia, wheel’s inertia, driving inertia. The main results of the present study show that polynomial chaos can be considered as an efficient tool for the uncertainty dispersions in the dynamic behavior study of a two-stage wind turbine gear system.