Fundamental winds are variations in the average wind speed in a wind field that persist throughout the operation of a wind turbine and significantly affect the transmission of rated power to the system. Simulation of the base wind model does not require consideration of the actual distribution of wind speeds. Typically, the fundamental wind speed is considered to be the average wind speed acting on the wind turbine and can be defined as a specific constant value as shown in Equation 1: V ¯ = a a  is  a  constant (1)

Random wind refers to the fluctuations in wind speed, denoted by the random wind speed V s , which can be described using a mathematical model [33] as shown in Equation 2: V s = V smax A − 1 , 1 cos ω v + ϕ v (2)

In the formula: V smax represents the maximum value of random wind, while A − 1 , 1 is a uniformly distributed random number ranging from −1 to 1; ω v denotes the average amplitude of wind speed fluctuations, typically ranging between 0.5 π and 2 π ; ϕ v is a uniformly distributed random quantity ranging from 0 to 2 π . Therefore, the comprehensive wind speed model can be represented as the superposition of basic wind and random wind as shown in Equation 3: V = V s + V ¯ = V s m a x A − 1 , 1 cos ω v + ϕ v + a (3)

The input and output torque of the wind turbine gear transmission system can be determined using aerodynamic principles and the Bates formula [32] as shown in Equations 4, 5: T i n = 1 2 ρ a i r π r y 2 V 3 C p / ω y (4) T o u t = T i n / Z s i m (5)

In the formula: Z s i m represents the total transmission ratio of the system; ρ a i r denotes the air density; r y stands for the radius of the wind turbine blade; V signifies the comprehensive wind speed, and C p indicates the wind energy utilization coefficient.

Tooth side clearance in gear manufacturing and assembly causes the meshing trajectory to change during operation. Repeated load effects impact meshing accuracy and gear lifespan. Thus, tooth side clearance is crucial in analyzing the nonlinear dynamic characteristics of gears.

Assuming that the expression for the relative displacement of gears along the meshing line is as shown in Equation 6: x av = r 1 θ 1 + r 2 θ 2 + e t (6)

In the formula: r 1 is the base circle radius of the driving wheel; r 2 is the base circle radius of the driven wheel; θ 1 is the rotation angle of the driving wheel; θ 2 is the rotation angle of the driven wheel. As illustrated in Figure 2, assuming b n represents the initial half tooth side clearance of the gear, the tooth side clearance function is [34] shown in Equation 7: f δ = δ − b n   δ > b n 0   δ ≤ b n δ + b n   δ < − b n (7)

During gear operation, factors like tooth friction can lead to power loss within the gear system. In a stable working state, well-lubricated gears can assume a uniform temperature field. Δ T is the change in temperature rise, and the meshing angle, indexing radius, and tooth thickness after temperature change are shown in Equations 8–12: cos ⁡ α ′ = a 0 ⁡ cos ⁡ α / d 1 + d 2 + Δ d 1 + Δ d 2 (8) r 1 ′ = r 1 ⁡ cos ⁡ α / ⁡ cos ⁡ α ′ (9) r 2 ′ = r 2 ⁡ cos ⁡ α / ⁡ cos ⁡ α ′ (10) S 1 ′ = S 1 r 1 ′ r 1 − 2 r 1 ′ i n v α ′ − i n v α (11) S 2 ′ = S 2 r 2 ′ r 2 − 2 r 2 ′ i n v α ′ − i n v α (12)

In the formula: d 1 and d 2 are the diameters of the indexing circles for the driving wheel and the driven wheel, respectively; temperature influences the indexing circle modifications of the two gears, designated as Δ d 1 and Δ d 2 , respectively; a 0 is the initial center distance; r 1 and r 2 denote the pitch radii of the driving wheel and the driven wheel, respectively; α is the meshing angle; tooth thickness of the driving wheel and driven wheel are denoted by S 1 and S 2 , respectively. The amount of change in gear backlash caused by temperature changes is shown in Equation 13: Δ b a = Δ T λ 1 ⁢ S 1 ′ + Δ T λ 2 ⁢ S c ′ = Δ T λ 1 ⁢ S 1 r 1 − 2 i n v α ′ − i n v α ⋅ r 1 ⁢ cos ⁡ α cos ⁡ α ′ − Δ T λ 2 ⁢ S 2 r 2 − 2 i n v α ′ − i n v α ⋅ r 2 ⁢ cos ⁡ α cos ⁡ α ′ (13)

In the formula: thermal expansion coefficients of the two gears are designated as μ 1 and μ 2 , respectively. This article believes that all gear materials are consistent. Therefore, λ = λ 1 = λ 2 , the change in the gear backlash caused by the change in the indexing circle is shown in Equation 14: Δ b b = Δ d 1 2 + Δ d 2 2 sin ⁡ α = Δ T λ a 0 ⁡ sin ⁡ α = Δ T λ m z 1 + z 2 sin ⁡ α / 2 (14)

