Linköpings universitet
Linköping University | Applied Thermodynamics and Fluid Mechanics (MVS)
Bacherlor’s thesis, 15 ECTS | CFD
Spring/ 2019 | LIU-IEI-TEK-G--19/01696--SE
CFD analysis of a glider aircraft
Using different RANS solvers and introducing improvements in
the design
David Perez Sancha
Supervisor : Roland Gårdhagen Examiner : Hossein Nadali Najafabadi
Copyright
The publishers will keep this document online on the Internet - or its possible replacement - for a period of 25 years starting from the date of publication barring exceptional circumstances.
The online availability of the document implies permanent permission for anyone to read, to down-load, or to print out single copies for his/hers own use and to use it unchanged for non-commercial research and educational purpose. Subsequent transfers of copyright cannot revoke this permission. All other uses of the document are conditional upon the consent of the copyright owner. The publisher has taken technical and administrative measures to assure authenticity, security and accessibility.
According to intellectual property law the author has the right to be mentioned when his/her work is accessed as described above and to be protected against infringement.
For additional information about the Linköping University Electronic Press and its procedures for publication and for assurance of document integrity, please refer to its www home page: http://www.ep.liu.se/.
Abstract
In this study, Computational Fluid Dynamics (CFD) simulations have been carried out in order to investigate and improve the performance of the Standard Cirrus glider, using different Navier-Stokes methods and solving the equations for the steady flow. The work has been divided in two parts:
First, a study is performed to test the quality of the transition model (γ-Reθ). The two di-mensional results of the glider´s airfoil are compared against the results from panel’s meth-ods and the open-source CFD codes: SU2 and OpenFoam. In addition, three dimensional glider´s models are simulated using the transition model with the purpose of creating a validated reference model of the glider’s performance in steady level flight.
The simulations are carried out in two dimensions for the outer wing airfoil for a 1.5 e+06 Reynolds number and in three dimensions for the Wing & Fuselage model and Tail & Fuse-lage model under a range of velocities. Both simulations are validated against experimental data.
In the second part of the study, the validated model is used to developed possible im-provements in the glider´s external geometry that could produce possible benefits in the performance and handling qualities of the glider.
Keywords:CFD, RANS, Standard Cirrus, glider, Fluent, SU2, OpenFOAM, SST, transition, γ-Reθ, winglet.
"Gliders, sailplanes,
they are wonderful flying machines.
It’s the closest you can come to being a bird.”
Acknowledgments
The author would like to thank T. Hansen for providing the glider geometry for this study. I would also like to express my profound gratitude to Linköping University for providing the computational resources and facilities to carry out this thesis.
Without the continuing support of my colleague and friend Luis G , I would probably not had been able to present such results.
I would especially like to express my gratitude to my thesis supervisor Roland Gårdha-gen, for the advice and support he provided during the course of this project. It has been a great rewarding experience to work with him, in which I have learned much more than I could have initially imagined.
I wish to thank my brother Carlos for being a model for me, encouraging me and taking the time to help me with the difficulties that have arisen through this venture.
Last, but certainly not least, I would like to thank my friends who have been there for me throughout this whole journey, and especially to my parents, for their constant love, their support and their unwavering belief that I can do anything I set my mind to.
Nomenclature
Abbreviations and Acronyms
Abbreviation Meaning
CFD Computational fluid dynamics
CAD Computer aided design
CG Center of gravity
CP Center of pressure
DES Detached Eddy Simulation
DNS Direct Numerical Simulation
k-ω SST Menter’s Shear Stress Transport (turb. model)
γ-Reθ Gamma Re-Theta Transition (turb. model)
SA Spalart–Allmaras (turb. model)
LES Large Eddy Simulation
TKE Turbulent Kinetic Energy
RANS Reynolds-Averaged Navier-Stokes
Latin Symbols
Symbol Description Units
A Reference area m2
AoA Angle of attack [deg]
CD Drag coefficient [´]
Cf Skin-friction coefficient [´]
CL Lift coefficient [´]
Cm Pitching moment coefficient [´]
e Oswald efficiency factor [´]
FLi f t Lift force [N]
FDrag Drag force [N]
g Gravity m/s2
p Pressure [Pa]
m Mass [Kg]
Reθt Momentum thickness Reynolds number [m]
y+ Dimensionless wall distance [´]
tturb Time the flow remains turbulent [s]
tlam Time the flow remains laminar [s]
Greek Symbols
Symbol Description Units
α Angle [degree]
ϕ Aerodinamic efficiency parameter [´]
ρ Density kg/m3
Θ Descent glide angle [deg]
µ Dynamic viscosity kgm´1s´1
e Dissipation rate m2s´3
µt Eddy or Turbulent viscosity kgm´1s´1
Symbol Description Units
τ Shear Stress [Pa]
κ Turbulent kinetic energy m2s´2
λ Wavelength / characteristic length scale [m]
Subscripts and superscripts
Abbreviation Meaning
Contents
Abstract iv
Acknowledgments vii
Contents viii
List of Figures x
List of Tables xii
1 Introduction 1
1.1 Motivation . . . 2
1.2 Background . . . 2
1.3 Aim . . . 3
1.4 The Standard Cirrus . . . 3
2 Theory 5 2.1 Laminar and turbulent flows . . . 5
2.2 Transition . . . 7 2.3 Separation . . . 8 2.4 Principles of Gliding . . . 9 2.5 Flight mechanics . . . 10 2.6 Drag Polar . . . 11 2.7 Improvements . . . 12 3 Method 14 3.1 Navier-Stokes Solver . . . 14 3.2 Turbulent scheme . . . 16
3.3 Two dimensional calculations . . . 18
3.4 Three dimensional calculations . . . 25
3.5 Procedure . . . 29
3.6 Original design improvements . . . 33
4 Results 35 4.1 Two dimensional mesh . . . 35
4.2 Two dimensional results . . . 36
4.3 Three dimensional results . . . 44
5 Discussion 52 5.1 Method . . . 52
5.2 Results . . . 53
5.3 Delimitations . . . 58
Bibliography 60
A Appendix 63
A.1 SU2 Configuration file . . . 63 A.2 OpenFOAM Configuration files . . . 71 A.3 Tables of performance results . . . 77
List of Figures
1.1 TKE spectrum, images extracted from [Arvidson2017] . . . . 2
1.2 Extracted from [standard] . . . 4
2.1 Turbulent flow properties, extracted from [Versteeg1995] . . . . 5
2.2 Typical eddie structure formation on smoke.. . . 6
2.3 Momentum-transport (left) and velocity profiles (right) of laminar & turbulent typical boundary layer. Images extracted from [Bertin] . . . 6
2.4 Boundary layer transition. . . 7
2.5 Schematic of the transition approach. . . 7
2.6 Boundary layer transition. . . 8
2.7 Extracted from: http://yorksoaring.com/index.php/what-is-soaring/how-gliders-fly 9 2.8 Diagram of forces acting on the the glider during steady flight . . . 10
2.9 Schematic of the induced drag process. . . 12
2.10 Total drag decomposition. . . 13
2.11 Examples of winglet mods for a Standard Cirrus glider . . . 13
3.1 Normalized airfoils. . . 18
3.2 Schematic of the boundary conditions applied through the domain. . . 20
3.3 CAD models used for the geometry . . . 25
3.4 Three dimensional mesh. . . 27
3.5 Domain with the correspondent boundary conditions. . . 28
3.6 Estimation of the needed AoA for each velocity. . . 29
3.7 Surfaces used in both models to calculate the aerodynamic coefficients. . . 29
3.8 Standard Cirrus glide polar. . . 30
3.9 Interpolation of the AoA for each velocity condition. . . 31
3.10 Winglet design parameters. . . 33
3.11 Renderings of the proposed winglet configuration. . . 34
4.1 Domain view and closed view of the meshes. Structured (left), Unstructured (right) . . . 35
4.2 Comparison of pressure coefficient distribution, [Tn2019]. . . . 36
4.3 Comparison of lift coefficient versus AoA (for structured mesh), [Tn2019],[Hansen2014]. . . . 36
4.4 Comparison of lift coefficient versus AoA (for unstructured mesh), [Tn2019],[Hansen2014]. . . . 37
4.5 Comparison of lift coefficient versus AoA (for unstructured mesh), [Tn2019],[Hansen2014]. . . . 37
4.6 Comparison of lift coefficient versus drag coefficient (for structured mesh), [Tn2019],[Hansen2014]. . . 37
4.7 Comparison of lift coefficient versus drag coefficient (for unstructured mesh), [Tn2019],[Hansen2014]. 38 4.8 Airfoil transition position.. . . 38
4.9 Performance comparison. . . 38
4.10 Turbulent kinetic energy prediction at α = 0 degrees for turbulent model (SST) . . . 39
4.11 Turbulent kinetic energy prediction at α = 0 degrees for transition model (γ-Reθ). . 39
4.12 Turbulent kinetic energy prediction at α = 8.05 degrees for turbulent model (SST) . 40 4.13 Turbulent kinetic energy prediction at α = 8.05 degrees for transition model (γ-Reθ). 40 4.14 Separation for high AoA (13.5 degrees). . . 41
4.15 Sensitivity of the Turbulent Kinetic Energy (TKE) in the airfoil surface for different AoA. . . 42
4.16 Sensitivity of the Skin Friction Coefficient in the airfoil surface for different AoA. . . 43
4.17 Sensitivity of the Turbulent Shear Stress in the airfoil surface for different AoA. . . 43
4.18 Non-dimensional TKE contour over both Wing & Tail models [95 Km/h] . . . 44
4.19 Comparison of the top side TKE transition for Wing models Left: Hansen´s result [Hansen2014], Right: Obtained [95 Km/h] . . . . 45
4.20 Comparison of the bottom side TKE transition for Wing models Left: Hansen´s result [Hansen2014], Right: Obtained [95 Km/h] . . . 45
4.21 Fuselage TKE transition for Wing models Top: Hansen´s result [Hansen2014], Bottom: Obtained [95 Km/h] . . . . 45
4.22 Tail section TKE transition for Tail model [95 Km/h] . . . 46
4.23 Pressure coefficient contour along the wingspan (plane X=2.1) [95 Km/h] Left: Initial geometry Right: Winglet addition. . . 46
4.24 Pressure distribution along the wingspan (plane X=2.1) [95 Km/h]. . . 46
4.25 Cpcoloured vortex core region based on a swirling strength of 0.0005 [90 Km/h] . . . 47
4.26 Eddie viscosity contour over the vortex swirling strength region and on three separated planes on the wing´s wake [90 Km/h] .. . . 47
4.27 Velocity streamlines in the wing-tip [90 Km/h] . . . 47
4.28 Lift coefficient versus AoA for the steady flight condition, [Ing2012]. . . . 48
4.29 Lift coefficient versus drag coefficient for the steady flight condition, [Ing2012]. . . . 48
4.30 Aerodynamic efficiency versus calibrated airspeed for the steady flight condition, [Ing2012]. . . . 49
4.31 Sink rate versus calibrated airspeed for the steady flight condition, [Ing2012]. . . . 49
4.32 Sink rate versus CAS (steady flight) with best gliding ratio, [Ing2012].. . . 49
4.33 Estimated drag polar for the different studied velocities. . . 50
4.34 Pressure and viscous drag contributions for the different studied velocities. . . 51
5.1 Zoom view from the pressure coefficient distributions, highlighting the critical points in the separation bubble. . . 53
5.2 Oil Flow Tests on a Standard Cirrus 60 with a wing twist of 0.75 deg [70 kn], extracted from [standard]. 55 5.3 Surfaces used in both models to calculate the aerodynamic coefficients. . . 58
List of Tables
3.1 Structured mesh independence study: α = 0 deg, γ-Reθtransition model. . . 19
3.2 Unstructured mesh independence study: α = 0 deg, SST turbulence model. . . . 19
3.3 Boundary conditions applied in the mesh. . . 20
3.4 Conditions applied at the inlet boundary. . . 20
3.5 Turbulent property values. . . 22
3.6 Lift coefficient validation (SST). . . 24
3.7 Mesh independence study Wing & fuselage model. . . 26
3.8 Mesh independence study Fuselage & Tail model. . . 26
3.9 Boundary conditions applied in the mesh. . . 27
3.10 Velocity range applied at the velocity inlet. . . 28
3.11 Basic data of Cirrus75 D-6607, extracted from [Ing2012] . . . . 31
3.12 Steady flight conditions, extracted from [Ing2012] . . . . 32
3.13 Lift coefficients for the total glider model. . . 32
3.14 Steady flight conditions, extracted from [Ing2012] . . . . 32
3.15 Reliability of citeIng2012 for steady flight conditions. . . 33
A.1 Model without tail lift contribution and steady AoA estimation. . . 77
A.2 Drag contributions of the complementary surfaces of both models for 0 and 3 de-grees AoA. . . 77
A.3 Steady flight AoA performance for the complete glider. . . 77
A.4 Model without tail lift contribution and steady AoA estimation. . . 78
A.5 Drag contributions of complementary surfaces of both models for 0 and 3 degrees AoA. . . 78
1
Introduction
When an engineer is tasked with designing a new product or improving an existing one, aero-dynamics usually play an important role in the engineering process. However, aerodynamic processes are not easily quantifiable during the concept phase and usually the only option for optimizing the design is to conduct physical tests on product prototypes. However, Com-putational Fluid Dynamics allows a detailed understanding of the flow field since the first stages, allowing to quantify the impact of possible modifications.
As CPUs have become more powerful in the last decades, most of the companies inter-ested in the field have decided to use CFD software as an alternative to the expensive ex-perimental testing that is often a barrier to new developments. However, its implementation requires from the adequate hardware and investing in commercial licenses which are usually not very affordable. Therefore, open source software provides a cheap alternative to carry out simulations without the cost of commercial software. Nonetheless, they required from a much more knowledgeable user to choose the suitable parameters for the case under study and the documentation available might be scarce.
This mentioned computer revolution has achieved the development of reliable turbulence models, which provide accurate simulations under a wide range of fully turbulent engineer-ing flows. Nonetheless, the effect of the laminar-turbulent transition is not often captured by the simulations. This is due to the fact that the mechanisms that develop the instabilities that induce the transition from a boundary laminar layer to turbulence are highly dependent of the case under study and therefore difficult to predict with general-purpose formulations, which often involve carrying out calibrations that allow applying the method to the specific case.
1.1. Motivation
1.1
Motivation
High performance gliders can hugely benefit from the possibilities that CFD offers to improve existing designs and optimise their performance and handling qualities. Therefore, current Computational Fluid Dynamics (CFD) has become an essential tool for predicting new de-signs since the completion of the design of the JS-1 sailplane, where the benefits of CFD were clearly demonstrated. In this study, it was stated that a single design iteration suppose 6 months of wind tunnel testing while only one day of CFD analysis, consequently, the equiva-lent wind tunnel time for a 3 month CFD design period would be 45 years (6 x 30 x 3 months) [29].
However, the flows around these aircrafts are complex and accurate tools are required in order to capture the existing transition from laminar to turbulent flow.
1.2
Background
The important development of technology and modelling capabilities has led simulations to be conducted on multiple millions of cells. Within the current computational landscape, it is possible to parallel processing and with world records such as the one achieved by ANSYS in partnership with the University of Stuttgart and Cray Inc, scaling ANSYS Fluent to over 172,000 computer cores, it has been proven to be feasible to simulate a case down to the characteristics of the smallest turbulent eddies without the need of modelling the physics.
However, despite the computational power available nowadays, these direct numerical simulations (DNS) often lack practicality due to their immense cost, as it requires solving the mean flow and all turbulent velocity fluctuations on a sufficiently refined grid, so it can resolve the Kolmogorov length scales at which the energy dissipation takes place, and with time steps sufficiently small to resolve the period of the fastest fluctuations. Nonetheless, Spalart [34] predicted that such simulations will become readily available for development towards the end of the century thanks to the progress of computers, whereas progress on turbulence modelling can stagnate.
The limiting factor in the efficiency of CFD simulations is usually the complexity of turbu-lence. Therefore, until technology is developed enough to allow for DNS to be widely used throughout the industry, the simplification and approximation of turbulence modelling has been a primary focus. The three main possibilities of simulating a modelled turbulence are: Reynolds Averaged Navier Stokes (RANS), Large Eddy Simulation (LES) and Detached Eddy Simulations (DES).
(a) Spectrum of the different ranges of turbu-lent kinetic energy.
E is the energy per wave number and κ is the wave number.
(b) Spectrum for turbulent kinetic energy over-laid by turbulence modelling techniques with ranges of resolved and modelled turbulence. Horizontal dashed grey lines: modelled turbu-lence; solid black lines: resolved turbulence.
1.3. Aim
RANS simulations time-average the fluid flow while solving for the Reynolds stresses, modelling the extra terms that appear in the time-averaged flow equations with a turbulence model.
LES models turbulence by ignoring all eddies smaller than the size of the smallest dynam-ically significant length-scale within the specific flow. This cut-off length is generally equal to the grid size and hence, only retains the scales within the inertial sub-range determined by the grid size.
The smaller and more numerous dissipative scales are modelled with a subgrid scale (SGS) model. This method requires space filtering the unsteady Navier-Stokes equations prior to the computations, passing the larger eddies and rejecting the smaller ones.
Finally, DES is the combination of the two previous techniques, whereas RANS is used to model the complicated near wall region, LES is used to solve the rest of the flow.
While LES and DES are commonly used in acadaemia due to the improved accuracy over RANS, especially in flows with massive separation, where RANS, steady or unsteady, is not able to predict the flow behaviour with acceptable tolerances. In most industries, its cost is a limiting factor and they opt to use the time-averaged option. However, Spalart [34] predicted that they will become readily available by mid century.
Flows with shallow or no separation appear to be within the reach of the current steady RANS methods or their finely calibrated derivatives. For such flows, transition prediction with generality, accuracy, and robustness may prove more challenging than simply turbu-lence prediction. This is the case of the present study, in which for the reduced time and computer cost, solving above the time-averaged steady equations goes beyond the scope of a modest study.
1.3
Aim
The study aims for two main objectives. First, the flow conditions around the Standard Cir-rus glider are simulated by solving the Reynolds-Averaged Navier-Stokes equations (RANS) and compared against experimental data in order to create a reference model. Second, the created model is used for investigating new wing geometry designs that may improve the aerodynamic performance.
The validation starts by analysing the performance of the FX 66-17 A II-182 airfoil which is located in the outer part of the Standard Cirrus wing. The two dimensional geometry is inves-tigated in detail by comparing the computation obtained from different CFD models against the experimental values from the low-turbulence pressure wind tunnel at NASA Langley [8]. Afterwards, the three dimensional glider model is simulated in steady flight for velocities between 90 Km/h to 160 Km/h. Finally, the conducted simulations are validated against the results from the flight tests performed with a Standard Cirrus at the Idaflieg summer meeting in 2011 [16].
1.4
The Standard Cirrus
The Standard Cirrus (Figure 1.2) is a 15-m span and no camber-changing flaps that was de-signed by Dipl.-Ing. Klaus Holighaus and was built at the Schempp-Hirth factory to compete in the Standard Class.
This specific glider competition class was introduced in the late fifties as an alternative to the increasingly heavy, difficult to fly and costly Open Class aircrafts of that time. Striv-ing for affordability and simplicity, the original standard class had a restricted span of 15 metres, fixed wing sections and maximum all-up mass of 525 kg. Moreover, retractable un-dercarriages, flight-disposable ballast, radios and lift-enhancing devices such as flaps were not allowed in this glider competition class.
1.4. The Standard Cirrus
The first flight took place in March 1969 and was equipped with an all-moving tailplane, air brakes on the upper surface of the wings and the possibility to carry up to 80 kg of water ballast to increase the flight performance.
The wing of the glider, as found in [35], is designed using two different airfoils, blending the root airfoil linearly into the airfoil that is used at the outer part of the wing. This outer airfoil is kept constant from the start of the aileron to the tip of the wing. The glider is widely known for its good handling qualities, large cockpit and ability to climb well in turbulent thermals. Therefore, despite the last Cirrus model was built in 1985, is still considered to be one of the best gliders for participating in club class competitions.
Some of the aircraft‘s most notable specifications are listed below [35]:
General characteristics
• Crew: One
• Capacity: 80 kg (176 lb) water ballast • Length: 6.35 m (20 ft 10 in)
• Wingspan: 15.00 m (49 ft 3 in) • Height: 1.32 m (4 ft 4 in) • Wing area: 10.0 m2 (108 ft2) • Aspect ratio: 22.5
• Wing profile: FX S-02-196 modified • Empty weight: ca. 215 kg (473 lb) • Gross weight: 390 kg (860 lb)
Performance
• Stall speed: 62 km/h (33 kn) • Never exceed speed: 220 km/h • Maximum speed: 220 km/h (140 mph) • Maximum glide ratio: 38.5:1 at 90
km/h (49 kn)
• Rate of sink: 0.6 m/s (120 ft/min)
2
Theory
2.1
Laminar and turbulent flows
Laminar flow is characterised by parallel streamlines with high diffusion of momentum be-tween particles, while the dissipation of momentum is very low, i.e. no momentum transport in the direction normal to the streamlines. Moreover, the kinetic energy of the flow is damp-ened by the viscosity of the fluid, however, when the viscosity of the fluid can no longer dampen all the disturbances, eddies (vortices) are formed and the flow turns turbulent. Tur-bulent flow is characterised by fluctuating quantities relative to the mean, as can be observed in Figure 2.1 (a).
(a) Typical velocity measurements in a turbulent flow
u(t) =U+u‘(t)
where U is the steady mean value and u‘ the fluctuating com-ponent.
(b) Energy spectrum of turbulence (Kolmogorov scales).
2.1. Laminar and turbulent flows
Figure 2.2:Typical eddie structure for-mation on smoke.
In addition, the vorticity of the flow field also fluctu-ates, which forms spinning structures called eddies. This fluid structures, such as the one depicted in Figure 2.2, rep-resent the chaotic and random swirling of a fluid.
The size of these eddies varies in the range of the tur-bulent region and increases with Reynolds number. The kinetic energy is handed down from large eddies to pro-gressively smaller ones in a process known as energy cas-cade, in which the largest eddies collapse into smaller and smaller eddies until the eddies are small enough for vis-cous forces to dissipate the energy. This is due to the fact that in this length and velocity scale, the viscous term is at least as important as the non-linear convective term, and therefore the viscous damping should be strong enough to smooth away any velocity gradient.
Figure 2.1 (b) gives the spectral energy E(κ)as a funtion of the wavenumber:
κ = 2 π / λ, where λ represents the characteristic length of the viscous eddies. The
dia-gram shows that the energy content peaks at the low wavenumbers, so the larger eddies are the most energetic. Nonetheless, the value of E(κ)rapidly decreases as the
wavenum-ber increases, so the smallest eddies have the lowest energy content. This means that the dissipation of energy is high and thus a persistent energy supply is necessary to maintain a turbulent flow. Particularly, the high diffusion of momentum in this flows, produces large wall shear stresses and thus increased drag compared to a laminar flow and therefore higher fuel consumption which is both expensive and increases the environmental impact.
The boundary layer characteristics are highly determined by the type of flow. The typical velocity profiles of laminar and turbulent boundary layers are presented in Figure 2.3, where can be observed how the turbulent flow increases velocity more rapidly away from the wall. The velocity gradient throughout the boundary layer gives rise to internal shear stresses which lead to friction known as skin-friction drag. This drag is predominant in streamlined flows where the majority of the body’s surface is aligned with the flow. As the velocity gradi-ent at the surface is greater for turbulgradi-ent than laminar flow, a streamlined body experiences more drag when the boundary layer flow over its surfaces is turbulent.
Figure 2.3:Momentum-transport (left) and velocity profiles (right) of laminar & turbulent typical boundary layer. Images extracted from [5]
2.2. Transition
2.2
Transition
Figure 2.4:Boundary layer transition. When the laminar boundary layer
develops in the streamwise direc-tion, it may be subjected to numer-ous disturbances. If these distur-bances are not damped, they might get amplified and force the transition to turbulent flow, as can be observed in Figure 2.6.
There are three fundamental stages through which laminar flow becomes turbulent: receptivity, dis-turbance growth and breakdown [32]. Receptivity controls how exter-nal perturbations enter the bound-ary layer and excite internal distur-bances. Instability how the latter grow, and breakdown when and where the flow first becomes turbu-lent.
The main factors responsible for the onset of the transition have been proposed to be the following transition mechanisms: natural transition, bypass transition, separation-induced transition and wake induced transition [9].
In external aerodynamic flows, with a very low mean-flow turbulence level (0.1% as-sumed for the 3D studied case), the transition mechanism is often the natural transition and typically occurs as a result of the Tollmien–Schlichting waves or cross-flow instability. The resulting exponential growth eventually results in the appearance of secondary instabilities with a three-dimensional nature and consequently a nonlinear breakdown to turbulence [24]. The T-S disturbances usually occur in the rear portion of the wing, decelerating the flow, whereas for the forward part of a wing, when the flow is accelerating in the streamwise di-rection, Crossflow disturbances are intensified by the favourable pressure gradient [32]. Even though both types of waves are present in typical swept-wings, transition is usually caused by either one, but not both, of these waves [30].
For a more complete review, the interested reader is referred to [14] or [30].
2.3. Separation
2.3
Separation
Over flat surfaces we can safely ignore any changes of pressure in the flow direction. Under these conditions, the boundary layer remains stable but grows in thickness in the flow direc-tion. This is, of course, an idealised scenario and in real-world applications, such as curved wings, the flow is most likely experiencing an adverse pressure gradient, which occurs when the static pressure increases in the direction of the flow. Under these conditions the bound-ary layer can become unstable and separate from the surface. The boundbound-ary layer separation induces a second type of drag, known as pressure drag. This type of drag is predominant for non-streamlined bodies, such as an aircraft wing at a high angle of attack.
The second, and more dramatic effect, of boundary layer separation in aircraft wings is aerodynamic stall. At relatively low angles of attack, which are typical for cruise conditions, the adverse pressure gradient acting on the top surface of the wing is benign and the bound-ary layer remains attached over the entire surface. However, if the angle of attack is increased, so does the pressure gradient and at some point the boundary layer will start to separate near the trailing edge of the wing, moving further upstream as the angle of attack is increased. If the airfoil is positioned at a sufficiently large angle of attack, separation will occur very close to the point of maximum thickness of the airfoil and a large wake will develop behind the point of separation.
Figure 2.6:Boundary layer transition. This wake redistributes the
flow over the rest of the airfoil and seriously reduces the produced lift in a condition known as aero-dynamic stall. Flows at moder-ate Reynolds number, typical in most gliders, are not able to sus-tain strong adverse pressure gra-dients and often separate in lam-inar flow regime. The turbulence developing inside the recirculation region enhances the momentum transport, and the flow reattaches, forming a laminar separation bubble.
A sketch of the typical structure of this process is shown in Figure 2.6. After the separation, a large part of the separated zone is characterized by a slow flow motion. This is named as dead-air region. The last part of the bubble presents a strong re-circulation vortex flow. Finally, a sudden pressure recovery leads to the reattachment of the flow.
Laminar separation bubbles are an important scale effect that must be understood by designers and researchers as their presence can cause deterioration in the performance of the airfoil.
2.4. Principles of Gliding
2.4
Principles of Gliding
There are three forces acting upon the glider during steady level flight. They are lift, gravity and drag, as is observed in Figure 2.7.
Figure 2.7: Extracted from: http://yorksoaring.com/index.php/what-is-soaring/ how-gliders-fly
Once the glider has been launched, towed, or winched, the thrust from the external source disappears. The glider is capable of maintaining the velocity needed for steady flight by con-verting the potential energy that it has accumulated into kinetic energy as it glides downward in an unaccelarated descent, trading height for distance. Nonetheless, the glider‘s lift, drag and glide ratio characteristics are governed solely by its construction, and are predetermined at takeoff. Without thrust, the only other characteristic that the pilot has control over (besides normal control surfaces) is the weight of the plane.
A heavier glider will sink faster than a light glider. However, the glide ratio is not affected by weight because while a heavier glider will sink faster, it will do so at a higher airspeed. The plane will come down faster, but will cover the same distance (at a higher speed) as a lighter glider with the same glide ratio and starting altitude. In order to help them fly faster, the Standard glider have tanks that can hold up to 60 kg of water.
Since gliders gain altitude by locating in a pocket of air that is rising faster than the glider is descending, the downsides of heavier sailplane includes reduced climbing rates in these lifting environments and, possibly, shorter flight duration if suitable lift cannot be located. To prevent this, the water ballast can be jettisoned at any time through dump valves, allowing the pilot to reduce the weight of the plane to increase climb rates, or to reduce speed as they come in for a landing.
2.5. Flight mechanics
2.5
Flight mechanics
The balance of forces (Figure 2.8) in a gliding steady flight meet the following equations: Li f t ¨ cos(θ) +Drag ¨ sin(θ) =W Li f t ¨ sin(θ) =Drag ¨ cos(θ) (2.1)
Or the equivalent:
Li f t=W ¨ cos(θ) =m ¨ g ¨ cos(θ) Drag=W ¨ sin(θ) (2.2)
Therefore, the lift force coefficient can be expressed as: CL= L
q8¨S
= m ¨ g ¨ cos(θ) q8¨S
(2.3) In addition, θ is usually a small value (around 1.6 for our glider), therefore cos(θ) « 1, assum-ing the followassum-ing simplified form:
CL = m ¨ g q8¨S
(2.4) The drag force coefficient can be expressed as:
CD= D q8¨S
(2.5) For the equations above, m is the mass of the glider, g is the gravitational constant and S is the reference area. The dynamic pressure q8is denoted:
q8= 1 2ρ8¨V
2
8 (2.6)
where ρ8is the density of air and V8is the free-stream velocity. Since the change in Reynolds number due to difference in density at different altitudes is small, the descent glide angle θ can be found from:
θ= 1
CL/CD (2.7)
Hence, the descent glide angle θ is only a function of the lift-to-drag ratio, CL/CD, and does not depend on altitude or wing loading. However, to achieve a given CL/CD at a certain altitude, the glider must fly at a specific velocity V8 called the equilibrium glide velocity. The value of V8is dependent on both altitude and wing loading. The rate of sink for each equilibrium velocity is given by:
h=V8¨sin(θ) (2.8)
2.6. Drag Polar
2.6
Drag Polar
A generic drag polar curve fits the expression CD = a + b CL + c CL2, but for the the following calculations, it is possible to simplify the expression by a parabolic function with great precision. CD=CD0+CDi =CD0+k ¨ C2L=CD0+ C2 L π Aϕ $ ’ ’ ’ ’ ’ & ’ ’ ’ ’ ’ %
CD0”Zero-lift drag coefficient CDi”Lift-induced drag coefficient k ” Unitary induced resistance parameter A ” Wing aspect ratio
ϕ ”Aerodynamic efficiency factor
, / / / / / . / / / / / -(2.9) As shown in eq. 2.9, the additional resistance that appears when the lift is not zero, re-ferred as induced drag, can be related to the square of the CLthrough the induced resistance parameter k.
In addition, the induced drag has two main contributions: one due to the variation of the friction on the surfaces as a consequence of the modification of the velocity field that the generation of lift induces, and another from the whirlpools of the marginal edge of the wing [27].
It is worth bearing in mind that some authors, when talking about resistance induced by the lift, they want to refer exclusively to the latter. An additional difference comes from the distribution of the lift between the wing and the horizontal stabilizer. If the lift distribution is forced to get the balance of the pitch moment, the consequent drag contribution is called balancing resistance. Therefore, it is referred as unbalanced polar (simplification consisting in supposing that the tail does not sustain), and of balanced polar (tail achieves the necessary lift to fly in pitching moments balance).
For the present study, as not having access to the information for balance, it will be used the unbalanced polar.
• The contribution due to the swirling wake is estimated from the value of CDi for a wing with elliptic lift corrected with the Oswald parameter (e):
CDi1= C2L
π Ae (2.10)
The Oswald factor usually has a value close to unity, e=1 corresponds to a wing with a elliptical lift distribution which minimizes the induced drag and to which we intend to get close.
• The viscous contributions due to the friction and pressure drag can be estimated from the polar of the two dimensional airfoil: cd =cd0+kac2l , by integrating the sum of the kac2l term along the span (the sum of cd0 has already been implicitly taken into account in the calculation of the CD0for the complete glider):
Finally, the induced drag coefficient can be subdivided as following:
CDi=kac2l + C 2 L π Ae = CL2 π Aϕ (2.11)
2.7. Improvements
2.7
Improvements
As the wing is driven through the air to develop the difference in air pressures which allows to generate a force known as lift, inevitably, another force is generated which is known as induced drag.
As can be noticed in Figure 2.9 (a), the higher pressure air on the lower surface of the airfoil located at the end of the wing, curves around the wing and fills in the lower pressure area on the upper surface. Therefore, much of the the lift is lost while the energy needed to produce the difference of pressures is still expended.
The result of this wasted energy is an unnecessary increase of drag force, which is spent producing wing-tip vortices such as those shown in the Figure 2.9 (b). Their formation is due the mixing of the streams of air from above and below the wing with the spanwise component of velocity going in opposite directions. When they combine along the trailing edge of the wing, vortices are formed which, when viewed from the rear, rotate clockwise from the left wing and counter clockwise from the right (represented in Figure 2.9 (a)). These vortices tend to move outwards towards the wing tip, joining together and forming within a short distance downstream, one large wing tip vortex at each of the wing tips, as it is represented in Figure 2.9 (b). Therefore, although they are often referred to as "wingtip vortices” because of the structures that line up fairly closely behind the wingtips, their presence is due the vorticity that feeds into the cores from the entire span of the wing trailing edge.
The more energy the glider requires to fly, the greater the required rate of descent is to sup-ply sufficient energy to convert into thrust to overcome that unnecessary drag. Consequently, as the energy that produces the vortices is wasted energy which translates into poorer per-formance, the purpose of improving the existing design is to convert as much of the energy as possible into useful lift and the necessary thrust. One way to attempt to reduce drag is to increase the aspect ratio, because the proportion of air which moves in this way is reduced and therefore more of it generates lift.
(a) Lateral air flow.
(b) Induced drag by wing tip vortices.
2.7. Improvements
Winglets
Wingtip devices, usually known as winglets, are also used to improve the efficiency of the glider. They are intended to reduce the aircraft’s drag by altering the airflow near the wingtips with an increase in the effective aspect ratio of a wing but without materially in-creasing the wingspan. Nonetheless, rather than being a simple fence that limits the span-wise flow, winglets try to create a flowfield that actively interacts with that of the main wing to reduce the total amount of spanwise flow. The sought effect is that the downwash (side-wash) produced by the winglet opposes the spanwise flow on the main wing, spreading out the influence of the tip vortex.
However, as their intention is to act like an endplate but with much less wetted area by carrying the proper aerodynamic loading, they also suffer the penalty of increasing the profile drag and consequently, the main designing goal is to obtain the largest reduction in induced drag for the smallest increase of profile drag. The crossover point of the trade-off between reducing the vortex induced drag against the penalty of additional profile, as stated by D. Maughmer in [22] is represented by the following equation:
∆DPROFILE=∆DINDUCED (2.12)
This simple expression indicates that the more the induced drag can be reduced for a given increase in profile drag, the higher will be the crossover point and therefore the more effective the winglet.
Figure 2.10:Total drag decomposition. The most notable drag benefits of winglets is
achieved for relatively low velocities, where the re-quired lift coefficient is high and therefore the in-duced drag is predominant, allowing the device to displace the concentrated vorticity away from the wing.
However, the profile drag penalty decreases with the velocity, as its contribution is mainly de-termined by the wetted area.
With the possible benefit and penalty being at different points in the flight regime, the optimum winglet geometry becomes fairly complicated and requires evaluating the changes in performance due to winglets over the entire flight envelope of the glider, because the global performance combines a low-speed, high-lift coefficient climb phase and a high-speed, low-lift coefficient cruise phase, both of relatively similar importance.
(a) Winglet mod by Steve Willits. (b) Blended winglets by Nick Gilbert.
3
Method
3.1
Navier-Stokes Solver
The numerical resolution of the Navier-Stokes equations in their incompressible form using the finite volume method, requires from powerful and robust codes. Although the possibil-ities are increasingly numerous for both commercial and open-source softwares, the choice should not be arbitrary. In the following sections, the selected options are explained.
Fluent
ANSYS® Fluent® software is a general-purpose (CFD) commercial package capable of mod-elling fluid flows, turbulence and other related physical phenomena. Its flow analysis ca-pabilities have been highly validated, obtaining relatively fast and accurate CFD results. In addition, it provides a good amount of tools for the analysis, with an interactive solver set-up, solution, and post-processing making it easy to pause a calculation, examine the results with the integrated post-processing capabilities, change any setting, and then continue the calculation within a single application. For the configuration of this software, it was followed the recommendations in the Fluent Theory Guide [1].
All the CFD simulations and numerical implementations involved in this current work within this solver have been performed in Fluent (version 19.2).
Open Source
The development of open source codes aims for usage on open source operating systems. Therefore, despite it is possible to compile versions for Windows of SU2 and OpenFOAM, in general, they run more stable and faster on open source operating systems and the versions available are more up to date.
Ubuntu has been selected as the open source operating system for this case, as most instal-lation guides and already compiled software downloads are made especially for Ubuntu and it is friendly for beginners.
For the open-source part, Ubuntu was installed on a laptop and all the simulations were carried out with a with 8GB RAM and an Intel® Core, i7-3630QM 2.4GHz processor.
3.1. Navier-Stokes Solver
SU2
The Stanford University Unstructured (SU2) is an open-source software tool written in C++ for the numerical solution of partial differential equations (PDE) and performing PDE con-strained optimization. As one of the primary applications for which this software has been developed is computational fluid dynamics, the performance of the solver coupled with the RANS SA and SST models had been validated against experimental data in previous papers: [10] , [26].
In the present study, its performance has been contrasted against the expectedly robust Fluent results and the experimental data.
The relatively few information available has encouraged the selection of this precise tool because of the innumerable possibilities of its future development as being an open-source software, which allows a wider range of users carry out CFD simulations at a low cost. Some of the features that made this precise tool especially interesting are the capability of monitoring lift and drag monitors and their corresponding contributions or its available design optimization tools with an option to do simple python scripting. Moreover, SU2 also allows for parallel computation and has CGNS mesh support. All the CFD simulations and numerical implementations involved in this current work within this solver have been performed in SU2 (version 6.2.0).
The post-processing was carried out in ParaView, an open-source, multi-platform data analysis and visualization application.
More detailed information about SU2 can be accessed at: https://su2code.github. io/
OpenFOAM
OpenFOAM (Open Field Operation and Manipulation) is an open source multi-physics modelling platform, written in C++ programming language that can perform numeri-cal simulations of fluid flow, combustion, electro-magnetics, heat transfer, etc. The code was under commercial license and owned by the company NABLA up until 2005, when it was released as open source under the GNU license. Nowadays, OpenFOAM is very well known for its fluid flow modelling capabilities using the finite volume method (FVM) and it is currently widely used for CFD simulations. Unlike commercial CFD packages, OpenFOAM provides users access to the source code as well as the flexibility to modify the code, which enables researchers to develop, implement and test various new turbulence models, boundary conditions, etc.
Due to this exclusive features, OpenFOAM is considered a very powerful tool for searchers in academia and widely adopted by engineering companies since its official re-lease in 2004.
All the CFD simulations and numerical implementations involved in this current work within this solver have been performed in OpenFOAM (version 1812).
A more detailed source of information about OpenFOAM can be accessed at: www. openfoam.com/or in the used user guide [13].
3.2. Turbulent scheme
3.2
Turbulent scheme
The flow under which the analysis is carried out is considered to be incompressible as the Mach number = 0.1 and with a relatively high Reynolds number (Re=1.5 ¨ 106).
Menter´s Shear Stress Transport (SST) model
The shear-stress transport (SST) k-ω model used, was developed by Menter in 1994 to deal with the strong freestream sensitivity of the k-ω turbulence model and improve the predic-tions of adverse pressure gradients. It blends the robust and accurate formulation of the k-ω model in the near-wall region with the free-stream independence of the k-e model in the far field. To achieve this, the standard k-ω model and the transformed k- e model are both mul-tiplied by a blending function and both models are added together. The blending function is designed to be one in the near-wall region, which activates the standard k-ω model, and zero away from the surface, which activates the transformed k-e model [1].
By this approach, the attractive near-wall performance of the k-ω can be used without expecting the potential errors resulting from the free stream sensitivity of that model in what is known as the baseline model (BSL).
The SST model incorporates a damped cross-diffusion derivative term in the ω equation and changes the turbulent viscosity definition to account for the transport of the turbulent shear stress, which can be interpreted as a variable Cµ, whereas Cµwas constant in the k-e
model.
These modifications are required to accurately capture the onset of separation under pres-sure gradients, making the SST k-ω model more accurate and reliable for a wider class of flows but from a completely turbulent approach.
γ-Reθ Transition model
The γ-Reθcorrelation based transition model is used in the simulations in order to capture the
laminar-turbulent transition process. Otherwise, fully turbulent flow simulations will most likely completely miss the separation. The model is based on two additional transport equa-tions, one for intermittency (γ ) and one for the transition onset criteria in terms of momentum thickness Reynolds number(Reθt). The model does not attempt to model the physics of the transition process but just implement a model based on experimental correlations.
The intermittency equation is used to trigger the transition process and it determines the percentage of time the flow is turbulent. Consequently, it produces a laminar flow when the intermittency is zero, a fully turbulent when the intermittency is 1 and a transition flow when 0< γ<1.
γ= tturb
tlam+tturb
(3.1) In addition, a second transport equation is solved in terms of the transition onset momentum-thickness Reynolds number (Reθt). This is necessary in order to capture the
non-local influence of the turbulence intensity, which changes due to the decay of the turbulence kinetic energy in the freestream, as well as due to changes in the freestream velocity outside the boundary layer. This new term is calculated outside the boundary layer, therefore in that region, the transport variable is forced to follow the value of Reθtprovided by the
experimen-tal correlation. However, a diffusion term is needed to diffuse it into the boundary layer. Basically, the transport equation takes this non-local empirical correlation, obtained by previous experimental measurements and transforms it into a local variable, Reθt, for the
flow that is aimed to be studied. This local variable can then be compared to the vorticity (or strain-rate) Reynolds number Reν, to determine where in the flow the transition criterion
3.2. Turbulent scheme
transport equation is activated which then produces turbulence. As a result of this mecha-nism, the strong variations of the turbulence intensity and the pressure gradient can be taken into account.
This additional transport equation is an essential part of the model as it ties the empirical correlation to the onset criteria in the intermittency equation and allows the models to be used in general geometries without interaction from the user.
Finally the Fluent Transition model is coupled with the SST model mentioned above, adding a total of four additional equations, two from the SST and two extra from the transi-tion modelling.
The implementation of this model is referred as "Transition SST Model" in the Fluent in-terface and as "Langtry-Menter k-omega Shear Stress Transport (SST)" in the OpenFOAM one.
Fluent by default, calculates three variables from the empirical correlations, based on the free-stream turbulence intensity (Tu=FSTI) and pressure gradients (λθ) [11]. The variables are
respectively: the transition momentum thickness Reynolds number Reθt (3.2), the length of
the transition zone Flength(3.3) and the point where the model is activated in order to match both Reθt and Flength, Reθc(3.4).
Reθt = f(Tu , λ) (3.2)
Flength= f(R˜eθt) (3.3)
Reθc= f(R˜eθt) (3.4)
These empirical correlations are provided by Langtry and Menter [19] (p. 741).
Nonetheless, the STAR-CCM+ results obtained in [15], were calibrated and imposed a free-stream edge for the simulations set to a distance of 50 mm.
Spalart Allmaras (SA) model
This one additional equation model has been developed mainly for aerodynamic flows. The formulation blends from a viscous sublayer formulation, to a logarithmic formulation based on the y+. In as such, no addition of highly non-linear damping functions for lami-nar/viscous sublayer modelling is in use.
It solves a modelled transport equation for the kinematic eddy viscosity without calcu-lating the length scale related to the shear layer thickness. The variable transported in the Spalart Allmaras model is ˜ν, which except from the regions that are affected by strong vis-cous effects such as the near wall region, is identical to the turbulent kinematic viscosity.
For all the turbulent models mentioned, the distance from the wall to boundary to the first cell centroid in the mesh needs to be controlled (by setting the appropiate first layer height) to ensure that the boundary layer is resolved by the turbulent model. This distance is defined by the y+value: y+ = yut ν where ut= c τw ρ (3.5)
where y is the normal distance from the wall to the first cell-centroid, utis the frictional velocity at the nearest wall, ν is the kinematic viscosity and τwis the wall shear stress.
To enable the models to converge, the y+values ought to be in the range between 0.1 to 1, with a relatively low growth rate and stream-wise mesh spacing in the transition area to capture the laminar separation bubble.
3.3. Two dimensional calculations
3.3
Two dimensional calculations
With the purpose of investigating the accuracy of the different turbulent models, the per-formance of the airfoil used on the outer part of the Standard Cirrus wing is investigated for two dimensions. Subsequently, the conducted simulations are validated by comparing the obtained values against the experimental data from the low-turbulence, pressure wind tunnel at NASA Langley [8].
Glider geometry
The Standard Cirrus wing has a linear transition from the root airfoil (FX S 02-196) to the outer airfoil (FX-66-17AII-182) at the inboard end of the aileron. Figure 3.1 shows the airfoil data normalized with the chord length.
Figure 3.1:Normalized airfoils.
This FX-66-17AII-182 airfoil was designed by Dr. F.X Wortmann at the University of Stuttgart and some of its main features are:
• Thickness: 18.2 % • Camber: 3.8 % • Efficiency: 37 • Max CL: 1.397 • Max CL angle: 15.0 • Max L/D: 31.105 • Max L/D angle: 6.0 • Stall angle: -0.5 • Zero-lift angle: -3.5
Mesh
The meshes were created using the Ansys meshing tool and were carefully designed to satisfy the y+ value required (0.1 < y+ < 1) with a fine transition near the airfoil to capture the laminar separation bubble. The high quality structured mesh was created with the purpose of getting accurate results as the computational cost was not restrictive. However, a second unstructured mesh was created in order to carry out the desired simulations on a limited power laptop.
With the purpose of fulfilling the above mentioned criteria, the structured mesh used a face meshing based on the number of divisions along the mesh faces. In addition, a bias factor was applied to obtain the required y+.
However, the unstructured mesh used several bodies of influence to refine the regions around the airfoil and its wake. Moreover, a inflation layer based on the first layer thickness was introduced to fulfil the required y+. It consisted of 56 layers with a 1.1 growth rate.
3.3. Two dimensional calculations
• Grid Independence Study
No matter how refined a grid is created to run the simulations, numerical errors will be introduced due to the discretization of a continuous domain in order to get a finite number of grid points. Therefore, it is important to be aware that the results may be af-fected by the composition of the grid and that the solution should become less sensitive to the grid spacing and approach the continuum solution as the mesh is refined. A grid study was therefore conducted for two different AoA, 0 degrees and 8.5 degrees which was considered restrictively enough, finally, a grid was developed under an iter-ative process.
Table 3.1: Structured mesh independence study: α = 0 deg, γ-Reθtransition model.
Mesh Elements Cl Cd Max speed / Uref
Coarse 237452 0.42000 0.007385 1.4417
Medium 510451 0.42431 0.007611 1.4419
Fine 882763 0.42287 0.007737 1.4426
Difference Coarse-Medium (%) 1.03 3.06 0.02
Difference Medium-Fine (%) -0.34 1.65 0.05
Structured mesh independence study: α = 11 deg, γ-Reθtransition model.
Mesh Elements Cl Cd Max speed / Uref
Coarse 247526 1.25443 0.03083 4.1969
Medium 510451 1.24383 0.0318183 4.1974
Fine 882763 1.24240 0.031945 4.1976
Difference Coarse-Medium (%) -0.85 3.21 0.01
Difference Medium-Fine (%) -0.11 0.40 0.00
Table 3.2: Unstructured mesh independence study: α = 0 deg, SST turbulence model.
Mesh Elements Cl Cd Max speed / Uref
Coarse 176366 0.39601 0.01231 1.4267
Medium 308755 0.39668 0.01221 1.4265
Fine 743868 0.39675 0.01212 1.4258
Difference Coarse-Medium (%) 0.17 -0.79 -0.01
Difference Medium-Fine (%) 0.02 -0.71 -0.05
Unstructured mesh independence study: α = 11 deg, SST turbulence model.
Mesh Elements Cl Cd Max speed / Uref
Coarse 175766 1.35711 0.0371737 2.0791
Medium 328128 1.35028 0.03622 2.0693
Fine 680230 1.34513 0.03591 2.0642
Difference Coarse-Medium (%) -0.50 -2.55 -0.47
3.3. Two dimensional calculations
Finally, for both cases, the medium meshes were employed for the other cases under study, because the relative error caused by the refinement was considered small enough as can be noticed in Tables 3.1 and 3.2. Therefore, a finer mesh would not provide more accuracy in the results, but only increase the computational cost.
Setup
To reproduce the flow conditions of the NASA wind tunnel, the turbulent intensity was set to 0.02% and a 10 turbulent viscosity ratio was used as was considered by Thomas Hansen in [15] based on [12].
The fluid simulated is incompressible air and the boundary conditions applied through the domain are shown in Figure 3.2 and described in Table 3.3.
Table 3.3: Boundary conditions applied in the mesh.
Boundary
Condition Description
Inlet Velocity inlet with a magnitude based on the Reynolds number (conditions in table 3.4) and the direction in the X-component Outlet Static pressure outlet (Gauge pressure = O Pa)
Upper Wall Free slip stationary wall (Shear stress = 0 Pa) Lower Wall Free slip stationary wall (Shear stress = 0 Pa) Airfoil Non slip wall
Interior Interior mesh
Table 3.4: Conditions applied at the inlet boundary. Reynolds number Density Viscosity Velocity
1.5 e+6 1.225 1.7894 e-5 V= Re¨µ
ρ¨L
Figure 3.2: Schematic of the boundary conditions applied through the domain. The velocity imposed at the velocity inlet is calculated based on the values from Table 3.4, where the chord of the airfoil is considered the characteristic length (L). These conditions can be considered representative of the glider´s steady flight. Finally, the airfoil geometry was rotated for the different AoA.
3.3. Two dimensional calculations
˝ SU2
Once downloaded and installed by building the SU2 suite from the source code on Linux, SU2 was ready to run simulations and design problems. In general, all SU2 execution occurs via command line arguments within a terminal, therefore, users can execute the individual C++ programs while specifying the problem parameters in the all-purpose configuration file.
SU2 consists of separate C++ modules whose sequential execution is automated within a Python framework. Each C++ module addresses a compute-intensive task in the pro-cess, and, in order to maintain high computational efficiency, encourage code reuse, and ease the integration of new features, the modules share a common set of classes and data structures within an object-oriented code architecture.
The module used in the study is SU2_CFD, a RANS CFD solver, which also includes a solver for the adjoint RANS equations. The solver requires two inputs: a mesh grid and a configuration file.
The unstructured mesh native format of this solver is a readable ASCII in a extension named .su2. However, the mesh was imported from Ansys meshing in order to facilitate the direct comparison of the results and take advantage of the possibilities of this tool for complex meshes. The mesh was exported in an unstructured, single-zone CGNS file which can be used directly within SU2 without the need for conversion to the native format.
As an unstructured code, SU2 requires information about both the node locations as well as their connectivity. The connectivity description provides information about the types of elements (triangle, rectangle, tetrahedron, hexahedral, etc.) that make up the volumes in the mesh and also which nodes make up each of those elements. Lastly, the boundaries of the mesh, or markers, are given names, or tags, and their connectivity is specified in a similar manner as the interior nodes.
The configuration file is a text file that contains the user’s options for the particular problem to be solved with the SU2 suite. The configuration model used in the study can be consulted in the Appendix A.1 at the end of this work.
The software currently supports .vtk and .plt output formats natively read by ParaView and Tecplot, respectively. It was chosen to use ParaView as it provides full functionality for data visualization and is freely available under an open source license. In addition, some SU2 results are also output to comma-separated value (.csv) files, which can be read by a number of software packages.
˝ OpenFOAM
The package was downloaded and a pre-compiled version of OpenFOAM for Linux was installed using Docker Hub, which enables the binaries compiled on a given Linux environment to be run on other platforms without any performance degradation. OpenFOAM does not have a generic solver applicable to all cases. Instead, users must choose a specific solver for a class of problems to solve. The solvers within the Open-FOAM distribution are numerous and each is focused on the physical models designed to solve. For the calculations in this work, SimpleFoam was used. It is a steady-state solver for incompressible, turbulent flow, using the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm, but also includes an option to use the SIMPLEC (Semi-Implicit Method for Pressure Linked Equations Consistent) algorithm.
3.3. Two dimensional calculations
The mesh was imported from Ansys using the OpenFOAM mesh conversion utility: fluentMeshToFoam, it converts a written in ASCII format (which is not the default option in Fluent) mesh file in .msh extension to a readable by OpenFOAM mesh.
In the conversion, two dimensional geometries are currently treated by defining a mesh in three dimensions, where the front and back planes are defined as an empty boundary patch type. Therefore, when reading the two dimensional Fluent mesh, the converter automatically extrudes the mesh in the third direction and adds the empty patch, nam-ing it as "frontAndBackPlanes".
All the turbulent initial properties must be defined for the initialization of the RAS modelling, the relevant boundary conditions applied are shown in Table 3.5, where the underlined values were calculated to match the flight conditions using the following equations:
For isotropic turbulence, the turbulence kinetic energy can be estimated by:
κ= 3
2 ¨(I|ure f|))
2 (3.6)
where I is the intensity, and ure f a reference velocity. The turbulence dissipation rate follows as:
e=
C0.75µ κ1.5
L (3.7)
and the turbulence specific dissipation rate as:
ω= κ 0.5 CµL
(3.8) where Cµis a constant equal to 0.09, and L a reference length scale.
The transition momentum thickness Reynolds number can be computed as: Reθ=1173.51 ´ 589.428Tu+0.2196 Tu2 i f Tu ď 1.3 (3.9) where Tu=100 ? 2/3κ |u8| =0.02 (3.10)
Table 3.5: Turbulent property values.
Freestream velocity U8 21.911 (Re=1.5 e6) m/s
Turbulence kinetic energy κ 2.88056 e-05 m2/s2
Turbulence dissipation e 5.1124 e-07 m/s3
Specific turbulence dissipation ω 0.1972 1/s
Turbulent intensity Tu 0.02 %
Kinematic viscosity ν 1.4607347 e-05 m2/s
Eddy viscosity ratio νt/ν 10
-Transition momentum thickness Reynolds number Reθ 1.7107 e+03 m
The different configuration files containing the selected parameters and conditions for the initialization used in the study, can be consulted in the Appendix A.2 at the end of this work.
3.3. Two dimensional calculations
Finally, the post-processing was carried out in paraFoam, the main post-processing tool provided with OpenFOAM, which is simply a reader module to run with the ParaView visualization application. Nonetheless, it could have been exported to the Ansys post-processing tools by using the following converters: first foamMeshToFluent converts the OpenFOAM mesh into Fluent format and writes it out as a .msh file; and afterwards, foamDataToFluent, converts the OpenFOAM results data into a .dat file readable by Flu-ent.
˝ XFLR5
The results from the two dimensional calculations are also compared to calculations performed using the panel method code XFLR5. This tool offers a user-friendly inter-face developed from the original XFOIL tool developed by the Massachusetts Institute of Technology. This update, translates the original code written in FORTRAN to the C++ language.
Although the limitations of this software are numerous, it has proved to give reasonable and consistent results for the design of sailplanes. Consequently, a Ncrit value of 12 was used to match the turbulence level of the wind tunnel.
Convergence
A converged solution was achieved by considering a drop in accuracy to the fourth decimal in all studied equation residuals as a valid stopping criterion. Moreover, the force coefficients: Cland Cdwere monitored, imposing an additional criterion to ensure a bounded accuracy of the fifth decimal for the last 50 iterations, forcing to achieve a constant trend. Additionally, a full multi grid (FMG) initialization was implemented to provide an optimal initialization, which resulted in a faster convergence.
3.3. Two dimensional calculations
Validation
With the intention of validating the numerical simulations, a comparison was carried out, taking as reference the investigation performed in the Langley low-turbulence pressure tun-nel to determine the low-speed aerodynamic characteristics of the Wortmann FX 66-17AII-182 airfoil in [8].
Specifically, the aerodynamic coefficients were compared to the obtained data from the SST turbulent models for all the different solvers, contrasting the differences among them.
Table 3.6: Lift coefficient validation (SST).
AoA
Cl
NASA Fluent Structured Fluent Unstructured
Value Value Relative error (%) Value Relative error (%)
0 0.418 0.3875 7.30 0.3967 5.10 2 0.619 0.6002 3.04 0.6130 0.97 4 0.86 0.8072 6.14 0.8260 3.96 6 1.08 1.0010 7.31 1.0311 4.53 9 1.395 1.3503 3.21 1.2858 7.83 AoA Cl
NASA SU2 OpenFoam
Value Value Relative error (%) Value Relative error (%)
0 0.418 0.4551 8.88 0.4291 2.66
2 0.619 0.6771 9.38 0.6770 9.37
4 0.86 0.9140 6.28 0.9203 7.01
6 1.08 1.1385 5.42 1.1568 7.11
9 1.395 1.4033 0.59 1.4905 6.85
Table 3.6 show relatively similar behaviours for the lift coefficients in the almost linear slope under low AoA , nonetheless, some differences are specially visible in different regions depending of the used solver.
3.4. Three dimensional calculations
3.4
Three dimensional calculations
The performance of the Standard Cirrus glider was studied simulating three dimensional models. The CFD simulations were carried out using the Fluent γ-Reθ transition model. When the mesh elements exceeded 15 million, the simulations were run on a supercomputer [33].
Glider geometry
Two CAD models were provided by Thomas Hansen from a previous study carried out in 2014 [15].
The first one consists of the wing, wing fairing and fuselage (Figure 3.3 (a)). It is created to simulate the lift and drag coefficients on the wing and fuselage. However, a different model is required to calculate the drag coefficient of the tail section. This second model is constructed with both the fuselage and tail section and has the elevator positioned at zero degrees AoA (Figure 3.3 (b)).
(a) Wing, wing-fairing & fuselage. (b) Fuselage & tail section
Figure 3.3: CAD models used for the geometry
Mesh
The discretization of both models was created using an unstructured mesh. The number of cells in the mesh was reduced by applying symmetry conditions in order to reduce the computational cost. The domain is constructed using half a sphere and half a cylinder from which the glider models were subtracted. The domain has an upstream velocity inlet in the sphere surface 50 meter from the glider surface. Downstream of the glider, a pressure outlet was positioned 70 meter (« 12 glider´s lengths) from the glider 3.5.
The three dimensional meshes were created using the Ansys meshing tool. They used several bodies of influence around the glider, the wing-tips and the wake region. The size of the elements in the glider surface were limited to 0.008 m and the growth ratio was set to be below 1.1 to get a smooth growth. In addition, an inflation layer based on the first layer height was applied to satisfy the y+value required (0.1 < y+ < 1). It consisted of 56 layers with a 1.1 growth rate, it can be visualized in Figures 3.4 (c) and (d).
