# Gust Response Analysis for Helicopter Rotors in the Hover and Forward Flights

Authors: Linpeng Wang; Yuting Dai; Chao Yang
Journal: Shock and Vibration
 (2017)

---

## Abstract

Dynamic load due to gust for helicopter rotors directly affects the structural stress and flight performance. In case of gust, it may cause the loss of trust force or the increase of deflection for rotors. In current work, an effective coupled aeroelastic model based on a medium-deflection beam theory and a nonlinear unsteady aerodynamic model in the time domain were constructed. Three types of gust in vertical direction were added in the model. The dynamic response and structural load for helicopter rotors under three types of gust were calculated, respectively. Results indicated that when rotors suffer a gust in hover at downward direction, the thrust force on rotor disk would decrease significantly when the gust amplitude increases, which should be paid attention in the design. Among the three gust types with the same gust strength, the maximum instantaneous shear force due to impulse shape gust is the largest. When the rotors suffer a gust in a forward flight, the shear force at the root of rotors would increase with the gust strength first but then it decreases. More attention should be paid to the decrease of thrust force and the increase of structural load in a forward flight.


---

## Body

## 1. Introduction

The response and dynamic load of helicopters to discrete gust are quite different from aircrafts, especially when unsteady aerodynamic force and elastic deformation of rotors make a coupling effect. The response of rotor under unsteady aerodynamic model is also different form steady or quasi-steady condition especially when rotor is under a large angle of attack where dynamic stall occurs. Moreover, the nonlinear aerodynamics couple with deformations of flap and twist of rotor blades. Because of the periodic inertial load and nonlinear aerodynamic load induced by rotation of rotor blades, the instantaneous dynamic stress of rotor blades undergoing gust would be significantly large [1]. This kind of periodic dynamic stress may degrade the performance of rotors which will affect design of stability and reliability for helicopter. In case of gust, it may cause the loss of trust force or the increase of deflection for rotors, which leads to the control instability or the blade failure. Hence, the dynamic load of rotors due to gust is considered as one of severe load cases in the design and certification processes of helicopter [2]. A helicopter is subjected to a wide variety of flight conditions with disparate objectives; these include efficient cruise and hover and high-performance maneuvering. The hover condition is a very important design consideration which is a limiting design point in terms of power requirements [3]. Accurate numerical prediction of rotor-blade aerodynamic parameters such as thrust requires accurate modeling of gust response [4, 5].Drees and Harvey [6] use single-rotor model to analyze the response due to different shapes of gust. They pointed out that the influences of gradual gust and gust shape on the response are quite strong. Arcidiacono et al. [7] applied the steady inflow theory to analyze the response due to step gust and sinusoidal gust on hinged rotors. They pointed out that the requirement of gust load alleviation factor in MIL-S-8698 is a little conservative. Elliott and Chopra [8] use finite element method and unsteady aerodynamic model based on a method using Fourier transform form Theodorsen’s lift to analyze response due to random gust.Yin and Xiang used the beam theory and quasi-steady aerodynamic model to investigate the aeroelastic response and hub load of a composite hingeless rotor in a forward flight [9]. Yasue et al. calculated the gust response by harmonic balance method [10]. To calculate aerodynamic forces, a sectional quasi-steady aerodynamic model was built and validated by a wind tunnel test for isolated rotors. Azuma and Saito applied the momentum theory to calculate the gust response with the consideration of flap bending motion [11], which is validated by a wind tunnel test for isolated rotors. This is the generally used structural dynamic model for helicopter rotors. Dai et al. [12] use large deformation beam theory to analyze the response of rotors under gust. It can be observed by comparison that medium-deflection beam theory is more suitable for gust analysis and has good accuracy to investigate the gust response for rotors.The rotor usually operates in the dynamic stall flow. The flow through rotor is complicated, especially near the blade root and the tip. In addition, nonlinear aerodynamic force and periodic inertial load may induce large deflections of blades. Hence, the steady and quasi-steady aerodynamic models are not accurate enough to be applied to dynamic gust load analysis [12]. Yeo and Johnson investigated the structural loads on several rotors by the CAMRAD II code [13]. Yeo and Johnson use CFD-CSD loosely coupling algorithm to develop a structural load analysis method for a rotor in a descending flight [13]. Referring to the unsteady aerodynamic model, Bir and Chopra constructed a unified rigid-elastic coupling model for rotor, with a famous medium-deflection beam theory [14], to calculate the gust response of helicopters in hover and forward flights. For reliability analysis, Maiti and Ammina calculated shear force for composite to analyze response due to gust [15].Generally speaking, for helicopter design, making an accurate prediction of dynamic load due to gust can reduce design redundancy and can give reference to the design of stability and reliability. The deformation of rotors due to the gust may change the load of rotor in hover and forward flight. The rigid model cannot take this effect into consideration. So a mathematical model of elastic rotors undergoing gust is constructed. The medium beam theory, which is based on the Hamilton principle was constructed for helicopter rotors. In this current work, the Pitt and Peters dynamic inflow model and Leishman Beddoes (L-B) dynamic stall model [16, 17] were used to calculate the induced velocity in the rotor disk and to calculate the aerodynamic force on the blades. After the dynamical differential equation of motion in the time domain was developed, it was solved by the widely used Newmark iteration algorithm. Finally, three shapes of vertical and multiscaled gust were investigated for elastic rotors at hover and forward flight.


## 2. Aeroelastic Model for Rotors

The framework for dynamic load analysis for rotors is shown in Figure1. The outputs of gust model and dynamic inflow model were described by varying the flow velocity components. They are added on the angle of attack calculation of aerodynamic model. First, the unsteady nonlinear aerodynamic model was constructed to calculate the aerodynamic force induced by blade motion and gust. Then the structural dynamic model was developed to calculate the displacement of rotor under the aerodynamic force and inertial force. The dynamic response and load can be calculated based on the coupling of structural and aerodynamic models in the time domain. Here, the blade deformation, such as flap deformation and twist angle, was regarded as a feedback to the aerodynamic model to calculate the local angle of attack in the next time step.Figure 1
Coupling strategy.


### 2.1. Aerodynamic Model

#### 2.1.1. Dynamic Inflow Model

In the forward flight, the Pitt and Peters model [18] was used as the dynamic inflow model to calculate the blade induced velocity. In the Pitt and Peters model, the distribution of induced velocity is not uniform on a rotor disk. The induced velocity can be expressed as a composition of several harmonics. For simplicity, the first order of harmonics was taken into consideration. So the nondimensional induced velocity can be expressed as follows:(1)v-=v-0+v-crRcos⁡ψ+v-srRsin⁡ψ,where v-0, v-c, and v-s are the nondimensional perturbed induced velocities and ψ is the azimuth angle.A distinguishing characteristic for unsteady aerodynamics is that the change in induced velocity lags behind the change in aerodynamic force and moment. The relationship of induced velocity with aerodynamic force and moment can be expressed as(2)ν˙-0ν˙-sν˙-c+Mm-1Lν-0ν-sν-c=Mm-1CTCLCM,where CT is the thrust coefficient; CL and CM are roll and pitch moment coefficients; Mm is the apparent mass matrix, and L is the inflow gain matrix. Their expressions are written in [18].The rotor shaft angle is defined asαs, advanced speed as Vf, rotation speed as Ω, and rotor radius as R. So the advanced ratio μ, the inflow ratio λ0, the steady included flow velocity v-1, wake skew angle at the disk αd, and mass flow parameter Vm can be expressed as follows [18]: (3)μ=Vfcos⁡αsΩR;λ0=Vfsin⁡αsΩR;ν-1=0.5CT;Vm=λ0+ν-1λ0+2ν-1+μ2λ0+ν-12+μ2;αd=arctan⁡λ0+ν-1μ.


