# Nonlinear Characteristics of Helicopter Rotor Blade Airfoils: An Analytical Evaluation

Authors: Constantin Rotaru
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Some results are presented  about the study of airloads of the helicopter rotor blades, the aerodynamic characteristics of airfoil sections, the physical features, and the techniques for modeling the unsteady effects found on airfoil operating under nominally attached flow conditions away from stall. The unsteady problem was approached on the basis of Theodorsen's theory, where the aerodynamic response (lift and pitching moment) is considered as a sum of  noncirculatory and circulatory parts. The noncirculatory or apparent mass accounts for the pressure forces required to accelerate the fluid in the vicinity of the airfoil. The apparent mass contributions to the forces and pitching moments, which are proportional to the instantaneous motion, are included as part of the quasi-steady result.


---

## Body

## 1. Introduction

The most important component of the helicopter is the main rotor for which there is a great deal of activity in developing new and improved mathematical models that predict the flow physics. A high tip speed gives the rotor a high level of stored rotational kinetic energy and reduces the rotor torque required for a given power, but there are two important factors that work against the use of a high tip speed: compressibility effects and noise.The additional effects of compressibility on the overall rotor profile power requirements, when the tip of the advancing blade approaches and exceeds the drag divergence Mach number, were estimated using blade element theory combined with the airfoil section characteristics [1]. A more detailed analysis of compressibility effects on the rotor must represent the actual nonlinear airfoil characteristics as functions of Mach number through stall at each blade element followed by numerical integration.The classical unsteady aerodynamic theories describing the observed behavior have formed the basis for many types of rotor analysis. The tools for the analysis of 2D, incompressible, and unsteady aerodynamic problems were extended to compressible flows, being a basis for developing linearized unsteady aerodynamic models applicable to compressible flows [2]. But, while the classical theories assume linearity in the airloads, the assumption of linearity can probably be justified for many of the problems encountered on the rotor, in practice. The advent of nonlinear methods based on CFD solutions to the Euler and Navier-Stokes equations has provided new results that justify and define the limits of the linear models and may give guidance in developing improved and more practical unsteady aerodynamic models for future use in helicopter rotor blade airloads prediction, aeroelastic analysis, and rotor design [3]. At the blade element level, the various sources of unsteady effects can be decomposed into perturbations to the local angle of attack and velocity field. At low angle of attack with fully attached flow, the various sources of unsteady effects manifest as moderate amplitude and phase variations relative to the quasi-steady airloads. At higher angles of attack when time-dependent flow separation from the airfoil may be involved, a phenomenon characterized by large overshoots in the values of the lift, drag, and pitching moments relative to the quasi-steady stall  values may occur [4].The helicopter rotor airfoil must assure a high maximum lift coefficient, a high drag divergence Mach number, a good lift-to-drag ratio over a wide range of Mach number, and a low pitching moment. At high angles of attack, the adverse pressure gradients produced on the upper surface of the airfoil result in a progressive increase in the thickness of the boundary layer and cause some deviation from the linear lift versus angle of attack behavior. On many airfoils, the onset of flow separation and stall occurs gradually with increasing angle of attack, but on some airfoils (those with sharp leading edges) the flow separation may occur suddenly.


## 2. Helicopter Rotor Blade Aerodynamics

