Verication of a Modelica Helicopter Rotor Model Using Blade Element Theory

(1)

Linköping University | Department of Management and Engineering Master’s thesis, 30 credits| MSc Aeronautical Engineering Autumn 2017| LIU-IEI-TEK-A--17/02958—SE

Verification of a

Modelica Helicopter

Rotor Model Using

Blade Element Theory

Author: Jorge Luis Lovaco Hernandez

Examiner: Ingo Staack

Supervisor: Magnus Sethson

Linköping University

(2)

Verification of a

Modelica Helicopter

Rotor Model Using

Blade Element Theory

Examiner: Ingo Staack

Supervisor: Magnus Sethson

(3)

Abstract

Helicopters have been valuable vehicles ever since their invention. Their capabilities for axial flight and hovering make them an outstanding resource. However, their complexity, directly related to their aerodynamics, makes them extremely hard to design. In today’s market competitivity resources must be optimized and accurate models are needed to obtain realizable designs.

The well known Blade Element Theory was used to model helicopter rotors using the Modelica based software SystemModeler. However, it remained unverified due to the lack of experimental data available.

The access to experimental data published by NASA motivated the comparison from the model to the measurements obtained during real testing to a scaled rotor. Some improvements were performed to the model obtaining unexpectedly accurate results for hover and axial flight. Two approaches based on the Blade Element Theory and related to Vortex Theory were followed: an infinite number of blades and a finite number of blades. Moreover, the model simulation speed was noticeably increased and prepared for the forward flight model development.

Nonetheless, even though the model was ready for forward flight simulations, further research is needed due to, again, the lack of experimental data available.

It is concluded from the present work that Wolfram’s SystemModeler can be used as a tool in early design phases of helicopters, even before CAD modeling and CFD due to its simplicity, speed, accuracy, and especially its capability for being used on simple desktop computers.

(4)

Acknowledgements

The present Master Thesis was carried out at Wolfram-Mathcore in Link¨oping under the supervision of Link¨oping University from January to June in 2017. I would like to acknowledge here all the people that gave me help at some point during these months.

I need to thank all the people in Wolfram-Mathcore their support, kindness and help during this time, with a special mention to Jan Brug˚ard for offering me this opportunity at the company, Markus Dahl for his endless patience and support as supervisor, and Leonardo Laguna for his valuable advices.

I also need to thank my supervisor, Magnus Sethson, for his help and sharing all his knowledge, and also my examiner, Ingo Staack, for his valuable comments and pointing out the right direction. I also need to mention Roland G˚ardhagen and Lars Johansson for replying my questions when I needed even though they were not directly involved.

(5)

Nomenclature

Abbreviations and Symbols

Abbreviation Meaning

A Area

A, A’ Axial Force, Axial Force per unit span

B Tip Loss Factor

c Chord

ca, CA Axial Coefficient per unit span, Axial Force Coefficient

cd, CD Drag Coefficient per unit span, Drag Coefficient

cf Friction Coefficient

clα Lift Slope

cl, CL Lift Coefficient per unit span, Lift Coefficient

cn, CN Normal Force Coefficient per unit span, Normal Force Coefficient

CP Pressure Coefficient

CQ Torque Coefficient

CT Thrust Coefficient

CFD Computational Fluid Dynamics

D Drag Force

F Force

F Prandtl’s Function

L, L’ Lift Force, Lift Force per unit span

LE Leading Edge

˙

m Mass Flow Rate

M Mach Number

M∞ Freestream Mach Number

N, N’ Normal Force, Normal Force per unit span

p Pressure p∞ Freestream Pressure P Power Q Torque r Dimensionless Radius R Radius Re Reynolds Number S Area T Thrust TE Trailing Edge U Resultant Velocity UT In Plane Velocity

UP Out of Plane Velocity

UR Radial Velocity

vi Induce Velocity

Vr Radial Velocity

Vθ Tangential Velocity

V∞ Freestream Velocity

Vtip Blade Tip Velocity

w Far Wake Flow Velocity

(6)

Symbols Meaning

α Angle of Attack

β Cone Angle

Γ Circulation

ξ Vorticity

θ Collective Pitch Angle λ Total Induced Inflow Ratio λh Hover Induced Inflow Ratio

λi Induced Inflow Ratio

µ Advance Ratio µ Viscosity ρ Density ρ∞ Freestream Density σ Solidity τ Shear Stress

φ Relative Inflow Angle

χ Wake Skew Angle

ψ Azimuthal Angle

ω Angular Velocity Ω Angular Velocity

(7)

(8)

List of Figures

1 Helicopter rotor characteristics. See Fig.1-4 in [Johnson, 1994] and

Fig.2-1 in[Leishman, 2006] respectively. . . 3

2 Aerodynamic forces on a lifting body (see figures(1.16) and (1.18) in [John D. Anderson, 2007]). a)Resultant forces b)Pressure and shear stress distribution . . . 7

3 Integration criteria for body surfaces (see [John D. Anderson, 2007]), adapted from figure(1.17) . . . 8

4 Circulation for a wing profile in a flow . . . 15

5 Compressibility corrections comparison . . . 18

6 Blade element velocity components (see [Leishman, 2006], adapted from figure(3.1) . . . 21

7 Rotor disk. Adapted from figure(3.4) - [Leishman, 2006] . . . 23

8 Forward flight model. Adapted from Figure(2.23) in [Leishman, 2006] 25 9 a)Modelica based block b)SystemModeler’s resultant animation . . . 29

10 Climbing disk with simple control system . . . 31

11 Tip path . . . 31

12 General model library . . . 32

13 NACA0015 Lift coefficient curve for Reynolds=43000 . . . 34

14 NACA0015 Lift slopes approximation criteria . . . 34

15 CFD Blade domain box . . . 36

16 NACA0015 Drag coefficient approximation for Reynolds=43000 . . . 37

17 Collective pitch at 15.3 degrees Ct comparison . . . 41

18 Collective pitch at 15.3 degrees Cp comparison . . . 41

19 Collective pitch at 16.9 degrees Ct comparison . . . 42

20 Collective pitch at 16.9 degrees Cp comparison . . . 43

21 Stalled blade streamlines, collective pitch at 15.3 degrees . . . 44

22 Reattached flow streamlines, collective pitch at 15.3 degrees . . . 44

23 Flow streamlines at the 75% blade length, collective pitch at 15.3 degrees . . . 45

24 Stalled blade streamlines, collective pitch at 16.9 degrees . . . 46

25 Reattached flow streamlines, collective pitch at 16.9 degrees . . . 46

26 Flow streamlines at the 75% blade length, collective pitch at 16.9 degrees . . . 47

27 Ct values depending on Lift-slope approximations at 15.3 degrees . . 48

28 Cp values depending on Lift-slope approximations at 15.3 degrees . . 48

29 Ct values depending on Lift-slope approximations at 16.9 degrees . . 49

30 Cp values depending on Lift-slope approximations at 16.9 degrees . . 49

31 Lift and Drag coefficient curves, Re ≤7000 . . . 59

32 Lift and Drag coefficient curves, 7000≤Re≤15000 . . . 59

(11)

42 Lift and Drag coefficient curves, 220000≤Re≤1e6 . . . 64

43 Lift and Drag coefficient curves, 1e6≤Re≤1e9 . . . 65

List of Tables

1 Calculation times for collective pitch at 15.3 degrees with Intel Core i7-4720HQ 2.60GHz . . . 39

2 Calculation times for collective pitch at 15.3 degrees with Intel Core i5-6500T 2,5GHz . . . 39

3 Calculation times for collective pitch at 16.9 degrees with Intel Core i7-4720HQ 2.60GHz . . . 40