#### 2.1.2. Dynamic Stall Model

The Leishman-Beddoes [19–21] dynamic stall model is a semiempirical model to calculate the unsteady aerodynamic force on an airfoil in the time domain. The flow characteristic in this model was divided into two stages: one is the attached flow and the other is the detached separation flow. They are described in detail in the following sections.(A) In the Attached Flow Stage. The circulatory normal force coefficient can be solved by the following equation:(4)CNCt=2π2VcβA1b1x1+A2b2x2,where β=1-Ma2 is the Prandtl-Glauert compressibility factor, Ma is the Mach number, V is the free-stream velocity, and c is the chord length. A1, A2, b1, and b2 are the parameters dependent on the airfoil profile. x1 and x2 can be solved by a state-space equation as follows:(5)x˙1x˙2=2Vββ2-b100-b2x1x2+11α3/4t,where(6)α3/4t=αt+qt2.α
3
/
4
(
t
) is the angle of attack at the 3/4 chord location, α(t) is angle of attack, and q(t) is a nondimensional pitch rate. That is, q=α˙c/V. The noncirculatory normal force coefficient can be calculated by(7)CNαIt=4Max˙3,where(8)x˙3=αt-1KαTIx3.T
I is the noncirculatory time constant. Its expression is TI=c/a, a is the sonic velocity, and Kα is expressed as(9)Kα=1-Ma+πβMa2A1b1+A2b2-1.The noncirculatory normal force coefficient caused by the pitch rate is written as(10)CNqIt=1Max˙4,where(11)x˙4=qt-1KqTIx4,KqMa=1-Ma+2πβMa2A1b1+A2b2-1.Hence, the total normal force coefficient can be got by the summation of them. That is,(12)CNt=CNCt+CNαIt+CNqIt.(B) In the Separation Stage. Here we define an effective angle of attack αf(t) to determine whether flow separation occurs on the airfoil:(13)αft=CN′tCNαMa,where CN′(t)=x9 and x9 can be solved via(14)x˙9=-x9TP+CNPtTP2Vc.HereTp is a time constant for leading edge pressure response.Thenαf is used to replace α in the following equation:(15)f=1-0.3exp⁡α-α1S1ifα≤α10.04+0.66exp⁡α-α1S1ifα>α1.The first-order lag in the separation point can be obtained as follows:(16)f′=1-0.3exp⁡αf-α1S1ifαf≤α10.04+0.66exp⁡αf-α1S1ifαf>α1,where S1 and S2 define the stall characteristics. While α1 represents the separate point corresponding to f=0.7, S1, S2 and α1 are determined empirically for each airfoil and they may vary with Mach number.The final separation pointf′′=x10 can be obtained via solving the following state-space equation:(17)x˙10=-x10Tf+f′Tf,where Tf is the time constant for separation point movement.When the flow separation occurs in the leading edge, the normal force coefficient is calculated by(18)CNf=CNαMa1+f′′2αEt,where the effective angle of attack αE(t) is defined as(19)αEt=β22VcA1b1x1+A2b2x2.In the stage of flow separation in the leading edge, the aerodynamic moment is(20)CMf=K0+K11-f′′+K2sin⁡πf′′mmCNC+CM0,where CM0 is the moment coefficient at zero angle of attack. The values of K0, K1, K2, andmm can be adjusted for different airfoils; typically for NACA0012,mm = 2.Here the critical normal force coefficientCN1 varies for different airfoils. The general case of dynamic stall involves the formation of a vortex near the leading edge of the airfoil which subsequently separates from the airfoil and it transports downstream.The critical condition used here isCN′(t)≥CN1. At this point, a catastrophic flow separation occurs and the accumulated circulation passes over the airfoil and it flows into the wake.The vortex lift can be represented by the following differential equations:(21)x˙11=-x11Tv+C˙vTv,CNvt=x11,where (22)Cv=CNC1-1+f′′24forτv≤2Tvl0forτv>2Tvl.Here(23)τv=τv+2×0.45VΔtc.T
v is the time constant for vortex lift. Tvl is the time constant for vortex traversing over the whole airfoil.The corresponding pitch momentCMv produced by the vortex lift component is given by(24)CMv=-CPvCNv,where(25)CPv=0.251-cos⁡πτvTvl.Then the final normal force coefficient can be calculated:(26)CNt=CNft+CNvt.And the aerodynamic moment coefficient is summated by(27)CMt=CMft+CMvt.


### 2.2. Structural Dynamic Model

The Hamilton principle is used to construct the structural dynamic equation for rotors [22–25]. The equation is expressed as(28)δΠ=∫t1t2δU-δT-δWdt=0,where U is the strain energy, T is the kinetic energy, and W is the virtual work of external forces.In the form of finite element method, the discrete equation is(29)δΠ=∫ψ1ψ2∑i=1NUMδUi-δTi-δWidψ=∫ψ1ψ2∑i=1NUMΔidψ=0,where NUM is the number of total beam elements for each blade.In each beam element, the axial deflectionu, the in-plane and out-plane deflections vss and w, and elastic torsion ϕ are taken into consideration:(30)usvssswsϕs=Hu0000H0000H0000Hϕqi,where H is the Hermite function; here,(31)HuT=Hu1Hu3Hu4Hu2=-4.5s3+9s2-5.5s+113.5s3-22.5s2+9s-13.5s3+18s2-4.5s4.5s3-4.5s2+sHT=H1H2H3H4=2s3-3s2+1lis3-2s2+s-2s3+3s2lis3-s2HϕT=Hϕ1Hϕ3Hϕ2=2s2-3s+1-4s2+4s2s2-s.s is independent variable; qi is defined as follows, which is shown in Figure 2: (32)qiT=u1u2u3u4vss1vss1′vss2vss2′w1w1′w2w2′ϕ1ϕ2ϕ3.Figure 2
Degree of freedom on element node.Then,Δi in the discrete form is deduced to the following form:(33)Δi=δqiTMq¨+Cq˙+Kq-Fi,M is the element mass matrix; C is the damping matrix; K is the stiffness matrix; F is the vector of generalized force.The strain and kinetic energy of each element are added together. Then the Hamilton equation is expressed as(34)δΠ=∫ψ1ψ2∑i=1NUMδqiTMq¨+Cq˙+Kq-Fidψ=0.


### 2.3. Gust Model

We assume that the gust field is invariant with space. It indicates that the gust velocities at the whole rotor disk are the same. According to the gust field, three typical vertical gust profiles are considered.(1) Impulse-Shape Gust. The form of impulse gust is shown in Figure 3(a). Here the vertical gust is chosen to investigate the structural dynamic response. The strength of gust is W with duration time oft.Figure 3
Gust models.
(a)
Impulse gust model
(b)
Slope-shape gust model
(c)
1
-
c
o
s shape gust model(2) Slope-Shape Gust. The form of slope-shape gust is shown in Figure 3(b). The time of slope length is t with the maximum strength ofW.(3)
1
-
c
o
s
Shape Gust. The form of 1-cos shape gust is shown in Figure 3 Gust models Figure 3(c). The maximum strength isW.


### 2.4. Aeroelastic Response Calculation Algorithm