The region of the rotor disk affected by compressibility effects is shown in Figure1 and is defined on the surface where the incident Mach number of the flow that is normal to the leading edge of the blade exceeds the drag divergence Mach number, Mdd. If MΩR is the hover tip Mach number, then the region of the disk affected by compressibility effects is defined by
(1)Mr,ψ=MΩR(r+μsinψ)≥Mdd.Figure 1
Helicopter rotor blade in forward flight.The angular or rotational speed of the rotor is denoted byΩ, the rotor radius by R, the advanced ratio by μ and the radial distance from the rotational axis by y. In the above equation the nondimensional quantities  are r=y/R and μ=V∞cosα/ΩR, where α is the angle between the forward velocity V∞ and the plane of the rotor [1, 2]. The azimuth angle for the onset drag divergence, ψ1, can be obtained by setting r=1, so that
(2)ψ1=arc sin([1μ(MddMΩR-1)])
and ψ2=180-ψ1.The increment in profile powerΔCP associated with this region on the disk is
(3)ΔCPσ=14π∫ψ1ψ2∫rdd1(r+μsinψ)3ΔCdrdrdψ,
where ΔCd is the extra drag on the blade section when it exceeds the drag divergence Mach number, Mdd, and σ is rotor solidity coefficient which represents the ratio of the blades area  to the rotor disk area. For the NACA 0012 airfoil, Prouty [4] suggests that this can be approximated by
(4)ΔCd(M)={12.5(M-0.74)3forM≥0.74,0otherwise.
The rotor limits may be determined by two conditions, one condition given by advancing blade compressibility effects and the other condition given by retreating blade stall. In either case the advancing blade operates at low angle of attack (AoA) but at high subsonic or transonic conditions, whereas the retreating blade operates at low Mach numbers and high lift coefficients.The aerodynamic behavior of airfoils in the high AoA regime is important for predicting the adverse effects produced in the reverse flow regime on the rotor. In the reverse flow region, the direction of the relative flow vector changes from the trailing edge toward the leading edge of the airfoil. While the fundamental process of the blade wake and tip vortex formation is similar to that found with a fixed wing, one difference with helicopter tip vortices is that they are curved, and so they experience a self-induced effect. Another complication with helicopter rotors is that the wakes and tip vortices from other blades can lie close to each other and to the plane of blade rotation, and so they have large induced effects on the blade lift distribution.If the wake is assumed to be undistorted in the tip path plane and no wake contraction occurs in the radial direction (Figure2), then the tip vortex trajectories are described by the equations
(5)x=Rcos(ψb-ψw)+Rμψw,y=Rsin(ψb-ψw),
where ψb is the position of the blade when the vortex was formed and ψw is the position of the vortex element relative to the blade.Figure 2
Tip vortex trajectory.One important parameter used in the description of unsteady aerodynamics and unsteady airfoil behavior is the reduced frequency,k, defined as k=ω·c/(2V), where ω is the angular frequency, c is the chord of the airfoil, and V is the flow velocity [3]. According to the dimensional analysis, the resultant force, F, on the airfoil chord c can be written in functional form as F/(ρV2c2)=f(Re,M,k). For k=0, the flow is steady, and for 0≤k≤0.05, the flow can be considered quasi-steady; that is, unsteady effects are generally small. Flows with characteristic reduced frequencies above 0.05 are considered unsteady [5, 6]. For a helicopter rotor in forward flight (Figure 3), the local sectional velocity, which appears in the denominator of the reduced frequency expression, is constantly changing.Figure 3
Helicopter main rotor.The approach to modeling of unsteady aerodynamic effects through an extension of steady, 2D thin airfoil theory gives a good level of analysis of the problem and provides considerable insight into the physics responsible for the underlying unsteady behavior. The Laplace’s equation for incompressible flow is eliptic; therefore, the unsteady aerodynamic theories cannot be obtained in a corresponding analytical form.


## 3. The Apparent Mass Tensor

The rate of change of the impulse vector, in general, is not in the direction of the acceleration of the body. The external forceFe applied to the body to translate it through the fluid has to be applied in a direction different from that of the acceleration of the body through the fluid [7]. Physical conditions that should be satisfied on given boundaries of the fluid (boundary conditions) depend on the assumptions made with regard to the nature of the fluid, more specifically on the nature of the differential equations that are assumed to govern the motion of the fluid. For a solid-fluid boundary, at each point of the solid-fluid surface, at every instant, the component normal to the surface of the relative velocity between the fluid and the solid must be vanish, V→·n→=0, where V→ represents the relative velocity and n→ the normal to the surface (Figure 4).Figure 4
Solid-fluid surface.If the surface is represented by a scalar function of position and time,F(r→,t)=0, then the total time rate of change is zero:
(6)D(F)Dt=∂(F)∂t+V→·grad(F)=0
on F(r→,t)=0.The fluid force acting on a rigid body of arbitrary shape translating with a velocityU→(t) is given by
(7)F→=-∬Spn→dS,
where S denotes the surface of the body and p is the pressure on the surface of the body. In general, the body may be translating, rotating, and deforming; consequently, the velocity U is a function of position on the surface and time. If the body is rigid and is in translatory motion, then U is a function of time but uniform over the surface of the body. The mathematical problem is to determine the external force F→e applied to the body to translate it through the fluid. According to Newton’s second law, we have
(8)ddt(mU→)=F→e+F→,
where m is the mass of the body. The above equation may be rewritten as
(9)F→e=ddt(mU→)-F→
or
(10)F→e=ddt(mU→+I→),
where I→ is the impulse applied on the fluid and -dI→/dt=F.→ The fluid force acting on the body is
(11)F→=∂∂t∯Sρϕn→dS-ρU→×∯Sn→×grad(ϕ)dS.
The integral
(12)I→C=∯Sn→×grad(ϕ)dS
is related to the circulation C around the body, and ϕ is the velocity potential.The unit vectore→ is normal to the cutting planes, the unit vector e→1 is tangent to the curve of intersection between the blade element surface and the cutting plane, and the unit vector e→2 is tangent to the blade element surface.The component of the vectorI→C in the direction e→ is
(13)e→·I→C=∯Se→·n→×q→dS,
and the vector n→dS corresponding to the surface element dS=dx·dy may be written as n→dS=dx→×dy→, where dx→=dxe→1, dy→=dye→2, and q=grad(ϕ).On the other hand,(14)n→dS×q→=(dx→×dy→)×q→=(q→·dx→)dy→-(q→·dy→)dx→,e→·n→×q→dS=e→(q→·dx→)dy→-e→(q→·dy→)dx→.
Since the unit vectors e→ and e→1are normal, it follows that e→·dx→=0 and
(15)e→·(n→×q→)dS=(q→·dx→)(e→·dy→).
The scalar product e→·dy→, is the normal distance between the cutting planes of the solid body (Figure 5). If we denote dh=e→·dy→, it follows that
(16)e→·∯Sn→×q→dS=∫h1h2(∮Cq→·dx→)dh=∫h1h2Γe(h)dh.
Here, h is the distance measured along the fixed direction e→, and
(17)Γe(h)=∮Cq→·dx→
is the circulation around the curve of intersection between the body surface and the cutting plane. The limits h1 and h2 denote the extremities of the body measured along the direction e→. It follows that for motions without circulation the force on the body is given by
(18)F→=∂∂t(∯Sρϕn→dS),
where the velocity potential ϕ is the solution of the system
(19)∇2ϕ=0,gradϕ·n→=∂ϕ∂n=U(t)nonS.
Since the equation and boundary condition for ϕ are linear, the solution could have the form
(20)ϕ=ϕ1+ϕ2+ϕ3,
where each of the functions ϕ1, ϕ2, and ϕ3 is a solution of the equations
(21)∇2ϕi=0,grad(ϕi)·n→=∂ϕi∂n=uinionS,
where i=1, 2, or 3.Figure 5
Rotor blade element.In the Cartesian coordinate system, the vectorsU→ and n→ have the expressions
(22)U→=u1i→+u2j→+u3k,→n→=n1i→+n2j→+n3k→.
Because time enters through ui, it is convenient to set ϕi=uiφi, so the system (21) takes the form
(23)∇2φi=0,grad(φi)·n→=∂φi∂n=nionS.
With these considerations the impulse I→ becomes
(24)-I→=∯Sρϕn→dS=∯Sρ(∑k=13ukφk)n→dS=∑k=13(∯Sρϕkn→dS)uk.
The components of the impulse I→=I1i→+I2j→+I3k→ are
(25)I1=j→·I→=∑k=13(-∯Sρφkn1dS)uk,I2=j→·I→=∑k=13(-∯Sρφkn2dS)uk,I3=k→·I→=∑k=13(-∯Sρφkn3dS)uk.
The surface integral in above equations may be written as follows:
(26)∯SρφknidS=∯Sρφk∂φi∂ndS.
According to Green’s theorem, if f1 and f2 are two harmonic functions, then
(27)∯Sf1∂f2∂ndS=∯Sf2∂f1∂ndS.
Introducing the symbol mki,
(28)mki=-∯Sρφk∂φi∂ndS
with mki=mik; the components of the impulse I→ are therefore given by
(29)Ii=∑k=13mikuk,i=1,2,3,
and the force applied externally to the body is
(30)F→e=ddtm(u1i→+u2j→+u3k→)+ddt(m11u1+m12u2+m13u3)i→+ddt(m11u1+m12u2+m13u3)j→+ddt(m11u1+m12u2+m13u3)k→
or
(31)F→e=[(m+m11)du1dt+m12du2dt+m13du3dt]i→+[m21du1dt+(m+m22)du2dt+m23du3dt]j→+[m31du1dt+m32du2dt+(m+m33)du3dt]k.→
The coefficients mik form a set of nine numbers which may be displayed as an array
(32)(m11m12m13m21m22m23m31m32m33)
and may be referred to as a virtual mass tensor or virtual masses that need to be added to the mass of the body in order to find the force that must be applied to translate it through the fluid. Introducing the symbol δik defined by δik=0 if i≠k and δik=1 for i=k, (30) may be rewritten as
(33)(Fe)i=∑k=13(mδik+mik)dukdt.
For any body, there are three perpendicular directions such that mik=0 for i≠k, so with respect to such axes, (33) becomes
(34)(Fe)i=(m+mii)duidt,i=1,2,3.
The sum (m+mii) represents the apparent mass for translation in the i-direction, and the corresponding mii is the additional apparent mass.


## 4. The Airloads on an Oscillating Airfoil

The oscillatory motion of the airfoil can be decomposed into contributions associated with angle of attack which is equivalent to a pure plunging motion (Figure6) and contributions associated with pitching (Figure 7).Figure 6
Plunge velocity.Figure 7
Pitch rate.A plunge velocityh˙ produces a uniform velocity perturbation w that is normal to the chord, w(x)=-h˙, and the pitch-rate term produces a linear variation in normal perturbation velocity.For a pitch rate imposed about an axis at “a” semichords from the midchord, then w(x)=-α˙(x-a), so that the induced chamber is a parabolic arc.The problem of finding the airloads on an oscillating airfoil was solved by Theodorsen, who gave a solution to the unsteady airloads on a 2D harmonically oscillated airfoil in inviscid, incompressible flow, with the assumption of small disturbances [6]. Both the airfoil and its shed wake were represented by a vortex sheet with the shed wake extending as a planar surface from the trailing edge downstream to infinity. The assumption of planar wake is justified if the angle of attack disturbances remain relatively small. As with the standard quasi-steady thin airfoil theory, the bound vorticity, γb, can sustain a pressure difference and, therefore, a lift force. The wake vorticity, γw, must be force-free with zero net pressure jump over the sheet. According to the Theodorsen’s theory, the solution for the loading γb on the airfoil surface under harmonic forcing conditions is obtained from integral equation
(35)w(x,t)=12π∫0cγb(x,t)x-x0dx+12π∫c∞γw(x,t)x-x0dx,
where w is the downwash on the airfoil surface. At the trailing edge, γb(c,t)=0, and the airfoil circulation Γ(t) is given by
(36)Γ(t)=∫0cγb(x,t)dx.
So long as the circulation about the airfoil is changing with respect to time, the circulation is continuously shed into the wake and will continuously affect the aerodynamic loads on the airfoil.For a general motion, where an airfoil of chordc=2b is undergoing a combination of pitching (α,α˙) and plunging (h) motion in a flow of steady velocity V, Theodorsen’s solution [8–10] for the lift coefficient and pitching moment coefficient corresponding to midchord, M1/2, is
(37)cl=πb[h¨V2+α˙V-bV2aα¨]+2π[h˙V+α+bα˙V(12-a)]C(k),cm1/2=π2[bah¨V2-b2V2(18+a2)α¨]+π(a+12)[h˙V+α+b(12-a)α˙V]C(k)-π2[(12-a)bα˙V],
where a is the pitch axis location relative to the midchord of the airfoil, measured in terms of semichord and C(k)=F(k)+iG(k) is the complex transfer function (known as Theodorsen’s function) which accounts for the effects of the shed wake on the unsteady airloads:
(38)C(k)=H1(2)(k)H1(2)(k)+i·H0(2)(k)=J1(J1+J0)+Y1(Y1-J0)(J1+Y0)2+(Y0-J1)2+iY1Y0+J1J0(J1+Y0)2+(Y0-J1)2
with J0,  J1,  Y0, and Y1 being Bessel functions of the first and second kinds, respectively (Figure 8).Figure 8
Bessel functions.The Hankel functions in above expression are(39)H0(2)=J0-i·Y0,H1(2)=J1-i·Y1.
The real and imaginary parts of C(k) function are plotted in Figure 9.Figure 9
Theodorsen’s function.It could be appreciated thatC(k) function serves to introduce an amplitude reduction and phase lag effect on the circulatory part of the lift response compared to the result obtained under quasi-steady conditions [11].This effect can be seen if a pure oscillatory variation in angle of attack is considered, that is,α=α-eiωt, so the circulatory part of the airfoil lift coefficient is given by
(40)cl=2πα-C(k)=2πα-[F(k)+iG(k)].
For k=0, the steady-state lift behavior is obtained; that is, cl is linearly  proportional to α. As k is increased, the lift plots develop into hysteresis loops, and these loops rotate such that the amplitude of the lift response (half of the peak-to-peak value) decreases with increasing reduced frequency.These loops are circumvented in a counterclockwise direction such that the lift is lower than the steady value, whenα is decreasing with time (i.e., there is a phase lag). For infinite reduced frequency, the circulatory part of the lift amplitude is half that at k=0, and there is no phase lag angle.


## 5. Pure Angle of Attack Oscillations

For a harmonic variation inα, that is, α=α-eiωt, the lift is
(41)L=2πρV2b[C(k)+12iωbV]α-eiωt
or, in terms of the lift coefficient, the result is
(42)cl=LρV2b=[2π(F+iG)+iπk]α-eiωt.
The term inside the square brackets can be considered the lift transfer function, which accounts for the difference between the unsteady and quasi-steady airloads.The first term inside the brackets is the circulatory term, and the second term is the apparent mass contribution, which is proportional to the reduced frequency and leads the forcing by a phase angle ofπ/2. The noncirculatory or apparent mass terms arise from the nonsteady term of the velocity potential and account for the pressure forces required to accelerate the fluid in the vicinity of the airfoil.The normalized lift amplitude is(43)cl2πα=(F+iG)+ik2.
The normalized lift amplitude and phase of lift for pure angle of attack oscillations are presented in Figures 10 and 11, where the significance of the apparent mass contribution to both the amplitude and phase can be appreciated.Figure 10
Normalized lift amplitude.Figure 11
Phase angle.At lower values of reduced frequency, the circulatory terms dominate the solution. At higher values of reduced frequency, the apparent mass forces dominate.


## 6. Pure Plunging Oscillations

For a harmonic plunging motion such as that contributed by blade flapping, the forcing ish=h-eiωt so that h˙=iωh-eiωt and h¨=-ω2h-eiωt. Substituting into the expression for the lift and solving for the lift coefficient, it gives
(44)cl=[2πk(iF-G)-πk2]h-beiωt.
The complete term inside the square brackets can be considered as the lift transfer function [12]. The circulatory part of the lift response leads the forcing displacement h by a phase angle of π/2. Also, the apparent mass force leads the circulatory part of the response by a phase angle of π/2 or the forcing by a phase angle of π. The corresponding pitching moment about midchord for this case is
(45)cm1/2=(π4)k2h-beiωt.


## 7. Pitching Oscillations

For harmonic pitch oscillations, additional terms involving pitch rateα˙ appear in the equations for the aerodynamic response. The forcing is given by α=α-eiωt and the pitch rate by α˙=iωα-eiωt. In this case, the lift coefficient is
(46)cl=2π[F(1+ik)+G(i-k)]α-eiωt+πk(i-k2)α-eiωt.
The lift amplitude initially decreases with increasing k because of the effects of the shed wake; and then, for k>0.5, it begins to increase, as the apparent mass forces begin to dominate the airloads. This is also shown by the phase angle, which exhibits an increasing lead for k>0.3.Von Karman and Sear analyzed the problem of a thin airfoil moving through a sinusoidal vertical gust field, where the gust can be considered as an upwash velocity that is uniformly convected by the free stream. The forcing function in this case is(47)wg(x,t)=sin(ωgt-ωgxV),
where ωg is the gust frequency. If the gust is referenced to the airfoil leading edge, then x=0 and wg(x,t) becomes wg(x,t)=sin(ωgt), and if the gust is referenced to the midchord, then x=b=c/2 and
(48)ωgxV=ωg(c/2)V=ωgc2V=kg.
Therefore, wg(t)=sin(ωgt)·cos(kg)-sin(kg)·cos(ωgt), which is equivalent to a phase shift. In this case, the lift coefficient can be written as
(49)cl=2π(w0V)S(kg),
where S(kg) is known as Sears function and the gust encounter frequency, kg, is given by
(50)kg=2πVλg,
and λg is the wavelength of the gust.In terms of Bessel functions, Sears’s function is given by(51)S(kg)=[J0(kg)-iJ1(kg)]C(kg)+iJ1(kg).
The terms of real and imaginary parts are
(52)RealS(kg)=F(kg)·J0(kg)+G(kg)·J1(kg),ImS(kg)=G(kg)·J0(kg)-F(kg)·J1(kg)+J1(kg).
If the gust is referenced to the leading edge of the airfoil, the function will be called S′ and can be written as
(53)RealS′(kg)=RealS(kg)·cos(kg)+ImS(kg)·sin(kg),ImS′(kg)=-RealS(kg)·sin(kg)+ImS(kg)·cos(kg).
The two results are plotted in Figure 12. The peculiar spiral shape of the S transfer function arises only when the gust front is referenced to the midchord of the airfoil.Figure 12
Sears function.The kinematics of the pitching and plunging airfoil of a typical blade element is the resultant of a combination of forcing from collective and cyclic blade pitch, twist angle, elastic torsion, blade flapping velocity, and elastic bending. At low angles of attack with fully attached flow, the various sources of unsteady effects manifest primarily as moderate amplitude and phase variations relative to the quasi-steady airloads. At higher angles of attack when time-dependent flow separation from the airfoil may be involved, the dynamic stall may occur. The amplitude and phase effects produced by the stalled airloads can lead to various aeroelastic problems on the helicopter rotor that may limit its performance. The need to control the aerodynamic forces on the rotor requires that the pitch of each blade be changed individually as the blades rotate about the shaft.


## 8. Summary

The first flap frequency of a helicopter rigid blade is about1.05Ω, then the reduced frequency at the 75% radius location is k=0.7c/R. For a blade aspect ratio R/c>10, the reduced frequency is in excess of 0.07, which is in the unsteady range. In the case of the first elastic torsion mode, which is typically about 3÷4Ω at the blade tip, the reduced frequency is in excess of 0.2. At these reduced frequencies, there is a significant amplitude and phasing introduced into the airloads by the effects of the unsteady aerodynamics.When a wing’s angle of attack is increased rapidly, it can momentarily generate a higher maximum lift coefficient than it could if the angle of attack was increased slowly. This overshoot can be related to the change in angle of attack during the time required for the air to travel one chord length. The dynamic overshoot is attributed to two effects (for the airfoils that stall first at the leading edge): the delay in the separation of the boundary layer and the momentary existence of a vortex shed at the leading edge after the boundary layer does separate. The delay in separation corresponds to the finite time required for the aft edge of the separation bubble to move forward to its bursting position. On the other hands, the airfoil can generate high lift as a result of a vortex that is shed at the leading edge at the instant of stall. The vortex travels back over the top of the airfoil carrying with it a low pressure wave that accounts for the very large lift coefficient. Airfoils that stall first at the trailing edge also exhibit a dynamic overshoot but considerably less than those airfoils that have leading edge stall.



---

*Source: 503858-2013-10-28.xml*