Simulation of atmospheric effects with the DLR code TAU
Full text
TESI DE MÀSTER Màster Títol Autor Tutor Intensificació Data Master’s degree in Numerical Methods in Engineering Simulation of atmospheric effects with the DLR code TAU Laia Alcaraz Capsada Dr. Ralf Heinrich i Dr. Roberto Flores Aeron`autica, Din`amica de Fluids Computational Juliol 2015
Universitat Polit` ecnica de Catalunya Master’s Thesis Simulation of atmospheric effects with the DLR code TAU Author: Laia Alcaraz Capsada Supervisors: Dr. Ralf Heinrich Dr. Roberto Flores A thesis submitted in fulfilment of the requirements for the Master’s degree in Numerical Methods in Engineering July 2015
Declaration of Authorship I, Laia Alcaraz Capsada, declare that this thesis titled, ”Simulation of atmospheric effects with the DLR code TAU” and the work presented in it are my own. I confirm that: This work was done wholly while in candidature for a research degree at the Universitat Polit`ecnica de Catalunya. Where I have consulted the published work of others, this is always clearly attributed. Where I have quoted from the work of others, the source is always given. With the exception of such quotations, this thesis is entirely my own work. I have acknowledged all main sources of help. Where the thesis is based on work done by myself jointly with others, I have made clear exactly what was done by others and what I have contributed myself. Laia Alcaraz Capsada. Barcelona, July 2015. i
UNIVERSITAT POLIT` ECNICA DE CATALUNYA Abstract Escola T`ecnica Superior d’Enginyers de Camins, Canals i Ports de Barcelona Simulation of atmospheric effects with the DLR code TAU by Laia Alcaraz Capsada The Computational Fluid Dynamics code TAU, developed at the Deutsches Zentrum f¨ur Luft- und Raumfahrt (German Aerospace Center) is used to model encounters between atmospheric effects (wind gusts and aircraft wake vortices) and aircraft. In the case of wind gusts, two different numerical methods have been used. The first one is the so-called Disturbance Velocity Approach, a simple method that captures the influence of a gust on an aircraft, but the influence of the aerodynamics of the aircraft on the gust properties is not captured. The second method is the Resolved Atmosphere Approach, which feeds a gust into the flow field by using an unsteady boundary condition. Compared to the Disturbance Velocity Approach, this method requires more computational effort, but it is more accurate: it captures the mutual interaction of the gust and the aircraft. In order to find the validity range for the Disturbance Velocity Approach, both methods are compared for two-dimensional as well as three-dimensional geometries. In the two-dimensional cases, there is an excellent agreement between the simple method and the highly accurate one for gust wavelengths greater than or equal to two reference chord lengths and an on-flow Mach number less than or equal to 0.78. In the three-dimensional cases and a vertical gust, an unexpected agreement between both methods is found for all the gusts tested. In the three-dimensional cases and a lateral gust, no validity range for the Disturbance Velocity Approach can be given, because the resolution of the grid should be improved. Finally, the Disturbance Velocity Approach has been implemented into the TAU code to enable the simulation of wake vortices encounter problems. The correctness of the implementation has been verified through different test cases. Key words: Computational Fluid Dynamics, wind gust, wake vortex, Disturbance Velocity Approach, Resolved Atmosphere Approach. ii
Acknowledgements Above all, I would like to acknowledge Dr. Ralf Heinrich, my supervisor at the Deutsches Zentrum f¨ur Luft- und Raumfahrt (DLR), for giving me the opportunity to do my Master’s thesis there and for the economic support. I thank him for his continuous patience and dedication, for teaching me so many things and solving all my doubts. I would also like to acknowledge Dr. Roberto Flores, my supervisor at the Universitat Polit`ecnica de Catalunya (UPC), for revising the text and giving me very useful advice. I am also grateful to many people at the DLR. First of all, to Britta Ernst, who assisted me during my first days there and taught me how to use TAU. Since then, she has always been willing to help me and answer my questions, and in fact, she provided one of the unstructured grids of this thesis. In addition, I acknowledge Ana Manso Jaume, who solved me many doubts about the grid generation software. The help of Alexander Kuzolap and Markus Widhalm was very useful for some of the plots of the thesis. I also thank the former one for many advices related to Linux and LaTeX, and for the nice working atmosphere. Moreover, I acknowledge my friend David Az´ocar Paredes, for answering many doubts, giving me advice, and for all the funny moments shared together. My stance at the DLR and in Braunschweig has been very pleasant thanks to a lot of people. Among others, I acknowledge Kaˇgan At¸cı, Christoph Bosshard, Andreas Michler, Eduard Frick, Britta Roll, Ana Pal, Verena Kuhlemann and Andrei Merle. It has been very nice to have lunch with them, to learn about the German culture and to improve my German communication skills. I am also grateful to the Barcelona staff from the Centre Internacional de M`etodes Num`erics en Enginyeria (CIMNE). In particular, to Dr. Jordi Pons-Prats and Dr. Antonia Larese de Tetto, who found the opportunity at the DLR and encouraged me to do my Master’s thesis there. In addition, I thank Lelia Zielonka, from the CIMNE secretariat, who has been really helpful with many administration details not only before moving to Braunschweig, but also during my whole stance there. And last, but by no means least, I would like to thank my family, who has always supported me since the beginning of my studies, and who encouraged me to do my thesis abroad. Finally, I really acknowledge my boyfriend for his patience and his unconditional support. iii
Contents Page Declaration of Authorship i Abstract ii Acknowledgements iii Contents iv List of Figures vi List of Tables ix List of Symbols xii 1 Introduction 1 1.1 Motivation and objectives of the thesis . . . . . . . . . . . . . . . . . . . . 1 1.2 Structureofthethesis ............................. 4 I Mathematical models 5 2 The DLR code TAU 6 2.1 The mathematical description of the flow . . . . . . . . . . . . . . . . . . 6 2.1.1 Reynolds Averaging and Favre (Mass) Averaging . . . . . . . . . . 9 2.1.2 The Favre- and Reynolds-Averaged Navier-Stokes equations . . . . 9 2.2 Spatial discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Primary grid and dual grid . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Centralscheme............................. 14 2.2.3 Upwindscheme............................. 15 2.3 Temporal discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 Global time-stepping . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Dual time-stepping (DTS) . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Initial and Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.1 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.2 Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5 Acceleration techniques: Multigrid . . . . . . . . . . . . . . . . . . . . . . 22 2.5.1 Agglomeration multigrid . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5.2 Full multigrid cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 iv
Contents v 3 The Chimera technique 25 3.1 Introduction and motivation . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 The Chimera technique in the DLR-TAU code . . . . . . . . . . . . . . . 27 3.2.1 Node identification . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2.2 Interpolation at the Chimera boundaries . . . . . . . . . . . . . . . 29 4 CFD modelling of atmospheric effects 33 4.1 Characterization of atmospheric effects . . . . . . . . . . . . . . . . . . . . 33 4.1.1 Atmospheric wind gusts . . . . . . . . . . . . . . . . . . . . . . . . 33 4.1.2 Aircraft wake vortex . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.2 Disturbance Velocity Approach (DVA) . . . . . . . . . . . . . . . . . . . . 40 4.2.1 Description of the method . . . . . . . . . . . . . . . . . . . . . . . 40 4.2.2 Advantages and disadvantages . . . . . . . . . . . . . . . . . . . . 41 4.3 Resolved Atmosphere Approach (RAA) . . . . . . . . . . . . . . . . . . . 41 4.3.1 Description of the method . . . . . . . . . . . . . . . . . . . . . . . 41 4.3.2 Advantages and disadvantages . . . . . . . . . . . . . . . . . . . . 42 II Numerical applications 44 5 Grid set-up 45 5.1 Griddensitystudy ............................... 45 5.1.1 Uniformgrid .............................. 46 5.1.2 Oversetgrid............................... 47 5.1.3 Conclusions............................... 50 5.2 2Dgrids..................................... 51 5.2.1 Grids for gusts encounters . . . . . . . . . . . . . . . . . . . . . . . 51 5.2.2 Grid for wake vortex encounters . . . . . . . . . . . . . . . . . . . 53 5.3 3Dgrids..................................... 54 6 Results and discussion 58 6.1 Atmospheric wind gusts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.1.1 2Dcomputations............................ 58 6.1.2 3Dcomputations............................ 64 6.2 Interaction with wake vortices . . . . . . . . . . . . . . . . . . . . . . . . . 76 7 Conclusions 84 7.1 Summary of the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 7.2 Futurework................................... 86 A Conservation laws for moving grids 88 B Conservation laws averaging 91 Bibliography 93
List of Figures Page 2.1 Primary grid (black lines) and dual grid (red lines) around a NACA0012 airfoil....................................... 13 2.2 Unsteady computation with dual time-stepping and Runge-Kutta. . . . . 19 2.3 Example of multigrid cycle. . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4 Generation of coarse grids by the agglomeration multigrid method in 2D. 24 3.1 Grids for the Chimera technique. . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Translation of Ωcg and Ωhwith respect to Ωbg. ............... 31 3.3 NACA0012 sample configuration with movable flap. . . . . . . . . . . . . 32 4.1 Wind data in flight and height directions. Taken from [18], modified. . . . 34 4.2 ”1 −cos ” gust velocity profile (FAR Part 25.341). . . . . . . . . . . . . . 35 4.3 Encounter between a SDM configuration and a ”1 −cos ” gust. . . . . . . 35 4.4 Possible vortex encounters and their consequences. Taken from [20]. . . . 36 4.5 Scheme of the tangential velocity and definition of vortex parameters. . . 37 4.6 Comparison of the tangential velocity vtpredicted by the Lamb-Oseen and the Burnham-Hallock models for different values of the core radius rc. 38 4.7 Up: vcontour field. Down: wcontourfield.................. 39 4.8 Gust travelling relative to a NACA0012 airfoil. . . . . . . . . . . . . . . . 40 4.9 Resolution of a lateral gust with the RAA. Boundaries where unsteady boundary conditions are applied (left and bottom) are dashed. . . . . . . 42 5.1 Geodesic (xg, yg, zg), TAU (x, y, z) and aircraft fixed (xb, yb, zb) coordinate systems for t= 0 and t > 0. Taken from [13]. . . . . . . . . . . . . . . . . 46 5.2 Comparison between the analytical gust profile and the numerical one, using a uniform Cartesian grid with 100 cells per λin the direction of movement..................................... 47 5.3 Overset grid for the grid density study. Ωbg is shown in blue, while Ωtg, which is moving, is plotted in red. A vertical ”1 −cos ” gust (in yellow) moves together with Ωtg............................. 48 5.4 Comparison between the analytical gust profile and the numerical one, using three overset grids with different resolutions in the xg−direction. . 50 5.5 Unstructured grid for a NACA0012 airfoil and a HTP (in red) and rectilinear background grid (in blue). Due to the effect of the holes, the nodes of the background grid around the airfoil and the HTP have been removed. 53 5.6 Grid around a NACA0012 airfoil. . . . . . . . . . . . . . . . . . . . . . . . 54 vi
List of Figures vii 5.7 Near field grid (in blue) and surface grid (in red) for one half of the LamAiRaircraft................................. 56 5.8 Surface grid for the complete LamAiR aircraft. . . . . . . . . . . . . . . . 56 5.9 Assembled grid and hole geometry for the complete aircraft. The grid for the aircraft, Ωbg and Ωtg are shown in red, blue and green, respectively. The hole geometry is plotted in black. . . . . . . . . . . . . . . . . . . . . 57 6.1 Lift coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with Ma∞= 0.28. The reference Reynolds number is Re∞= 6.53 ·106........................... 59 6.2 Lift coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with Ma∞= 0.78. The reference Reynolds number is Re∞= 18.20 ·106. ......................... 59 6.3 Interaction of a gust with λ/Cref = 1 and the NACA0012 airfoil with both the DVA and the RAA. . . . . . . . . . . . . . . . . . . . . . . . . . 60 6.4 Interaction of a gust with λ/Cref = 1 and the HTP with both the DVA andtheRAA................................... 60 6.5 Convergence history of the residual and the lift coefficient for the steady simulation with Ma∞= 0.28 and λ/Cref =1................. 62 6.6 Convergence history of the residual and the lift coefficient for the unsteady simulation with Ma∞= 0.28, λ/Cref = 1 and the DVA. . . . . . . . . . . 62 6.7 Zoom-inoffigure6.6. ............................. 63 6.8 Lift coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with Ma∞= 0.78. The reference Reynolds number is Re∞= 24.98 ·106. ......................... 65 6.9 Details of the interaction of vertical a gust and an aircraft with both the DVA and the RAA for λ/Cref = 1 and λ/Cref =2.............. 66 6.10 Convergence history of the residual and the lift coefficient for the steady simulation with Ma∞= 0.78 and λ/Cref =1................. 67 6.11 Convergence history of the residual and the lift coefficient for the unsteady simulation with Ma∞= 0.78, λ/Cref = 1 and the RAA. . . . . . . . . . . 68 6.12Zoom-inoffigure6.11.............................. 68 6.13 Lateral force coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with λ/Cref =1............ 69 6.14 Lateral force coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with λ/Cref =2............ 70 6.15 Lateral force coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with λ/Cref =4............ 70 6.16 Slice at xg=−29 m of the y−component of the velocity field at time step 1000, for λ/Cref =4andtheRAA.................... 71 6.17 Slice at zg= 0 m of the y−component of the velocity field at time step 1900, for λ/Cref =1andtheRAA....................... 72 6.18 Slice at zg= 0 m of the y−component of the velocity field at time step 1900, for λ/Cref =2andtheRAA....................... 72 6.19 Slice at zg= 0 m of the y−component of the velocity field at time step 1900, for λ/Cref =4andtheRAA....................... 73 6.20 Slice at zg= 0 m of the y−component of the velocity field at time step 1000, for λ/Cref =4andtheRAA....................... 74
List of Symbols xiv xεfringe point ∆ygrid spacing in y−direction ∆zgrid spacing in z−direction αangle of attack αkcoefficient of the Runge-Kutta scheme (in stage k) βkcoefficients for time integration scheme (k=−1,0,1,2) γratio of specific heat coefficients at constant pressure and volume γkcoefficients for time integration scheme (k= 0,1) Γ circulation δij Kronecker delta symbol ζtype of artificial dissipation ηcoordinate of a reference finite element λwavelenght µdynamic viscosity coefficient; coordinate of a reference finite element µTeddy viscosity νkinematic viscosity coefficient (= µ/ρ) ξcoordinate of a reference finite element; yaw angle ρdensity τij components of viscous stress tensor τF ij components of Favre-averaged Reynolds stress tensor τviscous stress tensor (normal and shear stresses) τFFavre-averaged Reynolds-stress tensor φroll angle Ω control volume; domain; grid; geometry Ωisubdomain Ωij components of rotation-rate tensor dΩ volume element ∂Ω boundary of a control volume; boundary of a domain; boundary of a grid; boundary of a geometry ∇Ugradient of scalar U=h∂U ∂x ,∂U ∂y ,∂U ∂z i Subscripts baircraft fixed coordinate system bg background grid crelated to convection cg Chimera grid ggeodesic coordinate system
List of Symbols xv ip interpolated Iindex of a control volume hgrid spacing; hole-cutting geometry Lleft mindex of control volume face; species max maximum min minimum ref reference Rright tg gust-transport grid vviscous part x, y, z components in the x−, y−, z−direction ∞at infinity (farfield) Superscripts bg background grid cg Chimera grid hnumerical solution mcurrent inner iteration n−2 two previous time levels n−1 previous time level ncurrent time level n+ 1 next time level ˜ Favre averaged mean value; virtual 00 fluctuating part of Favre decomposition ¯ Reynolds averaged mean value; complement (in set theory) 0fluctuating part of Reynolds decomposition ∗dual + improved solution (in multigrid)
Chapter 1 Introduction 1.1 Motivation and objectives of the thesis The importance of the air-transportation industry has grown considerably during the last decades, and at the same time, the necessity of making it safer and more efficient. Although the security of aircraft has been increased a lot, there are situations that can still be very dangerous and can have very serious consequences. One example is the encounter of an aircraft with an atmospheric wind gust, which is a very strong discrete wind pulse, usually with a random and sudden character. Wind gusts are responsible for the appearance of additional unsteady air-loads on the aircraft, and can lead to wind shear, which is an abrupt change in wind speed and/or direction, causing turbulence or a rapid increase or decrease in velocity. Wind shear is especially delicate if it is produced during approach and landing, and for instance, it was the reason for the crash of the Delta Airlines flight 191 in 1985. Therefore, it is really of interest to predict the additional loads that arise during an encounter with a gust. This is important for the layout of the flight control system and the control surfaces, which are used to control the aircraft’s direction in flight and the aircraft flight’s attitude, respectively. In addition, gust loading has to be taken into account for the design of the structure, since for transport aircraft, gusts are the largest source of fatigue loading for the major part of the structure. The encounter between an aircraft and vortices in the atmosphere is another potentially hazardous situation, yielding for instance to an imposed roll or yaw, to additional structural loads or to a loss of altitude. Obviously, this is especially dangerous during take-off and landing: when the aircraft is close to the ground, the pilot has only a limited altitude left to recover. One of the sources for the creation of vortices lies in 1
Chapter 1. Introduction 2 the planes themselves: the wake behind an aircraft generally consists of two coherent counter-rotating swirling flows of equal strength, proportional to the lift created by the aircraft. These flows, which are like horizontal tornadoes, are called wake vortex, and are a hazard to oncoming aircraft. In fact, some accidents have been reported due to wake vortex encounters, for example the crash of the American Airlines flight 587 in 2001. To avoid such encounters, a safety distance must be maintained between two aircraft. That is, there must be a waiting time of several minutes (depending on the size of the aircraft involved and other factors) before allowing another aircraft to take-off or land on the same runway. Nevertheless, the standard separations already limit the capacity of many airports, something really problematic in view of the strong growth of worldwide air traffic. Thus, modelling and understanding an encounter between an aircraft and the wake vortices can contribute to optimize the aircraft separation distances, hence increasing airport capacities while maintaining safety levels. In the early 1970’s, and parallel to the improvements of modern aircraft, started the history of the Computational Fluid Dynamics (CFD), which can be defined as a combination of physics, numerical mathematics and computer sciences employed to simulate fluid flows. The importance of this discipline has grown very fast, and thanks to the rapidly increasing speed of supercomputers and the improvement of numerical methodologies, CFD is nowadays employed in many fields: from aircraft, turbomachinery, or car and ship design to meteorology, astrophysics or oceanography. In the field of aeronautics, CFD is seen as a complementary tool to wind tunnel and flight tests, and it has become very important for aircraft design in industry and research. To date, several powerful CFD codes have been developed. One of them is the TAU code, developed at the Deutsches Zentrum f¨ur Luft- und Raumfahrt (DLR), the national aeronautics and space research centre of the Federal Republic of Germany. TAU is a finite volume software which can simulate viscous and inviscid flows about complex geometries from the low subsonic up to the hypersonic flow regime, employing hybrid unstructured grids [22]. The development of the TAU code, which started more than a decade ago, is still being carried out at the Institut f¨ur Aerodynamik und Str¨omungstechnik (Institute of Aerodynamics and Flow Technology), in particular at its numerical methods department, called Center for Computer Applications in AeroSpace Science and Engineering (C2A2S2E). Nowadays, the code is routinely used not only at the DLR, but also at many universities and in the European aeronautical industry. In particular, the aircraft manufacturer Airbus uses the TAU code for the aerodynamic design and the prediction of aerodynamic performance of aircraft. This thesis aims to use the TAU code (release 2014.2.0) to simulate encounters between aircraft and atmospheric effects, in particular wind gusts and wake vortex.
Chapter 1. Introduction 3 On the one hand, regarding gust encounters, two different approaches had been previously implemented in TAU, whose basic two-dimensional (2D) verification was done in [12]. This work has been taken as a starting point, and one of the objectives of the thesis is to improve it for 2D geometries and extend the verification to three-dimensional (3D) cases. Such simulations are very costly and take a lot of time, so this is why another of the intentions of the thesis is to obtain accurate results with the smallest possible computational cost. For doing so, attention must be paid during the grid generation process and also when using the solver module of TAU. On the other hand, the two approaches implemented to model gust-aircraft encounters can also be adapted to model encounters with wake vortex. Thus, the last purpose of the present work is to implement one of these methods in TAU, and to verify it through some equivalent test cases. To sum up, the main objectives of the thesis are: •To analyse the two gust modelling methods implemented in TAU for 2D and for 3D geometries. •To implement a numerical method to simulate wake vortex encounters in TAU. •To work with the smallest possible computational cost, but having sufficient accuracy. As it has been already said, the numerical simulations of the present thesis are very time consuming. In this context, it is essential to perform parallel simulations on a High-Performance-Computing system. In particular, the computations of the present thesis have been run on the C2A2S2E cluster, which has 13440 cores and uses the Son of Grid Engine resource management system. The TAU code, which is parallelized based on Domain Decomposition (DD) and on the Message Passing Interface (MPI), can be efficiently run in parallel. In fact, this an important feature of the software. Before starting a CFD computation, it is important to generate suitable meshes. TAU does not have a module for grid generation, so the grids for this thesis have been created with two external grid generation software: CENTAUR [1] and the Multiblock elliptic grid generation and Computer aided system (MegaCads). The former one is a hybrid grid generation package from the company CentaurSoft, and the second one, which is property of the DLR, generates structured multiblock grids. Another very important point is the postprocessing (i.e., the visualization) of data and results. For doing so, the professional software Tecplot [2], the IsoPlot99 software (property of the DLR) and the LaTeX package PGFPlots have been used.
Chapter 1. Introduction 4 1.2 Structure of the thesis This subsection provides an overview of the structure of the present thesis, which is divided into two different parts and has two appendices. Chapters 2, 3 and 4 belong to the first part, called ’Mathematical models’, which contains the theory of the models and the numerical methods used in the thesis. The second part (chapters 5, 6 and 7) is related to the numerical simulations of the thesis. Chapter 2 is devoted to the TAU code, and explains the most relevant theory behind the software. First of all, the governing equations of the flow are defined, and some information about the turbulence models implemented is given. After that, the governing equations are discretized both in space and time, something necessary to perform CFD simulations. The section devoted to the temporal discretization is very important, since the simulations of this thesis are mainly unsteady. In the next section, the initial and boundary conditions that must be prescribed are introduced. For the sake of simplicity, details on the procedures for solving the systems of equations have been avoided. Nevertheless, the last section of the chapter is devoted to a technique used in the solution procedure, which is very important for the reduction of the computational time. Chapter 3 is related to the Chimera technique, a DD method which has been crucial for the purposes of this thesis. After motivating the usage of this technique, its implementation in the TAU code is discussed. Finally, a description about how it allows to do simulations with moving grids is given. Chapter 4 describes the two atmospheric effects considered in this work: wind gusts and wake vortex. Within the first section, both disturbances are defined and characterized in detail. After that, the first method for modelling them, called Disturbance Velocity Approach (DVA), is described. Then, the same is done with the second method, which is the so-called Resolved Atmosphere Approach (RAA). In both cases, the advantages and the disadvantages of the methods are presented. The second part of the thesis begins with chapter 5, where the grids for all the computations are described. The first section aims to find a sufficient spatial resolution to achieve good results in the simulations. Then, and taking the results of the first section into account, the grids for the 2D computations are described, both for wind gusts and for wake vortex. Finally, the 3D grids are presented. Chapter 6 describes the numerical simulations done and presents their results, which are also discussed in detail. Finally, the last chapter of the thesis presents the conclusions of the main results achieved, and suggests some work to be done in the future.
Part I Mathematical models 5
Chapter 2 The DLR code TAU This chapter is devoted to the numerical methods implemented in the DLR code TAU [22], which is a second order unstructured finite volume solver for the compressible 3D Navier-Stokes equations. The DLR code TAU (from now on simply called TAU) is not only a code, but a software system which can be used for the prediction of viscous and inviscid flows about complex geometries (mainly complex aircraft-type configurations), from the low subsonic up to the hypersonic flow regime. TAU solves the so-called Unsteady Reynolds Averaged Navier-Stokes equations, which are introduced at the beginning of this chapter. After that, these equations are discretized both in space and time, and boundary and initial conditions are introduced. Finally, some information about the multigrid technique for convergence acceleration is also given. 2.1 The mathematical description of the flow The dynamical behaviour of a fluid is determined by the conservation of mass, momentum and energy. In the absence of body forces, the conservation of these quantities is considered in the following equations ∂ ∂t ZΩ ρdΩ + I∂Ω ρ(v·n)dS = 0 ∂ ∂t ZΩ ρvdΩ + I∂Ω ρv(v·n)dS =−I∂Ω pndS +I∂Ω (τ·n)dS ∂ ∂t ZΩ ρEdΩ + I∂Ω ρE (v·n)dS =I∂Ω k(∇T·n)dS +ZΩ ˙qhdΩ−I∂Ω p(v·n)dS +I∂Ω (τ·v)·ndS (2.1) 6
Chapter 2. The DLR code TAU 7 which are called the continuity equation, the momentum equation and the energy equation, respectively. In (2.1), Ω is an arbitrary control volume and ∂Ω its closed boundary, ρis the density, vthe velocity vector, na normal vector to ∂Ω, pthe static pressure, τthe so-called viscous stress tensor, Ethe total energy per unit mass, kthe thermal conductivity coefficient, Tthe temperature and ˙qhthe heat flux. To close equations (2.1), additional relations (i.e., constitutive equations) have to be considered. One of these relations involves the viscous stress tensor τ. If its components are τij = 2µSij −2 3µ∂vk ∂xk δij (2.2) where vkis the k-th component of v,µis the dynamic viscosity coefficient, δij is the Kronecker delta symbol and the components of the strain-rate tensor Sare Sij =1 2∂vi ∂xj +∂vj ∂xi−1 3 ∂vk ∂xk (2.3) the fluid is defined as Newtonian. The following equation for the thermodynamic variables must also be considered E=e+|v|2 2(2.4) with e=cvT=h−p ρ, h =cpT(2.5) where eis the internal energy per unit mass, his the enthalpy, and cvand cpare the specific-heat coefficients for constant volume and pressure processes, respectively. They are constant because the fluid is assumed to be calorically perfect. Finally, the last relation needed to close equations (2.1) is the equation of state, which in this case, is the perfect gas law p=ρRT (2.6) where R=cp−cvis the specific gas constant. Equations (2.1), together with the constitutive equations that have been just introduced are the so-called Navier-Stokes equations in the absence of body forces. Note, however, that the Navier-Stokes equations have been introduced in a fixed frame, that is, assuming that Ω is static. For the purposes of this thesis, it is of interest to consider that the boundary of the control volume is not fixed but moves with a certain velocity vb. As stated in [13], the most general expression for vbis vb=vtrans +vrot +vflex (2.7)
Chapter 2. The DLR code TAU 8 to take into account the velocity induced by a rigid translation, a rigid rotation, or a deformation of the control volume, respectively. In the case of this thesis, vflex =vrot = 0. Equations (2.1) can be rewritten to take into account vb ∂ ∂t ZΩ ρdΩ + I∂Ω ρ[(v−vb)·n]dS = 0 ∂ ∂t ZΩ ρvdΩ + I∂Ω ρv[(v−vb)·n]dS =−I∂Ω pndS +I∂Ω (τ·n)dS ∂ ∂t ZΩ ρEdΩ + I∂Ω ρE [(v−vb)·n]dS =I∂Ω k(∇T·n)dS +ZΩ ˙qhdΩ−I∂Ω p(v·n)dS +I∂Ω (τ·v)·ndS (2.8) and are defined as the Navier-Stokes equations for an arbitrary control volume (from now on, simply called Navier-Stokes equations). As said in [5], it is possible to write these equations as ∂ ∂t ZΩ WdΩ + I∂Ω (Fc−Fv)dS =0(2.9) where Wis the vector of the conservative variables (the unknowns), Fcis the vector of convective fluxes and Fvis the vector of viscous fluxes. The definition of each of these terms can be found in the appendix A. It must be highlighted that, in a moving frame, not only (2.9) but also the so-called Geometric Conservation Law (GCL) [26] must be satisfied everywhere in the domain during the whole computation. The integral form of the GCL is ∂ ∂t ZΩ dΩ−I∂Ω vb·ndS = 0 (2.10) If this equation did not hold, the solution would be perturbed by numerical errors induced by the presence of vb. If vflex 6=0(the grid is deformed), the GCL should be discretized and solved using the same scheme applied to the Navier-Stokes equations, see [5]. For high Reynolds numbers Re, the flow usually becomes turbulent: the solution is 3D, time-dependent and also includes a large range of time- and length-scales. Since the current computer power resources are not sufficient to take all scales into account, a suitable approach is to average the Navier-Stokes equations and decompose the flow variable in mean (or averaged) values and turbulent fluctuations. It is advisable to perform the so-called Reynolds averaging for the density and pressure, while the mean values of the rest of the flow variables (for instance, the velocity v) shall be computed
Chapter 2. The DLR code TAU 15 •Matrix dissipation: Despite its simplicity, scalar dissipation may result in a lack of accuracy of the results. However, by considering a matrix (the Jacobian of the convective flux), the dissipation terms can be weighted in a proper way. As a result, the accuracy of the central scheme is improved. 2.2.3 Upwind scheme Several upwind schemes are implemented in TAU. A common feature of them is the consideration of the transport of information across the cell faces in the formulation of the convective fluxes. Compared to a central scheme with scalar artificial dissipation, these methods show an improved resolution of shocks, so they are of importance specially at the supersonic or hypersonic flow regime. Nevertheless, the numerical effort is increased, in particular in the case of unstructured grids. Several options are available in TAU, for instance the Advection Upstream Splitting Method combining flux Difference and Vector splitting (AUSMDV), the Roe scheme or the Van Leer’s scheme, either of first or second order spacial accuracy. The choice of the scheme is a decision of the user and influences the results: for example, the Van Leer’s scheme usually leads to results of lower accuracy because of its diffusive character, but if one aims to improve the convergence of the simulation, it might be a good option to perform the first iterations of the computation with this scheme [5, 15]. 2.3 Temporal discretization As said previously, the approximation of the governing equations follows the method of lines, which decouples the discretization of space and time. Thus, after dealing with spatial discretization in the previous section, (2.28) shall now be discretized in time. Before defining the time discretization techniques, two important concepts are introduced: •Global time-stepping: Consists in using the same value of the time step ∆tfor the computations in all the cells of the dual grid. On the one hand, global timestepping is easy to implement, but on the other hand, the value of ∆tis restricted by the grid geometry and the characteristics of the governing equations. In practice, this means that the scheme is only stable if the so-called Courant- Friedrichs-Lewy (CFL) condition is fulfilled, and this is a strong limitation for the computation. More information about schemes using global time-stepping is given in subsection 2.3.1.
Chapter 2. The DLR code TAU 16 •Local time-stepping: Consists in prescribing a different value of ∆tat each of the dual cells. The value of ∆tis chosen in order to enforce a stable computation at each cell. Local time-stepping can be used with explicit schemes for steady computations, where global time-stepping is extremely inefficient because of the restrictions of the CFL condition. Therefore, local time-stepping is also regarded as a convergence acceleration technique, see section 2.5. After introducing the two ways of defining ∆t, the temporal discretization of (2.28) can be done through the two schemes available in TAU: global time-stepping and dual time-stepping. 2.3.1 Global time-stepping In TAU, a Runge-Kutta multistage explicit scheme has been implemented using global time-stepping. In this case, the vector of averaged flow variables ¯ WIis computed through ¯ W(0) I=¯ Wn I ¯ W(1) I=¯ W(0) I−α1∆tI ΩI ¯ R(0) I ¯ W(2) I=¯ W(0) I−α2∆tI ΩI ¯ R(1) I . . . ¯ Wn+1 I=¯ W(m) I=¯ W(0) I−αm∆tI ΩI ¯ R(m−1) I (2.30) where αkare the stage coefficients, and ¯ R(k) Iis the residual evaluated with the solution ¯ W(k) Iof the k-th stage. Superscripts nand n+ 1 denote data at time tand at t+ ∆t, respectively, while mis the number of steps (stages) in which the solution advances in time. As it is shown in (2.30), global time-stepping explicit schemes use known data at time level nto compute the solution at the next time level, n+ 1. This makes such methods easy and fast to implement. However, as said before, its main disadvantage is that the CFL condition has to be fulfilled, which in an explicit time scheme, implies that ∆t should be equal or smaller than the time required to transport information across the stencil of the spatial discretization scheme. Otherwise, the method is not stable. The fulfilment of the CFL condition for ∆tdepends on the governing equations to be solved (Euler, Navier-Stokes...), on the dimension of the problem (one-dimensional (1D), 2D or 3D) and on the grid used. A detailed estimation of ∆tfor several cases can be found in [5].
Chapter 2. The DLR code TAU 17 2.3.2 Dual time-stepping (DTS) This method was defined with the aim of avoiding the limitations that methods with global time-stepping have due to the CFL condition. In fact, it combines local and global time-stepping, resulting in a scheme that is suitable for time-accurate computations. Hence, it has been used in the unsteady simulations of the present work. As stated in [13], the temporal discretization of (2.28) can be carried out using Backward Difference Formulas (BDF). This yields to the following implicit equation β1¯ WΩn+1 I+β0¯ WΩn I+β−1¯ WΩn−1 I+β−2¯ WΩn−2 I ∆t +γ1¯ Rn+1 I+γ0¯ Rn I=0 (2.31) where the superscript n−1 denotes the previous time level (t−∆t), and n−2 indicates two previous time levels (t−2∆t). The value of the coefficients βand γinfluences the accuracy of the BDF, which can be of first, second or third order. Table 2.1 presents the coefficients for each one of the schemes: Scheme Order β1β0β−1β−2γ1γ0 BDF1 1 1 -1 0 0 1 0 BDF2 2 3/2 -2 1/2 0 1 0 BDF3 3 11/6 -3 3/2 -1/3 1 0 Table 2.1: Coefficients for the BDF integration schemes. Due to the presence of the residual ¯ Rn+1 I=¯ R¯ Wn+1 I, (2.31) is a non-linear equation in the unknowns ¯ Wn+1 I. Hence, it has to be solved iteratively. In TAU, this is done through the so-called Dual Time Stepping (DTS) method [16], which is very popular for unsteady flows. The DTS scheme suggests solving (2.31) by defining, at each physical time step, a pseudo steady state-problem with a fictitious dual time t∗. In this state, the following equation holds dΩ¯ W∗I dt∗+¯ R∗ I=0(2.32) where ¯ W∗ Iis an approximation to ¯ Wn+1 I, the unsteady residual ¯ R∗ Iis defined as ¯ R∗ I=β1¯ W∗ IΩn+1 I+β0¯ WΩn I+β−1¯ WΩn−1 I+β−2¯ WΩn−2 I ∆t +γ1¯ R¯ W∗ I+γ0¯ Rn I (2.33)
Chapter 2. The DLR code TAU 18 and the dual time t∗satisfies lim t∗→∞ d(Ω¯ W∗)I dt∗=0 ¯ R∗ I=0 ¯ W∗ I=¯ Wn+1 I (2.34) Due to the definition of t∗, (2.31) is fulfilled when (2.32) is evaluated at the so-called pseudo-steady state (t∗→ ∞). Therefore, (2.32) is solved (that is, ¯ W∗ Iis computed) by advancing in dual time until the pseudo-steady state is reached and (2.34) holds. In that case, ¯ W∗ I=¯ Wn+1 I, i.e, the solution of (2.31) is found. The advantage of the DTS scheme is that the global time step ∆tcan be chosen based solely on flow physics and is not restricted by the CFL condition. Equation (2.32) can be solved with the multistage Runge-Kutta method. The idea is to advance in t∗, and at each dual time step ∆t∗, perform Runge-Kutta iterations (see subsection 2.3.1 for the numerical scheme) until the pseudo-steady state is reached and ¯ W∗ Iapproximates ¯ Wn+1 Iwith sufficient accuracy (in practice, this is achieved when ¯ R∗ I=¯ R∗ I,min, where ¯ R∗ I,min is the residual reduced by at least two or three orders of magnitude). Once the pseudo-steady state has been reached, the next physical time step ∆tis conducted. The computation ends when Nmax, the maximum number of physical time steps, has been reached. This procedure is shown in the flow chart of figure 2.2, where mand ncontrol the iterations in dual and physical time, respectively. The time-marching process starts with ¯ W∗ Imas an input, which can be ¯ Wn I(the vector of averaged flow variables at the previous physical time level) or can be extrapolated from previous time steps, as said in [5, 13]. Equation (2.32) can also be solved using a Lower-Upper Symmetric Gauss-Seidel (LUSGS) scheme, which is implicit and therefore usually more robust. However, the explicit Runge-Kutta multistage method combined with acceleration techniques has a great popularity. These techniques are presented in section 2.5.
Chapter 2. The DLR code TAU 19 Start ¯ W∗ Im Runge-Kutta scheme ¯ R∗ I≤¯ R∗ I,min? Advance in dual time, m=m+ 1 t∗=t∗+ ∆t∗ Advance in physical time, n=n+ 1 t=t+ ∆t Pseudo-stady state, ¯ W∗ Im=¯ Wn+1 I n=Nmax? Stop m=n= 0 yes yes no no Figure 2.2: Unsteady computation with dual time-stepping and Runge-Kutta.
Chapter 2. The DLR code TAU 20 2.4 Initial and Boundary conditions The aim of this section is to deal with the initial and boundary conditions that must be specified when solving (2.28) and computing ¯ W. 2.4.1 Boundary conditions The correct implementation of boundary conditions is one of the crucial points of the solver, since they influence the accuracy of the solution, and also the robustness and convergence speed of the numerical scheme. Two different type of boundaries can be distinguished: •Artificial boundaries. Since the computational domain is only a certain part of the (real) physical domain, the last one has to be truncated, and this results in the creation of artificial boundaries like the symmetry plane or the farfield. •Surface boundaries. These boundaries, which have a real physical meaning, have to be considered when the surface of a body is exposed to the flow. For the simulations of this thesis, some of the artificial boundaries that can be considered in TAU have been used. In the following, they are generally described (details can be found in [15]) with the exception of the so-called Chimera boundaries, to which a complete chapter is devoted (see chapter 3 for the details of the Chimera technique). Finally, the surface boundaries used are also introduced. Farfield This type of boundary condition is of importance when simulating external flows past aerodynamic configurations. The farfield boundaries have to be implemented such that there are not notable effects on the flow solution as compared to the infinite domain. In addition, any outgoing information must not be reflected back into the flow field. As said before, the hyperbolic character of the system of equations (2.28) influences how boundary conditions are applied, both in the inflow and the outflow boundaries. In the governing equations, information travels along the characteristics, and the sign of the eigenvalues of the Jacobian of the convective flux determines if it travels out or into the computational domain. According to [5], the number of variables to be imposed from outside the boundary should be equal to the number of incoming characteristics, and the remaining ones should be determined from the solution inside the domain.
Chapter 2. The DLR code TAU 21 Since all the simulations in this thesis have been performed in a subsonic regime (onflow Mach number lower than 0.8), the treatment of farfield boundary conditions will be restricted to this particular case. It turns out that, in a 3D domain, a subsonic inflow boundary has four incoming characteristics and one outgoing, while the situation for a subsonic outflow boundary is just the opposite. Therefore, four variables have to be specified for the inflow boundary, which in TAU are the three components of the velocity and the temperature [15, 28] v=v∞, T =T∞(2.35) where the subscript ∞references data at the farfield. The remaining variable, ρ, gets a value extrapolated from the first inner grid point of the domain. On the other hand, the following is prescribed at an outflow boundary ρ=ρ∞(2.36) while the values for vand Tare extrapolated from the first inner grid point of the domain. After prescribing the flow conditions at the farfield boundaries, the fluxes crossing the boundary faces have to be determined by solving a Riemann problem. Symmetry plane This type of boundary condition has to be considered when the flow is symmetrical with respect to a plane. In that case, the velocity perpendicular to the symmetric boundary must be zero, which is equivalent to enforcing zero flux across this boundary. In particular, the following shall be imposed in a symmetry plane v·n= 0,(∇ρ)·n= 0,(∇T)·n= 0 (2.37) Viscous wall This type of surface boundary, which has a physical meaning, has to be considered when a viscous fluid passes a solid wall. In this case, a no-slip boundary condition has to be implemented, that is, v=0at the surface (2.38) Moreover, the following is prescribed for the density and the temperature (∇ρ)·n= 0,(∇T)·n= 0 (2.39) When creating the primary grid, it is advisable to discretize the viscous region near the surface of the geometry (the viscous wall) with prismatic elements, which are well-suited
Chapter 2. The DLR code TAU 22 to capture boundary layers. Inviscid/Euler wall This type of surface boundary, which has a physical meaning, neglects all the viscous effects of the flow, and considers v·n= 0 at the surface,(∇ρ)·n= 0,(∇T)·n= 0 (2.40) 2.4.2 Initial conditions Before the solution process starts, the initial flow conditions have to prescribed. They determine the state of the flow at the first iteration (t= 0) in a steady (unsteady) simulation. It is important to choose appropriate initial conditions, since they can influence in how fast the final solution is found, and also in the behaviour of the numerical scheme. In the present work, the values defined at the farfield are taken as initial values in the whole flow field. In the case of an unsteady simulation, these values are initialized from a steady solution that has to be computed before. 2.5 Acceleration techniques: Multigrid As stated in [15], different techniques have been implemented in TAU to accelerate the convergence of the solution, and therefore reduce the computational time. Among them, the so-called local time-stepping (see section 2.3 for more details) and multigrid methods have been really useful for the numerical simulations. Given the dual grid of the problem, the multigrid technique consists in solving the governing equations on a series of successively coarser grids. By doing so, the following holds: •Since the control volumes are larger in a coarser grid, the value of the local time step can also be increased, resulting in an increase of the convergence speed. •The high-frequency component of the error is reduced during the computation (due to the properties of the time-stepping and numerical schemes), but this does not happen with the low-frequency components. Nevertheless, the low-frequency components of the error in the finer grid can be seen as high-frequency ones in the successively coarser grids. Therefore, they can be reduced. If this scheme is followed in a sequence of grids, the error is successively damped, and as a consequence, the convergence speed is increased.
Chapter 2. The DLR code TAU 23 2.5.1 Agglomeration multigrid Figure 2.4 shows the coarser grid levels generated from the dual grid of figure 2.1. These grids are generated through the so-called agglomeration multigrid method: given a fine dual grid, the next coarser level is generated by fusing the control volumes with their neighbours, so the dual cells of the successively coarser grids are larger and irregularly shaped. More information about the agglomeration multigrid can be found in [5]. 2.5.2 Full multigrid cycle The aim of this section is to describe how a computation with the multigrid technique is done, taking the grids of figures 2.1 and 2.4 as an example. For consistency, data in the finest grid has the subscript h(in reference to its grid spacing), while data from the coarser grids from figure 2.4 have the subscripts 2h, 4hand 8h, respectively. The so-called full multigrid cycle has been implemented in TAU to accelerate the convergence speed. It begins by computing the solution of (2.28) in one of the grid levels (h, 2h, 4hor 8h). If the computation has started in a level different from h, it is interpolated to the consecutive finer levels, until the finest one (level h) is reached. By doing that, a first solution on the finest grid is obtained efficiently. Once this has been done, the solution and the residual ¯ Wh,¯ Rhare transferred to the next coarser level by means of interpolation, and ¯ W2h,¯ R2hare computed. This procedure is repeated until the coarsest grid (level 8h) is reached. After computing ¯ W8h,¯ R8h, the so-called coarse grid correction is calculated [5]. Then, this correction is interpolated to the next finer level (4h), and thanks to it, data in that level is improved, that is, a new solution ¯ W+ 4hand a new residual ¯ R+ 4hare computed. This procedure is repeated successively until the finest grid level is reached. As a consequence, low-frequency error components of the solution and the residual are smoothed. ¯ Wh,¯ Rh ¯ W2h,¯ R2h ¯ W4h,¯ R4h ¯ W8h,¯ R8h ¯ W+ 4h,¯ R+ 4h ¯ W+ 2h,¯ R+ 2h ¯ W+ h,¯ R+ h 8h 4h 2h h Figure 2.3: Example of multigrid cycle.
Chapter 2. The DLR code TAU 24 (a) Second multigrid level, 2h. (b) Third multigrid level, 4h. (c) Fourth multigrid level, 8h. Figure 2.4: Generation of coarse grids by the agglomeration multigrid method in 2D.
Chapter 3. The Chimera technique 31 3. Given one of the cells of the box, check if the fringe point xεwhere the flow values will be interpolated is located inside the planar element to which the donor cell is related. If so, solve the following system of equations using a Newton method X i fi(ξ, η, µ)xi=xε(3.5) being xia corner of the planar element, and ξ, η, µ the coordinates of a reference finite element. The weight of the interpolation fiis 0 if xεis outside the planar element, and has a value between 0 and 1 if xεis inside the planar element. When performing an unsteady computation with the DTS scheme of (2.31), the solutions ¯ WΩn−1 Iand ¯ WΩn−2 Iof the two previous time steps are required. Nevertheless, the hole moves together with the grid to which is associated (in the present case, Ωcg), so a point in Ωbg \¯ Ωhcan become active. Against this backdrop, data at old time steps are reconstructed by interpolation of the respective values on Ωcg [21]. As an example, the following figures show two types of relative motion between grids. In the first one, there is a translation of one of the grids, while in the second one, one of the grids is rotated such that the flap of a NACA0012 airfoil is deflected. Note that, in both cases, the hole associated to the moving grid moves together with it. (a) Initial configuration, t= 0. (b) Final configuration, t > 0. Figure 3.2: Translation of Ωcg and Ωhwith respect to Ωbg.
Chapter 3. The Chimera technique 32 (a) Initial configuration, t= 0. (b) Final configuration (t > 0), flap deflected 20◦. Figure 3.3: NACA0012 sample configuration with movable flap.
Chapter 4 CFD modelling of atmospheric effects This chapter is devoted to the numerical methods implemented in TAU to enable encounters between aircraft and atmospheric effects (also called disturbances). In the first section, the particular atmospheric turbulences considered in this thesis are defined and characterized. After that, the second and the third sections present the two numerical methods used in the simulations of chapter 6: the Disturbance Velocity Approach and the Resolved Atmosphere Approach. 4.1 Characterization of atmospheric effects The two examples of atmospheric effects considered in the present work are atmospheric wind gusts and aircraft wake vortex. 4.1.1 Atmospheric wind gusts Among others, pressure differences in the air-layer around the Earth, together with solar radiation, the Coriolis effect or the presence of oceans and mountains, lead to a global air circulation in the atmosphere which changes in time and in space (see figure 4.1). Atmospheric motion can be classified into different scales, since it comprises phenomena of different magnitudes. If the smallest time-scales are considered (minutes and seconds), the variations in wind speed and direction are caused by wind gusts and turbulence [27]. A wind gust (from now on, simply called gust) is a strong, discrete wind pulse. Caused by a combination of thermal effects and friction with the Earth’s surface, the suddenness 33
Chapter 4. CFD modelling of atmospheric effects 34 and the random character of a gust can lead to very dangerous situations when interacting with an aircraft: gusts are responsible for the appearance of additional unsteady air-loads, which can have serious consequences for the aircraft structure, stability and control. x[m] z[m] w[m/s] (a) z−component of the velocity field. x[m] w[m/s] (b) z−component of the velocity field obtained at a height of 23 m, which corresponds to the lowest horizontal black line in figure 4.1(a). The single gusts obtained are marked in red. Crosses indicate the starting and ending points of each gust. Figure 4.1: Wind data in flight and height directions. Taken from [18], modified. TAU can simulate encounters with vertical and lateral gusts, whose shape can correspond to a sharp edge gust (modelled by a Heaviside step function), to the one considered in [18], or to the ”1 −cos ” law of the Federal Aviation Regulations (FAR) Part 25.341 [8]. For example, a ”1 −cos ” vertical gust is defined as w(x) = 1 2w01−cos 2πx λ (4.1) where wis the z−velocity component (in Cartesian coordinates), w0is the amplitude of the gust and λis its wavelenght.
Chapter 4. CFD modelling of atmospheric effects 35 −0.5−0.25 0 0.25 0.5 0.75 1 1.25 1.5 −0.2 0 0.2 0.4 0.6 0.8 1 λ w0 x/λ [-] w/w0[-] Figure 4.2: ”1 −cos ” gust velocity profile (FAR Part 25.341). The ”1 −cos ” gust shape is mandatory for the aircraft certification process, and has been used in aircraft design processes and loads calculations in the last 60 years [18]. Thus, it is the one considered in this thesis, where gusts move relative to the aircraft at a translational velocity in the x−direction u∞=a∞Ma∞=rγp∞ ρ∞ Ma∞(4.2) being athe speed of sound, Ma the Mach number and γ= 1.4. Figure 4.3 shows two examples of gust encounters with the Standard Dynamic Model (SDM) configuration [29], which is a generic fighter configuration similar to an F16. (a) Encounter with a vertical gust. (b) Encounter with a lateral gust. Figure 4.3: Encounter between a SDM configuration and a ”1 −cos ” gust.
Chapter 4. CFD modelling of atmospheric effects 36 4.1.2 Aircraft wake vortex Aircraft wake vortex are the second type of atmospheric effect considered in this thesis. As stated in [10], the wake flow behind an aircraft is a natural consequence of the existence of aerodynamic lift. On the one hand, when an aircraft is flying, an upward motion, called ”upwash”, prevails in the region beyond both wing tips. On the other hand, there is a stronger downward motion just behind the trailing edge of each of the wings, which is known as ”downwash”. These two motions are shown in the following figure: Figure 4.4: Possible vortex encounters and their consequences. Taken from [20]. In the near field of the wing, small vortices emerge from the so-called vortex sheet at the wing tips and at the edges of the landing flaps [10], in a process that is governed by several physical phenomena: boundary-layer separation, initiation of vortex stabilities, etc. The vortex sheet experiences a roll-up and the co-rotating vortices merge, resulting in a wake behind an aircraft which consists of two coherent counter-rotating swirling flows: the aircraft wake vortex. Each of these flows, which are like horizontal tornadoes of about equal strength, has its own strong circulation, which is proportional to the weight of the aircraft. Once the aircraft has passed, the two wake vortex dissipate slowly in the atmosphere. Nevertheless, the wingtip vortices cause the so-called wake turbulence, which is hazardous to other aircraft. The situation is specially delicate during take-off and landing, where strong wake vortex are created because aircraft fly at high Angle of Attack (AoA). In this case, an encounter between a light aircraft and the wake turbulence of a heavier one can be really dangerous, see figure 4.4 for some of the consequences.
Chapter 4. CFD modelling of atmospheric effects 37 Consider a system of coordinates with origin in the centre of gravity of an aircraft, and x,yand zdefining the axes of flight, span and height, respectively. In this context, several parameters can be defined to describe the wake vortex behind an aircraft [10]: •Tangential velocity vt=√v2+w2, which is the most obvious feature of a swirling flow. •Circulation Γ, related to the amount of vortex strength Γ = Γ (r) = Ic vtds (4.3) where cis a closed curve, and r=py2+z2is the distance from the vortex centre. •Core radius rc, which, once the roll-up is completed, defines the value of rcorresponding to the maximum tangential velocity vt,max. •Distance between both vortices b, corresponding to a 70% or 75% of the wingspan of the preceding aircraft. b rcrc left vortex right vortex y vt Figure 4.5: Scheme of the tangential velocity and definition of vortex parameters. The tangential velocity vtis modelled independently for each of the vortices through an analytical model. In the present work, the following ones have been considered: •Lamb-Oseen model [19]: vt=Γ 2πr h1−exp −1.2544(r/rc)2i (4.4)
Chapter 4. CFD modelling of atmospheric effects 38 •Burnham-Hallock model [7]: vt=Γ 2π r r2 c+r2(4.5) Figure 4.6 shows a comparison between both models, considering only one vortex. The plots have been done with Γ = 158 m2/s and b= 21.5 m, corresponding to a VFW- Fokker 614 jetliner, see [9]. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 2 4 6 8 10 12 14 16 18 r/b [-] vt[m/s] Lamb-Oseen, rc= 1.075 m Burnham-Hallock, rc= 1.075 m Burnham-Hallock, rc= 0.75 m Figure 4.6: Comparison of the tangential velocity vtpredicted by the Lamb-Oseen and the Burnham-Hallock models for different values of the core radius rc. It is clear that, for the same core radius rc, the value of the maximum tangential velocity vt,max is different. The opposite holds for its position, rmax/b. This is detailed in the table below: Model rc[m] vt,max [m/s] rmax/b [-] Lamb-Oseen 1.075 16.72 0.050 Burnham-Hallock 1.075 11.70 0.050 Burnham-Hallock 0.75 16.76 0.035 Table 4.1: Maximum value of the tangential velocity profiles of figure 4.6. For the same core radius rc, the value of vt,max predicted by the Lamb-Oseen model is around 40 % greater than the value of vt,max predicted by the Burnham-Hallock model. However, the position of the maximum is the same, rmax/b = 0.050. The maximum value of vt,max predicted by the Burnham-Hallock model with rc= 0.75 m is really similar to the Lamb-Oseen one, but rmax/b diminishes.
Chapter 4. CFD modelling of atmospheric effects 39 The wake vortex are modelled by a superposition of two ideal vortices with an opposite circulation (the left vortex rotates clockwise, while the right one rotates counterclockwise), see figures 4.4 and 4.5. Defining (xvortex,1, yvortex,1, zvortex,1) and (xvortex,2, yvortex,2, zvortex,2) as the Cartesian coordinates from the left vortex and the right vortex, respectively, the following radii can be computed [9] r1=q(y−yvortex,1)2+ (z−zvortex,1)2(4.6) r2=q(y−yvortex,2)2+ (z−zvortex,2)2(4.7) Then, the superposition of the left and the right vortices yields to a velocity field that, in Cartesian coordinates, has the following components in span and height direction v=vt(r1)z−zvortex,1 r1−vt(r2)z−zvortex,2 r2(4.8) w=−vt(r1)y−yvortex,1 r1+vt(r2)y−yvortex,2 r2(4.9) where, as it has been said before, the tangential velocity vtis given by (4.4) or (4.5). Plots of vand w, using Γ = 300 m2/s, b= 21.5 m and rc= 0.75 m are shown in figure 4.7 for the Burhnam-Hallock model. X Y Z w [m/s] 22 18 14 10 6 2 -2 -6 -10 -14 -18 -22 X Y Z v [m/s] 22 18 14 10 6 2 -2 -6 -10 -14 -18 -22 Figure 4.7: Up: vcontour field. Down: wcontour field.
Chapter 4. CFD modelling of atmospheric effects 40 4.2 Disturbance Velocity Approach (DVA) The Disturbance Velocity Approach (DVA) is the first method used to model encounters between aircraft and atmospheric effects (gusts and wake vortex). 4.2.1 Description of the method The approach followed to implement the DVA is very similar to the one done in section 2.1 when introducing the velocity vbof a control volume Ω. In (2.8), vbslightly alters the flux balance, yielding, for instance, to the following continuity equation ∂ ∂t ZΩ ρdΩ + I∂Ω ρ[(v−vb)·n]dS = 0 (4.10) where the convection across the cell interface of the control volume is altered from v(in the original continuity equation (2.1)) to v−vb. This balance also changes if the velocity induced by an atmospheric effect viis taken into account. Then, (4.10) becomes [12] ∂ ∂t ZΩ ρdΩ + I∂Ω ρ[(v−vb−vi)·n]dS = 0 (4.11) The same can be done for the momentum and the energy equations of (2.8). The induced velocity viis a function of space and time, and in TAU it can be prescribed either for vertical and lateral gusts or wake vortex. In this last case, the normal vector to the plane containing the vortices can have any orientation in space. Figure 4.8 shows a vertical gust of wavelength λthat moves relative to a NACA0012 airfoil with a translational speed u∞given by (4.2). If the gust is still far from the leading edge, the boundary of the control volume (in gray) has vi=0. However, when the gust is just beneath the airfoil, the induced velocity in the z−direction is equal to the amplitude of the gust. As said in [13], the local effect of the gust is the same, as if the airfoil is moving with a velocity (0,0,−w0) downward. u∞ vi=0 λ z x (a) Initial situation (t= 0). u∞ vi= (0,0,−w0) λ z x (b) Gust beneath the airfoil (t > 0). Figure 4.8: Gust travelling relative to a NACA0012 airfoil.
Chapter 5. Grid set-up 47 −0.500.511.52 0 0.02 0.04 0.06 0.08 0.1 xg[m] w/u∞[-] analytical 100 cells Figure 5.2: Comparison between the analytical gust profile and the numerical one, using a uniform Cartesian grid with 100 cells per λin the direction of movement. The error observed is measured using an error defined in the L2norm eL2=sZΩ (wh−w)2dΩ.ZΩ w2dΩ (5.1) where Ω is the computational domain, whis the numerical data, and wis the analytical gust profile, given by (4.1). In the case of figure 5.2, the error is eL2= 15.06%. Therefore, an uniform Cartesian grid with ∆xg= 0.01 m (100 cells per λ) yields to an unacceptably large error. The value of eL2can be reduced by refining the grid in the direction of the gust’s movement. However, this implies much more expensive simulations, which is especially critical for the 3D cases. The alternative approach suggested in [12] provides much better results within an acceptable computational time, as it is explained in the following section. 5.1.2 Overset grid Motivated by the results of figure 5.2, a new density study is made. In this case, the primary grid is assembled from the following two component grids: •A background grid, Ωbg. This grid is rectilinear, i.e., its spacings in both xg−and zg−directions are not constant. This is done in order to reduce the number of nodes.
Chapter 5. Grid set-up 48 •A Chimera grid, which, following the nomenclature of [12], is called gust-transport grid, Ωtg. This Cartesian grid, whose length in xg−direction is 2λand overlaps with Ωbg, is used to transport the gust through the domain. At the beginning of the computation (t= 0), the gust is just in front of the computational domain. For t > 0, the gust is fed into the flow field through unsteady farfield boundary conditions at the left and bottom boundaries. After that, the gust advances in time, moving to the right of the domain. Once the gust reaches the center of Ωtg, the gust-transport grid starts to move, together with the gust, with the convection velocity of the flow. In the moving grid, the velocity of the boundary of a control volume is vb= (−u∞,0,0) (in geodesic coordinates). Therefore, equations (2.8), presented in chapter 2, have to be taken into account. To create such an overset grid and enable the movement of Ωtg, it is necessary to use the Chimera technique, which is described in detail in chapter 3. xg[m] zg[m] Figure 5.3: Overset grid for the grid density study. Ωbg is shown in blue, while Ωtg, which is moving, is plotted in red. A vertical ”1 −cos ” gust (in yellow) moves together with Ωtg. Once the gust is centred in Ωtg, it is placed in the same nodes (the ones belonging to Ωtg) until the end of the computation. This implies that ∂w ∂t = 0 at these nodes, and the properties of the gust are expected to be conserved. Note that this is contrary to what happens for the single grid of subsection 5.1.1. In that case, ∂w ∂t 6= 0 at the nodes of the
Chapter 5. Grid set-up 49 grid, because the gust moves through the domain, so the properties of the gust are not conserved (see figure 5.2). As it has been said in chapter 3, the Chimera technique needs a hole-cutting geometry, Ωh, to define the Chimera boundaries where communication between the grids is done. Furthermore, the governing equations are solved twice in the region where both grids overlap, so it is usually more expensive compared to standard computations. In the grid density study of [12], three Cartesian component grids are used: Ωbg, Ωtg, and a fine grid at the left part of the domain, used to feed the gust into the domain with sufficient resolution. With the aim of reducing the computational cost of using three component grids, only Ωbg and Ωtg have been used in the present work, but Ωbg is rectilinear instead of Cartesian (see figure 5.3). Hence, Ωbg has enough spatial resolution to feed in the gust successfully, and once the gust is centred in Ωtg, it becomes coarser in order to reduce the number of nodes (because the gust is already placed in Ωtg, not in Ωbg). This type of grid has been created with the Multiblock elliptic grid generation and Computer aided system (MegaCads), which is a grid generation software property of the DLR suitable for creating block-structured meshes. The grid density study has been done with the overset grid just described for three different gust resolutions: 25 cells per λ, 50 cells per λand 100 cells per λ, where λ= 1 m. Hence, several grids Ωbg and Ωtg have to be generated for each one of the resolutions. Some details of the grids are given in the table below: Resolution ∆xgfor Ωtg [m] Nodes (assembled grid) 25 cells per λ0.04 147204 50 cells per λ0.02 161304 100 cells per λ0.01 189504 Table 5.1: Overset grids for the density study. The simulation described in subsection 5.1.1 has been done with the three overset grids, and their results are compared to the analytical gust profile in figure 5.4. It can be seen, for a resolution of 25 cells per λ, that the maximum value of w/u∞is lower than the analytical one. Besides, deviations from the analytical solution are specially visible upstream of the gust. As expected, higher resolutions offer better results, as it is shown by the values of eL2in table 5.2. Although a resolution of 50 cells per λyields very good results, using 100 cells per λis the best option: the amplitude of the numerical gust is conserved during the transportation, and only minor deviations are visible upstream and downstream of the gust, so the value of eL2is acceptable.
Chapter 5. Grid set-up 50 −0.500.511.5 0 0.02 0.04 0.06 0.08 0.1 xg[m] w/u∞[-] analytical 100 cells 50 cells 25 cells Figure 5.4: Comparison between the analytical gust profile and the numerical one, using three overset grids with different resolutions in the xg−direction. Resolution eL2[%] 25 cells per λ5.31 50 cells per λ1.52 100 cells per λ1.39 Table 5.2: Error in the L2norm for the overset grid. 5.1.3 Conclusions With the aim of finding a sufficient grid resolution in the xg−direction for both 2D and 3D simulations, two different grid density studies have been performed. The first one considers a Cartesian uniform grid with a resolution of 100 cells per λ, while for the second one, an overset grid has been used. The overset grid consists of a background grid and a gust-transport grid, which are assembled using the Chimera technique. Figures 5.2 and 5.4 compare the numerical gust and the analytical one for the two types of grids described. It is clear that the overset grid offers much better results than the uniform grid, where the error in the L2norm is high (eL2= 15.06 %). This is because: •∂w ∂t = 0 in the case of the overset grid. Once the gust is centred in Ωtg and its movement starts, the gust is always placed in the same nodes (the ones belonging to Ωtg), so the properties of the gust are conserved. •∂w ∂t 6= 0 in the case of an uniform grid. Thus, the properties of the gust are not conserved as the gust moves through the domain.
Chapter 5. Grid set-up 51 In the case of an uniform grid, the properties of the gust can be conserved if ∆xgis sufficiently small. Unfortunately, this implies very expensive computations, so using an overset grid is more advisable. In this case, the gust is transported with higher fidelity, and even with a resolution of 25 cells per λ, the value of eL2is clearly better than the one from an uniform grid, for a significant lower computational cost. Motivated by the results obtained in the grid density studies, the grids used for gust simulations, which are described in the following sections, contain a moving grid with a resolution of 100 cells per λin the direction of the gust movement. 5.2 2D grids In this section, all details related to the 2D grids are given. In the first subsection, the grids for the simulations involving atmospheric gusts are described. The same is done in the last subsection for the grid to simulate wake vortex encounters. 5.2.1 Grids for gusts encounters The gust simulations of chapter 6 aim to set the range of validity of the DVA by comparing its results with the RAA ones. One of the critical points of the RAA is the requirement of a high spatial resolution to transport the gust without too much numerical diffusion. According to the grid density study of the previous section, an overset grid with a resolution of 100 cells per λis sufficient. A lower resolution might be sufficient for the DVA, but for the sake of consistency, the same grid has to be used to compare it to the RAA. The DVA and the RAA are compared using a configuration consisting in a symmetrical NACA0012 airfoil and a Horizontal Tail Plane (HTP), also called horizontal stabilizer. Both airfoils are at α= 0◦, where αis the AoA. The reason for choosing such a configuration is the expectation, that the aerodynamics of the airfoil will influence the properties of the gust, which afterwards interacts with the HTP [12]. As said in chapter 4, the DVA cannot capture the effect of the aerodynamics of the airfoil on the gust, but this can be done with the RAA. Hence, prediction errors for the DVA are expected, especially for gusts of short wavelength. In order to check this, an encounter between a ”1 −cos ” vertical gust and the NACA0012-HTP configuration has to be computed for different wavelengths.
Chapter 5. Grid set-up 52 The primary grid is created with the Chimera technique, and it is assembled from the following component grids: •A rectilinear background grid Ωbg =xg×zg∈[−20,20] ×[−20,20] m2, which has a growing value of ∆xg. The finest zone (∆xg=λ/100) has a length of λ/2 in the xg−direction, and, as explained in subsection 5.1.2, is used to feed the gust into the domain with sufficient spatial resolution. After that, and within a region of length 2λ, the grid spacing grows, according to a Poisson distribution, until ∆xg= 0.1 m. With the aim of saving nodes and reducing the computational cost, the grid is also rectilinear in the direction of height, with a higher resolution for zg∈[−3,3] m (close to the airfoils). This grid has been created with MegaCads. •A Cartesian gust-transport grid Ωtg, which has a size of 2λ×40 and a constant spacing ∆xg=λ/100. As said in subsection 5.1.2, this grid, which has been created with MegaCads, moves during the computation thanks to the Chimera technique. •An unstructured mesh containing the near field of the NACA0012 airfoil and the HTP, with a size of [−5.5,2]×[−2,2] m2. The distance between the inflow farfield boundary and the airfoil is 20 Cref , being Cref = 1 m the reference chord length of the airfoil. This grid has been created with the grid generation software CENTAUR [1], and has prismatic elements close to the viscous wall to properly capture the boundary layer of the flow. In order to check different gust wavelengths, simulations with λ/Cref = 1, λ/Cref = 2 and λ/Cref = 4 are made. Since many features of Ωbg and Ωtg depend on λ, different grids have to be created for each value of λ. Table 5.3 presents some details of the different Ωtg grids. λ/Cref = 1 [-] λ/Cref = 2 [-] λ/Cref = 4 [-] ∆xg[m] 0.01 0.02 0.04 Size [m2] 2 ×40 4 ×40 8 ×40 Table 5.3: Features of Ωtg depending on λ/Cref . When using the Chimera technique, the hole-cutting geometries remove some nodes of the grids and create the Chimera boundaries, see chapter 3. In the present case, there are three 2D holes: •The first hole is associated to Ωtg, the gust-transport grid (with movement). During the simulation, this hole removes some nodes of Ωbg. Therefore, before reaching
Chapter 5. Grid set-up 53 the airfoil, the gust is placed within the nodes of Ωtg, which has a sufficient spatial resolution. Since this hole does not cut the near field grid, its nodes are not removed, but this does not affect the accuracy of the results (the near field grid has sufficient spatial resolution). •The second hole is associated to the NACA0012 airfoil, and removes nodes of Ωbg and Ωtg. Thus, the near field of the airfoil is always discretized with the unstructured mesh that has also elements for resolving the boundary layer. •The third hole is associated to the HTP, and removes nodes of Ωbg and Ωtg. It works like the hole for the NACA0012 airfoil. In TAU, the user can select in which grids nodes should be removed, and also the grid, to which a hole is associated. This is specially of importance in unsteady computations with moving grids, because the hole can move together with the grid to which it is associated. In the present case, this is done for the hole associated to Ωtg. xg[m] zg[m] Figure 5.5: Unstructured grid for a NACA0012 airfoil and a HTP (in red) and rectilinear background grid (in blue). Due to the effect of the holes, the nodes of the background grid around the airfoil and the HTP have been removed. 5.2.2 Grid for wake vortex encounters One of the objectives of this thesis is to implement the DVA to enable the simulation of encounters with wake vortices. The implementation has been done in the most general
Chapter 5. Grid set-up 54 way, allowing the normal vector to the plane containing the vortices to have any orientation in space. The code has to be verified through some equivalent test cases (see section 6.2 for more details), which, for the sake of simplicity, are inviscid and consider a NACA0012 airfoil with α= 0◦. The primary grid for the airfoil has been created with CENTAUR and it is shown in figure 5.6. Since the simulations are inviscid, the inviscid wall boundary condition has to be prescribed (see subsection 2.4.1). Thus, the grid does not have any prismatic elements near the wall. x[m] z[m] Figure 5.6: Grid around a NACA0012 airfoil. 5.3 3D grids In subsection 5.2.1, a configuration involving a NACA0012 airfoil and an HTP is suggested for comparing the DVA and RAA. Nevertheless, more realistic results are expected with a 3D simulation involving a complete aircraft, where before reaching the tail of the aircraft, the gust interacts with the fuselage, and also with its wings. As a consequence, more discrepancies between the DVA and RAA may arise. In this case, the following simulations are considered: •An encounter between an aircraft (α= 0◦) and a ”1 −cos ” vertical gust with λ/Cref = 1, λ/Cref = 2 and λ/Cref = 4, where Cref = 4 m. Such simulations have a clear advantage: there is a symmetry with respect to the yg= 0 m plane
Chapter 5. Grid set-up 55 (where ygis the span axis in geodesic coordinates), so only one half of the aircraft has to be considered. This reduces the number of mesh nodes and hence the computational cost of the simulation by a factor of 2. •An encounter between an aircraft (α= 0◦) and a ”1 −cos ” lateral gust with λ/Cref = 1, λ/Cref = 2 and λ/Cref = 4. In this case, the situation is not symmetric with respect to the yg= 0 plane, so the complete aircraft must be considered. When creating the 3D primary grids, the conclusions of the 2D grid density study done in section 5.1 can be also considered to be valid. Thus, the approach followed for the 2D grids in subsection 5.2.1 (that is, overset grid and Ωtg with a resolution of 100 cells per λ) is also used for the 3D applications. The primary grid is assembled from: •A rectilinear background grid Ωbg =xg×yg×zg∈[−40,70] ×[−40,40] ×[−40,45] m3for the half configuration. In the case of the complete configuration, yg∈ [−80,80] m. The resolution in span and height direction is increased close to the near field of the aircraft, according to a Poisson distribution, and the finest spacings in span and height direction are ∆yg= ∆zg= 0.25 m. As in the 2D case, the spacing also changes in the direction of flight in order to feed in the gust properly: the finest spacing is equal to λ/100, while the coarsest one is equal to 0.25 m. •A Cartesian gust-transport grid Ωtg =xg×yg×zg∈[0,2λ]×[−40,40]×[−40,45] m3 for the half configuration. In the case of the complete configuration, yg∈[−80,80] m. This grid has a spacing ∆xg=λ/100, and moves with the convection velocity u∞of the gust to transport it through the domain. •An unstructured grid for the surface and the near field of the aircraft. The aircraft chosen is the one from the DLR project Laminar Aircraft Research (LamAiR), see [24] for an overview of the work. The grid has been created with CENTAUR, and has prismatic elements close to the viscous wall to capture the boundary layer of the flow. Figure 5.7 shows both the surface and the near field grids for one half of the aircraft, while figure 5.8 shows the complete aircraft with its surface grid. As in the 2D case, the features of the grids depend on λ, so several Ωbg and Ωtg grids have to be created. The number of nodes for the final assembled grid is shown in table 5.4 for each one of the cases. As it can be seen, the number of nodes is much higher compared to the 2D cases, in particular for the complete aircraft. Thus, these simulations are the most expensive ones that have been considered in this thesis.
Chapter 5. Grid set-up 56 λ/Cref = 1 [-] λ/Cref = 2 [-] λ/Cref = 4 [-] Half aircraft 24.68 ·10624.45 ·10623.49 ·106 Complete aircraft 49.16 ·10651.72 ·10650.21 ·106 Table 5.4: Number of nodes of the 3D assembled grids. Figure 5.7: Near field grid (in blue) and surface grid (in red) for one half of the LamAiR aircraft. Figure 5.8: Surface grid for the complete LamAiR aircraft.
Chapter 6. Numerical results 63 38000 38100 38200 38300 38400 38500 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 0.01 0.02 0.03 0.04 0.05 MG-cycles Residual [-] CL[-] Physical time step Pseudo-steady state Figure 6.7: Zoom-in of figure 6.6. As it has been already said, figure 6.5 shows the convergence history of the steady simulation, used to find the initial flow conditions for the following unsteady computation. A central scheme with matrix dissipation and local time-stepping with a Runge-Kutta scheme has been used. Notice that, after 1500 MG-cycles, the residual has decreased some orders of magnitude (it is of order 10−7), and CLhas a constant value (the lift is stable). This is the criterion considered to ensure that the computation has converged. The basic numerical parameters of the unsteady simulation are the same as in the steady one. Nevertheless, the DTS scheme has been used to allow a time accurate flow solution. The physical time step size ∆tis computed as ∆t=tref 100 =Cref 100u∞ (6.3) In particular, 4000 time steps are necessary for the simulation of figure 6.6. As the flow chart of figure 2.2 shows, several iterations in dual time (now MG-cycles, because the Multigrid technique is used) are performed within each physical time step. As described in chapter 2, the solution in dual time has to converge, within each physical time step, to the so called pseudo-steady state, which is reached, in practice, if the residual of the computation decreases some orders of magnitude and CLhas an horizontal tangent. A zoom-in of figure 6.6 is presented in figure 6.7. The plot clearly shows the behaviour of the residual and CLas the simulation advances in physical time. It can also be seen that about 50 MG-cycles in dual time are necessary to reach the pseudo-steady state,
Chapter 6. Numerical results 64 where the residual has decreased some orders of magnitude and CLhas a constant value (horizontal tangent). In order to reduce the computational effort, the number of MG-cycles that are performed within each time step changes depending on the requirements of the computation. In particular, it changes from a minimum of 5 to a maximum of 50 (shown in figure 6.7). For example, 50 MG-cycles are necessary to reach the pseudo-steady state when the gust is interacting with the airfoil, but when it is still far away from its leading edge, 5 MG- cycles are sufficient. The number of MG-cycles also depends on the numerical method chosen to model the gust: compared to the DVA, the RAA needs more MG-cycles at the beginning of the computation, to feed the gust into the domain with sufficient resolution. Therefore, the RAA has greater computational cost compared to the DVA. 6.1.2 3D computations If a gust interacts with a realistic 3D configuration including the fuselage and the tail of an aircraft, more discrepancies between the DVA and the RAA may arise compared to the 2D case: the gust interacts three times with the aerodynamics of the aircraft instead of two. As a consequence, the range of validity of the DVA set in 2D may change. In the following, encounters between a vertical and a lateral gust with λ/Cref = 1, λ/Cref = 2 and λ/Cref = 4 with the LamAiR aircraft (with Cref = 4 m, presented in section 5.3) are simulated with both the DVA and the RAA. The vertical gust influences the aerodynamics of the fuselage, the wings and horizontal tail, whereas the lateral one mainly influences the fuselage and the vertical tail. The computations done have a large computational effort, specially in the case of a lateral gust, since a grid for the complete aircraft has to be used. Hence, to avoid many expensive calculus, only simulations with Ma∞= 0.78 have been done, because it corresponds to the biggest discrepancies between the DVA and the RAA in 2D (see table 6.1). 6.1.2.1 Interaction with a vertical gust Figure 6.8 shows the history of the lift coefficient CLas a function of the dimensionless time t/tref . At the beginning of the simulation, the value of CLis constant, because the gust is still far away from the aircraft. Then, the gust interacts with the fuselage of the aircraft, as the first growth in CLshows, and after that, there is an interaction with the wing, yielding to the major growth of the lift. Finally, the gust interacts with the tail of the aircraft, which corresponds to the last growth of the lift.
Chapter 6. Numerical results 65 7 9 11 13 15 17 19 21 23 25 27 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 t/tref [-] CL[-] λ/Cref = 1, RAA λ/Cref = 1, DVA λ/Cref = 2, RAA λ/Cref = 2, DVA λ/Cref = 4, RAA λ/Cref = 4, DVA Figure 6.8: Lift coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with Ma∞= 0.78. The reference Reynolds number is Re∞= 24.98 ·106. As in the 2D case, the maximum CLcan be related to the maximum load, and it is used to evaluate the discrepancies between the DVA and the RAA through eCL,max, defined in (6.2). Table 6.2 shows the value of eCL,max corresponding to data of figure 6.8: λ/Cref [-] eCL,max [%] 1 0.90 2 1.02 4 0.41 Table 6.2: eCL,max values for the 3D computations involving a vertical gust. Compared to the errors of table 6.1, the ones from table 6.2 are surprisingly small. There is a very good agreement between the DVA and the RAA in terms of the maximum load, independently of the wavelenght of the wind gust. Presumably, there are lift cancellations due to the global character of this magnitude, so the values of eCL,max are smaller than expected. Lift cancellations are also assumed to be the reason why eCL,max for λ/Cref = 1 is even smaller than eCL,max for λ/Cref = 2, which is unexpected and contrary to the results of the 2D simulations. Figure 6.9 shows zoom-ins of figure 6.8 for both λ/Cref = 1 and λ/Cref = 2, the two wavelengths that present more differences between the DVA and the RAA.
Chapter 6. Numerical results 66 10 10.5 11 11.5 12 12.5 13 0.32 0.33 0.34 0.35 t/tref [-] CL[-] RAA DVA (a) Interaction with the fuselage, λ/Cref = 1. 10 10.5 11 11.5 12 12.5 13 0.33 0.34 0.35 0.36 t/tref [-] CL[-] RAA DVA (b) Interaction with the fuselage, λ/Cref = 2. 14 15 16 17 18 19 0.3 0.35 0.4 0.45 0.5 t/tref [-] CL[-] RAA DVA (c) Interaction with the wing, λ/Cref = 1. 14 15 16 17 18 19 0.4 0.45 0.5 0.55 0.6 t/tref [-] CL[-] RAA DVA (d) Interaction with the wing, λ/Cref = 2. 19 19.5 20 20.5 21 21.5 22 0.34 0.35 0.36 0.37 0.38 0.39 0.4 t/tref [-] CL[-] RAA DVA (e) Interaction with the tail, λ/Cref = 1. 19 19.5 20 20.5 21 21.5 22 0.37 0.38 0.39 0.4 0.41 0.42 0.43 t/tref [-] CL[-] RAA DVA (f) Interaction with the tail, λ/Cref = 2. Figure 6.9: Details of the interaction of vertical a gust and an aircraft with both the DVA and the RAA for λ/Cref = 1 and λ/Cref = 2. Note that the agreement between both methods is surprisingly well for the interaction with the wing of the plane, and also for the interaction with the tail. This is something unexpected, because before reaching the tail, the gust experiences two interactions with
Chapter 6. Numerical results 67 the aircraft’s aerodynamics (with the fuselage and with the wing), so in principle, there should be larger discrepancies between the DVA and RAA. As said before, lift cancellations are supposed to be the reason for the good agreement of both methods. Convergence of the solutions has been checked in terms of the residual and CL. The convergence history of the steady simulation for λ/Cref = 1 is shown in figure 6.10, while figure 6.11 presents the convergence of the unsteady computation for λ/Cref = 1 and the RAA. A zoom-in of it is shown in figure 6.12. Similar plots have been obtained for the rest of the computations. In order to improve the convergence properties of the steady simulation, 500 MG-cycles using an upwind scheme with AUSMDV have been first performed, followed by 9000 MG-cycles using a central scheme with scalar dissipation. In both cases, local timestepping with the Runge-Kutta scheme has been used. It can be seen that the residual decreases some orders of magnitude, and that after 6000 MG-cycles, it is of order 10−6. The lift coefficient has also a constant value, so the convergence criterion is fulfilled. The unsteady simulation has been performed using a central scheme with scalar dissipation and the DTS method with 2750 physical time steps, with ∆tdefined in (6.3). Convergence within each time step is achieved if the residual has decreased some orders of magnitude and the CLhas an horizontal tangent, which can be seen in figure 6.12. As in the 2D unsteady simulations, the number of MG-cycles in dual time that are performed within each time step changes depending on the necessities of the simulation. 0 2000 4000 6000 8000 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 MG-cycles Residual [-] CL[-] Figure 6.10: Convergence history of the residual and the lift coefficient for the steady simulation with Ma∞= 0.78 and λ/Cref = 1.
Chapter 6. Numerical results 68 0 20000 40000 60000 80000 10-7 10-6 10-5 10-4 10-3 10-2 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 MG-cycles Residual [-] CL[-] Figure 6.11: Convergence history of the residual and the lift coefficient for the unsteady simulation with Ma∞= 0.78, λ/Cref = 1 and the RAA. 35000 35100 35200 35300 35400 35500 10-6 10-5 10-4 10-3 10-2 0.385 0.39 0.395 0.4 0.405 MG-cycles Residual [-] CL[-] Figure 6.12: Zoom-in of figure 6.11.
Chapter 6. Numerical results 69 6.1.2.2 Interaction with a lateral gust For the interaction with a lateral gust, the DVA and the RAA are compared for λ/Cref = 1, λ/Cref = 2 and λ/Cref = 4 in terms of the lateral force coefficient CFy. This coefficient is defined as CFy=2Fy ρ∞u2 ∞A(6.4) where Fyis the lateral force (in span direction) that acts on the aircraft. Figures 6.13, 6.14 and 6.15 show the time history of CFyas a function of the dimensionless time t/tref for an encounter with a gust of λ/Cref = 1, λ/Cref = 2 and λ/Cref = 4, respectively. For all the cases, the on-flow Mach number is Ma∞= 0.78 and the on-flow Reynolds number is Re∞= 24.98 ·106. The value of CFyis equal to zero at the beginning of the simulation, because the lateral gust is still far away from the aircraft and it does not influence the flow around it. The first decrease of CFycorresponds to the interaction with the fuselage of the aircraft. After that, the gust reaches the wings. It can be seen that the change of CFyduring the interaction between the wings and the gust really depends on the value of λ: for example, for λ/Cref = 4, there is a small change on the value of CFy. Finally, the gust interacts with the Vertical Tail Plane (VTP), yielding to the major decrease of CFy. 9 11 13 15 17 19 21 23 25 27 −0.04 −0.03 −0.02 −0.01 0 0.01 t/tref [-] CFy[-] RAA, λ/Cref = 1 DVA, λ/Cref = 1 Figure 6.13: Lateral force coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with λ/Cref = 1.
Chapter 6. Numerical results 70 9 11 13 15 17 19 21 23 25 27 −0.04 −0.03 −0.02 −0.01 0 0.01 t/tref [-] CFy[-] RAA, λ/Cref = 2 DVA, λ/Cref = 2 Figure 6.14: Lateral force coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with λ/Cref = 2. 9 11 13 15 17 19 21 23 25 27 −0.04 −0.03 −0.02 −0.01 0 0.01 t/tref [-] CFy[-] RAA, λ/Cref = 4 DVA, λ/Cref = 4 Figure 6.15: Lateral force coefficient history as a function of the dimensionless time predicted by the DVA and the RAA, with λ/Cref = 4.
Chapter 6. Numerical results 71 The prediction error of the DVA is defined as eCFy,max =CFy,max,RAA −CFy,max,DV A CFy,max,RAA(6.5) being CFy,max,DV A and CFy,max,RAA the maximum lateral force coefficient computed with the DVA and the RAA, respectively, which correspond to the interaction with the VTP. Table 6.3 presents the value of this error for the three wavelenghts considered: λ/Cref [-] eCFy,max [%] 1 2.47 2 6.40 4 9.73 Table 6.3: eCFy,max values for the 3D computations involving a lateral gust. As it can be seen, the larger λis, the greater is the value of eCFY,max. Note, however, that this is totally opposite to what is predicted for the interaction with the fuselage: figures 6.13, 6.14 and 6.15 show that the greatest differences between the DVA and the RAA correspond to λ/Cref = 1. This change of behaviour is assumed to be due to the two vortices that emerge from the tips of the wings. Hence, when the gust interacts with the wings, it is perturbed due to the presence of the vortices, and this lasts until the end of the computation. The presence of the vortices is shown in the following figure: Figure 6.16: Slice at xg=−29 m of the y−component of the velocity field at time step 1000, for λ/Cref = 4 and the RAA.
Chapter 6. Numerical results 72 The effect of the vortices on the gust can be clearly seen in the following figures, which show the y−component of the velocity field at time step 1900, for all the values of λ and the RAA: Figure 6.17: Slice at zg= 0 m of the y−component of the velocity field at time step 1900, for λ/Cref = 1 and the RAA. Figure 6.18: Slice at zg= 0 m of the y−component of the velocity field at time step 1900, for λ/Cref = 2 and the RAA.
Chapter 6. Numerical results 79 Both, a steady and an unsteady simulation, have been done for each of the test cases just described. These computations have the following features in common: •The initial distance between the leading edge of the airfoil and the right vortex is 33 m, i.e., there is no influence on the flow around the airfoil. •The cores of the two vortices are separated by b= 21.5 m, corresponding to a VFW-Fokker 614 jetliner. •The circulation is Γ = 300 m2/s. •The core radius is rc= 0.75 m. •The reference Mach and Reynolds numbers are Ma∞= 0.5 and Re∞= 11.06·106, respectively. Figures 6.29, 6.30 and 6.31 show the history of the CLas a function of t/tref for the different pairs of equivalent test cases. As expected, each equivalent pair leads to the same results. For the sake of simplicity, the CLtime history corresponding to the first test case (figure 6.29) is discussed in what follows. Nevertheless, equivalent analysis can be done for the rest of the test cases. 10 20 30 40 50 60 70 80 −1 −0.75 −0.5 −0.25 0 0.25 0.5 0.75 1 t/tref [-] CL[-] test case 1 test case 2 Figure 6.29: Time history of the lift coefficient for test cases 1 and 2.
Chapter 6. Numerical results 80 10 20 30 40 50 60 70 80 −1 −0.75 −0.5 −0.25 0 0.25 0.5 0.75 1 t/tref [-] CL[-] test case 1 test case 3 Figure 6.30: Time history of the lift coefficient for test cases 1 and 3. 10 20 30 40 50 60 70 80 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 t/tref [-] CL[-] test case 4 test case 5 Figure 6.31: Time history of the lift coefficient for test cases 4 and 5. Since the airfoil is symmetric and has a zero AoA, the variations of the lift coefficient are purely created due to the interaction with the vortices. Thus, CL= 0 for t/tref ≤15 and t/tref ≥75: the flow around the airfoil is not perturbed by the wake vortices, which are far away from it. For t/tref >15, the vortices advance in the domain. When the right one starts to disturb the flow, there is a lift growth, because the vortex rotates counter-clockwise. The closer is the airfoil to rc, the stronger is the velocity field of the vortex, as it is shown in figure 6.32(a) for the local z−component of the velocity (plotted in blue). As the right vortex advances in time, stronger loads act on the airfoil, and
Chapter 6. Numerical results 81 there is a progressive growth of the lift. The maximum CLcorresponds to the maximum value of the tangential velocity at rc(see figure 4.5). xg zg (a) Counter-clockwise flow, before the interaction with the core of the right vortex. xg zg (b) Clockwise flow, after the interaction with the core of the right vortex. Figure 6.32: Local z−component of the velocity field (in blue) of the right vortex. The situation corresponds to the first test case, and it is not drawn at scale. Just after interacting with the core of the right vortex, the direction of the flow changes, as it is shown for the tangential velocity in figure 4.5. The situation resembles the one from figure 6.32(b), which results in a strong lift decrease. Once the right vortex is away from the airfoil, its influence on the flow decreases. Hence, the loads applied on the airfoil also decrease, and the value of CLgrows for t/tref ∈(35,45]. From t/tref = 46, the flow starts to be influenced by the left vortex, which rotates clockwise (see figure 4.5). Thus, the lift decreases, and after the airfoil interacts with the core of the vortex, the direction of the flow changes and the lift grows. From t/tref = 56, the influence of the left vortex on the airfoil decreases, and the initial lift is recovered, when both wake vortices are far away from the airfoil. Notice the difference between figures 6.29 and 6.31: the values of CLpredicted with the Lamb-Oseen model are greater than the ones predicted with the Burnham-Hallock model. This agrees with figure 4.6, where for the same value of rc, the maximum tangential velocity prescribed by the Lamb-Oseen model is greater than the one prescribed by the Burnham-Hallock model. Hence, the greater the maximum tangential velocity, the greater is the maximum lift. The convergence of the solution is analysed for the first test case in figures 6.33, 6.34 and 6.35. Similar results have been obtained for the rest of test cases. In the case of the steady simulation (figure 6.33), the residual is of order 10−14 after 225 MG-cycles, while the value of CLis constant. A central scheme with matrix dissipation is used to achieve convergence without much computational effort both for the steady and the unsteady simulations.
Chapter 6. Numerical results 82 In the unsteady simulation, 20 MG-cycles in pseudo-time have been performed within each physical time step ∆t, where ∆t= 0.0001 s. Figure 6.34 shows that the residual decreases some orders of magnitude for each physical time step. This is sufficient to reach the pseudo-steady state, as indicated by the behaviour of CLin figure 6.35. 0 25 50 75 100 125 150 175 200 225 10-14 10-12 10-10 10-8 10-6 10-4 10-2 100 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 MG-cycles Residual [-] CL[-] Figure 6.33: Convergence history of the residual and the lift coefficient for the steady simulation of test case 1. 0 20000 40000 60000 80000 100000 10-9 10-7 10-5 10-3 10-1 101 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 MG-cycles Residual [-] CL[-] Figure 6.34: Convergence history of the residual and the lift coefficient for the unsteady simulation of test case 1.
Chapter 6. Numerical results 83 32100 32150 32200 32250 32300 32350 10-6 10-5 10-4 10-3 10-2 10-1 100 101 102 0.17 0.172 0.174 0.176 0.178 0.18 MG-cycles Residual [-] CL[-] Figure 6.35: Zoom-in of figure 6.34.
Chapter 7 Conclusions 7.1 Summary of the results This section presents an overview of the main results achieved in the numerical computations of this thesis, where the accuracy of the DVA has been compared to the RAA for wind gusts, and the implementation of the DVA to simulate wake vortex encounters has been verified. As explained in chapter 6, this has been done through encounters between airfoils or aircraft and atmospheric disturbances. In particular, the following simulations have been performed: •Vertical gust encounters with a 2D NACA0012-HTP configuration. •Vertical gust encounters with the 3D LamAiR aircraft. •Lateral gust encounters with the 3D LamAiR aircraft. •Wake vortex encounters with a 2D NACA0012 airfoil. As stated in chapter 1, one of the objectives of this thesis is to reduce the computational cost of the simulations. This has been mainly done by varying the number of MG-cycles depending on the requirements of the computation and by using the Chimera technique to set up the grids and enable the movement of a high-resolution grid for a gust. Vertical gust encounters with a 2D NACA0012-HTP configuration The work done in [12] has been taken as a starting point for the 2D comparison of the DVA and RAA. However, the grid density study and the grid set-up have been improved, which has resulted in a reduction of the computational cost. In order to check the influence of compressibility, two different on-flow Mach numbers have been 84
Chapter 7. Conclusions 85 selected: for Ma∞= 0.28, there is a nearly incompressible flow, while compressibility effects are expected for Ma∞= 0.78. The two gust approaches have been compared in terms of the maximum lift coefficient for both Mach numbers and three different gust wavelengths, through the prediction error of the DVA. On the one hand, the effect of compressibility results in greater discrepancies between the DVA and the RAA for Ma∞= 0.78, independently of the wavelength of the gust. On the other hand, the greatest differences between both methods correspond to the smallest wavelength, independent of the Mach number. Therefore, a safe range of validity of the DVA corresponds to gust wavelengths greater than or equal to two reference chord lengths and an on-flow Mach number less than or equal to 0.78. These results agree with the ones from [12]. Vertical gust encounters with the 3D LamAiR aircraft In order to study the influence of the fuselage in the comparison of the DVA and the RAA, 3D computations have been made with Ma∞= 0.78 and a vertical gust. The number of nodes of the simulation is reduced by only considering one half of the aircraft, so the computation is less costly. The two gust modelling methods are compared in terms of the maximum lift coefficient, which is related to the interaction with the wing. The results obtained do not agree with the ones from the 2D simulations, which is totally unexpected. In particular, the DVA results are acceptable for all the wavelengths tested, since the DVA prediction errors are surprisingly small. In addition, the agreement between the DVA and the RAA is better for a gust wavelength equal to two reference chord lengths than to one reference chord length. Presumably, there are lift cancellations due to the global character of this magnitude that cause these unexpected results. Lateral gust encounters with the 3D LamAiR aircraft Simulations involving the complete aircraft and a lateral gust have been performed with Ma∞= 0.78. In this case, the DVA and the RAA have been compared in terms of the maximum lateral force coefficient, which is related to the interaction with the VTP. The wake vortices that emerge from the tips of the wings influence the flow, so the gust is perturbed by the presence of the vortices (from the interaction with the wings until the end of the computation). On the one hand, when the gust has not reached the wings, the discrepancies between the DVA and the RAA are greater for small wavelengths. On the other hand, during the interaction with the VTP, the most important discrepancies between both methods are found for the largest wavelength (four reference chord lengths). This change of behaviour is presumably associated to the presence of the
Chapter 7. Conclusions 86 vortices, whose effect on the gust is more important for large wavelengths. Nevertheless, the computations can also be influenced by the resolution of the near field grid, which is maybe insufficient. Therefore, no range of validity for the DVA and a lateral gust can be given at the moment, because the resolution of the grid should be improved. Wake vortex encounters with a 2D NACA0012 airfoil The DVA has been implemented in TAU to enable encounters between an aircraft and the wake vortices from another aircraft. The implementation has been done in the most general way, such that encounters between aircraft flying in different directions can be simulated. To verify the correctness of the implementation, some equivalent test cases have been defined. To reduce the computational cost, the simulations consider an inviscid flow and a NACA0012 airfoil. The time history of the lift coefficient that has been obtained for each pair of equivalent test cases is the same, so the implementation can be considered successful. 7.2 Future work The work that has been done in this thesis can be extended and improved in the future. For doing so, it is worth to take the following suggestions into account: •Study the global distribution of the lift force in the span direction. This is necessary to ensure that the small values of eCL,max obtained for the interaction between the LamAiR aircraft and a vertical gust are really caused by lift cancellations. In addition, it might be recommended to compute the prediction error of the DVA in terms of another magnitude, like for instance the wing root bending moment. •Improve the resolution of the near field grid for the LamAiR aircraft. As said in chapter 5, the lateral gust simulations consider only two hole geometries, so during the interaction with the aircraft, the gust is not contained in the nodes of the moving grid but in the nodes of the near field grid. Thus, it could happen that not only the wake vortices but also the lower resolution of the near field grid have an influence on the numerical results. In order to discard this, the resolution of the near field grid should be improved. •Do more realistic simulations by considering not only the influence of the gust and the aircraft’s aerodynamics, but also the coupling to other relevant disciplines. In particular, the next logical step is to consider the coupling to flight mechanics and also the fluid-structure interaction.
Chapter 7. Conclusions 87 •Implement the so-called Source Gust Model, which is a third method for gust modelling developed in [17]. Like the RAA, this method allows a mutual interaction of gust and aircraft, but at a lower computational cost. Therefore, the main disadvantage of the RAA (very costly simulations) could be avoided. •Implement the RAA to allow encounters with wake vortex turbulence, by prescribing unsteady boundary conditions for the velocity components, the pressure and the density at the farfield. As in the case of gust encounters, the RAA has to be compared to the DVA, to set the range of validity of the last approach.
Appendix A Conservation laws for moving grids If a control volume Ω moves with a certain velocity vb, the conservation laws of mass, momentum and energy can be written in integral form as ∂ ∂t ZΩ ρdΩ + I∂Ω ρ[(v−vb)·n]dS = 0 ∂ ∂t ZΩ ρvdΩ + I∂Ω ρv[(v−vb)·n]dS =−I∂Ω pndS +I∂Ω (τ·n)dS ∂ ∂t ZΩ ρEdΩ + I∂Ω ρE [(v−vb)·n]dS =I∂Ω k(∇T·n)dS +ZΩ ˙qhdΩ−I∂Ω p(v·n)dS +I∂Ω (τ·v)·ndS (A.1) Equations (A.1) can be rewritten in a more compact way ∂ ∂t ZΩ WdΩ + I∂Ω (Fc−Fv)dS =0(A.2) where Wis the vector of the conservative variables, which has the following components W= ρ ρv ρE = ρ ρu ρv ρw ρE (A.3) 88
Bibliography 95 [23] T. Schwarz. Ein blokstrukturiertes Verfahren zur Simulation der Umstr¨omung komplexer Konfigurationen. Technical Report DLR-FB 2005-20, Deutsches Zentrum f¨ur Luft- und Raumfahrt, Braunschweig, Germany, 2005. [24] A. Seitz, M. Kruse, T. Wunderlich, J. Bold, and L. Heinrich. The DLR Project LamAiR: Design of a NLF Forward Swept Wing for Short and Medium Range Transport Application. 29th AIAA Aplied Aerodynamics Conference, Honolulu, Hawaii, United States of America, 2011. [25] P. R. Spalart and S. R. Allmaras. A One-Equation Turbulence Model for Aerodynamic Flows. Recherche Aerospatiale, No. 1, pp. 5-21, 1994. [26] P. D. Thomas and C. K. Lombard. Geometric conservation law and its application to flow computations on moving grids. AIAA Journal, Vol. 17, No. 10, pp. 1030–1037, DOI 10.2514/3.61273, 1979. [27] V. Vandeucarter. CFD simulation of atmospheric wind gusts. A CFD simulation of three different gust models and its first effects on a NACA4415 airfoil. Master’s thesis, Delft University of Technology, The Netherlands, 2011. [28] C. Wolf. A Chimera Simulation Method and Detached Eddy Simulation for Vortex- Airfoil Interactions. PhD thesis, Georg-August University, G¨ottingen, Germany, 2010. [29] S. Zan and et al. Wing and Fin Buffet on the Standard Dynamic Model. RTO Technical Report, RTO-TR-26 pp. 361-381, 2000.