To calculate the aeroelastic response or gust response, the aerodynamic model needs to be coupled with the structural one. In the dynamic inflow model,v-0, v-c, and v-s can be obtained by solving the state-space equation (2). Thus the total induced velocity can be obtained by the summation of v-0, v-c, and v-s It is pointed out that in the coupled model we recover the induced speed to a dimensional one. That is, v=v-ΩR.The geometric relationship of blade element is shown in Figure4. The inflow angle can be calculated via the following equation:(35)βin=arctan⁡ν-vhΩr+Vfcos⁡Ωt-2jπ/N.vh is the flap speed of blade, which can be calculated via solving the structural dynamic model. N is number of blades andj is the corresponding number of each blade. In this paper, the rotor has a total of four blades. They are marked as 1,2,3,4, respectively.Figure 4
Geometric relationship of blade element.Because the flap and twist are the main deformations in the aeroelastic response, the lag deflection is too small, so it is neglected.The angle of attackα is composed of (36)α=θ0-βin+ϕ,where ϕ is twist deformation of a blade which can be calculated in the structural dynamic model and θ0 is the control angle, which is expressed as(37)θ0=θ75+θtwxR-0.75+θ1ccos⁡ψ+θ1ssin⁡ψ.θ
75 is the collective pitch, θtw is the pretwist angle, θ1c is the transverse cycle control, and θ1s is the longitudinal cycle control.The Process of Calculation. All variables need to be initialized first. Then the induced velocity is calculated. Afterwards, the angle of attack is updated as the input of aerodynamic model. The normal force coefficient and moment coefficient CN and CM can be calculated by the L-B aerodynamic model.The normal force and moment are calculated as follows:(38)FN=0.5ρVl2CNcl,FC=0.5ρVl2CCcl,MM=0.5ρVl2CMc2l,where ρ is the air density, Vl is the local speed, and l is length of each element.The lift and dragon on each element are calculated as follows:(39)FL=FNcos⁡α-FCsin⁡α,FD=FNsin⁡α+FCcos⁡α,FZ=FLcos⁡βin-FDsin⁡βin.The total force of each rotor is calculated by adding all the forces of elements. Forces on the whole rotor disc can be calculated via the following equations:(40)FTT_all=∑0N∑0NUMFZ,FLL_all=∑0N∑0NUMFZrcos⁡Ωt-2jπN,FMM_all=∑0N∑0NUMFZrsin⁡Ωt-2jπN.The following forces are nondimensionalized to be regarded as input in the dynamic inflow model:(41)C_T=FTT_allρπR2ΩR2,C_L=FLL_allρπR3ΩR2,C_M=FMM_allρπR3ΩR2.Then the induced velocity over the whole rotor disk can be calculated again via the following equation:(42)vji=ν-jiΩR=ν-0+ν-criRcos⁡Ωt-2jπN+ν-sriRsin⁡Ωt-2jπNΩR.All variables are calculated again to start a new iteration.C language was used to code the model; the whole calculation flowchart is shown in Figure 5.Figure 5
Flowchart of coupled strategy.When the effect of gust is considered in the model, the effect is added in the calculation of inflow angle. In this case, the inflow angle is modified as follows:(43)βin=arctan⁡ν-vh+vgustΩr+Vfcos⁡Ωt-2jπ/N,where vgust is the speed of downward gust.


## 2.1. Aerodynamic Model

### 2.1.1. Dynamic Inflow Model

In the forward flight, the Pitt and Peters model [18] was used as the dynamic inflow model to calculate the blade induced velocity. In the Pitt and Peters model, the distribution of induced velocity is not uniform on a rotor disk. The induced velocity can be expressed as a composition of several harmonics. For simplicity, the first order of harmonics was taken into consideration. So the nondimensional induced velocity can be expressed as follows:(1)v-=v-0+v-crRcos⁡ψ+v-srRsin⁡ψ,where v-0, v-c, and v-s are the nondimensional perturbed induced velocities and ψ is the azimuth angle.A distinguishing characteristic for unsteady aerodynamics is that the change in induced velocity lags behind the change in aerodynamic force and moment. The relationship of induced velocity with aerodynamic force and moment can be expressed as(2)ν˙-0ν˙-sν˙-c+Mm-1Lν-0ν-sν-c=Mm-1CTCLCM,where CT is the thrust coefficient; CL and CM are roll and pitch moment coefficients; Mm is the apparent mass matrix, and L is the inflow gain matrix. Their expressions are written in [18].The rotor shaft angle is defined asαs, advanced speed as Vf, rotation speed as Ω, and rotor radius as R. So the advanced ratio μ, the inflow ratio λ0, the steady included flow velocity v-1, wake skew angle at the disk αd, and mass flow parameter Vm can be expressed as follows [18]: (3)μ=Vfcos⁡αsΩR;λ0=Vfsin⁡αsΩR;ν-1=0.5CT;Vm=λ0+ν-1λ0+2ν-1+μ2λ0+ν-12+μ2;αd=arctan⁡λ0+ν-1μ.


### 2.1.2. Dynamic Stall Model

The Leishman-Beddoes [19–21] dynamic stall model is a semiempirical model to calculate the unsteady aerodynamic force on an airfoil in the time domain. The flow characteristic in this model was divided into two stages: one is the attached flow and the other is the detached separation flow. They are described in detail in the following sections.(A) In the Attached Flow Stage. The circulatory normal force coefficient can be solved by the following equation:(4)CNCt=2π2VcβA1b1x1+A2b2x2,where β=1-Ma2 is the Prandtl-Glauert compressibility factor, Ma is the Mach number, V is the free-stream velocity, and c is the chord length. A1, A2, b1, and b2 are the parameters dependent on the airfoil profile. x1 and x2 can be solved by a state-space equation as follows:(5)x˙1x˙2=2Vββ2-b100-b2x1x2+11α3/4t,where(6)α3/4t=αt+qt2.α
3
/
4
(
t
) is the angle of attack at the 3/4 chord location, α(t) is angle of attack, and q(t) is a nondimensional pitch rate. That is, q=α˙c/V. The noncirculatory normal force coefficient can be calculated by(7)CNαIt=4Max˙3,where(8)x˙3=αt-1KαTIx3.T
I is the noncirculatory time constant. Its expression is TI=c/a, a is the sonic velocity, and Kα is expressed as(9)Kα=1-Ma+πβMa2A1b1+A2b2-1.The noncirculatory normal force coefficient caused by the pitch rate is written as(10)CNqIt=1Max˙4,where(11)x˙4=qt-1KqTIx4,KqMa=1-Ma+2πβMa2A1b1+A2b2-1.Hence, the total normal force coefficient can be got by the summation of them. That is,(12)CNt=CNCt+CNαIt+CNqIt.(B) In the Separation Stage. Here we define an effective angle of attack αf(t) to determine whether flow separation occurs on the airfoil:(13)αft=CN′tCNαMa,where CN′(t)=x9 and x9 can be solved via(14)x˙9=-x9TP+CNPtTP2Vc.HereTp is a time constant for leading edge pressure response.Thenαf is used to replace α in the following equation:(15)f=1-0.3exp⁡α-α1S1ifα≤α10.04+0.66exp⁡α-α1S1ifα>α1.The first-order lag in the separation point can be obtained as follows:(16)f′=1-0.3exp⁡αf-α1S1ifαf≤α10.04+0.66exp⁡αf-α1S1ifαf>α1,where S1 and S2 define the stall characteristics. While α1 represents the separate point corresponding to f=0.7, S1, S2 and α1 are determined empirically for each airfoil and they may vary with Mach number.The final separation pointf′′=x10 can be obtained via solving the following state-space equation:(17)x˙10=-x10Tf+f′Tf,where Tf is the time constant for separation point movement.When the flow separation occurs in the leading edge, the normal force coefficient is calculated by(18)CNf=CNαMa1+f′′2αEt,where the effective angle of attack αE(t) is defined as(19)αEt=β22VcA1b1x1+A2b2x2.In the stage of flow separation in the leading edge, the aerodynamic moment is(20)CMf=K0+K11-f′′+K2sin⁡πf′′mmCNC+CM0,where CM0 is the moment coefficient at zero angle of attack. The values of K0, K1, K2, andmm can be adjusted for different airfoils; typically for NACA0012,mm = 2.Here the critical normal force coefficientCN1 varies for different airfoils. The general case of dynamic stall involves the formation of a vortex near the leading edge of the airfoil which subsequently separates from the airfoil and it transports downstream.The critical condition used here isCN′(t)≥CN1. At this point, a catastrophic flow separation occurs and the accumulated circulation passes over the airfoil and it flows into the wake.The vortex lift can be represented by the following differential equations:(21)x˙11=-x11Tv+C˙vTv,CNvt=x11,where (22)Cv=CNC1-1+f′′24forτv≤2Tvl0forτv>2Tvl.Here(23)τv=τv+2×0.45VΔtc.T
v is the time constant for vortex lift. Tvl is the time constant for vortex traversing over the whole airfoil.The corresponding pitch momentCMv produced by the vortex lift component is given by(24)CMv=-CPvCNv,where(25)CPv=0.251-cos⁡πτvTvl.Then the final normal force coefficient can be calculated:(26)CNt=CNft+CNvt.And the aerodynamic moment coefficient is summated by(27)CMt=CMft+CMvt.