In summary, the change in clearance caused by temperature in gears is Δ b = Δ b a + Δ b b , and the temperature change causes the tooth side clearance to become b ′ = b n − Δ b . The tooth side clearance function under the influence of temperature is shown in Equation 15: f δ = δ − b ′   δ > b ′ 0   δ ≤ b ′ δ + b ′   δ < − b ′ (15)

In the absence of lubricating oil during gear operation, teeth come in direct contact, leading to heightened friction and wear. Introducing lubricating oil creates a film between contact surfaces, preventing direct gear contact and significantly reducing wear. Temperature fluctuations affect the lubrication efficacy of oil by changing its viscosity, a key parameter of oil performance. According to Sutherland equation, the expression between lubricant viscosity and temperature change can be derived [35] shown in Equation 16: η Δ T = η 0 e c 1 1 Δ T + 273 − 1 293 + c 2 1 Δ T + 273 − 1 293 2 (16)

In the formula: η Δ T is the viscosity of the lubricating oil (under temperature changes); η 0 is the viscosity of the lubricating oil at 20°C.

A simplified diagram in Figure 3A illustrates the engagement between two gears. Point c marks their meshing interface, where the velocity can be mathematically expressed as Equation 17: V o 1 = ω 1 L o 1 c V o 2 = ω 2 L o 2 c (17)

The sliding speed of meshing point c is: V c = V o 1 ⁡ sin ⁡ ψ 1 + V o 2 ⁡ sin ⁡ ψ 2 (18)

Among them, ψ 1 = ω 1 t + θ 1 , ψ 2 = ω 2 t + θ 2 , ω 1 signifies the angular displacement of the driving wheel, whereas ω 2 corresponds to that of the driven wheel, respectively; θ 1 and θ 2 denote the minute displacements experienced by the driving wheel and the driven wheel, respectively. Observing the geometric relationships depicted in the figure, can deduce the friction arms associated with the two gears: L A c = L o 1 c sin ⁡ ψ 1 L B c = L o 2 c sin ⁡ ψ 2 (19)

Utilizing the comprehensive Equations 18, 19, the sliding velocity at the meshing point c can be using mathematical Equation 20: V c = ω 1 L A c + ω 2 L B c (20)

The friction force arm can be written as shown in Equation 21: L A c = r 1 + r 2 tan ⁡ α − r 1 2 + r 2 2 + r 1 ω 1 t L B c = r 1 2 + r 2 2 − r 1 ω 1 t (21)

The frictional force and torque experienced by the two gears can be calculated as follows, taking into account the relevant factors and Equations 22, 23: F f = μ γ c F n (22) T f 1 = F f L A c , T f 2 = F f L B c (23)

In the formula: μ is the friction coefficient; F n gear meshing force; γ c is the directional coefficient. Equation 24 is the directional coefficient function. γ c = 1   V c > o 0   V c = o − 1   V c < 0 (24)

The friction coefficient on gear surfaces varies during meshing due to the presence of lubricating oil. Figure 3B depicts the instantaneous curvature radius of a pair of meshing gears. Consequently, the friction coefficient of the tooth surface can be determined based on this geometric characteristic and the properties of the lubricating oil present during the meshing process [36]; as shown in Equation 25: μ = 0.872 b a L ⁡ ln η \ / η 0 τ η F i + 1.74 τ η 2 b a L F i ln ⁡ 1.2 τ η h η 2 k η η 0 1 + 9.6 ζ 1 / 2 (25)

Among them: ζ = 4 k l π h c / 1 / R 1 + 1 / R 2 ⁢ 4 F n 1 − ν 2 E k c 1 / R 1 + 1 / R 2 ρ c 2 b a L R 1 ω 1 sin ⁡ α + 2 l / r a + R 2 ω 2 sin ⁡ α + 2 l / r b 1 / 2 (26)

In the formula: τ η represents the shear stress experienced by the lubricating oil; L is the width of the action line on the tooth surface; b a is the contact width of the tooth surface; R 1 and R 2 represent the curvature radii of the driving wheel and the driven wheel, respectively; k c and k η signify the thermal conductivity of the gear material and the lubricating oil, respectively; h η is the thickness of the contact oil film; E is the elastic modulus; ν is the Poisson’s ratio; The density and specific heat capacity of the gear material are designated by ρ and c , respectively.

During meshing, gears can be simplified as a cantilever beam consisting of a combination of trapezoidal and rectangular shapes. Thus, when shifting the meshing point from the gear crest to its root, the meshing stiffness undergoes periodic variations: k = F N / Δ δ (27)

In the formula: F N is the meshing force; Δ δ represents elastic deformation. Expanding Equation 27 into a Fourier series yields the time-domain expression as follows Equation 28: k t = k m + ∑ i = 1 ∞ k 2 i ⁡ cos n ω i t + k 2 i + 1 ⁡ cos n ω i t (28)