4 Calculation times for collective pitch at 16.9 degrees with Intel Core i5-6500T 2,5GHz . . . 40 5 Estimated Inflow Variations. Adapted from Table(3.1)[Leishman, 2006] 58 6 Thrust coefficient results comparison. Collective pitch at 15.3 degrees 65 7 Power coefficient results comparison. Collective pitch at 15.3 degrees 66 8 Thrust coefficient results comparison. Collective pitch at 16.9 degrees 66 9 Power coefficient results comparison. Collective pitch at 16.9 degrees 67

(12)

1 Introduction

Rotary wing aircrafts had a parallel development to fixed wing aircrafts at the start of the 1900s. However, after Ludwig Prandtl published his Lifting-Line Theory in 1918 developed with the help of the Vortex Theory, the calculations of the forces related to fixed wings started to become really accurate. Years later, Hermann Glauert formalized in 1935 the Momentum Theory developed by Rankine in 1865 and connected it to the Blade Element Theory suggested by Drzewiecki in 1892. The combination of these latter two theories gave the first accurate results for rotary wing aircrafts.

(a) Articulated rotor hub

(b) Rotor velocity distribution

Figure 1: Helicopter rotor characteristics. See Fig.1-4 in [Johnson, 1994] and Fig.2-1 in[Leishman, 2006] respectively.

The work of both Theodore Theodorsen and Sydney Goldstein gave an impulse to the theory related to rotating lifting bodies, although their work was more oriented to aircraft propellers. During the upcoming years, the application of the Vortex

(13)

Theory to helicopter rotors started to give good results; however, the mathematical models obtained depended on experimental observations beside complicated and long calculations involving complete elliptic integrals (see [Walter Castles, 1952]) and the analysis of the harmonics related to flapping (see [Robin B. Gray, 1954]). Figure(1a) shows the structure of an articulated rotor hub similar to the originally developed by Juan De la Cierva, whilst figure(1b) shows the velocity profiles for blades during hover and forward flight.

Helicopters are complicated to build due the difficulties and uncertainties that need to be faced. Modeling their main important component, being it the rotor, is the main task. Their capability for axial flight has been proven extremely useful for society and therefore improving our understanding of them is crucial. However, due to their complexity, it is always hard to find new ideas or design and moreover experimental data.

1.1 Background

Modeling rotary wing aircrafts nowadays is still far more complex than fixed wing aircrafts due to the flow properties. The tools available for these latter ones, such as the mathematical models developed or CFD analysis, can be rather simple for first estimations. The flow properties involved in rotary wing aircrafts are far more complex and hard to approximate since the flow approximations involve complex mathematics and several unknowns.

Some software tools available for rotary wing aircrafts are ANSYS-Fluent and ANSYS-CFX, oriented to flow simulations or CFD; the CAMRAD software devel-oped by NASA focused on blade-vortex iterations and rotor wake modeling (see [Johnson, 2012]); and CAMRAD II developed by Johnson Aeronautics with im-proved wake modeling and CFD capabilities.

However, a new approach was done in 2015 using Wolfram’s Modelica based SystemModeler by M. Dahl (see [Dahl, 2015]). A model developed from scratch was created in the way of blocks that can be connected to perform different simulations, obtaining not only numerical results but visual and realistic behavior to actually see how a rotor or helicopter will perform. The model created included complete helicopter aerodynamics, flight mechanics and control system. The Blade Element Momentum Theory was the chosen mathematical model for the calculations related to the rotor aerodynamics. Nonetheless, the lack of real testing data available hin-dered the validation of this outstanding model and the most certainly troublesome and important variable involved in rotary wing aerodynamics, the induced inflow.

1.2 Aim and Contribution

The aim of this present thesis will be the validation of the developed model by M.Dahl (see [Dahl, 2015]), based on Blade Element Theory and the capabilities of Wolfram’s SystemModeler software as a strong and proper tool for helicopter simulations and estimations. The research questions that this work will try to answer will be:

• Is the Blade Element Theory a suitable mathematical model compared to other theories?

(14)

• Is it possible to obtain a rotor modeling tool for both real sized and scaled rotors?

• How noticeable are compressibility effects?

• Is it possible to obtain good accuracy and simulation speed at the same time?

Contribution

The contributions to the model were:

• Relate Blade Element Theory and Vortex Theory. • Solve Lift - Slope uncertainties.

• Solid Stall and Drag criteria.

1.3 Literature Survey

With the development of Blade Element Theory and Vortex theories, several ap-proaches were done to obtain a general equation for the induced inflow of rotary wing aircrafts. This goal is, nevertheless, rather complicated due to the differ-ent flight properties of the rotor and the complexity of the mathematics involved. Therefore, the approximations recommended specifically for rotary wing aircrafts in Principles of Helicopter Aerodynamics ([Leishman, 2006]), Rotorcraft Aeromechanics ([Johnson, 2013]), Helicopter Theory ([Johnson, 1994]), The elements of Aerofoil and Airscrew Theory ([Glauert, 1983]) were used; whilst for general aerodynamics and vortex theory, Fundamentals of Aerody-namics ([John D. Anderson, 2007]) and Introduction to Vortex Theory ([Lugt, 1996]) were used respectively alongside many others.

As mentioned previously, experimental data was needed to validate the model. For this matter, the article Comparisons of Predicted and Measured Rotor Performance in Vertical Climb and Descent ([Felker and McKillip, 1994]) was used. The validity of this publication as a proper source was guaranteed by the NASA agency and the already mentioned literature which uses this article as a reliable source.

1.4 Delimitations

The main delimitation for the present work was the lack of experimental data. Sev-eral companies and universities were were contacted without any results whatsoever. Thus, more experimental data for different rotor configurations would be the perfect way to achieve a wide validation for the model, although the methodology followed should be valid, as it will be seen, for almost every rotor configuration.

Another limitation is, that the article used ([Felker and McKillip, 1994]), only has experimental data from a rotor hovering and in axial flight, which means that no data were available for forward flight. Nevertheless, as it will be mentioned, forward flight estimations are based in the calculations made for a hovering rotor hence the model will be a proper starting point for future development.

(15)

1.5 Outline

The present work outline is: • Chapter 1: Introduction.

• Chapter 2: Theory. All theories applied in the present work are described in this chapter.

• Chapter 3: Method. The assumptions and developed method of the present work is described here.

• Chapter 4: Results. Values obtained during the different simulations are shown.

• Chapter 5: Discussion. The results are evaluated here, with respect to their proximity to the measurements.

• Chapter 6: Conclusion. Answers to the questions formulated in the aim. • Chapter 7: Perspectives. Plausible future work suggestions are given here. • Appendix: Appendix containing useful equations and a complete list of

(16)

2 Theory

2.1 Aerodynamic Forces

Flow along bodies is hard to predict and involves several considerations and long and though calculations, which nowadays are done with the help of computers through CFD simulations. Nevertheless, the aerodynamic forces on bodies are the result of the of the pressure distribution and the shear stress distribution over the body surface. In this section a short demonstration will be carried out to demon-strate that ”the sources of the aerodynamic lift, drag, and moments on a body are the pressure (p) and shear stress (τ ) distributions integrated over the body”(from [John D. Anderson, 2007]).

Figure 2: Aerodynamic forces on a lifting body (see figures(1.16) and (1.18) in [John D. Anderson, 2007]). a)Resultant forces b)Pressure and shear stress distribu-tion

As it is seen in figure(2), the resultant force, R, of a body with chord c (distance in a straight line from the leading edge to the trailing edge, see figure(3)) and planform area S in a freestream, V∞,is defined by the sum of its components, being

them:

• Lift force, L, defined as the component of the resultant force perpendicular to the freestream. The equation for the lift force is:

L = 1 2ρ∞V

2

∞SCL (1)

• Drag force, D, defined as the component of the resultant force parallel to the freestream. The equation for drag is:

D = 1 2ρ∞V

2

∞SCD (2)