## 2.1.1. Dynamic Inflow Model

In the forward flight, the Pitt and Peters model [18] was used as the dynamic inflow model to calculate the blade induced velocity. In the Pitt and Peters model, the distribution of induced velocity is not uniform on a rotor disk. The induced velocity can be expressed as a composition of several harmonics. For simplicity, the first order of harmonics was taken into consideration. So the nondimensional induced velocity can be expressed as follows:(1)v-=v-0+v-crRcos⁡ψ+v-srRsin⁡ψ,where v-0, v-c, and v-s are the nondimensional perturbed induced velocities and ψ is the azimuth angle.A distinguishing characteristic for unsteady aerodynamics is that the change in induced velocity lags behind the change in aerodynamic force and moment. The relationship of induced velocity with aerodynamic force and moment can be expressed as(2)ν˙-0ν˙-sν˙-c+Mm-1Lν-0ν-sν-c=Mm-1CTCLCM,where CT is the thrust coefficient; CL and CM are roll and pitch moment coefficients; Mm is the apparent mass matrix, and L is the inflow gain matrix. Their expressions are written in [18].The rotor shaft angle is defined asαs, advanced speed as Vf, rotation speed as Ω, and rotor radius as R. So the advanced ratio μ, the inflow ratio λ0, the steady included flow velocity v-1, wake skew angle at the disk αd, and mass flow parameter Vm can be expressed as follows [18]: (3)μ=Vfcos⁡αsΩR;λ0=Vfsin⁡αsΩR;ν-1=0.5CT;Vm=λ0+ν-1λ0+2ν-1+μ2λ0+ν-12+μ2;αd=arctan⁡λ0+ν-1μ.


## 2.1.2. Dynamic Stall Model

The Leishman-Beddoes [19–21] dynamic stall model is a semiempirical model to calculate the unsteady aerodynamic force on an airfoil in the time domain. The flow characteristic in this model was divided into two stages: one is the attached flow and the other is the detached separation flow. They are described in detail in the following sections.(A) In the Attached Flow Stage. The circulatory normal force coefficient can be solved by the following equation:(4)CNCt=2π2VcβA1b1x1+A2b2x2,where β=1-Ma2 is the Prandtl-Glauert compressibility factor, Ma is the Mach number, V is the free-stream velocity, and c is the chord length. A1, A2, b1, and b2 are the parameters dependent on the airfoil profile. x1 and x2 can be solved by a state-space equation as follows:(5)x˙1x˙2=2Vββ2-b100-b2x1x2+11α3/4t,where(6)α3/4t=αt+qt2.α
3
/
4
(
t
) is the angle of attack at the 3/4 chord location, α(t) is angle of attack, and q(t) is a nondimensional pitch rate. That is, q=α˙c/V. The noncirculatory normal force coefficient can be calculated by(7)CNαIt=4Max˙3,where(8)x˙3=αt-1KαTIx3.T
I is the noncirculatory time constant. Its expression is TI=c/a, a is the sonic velocity, and Kα is expressed as(9)Kα=1-Ma+πβMa2A1b1+A2b2-1.The noncirculatory normal force coefficient caused by the pitch rate is written as(10)CNqIt=1Max˙4,where(11)x˙4=qt-1KqTIx4,KqMa=1-Ma+2πβMa2A1b1+A2b2-1.Hence, the total normal force coefficient can be got by the summation of them. That is,(12)CNt=CNCt+CNαIt+CNqIt.(B) In the Separation Stage. Here we define an effective angle of attack αf(t) to determine whether flow separation occurs on the airfoil:(13)αft=CN′tCNαMa,where CN′(t)=x9 and x9 can be solved via(14)x˙9=-x9TP+CNPtTP2Vc.HereTp is a time constant for leading edge pressure response.Thenαf is used to replace α in the following equation:(15)f=1-0.3exp⁡α-α1S1ifα≤α10.04+0.66exp⁡α-α1S1ifα>α1.The first-order lag in the separation point can be obtained as follows:(16)f′=1-0.3exp⁡αf-α1S1ifαf≤α10.04+0.66exp⁡αf-α1S1ifαf>α1,where S1 and S2 define the stall characteristics. While α1 represents the separate point corresponding to f=0.7, S1, S2 and α1 are determined empirically for each airfoil and they may vary with Mach number.The final separation pointf′′=x10 can be obtained via solving the following state-space equation:(17)x˙10=-x10Tf+f′Tf,where Tf is the time constant for separation point movement.When the flow separation occurs in the leading edge, the normal force coefficient is calculated by(18)CNf=CNαMa1+f′′2αEt,where the effective angle of attack αE(t) is defined as(19)αEt=β22VcA1b1x1+A2b2x2.In the stage of flow separation in the leading edge, the aerodynamic moment is(20)CMf=K0+K11-f′′+K2sin⁡πf′′mmCNC+CM0,where CM0 is the moment coefficient at zero angle of attack. The values of K0, K1, K2, andmm can be adjusted for different airfoils; typically for NACA0012,mm = 2.Here the critical normal force coefficientCN1 varies for different airfoils. The general case of dynamic stall involves the formation of a vortex near the leading edge of the airfoil which subsequently separates from the airfoil and it transports downstream.The critical condition used here isCN′(t)≥CN1. At this point, a catastrophic flow separation occurs and the accumulated circulation passes over the airfoil and it flows into the wake.The vortex lift can be represented by the following differential equations:(21)x˙11=-x11Tv+C˙vTv,CNvt=x11,where (22)Cv=CNC1-1+f′′24forτv≤2Tvl0forτv>2Tvl.Here(23)τv=τv+2×0.45VΔtc.T
v is the time constant for vortex lift. Tvl is the time constant for vortex traversing over the whole airfoil.The corresponding pitch momentCMv produced by the vortex lift component is given by(24)CMv=-CPvCNv,where(25)CPv=0.251-cos⁡πτvTvl.Then the final normal force coefficient can be calculated:(26)CNt=CNft+CNvt.And the aerodynamic moment coefficient is summated by(27)CMt=CMft+CMvt.


## 2.2. Structural Dynamic Model