In the formula: k m is the average meshing stiffness; k 2 i and k 2 i + 1 are Fourier coefficients, respectively; ω i is the meshing frequency. Due to the small difference in high-order meshing errors, only the first order is considered as shown in Equation 29: k t = k m + k A ⁡ cos ω i t + ϕ (29)

In the formula: k A represents the amplitude of the variation in meshing stiffness, while ϕ denotes the initial phase angle. Figure 4 shows the meshing stiffness curves of each meshing pair in the gear system.

The meshing damping of gear pairs can be obtained through empirical formulas as shown in Equation 30: c m = 2 ξ k 1 / M 1 + 1 / M 2 (30)

In the formula: Respectively, M 1 and M 2 signify the equivalent masses of the main and driven wheels; ξ is the meshing damping ratio; k is the meshing stiffness of two gears. During gear meshing, oil film damping arises between the teeth due to the lubrication provided by the oil present. The coupling of gear meshing damping c m and oil film damping c η results in the formation of an equivalent damping, denoted as c 0 . Its expression is Equation 31: c 0 = c η c m c η + c m (31)

The Reynolds equation can be expressed as: ∂ ∂ ξ h η 3 ∂ p ∂ ξ = 12 η ∂ h η ∂ t (32)

Solve Equation 32 to derive the expression Equation 33 for the oil film damping force: F 0 = c 0 ∂ h η ∂ t = − 8 a 3 η l h η 3 ∂ h η ∂ t (33)

In the formula: c 0 = − 8 a 3 η l / h η 3 ; a = 8 F N R 1 − ν 2 / B E π ; l represents the oil film thickness along the contact surface; ∂ h η / ∂ t is the change in lubricating oil; R denotes the comprehensive curvature radius of the gear; B is the tooth width.

The representation of meshing error can be achieved through the utilization of Fourier series as shown in Equation 34: e t = e m + ∑ i = 1 ∞ e 2 i ⁡ cos n ω i t + e 2 i + 1 ⁡ cos n ω i t (34)

In the formula: e m signifies the mean value of meshing error. Due to the small difference in the order of high-order meshing errors, only the first order needs to be considered, e m = 0 , Simplify the meshing error as follows Equation 35: e t = e a ⁡ cos ω i t + φ 0 (35)

In the formula: e a is the amplitude of the error.

Create a simplified pure torsional model for a wind power planetary-face gear system utilizing the lumped mass approach, depicted in Figure 5, with the following assumptions: All gears are installed as standard; Assuming rigid body behavior for all gears, the interacting pair is modeled as a spring-damping system for simplicity; The friction and lubrication properties are considered uniform across all system components; Lateral vibration displacement of gears is disregarded in the model; Without considering the stiffness of the shaft, the sun gear s and the face gear m 1 are considered as one component, therefore the angular acceleration of the two gears is considered equal.

The degrees of freedom of the planetary-facing gear system were first determined, and for the system studied in this paper, only the torsional direction displacement was investigated. Each gear assembly is reduced to a centralized mass point with the mass distribution concentrated by the center of mass. For the rotational degrees of freedom, the gears are considered as equivalent rotational inertia in order to simplify to a mass-spring-damping system with single or multiple degrees of freedom, and secondly, each gear meshing pair is considered as a spring and damping unit. The springs represent the effect of the gears on the meshing stiffness, while the damping represents the energy dissipation characteristics of the system, and finally the dynamical equations are developed according to Newton’s second method. In this paper, only the torsional degrees of freedom are considered. The elastic deformation of each gear along the meshing line direction is: δ s p i = r s θ s + r p i θ p i − r p i + r s r c r c θ c − e s p i t (36) δ r p i = − r p i θ p i − r r − r p i r c r c θ c − e s p i t (37) δ m 1 z i = r z i θ z i − r m 1 θ m 1 − e m 1 z i t (38) δ m 2 z i = r z i θ z i − r m 2 θ m 2 − e m 2 z i t (39)