Normal and axial forces are related to lift and drag forces geometrically through the angle of attack α (as can be seen in figure(2a). The equations to express this relation can be easily written such as:

L = N cos α − A sin α (3)

D = N sin α + A cos α (4)

The resultant aerodynamic force can also be defined as the sum of the compo-nents parallel (axial force, A) and perpendicular (normal force, N ) to the chord line:

(17)

• Normal force, N, defined as the component of the resultant force perpendicular to the chord line. The integration criteria can be seen in figure(3) with the resultant equation for the normal force per unit span ((N’ )) being:

N0= − Z T E LE (pucos θ + τusin θ)dsu+ Z T E LE (plcos θ − τlsin θ)dsl (5)

• Axial force, A, defined as the component of the resultant force parallel to the chord line. The integration criteria can be seen in figure(3) with the resultant equation for the axial force per unit span (A’ ) being:

A0 = Z T E LE (−pusin θ + τucos θ)dsu+ Z T E LE (plsin θ − τlcos θ)dsl (6)

Figure 3: Integration criteria for body surfaces (see [John D. Anderson, 2007]), adapted from figure(1.17)

2.1.1 Aerodynamic Dimensionless Coefficients

Dimensionless coefficients that are being used throughout this thesis are:

• Lift coefficient for a complete three-dimensional body and per unit span are respectively: CL= L 1 2ρ∞V∞2S (7) cl = L0 1 2ρ∞V∞2c (8)

(18)

• Drag coefficient for a complete three-dimensional body and per unit span are respectively: CD = D 1 2ρ∞V∞2S (9) cd= D0 1 2ρ∞V∞2c (10) • Defining here pressure as ”the normal force per unit are exerted on a surface due to the time rate of change of momentum of the gas molecules impacting on that surface” (from [John D. Anderson, 2007], Chapter 1), the pressure coefficient is expressed as:

Cp =

p − p∞ 1 2ρ∞V∞2

(11) Pressure coefficient can be expressed in another useful form using Bernoulli’s equation as follows: p∞+ 1 2ρV 2 ∞= p + 1 2ρV 2 ₍₁₂₎ Cp= 1 − ( V V ∞) 2 ₍₁₃₎

• Defining here the shear stress (τ ) as ”the limiting form of the frictional force per unit area” (from [John D. Anderson, 2007], Chapter 1), the skin friction coefficient is expressed as:

cf =

τ

1 2ρ∞V∞2

(14) • Normal force coefficient for a complete three-dimensional body and per unit

span are respectively:

CN = N 1 2ρ∞V∞2S (15) cn= 1 c[ Z c 0 (Cp,l− Cp,u)dx + Z c 0 (cf,u dyu dx + cf,l dyl dx)dx] (16)

• Axial force coefficient for a complete three-dimensional body and per unit span are respectively: CA= A 1 2ρ∞V∞2S (17) ca= 1 c[ Z c 0 (Cp,u dyu dx − Cp,l dyl dx)dx + Z c 0 (cf,u+ cf,l)dx] (18)

Note that equations(5)-(6) are related to equations(16)-(18) by a change of vari-ables: dx = ds·cos θ; dy = -(ds·sin θ); S = c(1). For a more detailed explanation see [John D. Anderson, 2007].

Finally, these coefficients can be expressed through the same geometrical relation using the angle of attack α as was shown previously:

cl= cncos α − casin α (19)

cd= cnsin α + cacos α (20)

This set of equations relates lift and drag coefficients to the pressure coefficient, which is an important result for the upcoming sections.

(19)

2.2 Dimensionless similarity parameters

Reynolds and Mach numbers are probably the most well known parameters in aero-dynamics. They can be used for model scaling and therefore it is interesting to introduce what they represent.

2.2.1 Reynolds number

The Reynolds number represents a ”ratio between inertial and viscous forces” (from [Leishman, 2006], Chapter 7). It is expressed as:

Re = ρV∞c

µ ≡

Inertialf orce

V iscousf orce (21)

Being ρ the fluid density, V∞ the flow velocity, c the characteristic length and µ the

fluid viscosity.

Usually Reynolds numbers with values lower than 106 represent flow where vis-cous forces prevail over inertial forces. In flow with Reynolds number values higher than 106 the inertial forces prevail and the approximation of inviscid flow can be used. For a deeper explanation of flow behavior related to Reynolds number, see Chapter 14 in [Lugt, 1996].

2.2.2 Mach number

The Mach number represents a ”ratio of inertia forces in a fluid to forces resulting from compressibility” (from [Leishman, 2006], Chapter 7).

M∞=

V∞

a (22)

Being a the local speed of sound. A low ratio implies that the flow changes its velocity and pressure gradually in a way that compressibility effects are negligible. A high ratio implies that compressibility effects are noticeable and, if the critical Mach number is reached on a surface (this depends on numerous factors such as geometry or angle of attack), shock waves can start developing on the surface. The development of shock waves is usually quite noticeable at transonic local Mach num-bers (0.7<M<1.2) and when the free stream reaches unity (M∞ = 1), the sound

barrier is reached and singularities start to happen. For further information see [Leishman, 2006], Chapter 7; or [John D. Anderson, 2007], Chapter 7, 8, 9, 11 & 12. Mach number regimes are: Subsonic when M<1, Sonic when M=1, Supersonic when M>1 and Hypersonic when M> 5.

2.2.3 Scaling and similarity

Full scale models are expensive to build and extremely resource consuming, for these reasons flow similarity is widely used for scale testing. Through dimensional analysis it can be demonstrated that flows are dynamically similar (see Chapter 1 in [John D. Anderson, 2007]) when:

• When the solids that are subject to the flow are geometrically proportional. • Similarity parameters have equal value.

(20)

A further explanation and demonstration can be found in both [Leishman, 2006] and [John D. Anderson, 2007] but, as can be easily noticeable with equations(21) -(22), the usual way to achieve the conditions previously asserted is by an increment in the flow speed to achieve the same Reynolds number(Re), which implies that compressibility variations can become rather considerable and by no means negligible (M>0.1).

2.2.4 Introduction to Boundary Layers

In 1904 Ludwig Prandtl introduced a revolutionary concept that became a break-through in the study of fluid dynamics. He proposed a concept known as boundary layer which states that there is a small region close to the surfaces of bodies in a flow where the velocity profile of the flow rises from 0 at the surface of the body, to a velocity of 99% of the free stream. Inside this boundary layer the effects of friction are quite noticeable and have a strong effect on the flow. Boundary layers are classified as:

• Laminar: In which the flow behaves neatly, so the velocity profile can be approximated as parallel layers of fluid along the surface. Up to Re=170, see [Lugt, 1996].

• Transitional: The flow is hard to model and calculate since the it is changing from the laminar with smooth behavior and linear properties to the more disordered and last kind of boundary layer, known as turbulent. Increasing irregularities and three-dimensional effects become noticeable. The Reynolds number range is 170<Re<4.5·105, see [Lugt, 1996].

• Turbulent: The flow does not follow parallel path in opposition to the ordered laminar flow, mixing continuously and also giving a velocity profile larger and more similar to the free stream flow. A flow is turbulent for Reynolds numbers bigger than 4.5·105.

The equations for each kind of boundary layer are reductions of the Navier-Stokes equations (see A.1). All of them are extensively explained and can be found in [John D. Anderson, 2007], Chapters 17 to 19.

It will be noted here that the characteristics of boundary layers are related to lift and drag:

• In the case of laminar flow, the local skin friction coefficient for incompressible flow over a flat plate is:

cf =

0.664 √

Rex

(23) • In the case of turbulent flow (usually Reynolds numbers bigger than 5x105_),

the local skin friction coefficient for incompressible flow over a flat plate is: cf =

0.02296 Re0.139

x

(24) These expressions can be found in [John D. Anderson, 2007], Chapters 18 & 19. As can be seen, these latter equations, together with equations(16)(18)(19)(20) show that flow properties, lift and drag values are related to the Reynolds number and its meaning.

(21)

2.3 Introduction to Vortex Theory

The Vortex Theory helped to develop the knowledge of fluid dynamics and can explain accurately the behavior of fluids. This section will explain the basic concepts of it since they will become handy at some points to explain different approaches that will be done and might not be obvious at a first glance.

2.3.1 Vorticity

Fluids are usually moving in the three dimensions, this results in a resultant angular velocity of its elements. A proper demonstration of this resultant velocity can be found in [John D. Anderson, 2007], being it:

~ ω = ωx~i + ωy~j + ωz~k (25) ~ ω = 1 2[( ∂w ∂y − ∂v ∂z)~i + ( ∂u ∂z − ∂w ∂x)~j + ( ∂v ∂x − ∂u ∂y)~k] (26)

This allows the definition of vorticity, ξ, as: ~ ξ = 2~ω = (∂w ∂y − ∂v ∂z)~i + ( ∂u ∂z − ∂w ∂x)~j + ( ∂v ∂x − ∂u ∂y)~k (27)

Thus, vorticity is defined as the curl of the velocity (see Chapter 5 in [Lugt, 1996] and Chapter 2 in [John D. Anderson, 2007], eq (2.127) to eq.(2.129)):

~

ξ = ∇ × ~V (28)

2.3.2 Circulation

Circulation of a vector field (~V ) is defined as ”the line integral of the tangential ve-locity component taken round a closed curve C” (from [Glauert, 1983], [Lugt, 1996] and [John D. Anderson, 2007]). Notice here that albeit circulation is usually ex-plained with the example of flow round cylinders, it does not imply a circular path of the flow. Γ ≡ − I C ~ V · d~s (29)

The minus sign in equation(29) appears due to the clockwise sense taken by con-vention in aerodynamics in opposition to the counterclockwise sense that is used by convention for line integrals.

Stokes’ theorem relates the line integral of ~V over C to the surface integral of ~V over S and thus vorticity (ξ) to circulation (Γ) can be related through equation(30):

Γ ≡ − I C ~ V · d~s = − Z Z S (∇ × ~V ) · d ~S (30)

Let the definition of streamline be a ”curve whose tangent at any point is in the direction of the velocity vector at that point” (from Chapter 2, point 2.11 in[John D. Anderson, 2007]), or in equation form:

(22)

If all the streamlines of a flow were concentric about a point so the radial velocity is zero (Vr) and the tangential velocity Vθ = constant/r, it is possible to demonstrate

that, for any streamline at a distance r from the center, the circulation is: Γ ≡ − I C ~ V · d~s = Constant r = −Vθ(2πr) (32) and thus: Constant = − Γ 2π (33)

Circulation, Γ, is usually called strength of a vortex flow. By convention, positive strength is defined in the clockwise direction. Notice the singularity at r=0 in the flow field. See Chapter 3 in [John D. Anderson, 2007] and Chapters 5&6 in [Lugt, 1996]. As a historical note, this same kind of procedure was proposed by R.P. Feynman to account quantized vortices in superfluid helium.

2.4 Two-Dimensional Thin Airfoil theory

The equations and definitions included so far give the tools for the theory that tries to approximate the lift generated by geometries in a flow stream. Historically this theory helped to develop the design and construction of sophisticated wings and lifting surfaces.

2.4.1 Flow around a rotating cylinder

It is needed to define first the velocity field for a nonlifting flows over circular sta-tionary cylinders. The demonstration is easy to follow but still quite long, so it is not included here but can be found in Chapter 3 in [John D. Anderson, 2007] and Chapter 6 in [Lugt, 1996]. Radial and tangential velocities are, respectively, in this case: Vr= (1 − R2 r2)V∞cos θ (34) Vθ = −(1 + R2 r2)V∞sin θ (35)

Equaling these latter equations to zero the stagnation points can be found.

If the case of a cylinder spinning about its axis is considered now, the flow gets more complex, but the tools to describe what happens are already available. The nonlifting flow mentioned previously can be combined (due to the flow linearity) with the circulation introduced in section 2.3.2 and thus equations (34),(35),(32) give the resultant velocity field at a point of distance r from the center:

Vr= (1 − R2 r2)V∞cos θ (36) Vθ = −(1 + R2 r2 )V∞sin θ − Γ 2πr (37)

If the previous definition of pressure coefficient in equation(13) is used combined with equation (37) for points located on the cylinder surface (r = R), the pressure

(23)

coefficient is then related to both the free stream velocity (V∞) and the circulation (Γ) as follows: Cp = 1 − ( V V∞ )2 = 1 − (−2 sin θ − Γ 2πRV∞ )2 = 1 − [4 sin2θ + 2Γ sin θ πRV∞ + ( Γ 2πRV∞ )2] (38)

See Chapter 3 in [John D. Anderson, 2007] for a more detailed demonstration and full explanation.

If the usual assumptions to model flows are done now (at high Reynolds numbers flows are usually considered inviscid and incompressible), considering cf=0, equation

16 becomes simplified and related to lift as follows: cn= 1 c Z c 0 (Cp,l− Cp,u)dx = cl (39)

If this latter equation is converted to polar coordinates (x=R·cosθ; dx=-R·sinθ·dθ), with the relation between radius and chord being R=c/2 and definite integral prop-erties, and that eq(38) defines the pressure coefficient for both lower (Cp,l) and upper

(Cp,u) sides,it becomes:

cl = 1 c Z c 0 (Cp,l− Cp,u)dx = 1 c Z c 0 Cp,ldx − 1 c Z c 0 Cp,udx = −R c Z 2π π Cp,lsin θdθ + R c Z 0 π Cp,usin θdθ = − c 2c Z 2π π Cp,lsin θdθ − c 2c Z π 0 Cp,usin θdθ = −1 2( Z 2π π Cp,lsin θdθ + Z π 0 Cp,usin θdθ) = −1 2 Z 2π 0 Cpsin θdθ (40)

Introducing equation(38) into equation(40) the result is: cl=

1 2

Z 2π

0

(− sin θ + 4 sin3θ +2Γ sin

2_θ

πRV∞

+ ( Γ

2πRV∞

)2sin θ)dθ (41) and after integration:

cl= 1 2(cos θ] 2π 0 +4( −3 cos θ 4 + cos 3θ 12 )] 2π 0 + 2Γ πRV∞ (θ 2− sin 2θ 4 )] 2π 0 −( Γ 2πRV∞ )2cos θ]2π₀ ) (42) or simply: cl= Γ RV∞ (43)