The Hamilton principle is used to construct the structural dynamic equation for rotors [22–25]. The equation is expressed as(28)δΠ=∫t1t2δU-δT-δWdt=0,where U is the strain energy, T is the kinetic energy, and W is the virtual work of external forces.In the form of finite element method, the discrete equation is(29)δΠ=∫ψ1ψ2∑i=1NUMδUi-δTi-δWidψ=∫ψ1ψ2∑i=1NUMΔidψ=0,where NUM is the number of total beam elements for each blade.In each beam element, the axial deflectionu, the in-plane and out-plane deflections vss and w, and elastic torsion ϕ are taken into consideration:(30)usvssswsϕs=Hu0000H0000H0000Hϕqi,where H is the Hermite function; here,(31)HuT=Hu1Hu3Hu4Hu2=-4.5s3+9s2-5.5s+113.5s3-22.5s2+9s-13.5s3+18s2-4.5s4.5s3-4.5s2+sHT=H1H2H3H4=2s3-3s2+1lis3-2s2+s-2s3+3s2lis3-s2HϕT=Hϕ1Hϕ3Hϕ2=2s2-3s+1-4s2+4s2s2-s.s is independent variable; qi is defined as follows, which is shown in Figure 2: (32)qiT=u1u2u3u4vss1vss1′vss2vss2′w1w1′w2w2′ϕ1ϕ2ϕ3.Figure 2
Degree of freedom on element node.Then,Δi in the discrete form is deduced to the following form:(33)Δi=δqiTMq¨+Cq˙+Kq-Fi,M is the element mass matrix; C is the damping matrix; K is the stiffness matrix; F is the vector of generalized force.The strain and kinetic energy of each element are added together. Then the Hamilton equation is expressed as(34)δΠ=∫ψ1ψ2∑i=1NUMδqiTMq¨+Cq˙+Kq-Fidψ=0.


## 2.3. Gust Model

We assume that the gust field is invariant with space. It indicates that the gust velocities at the whole rotor disk are the same. According to the gust field, three typical vertical gust profiles are considered.(1) Impulse-Shape Gust. The form of impulse gust is shown in Figure 3(a). Here the vertical gust is chosen to investigate the structural dynamic response. The strength of gust is W with duration time oft.Figure 3
Gust models.
(a)
Impulse gust model
(b)
Slope-shape gust model
(c)
1
-
c
o
s shape gust model(2) Slope-Shape Gust. The form of slope-shape gust is shown in Figure 3(b). The time of slope length is t with the maximum strength ofW.(3)
1
-
c
o
s
Shape Gust. The form of 1-cos shape gust is shown in Figure 3 Gust models Figure 3(c). The maximum strength isW.


## 2.4. Aeroelastic Response Calculation Algorithm

To calculate the aeroelastic response or gust response, the aerodynamic model needs to be coupled with the structural one. In the dynamic inflow model,v-0, v-c, and v-s can be obtained by solving the state-space equation (2). Thus the total induced velocity can be obtained by the summation of v-0, v-c, and v-s It is pointed out that in the coupled model we recover the induced speed to a dimensional one. That is, v=v-ΩR.The geometric relationship of blade element is shown in Figure4. The inflow angle can be calculated via the following equation:(35)βin=arctan⁡ν-vhΩr+Vfcos⁡Ωt-2jπ/N.vh is the flap speed of blade, which can be calculated via solving the structural dynamic model. N is number of blades andj is the corresponding number of each blade. In this paper, the rotor has a total of four blades. They are marked as 1,2,3,4, respectively.Figure 4
Geometric relationship of blade element.Because the flap and twist are the main deformations in the aeroelastic response, the lag deflection is too small, so it is neglected.The angle of attackα is composed of (36)α=θ0-βin+ϕ,where ϕ is twist deformation of a blade which can be calculated in the structural dynamic model and θ0 is the control angle, which is expressed as(37)θ0=θ75+θtwxR-0.75+θ1ccos⁡ψ+θ1ssin⁡ψ.θ
75 is the collective pitch, θtw is the pretwist angle, θ1c is the transverse cycle control, and θ1s is the longitudinal cycle control.The Process of Calculation. All variables need to be initialized first. Then the induced velocity is calculated. Afterwards, the angle of attack is updated as the input of aerodynamic model. The normal force coefficient and moment coefficient CN and CM can be calculated by the L-B aerodynamic model.The normal force and moment are calculated as follows:(38)FN=0.5ρVl2CNcl,FC=0.5ρVl2CCcl,MM=0.5ρVl2CMc2l,where ρ is the air density, Vl is the local speed, and l is length of each element.The lift and dragon on each element are calculated as follows:(39)FL=FNcos⁡α-FCsin⁡α,FD=FNsin⁡α+FCcos⁡α,FZ=FLcos⁡βin-FDsin⁡βin.The total force of each rotor is calculated by adding all the forces of elements. Forces on the whole rotor disc can be calculated via the following equations:(40)FTT_all=∑0N∑0NUMFZ,FLL_all=∑0N∑0NUMFZrcos⁡Ωt-2jπN,FMM_all=∑0N∑0NUMFZrsin⁡Ωt-2jπN.The following forces are nondimensionalized to be regarded as input in the dynamic inflow model:(41)C_T=FTT_allρπR2ΩR2,C_L=FLL_allρπR3ΩR2,C_M=FMM_allρπR3ΩR2.Then the induced velocity over the whole rotor disk can be calculated again via the following equation:(42)vji=ν-jiΩR=ν-0+ν-criRcos⁡Ωt-2jπN+ν-sriRsin⁡Ωt-2jπNΩR.All variables are calculated again to start a new iteration.C language was used to code the model; the whole calculation flowchart is shown in Figure 5.Figure 5
Flowchart of coupled strategy.When the effect of gust is considered in the model, the effect is added in the calculation of inflow angle. In this case, the inflow angle is modified as follows:(43)βin=arctan⁡ν-vh+vgustΩr+Vfcos⁡Ωt-2jπ/N,where vgust is the speed of downward gust.


## 3. Model Validation

To make sure the model has good accuracy, which can be used to solve the problem of aeroelastic response due to gust, validation is conducted in this section. It includes the validation of aerodynamic model, the validation of structural model, and the validation of the whole aeroelastic coupled model.


### 3.1. Aerodynamic Model Validation

To validate the aerodynamic model, the result of this paper is compared with ones of Leishman and Crouse in paper [19]. The results at the condition of stall onset (angle of attack is α=5°+10°sin⁡ωt) and at the condition of deep dynamic stall (angle of attack is α=15°+10°sin⁡ωt) are compared. The compared results are shown in Figures 6 and 7.Figure 6
Coefficients of unsteady normal force and pitching moment during stall onset.
(a)
Coefficient of normal force
(b)
Coefficient of pitching momentFigure 7
Coefficients of unsteady normal force and pitching moment during deep dynamic stall.
(a)
Coefficient of normal force
(b)
Coefficient of pitching momentFrom the compared figures (Figures6 and 7), the aerodynamic model has good accuracy, which can be applied to calculate the nonlinear unsteady aerodynamics caused by large angles of attack.


### 3.2. Structural Model Validation

For the structural model validation, a structural model is chosen form paper [25]. The model parameters are shown in Table 1. The modal frequencies are calculated based on the medium-deflection beam theory and they are compared with the results in this paper, shown in Table 2. Form Table 2, it can be seen that the frequency calculation error is within 5%. Hence, the structural model has a good accuracy. The nondimensional modal shapes with tip deformation of first three flaps and the first two twists are shown in Figure 8.Table 1
Model parameters of the blade.
Element
l
/
R
m
/
m
0
E
I
y
/
m
0
Ω
2
R
4
E
I
z
/
m
0
Ω
2
R
4
G
J
/
m
0
Ω
2
R
4
E
A
/
m
0
Ω
2
R
4
m
k
m
1
2
/
m
0
R
2
m
k
m
2
2
/
m
0
R
2
1 to 10
0.1
1
0.008345
0.023198
0.003822
378.1
0.0001
0.0004Table 2
Modal frequency.
Rotating speed: 383 RPM
This paper (Hz)
This paper (per/rev)
Paper [25] (per/rev)
Error (%)
Lag 1
4.79
0.75
0.747
0.5
Flap 1
7.33
1.15
1.146
0.9
Flap 2
21.7
3.44
3.389
1.5
Lag 2
27.1
4.24
4.315
−1.7
Twist 1
28.2
4.42
4.59
−3.8
Flap 3
46.2
7.24
7.416
−2.4
Twist 2
84.0
13.17
13.6
−3.1Figure 8
The modal shapes of flaps and twists.
(a)
Modal shape of flap
(b)
Modal shape of twist