Establish system dynamics equations using the Newton Euler method: I c ⁢ θ ¨ c + ∑ i = 1 n k r p i ⁢ f ⁡ δ r p i + c r p i ⁢ δ ˙ r p i ⁢ r c + ∑ i = 1 n k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r c = T i n I p i ⁢ θ ¨ p i + k r p i ⁢ f ⁡ δ r p i + c r p i ⁢ δ ˙ r p i ⁢ r p i − k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r p i = − T f r p i − T f s p i I r ⁢ θ ¨ r + ∑ i = 1 n k r p i ⁢ f ⁡ δ r p i + c r p i ⁢ δ ˙ r p i ⁢ r r = ∑ i = 1 n − T f r p i I s ⁢ θ ¨ s − ∑ i = 1 n k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r s = ∑ i = 1 n − T f s p i I m 1 ⁢ θ ¨ m 1 + ∑ i = 1 n k z i m 1 ⁢ f ⁡ δ z i m 1 + c z i m 1 ⁢ δ ˙ z i m 1 ⁢ r m 1 = − ∑ i = 1 m T f z i m 1 I z i ⁢ θ ¨ z i − k z i m 1 ⁢ f ⁡ δ z i m 1 + c z i m 1 ⁢ δ ˙ z i m 1 ⁢ r z i + k z i m 2 ⁢ f ⁡ δ z i m 2 + c z i m 2 ⁢ δ ˙ z i m 2 ⁢ r z i = − ∑ i = 1 m T f z i m 1 − ∑ i = 1 m T f z i m 2 − T o u t I m 2 ⁢ θ ¨ m 2 − ∑ i = 1 n k z i m 2 ⁢ f ⁡ δ z i m 2 + c z i m 2 ⁢ δ ˙ z i m 2 ⁢ r m 2 = − ∑ i = 1 m T f z i m 2 − T o u t ′ (40)

In the formula: I c . I p i , I r , I m 1 , I m 2 and I z i are the moment of inertia of the planetary carrier, planetary gear, inner ring gear, face gear m 1 , face gear m 2 , and spur gear, respectively; k r p i , k s p i , k z i m 1 and k z i m 2 are the meshing stiffness of the inner ring gear planetary gear, sun gear planetary gear, spur gear and face gear m 1 , and spur gear and face gear m 2 , respectively; c r p i , c s p i , c z i m 1 , c z i m 2 refers to the meshing damping of inner ring gear planetary gear, sun gear planetary gear, spur gear and face gear m 1 , and spur gear and face gear m 2 ; T f r p i , T f s p i , T f z i m 1 and T f z i m 2 are the friction torques of the meshing pairs of inner ring gear planetary gear, sun gear planetary gear, spur gear and face gear m 1 , and spur gear and face gear m 2 , respectively; T i n is the input torque of the system; T o u t and T o u t ′ are both the output torque of the system; r c , r p i , r r , r s , r m 1 and r m 2 are the base circle radii of planetary carrier, planetary gear, inner ring gear, sun gear, face gear m 1 , and face gear m 2 , respectively; n and m represent the number of planetary gears and spur gears, respectively. By combining Equations 36–40, we can obtain: δ ¨ s p i = ∑ i = 1 n k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r s I s − ∑ i = 1 n T f s p i I s − k r p i ⁢ f ⁡ δ r p i + c r p i ⁢ δ ˙ r p i ⁢ r p i I p i + k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r p i I p i − T f r p i I p i − T f s p i I p i − T i n I c + ∑ i = 1 n k r p i ⁢ f ⁡ δ r p i + c r p i ⁢ δ ˙ r p i ⁢ r c I c − ∑ i = 1 n k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r c I c − e ¨ s p i ⁢ t δ ¨ r p i = k r p i ⁢ f ⁡ δ r p i + c r p i ⁢ δ ˙ r p i ⁢ r p i I p i − k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r p i I p i + T f r p i I p i + T f s p i I p i − T i n I c + ∑ i = 1 n k r p i ⁢ f ⁡ δ r p i + c r p i ⁢ δ ˙ r p i ⁢ r c I c − ∑ i = 1 n k s p i ⁢ f ⁡ δ s p i + c s p i ⁢ δ ˙ s p i ⁢ r c I c − e ¨ s p i ⁢ t δ ¨ m 1 z i = k z i m 1 ⁢ f ⁡ δ z i m 1 + c z i m 1 ⁢ δ ˙ z i m 1 ⁢ r z i I z i − k z i m 2 ⁢ f ⁡ δ z i m 2 + c z i m 2 ⁢ δ ˙ z i m 2 ⁢ r z i I z i − ∑ i = 1 m T f z i m 1 I z i − ∑ i = 1 m T f z i m 2 I z i − T o u t I z i + ∑ i = 1 n k z i m 1 ⁢ f ⁡ δ z i m 1 + c z i m 1 ⁢ δ ˙ z i m 1 ⁢ r m 1 I m 1 + ∑ i = 1 m T f z i m 1 I m 1 − e ¨ m 1 z i ⁢ t δ ¨ m 2 z i = k z i m 1 ⁢ f ⁡ δ z i m 1 + c z i m 1 ⁢ δ ˙ z i m 1 ⁢ r z i I z i − k z i m 2 ⁢ f ⁡ δ z i m 2 + c z i m 2 ⁢ δ ˙ z i m 2 ⁢ r z i I z i − ∑ i = 1 m T f z i m 1 I z i − ∑ i = 1 m T f z i m 2 I z i − T o u t I z i − ∑ i = 1 n k z i m 2 ⁢ f ⁡ δ z i m 2 + c z i m 2 ⁢ δ ˙ z i m 2 ⁢ r m 2 I m 2 + ∑ i = 1 m T f z i m 2 I m 2 + T o u t ′ I m 2 − e ¨ m 2 z i ⁢ t (41)