(24)

This latter equation introduced in the equation for lift per unit span (denoted as L’ ) in the case of a cylinder [S=2R(1)] gives the resultant:

L0 = 1 2ρ∞V 2 ∞Scl= 1 2ρ∞V 2 ∞2R Γ RV∞ (44) or in the simpler form in which is well known as the Kutta-Joukowski theorem :

L0 = ρ∞V∞Γ (45)

This theorem states that circulation is directly related to the lift force and is one of the basis for modern aviation. Note here that this demonstration implies that the drag for a rotating cylinder is zero due to the inviscid flow assumption.

2.4.2 Two-Dimensional Wing Theory

This section will relate the Kutta-Joukowski theorem to the lift generated by thin wings. The demonstration, albeit short, requires knowledge of the Vortex Theory and several theorems related to it, thus we refer to Chapter 6 in [Lugt, 1996].

As was previously shown, circulation generates lift. Lanchester-Prandtl used this relation to create a model using the concept of bound vortex. Figure (4) shows an approximation of circulation due to a flow for a wing airfoil.

Figure 4: Circulation for a wing profile in a flow

If the airfoil profile is approximated to a flat plate of width c at an angle of attack α, after expressing its coordinates in the complex plane, expressing the flat plate as a circle with a radius r=c/4 and several mathematical operations and the use of Blasiu’s Theorem, the Kutta Joukowski equation for wings is obtained:

X − iY = iρ∞ΓV∞e−iα (46)

Being the real part, X, the drag generated by the flat plate while the imaginary part would be the lift force L’=ρ∞ΓV∞. This simple equation states that the drag is zero

and the lift is completely independent of the shape. Moreover, equation(46) at zero angle of attack gives exactly the same result that the cylinder case, equation(45). See Chapter 6, point 6.4, and equations(6.51) to (6.60) in [Lugt, 1996].

If the relation stated previously for a flat plate of width c and a cylinder of radius r=c/4 is kept, the circulation is:

(25)

Thus, with the previous relation stated for lift, introducing the latter equation in it:

L0 = πρ∞V∞2c sin α (48)

Using here the equation8 to define the lift coefficient:

cl=

πρ∞V∞2c sin α 1

2ρ∞V∞2c

(49)

Resulting in the well know approximation for small angles of attack for the lift coef-ficient of 2πsinα, or 2πα. It has been demonstrated that this approximation follows experimental data accurately at high Reynolds numbers (Re>106) and small angles of attack (see figure(14.22) in [Lugt, 1996] or figure(4.25) in [John D. Anderson, 2007]).

2.4.3 Airfoil Lift Slopes

As demonstrated in the previous point, the lift coefficients of lifting surfaces or thin airfoils change at a certain rate, following the definition of gradient or, as they are usually called, a lift-slope:

clα =

dcl

dα (50)

As mentioned previously, it can be approximated to a linear lift-slope of 2π. However this has been stated under the assumption of inviscid flow, which implies a high Reynolds number and it still is unrealistic. These lift slopes are characteristic for each airfoil profile, vary depending of the characteristics of the flow (laminar or turbulent) and are obtained through experimental data or approximated by refined computer simulations. In the Appendix(A.1) several curves of lift can be seen; also in [Abbott and Doenhoff, 1976] a extensive explanation of airfoil design, approximation and experimental lift-slopes can be found.

2.4.4 Stall

An airfoil is considered to be stalled when, due to viscous forces, the flow gets sep-arated from the lifting surface leading then to a sudden decrease in the resultant lift force. This effect can be observed in the characteristic lift curves of airfoils and always happens after the angle of attack where the lift coefficient reaches its maxi-mum. Beyond the point of maximum lift coefficient, further increasing of the angle of attack makes the flow incapable of following the surface and a reverse flow starts to happen. Viscous forces prevail over inertial forces when this reverse flow appears. A sudden increase in drag happens once stall begins. See [John D. Anderson, 2007]. The growth rate, defined using equation(50), does not always follow this linear growth rate. The characteristic lift-slopes for each airfoil differs until the Reynolds number reaches values over 1x106, when the growth becomes linear, as has been experimentally demonstrated. See figures from (31) to (43) in Appendix(A.1); also [Abbott and Doenhoff, 1976] and [Pope, 2009] can be consulted for further informa-tion and experimental lift-coefficient curves.

(26)

2.4.5 Drag

To describe drag, equation(20) needs to be referred here. In order to obtain the mathematical models for thin airfoils described in previous points this chapter, the flow was assumed inviscid and therefore the friction coefficient neglected. This as-sumption leads to a drag coefficient of zero, whether the cylinder or flat plate ap-proaches are followed (see sections 2.4.1 and 2.4.2; also [John D. Anderson, 2007] and [Lugt, 1996]).

This result is obviously unreal, and it has been demonstrated that, even though at an angle of attack where a profile has a resultant lift coefficient of zero (usually at zero degrees angle of attack), the drag coefficient is, indeed, not zero (and usually denoted as cd0). This means the viscosity is present and leads to a sear stress on

the surface. See equations(20) and (18).

Also, as mentioned previously, when a profile reaches a certain angle and starts stalling, a sudden increase in drag happens. This is due to the separated flow generating what it is called pressure drag (see Chapter 4 in [John D. Anderson, 2007] for further information).

2.4.6 Compressibility corrections

The usual assumption of incompressible flow is applicable or rather true below Mach numbers below M∞ < 0.1. For values beyond that point, the

assump-tion of incompressibility becomes weak and hence the results lose accuracy. Lud-wig Prandtl in 1922 and Hermann Glauert in 1927 started using a compressibil-ity correction that today is known as the Prandtl-Glauer rule (see Chapter 11 in [John D. Anderson, 2007]): Cp = Cp,0 p|1 − M2 ∞| (51) Equation(51) relates compressible and incompressible pressure distribution over an airfoil. If the incompressible distribution is known, it can be corrected through this latter equation and obtain the compressible distribution directly. The whole demonstration of this equation can be found in the original publication by H.Glauert (see [Glauert, 1927]) or in Chapter 11 in [John D. Anderson, 2007].