### 3.3. Trim Validation

The wind tunnel trim state is used as the basis for analysis. The properties of vehicle and blade are shown in Table3.Table 3
Blade and vehicle (box-beam) properties.
Number of blades
4
Radius (m)
4.9377
Hover tip speed (m/s)
198.12
Airfoil
NACA0015
Hub length (xhub/R)
0.04
Solidity
0.1
Lock number
6.34
Mass per unit length (kg/m)
0.135
Aerodynamic root cutout (xroot/R)
0.1
CG below hub (h/R)
0.2
Flat plate area (f/πR2)
0.01
C
T
/
σ
0.07
Precone
0
Tail rotor radius (m)
0.9876
Tail rotor solidity
0.15
Tail rotor location (xtr/R)
1.2
Tail rotor above CG (htr/R)
0.2
Horizontal tail location (xht/R)
0.95
Horizontal tail planform area (Sht/πR2)
0.011For forward flight, the variables that need to be trimmed are three manipulating variables,θ75, θ1c, and θ1s, two attitude angles, αs and ϕs, and θtail. The trim equations for isolated rotors in wind tunnel state are shown as follows:(44)Fx=Fy=Fz=0,Mx=My=Mz=0.The forces and moments due to unsteady aerodynamic model are used to update the vehicle equilibrium equations at each step of iterative trim solution [25, 26]. The flowchart of trim strategy is shown in Figure 9.Figure 9
Flowchart of trim strategy.Trim results of current paper are shown in Figure10. From the results, it can be seen that although the same properties were used to calculate the trim equation, there is still a little difference between current work and the results of [25]. This difference is mainly caused by the difference between the unsteady and steady aerodynamics. The unsteady aerodynamic model is closer to the real situation of rotors.Figure 10
Trim results of control angle in forward flight.


### 3.4. Aeroelastic Model Validation

Based on the trimmed result, a medium advanced ratio of 0.2 of rotor in Section3.3 is chosen to calculate the aeroelastic response. The rotating speed is also 383 RPM. The control angle is θ0=7.1°+1.2°cos⁡ψ-3.1°sin⁡ψ. The rotor for validation is hingeless rotor.The compared result is shown in Figure11. The result of this paper coincides with the one of paper [25], so the coupled strategy of this paper is correct which can be used to investigate the aeroelastic response due to gust.Figure 11
Deformations of flap and twist at the advanced ratio of 0.2.
(a)
Deformations of flap
(b)
Deformations of twist


## 3.1. Aerodynamic Model Validation

To validate the aerodynamic model, the result of this paper is compared with ones of Leishman and Crouse in paper [19]. The results at the condition of stall onset (angle of attack is α=5°+10°sin⁡ωt) and at the condition of deep dynamic stall (angle of attack is α=15°+10°sin⁡ωt) are compared. The compared results are shown in Figures 6 and 7.Figure 6
Coefficients of unsteady normal force and pitching moment during stall onset.
(a)
Coefficient of normal force
(b)
Coefficient of pitching momentFigure 7
Coefficients of unsteady normal force and pitching moment during deep dynamic stall.
(a)
Coefficient of normal force
(b)
Coefficient of pitching momentFrom the compared figures (Figures6 and 7), the aerodynamic model has good accuracy, which can be applied to calculate the nonlinear unsteady aerodynamics caused by large angles of attack.


## 3.2. Structural Model Validation

For the structural model validation, a structural model is chosen form paper [25]. The model parameters are shown in Table 1. The modal frequencies are calculated based on the medium-deflection beam theory and they are compared with the results in this paper, shown in Table 2. Form Table 2, it can be seen that the frequency calculation error is within 5%. Hence, the structural model has a good accuracy. The nondimensional modal shapes with tip deformation of first three flaps and the first two twists are shown in Figure 8.Table 1
Model parameters of the blade.
Element
l
/
R
m
/
m
0
E
I
y
/
m
0
Ω
2
R
4
E
I
z
/
m
0
Ω
2
R
4
G
J
/
m
0
Ω
2
R
4
E
A
/
m
0
Ω
2
R
4
m
k
m
1
2
/
m
0
R
2
m
k
m
2
2
/
m
0
R
2
1 to 10
0.1
1
0.008345
0.023198
0.003822
378.1
0.0001
0.0004Table 2
Modal frequency.
Rotating speed: 383 RPM
This paper (Hz)
This paper (per/rev)
Paper [25] (per/rev)
Error (%)
Lag 1
4.79
0.75
0.747
0.5
Flap 1
7.33
1.15
1.146
0.9
Flap 2
21.7
3.44
3.389
1.5
Lag 2
27.1
4.24
4.315
−1.7
Twist 1
28.2
4.42
4.59
−3.8
Flap 3
46.2
7.24
7.416
−2.4
Twist 2
84.0
13.17
13.6
−3.1Figure 8
The modal shapes of flaps and twists.
(a)
Modal shape of flap
(b)
Modal shape of twist


## 3.3. Trim Validation

The wind tunnel trim state is used as the basis for analysis. The properties of vehicle and blade are shown in Table3.Table 3
Blade and vehicle (box-beam) properties.
Number of blades
4
Radius (m)
4.9377
Hover tip speed (m/s)
198.12
Airfoil
NACA0015
Hub length (xhub/R)
0.04
Solidity
0.1
Lock number
6.34
Mass per unit length (kg/m)
0.135
Aerodynamic root cutout (xroot/R)
0.1
CG below hub (h/R)
0.2
Flat plate area (f/πR2)
0.01
C
T
/
σ
0.07
Precone
0
Tail rotor radius (m)
0.9876
Tail rotor solidity
0.15
Tail rotor location (xtr/R)
1.2
Tail rotor above CG (htr/R)
0.2
Horizontal tail location (xht/R)
0.95
Horizontal tail planform area (Sht/πR2)
0.011For forward flight, the variables that need to be trimmed are three manipulating variables,θ75, θ1c, and θ1s, two attitude angles, αs and ϕs, and θtail. The trim equations for isolated rotors in wind tunnel state are shown as follows:(44)Fx=Fy=Fz=0,Mx=My=Mz=0.The forces and moments due to unsteady aerodynamic model are used to update the vehicle equilibrium equations at each step of iterative trim solution [25, 26]. The flowchart of trim strategy is shown in Figure 9.Figure 9
Flowchart of trim strategy.Trim results of current paper are shown in Figure10. From the results, it can be seen that although the same properties were used to calculate the trim equation, there is still a little difference between current work and the results of [25]. This difference is mainly caused by the difference between the unsteady and steady aerodynamics. The unsteady aerodynamic model is closer to the real situation of rotors.Figure 10
Trim results of control angle in forward flight.


## 3.4. Aeroelastic Model Validation