In order to facilitate numerical solutions and improve the stability of the calculation process, it is necessary to dimensionless the dynamic equations. Define dimensionless time as t ¯ = ω d t , Where ω d is the nominal angular frequency of the system when dimensionless, ω d 1 = k s p i m s + m p i / m s m p i 1 / 2 , ω d 2 = k m z i m m + m z i / m m m z i 1 / 2 , and introduce displacement scale b c , The non-dimensional displacement, non-dimensional velocity, non-dimensional acceleration, and non-dimensional frequency of the system can be expressed as: X = δ / b c , X ˙ = δ ˙ / b c ω d , X ¨ = δ ¨ / b c ω d 2 , Δ = ω / ω d . By nondimensionalizing Equation 41, we can obtain Equations 42: X ¨ s p i = − 1 ω d 1 2 b c r s ∑ i = 1 n T f s p i + ∑ i = 1 n k s p i f δ s p i ω d 1 2 b c + c s p i δ ˙ s p i ω d 1 b c ⁢ r s − T i n ω d 1 2 b c r c − T f r p i ω d 1 2 b c r p i − T f s p i ω d 1 2 b c r p i + k s p i f δ s p i ω d 1 2 b c + c s p i δ ˙ s p i ω d 1 b c ⁢ r p i − k r p i f δ r p i ω d 1 2 b c + c r p i δ ˙ r p i ω d 1 b c ⁢ r p i + ∑ i = 1 n k r p i f δ r p i ω d 1 2 b c + c r p i δ ˙ r p i ω d 1 b c ⁢ r c + ∑ i = 1 n k s p i f δ s p i ω d 1 2 b c + c s p i δ ˙ s p i ω d 1 b c ⁢ r c − e ¨ s p i t ω d 1 2 b c X ¨ r p i = − T i n ω d 1 2 b c r c + T f r p i ω d 1 2 b c r p i + T f s p i ω d 1 2 b c r p i − k s p i f δ s p i ω d 1 2 b c + c s p i δ ˙ s p i ω d 1 b c ⁢ r p i + k r p i f δ r p i ω d 1 2 b c + c r p i δ ˙ r p i ω d 1 b c ⁢ r p i + ∑ i = 1 n k r p i f δ r p i ω d 1 2 b c + c r p i δ ˙ r p i ω d 1 b c ⁢ r c + ∑ i = 1 n k s p i f δ s p i ω d 1 2 b c + c s p i δ ˙ s p i ω d 1 b c ⁢ r c − e ¨ s p i t ω d 1 2 b c X ¨ m 1 z i = − 1 ω d 2 2 b c r z i ∑ i = 1 2 T f z i m 1 − 1 ω d 2 2 b c r z i ∑ i = 1 2 T f z i m 2 + T o u t ω d 2 2 b c r z i − k z i m 2 f δ z i m 2 ω d 2 2 b c + c z i m 2 δ ˙ z i m 2 ω d 2 b c ⁢ r z i + k z i m 1 f δ z i m 1 ω d 2 2 b c + c z i m 1 δ ˙ z i m 1 ω d 2 b c ⁢ r z i + 1 ω d 2 2 b c r m 1 ∑ i = 1 2 T f z i m 1 + ∑ i = 1 n k z i m 1 f δ z i m 1 ω d 2 2 b c + c z i m 1 δ ˙ z i m 1 ω d 2 b c ⁢ r m 1 − e ¨ m 1 z i t ω d 2 2 b c X ¨ m 2 z i = − 1 ω d 2 2 b c r z i ∑ i = 1 2 T f z i m 1 − 1 ω d 2 2 b c r z i ∑ i = 1 2 T f z i m 2 + T o u t ω d 2 2 b c r z i − k z i m 2 f δ z i m 2 ω d 2 2 b c + c z i m 2 δ ˙ z i m 2 ω d 2 b c ⁢ r z i + k z i m 1 f δ z i m 1 ω d 2 2 b c + c z i m 1 δ ˙ z i m 1 ω d 2 b c ⁢ r z i − ∑ i = 1 n k z i m 2 f δ z i m 2 ω d 2 2 b c + c z i m 2 δ ˙ z i m 2 ω d 2 b c ⁢ r m 2 + 1 ω d 2 2 b c r m 2 ∑ i = 1 2 T f z i m 2 + T o u t ′ ω d 2 2 b c r m 2 − e ¨ m 2 z i t ω d 2 2 b c ) (42)

Wind speed directly affects the performance of the entire power generation system during wind turbine operation. Establishing a random wind speed model is essential for simulating the nonlinear dynamic response of wind speed in the gear transmission system of the turbine. Figure 6 illustrates the dynamic simulation of random wind speed under average wind speeds of 5 m/s, 10 m/s, and 15 m/s.