If the previous equations(16),(18) and (19) are revised at this point, it can be seen that as a matter of fact, the same compressibility correction can be used for the lift coefficient since it is directly related to the pressure coefficient (see Chapter 11 in [John D. Anderson, 2007], thus:

cl=

cl,0

p|1 − M2 ∞|

(52)

Equations(53) and (54) are more recent compressibility corrections. Being them, respectively, the Karman-Tsien rule and the Laitone’s rule:

Cp = Cp,0 p|1 − M2 ∞| + [M∞2/(1 +p|1 − M∞2 |)]Cp,0/2 (53) Cp = Cp,0 p|1 − M2 ∞| + (M∞2 [1 + [(γ − 1)/2]M∞2 ]/2p|1 − M∞2|)Cp,0 (54)

(27)

Compressibility corrections, nonetheless, present some issues. In figure(5) can be seen that when M∞= 1 the Prandtl-Glauert rule has a singularity and the Laitone’s

rule tend to zero. Moreover, albeit they are valid for subsonic and supersonic flows (up to M∞= 5 approximately), they do not give an accurate correction for transonic

flow since it is highly nonlinear (see [Vos and Farokhi, 2015] and Chapter 11 in [John D. Anderson, 2007]). 1 2 3 4 5 0.0 0.5 1.0 1.5 2.0 Cp0 1-M 2_ Cp0 1-M 2_ + M2Cp0 1 1M 2_{ 2} Cp0 1-M 2_ + Cp0M2_11 2(1.41) M 2_ 2 1M 2_

Figure 5: Compressibility corrections comparison

2.5 Momentum theory

Flow through a helicopter rotor is quite complex to describe, however, they way to approach was using a surface as a control volume that wraps the rotor, its wake downstream and the upstream flow. With the assumptions of quasi-steadiness, one-dimensional incompressible and inviscid flow, the following set of equations can be used:

• A reduced version of the continuity equation, implying that the mass flow coming into the control volume is equal to the mas flow leaving the control

volume: _{Z Z}

S

ρ~V · d ~S (55)

• The governing equation of fluid momentum, implying that the net force on the fluid equals the fluid momentum change rate:

~ F = Z Z S pd ~S + Z Z S (ρ~V · d ~S)~V (56)

• The governing equation of conservation of energy, implying that the increase on the kinetic energy of the fluid depends on the work done by the rotor:

W = Z Z

S

1

2(ρ~V · d ~S)~V (57)

These assumptions and equations are written following [Leishman, 2006] and represent the basis of the Momentum Theory for rotary wing aircrafts. For the his-torical development of the Momentum Theory see [Glauert, 1983]. For the complete

(28)

Continuity equation, Momentum equation and Energy equation development, see [John D. Anderson, 2007]. The complete form of the equations can also be found in A.1

It is possible to continue henceforth with the assumptions and the equations obtained to obtained the different values related to the rotor.

The mass flow rate through the control volume using equation(55) is: ˙ m = Z Z ∞ ρ~V · d ~S = Z Z 2 ρ~V · d ~S (58) or simply ˙ m = ρA∞w = ρA2vi (59)

Being w the flow velocity at the far wake of the rotor and vi the induced velocity

at the rotor disk. This latter value is probably the most important factor for rotor design and specially hard to predict.

If now equation(56) is used to analyze a helicopter hovering, then it becomes: − ~F = T = Z Z ∞ ρ(~V · d ~S)~V − Z Z 0 ρ(~V · d ~S)~V (60)

or, due to the null velocity far upstream, T =

Z Z

∞

ρ(~V · d ~S)~V = ˙mw (61)

Being T the thrust produced by the rotor, equation (61) gives the relationship between thrust, mass flow and flow velocity.

Now, using (57) for the work produced by the rotor, it is possible to find an additional equation: T vi= Z Z ∞ 1 2(ρ~V · d ~S)~V 2₋ Z Z 0 1 2(ρ~V · d ~S)~V 2 ₍₆₂₎

Being the second term of the right-hand side of the equation zero for a hovering rotor, the reduced expression becomes:

T vi = Z Z ∞ 1 2(ρ~V · d ~S)~V 2 ₌ 1 2mw˙ 2 ₍₆₃₎ or simply vi = 1 2w (64)

Equation (64) represents an extremely important relation between the flow velocity at the far wake and the induced velocity at the rotor disk. Several approaches have been done to calculate both of them throughout the history and together with equation (61) are the basis of the Momentum Theory for rotary wing aircrafts.

(29)

Now it is possible to obtain values for a hovering helicopter such as the thrust (T ):

T = ˙mw = ˙m(2vi) = 2(ρAvi)vi = 2ρAv2i (65)

or the power required:

P = T vi = T

s T

2Aρ (66)

From now on, the induced velocity at the rotor, vi, will be differenced from the

induced velocity for a hovering rotor vh.

2.5.1 Useful Dimensionless Coefficients

Usually the different variables are nondimensonalized using the blade tip speed of the rotor, Vtip, thus, the dimensional analysis gives the following dimensionless

coefficients (see [Leishman, 2006]):

• Thrust coefficient is nondimensionalized: CT = T ρAV2 tip = T ρAΩ2_R2 (67)

Being Ω the rotational velocity of the rotor and R the rotor radius.

• The induced velocity has a nondimensionalized form as well and it is known as the induced inflow ratio, λ, being for a hovering rotor:

λh = vh ΩR = 1 ΩR s T 2ρA = r CT 2 (68)

Being the general case for all flight modes, with Vc the climbing speed of the

helicopter:

λ = vi+ Vc

ΩR (69)

It is obvious that this nondimensionalization assumes that the induced velocity is distributed equally all over the rotor disk.

• Power coefficient is defined then such as: CP = P ρAV3 tip = P ρAΩ3_R3 (70)

or according to equations(66) and (68)

CP = T vi ρAΩ3_R3 = C_T3/2 √ 2 = (71)

These coefficients (CT and CP) will be used throughout this thesis and follow the

US definition. They can be nondimensionalized using a different approach (usually Europe, Britain and Russia), resulting in:

CT =

T

1

2ρAΩ2R2

(30)

CP =

P

1

2ρAΩ3R3

(73) This section 2.5 follows the thorough theoretical analysis and equations that can be found in the chapter 2 in [Leishman, 2006], a longer description can be found there.

2.6 Blade element theory

The Blade Element Theory was developed in an attempt to obtain thrust gen-erated by propellers and helicopter rotors. It has, however, some important issues.

Figure 6: Blade element velocity components (see [Leishman, 2006], adapted from figure(3.1)

.

The idea of the Blade Element Theory is to divide each blade in small elements and calculate the lift generated by each element in order to obtain the lift force profile generated by the blade. Now, if the same nomenclature is used as in figure(6), some parameters can be defined (see [Leishman, 2006]):

(31)

• UP, as the local velocity perpendicular to the rotor at the blade element.

• U, as the resultant velocity at the blade element. • α, as the angle of attack at the blade element • θ, pitch angle at the blade element

• φ, relative inflow angle or induced angle of attack at the blade element Now some geometrical relations can be established as can be seen in figure(6):

U = q U_T2 + U_P2 (74) φ = tan−1(UP UT ) (75) α = θ − φ (76) UP = Vc+ vi (77)

With these equations and equation(69) the main issue can be noticed. Through some calculation, these values can be obtained but one, which is the main issue for all the theories related to the calculation of the thrust generated by a rotor. This value is the induced flow, vi. This value influences the resultant angle of attack,

hence the lift coefficient of the blade element and thus the lift generated and so on. The problem is that Momentum Theory and Blade Element Theory can only offer rough and completely linear approximations, which represent a limitation for rotor design.

To solve this problem, a combination of the Blade Element Theory and the Momentum Theory was proposed (see [Knight and Hefner, 1937], [Glauert, 1983] and [Leishman, 2006]).

To define the equations used by the combined Blade Element-Momentum The-ory, or BEMT, the concept of solidity (σ) needs to be defined. It represents, as H.Glauert defined it, the ratio of the rotor blade area to the rotor disk area:

σ = Nbc

πR (78)

Being Nb the number of blades, c the blade chord, the position y as in figure(7)

and R the rotor radius. Now, if we proceed in the same way as in section(2.5) for Momentum Theory :

(32)

Figure 7: Rotor disk. Adapted from figure(3.4) - [Leishman, 2006]

If the latter equation is expressed as a nondimensional parameter: dCT = dT ρ(πR2_)(ΩR)2 = 2ρ(Vc+ vi)vidA ρ(πR2_)(ΩR)2 = 4πρ(Vc+ vi)viydy ρ(πR2_)(ΩR)2 (80)

If the definition of inflow ratio, λ, is used now, this latter equation can be expressed simply as:

dCT = 4λλirdr (81)

Using the relation between power and thrust (see eq.66), the induced power coeffi-cient is:

dCPi = λdCT = 4λ

2_λ

irdr (82)

If now the definitions of lift and drag per unit span are used: dL = 1 2ρcU 2_c ldy; dD = 1 2ρcU 2_c ddy (83)

As can be seen in figure(6), the normal and axial forces are geometricaly related again to the lift and drag forces, in this case through φ:

dFz= dL cos φ − dD sin φ; dFx = dL sin φ + dD cos φ (84)

Thus:

dT = NbdFz; dQ = NbdFxy; dP = NbdFxΩy (85)

Now, since UP is way smaller than UT, it can be assumed that U≈ UT and small

angle approximations (φ ≈ UP/UT), thus:

dT = NbdL (86)

dQ = Nb(φdL + dD)y (87)

dP = Nb(φdL + dD)Ωy (88)

With these results, we can obtain new expressions for the nondimensional parame-ters: λ = Vc+ vi ΩR = Vc+ vi Ωy ( Ωy ΩR) = UP UT (y R) = φr (89)

(33)

dCT = 1 2σclr 2_dr ₍₉₀₎ dCP = 1 2σ(φcl+ cd)r 3_dr ₍₉₁₎

If now equation(50) is used, the lift coefficient can be related to the local lift coefficient through the lift-slope and the different angles:

cl = clα(θ − α − φ) (92)

Combining α and θ in this latter expression (α ≈ 0 since this analysis is suitable for hover and axial flight), equation(90) and (89) can be written now:

dCT =

σclα

2 (θr

2_{− λr)dr} ₍₉₃₎

If equations(93) and (81) are combined, an expression for the inflow ratio can be obtained: λ(r, λc) = r (σclα 16 − λc 2 ) 2₊σclα 8 θr − ( σclα 16 − λc 2) (94)

This latter equation is the basis of the BEMT and represents a powerful tool for helicopter rotors development. It is easy to include a twist ratio for the blades and gives the approximated inflow due to circulation on each element of the blade. 2.6.1 Tip Loss

Approximately in 1919, Ludwig Prandtl calculated the approximated losses at the tip. Using Vortex Theory, he approximated the lift reduction to a factor (Prandtl’s F function):

F = 2 πcos

−1_{e(r − 1)N/2λ} ₍₉₅₎

This factor is related to the blade circulation (Γ) and it is well explained and demon-strated in Chapter 3,[Johnson, 2013]. To calculated this factor F is necessary to iterate with the inflow ratio λ, however, it is usually approximated to easier expres-sions, such as the one proposed by Gessow & Myers in 1952 based on the blades geometry and that was developed empirically:

B = 1 − c

2R (96)

Several other approximations are described in [Leishman, 2006] and [Johnson, 2013].

2.7 Forward flight

Forward flight is approximated using Momentum Theory and was proposed by H.Glauert (see [Glauert, 1983] and [Leishman, 2006]). However, Glauert himself commented for the proposal that the assumption is weak although the mathemat-ical model described all the flight modes and was accurate for hover, climb and descent flight. The assumptions made by Glauert can be seen in figure(8)

(34)

Figure 8: Forward flight model. Adapted from Figure(2.23) in [Leishman, 2006]

If the concept of advance ratio, µ is defined here as: µ = V∞cos α

ΩR (97)

The advance ratio can be related to the inflow ratio with the expression: λ = V∞sin α

ΩR +

vi

ΩR = µ tan α + λi (98)

With λi being in this case:

λi =

λ2_h p

µ2_{+ λ}2 (99)

To obtain this inflow, numerical iterations are needed. 2.7.1 General Inflow Equation

The general equation that describes all the flight modes for helicopters using the Momentum Theory is (see [Leishman, 2006]:

λ = (V∞cos α ΩR ) tan α + λ2_h q (V∞cos α ΩR )2+ λ2 +V∞sin α ΩR (100)

Due to the complexity of the inflow and its dependence to velocity, numerous cor-rections were proposed for the inflow model. See table (5) in A.3.

2.7.2 Flapping

For a complete model of the rotor using the Blade Element Theory, it needs to be mentioned the flapping effect together with the forward flight since it is when it becomes more noticeable. If we define:

(35)

• β as the cone angle, due to the lift forces that make the blade to displace upwards.

• y ˙β as the flapping rate induced velocity.

• ψ as the azimuthal angle that is used to locate the position of the blade (see figure(1b)

If the flapping effect is included to the model, the velocities described for the Blade Element Theory become:

UT(y, ψ) = Ωy + V∞sin ψ = Ωy + µΩR sin ψ (101)

UP(y, ψ) = (λc+ λi)ΩR + y ˙β(ψ) + µΩRβ(ψ) cos ψ (102)

(36)

3 Method

In this section the different decisions taken will be justified with the theoretical definitions described in section 2 as well as the different assumptions and inputs. The differences to the previous approach (see [Dahl, 2015]) will be explained and the suitability of the different approach for a faster modular based Modelica library. Moreover this section will help to motivate the election of this Modelica based ap-proach over the more resource consuming CFD simulations. The outline of the method will be:

• Boundary Conditions

• SystemModeler Motivation and Implementation • Modeling with Reynolds and Mach numbers • NACA0015 Airfoil properties

3.1 Experiment Setup Geometry

The way to validate a model is to compare with real test data, so the difference be-tween the predictions and reality can be evaluated. For the model developed in this thesis, based on the Blade Element Momentum Theory, data obtained from hovering rotors, and in axial flight was needed since these cases are well studied, with the the-oretical equations properly developed and serve as the results used to approximate forward flight values, validate the rotor design, blade configuration, airfoil choice, etc. For this purpose, a NASA report was used (see [Felker and McKillip, 1994]) where an experiment was carried out comparing several measurements with differ-ent theoretical approaches. In this report, all the necessary information is given, so the experiment can be reproduced using SystemModeler. In further sections the results will be shown and compared, whereas the suitability of the model will be deeply discussed.

The parameters from the report used are: • 4-bladed rotor • Radius = 1.22[m] • NACA0015 airfoil • Root Cutout = 10% • Blade chord = 0.0635 • Solidity (σ) = 0.0663

• Blade twist rate = -8 degrees/radius • Tip speed = 55[m/s]

(37)

3.2 Mathematical Model

In the theory chapter, the Blade Element Momentum Theory was thoroughly de-scribed, whilst the Vortex Theory is presented as a short introduction. Nevertheless, the Vortex Theory is described as the basis of several theories related to lift calcu-lation. Equation (94), obtained using the BEMT approach and shown once again below, was also obtained by M. Knight and R.A. Hefner in 1937 using theVortex Theory by means of Prandtl’s flow potential at a point P due to vortex elements (see [Knight and Hefner, 1937] and [Prandtl, 1923] for the complete explanation).

λ(r, λc) = r (σclα 16 − λc 2 ) 2₊σclα 8 θr − ( σclα 16 − λc 2)

Thus, the equation(94), obtained by means of two different theories, can hence-forth be used as a mathematical approximation for the rotor inflow ratio and as the basis for rotor modeling estimates. To avoid negative values of the square root for negative angles of attack, which are likely to happen with twisted blades, absolute values are taken in one of the terms:

|σclα

8 θr|

The square root is multiplied then +1 or -1 depending on the resultant angle of the sum of the collective pitch angle plus the resultant angle of the twist ratio at the blade element:

Collective − P itch − Angle + T wist − Ratio ∗ r |Collective − P itch − Angle + T wist − Ratio ∗ r|

This small modification removes all the problems related to negative values inside the square root and also provide the model with the capability to simulate every angle.

3.3 SystemModeler Motivation and Implementation

SystemModeler is a Modelica based software, whose properties allow the user to de-velop libraries capable of simulating real systems through the inclusion of predefined or user-defined components. These will interact together and approach the model as close as possible to the way the system will behave in reality. This capability allow the users to create and modify complete libraries to simulate several different situations, from easy and simple approaches to complex and complete situations. In the present case, SystemModeler features allow the user to simulate how a simple rotating rotor will behave, which is the present case, and obtain important values such as the resultant thrust or the inflow. However, it is possible to go even further, a complete model can be created (see [Dahl, 2015]) and simulate a complete heli-copter. From a situation starting on the ground, axial flight to gain altitude, hover, forward flight or descent. Moreover, it can be added a complete control system, complete flight mechanics models and so on.

Other tools available mentioned before provide more specific results that can be used to confirm the simulations or find out possible issues or unexpected interactions, for instance, the CAMRAD software mentioned before that is used to confirm the

(38)

wake of the rotor and the lifting surface. A special mention for the different CFD analysis needs to be done here. It is well known that they give accurate results, show the flow behavior and resultant forces. Nonetheless, the computational costs and the complexity of these simulations, especially for rotary wing aircrafts is extremely big and should be used as a verification of the designs, not as a testing tool. 3.3.1 Previous Work

A short description of the approaches followed by M. Dahl will be described here shortly. For a complete description see [Dahl, 2015]. As can be seen in figure (9), the block consist in four blades that receive the airfoil profile characteristics, a revolute that makes then rotate at a certain velocity about the z-axis, an input for the wind and finally the air characteristics depending on the altitude. Each blade is divided in smaller blade elements and the resultant lift and drag forces are calculated for each one of them together with the inflow. The approach included flapping of the blades as well, allowing them to move three dimensionally but this will not be described here. The equations used for the aerodynamic forces were (see equations 3.165 and

(a) Articulated rotor hub

(b) Rotor velocity distribution

Figure 9: a)Modelica based block b)SystemModeler’s resultant animation

(39)

• Lift force: dL = 1 2ρ · Cl· c · U 2_{· dy} • Drag force: dD = 1 2ρ · Cd· c · U 2_{· dy}

• The induced velocity was calculated using momentum theory. the variable Trev that expresses the time for a blade to rotate to the position were the

previous blade was. This approach was similar to [Knight and Hefner, 1937]. The equation used was:

vi=

−lρTrevUTVc±

q

2lρTrevUT + l2ρ2Trev2 UT2Vc2

2lρTrevUT

This approach, although simple and rather easy to implement had the identified downsides of the Momentum theory and did no include the angle of attack variation due to the induced velocity, therefore it needed to be improved. Another identified issue was a delay caused by the flapping calculations that made the model rather heavy in therms of computational costs.

3.3.2 Early Approach

A class or block in the simplest form possible, being it, the set of equations and parameters needed to obtain the power coefficient (CP), the thrust coefficient (CT)

and the inflow ratio for hover and axial flight, being respectively the equations (91), (93) and (94). dCP = 1 2σ(φcl+ cd)r 3_dr dCT = σclα 2 (θr 2_{− λr)dr} λ(r, λc) = r (σclα 16 − λc 2 ) 2₊σclα 8 θr − ( σclα 16 − λc 2)

The input parameters from the user for this case were: radius of the rotor, rota-tional velocity of the rotor, number of blade elements, blade chord, the root cutout, solidity, twist rate, collective pitch angle. The block included a CombiTable (see [Wolfram, 2017]) component which interpolates the properties of the NACA0015 (lift-slopes, stall angles, drag approximations) depending on the Reynolds number. These properties are loaded into the CombiTable by a user defined file (a simple *.txt file), which means that any kind of airfoil profile can be included as long as all its properties are defined.

(40)

3.3.3 Computer Controlled Model

Once the first and simple approach was concluded and tested, a more elaborated con-figuration was created. The previous block was reused, only that the rotational ve-locity and the collective pitch were removed as user-defined parameters and changed to block inputs. This model basically consists in a rotating cylinder (bodyCylinder1 in figure(10)) that represents the rotor with a weight equal to the thrust gener-ated when hovering for the angles in [Felker and McKillip, 1994]. A simple control system was created: a sensor measuring its absolute position (absolutePosition in figure(10)) was fed to a parameter controller, called LimPID (see [Wolfram, 2017] for a complete description), capable of changing the pitch angle to obtain enough thrust to make the disk climb up to a certain height and keep it hovering.

Figure 10: Climbing disk with simple control system

Figure 11: Tip path

Figure (11) shows the tip path from the disk. The rod on the disk represents the azimuth angle and as can be seen it starts at 0 degrees (aligned with x-axis). Although similar to vortex wakes for ideal cases, these path lines follow only the tip trajectory. For proper vortex wakes and their real behavior see, for example, [Lugt, 1996] (Chapter 4, ”Decay of a Rankine vortex”).

(41)

3.3.4 Complete model

Again, using SystemModeler’s scalability, the original library was modified again so the system will feed the collective pitch angle, the rotational velocity and a tilt angle (or angle angle of attack of the whole rotor, being zero the angle at which the rotor is parallel to the ground). It also has the general inflow equation(100) implemented together with the BEMT, so the results obtained with this latter theory can be used to iterate and obtain the forces from each blade separately. The way this was achieved was using the linear inflow assumption using equation (102):

UP(y, ψ) = (λc+ λi)ΩR + y ˙β(ψ) + µΩRβ(ψ) cos ψ

The right hand side of this latter equation includes the flapping terms based in β and ˙β . They are both included but, unfortunately, no available data was found and hence remained unverified.

The class was used to obtain the results of axial flight and hover from the blades with the uniform inflow assumption (similarly to the CAMRAD software); the disk was only allowed to move parallel to the ground this time, since it was also used to observe the behavior of the model in forward flight, although the lack of available experimental data for the forward flight mode stopped the possibility of further comparisons.

Figure 12: General model library

Figure (12) shows how the complete model can be achieved with a rather sim-ple block, capable of simulate all the different flying modes. It can also be easily manipulated to connect it to other blocks, create a more complex model, etc.

(42)

3.4 Scaling with Reynolds and Mach numbers

The Reynolds number is widely used as a similarity parameter between scaled mod-els. This property is used in the present model so any kind of geometry can be defined and SystemModeler will load the airfoil properties from a file defined by the user. For the present model, a range of Reynolds number from Re=2000 to Re=109 was included in the file (the complete list of values can be seen in A.4). SystemModeler will extrapolate the results by itself hence the more detailed the file, the more realistic the model will be.

On the other hand, the Reynolds number similarity has a serious drawback: to achieve the same Reynolds number with smaller geometries, velocity needs to be increased. The usually negligible compressibility effects become quite noticeable. For this matters, the Mach number becomes quite relevant and compressibility cor-rections that depend on it need to be done.

There are several compressibility corrections, being widely known the Prandtl-Glauert rule as well as the more refined ones by Karman-Tsien rule and Laitone’s rule. The chosen one for the present model is the Karman-Tsien rule, equation(94), not only for its proximity to experimental results that ensure values corrected with good accuracy but also for its obvious avoidance of singularities at the value of M=1.

Cp =

Cp,0

p|1 − M2

∞| + [M∞2/(1 +p|1 − M∞2 |)]Cp,0/2

To demonstrate the validity of having a compressibility correction included in the model, it was tested out with both corrected and uncorrected properties due to compressibility and compared with the available experimental data as it will be seen in upcoming sections.

3.5 NACA0015 Airfoil properties

Each airfoil has its own properties due to its geometry. The NACA agency has worked on the development of profiles and the evaluation of their properties for years. One of the most well known profiles is the NACA0015 profile. Its properties are widely known and have been studied for years. The complete collection of lift coefficient curves used has been included in A.4. It must be mentioned here that the first helicopter blades were developed using the NACA0012 profile being it still the most common profile for rotors.

3.5.1 Lift Slope

The lift slope at small angles of attack can be approximated accurately to 2π per radian when the Reynolds number is sufficiently high. To develop a general model, this approximation does not suffice, since it does not consider the whole range of flow behaviors such as laminar, transient or turbulent. Using the computed lift coefficient curves (see A.4), and with figure 13 as an example, the different lift slopes where fed to SystemModeler in this way:

(43)

Reynolds=43000 -10 -5 5 10 AoA_° -0.5 0.5 Lift_{- Coefficient}

Figure 13: NACA0015 Lift coefficient curve for Reynolds=43000

• Depending on the flow properties, the lift curve can change the slope, therefore, the first step was to check whether there is an angle at which the slope changes. In the figure (13) example, it was approximated to 3 degrees, considering its starting point 0 since the instability from, approximately 0 to 1 degrees can be neglected without making a big difference.

• The second step was to calculate the lift slopes of the obtained curves at each Reynolds number. As can be seen in A.4 and figure (13), some of them have non-linear curves that are really hard to approximate, especially at low Reynolds numbers. In this sense, the approximation depends on personal criteria and it is hard to achieve good accuracy. Following with the with the example of figure (13), two lift slopes can be calculated and included, from 0 to 3 degrees and then from 3 to the stall angle, being in this case approximately 9 degrees.

Verication of a Modelica Helicopter Rotor Model Using Blade Element Theory

Verification of a

Modelica Helicopter

Rotor Model Using

Blade Element Theory

Author: Jorge Luis Lovaco Hernandez

Examiner: Ingo Staack

Supervisor: Magnus Sethson

Verification of a

Modelica Helicopter

Rotor Model Using

Blade Element Theory

Examiner: Ingo Staack

Supervisor: Magnus Sethson

Abstract

Acknowledgements

Nomenclature

Contents

List of Figures

List of Tables

1

Introduction

Contribution

2

Theory

3

Method