Based on the trimmed result, a medium advanced ratio of 0.2 of rotor in Section3.3 is chosen to calculate the aeroelastic response. The rotating speed is also 383 RPM. The control angle is θ0=7.1°+1.2°cos⁡ψ-3.1°sin⁡ψ. The rotor for validation is hingeless rotor.The compared result is shown in Figure11. The result of this paper coincides with the one of paper [25], so the coupled strategy of this paper is correct which can be used to investigate the aeroelastic response due to gust.Figure 11
Deformations of flap and twist at the advanced ratio of 0.2.
(a)
Deformations of flap
(b)
Deformations of twist


## 4. Gust Response Analysis

Dynamic gust load of an isolated rotor is investigated with the validated aeroelastic model, coupled with the gust model. The blades are rotating at a speed of 383 RPM. Three different profiles of gust models are shown in the Figure12. They are impulse gust, sinusoidal gust, and the slope gust. The time endurance is also indicated in Figure 12. The gust amplitude is 8 m/s in all the three gust types. Gust responses in hover and forward flights with the advanced ratios of 0.2 and 0.35 are calculated, respectively. Since the equilibrium state which will be used as the initial condition is very important for gust analysis, steady rotor forces and moments are used to update the vehicle equilibrium equations at each step of iterative trim solution [25, 26]. The results of equilibrium state are used in the current work. At the advanced ratio of 0.2, the control angle is θ0=7.1°+1.2°cos⁡ψ-3.1°sin⁡ψ; CT/σ=0.07. At the advanced ratio of 0.35, the control angle is θ0=10.1°+2.1°cos⁡ψ-8°sin⁡ψ; CT/σ=0.07. The time history of the dynamic load and response is calculated in a coupling way shown in Section 2.4.Figure 12
Three different profiles of gust models.


### 4.1. Gust Response in Hover

Figure13 shows the shear force at the root of rotor due to three different gust models. They are also compared with the shear force without gust. Figure 13 compares the thrust coefficient on the rotor disk due to different shape of gust. From Figure 13, when there is no gust in hover, the shear force keeps a fixed value. When the gust models were added in the calculation, the shear force will decrease most as the gust velocity reaches its maximum value. It also can be observed that at the beginning and end of impulse gust the value of shear force is larger than the shear force of the other two gust models due to the instantaneous shock effect. From Figure 14, when the downward gust acts on the rotors, the thrust coefficient on the rotor disk would decrease. Especially when the gust amplitude reaches its maximum value, the thrust coefficient decreases quarter of the value which has no gust. This is mainly because of the decrease of angle of attack. When the rotors are in the hover, more attention should be paid to the decrease of thrust force due to downward gust. Insufficient rotor thrust may lead to helicopter crash. Figure 15 illustrates the tip deflection of rotor due to different shape of gust; it can be seen that the tip deflection will increase when the gust is in low velocity; when the gust reaches its maximum value the tip deflection will decrease.Figure 13
Shear force at the root of rotor.Figure 14
Thrust coefficientC_T on the rotor disk.Figure 15
Tip deflection of rotor due to different shape gust.In order to investigate the effect of gust frequency on the rotor trust force, a multiscaled1-cos shape gust is acted on the rotors. Three segments of this multifrequency gust were investigated as shown in Figure 16. Frequency of first segment, which is acted at the time of 1 s, is 6.3 Hz. It is the same with the rotor rotational frequency. The frequency of the second one is 7.33 Hz which is the same with the frequency of the first flap mode. The frequency of the last one is 0.64 Hz, which is a very low frequency in order to compare with the above two gust segments. Figure 17 shows the tip deflection of rotor due to multiscale gust. Figure 16 shows the variation of thrust coefficient on rotor disk due to multiscale gust. From Figures 17 and 18, it can be observed that the thrust coefficient of disk reaches its same minimum value when three-segment gust reaches its maximum strength. The tip deflection of rotor at high gust frequency is larger than the one at low frequency.Figure 16
Multifrequency sinusoidal gust.Figure 17
Tip deflection of rotor due to multifrequency sinusoidal gust in hover.Figure 18
Thrust coefficient of disk for multiscale gust.


### 4.2. Gust Response at Forward Flight

The rotor tip deflection at the advanced ratio of 0.2 is shown in Figure19. It can be observed that when gust acts on rotor, the flap deflection will increase. Figure 20 shows the shear force at the root of rotors at the advanced ratio of 0.2 due to three kinds of gusts whose maximum strengths are 8 m/s. Figure 21 shows that the thrust coefficient on rotor disk varies when the gust acts on the rotor. It can be observed that when gust acts on rotor the flap deflection will increase. From the curves of 1-cos shape gust and slope gust in Figure 20, it can be observed that when the gust strength increases the shear force at the root will increase at first. However, as the strength of gust increases sequentially, the shear force would decrease when the gust strength reaches its maximum value. This is mainly because when the gust strength is not large enough, the gust induced velocity will be balanced by flap velocity. When the gust strength is large enough, as the increase of gust strength, it will lead to decrease for angle of attack, which will decrease shear force. It can be seen that when the gust was added on the rotors at the advanced ratio of 0.2, the thrust coefficient would decrease. Angle of attack for tip element due to 1-cos gust at the advanced ratio of 0.2 is shown in Figure 22. It can be observed that when the helicopter flights at this advanced ratio, the rotor is at aerodynamic stall onset state. When the gust was acted on rotor, the angle of attack would decrease which would decrease the thrust force of rotor disk.Figure 19
Tip deflection of rotors at the advanced ratio of 0.2.Figure 20
Shear force at root of rotors for advanced ratio of 0.2.Figure 21
Thrust coefficientC_T of disk at the advanced ratio of 0.2.Figure 22
Angle of attack at the tip element due to1-cos gust at the advanced ratio of 0.2.Figure23 shows the frequency response of the flap displacement at the tip due to three gust models. It can be seen that the flap displacement increases obviously when slope gust acts on it. At the rotating speed of 383 RPM, It can be observed that the peak response frequency is approximately 6.38 Hz, which is the rigid rotating frequency. From Figure 23, it is obviously seen that all the peak frequencies of three shapes of gust are nearly the same as the rigid rotating frequency. There is also a peak response frequency which is four times of rigid rotating frequency. The response also can be observed form Figure 21 when no gust was added. It can also be observed that when the gust was added on the rotors, the amplitude of response for the four times rotating frequency would decrease.Figure 23
Frequency response of the flap displacement at the tip due to gusts with three different shapes.Figure24 shows the rotor tip deflection at the advanced ratio of 0.35. Figure 25 shows the thrust coefficient C_T on the rotor disk. Figure 26 shows angle of attack for tip element under 1-cos gust. From Figure 24, it can be observed that the flap deflection of rotors at a higher advanced ratio is larger than lower one. Thrust coefficient of disk of advanced ratio of 0.35 is different from the one at the advanced ratio of 0.2. When the gust acts on the rotor at the advanced ratio of 0.35, the rotor has a large flap velocity. When the gust velocity increases, it may be cancelled out by the flap velocity first. Therefore, the angle of attack may increase first. However, when the velocity of gust increases large enough, the angle of attack decreases afterwards. Since the aerodynamic force of rotor at the advanced ratio of 0.35 is in the state of deep dynamic stall, when the gust strength reaches its maximum value, the angle of attack reaches its minimum value. This is because the angle of attack for rotor due to no gust model is larger than the angle of attack for rotor due to models with three kinds of gust. When the rotor is in deep dynamic stall, the increase of angle will decrease the normal force coefficient which will affect thrust coefficient directly. During the forward flight, especially at a medium advanced ratio, more attention should be paid to both the increase of shear force at root of rotors and decrease of thrust coefficient of disk, because a large shear force may lead to structural damage and a significant decrease of thrust force may lead to crash.Figure 24
Tip deflection at the advanced ratio of 0.35.Figure 25
Thrust coefficientC_T of disk at the advanced ratio of 0.35.Figure 26
Angle of attack for tip element under1-cos gust at the advanced ratio of 0.35.