In the random wind speed model, higher average wind speeds cause greater fluctuations and amplitudes in the wind speed curve. Using aerodynamic principles, the input torque curve of the wind turbine gear system can be modeled and simulated, as shown in Figure 7. To maintain the randomness of gear system excitation from the wind turbine input, the torque data from the first 10 seconds is used as the system’s input excitation. As average windspeed increases, the torque input from the wind turbine system also rises. Torque variations and torque wave momentum impact the gear system’s dynamic characteristics. Thus, analyzing the impact of average wind speed on the nonlinear dynamic characteristics of wind power gears is highly significant.

Analyzing the nonlinear dynamic characteristics of the system involves setting the system temperature rise to 60°C, with a time-varying mesh stiffness coefficient of 0.3, mesh damping ratio of 0.07, mesh error amplitude of 0.35, initial tooth side clearance of 30 microns, dimensionless excitation frequency of 1 [39]. The analysis spans average wind speeds from 0 m/s to 20 m/s.

Figure 8 shows the bifurcation of the gear system as the mean wind speed varies. The bifurcation diagram shows a clear bifurcation in the range of mean wind speed from 0 m/s to 17.5 m/s. Figure 9 shows the maximum Lyapunov exponent of the system as a function of the mean wind speed. In the above wind speed range, the maximum Lyapunov exponent stays below 0, indicating the presence of periodic motion in the system. When the mean wind speed is between 0 m/s and 5 m/s, the system shows a single-period motion. When the mean wind speed is between 5 m/s and 12 m/s, the system is in double periodic motion. Between 12 m/s and 17.5 m/s, the system is in a four-cycle state. When the mean wind speed increases between 17.5 m/s and 25 m/s, the bifurcation diagram does not show any obvious bifurcation and the maximum Lyapunov exponent of the system is constant greater than 0, which suggests that there is unstable chaotic motion in this range.

To further verify the global bifurcation diagram of the gear system with varying average wind speed, we analyzed phase trajectory diagrams, time-domain plots, Poincaré cross-sections, FFT spectrum diagrams 3D phase diagram and time–frequency diagram of the gear system at average wind speeds of 2 m/s, 10 m/s, 15 m/s, and 23 m/s.

Figure 10 depicts the time-domain plot, phase trajectory, FFT spectrum, and Poincaré cross-section of the system at an average wind speed of 2 m/s. The time domain plot exhibits clear periodicity, while the phase trajectory forms a closed curve encircling a circle. A peak is evident in the FFT spectrum, and a discrete point appears in the Poincaré cross-section, indicating that the gear system is in a single-cycle motion state at this time; The illustration in Figure 11 presents various diagrams including the time domain, phase trajectory, FFT spectrum, and Poincaré cross-section, all depicting the system’s behavior at an average wind velocity of 10 m/s. The time-domain graph displays a curve featuring two distinct periods, paralleled by the phase trajectory with two enclosed loops. The FFT spectrum exhibits two peaks, and the Poincaré cross-section identifies two separate points, suggesting the gear system has transitioned into a two-cycle motion state at present; Figure 12 shows the time domain diagram, phase trajectory diagram, FFT spectrum diagram, and Poincaré cross-section diagram of the system at an average wind speed of 15 m/s, indicating that the gear system has entered a four period motion state at this time; At an average wind speed of 23 m/s, the gear system demonstrates an unstable chaotic motion. Figure 13 illustrates the system’s time-domain, phase trajectory, FFT spectrum, and Poincaré cross-section at this velocity. Figures 14, 15 present time-frequency and 3D phase diagrams, respectively, showcasing the system’s behavior across varying wind speeds. It is evident that with rising wind speeds, the system progressively shifts from stability to instability.

The study reveals that the average wind speed significantly influences system characteristics. As it rises, the system transitions from periodic to chaotic motion, consequently, this wind power gear system is optimally suited for regions with an annual average wind speed below 17.5 m/s.

Variations in the radius of the wind turbine for a constant mean wind speed can change the aerodynamic performance and thus affect the input torque of the geared system. Therefore, wind turbine dimensions are valuable variables for exploring the nonlinear dynamics of the gear system. A time-varying meshing stiffness factor of 0.3, a damping ratio of 0.07, a meshing error magnitude of 0.35, an initial tooth side clearance of 30 μm, a dimensionless excitation frequency of 1, and an operating environment with an average wind speed of 10 m/s are set for a system temperature rise of 60°C. The nonlinear dynamic behavior of the system is analyzed for a range of wind turbine radii from 0 m to 30 m. The system is then subjected to a time-varying meshing stiffness coefficient of 0.3.

