Abstract
The design optimization of wind turbines and their subsystems will make them competitive as an ideal alternative for energy. This paper proposed a design procedure for one of these subsystems, which is the Wind Turbine Drive-Train (WTDT). The design of the WTDT is based on the load assumptions and considered as the most significant parameter for increasing the efficiency of energy generation. In industry, these loads are supplemented by expert assumptions and manipulated to design the transmission elements. In contrary, in this work, the multibody system approach is used to estimate the static as well as dynamic loads based on the Lagrange multipliers. Lagrange multipliers are numerical parameters associated with the holonomic and nonholonomic constraints assigned in the drive-train model. The proposed scheme includes computational manipulations of kinematic constraints, mapping the generalized forces into Cartesian respective, and enactment of velocity-based constrains. Based on the dynamic model and the obtained forces, the design process of a planetary stage of WTDT is implemented with trade-off’s optimization in terms of gearing parameters. A wind turbine of 1.4 megawatts is introduced as an evaluation study of the proposed procedure, in which the main advantage is the systematic nature of designing complex systems in motion.
1. Introduction
Wind turbine is a mechanical system utilized to convert wind to electrical energy through high-speed generator. The size of commercial wind turbines has increased dramatically in the last 25 years from approximately rated power of 50 [KW] and a rotor diameter of 10 to 15 [m] up to 8 [MW] with a rotor diameter of more than 140 [m] [1]. The mechanical power generated via turbine blades is characterized by high torque and low rotational speed of the rotor, between 5 and 22 [rpm]. Thus, the next step is necessarily the design of power conversion system or drive-train that produces standard electricity from the shaft power developed by the rotor. Most of wind turbines in operation field have the traditional Danish design; that is, three blades are rigidly fixed to the rotor, which is indirectly coupled with an electric generator through a gearbox. Figure 1 shows the nacelle arrangement of such turbines including a high-speed induction generator, typically operated around 1000 to 2000 [rpm], in a clearly modular arrangement.
The design optimization of such a turbine and its subsystems will make future products more attractive compared with fossil and nuclear power plants. Instead of using the torsional models of wind turbine to simulate and design the various components of wind turbines, it is proposed to implement the multibody system dynamics approach [2, 3], to develop a consistent model to describe correctly the dynamic behavior of wind turbine blades and of the drive-train as well.
Under the umbrella of the multibody system dynamics, efficient procedure was developed for mapping airfoils geometry of wind turbine blades into flexible models using absolute nodal coordinate formulation (ANCF) [4, 5]. This procedure is extended to involve the complete aerodynamical and structural behavior of blades with the nonuniform and twisted configuration [6]. The effectiveness of using ANCF in modeling large-size wind turbine blades is examined through dynamic/static simulation procedure, and a method for designing such blades within DFD process is introduced successfully for practical wind conditions [7].
Although the cost of the rotor blades is approximately 20% of the total cost of wind turbine, it has a working lifetime of twenty years. In contrast, the gearboxes are noticed to fail ordinarily within operational time interval of five years, and need substitution. This is a costly task, since the substitution of a gearboxes accounts for about 10% of the construction and installation cost of the wind turbine. That will negatively affect the estimated income from a wind turbine.
Particular design procedures with the auxiliary function of computer simulation on wind turbine gearbox design have been discussed in literature [8]. Instead of the torsional-based models of the wind turbine drive-train, which were studied in many literatures [9–12], the multibody model enables obtaining reaction forces as well as reaction moments between the meshing gears and between gears and ground-structure.
This paper proposes a design procedure for wind turbine drive-train (WTDT) based on multibody system approach. This scheme can be integrated with the previously introduced procedures for modeling and design of wind blades [7]. The current work, which considers only rigid bodies, presents a computational scheme for designing the wind turbine drive-train based on gear tooth stresses. The flexibility effects are ignored assuming that they have little effect in estimating the maximum forces proposed for design process. Static as well as dynamic loads are calculated based on the Lagrange multipliers, which are numerical parameters associated with the holonomic and nonholonomic constraints assigned for the drive-train model. The numerical manipulations of kinematic stabilization, mapping the generalized forces into Cartesian respective, and implementation of velocity-based constrains into quadratic terms are carried out, and the corresponding subroutines are constructed. Based on the dynamic model and the obtained forces, the design process of a standard gearbox with (at least) one planetary stage is demonstrated with trade-offs optimization in terms of gearing parameters.
This paper is organized as follows: Section 2 presents the general form of equations of motion for multibody system and describes the relative motion between the interconnected bodies, that is, whether holonomic or nonholonomic type of constraints. Section 3 illustrates the proposed procedure and introduces the various design modules. Furthermore, two subsections are included to present the method of calculating the generalized aerodynamic forces and also the necessary equations for stress analysis. In Sections 4 and 5, a case study of practical wind turbine is manipulated using the proposed procedure. Section 6 summarizes the concluding remarks and conclusion.
2. Multibody Dynamics
In contrast to the blade model, which is definitely defined as flexible body, the gearbox model is constructed by considering only rigid bodies for its components. It can be shown that rigid body can be considered as a special case of the floating frame of reference formulation [13, 14], which can be used in estimating the stress along certain elements after carrying out the simulation. In this case, the system of interconnected bodies, the overall equation of motion can be written in a matrix form aswhere is the mass matrix and are the generalized coordinates and expressed as , such that are the Cartesian coordinates of the body frame origin and is the vector of Euler parameters that represent the rotational coordinates [15]. is the vector of Lagrange multipliers and is the vector of generalized forces. The vector absorbs all terms that are quadratic in the velocity; it can be expressed as [2]The matrix is the kinematic constraints Jacobian matrix, where the constraint equations can be presented assuch that and are the holonomic and nonholonomic constraint equations, respectively. If the nonholonomic constraints are linear in the generalized velocities, they could be expressed as follows [16]:where the matrix and vector depend on the system coordinates and may depend on time as well. These types of constraints restrict the generalized velocities and accordingly restrict the generalized acceleration but do not, however, apply any constraint on the generalized coordinates. Consequently, the number of independent coordinates is greater than the number of independent velocity or acceleration values. It should be mentioned here that the nonholonomic constraint is best suited to describe the relative motion between meshing gears; see Figure 2.
The Jacobian of constrained system corresponding to both holonomic and nonholonomic types of constraints can be expressed aswhere is the Jacobian matrix of the holonomic constraints and is the Jacobian matrix of the nonholonomic constraints. The time derivative of the constraints function can be written asIt therefore can be concluded thatThus, (2) can be rewritten asIt should be pointed out that the nonholonomic constraint equations given in (4) must be nonintegrable and cannot be reduced to an integrable form like the ones of geometric constraints of other joints.
In this case, the case of multibody system involving holonomic and nonholonomic constraints, the equations of motion can be expressed aswhere the vector includes all external and elastic forces, are the generalized reaction forces due to holonomic constraints, and are the generalized forces due to the nonholonomic constraints. These forces can be written as follows:
2.1. Constraints Forces
In this section, a procedure for determining the reaction forces associated with the generalized coordinates of constrained bodies is presented. These reactions are mapped into the Cartesian coordinates in the local and global domain as well. In spatial analysis, joints that connect two adjacent bodies have forces and moments at the joint definition points defined in the Cartesian space. These constraint forces significantly affect the dynamic stresses because these forces are heavily influenced by the impulsive contact forces. In the next subsections, the vector of Lagrange multipliers is evaluated corresponding to specified holonomic and nonholonomic constraints, and, accordingly, the reaction forces equations are derived in the form suitable for design process.
2.1.1. Holonomic Constraints Forces
The generalized constraint forces acting on rigid body , as a result of certain joint with body , can be written in terms of the Lagrange multipliers associated with the joint asThe virtual work of these generalized constraint forces can then be expressed asIn order to calculate the joint reaction forces associated with the Cartesian coordinates, let us assume that and be, respectively, the vectors of actual joint reaction forces and moments acting on body as a result of the certain (revolute, prismatic, etc.) joint with body . The virtual work of these reaction forces and moments can be written aswhere (see Figure 2) is the position vector of the point of application (say point ) of the force . is the absolute rotation vector of the body local coordinate system. The vector can be expressed as and therefore the vectors and can be expressed in terms of the generalized coordinates of body aswhere is the local position vector of the point of application of the force , defined in the body coordinate system, and virtual rotation is defined as of the body local coordinate system. Substituting (14) into (13), one obtains
Comparing (12) with (15), one obtainsSubstituting the value of of (16) into (17), one can conclude thatEquations (16) and (18) represent the necessary forms to evaluate the reaction forces and reaction moments for certain kinematic constraints.
2.1.2. Nonholonomic Constraints Forces
The virtual work of these generalized (nonholonomic) constraint forces, , can then be written asIn order to map , into the Cartesian coordinates, let us assume that is the vector of forces acting on body ; the virtual work of these forces can be written aswhere is the position vector of the tangential point, at which the nonholonomic constraints are formulated. Recall that the virtual displacement can be presented asOne obtainsComparing (19) with (22) yieldsThe following equation can be obtained:Thus, the resulting moment vector in the global frame can be concluded asIn the local coordinate system, one can obtain the following equations:Equations (23) and (27) represent the resulting forces, while (26) and (27) represent the forms of the resulting moments due to nonholonomic constraints in global and local frames, respectively.
2.1.3. Position and Velocity Stabilization
After each integration step, some constraint drifting at the position and velocity levels may take place. Therefore after each integration step relying on acceleration, it is necessary to move positions and velocities back to their manifolds. This corrections process is called poststabilization procedure and can be implemented using Newton-Raphson method [17]. The poststabilization of rotational multibody system subjected to nonholonomic constraints has been presented recently in [18].
3. DFD Procedure
Dynamics for Design (DFD) is the integration of recent advances in system dynamics, including nonlinearities, vibration analysis, and multibody systems with current design methodologies. The goal of integrating these subjects in one procedure is to obtain an efficient design cycle in terms of improved system reliability and consistent behavior across operating environments.
WTDT components include the rotor blades, gearbox (one, two, or three stages), and electrical generator. The purpose of the gearbox is to transfer the low rotational speed of the rotor to relatively high speed required for the generator. Consequently, the design process starts by assigning the preliminary size of the blade structure and from the other side by calculating the pressure distribution along the blade [7] and thus obtains the induced lift and drag forces. Mapping these forces from the blade coordinate system (frame) into the rotor frame enables the designer to calculate the generalized forces applied along the blade structure. Introducing the size of the blade as well as the generalized forces into the flexible multibody dynamics module based on continuum model of blade enables calculation of the induced torque and the associated rotational speed of the rotor blades. As shown in Figure 3, based on the outputs of blade dynamics module, the designer can assign the preliminary configuration and size of the required gearbox. It should be mentioned that small wind turbines are equipped with parallel-shaft gear trains, which are commercially available from numerous manufacturers. On the other hand, in larger wind turbines, the planetary design (epicyclic configuration) is definitely dominating. For outputs of several megawatts, two or three planetary stages are used within one or two parallel stages [19]. This step helps in calculating the geometrical dimensions and inertia parameters of the gearbox as multibody system.
Based on the multibody system approach, described in Section 2, and according to the preliminary configuration selected for the gearbox, the designer can calculate the reaction forces (between gears and ground) and contact forces between the gears. This step will be presented clearly during practical case study, presented in Section 4.
Once the reactions and contact forces are calculated, these forces can be mapped into components in the normal, tangential, and axial directions. Accordingly, the normal as well as shear stresses can be evaluated. Checking the stresses of the gears is carried out according to the American Gear Manufacturers Association (AGMA) [20], which has the responsible authority for the dissemination of standards pertaining to the design and analysis of gearing [20]. Resistance to tooth breakage is normally dependent upon the bending stress occurring at the contact area of mating teeth. These stresses are the key start of the design process by which trade-offs optimization can be carried out between material selection and the related dimensions of designed components.
In this section, two subsections are introduced in order to complete the picture in its general form. The first subsection presents the necessary equations of estimating the generalized forces induced due to the aerodynamic loads. The other subsection presents the necessary equations to calculate and examine the stresses according to AGMA code.
3.1. Aerodynamic Loads
The aerodynamic loads resulting from the attacking wind stream consist of two main components. The first component is called drag force, which comes parallel to the relative wind velocity. The second component is called lift force, and its direction is perpendicular to relative wind velocity. Based on the numeric values of the angle of attack , as illustrated in Figure 4, Reynolds and Mach number, the blade coefficients can be identified. These coefficients are the pressure distribution across the blade, , lift coefficient, , and the drag coefficient, , and can be estimated for specific geometric shape of turbine blade [7, 21].
As shown in Figure 4, the relative wind velocity is always greater than the upstream velocity at the rotor plane by the vectorial sum of tangential velocity ; that is, . This tangential velocity is induced from the rotor’s angular speed and thus depends on the local distance along the -axis on the blade frame; that is, , and, consequently, there is a relative velocity for each cross section of the blade; that is, . The equations of lift and drag forces can be written in following forms [22]:where is the air density, is the projected area, and is the dynamic pressure. These forces are defined with respect to the blade frame; accordingly, they should be mapped into rotor frame, which can be carried out by the following transformation:Thus, the generalized forces can be estimated by accumulating these terms along the blade span as follows [23]:
These forces are introduced to the blade structure module in order to estimate the induced torque and the associated rotational speed of the rotor blades.
3.2. Stress Analysis
The driver gear transmits motion to driven gear at contact point at pressure angel with respect to the contact line. This force can be calculated by obtaining the nonholonomic Lagrange multipliers, which is determined using the rigid body simulation of the model. It can be analyzed into two components, tangential and normal (radial) components. According to AGMA standard, the tangential component is the source of bending and contact stresses acting on gear tooth and these stresses are used to establish design or check safety for realistic design.
Under the assumptions stated in [24] and according to the flowchart shown in Figure 3, the stress analysis can be carried out in the postprocessing for design stage. The bending and contact stresses are estimated according to AGMA equations [20, 25]. The estimated stresses are compared with allowable stresses that gears materials can withstand; this ratio is called the factor of safety, , and it must be greater than one; otherwise the dimensions of gears should be modified and/or other material can be selected. The basic equations and terminology of bending and contact stresses according to AGMA are stated in the following paragraph.
The bending stress can be estimated aswhere is tangential force component, is the overload factor, is the dynamic factor, is the size factor, is the load distribution factor, and is the rim thickness factor. Also, is facewidth, is module of the gear, and is geometry factor on bending. Equation (32) is used to calculate actual bending stress acting on gear teeth. This value of stress is compared with allowable bending stress. The allowable bending stress can be estimated bywhere is the allowable bending stress (from materials tables), is the bending stress cycle factor, is the reliability factor, and is the temperature factor.
Furthermore, the actual contact stress can be estimated as follows:where is the elastic coefficient factor, is the pitch diameter of gear, and is the geometry factor on contact. The other values are defined previously in the bending stress equation. The allowable contact stress can be estimated aswhere is the allowable contact stress (from gears materials tables), is the contact stress cycle factor, and is the hardness ratio factor. By using (32)–(35), the safety factors can be estimated for bending and contact stresses as follows:
4. Case Study of Practical Wind Turbine
A case study of 1.4 MW nominal power is introduced in this section in order to evaluate the proposed procedure for WTDT design. The turbine is selected with the specifications listed in Table 1.
The nominal power is attained at wind speed of 10 [m/s], by which cut-in/cut-out speeds can be proposed as 4 [m/s] and 20 [m/s], respectively. By considering a constant wind speed of 10 [m/s], using the blade dynamics module (see Figure 3), the steady state value of the exerted torque on the rotor-blade is estimated as 1 [MN·m] (see Figure 5), and the associated rotational speed is found to be 14.2 [rpm] (see Figure 6). The high oscillatory nature shown in Figures 5 and 6 results from the large flexibility of blades structure. It is presumable that this oscillation will be damped out when the rotor hub is attached to the gearbox; this is because of the damping effect between rotating components and the wall structure [26].
4.1. Preliminary Gearbox Design
Planetary gearing systems introduce higher ratios of power transmission than parallel axis gears and are able to provide high-speed ratios within a compact volume. However, the inaccessibility of their components and high loads on the shaft bearings are considered one of the most important difficulties of maintenance and operations. The standard gearbox configurations include one or two planetary stages followed by one or two parallel stages. As is common in the horizontal wind turbines, the ring gear would be connected to the bedplate, and the planet carrier would be connected to the rotor hub, while the sun gear would be connected to the output shaft, which is connected to the next stage.
In the case study under consideration, standard gearbox configuration is selected such that one planetary stage and two parallel stages are connected in series. The overall speed ratio should be in range of 85 to 100 in order to attain the nominal speed of the electric generator. The preliminary size and specifications can be selected as listed in Table 2. In the planetary stage, the center distance incremental factor is held at zero, and thus the geometric relation of is valid, where , , and are the radii of ring, sun, and planet gears, respectively.
4.2. Multibody System and Constraints
As shown in Table 2, the three stages can be constructed using thirteen bodies; the ring gear is held to be fixed and considered as the wall structure (ground). The planetary stage consists of nine bodies (see Figure 7); three pins are considered to connect the carrier body with planet gears. The multibody system of the planetary stage can be constructed with the arrangement listed in Table 3. The total number of generalized coordinates is 63; the dependency between them can be determined by defining the constraints function according to the type of joints.
The holonomic constraints can be considered as listed in Table 4; the total number of holonomic constraints is 58; therefore, the system has 5 degrees of freedom. It should be noted here that every revolute joint is accompanied with significant damping proportion, between the carrier, sun, and the ground and between the pins and planet gears.
The nonholonomic constraints result from the equal tangential velocity of the meshing gears at their contact points and therefore zero relative velocity between them. Thus, the following equations can be applied:where is the center distance and is the radius of the ring gear. By considering no profile shifting, that is, the center distance incremental factor is zero, and using the definition of gear module, that is, , where is the pitch diameter and is the number of teeth, thus (38) can be rewritten as follows:Two more constraints equations can be added for the other two planet gears, to relate the rotational speed of the carrier, planet, and sun gears, such thatThe total number of constraints, the holonomic and nonholonomic constraints, became 62; thus the system has only one degree of freedom. By applying the exerted torque (see Figure 5) on the carrier body along with the rotational axis (local -axis as shown in Figure 7), the dynamic simulation of the planetary stage system shows that the system angular velocities obey the required ratios. As shown in Figure 8, the angular velocity of the carrier body is 14.2 [rpm] and the angular velocity of the sun gear is 85.2 [rpm] with speed ratio of 6.
The simulation results of the multibody system represented by (1) include the states vector (generalized acceleration and velocities) plus a vector of 62 Lagrange multipliers. Some of these multipliers are plotted in Figures 9 and 10. The number of Lagrange multipliers associated with holonomic constraints is 58 multipliers; among them, there are 9 trivial multipliers. These trivial multipliers represent the Euler parameters constraints , for each body. The other 49 multipliers can be used to calculate the reaction forces and moments by using (16) and (18).
Figure 9 shows the nonzero Lagrange multipliers associated with the first joint, that is, the rigid joint between the ring gear and the ground. The first and second Lagrange multipliers are zero, because they are related to the forces along - and -axes; however, the third multiplier is related to the force along the -axis, which is nonzero value at all cases because of the gear weight. The other multipliers are related to the reaction moments and associated with two Euler parameters and . The reaction forces and moments of the ring gear are plotted in Figures 11 and 12; note that the values are considered with respect to the global coordinate system (see Figure 7).
Figure 10 shows the Lagrange multipliers associated with the nonholonomic constraints, that is, (39)-(40). Since these multipliers are related to the angular velocities of the rotating gears and the ring, the resulting forces estimated by (23) are zeros; however, (26) and (28) give the driving torques of the rotating gears. Thus, the driving forces between gears can be obtained for the planet, sun, and also the ring gears. The tangential force between the sun and planet gear is plotted in Figure 13, with a steady state value of 412 [KN]. It should be noted that arbitrary values are selected for the pressure angle, module, center distance, and the facewidth of the matching gears. These values are taken as follows: , [mm], [mm], and [mm]. These numerical values are used to estimate the inertia properties of Table 3 by aid of CAD software. For the preliminary design and according to the tangential force estimated, it is found that the safety factor on contact between the sun and planet gears is and between the planet and ring gear is Also, the minimum safety factor on bending at the gear root is 2.1; that is, = 2.1, and at the gear flank it is 1.66; that is, = 1.66. These factors can be accepted as working design for the current gearing configuration.
5. Design Procedure
The design procedure assumes zero profile shifting coefficients for the planetary stage. These coefficients can be selected in order to obey the geometrical boundaries as well as to achieve some optimum values for the sliding velocities and safety factors. As mentioned earlier, the multibody model constructed in this paper concerns only the case of perfect center distances and, therefore, these coefficients are excluded from the optimization process. Based on the acceptable ranges of gearing ratios, center distances, and the gearing normal modules, 39 workable designs are found and listed in Table 5.
Figure 15 shows the workable designs which satisfy a center distance between 420 [mm] and 600 [mm]. Figure 14 shows the workable designs that satisfy the required range of speed ratios. As stated in Section 4.1, the speed ratio of the gearbox is accepted for values between 85 and 100; therefore, the speed ratio of the planetary stage should be accepted in the range of 6 to 6.7. By the use of the multibody model (see Figure 3), the tangential force is calculated for all workable designs listed in Table 5. The facewidth is reestimated according to AGMA; that is, = (3~5) [20]. The inertia properties are updated and then the contact and bending stresses and the corresponding safety factors are calculated. The accepted workable designs based on the safety factor on contact between the sun and planet gear and between the planet and ring gears are plotted in Figures 16 and 17. Similarly, the accepted workable designs based on the minimum safety factor on bending at the gears roots and flanks are illustrated in Figures 18 and 19. The center distances, speed ratios, and safety factors are plotted against the normal module, which comes in values between and 15 [mm].
The trade-offs optimization is a logical operation (AND/OR/NOT) among all workable designs. This operation is based on the required center distance , speed ratio , and acceptable safety factors . Practically, it can be carried out using the results shown in Figures 14–19. The optimal design, , is computed according to the following formula:
For instance, if the required center distance is 600 [mm], speed ratio is , and normal module is 10 [mm], (41) gives . The design number is the only design that satisfies all criteria, and the safety factors are greater than (see Table 5).
If minimum center distance is required, with acceptable safety factors , whatever the values of the module and speed ratio, the solution is . This solution gives a center distance of 420 [mm] with a speed ratio of , which lies in the required range, and module of 10 [mm]. In Table 5, design number 31 refers to an optimal design for compact size, precise speed ratio, and acceptable factors of safety.
The previous systematic procedure can be applied for designing the other parallel axis stages with less computational effort than the planetary stage.
6. Summary and Conclusions
Dynamics for Design (DFD) is the integration of recent advances in system dynamics including nonlinearities, vibration analysis, and multibody systems with current design methodologies. This paper proposed a design procedure for the wind turbine drive-train (WTDT). The design procedure is based on the multibody system (MBS) approach, which is used to estimate static as well as dynamic loads based on the Lagrange multipliers. The paper presented the necessary derivation of Lagrange multipliers associated with holonomic and nonholonomic constraints assigned for the drive-train model. Based on the dynamic model and the obtained forces, the design process of a planetary gear stage is implemented with trade-offs optimization in terms of gearing parameters. A wind turbine of 1.4 Megawatts nominal power is introduced as an evaluation study of the proposed procedure, in which the main advantage is the systematic nature of designing such complex system in motion. Finally, the nonholonomic constraints of the engaged gears presented in this paper assume zero profile shifting coefficients between them. This issue will be considered in future work to obey the allowable stresses and the geometrical boundaries of the drive-train as well.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.