## 4.1. Gust Response in Hover

Figure13 shows the shear force at the root of rotor due to three different gust models. They are also compared with the shear force without gust. Figure 13 compares the thrust coefficient on the rotor disk due to different shape of gust. From Figure 13, when there is no gust in hover, the shear force keeps a fixed value. When the gust models were added in the calculation, the shear force will decrease most as the gust velocity reaches its maximum value. It also can be observed that at the beginning and end of impulse gust the value of shear force is larger than the shear force of the other two gust models due to the instantaneous shock effect. From Figure 14, when the downward gust acts on the rotors, the thrust coefficient on the rotor disk would decrease. Especially when the gust amplitude reaches its maximum value, the thrust coefficient decreases quarter of the value which has no gust. This is mainly because of the decrease of angle of attack. When the rotors are in the hover, more attention should be paid to the decrease of thrust force due to downward gust. Insufficient rotor thrust may lead to helicopter crash. Figure 15 illustrates the tip deflection of rotor due to different shape of gust; it can be seen that the tip deflection will increase when the gust is in low velocity; when the gust reaches its maximum value the tip deflection will decrease.Figure 13
Shear force at the root of rotor.Figure 14
Thrust coefficientC_T on the rotor disk.Figure 15
Tip deflection of rotor due to different shape gust.In order to investigate the effect of gust frequency on the rotor trust force, a multiscaled1-cos shape gust is acted on the rotors. Three segments of this multifrequency gust were investigated as shown in Figure 16. Frequency of first segment, which is acted at the time of 1 s, is 6.3 Hz. It is the same with the rotor rotational frequency. The frequency of the second one is 7.33 Hz which is the same with the frequency of the first flap mode. The frequency of the last one is 0.64 Hz, which is a very low frequency in order to compare with the above two gust segments. Figure 17 shows the tip deflection of rotor due to multiscale gust. Figure 16 shows the variation of thrust coefficient on rotor disk due to multiscale gust. From Figures 17 and 18, it can be observed that the thrust coefficient of disk reaches its same minimum value when three-segment gust reaches its maximum strength. The tip deflection of rotor at high gust frequency is larger than the one at low frequency.Figure 16
Multifrequency sinusoidal gust.Figure 17
Tip deflection of rotor due to multifrequency sinusoidal gust in hover.Figure 18
Thrust coefficient of disk for multiscale gust.


## 4.2. Gust Response at Forward Flight

The rotor tip deflection at the advanced ratio of 0.2 is shown in Figure19. It can be observed that when gust acts on rotor, the flap deflection will increase. Figure 20 shows the shear force at the root of rotors at the advanced ratio of 0.2 due to three kinds of gusts whose maximum strengths are 8 m/s. Figure 21 shows that the thrust coefficient on rotor disk varies when the gust acts on the rotor. It can be observed that when gust acts on rotor the flap deflection will increase. From the curves of 1-cos shape gust and slope gust in Figure 20, it can be observed that when the gust strength increases the shear force at the root will increase at first. However, as the strength of gust increases sequentially, the shear force would decrease when the gust strength reaches its maximum value. This is mainly because when the gust strength is not large enough, the gust induced velocity will be balanced by flap velocity. When the gust strength is large enough, as the increase of gust strength, it will lead to decrease for angle of attack, which will decrease shear force. It can be seen that when the gust was added on the rotors at the advanced ratio of 0.2, the thrust coefficient would decrease. Angle of attack for tip element due to 1-cos gust at the advanced ratio of 0.2 is shown in Figure 22. It can be observed that when the helicopter flights at this advanced ratio, the rotor is at aerodynamic stall onset state. When the gust was acted on rotor, the angle of attack would decrease which would decrease the thrust force of rotor disk.Figure 19
Tip deflection of rotors at the advanced ratio of 0.2.Figure 20
Shear force at root of rotors for advanced ratio of 0.2.Figure 21
Thrust coefficientC_T of disk at the advanced ratio of 0.2.Figure 22
Angle of attack at the tip element due to1-cos gust at the advanced ratio of 0.2.Figure23 shows the frequency response of the flap displacement at the tip due to three gust models. It can be seen that the flap displacement increases obviously when slope gust acts on it. At the rotating speed of 383 RPM, It can be observed that the peak response frequency is approximately 6.38 Hz, which is the rigid rotating frequency. From Figure 23, it is obviously seen that all the peak frequencies of three shapes of gust are nearly the same as the rigid rotating frequency. There is also a peak response frequency which is four times of rigid rotating frequency. The response also can be observed form Figure 21 when no gust was added. It can also be observed that when the gust was added on the rotors, the amplitude of response for the four times rotating frequency would decrease.Figure 23
Frequency response of the flap displacement at the tip due to gusts with three different shapes.Figure24 shows the rotor tip deflection at the advanced ratio of 0.35. Figure 25 shows the thrust coefficient C_T on the rotor disk. Figure 26 shows angle of attack for tip element under 1-cos gust. From Figure 24, it can be observed that the flap deflection of rotors at a higher advanced ratio is larger than lower one. Thrust coefficient of disk of advanced ratio of 0.35 is different from the one at the advanced ratio of 0.2. When the gust acts on the rotor at the advanced ratio of 0.35, the rotor has a large flap velocity. When the gust velocity increases, it may be cancelled out by the flap velocity first. Therefore, the angle of attack may increase first. However, when the velocity of gust increases large enough, the angle of attack decreases afterwards. Since the aerodynamic force of rotor at the advanced ratio of 0.35 is in the state of deep dynamic stall, when the gust strength reaches its maximum value, the angle of attack reaches its minimum value. This is because the angle of attack for rotor due to no gust model is larger than the angle of attack for rotor due to models with three kinds of gust. When the rotor is in deep dynamic stall, the increase of angle will decrease the normal force coefficient which will affect thrust coefficient directly. During the forward flight, especially at a medium advanced ratio, more attention should be paid to both the increase of shear force at root of rotors and decrease of thrust coefficient of disk, because a large shear force may lead to structural damage and a significant decrease of thrust force may lead to crash.Figure 24
Tip deflection at the advanced ratio of 0.35.Figure 25
Thrust coefficientC_T of disk at the advanced ratio of 0.35.Figure 26
Angle of attack for tip element under1-cos gust at the advanced ratio of 0.35.


## 5. Conclusions

A structural dynamic model and a consequent aeroelastic model for an isolated helicopter rotor are constructed. The dynamic load induced by impulse gust, slope-shape gust and1-cos shape gust is investigated, respectively. We conclude the following:(1)
The structural dynamic model based on the principle of Hamilton is accurate and is validated for dynamic analysis of rotors. The aerodynamic model is constructed in this paper and is validated. The result of validation shows that the model is accurate to calculate unsteady aerodynamics.(2)
The coupled strategy introduced in this paper is validated; the result of validation shows that coupled strategy is reasonable and the aeroelastic model is accurate to investigate gust response.(3)
In the hover flight, the thrust force may decrease when the gust strength increases. A lower rotor thrust force may lead to helicopter crash. The shear force would also decrease as the gust strength increases. More attention should be paid to the decrease of thrust force in hover.(4)
At the forward flight, the shear force may increase as the gust strength increases. The thrust force may also decrease when the gust acts on rotor, which is similar to the hover flight condition. More attention should be paid to the structural shear force, which increases with the gust strength at the forward flight.(5)
The dynamic load due to gust is periodic with the harmonics of rotational speed of rotors. The peak response is found at the first and the fourth harmonics of rigid rotational frequency. The response is not significant at the elastic modal frequencies.



---

*Source: 8986217-2017-06-19.xml*