Figure 16 depicts the global bifurcation diagram of the gear system with the variation of wind turbine radius. Within the range of 0–12 m, there is a noticeable chaotic phenomenon in the graph. Figure 17 shows the maximum Lyapunov exponent as a function of radius; The index is always greater than 0 within the range of 0–12 m. When the radius is 14–16 m, the bifurcation diagram shows a three-period motion, and the maximum Lyapunov exponent is also negative. At a distance of 16–17.5 m, the bifurcation diagram shows a two-period motion with a constant negative exponent. Finally, stable single period motion occurs at 17.5–30 m, with a maximum negative Lyapunov exponent.

To validate the global bifurcation diagram of the gear system across varying wind turbine radii, phase trajectory, time domain, Poincaré section, FFT spectrum diagrams, 3D phase diagram and time–frequency diagram were analyzed specifically for radii of 5 m, 13 m, 15 m, 17 m, and 23 m.

Supplementary Figure S1 shows the time-domain, phase trajectory, FFT spectrum, and Poincare section of a 5 m wind turbine radius, indicating chaotic motion. Supplementary Figure S2 shows that when the radius of the wind turbine is 13 m, the system is in a six-cycle motion state. As the radius increases, Supplementary Figure S3 shows that the wind turbine system is in a three-cycle motion state at a radius of 15 m, Supplementary Figure S4 depicts the transition of the system to a two-cycle motion state at a radius of 17 m, and Supplementary Figure S5 depicts a stable single cycle motion at a radius of 23 m. Supplementary Figures S6, S7 show the time-frequency and 3D phase diagrams at different radii, respectively. It can be seen that increasing the radius of the wind turbine will cause the gear system to stabilize at a constant wind speed.

The findings indicate that the wind turbine radius significantly influences system behavior. As the radius increases, the system transitions from chaotic motion to periodic states (six-period, three-period, two-period, and single-period). The global bifurcation diagram illustrates increasing system stability with greater turbine radius. Systems with radius less than 12 m exhibit chaotic behavior, increasing gear collision severity and accelerating wear. Optimal operating speeds should be chosen to avoid prolonged exposure to chaotic conditions, recommending a wind turbine radius exceeding 12 m for system stability.

Supplementary Figure S8 shows the relationship between the maximum Lyapunov exponent of the gear system and the radius and average wind speed of the wind turbine. In the case of a small radius of a wind turbine, regardless of changes in wind speed, the gear system will maintain a chaotic state, indicating that the radius of the wind turbine is the main factor affecting system chaos in the same working environment. When the radius of the wind turbine is large and the wind speed is low, the gear system tends to stabilize.

Due to the complex operating conditions of most wind power gear systems, factors such as sunlight exposure and power loss during high-speed operation can lead to increased gear surface temperatures. Additionally, environmental temperatures vary across different regions, necessitating an analysis of temperature rise effects on the system. The gear system parameters include a time-varying meshing stiffness coefficient of 0.3, meshing damping ratio of 0.07, meshing error amplitude of 0.35, initial tooth side clearance of 30 microns, dimensionless excitation frequency of 1. The average wind speed in the region is set to 10 m/s.

Supplementary Figure S9 shows the bifurcation diagram of the gear system with temperature rise. Supplementary Figure S10 shows the maximum Lyapunov exponent of the corresponding gear system. From 0°C to 55°C, the bifurcation diagram indicates that the gear system is in unstable chaotic motion, and the maximum Lyapunov exponent is always positive. When the temperature rise range is from 55°C to 61°C, the bifurcation diagram indicates that the gear system is in double periodic motion, and the Lyapunov exponent remains negative. The temperature rise range is from 61°C to 100°C, and the bifurcation diagram indicates that the gear system is in single cycle motion. At this stage, the maximum Lyapunov exponent remains negative.

In order to accurately illustrate the motion characteristics, the phase trajectory diagram, time-domain diagram, Poincare section diagram, FFT spectrum diagram, 3D phase diagram, and time-frequency diagram of the gear system were analyzed at temperatures of 30°C, 59°C, and 90°C.

Supplementary Figure S11 displays the system’s time-domain, phase trajectory, FFT spectrum, and Poincaré cross-section diagrams with a 30°C rise in gear temperature, indicating unstable chaotic motion. Supplementary Figure S12 illustrates these diagrams with a 59°C temperature rise, showing the system in a two-cycle motion state. Supplementary Figure S13 depicts these diagrams with a 90°C temperature rise, indicating the system is in single-cycle motion. Supplementary Figure S14 shows the time-frequency diagram of the system at different temperature rises, indicating that as the system temperature increases, it transitions from chaotic motion to single period stable motion. Supplementary Figure S15 shows the 3D phase diagram of the system at different temperature rises, which can further illustrate its dynamic behavior.

The results indicate that as the temperature of the gear system increases, it transitions from chaotic motion to periodic motion, experiencing chaotic, bi periodic, and single periodic states. When the temperature rises to 55°C, the system is in a long-term chaotic state, and gear collisions increase.

