Abstract
Dynamic load of helicopter rotors due to gust directly affects the structural stress and flight performance for helicopters. Based on a large deflection beam theory, an aeroelastic model for isolated helicopter rotors in the time domain is constructed. The dynamic response and structural load for a rotor under the impulse gust and slope-shape gust are calculated, respectively. First, a nonlinear Euler beam model with 36 degrees-of-freedoms per element is applied to depict the structural dynamics for an isolated rotor. The generalized dynamic wake model and Leishman-Beddoes dynamic stall model are applied to calculate the nonlinear unsteady aerodynamic forces on rotors. Then, we transformed the differential aeroelastic governing equation to an algebraic one. Hence, the widely used Newton-Raphson iteration algorithm is employed to simulate the dynamic gust load. An isolated helicopter rotor with four blades is studied to validate the structural model and the aeroelastic model. The modal frequencies based on the Euler beam model agree well with published ones by CAMRAD. The flap deflection due to impulse gust with the speed of 2m/s increases twice to the one without gust. In this numerical example, results indicate that the bending moment at the blade root is alleviated due to elastic effect.
1. Introduction
The response and dynamic load of a rotorcraft to discrete gusts are quite different from those of an aircraft. Because of the periodic inertial load and unsteady aerodynamic load induced by the rotation of rotor, the instantaneous dynamic stress undergoing gust would be significantly large. This kind of periodic dynamic stress leads to severe vibratory hub loads, which limits the operational envelope of helicopter rotors. Hence, the dynamic load of rotors to gust is considered as one of severe load cases in the design of helicopter and in the certification processes.
The dynamic load of hinged rotors by step gust and sinusoidal gust is analyzed by Arcidiacono et al. [1], with the steady inflow theory. They pointed out the requirement of gust load alleviation factor in MIL-S-8698 is a little conservative. Yasue et al. calculated the gust response by a harmonic balance method [2]. The sectional quasi-steady aerodynamic model was used to calculate aerodynamic forces. In addition, a wind tunnel test for isolated rotors was carried out to validate the model. Azuma and Saito applied the Momentum Theory to calculate the gust response with the consideration of flap-torsion motions [3], which is validated by a wind tunnel test for isolated rotors. Yin and Xiang investigated the aeroelastic response and hub loads of a composite hingeless rotor in a forward flight [4], with a 21 degrees-of-freedom (DOFs) beam theory and the quasi-steady aerodynamic model. The 21 DOFs beam theory is the generally used structural dynamic model for helicopter rotors.
The air flow through rotors is unsteady and complicated, especially near the blade root and the tip in the retreating side of the rotor disk. In addition, nonlinear aerodynamic force and periodic inertial load may induce large deflections of blades. Hence, the steady as well as quasi-steady aerodynamic model is not accurate enough to be applied to dynamic gust load analysis. Referring to the unsteady aerodynamic model, Bir and Chopra constructed a unified aeroelastic coupling model for rotors, with a famous medium deflection beam theory [5], to calculate the gust response of helicopters in hover and forward flights. Yeo and Johnson investigated the structural vibratory loads on several rotors by the CAMRAD II code [6]. Wang and Gao applied a small deflection beam theory and a medium deflection beam theory, to investigate the aeroelastic stability and structural load for helicopters [7]. A three DOFs rigid model coupled with a beam model of a small rotation angle was constructed to predict the structural load for rotors [8]. With the CFD-CSD loosely coupling algorithm, Park et al. developed a structural load analysis method for a rotor in descending flight [9], which is accurate in the modeling of both structural dynamics and aerodynamics.
Generally speaking, the aspect ratio of a helicopter blade is large enough to employ the beam theory for most aeroelastic analysis, while for gust response, the deformation of elastic blade due to gust may be very large. Hence, a large deflection beam model is developed to investigate suddenly gust-induced response. In this current work, the generalized dynamic wake (GDW) model [10] is used to calculate the induced velocity in the rotor disk, and the Leishman-Beddoes (L-B) dynamic stall model [11, 12] is to calculate the aerodynamic force on the blades. For the structural model, a 36 DOFs Euler beam theory is constructed for helicopter rotors [13, 14]. After the algebraic equation of dynamic motion in the time domain is developed, it is solved by the widely used Newton-Raphson iteration algorithm. Consequently, the dynamic load due to gust is investigated for elastic rotors.
2. Aeroelastic Model for Helicopter Rotors
The framework for dynamic load analysis for rotors is shown in Figure 1. The gust model, described by varied gust velocity, is added on the aerodynamic model. First, the unsteady and nonlinear aerodynamic model is constructed to calculate the aerodynamic force on the rotor disk. The dynamic pressure is induced by both rotor rotation and the gust velocity around the rotor disk. Then the structural dynamic model is developed to calculate the elastic displacement and velocity of a blade. The dynamic response and load can be calculated based on the coupling of structural and aerodynamic models in the time domain.
2.1. Structural Dynamic Model
A rotor is composed of several blades which are hinged to a hub. It is supposed that the rotor blades rotate around the hub with a fixed speed of . The global coordinate system is located at the hub, and its -axis is along one blade. It is also rotating along the -axis. Besides the global coordinate system, the local coordinated system is originated at the elastic center of each beam element, whose -axis is along the chord, named axis. -axis is perpendicular to -axis. It is also in the same cross section, named axis. The structural dynamics of each rotor blade is expressed as Euler beam elements [15, 16]. The displacement and rotation angles in the global coordinate system are used to describe the motion of each cross section. In addition, the nodal force and moment in the local coordinate system are also used as nodal variables. Hence, the nodal variables on each cross section are , where is the location of each node, is the velocity, is the rotation angle of each section, and is the angular velocity. and are the moments and forces on each cross section under the local coordinate system. On each section, the number of nodal variables is 18. Hence, there are total 36 nodal variables to describe the elastic motion for each beam element. The Euler beam element with nodal forces is shown in Figure 2.
For each beam element, the force-displacement relationship, the moment curvature relationship, the force-balance function, and displacement constraints between nodal variables are written by [15, 16]:where (1) and (2) are the constitution equations, relating the displacement ratio and curvature with the structural force and moment. Equation (3) is the force and moment-balance equations. Equation (4) is the constraints of displacement and velocity variables on each beam element. is a transformation matrix from the local beam-element coordinate system to the global coordinate system. In all equations, subscript denotes an average between two nodes of a beam element. 0 represents the initial value before deformation. is the curvature matrix. is the averaged stiffness matrix. and are the applied distributed force and moment, respectively. They are generated by the inertial and gravity loads. and are the applied concentrated force and moment on a beam node, respectively. In this current work, it is given by the aerodynamic load and joint load. Though the aerodynamic force is a distributed one, it is equivalently acting on the element node. is the unit-impulse function. Because the coefficients and are related with beam-element variables , the above equations are nonlinear differential equations. In addition, is not supposed to be small in the structural dynamic model. Therefore, this Euler beam model can be applied to rotation blades with small to large deflections.
Usually a helicopter rotor is composed of three or four blades. They are connected together by a hinge to the hub. Therefore, their relationships are described by a joint element. Under the global coordinate system, the rigid joint is supposed to connect two beam elements while keeping its length and orientation. Each joint element has 12 DOFs. It is denoted as . The variables of a joint are described by 12 equations. Those are [14]where stands for the tensor product operator. The subscripts 1 and 2 denote the left and right stations of a joint. and represent the applied moments and forces on the joint. Equations (5) and (6) are the displacement constraints for each joint. They mean the length and orientation of a joint element remain their initial values. Equations (7) and (8) are the force and moment-balance relationships. Since the equations of joint and beam element are expressed in the rotational global coordinate system, an additional rotation motion should be defined in an inertial coordinate system. Therefore, additional 18 global variables are introduced to describe the global rigid motion of rotors. More details of these equations are seen in [15].
Until now, the equations of beam element, joint element, and the global motion constitute the whole structural dynamic models in the time domain. The global boundary conditions here are set on the angular velocity and the angular acceleration. That is, , and . It means the rotor blades rotate at a constant angular velocity.
2.2. Aerodynamic Model
A suitable aerodynamic model should be constructed to depict the two characteristics of airflow through blades. One is unsteadiness due to blade rotating and deformation; the other is the nonlinearity of forces with angles of attack. Driving by these, two basic models are included in the aerodynamic model. One is the wake flow model; the other is the stall flow aerodynamic model. The generalized dynamic wake (GDW) flow model is used to calculate the induced velocity through the disk, considering the lag of flow wake [10]. The induced flow velocity varies with azimuth angles and blade radial locations. In addition, after calculating the induced velocity and gust velocity, the L-B Semi-Empirical model for dynamic stall is used to calculate the aerodynamic force on each blade section [11].
A basic issue to calculate the aerodynamic force is to calculate the local velocity through blades. Neglecting the spanwidth flow along the blade, the local velocity on each airfoil section depends on the summation of gust velocity, the induced wake velocity, and the velocity due to structural motion. That is given byIn this equation, the gust velocity is given by the gust model, which is mentioned in the following section. The induced velocity is calculated by the GDW method. The velocity of airfoil’s motion is calculated by the rotation speed and the structural dynamic model which is a coupling variable between the aerodynamic model and the structural dynamic model.
Afterwards, the aerodynamic forces are calculated by the L-B model. In this aerodynamic model, the aerodynamic forces are separated as three parts. Those are the forces due to separated flow in the trailing-edge, the vortex lift, and the noncirculatory lift part. They are calculated by each aerodynamic model, respectively. Consequently, the aerodynamic force and moment can be calculated by the summation of these three aerodynamic parts. The normal force , chord-wise force , and moment across the aerodynamic center are written aswhere is the chord length, is the above local flow velocity, and , , and are the normal force coefficient, chord-wise force coefficient, and the pitch moment coefficient, respectively. These aerodynamic forces can be transformed to the global coordinated system. Then the concentrated force on the structural dynamic model, indicated in (3), can be obtained. That is,
This generalized dynamic wake model and L-B stall aerodynamic model are available in AeroDyn [12] which is a turbine aerodynamic solver code. In each time step, the concentrated aerodynamic forces are the inputs of structural dynamic model. Vice versa, the blade motion velocity from the structural dynamic model is used to calculate the aerodynamic inflow velocity and forces. This coupling framework has been illustrated in Figure 1.
2.3. Gust Model
Discrete gust model is applied here. The gust field is assumed to be invariant with space. It indicates that the gust velocities at the whole rotor disk are the same. Two typical discrete gust models are considered. One is the impulse gust model; the other is the slope-shape gust model.(1)Impulse shape gust model: the form of impulse gust is shown in Figure 3(a). Here the vertical gust is chosen to investigate the structural dynamic load response. The strength of gust is with duration time of .(2)Slope-shape gust model: the form of slope-shape gust is shown in Figure 3(b). The time of slope length is with the maximum strength of .
(a) Impulse gust model
(b) Slope-shape gust model
3. Algorithm for Gust Load Calculation
The structural dynamic model is coupled with the aerodynamic model and the gust model in the time domain. Since the structural dynamic and the aerodynamic models are nonlinear ones, the coupling aeroelastic model is also expressed as a nonlinear differential equation in the time domain. Here, we introduce an extra system variable to represent . Then the differential equation is transformed to an algebraic one. It is written as
This equation is solved in two stages in sequence. In the first stage, the algebraic equation is solved by the Newton-Raphson iteration at each time step. For example, in the iteration loop at time step , the system variable in the equation is denoted as . Equation (12) can be written as where is the Jacobi matrix of to variable . is the Jacobi matrix of to the derivative of variables . In order to calculate the increment of in (13), we should deduce the expression of . Hence, a forward second order time differentiation scheme in the second stage is used to calculate the derivative of in the time domain, which yieldswhere , and are the differentiation coefficients. They are written asFrom (14), in the first initial calculations, we can also calculate that . Then (13) can be updated and solved. Afterwards, is updated until the change is less than tolerance:After the converged solution at time step is obtained.
Since the system variable includes structural velocity, the solution of this coupled aeroelastic model will approach to a quasi-steady state at each time step . This is a significant difference from the most Euler beam model.
The flowchart of dynamic load analysis for helicopter rotors due to gust is illustrated in Figure 4.
4. Model Validation
First, the modal analysis algorithm of the developed large-deflection beam theory is validated by a Mach Scale rotor in [17]. Afterwards, the aeroelastic response analysis algorithm is validated by another composite rotor in [18]. In order to validate the two algorithms based on the large deflection beam theory, both the results of the well-known medium deflection beam theory [5] with 15 DOFs and the published results in papers are presented and compared in this section.
4.1. Validation of Structural Dynamic Model Based on the Large Deflection Beam Theory
The modal frequencies are calculated by the large deflection beam theory. Their modal frequencies are also calculated by the medium deflection beam theory and compared with the published results. The rotor used here is called “Mach Scale.” It is a cantilever hingeless rotor with identical fiber-reinforced composite root flexures, designed and tested by U.S. Army Aeroflight dynamics Directorate [17]. This rotor has four blades with the outer shape shown in Figure 5(a). The blade radius is 1.143 m. The pretwist angle of this rotor is zero with solidity of 0.096. The tested structural parameters are shown in Table 1. With the geometric and structural parameters, its beam element model is constructed, shown in Figure 5(b). There are six elements on each blade.
![(a) Rotor blade [17]](https://static.hindawi.com/articles/sv/volume-2016/5727028/figures/5727028.fig.005a.jpg)
(a) Rotor blade [17]
(b) Rotor beam model
In Table 1, and represent the bending stiffness in the flap direction and in the lag direction, respectively. is the torsional stiffness and is the axial stiffness. represents the location of mass center in the sectional coordinate system. is the location of elastic axis in the sectional coordinate system. From this table, it can be seen that equals zero at all the cross sections. It indicates that the origin of the cross-sectional coordinate system is set on the elastic axis.
The blades are attached to a hingeless rotor hub. When the blades rotate around the hub, the structural dynamics is certainly different from the nonrotating state, which names dynamic stiffening effect. The large deflection Euler beam element is used to construct the finite element model for the rotor. Then modal analysis is conducted for both the rotating and nonrotating state. The frequencies calculated by large-deflection beam theory and medium-deflection beam theory are shown in Table 2. The modal shapes for the first two elastic modes are shown in Figure 6. From this table, the modal frequencies by the large deflection Euler beam model agree well with published results and the frequencies by medium deflection beam theory. In detail, in the rotation condition, the frequencies of the medium deflection beam model are almost the same as the published ones. From another aspect, in the nonrotation condition, the frequencies by the large deflection beam model agree quite well the experimental ones. Moreover, in the rotation condition, the errors between the modal frequencies by the large deflection beam theory and the published ones are less than 5%. It indicates that the structural dynamic model based on the large deflection beam theory is accurate enough. Though the modal analysis is not applied to the gust load analysis, the nonlinear structural dynamic model indicated in Section 2.1 is used to construct the aeroelastic model in the time domain.
For this “Mach Scale” rotor, the modal frequencies increase much with the rotational frequency. The natural frequencies with different rotational speeds are shown in Figure 7. It is observed that the flap frequency increases quickly with rotational speed while the lag frequency increases little.
4.2. Validation of Aeroelastic Response Algorithm by the Large Deflection Beam Theory
We programmed the aeroelastic response by both the large deflection beam model and the medium deflection beam model. However, we do not have the detailed and appropriate numerical examples to validate the aeroelastic algorithm by the large deflection beam theory. Hence, in order to validate it, the aeroelastic response by the medium deflection beam theory is first compared with published ones. Then the results by the large deflection beam theory (denoted as Model 1) are validated through the medium deflection beam theory (denoted as Model 2).
The parameters of this rotor are found in [18]. It is a composite rotor with four blades. The modal frequencies and aeroelastic behavior are calculated by the UMARC program, which is a widely accepted program in the helicopter aeroelastic analysis [18].
For this model, the cyclic control angle is . And the advanced ratio is 0.2. It is calculated by Model 2, the medium deflection beam theory. Figure 8 shows that the calculated blade tip deflections of Model 2 have a good agreement with the results of UMARC. It is indicated that Model 2 has a good accuracy to depict the aeroelastic behavior of a rotor.
Afterwards, both Model 1 and Model 2 are employed to calculate the aeroelastic response of the “Mach Scale” rotor in Section 4.1. In the compared calculation, the composite blades are rotating at a speed of 500 rpm. The forward flight speed is 30 m/s. The time step is 0.001 s in the Newton-Raphson algorithm of Model 1, and the tolerance is 10−5 at each time step. The calculation condition of Model 2 was set the same as the one of Model 1.
Figure 9 shows that the flap deformations at the blade tip of Model 1 and Model 2 are nearly the same. Since Model 1 applies the large deflection beam theory in the structure model, the deformation is a little bit larger than Model 2. The compared results indicate that Model 1 with the large deflection beam theory has a good accuracy to calculate the aeroelastic response. It is then applied to gust response analysis in the following section.
5. Gust Load Analysis Results
After the aeroelastic response of Model 1 is validated, it is applied for dynamic gust load calculation, coupled with the two gust models in Section 2.3. The rotor used here is also the “Mach Scale” helicopter model, shown in Section 4.1. The blades rotate at the speed of 500 rpm. The forward flight speed is 30 m/s. First, an impulse gust is acted on the rotor with a vertical velocity of 2 m/s. The starting time of impulse gust is s. The duration time of gust is s. The time step is 0.001 s in the Newton-Raphson algorithm, and the tolerance is 10−5 at each time step. We calculated the dynamic load according to the blade deflections and stiffness on each cross section.
Figure 10 shows the flap displacements at the blade tip due to impulse gust and slope-shape gust, respectively. From Figure 10(a), the flap displacement increases nearly 100% when impulse gust acts on it. Afterwards, when the gust disappears, the flap displacement recovers to a quasi-steady state. In Figure 10(b), induced by the slope-shape gust, the displacement magnitude grows with the gust amplitude and then it keeps a periodic motion with larger deformations. A direct reason is that the local blade speed and forces increase with gust. It is observed that the slope of flap displacement is similar as the gust shape shown in Figure 3. The maximum flap displacement reaches almost 2% of its span length. This deformation is not large. By the Fourier transformation, the frequency response of flap displacement due to impulse gust is indicated in Figure 11. From this figure, it is obviously seen that both the peak frequencies of lag displacement and flap displacement are 8.49 Hz. It is the first harmonic frequency of the rotational speed. The magnitude at the second harmonic frequency is much smaller. The third harmonic peak vanishes in this case. From the frequency response of the lag displacement, we can also observe a peak magnitude at the frequency of 14.99 Hz, which locates at the elastic lag mode. Hence, the instantaneous gust load is not only affected by the rotational motion but also affected by the elastic motions.
(a) Response due to impulse gust
(b) Response due to slope-shape gust
(a) Lag displacement
(b) Flap displacement
Figure 12 shows the bending moment at the blade root, due to the impulse gust. Compared with Figure 10(a), the maximum of bending moment is also increased due to gust during the time of 0.5 s to 0.6 s. We also calculated the dynamic load of a rigid blade with only a flexure in the root. The blades are assumed to be rigid by setting the blade stiffness coefficient to be 1015. The bending moment of the elastic blades is compared with the one of rigid blades, shown in Figure 13. From the comparisons of both the maximum bending moment and the mean bending moment, it can be seen that the bending moment along the blade spanwidth location is smaller than the one of rigid case. From Figure 14, we can find the reason that the shear force at the root location is larger than the rigid blade. With elastic deformation, the shear forces on the blades redistribute. It is larger in the root and smaller in the tip. Hence, to some extent, the bending moment is alleviated by the elastic deformation, which would be a benefit for structural strength.
(a) Maximum of the bending moment
(b) Mean of the bending moment
6. Conclusions
A structural dynamic model and a consequent aeroelastic model for isolated helicopter rotors are constructed. The dynamic load induced by the impulse gust and slope-shape gust is investigated. A “Mach Scale” helicopter rotor is used to validate the structural dynamic model and the aeroelastic model. From this work, we have the following conclusions:(1)The structural dynamic model based on the Euler beam theory with 36 DOFs per element is accurate and validated for rotor dynamic analysis. With the aeroelastic response comparison of UMARC result, the aeroelastic algorithm based on the large deflection beam theory has a good accuracy.(2)The flap displacement at the blade tip increases nearly 100% when impulse gust acts on it. Hence, the instantaneous gust load is not only affected by the blade rotational motion but also affected by the flap and lag elastic motions, which should be paid attention to in the dynamic gust load prediction for helicopter rotors.(3)The dynamic gust load will be redistributed with suitable elastic properties and may be alleviated, compared with the load of rigid blades.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (no. 11302011) and the Research Fund for the Doctoral Program of Higher Education of China (no. 20131102120051).