Wind power gear system in the process of operation, usually add lubricant, lubricant will be added to form an oil film on the gear surface, which can not only avoid direct contact with the gear, but also reduce the friction of the gear surface in the process of work, slowing down the wear of the tooth surface. The lubrication performance of the lubricant is mainly reflected in the lubricant viscosity, so the effect of lubricant viscosity on the nonlinear dynamic characteristics of the gear system can be analyzed. The time-varying meshing stiffness coefficient of the gear system is set to be 0.3, the meshing damping ratio to be 0.07, the magnitude of the meshing error to be 0.35, the initial tooth flank clearance to be 30 μm, the dimensionless excitation frequency to be 1, the temperature rise of the gears to be 60°C, and the average wind speed of the wind machine’s working area to be set to be 10 m/s.

Supplementary Figure S16 displays the global bifurcation diagram of the gear system as the viscosity of the lubricating oil varies from 0.01 Pa·s to 0.1 Pa·s. Supplementary Figure S17 shows the maximum Lyapunov exponent of the gear system with changing viscosity of the lubricating oil. For viscosity values between 0.01 Pa·s and 0.016 Pa·s, the bifurcation diagram indicates the gear system is in unstable chaotic motion, with a consistently positive maximum Lyapunov exponent. As viscosity increases to between 0.016 Pa·s and 0.022 Pa·s, the system exhibits two-period motion, with a negative maximum Lyapunov exponent. When viscosity ranges from 0.022 Pa·s to 0.1 Pa·s, the system shows single-cycle motion, maintaining a negative maximum Lyapunov exponent throughout this range.

To accurately determine the motion characteristics, the gear system’s phase trajectory diagram, time domain chart, Poincaré section diagram, FFT spectrum diagram, 3D phase diagram and time–frequency diagram was analyzed at lubricating oil viscosities of 0.01 Pa·s, 0.02 Pa·s, and 0.05 Pa·s.

Supplementary Figure S18 depicts the time-domain diagram, phase trajectory, FFT spectrum, and Poincaré cross-section of a gear system lubricated with oil viscosity of 0.01 Pa·s, showing the system in an unstable chaotic motion state. Supplementary Figure S19 illustrates the time-domain diagram, phase trajectory diagram, FFT spectrum diagram, and Poincaré cross-section diagram of the system at 0.02 Pa·s viscosity, indicating the system has entered a two-cycle motion state. Supplementary Figure S20 shows the time-domain diagram, phase trajectory diagram, FFT spectrum diagram, and Poincaré cross-section diagram of the system at 0.05 Pa·s viscosity, where the system has entered a single-cycle motion state. Supplementary Figure S21 shows the time-frequency diagram of the system. It can be seen that as the viscosity of the lubricating oil increases, the system gradually becomes a single cycle motion and tends to stabilize. Under chaotic motion, the time-frequency diagram is a discrete spectrum with a certain width. Supplementary Figure S22 shows a three-dimensional phase diagram of the system when the viscosity of the lubricating oil changes. Under other unchanged conditions, using high viscosity lubricating oil can improve the stability of the system.

The results indicate that increasing the viscosity of lubricating oil will gradually transition the gear system from chaotic motion to periodic motion. It is recommended to choose lubricating oil with a viscosity higher than 0.016 Pa·s to reduce gear collision and prevent gear failure.

According to Equation 16, there is a relationship between the viscosity of lubricating oil and temperature rise. An increase in gear temperature rise will cause a decrease in the viscosity of lubricating oil. Therefore, it is necessary to consider the effects of temperature rise and lubricating oil viscosity on system characteristics simultaneously.

Supplementary Figure S23 shows the maximum Lyapunov exponent of the gear system as the viscosity of the lubricating oil and the temperature rise of the gear change simultaneously. It can be seen that when the viscosity of the lubricating oil is low ( η < 0.016 Pa·s), regardless of how the system temperature rises, the maximum Lyapunov value is always positive, indicating that the system has been in an unstable chaotic motion; When the viscosity of the lubricating oil is small (0.016 Pa·s < η < 0.05 Pa·s) and the temperature rise is low, the maximum Lyapunov value is less than 0. When the temperature rise is high, the maximum Lyapunov value is greater than 0. This indicates that in this range of lubricating oil, an increase in temperature rise will cause a decrease in the viscosity of the lubricating oil, making the system unstable; When the viscosity of the lubricating oil is high ( η > 0.05 Pa·s) and the temperature rise is high, although the viscosity of the lubricating oil will decrease, the expansion of the gear material causes a decrease in the backlash between the teeth, resulting in a negative maximum Lyapunov value of the system and a stable state of the gear system.

The results suggest that at low viscosity, the primary cause of gear chaos is the lubrication performance of the oil. At high viscosity, the main factor is the change in tooth side clearance due to gear temperature rise.