Study of unstructured finite volume methods for the solution of the Euler equations

Study of unstructured finite volume methods for the solution of the Euler equations Master Thesis Author: Iván Padilla Montero Supervisor: Roberto Maurice Flores Le Roux Document: Report & Appendix Call: July 2016 Master’s Degree in Aerospace Engineering Universitat Politècnica de Catalunya Terrassa - June 17, 2016 Abstract This work deals with the numerical solution of inviscid compressible flows by means of the Euler equations. It focuses on the description of an unstructured finite volume method for these equations and its numerical application to solve external, two-dimensional steady problems. On first place, the standard formulation of the Euler equations is presented, reviewing the most important properties that characterize their mathematical behavior. The hyperbolic nature of the system is discussed, emphasizing the fundamental importance of taking into account the propagation of information in the flow field in order to obtain physically meaningful solutions, which also leads to a description of how the boundary conditions should be treated to avoid undesirable behaviors. To complete this presentation, a dimensionless form of the equations is derived, which provides substantial advantages to the numerical solution. The attention is then focused on the unstructured finite volume formulation, which is based on a central approximation of the fluxes at the volume interfaces. According to the need of properly accounting for the propagation of characteristic variables, the requirement to add artificial dissipation terms to the central discretization is justified. Then, two classical forms of artificial dissipation are defined, namely, the first-order upwind scheme and the JamesonSchmidt-Turkel high-order model, detailing how to adapt the formulation of the dissipation terms to an unstructured mesh. Eventually, the time integration of the spatially discretized equations is assessed. With the objective of performing a practical implementation of the theoretical concepts studied, the development of a numerical solver is presented next, briefly describing the program structure and characteristics. After that, five different test cases are solved with the purpose of validating the code, consisting on two transonic flows around a NACA0012 airfoil and three supersonic examples, respectively around a NACA0012 airfoil, a double wedge airfoil and circular cylinder. The results obtained for each case are then analyzed and compared against reference solutions, showing an overall satisfactory performance of the solver developed. Resumen Este trabajo trata sobre la solución numérica de flujos compresibles no viscosos mediante las ecuaciones de Euler. La atención se centra en la descripción de un método de volúmenes finitos no estructurado para estas ecuaciones y en su aplicación numérica para la solución de flujos bidimensionales, externos y estacionarios. En primer lugar se presenta la formulación estándar de las ecuaciones de Euler, revisando las propiedades más importantes que caracterizan su comportamiento matemático. Se describe la naturaleza hiperbólica del sistema, haciendo hincapié en la importancia fundamental de tener en cuenta la propagación de información en el fluido para poder obtener soluciones físicas correctas, lo que lleva directamente a especificar como deben tratarse las condiciones de contorno a fin de evitar posibles comportamientos indeseados. Esta presentación se completa exponiendo la utilización de una forma adimensional de las ecuaciones, lo que aporta ventajas sustanciales a la solución numérica. A continuación, el documento se centra en la explicación del método de volúmenes finitos no estructurado escogido, que se basa en una aproximación central de los flujos en la entrefase de los volúmenes. De acuerdo con la necesidad de tener en cuenta la propagación de las variables características en el fluido, se procede a justificar el requerimiento de introducir difusión artificial en la discretización producida por el esquema central. Se prosigue entonces con la descripción de dos métodos clásicos de difusión artificial, que son el esquema aguas arriba de primer orden y el modelo de alto orden de Jameson-Schmidt-Turkel, detallando como adaptar la formulación de los términos disipativos a la formulación no estructurada. Finalmente se discute la integración temporal de la discretización espacial de las ecuaciones. Con el objetivo de realizar una aplicación práctica directa de los conceptos teóricos estudiados, a continuación se presenta el desarrollo de un código de simulación numérica, describiendo en términos generales la estructura del programa y sus características. Con el propósito de validar el código, se pasa entonces a la solución de cinco ejemplos de aplicación, que consisten en dos flujos transónicos alrededor de un perfil NACA0012 y tres casos supersónicos, respectivamente alrededor de un perfil NACA0012, un perfil romboidal y un cilindro circular. Los resultados obtenidos en cada caso se comparan con una solución de referencia, demostrando un comportamiento general satisfactorio del código desarrollado. Contents 1 INTRODUCTION ........................................................................................................... 1 1.1 Motivation ............................................................................................................... 1 1.2 Objectives and scope ............................................................................................. 2 2 THE EULER EQUATIONS ............................................................................................ 3 2.1 Differential conservation form of the Euler equations ............................................. 4 2.2 Some mathematical properties of the Euler equations ........................................... 5 2.2.1 Quasi-linear form ............................................................................................. 5 2.2.2 Eigenvalues of the system .............................................................................. 6 2.2.3 One-dimensional linearized characteristic formulation .................................... 6 2.3 Physical boundary conditions for external inviscid flow .......................................... 8 2.4 Nondimensionalization ........................................................................................... 9 3 UNSTRUCTURED FINITE VOLUME METHOD FOR THE EULER EQUATIONS...... 11 3.1 Spatial discretization ............................................................................................ 12 3.1.1 Central schemes and the need for artificial dissipation ................................. 14 3.2 Artificial dissipation ............................................................................................... 16 3.2.1 First-order artificial dissipation ...................................................................... 19 3.2.2 Jameson-Schmidt-Turkel high-order artificial dissipation .............................. 22 3.3 Numerical treatment of boundary conditions ........................................................ 25 3.3.1 Body surface boundary condition .................................................................. 25 3.3.2 Far-field boundary condition .......................................................................... 27 3.4 Time integration .................................................................................................... 28 3.4.1 Calculation of the time step ........................................................................... 29 3.4.2 A relaxed update procedure to promote the positivity of thermodynamic variables ...................................................................................................................... 30 4 NUMERICAL IMPLEMENTATION AND TEST CASES .............................................. 31 4.1 Description of the solver developed ..................................................................... 31 4.2 Test cases ............................................................................................................ 33 4.2.1 Calculation of nondimensional quantities for the analysis of results ............. 34 4.2.2 Transonic flow around a NACA0012 airfoil: case 1 ....................................... 36 4.2.3 Transonic flow around a NACA0012 airfoil: case 2 ....................................... 38 4.2.4 Supersonic flow around a NACA0012 airfoil ................................................. 38 4.2.5 Supersonic flow past a double wedge airfoil ................................................. 44 4.2.6 Supersonic flow past a circular cylinder ........................................................ 47 5 CONCLUSIONS AND FUTURE WORK ...................................................................... 50 Bibliography ........................................................................................................................ 51 APPENDIX: Practical calculation of geometrical quantities for triangular cells .................. 53 List of figures Figure 3.1 Schematic of an unstructured triangular mesh. ................................................ 13 Figure 3.2 Example of strong unphysical oscillations near a shock wave on a 10 degree compression corner. This solution was obtained using a central finite difference scheme with no artificial dissipation. ................................................................................................ 14 Figure 3.3 Example geometry for reconstruction of the fourth difference stencil. The dummy nodes are located in the line joining the centroids of cells i and j. ......................... 23 Figure 4.1 Diagram showing the different structural blocks and flow of the developed code. ................................................................................................................................... 32 Figure 4.2 Different views of the mesh used for the NACA0012 test cases. ..................... 36 Figure 4.3 Pressure coefficient and Mach number results for case 1. ............................... 40 Figure 4.4 Pressure coefficient and Mach number results for case 2. ............................... 41 Figure 4.5 Pressure coefficient and Mach number results for case 3. ............................... 42 Figure 4.6 Convergence history and entropy generation for cases 1 (top), 2 (middle) and 3 (bottom). ............................................................................................................................. 43 Figure 4.7 Pressure coefficient and Mach number results for case 4. ............................... 45 Figure 4.8 Convergence history and entropy change for case 4. ...................................... 46 Figure 4.9 Details of the mesh used for the numerical solution of case 4. ........................ 46 Figure 4.10 Pressure coefficient and Mach number results for case 5. ............................. 48 Figure 4.11 Convergence history and entropy change for case 5. .................................... 49 Figure 4.12 Detail views of the mesh for case 5. ............................................................... 49 List of tables Table 4.1 Summary of the different test cases considered. ............................................... 33 Table 4.2 Comparison of the aerodynamic force and moment coefficients for case 1. ...... 37 Table 4.3 Summary of the aerodynamic force and moment coefficients for case 2. .......... 38 Table 4.4 Comparison of lift, drag and moment coefficients for case 3. ............................ 39 Table 4.5 Comparison of aerodynamic force coefficients for case 4. ................................. 44 Iván Padilla Montero 1 1 INTRODUCTION 1.1 Motivation During the last decades, the rapid increase in the computational power available has promoted the development of numerical analysis techniques in many different engineering fields. Due to the technological challenges that characterize the design and construction of aircraft, and the important costs associated with wind tunnel testing, aerospace engineering is one of the areas where numerical methods show a wide range of applications. At present, computational solid mechanics and computational fluid dynamics constitute two fundamental pillars of any serious design process in the aerospace industry. It is necessary, then, for the modern aerospace engineer to develop a solid understanding of the foundations upon which these disciplines build-up. With the numerical solution of complete flow fields around complex geometries becoming a routine calculation nowadays, there is a strong necessity to develop robust and efficient programs that allow the obtention of accurate results in a short time. However, even with the most recent developments in computing performance, the solution of the full system of Navier-Stokes equations is still an expensive task from the computational point of view, especially if the turbulence effects are to be modelled with high fidelity. This situation has motivated the search for alternate strategies, which focus on simplified flow descriptions that, if applied with care, can capture the relevant physics of the complete model at a lower computational cost. The best simplified model for the solution of actual aircraft configurations is that given by the Euler equations, which allow the obtention of realistic results in the transonic and low to moderate supersonic regimes when flow separation is not important (Anderson, 2011). These equations, which govern the behavior of inviscid flows, offer a very attractive alternative especially in the case of preliminary design stages. They are able to capture the compressible flow phenomena associated to the presence of discontinuities in the flow field, such as shock waves, which shows to be one of the most important design drivers for any machine flying at or near the supersonic regime. Although the inviscid flow model is obviously not of universal validity, the importance of its accurate numerical solution also resides in the dominating convective character of the Navier-Stokes equations at high Reynolds numbers, which is the case of the immense majority of flow situations encountered in practice. Since the Euler equations retain the convective properties of the general formulation, almost all of the methods developed for the Euler system are also valid for the Navier-Stokes equations. As a result, another advantage of developing solvers for the Euler equations is that they serve as the base for possible extensions to the complete, viscous model. Iván Padilla Montero 2 1.2 Objectives and scope With the previous ideas in mind, the objective of this work is the numerical solution of compressible inviscid flows by means of the Euler equations. This has been the subject of extensive research for several years, and many different techniques have been developed and can be found in the literature; see for example (Hirsch, 1990). Here, the aim is to learn the fundamental concepts that lie behind the practical implementation of inviscid flow solvers for the aerospace industry. Accordingly, the scope is focused on the solution of twodimensional, external, steady inviscid flows, which retain all the important properties of the numerical calculations carried out in practice. To achieve the proposed goal, a numerical solver has been developed. Due to its natural connection with conservation laws, a finite volume method has been chosen for the discretization of the governing equations, based on an unstructured mesh strategy. Among the different numerical schemes available, a well-established central scheme has been adopted, which relies on a combination of simple artificial dissipation terms and a multistage time integration method to offer a robust solution to most of the flow situations under consideration (Jameson, Schmidt, & Turkel, 1981; Hirsch, 1990). The use of an unstructured approach has nowadays the enormous advantage of allowing an almost automatic mesh generation on arbitrary geometries. This is an essential aspect of any practical application of computational fluid dynamics to real engineering problems, and although the scope of the numerical solutions treated here is not that big, it is important to follow this philosophy when looking into the future. Furthermore, the possibility to perform local refinements in a certain region without affecting the rest of the domain opens the way to perform flexible mesh adaptation techniques that can optimize the number of grid points for a given level of accuracy. Five different examples of application have been carried out in order to investigate the performance of the code developed, which cover the transonic and supersonic regimes as well as different body geometries. In all the cases considered, the results obtained are compared against reference solutions. Iván Padilla Montero 9 information in the flow field is different depending on whether the flow is locally subsonic or supersonic. Then, in order to properly account for the correct transport of characteristic variables across the far-field boundary, there is a need to differentiate between two cases. On one side, if the normal freestream velocity component is supersonic, all the characteristics move in the same direction (all positive eigenvalues), meaning that all the flow variables have to be prescribed at the inlet portion the far-field boundary and none of them have to be imposed at the outlet portion. On the other side, for a subsonic normal velocity there is a negative eigenvalue, telling us that one component of the solution moves upstream. In this case, for a two-dimensional flow, only three conditions shall be prescribed at the inlet portion and one at the outlet. Finally, it is also important to take into account that the far-field boundary condition should be imposed on the characteristic variables only, not on the primitive or conservative ones. If the characteristic variables are not properly prescribed at the boundary, their associated waves may not be able to leave or enter the computational domain in the correct way, causing unwanted reflections which decrease the numerical quality of the solution. 2.4 Nondimensionalization The use of the equations in nondimensional form is very convenient from the numerical perspective. First, the use of normalized quantities contributes to the numerical quality of the solution by enhancing the conditioning of the system of the equations. Second, the use of dimensionless variables allows us to work with the similarity parameters of the flow, which reduces the number of variables involved in the calculations and yields more general results. A useful set of dimensionless variables for the Euler equations is 𝑡=𝑡𝑐∞ 𝐿 , 𝑥=𝑥𝐿 , 𝑦=𝑦𝐿 , 𝑢=𝑢 𝑐∞ , 𝑣=𝑣 𝑐∞ 𝑝= 𝑝 𝜌∞𝑐∞ 2 , 𝑒0=𝑒0 𝑐∞ 2 , 𝑇=𝑇 𝑇∞ (2.23) where 𝐿 is a characteristic length of the problem under study (for example the chord of an airfoil) and 𝑐∞, 𝜌∞ and 𝑇∞ respectively denote the freestream speed of sound, density and temperature of the fluid. For a calorically perfect gas, the speed of sound satisfies the following relationship: 𝑐=√𝛾𝑅𝑇. With these variables, the system (2.1) simply transforms in 𝜕𝐔 𝜕𝑡+𝜕𝐅 𝜕𝑥+𝜕𝐆 𝜕𝑦=𝟎 (2.24) Iván Padilla Montero 10 and the respective vectors become 𝐔=[𝜌 𝜌𝑢 𝜌𝑣 𝜌𝑒0], 𝐅= [ 𝜌𝑢 𝜌𝑢2+𝑝 𝜌𝑢𝑣 𝜌𝑢ℎ0 ] , 𝐆= [ 𝜌𝑣 𝜌𝑢𝑣 𝜌𝑣2+𝑝 𝜌𝑣ℎ0 ] (2.25) with 𝑒0=𝑒+𝑉2 2 , ℎ0=𝑒0+𝑝𝜌 (2.26) Then, the thermodynamic relationships for a calorically perfect gas change as 𝑝=𝜌𝑇 𝛾 , 𝑒 = 1 𝛾(𝛾−1)𝑇 (2.27) Similarly, the freestream values of the primitive variables result in the following expressions 𝜌∞=1 , 𝑢∞=𝑀∞ 𝑥 , 𝑣∞=𝑀∞ 𝑦 , 𝑝∞=1𝛾 , 𝑇∞=1 (2.28) which conveniently transform the freestream vector of conservative variables in 𝐔∞= [ 1 𝑀∞ 𝑥 𝑀∞ 𝑦 1 𝛾(𝛾−1)+𝑀∞ 2 2 ] (2.29) where 𝑀∞ 𝑥 and 𝑀∞ 𝑦 denote the components of the freestream Mach number along each spatial direction, given by the angle of attack 𝛼 as 𝑀∞ 𝑥=𝑀∞cos𝛼 , 𝑀∞ 𝑦=𝑀∞sin𝛼 (2.30) Observing equations (2.28) and (2.29), it is easy to notice that the nondimensional solution of the Euler equations only depends on the freestream Mach number (including the angle of attack), the ratio of specific heats and the shape (but not the size) of the body under consideration. These are precisely the similarity parameters for an inviscid flow (Anderson, 2011), which clearly shows the benefits of working with the equations in nondimensional form. From now on, for the sake of clarity, the tilde will be omitted in the equations, but it will be implicitly assumed that the dimensionless form of the equations is being used unless otherwise specified. Iván Padilla Montero 11 3 UNSTRUCTURED FINITE VOLUME METHOD FOR THE EULER EQUATIONS The finite volume method is advantageous for the discretization of conservation laws such as the Euler equations due to its direct connection to the physical flow properties. It is based upon the discretization of the integral form of the conservation equations, as opposed to the finite difference method, which deals with the differential form, and is therefore more general and fundamentally appropriate for the solution of external compressible flows. Assuming a control volume Ω fixed in space, enclosed by a surface Γ, the integral conservation form of the two-dimensional Euler equations can be written as ∫𝜕𝐔 𝜕𝑡 Ω𝑑Ω+∮𝐅𝑛 Γ𝑑Γ=𝟎 (3.1) where 𝐅𝑛 is the flux vector across the boundary of the control volume, given by the unit normal vector 𝐧 pointing outwards from 𝑑Γ 𝐅𝑛=𝑛𝑥𝐅+𝑛𝑦𝐆 (3.2) with 𝑛𝑥 and 𝑛𝑦 being the components of 𝐧. It is important to note that with this formulation, the order of the equations has been decreased by one, which reduces the continuity requirements of the flow variables. This is, as commented before, the reason why this form of the Euler equations allows for mathematical discontinuities in the flow field. The system of integral conservation laws given by equation (3.1) can be either discretized both in space and time simultaneously, or by a spatial discretization with independent time integration. The later approach converts the system of partial differential equations into a system of ordinary differential equations in time, which are then solved using any suitable time integration method. This offers more flexibility regarding the stability and accuracy properties of the time-marching scheme, and is also the approach taken in this work. It may be somewhat surprising at first that in order to obtain the numerical solution to a steady flow, the time dependent Euler equations are considered. There is, however, a strong mathematical reason to justify this choice. We have stated previously that the complete system of Euler equations is hyperbolic, no matter the type of inviscid flow being computed. Nevertheless, when the time derivatives are removed from the system, the mathematical behavior of the equations changes, and it can actually be proved that the system of steady Euler equations exhibits a mixed elliptic-hyperbolic nature (Anderson, 1995). Moreover, this mixed behavior is associated with the local flow regime, namely, the equations are elliptic in the subsonic regions and hyperbolic in the supersonic ones. This situation poses enormous numerical difficulties, since any steady technique that is suitable for the solution of the Iván Padilla Montero 12 subsonic region is usually not valid for the supersonic counterpart, and vice versa. This is the reason why the unsteady Euler equations are also used for the solution of steady flows. It is the steady-state flow field what we want, and the time-dependent approach is simply a means to that end. 3.1 Spatial discretization The formulation given by (3.1) expresses that the variation of the conservative flow variables inside the control volume only depends on the net balance of the flux vectors across its boundary. This implies that for an arbitrary division of the domain into smaller subdomains, we can write the integral conservation equations for each subdomain and recover the global conservation law by simply adding up the contribution of each one. This is the basis of the spatial discretization by means of the finite volume method, i.e., the subdivision of the domain Ω into a series of finite control volumes, also known as cells, and the application of the conservation statement to each one of them. In order to obtain the semi-discrete form of the integral equations, the volume integral in (3.1) is usually replaced by the average value of the vector of conservative variables over the cell, which for a cell with area 𝐴 is defined as (the finite volume is two-dimensional) 𝐔=1𝐴∫𝐔 𝐴𝑑𝐴 (3.3) Regarding the discretization of the flux term, the surface integral can be replaced by the sum over all the bounding faces of the cell, so that the fluxes are assumed constant along each face. This proves to be a second-order approximation, which is the recommended accuracy for the majority of CFD applications. Thus, such an approximation is acceptable for our purpose. Under the previous considerations, for a given generic cell 𝑖 with area 𝐴𝑖, the spatial discretization of equation (3.1) may be expressed as 𝐴𝑖𝜕𝐔𝑖 𝜕𝑡+∑𝐅𝑛𝑒 𝑁𝑒 𝑒=1 𝑙𝑒=𝟎 (3.4) where 𝐔𝑖 is the average vector of conservative variables over the cell, 𝑁𝑒 is the number of faces of the cell, 𝐅𝑛𝑒 denotes the numerical flux vector across the cell face 𝑒, and 𝑙𝑒 is the length of face 𝑒. As before, the numerical flux across the face is determined by the local unit outward normal 𝐧𝑒 𝐅𝑛𝑒=𝑛𝑥𝑒𝐅𝑒+𝑛𝑦𝑒𝐆𝑒 (3.5) Iván Padilla Montero 13 where 𝐅𝑒 and 𝐆𝑒 denote the numerical flux vectors along each spatial component evaluated at the face. This is a quite general formulation of the finite volume method applied to the two-dimensional Euler equations. In practice, when the numerical results obtained from a discretization like the one already presented are to be analyzed, we need to assign the cell-averaged values to a mesh point, for example the center (centroid) of the cell. Taking this into account, it can be shown, see (Lomax, Pulliam, & Zingg, 2001), that the cell-averaged values and the values at the center of the cell only differ by a term of second-order. This means that the volume integral of the original conservation equation can be approximated as the value of the vector of conservative variables at the center of cell, once again with second-order accuracy. As a result, 𝐔𝑖 may be substituted by 𝐔𝑖 in the discrete equation, denoting the conservative variables at the centroid of the cell. At this point, one has to select the type of cell to subdivide the computational domain and choose how to approximate the fluxes at the cell faces. As can be expected, there exist many different options, and a general discussion can be found in (Hirsch, 2007). As announced previously, in this study, an unstructured discretization based on triangular cells has been considered (see Figure 3.1). Triangular cells are the simplest two-dimensional control volumes, and are widely used since they allow a lot of flexibility to discretize the majority of geometries usually encountered (Löhner, 2008). During practical implementation, the calculation of cell areas, centroids, face lengths and face normals is required. The mathematical formulation of such quantities for triangular cells can be found in the Appendix. The evaluation of fluxes at the cell faces is one of the key aspects of any numerical scheme based on the finite volume method. For the Euler equations, we can distinguish essentially between two families: central and upwind discretization schemes. Both of them have been widely applied to the numerical solution of high-speed inviscid flows with satisfactory results (Hirsch, 1990). From their pure definition, upwind schemes are designed to numerically account for the direction of propagation of information in the flow field, whereas central cell 𝑖 cell 𝑗 face 𝑖𝑗 cell centroid Figure 3.1 Schematic of an unstructured triangular mesh. Iván Padilla Montero 14 schemes are based on centered discretizations which can draw numerical information from outside the correct domain of dependence of a given grid point. In this sense, upwind schemes obey more properly the physics of the flow, and can be viewed as a more natural discretization of the Euler equations. However, as will be discussed in the next section, in view of the limitations of pure central discretizations, different corrections have been developed for central schemes, which modify them to mimic the behavior of upwind methods. Here, a second-order central discretization has been considered. 3.1.1 Central schemes and the need for artificial dissipation Focusing on a cell-centered approach, which is the most common practice in CFD, and assuming a piecewise constant approximation, that is, a constant value of the fluid variables inside a given cell, the approximation of numerical fluxes at the interface by the secondorder central scheme can be formulated for an unstructured mesh as 𝐅𝑖𝑗=12(𝐅𝑖+𝐅𝑗) 𝐆𝑖𝑗=12(𝐆𝑖+𝐆𝑗) (3.6) where 𝐅𝑖𝑗 and 𝐆𝑖𝑗 are the spatial components of the numerical flux at the face shared by the cells 𝑖 and 𝑗 (see Figure 3.1), and 𝐅𝑖, 𝐅𝑗, 𝐆𝑖 and 𝐆𝑗 respectively denote the fluxes associated to each one of the cells, namely 𝐅𝑖=𝐅(𝐔𝑖), 𝐅𝑗=𝐅(𝐔𝑗), and the same for the 𝑦 components. Figure 3.2 Example of strong unphysical oscillations near a shock wave on a 10 degree compression corner. This solution was obtained using a central finite difference scheme with no artificial dissipation. Iván Padilla Montero 15 Then, the central discretization approximates the fluxes at a given interface as the simple arithmetic average between the fluxes of the cells sharing that face. This estimation does not take into account the direction of the transport of the characteristic variables in the solution, something which is the source of significant numerical problems, especially in the presence of discontinuities, as described below. It is known that in processes governed by nonlinear equations such as the Euler and NavierStokes systems, there can be a continuous production of high-frequency components in the solution. These are the responsible, for example, of the production of shock waves. In real flows, the production of these high-frequency modes is actually limited by viscosity. However, when dealing with the Euler equations, there is no such limitation. Then, if a nondissipative scheme is used to discretize the equations in the presence of sharp gradients, the numerical errors associated to the discretization usually introduce severe oscillations in the vicinity of discontinuities (see Figure 3.2) and, in some cases, can produce unphysical solutions due to a violation of the entropy condition. As a consequence, the numerical scheme in use must contain some form of numerical dissipation in order to be able to deal with this phenomenon. Upwind schemes, on the other hand, respect the correct propagation of the flow characteristics and do not present these problems. In fact, they show an intrinsic dissipation. Unfortunately, it can be found that the second-order central approximation to a first derivative is nondissipative, see for instance (Lomax, Pulliam, & Zingg, 2001). This means that if the central scheme presented above is to be used to obtain the numerical solution of a high-speed inviscid flow by means of the Euler equations, some numerical dissipation has to be added to the solution. Otherwise, strong oscillations and incorrect, unphysical solutions may be obtained. This explicitly added dissipation is usually referred to as artificial dissipation, artificial diffusion or artificial viscosity, and can also be viewed as a means of stabilizing the numerical solution. Different forms of artificial dissipation have been developed for central schemes (Hirsch, 1990), which are usually based upon the properties of upwind schemes. When simplified equations that retain most of the properties of the original system are discretized using upwind-type discretizations, a modified partial differential equation can be obtained in which the inherent dissipative terms that the upwinding procedure generates appear explicitly. These terms can then be used to design an artificial dissipation model such as to correct the behavior of the central schemes. It can be stated, then, that the artificial dissipation terms introduce an upwind-like correction to the central schemes such as to remove the non-physical effects arising from the central discretization of wave propagation phenomena. In this sense, the terms that are added are not as “artificial” as it may seem, but are strongly connected to the nature of the flow. This is the reason why the use of a central discretization plus dissipation terms can usually be considered equivalent to the use of an upwind technique. The artificial dissipation models that have been considered in this analysis are described in the following section. Iván Padilla Montero 16 3.2 Artificial dissipation In order to add artificial dissipation to the finite volume discretization, the central approximation of the interface fluxes given by (3.6) is generally modified as 𝐅𝑖𝑗=12(𝐅𝑖+𝐅𝑗)−𝐃𝑥𝑖𝑗 𝐆𝑖𝑗=12(𝐆𝑖+𝐆𝑗)−𝐃𝑦𝑖𝑗 (3.7) with 𝐃𝑥𝑖𝑗 and 𝐃𝑦𝑖𝑗 being the artificial dissipation terms associated to the fluxes 𝐅𝑖𝑗 and 𝐆𝑖𝑗. For the purpose of properly understanding the nature of the artificial dissipation terms commonly used in practice, it is important to review the basic philosophy of upwind schemes. Recalling the analysis performed in section 2.2.3, and taking into account the homogeneous property of the Euler equations, a characteristic flux vector can be defined as (Flores, Ortega, & Oñate, 2011) 𝐟=𝚲𝛙 (3.8) which allows rewriting the one-dimensional characteristic formulation into 𝜕𝛙 𝜕𝑡+𝜕𝐟 𝜕𝑥=𝟎 (3.9) Let us assume for the time being that (3.9) is the system of governing equations that we want to solve using the finite volume method. After applying a spatial discretization equivalent to that given by (3.4), we need to decide how to approximate the characteristic fluxes at the cell interfaces. In this one-dimensional case, the cells become segments and the faces reduce to the midpoint between the nodes defining each segment. We now choose to use an upwind scheme, for simplicity the first-order one. Assuming a segment defined by the nodes 𝑖 and 𝑗, and that the positive direction of the 𝑥 axis (from left to right) is from 𝑖 to 𝑗, the first-order upwind approximation for each characteristic component 𝑘 can be written as 𝑓𝑘𝑖𝑗=𝑓𝑘𝑖 if 𝜆𝑘>0 𝑓𝑘𝑖𝑗=𝑓𝑘𝑗 if 𝜆𝑘<0 (3.10) for 𝑘=1,2,3. As can be seen, the previous expression takes into account the direction of propagation of the wave components in the flow, as given by the sign of the corresponding eigenvalues, which are assumed piecewise constant between 𝑖 and 𝑗 due to the approximation previously considered. If the wave is propagating from left to right (𝜆𝑘>0), Iván Padilla Montero 17 the value of the interface flux is assumed to be that of node 𝑖, that is, upstream of the midpoint 𝑖𝑗. On the contrary, if the wave travels from right to left (𝜆𝑘<0), the assigned flux value is the one of node 𝑗, once again upstream of the midpoint. It is clear, then, which is the fundamental idea of upwinding. Using expression (3.8), equation (3.10) can be transformed in 𝑓𝑘𝑖𝑗=𝜆𝑘𝜓𝑘𝑖 if 𝜆𝑘>0 𝑓𝑘𝑖𝑗=𝜆𝑘𝜓𝑘𝑗 if 𝜆𝑘<0 (3.11) Now, to see the equivalence between the central scheme with artificial dissipation and the upwind scheme, we can also write the upwind approximation in (3.11) as 𝑓𝑘𝑖𝑗=12(𝑓𝑘𝑖+𝑓𝑘𝑗)−12𝜆𝑘(𝜓𝑘𝑗−𝜓𝑘𝑖) if 𝜆𝑘>0 𝑓𝑘𝑖𝑗=12(𝑓𝑘𝑖+𝑓𝑘𝑗)+12𝜆𝑘(𝜓𝑘𝑗−𝜓𝑘𝑖) if 𝜆𝑘<0 (3.12) which is the same as the second-order central approximation plus an additional term which takes into account the direction of propagation of information. Both cases in (3.12) can be combined in a single expression by making use of the absolute value of the eigenvalues (Hirsch, 1990; Lomax, Pulliam, & Zingg, 2001; Lyra & Morgan, 2000) 𝑓𝑘𝑖𝑗=12(𝑓𝑘𝑖+𝑓𝑘𝑗)−12|𝜆𝑘|(𝜓𝑘𝑗−𝜓𝑘𝑖) (3.13) Then, defining |𝚲| as the diagonal matrix which contains the absolute value of the eigenvalues in its diagonal terms, the previous equation can be directly expressed in terms of the vectors of characteristic variables and fluxes, that is 𝐟𝑖𝑗=12(𝐟𝑖+𝐟𝑗)−12|𝚲|(𝛙𝑗−𝛙𝑖) (3.14) This constitutes the first-order upwind approximation to the interface characteristic fluxes for the formulation given by equation (3.9). Focusing again on the conservative form, the relationship (2.19) can be used to change from characteristic to conservative variables, obtaining 𝐅𝑖𝑗=12(𝐅𝑖+𝐅𝑗)−12𝐑|𝚲|𝐑−1(𝐔𝑗−𝐔𝑖) (3.15) where the product 𝐑|𝚲|𝐑−1=|𝐀| (3.16) Iván Padilla Montero 18 is defined as the positive flux Jacobian. It is obtained by taking the absolute value of all the eigenvalues of the flux Jacobian 𝐀. With this definition, the upwind approximation to the interface flux for the one-dimensional Euler equations becomes 𝐅𝑖𝑗=12(𝐅𝑖+𝐅𝑗)−12|𝐀|(𝐔𝑗−𝐔𝑖) (3.17) This is an important result, since it states that all the information necessary to properly account for the propagation of characteristic quantities in the flow is contained inside the positive flux Jacobian. Then, |𝐀| plays a key role in any upwind discretization. The extension to multiple dimensions is straightforward, see for instance (Hirsch, 1990), allowing us to write 𝐆𝑖𝑗=12(𝐆𝑖+𝐆𝑗)−12|𝐁|(𝐔𝑗−𝐔𝑖) (3.18) Now compare equation (3.7) with (3.17) and (3.18). It is clear that choosing 𝐃𝑥𝑖𝑗=12|𝐀|(𝐔𝑗−𝐔𝑖) 𝐃𝑦𝑖𝑗=12|𝐁|(𝐔𝑗−𝐔𝑖) (3.19) converts the second-order central scheme into the first-order upwind scheme. Hence, this is a form of artificial dissipation, actually, one of the simplest. However, it is important to bear in mind that in the process, the order of the central scheme is being reduced to first-order, which is usually too low for the desirable accuracy of the solution. In order to maintain second-order accuracy, higher-order artificial dissipation terms have to be used, which, once again, should be based on the foundations of high-order upwind discretizations. Taking advantage of the definition of the Jacobian along an arbitrary direction, as given by equation (2.11), the components of the dissipation term associated to the face fluxes 𝐅𝑖𝑗 and 𝐆𝑖𝑗 can be grouped as 𝐃𝑛𝑖𝑗=𝑛𝑥𝑖𝑗𝐃𝑥𝑖𝑗+𝑛𝑦𝑖𝑗𝐃𝑦𝑖𝑗=12|𝐀𝑛|(𝐔𝑗−𝐔𝑖) (3.20) where 𝑛𝑥𝑖𝑗 and 𝑛𝑥𝑖𝑗 are the components of the unit vector 𝐧𝑖𝑗 normal to the face shared by the cells 𝑖 and 𝑗. Recalling that the normal face flux is then 𝐅𝑛𝑖𝑗=𝑛𝑥𝑖𝑗𝐅𝑖𝑗+𝑛𝑦𝑖𝑗𝐆𝑖𝑗 (3.21) the two-dimensional first-order upwind approximation can be expressed as follows Iván Padilla Montero 25 it to be equivalent to the low-order artificial dissipation scheme presented previously. This nonlinear switching behavior serves to comply with the restrictions imposed by Godunov’s theorem, thus allowing a sharp capturing of discontinuities with no oscillations. The reconstruction method chosen here has the advantage of being computationally efficient since the coordinates of the dummy nodes do not change with time for a fixed mesh. Hence, during practical implementation, the procedure to find the cells 𝑖𝑖 and 𝑗𝑗 in which each dummy node falls only has to be executed once, namely, at the beginning of the numerical solution. Note also that in the case of boundary cells, the dummy nodes may fall outside the computational domain. In order to handle this situation, the conservative variables for a point located outside the domain are extrapolated from the interior cells following a constant difference extrapolation, that is 𝐔𝑖𝑖=2𝐔𝑖−𝐔𝑗 𝐔𝑗𝑗=2𝐔𝑗−𝐔𝑖 (3.45) The pressures 𝑝𝑖𝑖 or 𝑝𝑗𝑗 to be used in the sensor can then be directly obtained from these extrapolated conservative variables. The artificial dissipation scheme presented in this section, which can also be implemented in both a scalar or matricial manner, is often also referred to as JST artificial dissipation scheme, or Jameson’s artificial dissipation. It has been widely applied to the solution of compressible flows with good results (Hirsch, 1990; Swanson & Turkel, 1990). 3.3 Numerical treatment of boundary conditions The necessary physical boundary conditions for external inviscid flow were described in section 2.3. In practice, their numerical implementation should be done with care in order to obtain a meaningful numerical solution. The approaches that have been adopted in this study are described next. 3.3.1 Body surface boundary condition The body surface is assumed to be a nonporous solid wall, so the normal velocity component must be zero at any point in the surface in order to satisfy the flow tangency condition. This boundary condition is introduced naturally in the finite volume method by forcing the flux vector across an face 𝑤 located at the wall to be 𝐅𝑛𝑤=[ 0 𝑛𝑥𝑤𝑝𝑤 𝑛𝑦𝑤𝑝𝑤 0] (3.46) Iván Padilla Montero 26 where 𝑛𝑥𝑤 and 𝑛𝑦𝑤 are the components of the face unit normal vector 𝐧𝑤, pointing into the wall, and 𝑝𝑤 is the pressure at the wall. The boundary face 𝑤 belongs to a single cell 𝑖𝑏 located at the boundary, which also has two interior faces shared with other cells of the mesh. Then, it is important to recognize that from the discretization given by equation (3.4), the value of the flow variables in the boundary cell, 𝐔𝑖𝑏, will be a result from the contribution of each of the three faces belonging to the cell, and as a result the boundary condition is only enforced in a weak form. Actually, the fluid at the centroid of the cell 𝑖𝑏 does not need to have a tangent velocity to the wall because it is not located at the wall. For consistency with the piecewise constant approximation, the fluid variables at the wall are assumed to be the same as that of the cell 𝑖𝑏, so that (3.46) becomes 𝐅𝑛𝑤=[ 0 𝑛𝑥𝑤𝑝𝑖𝑏 𝑛𝑦𝑤𝑝𝑖𝑏 0] (3.47) Other methods such as extrapolation or the solution of characteristic compatibility relations can be used in order to estimate the flow at the wall, see (Hirsch, 1990) for a general description. Note also that for the calculation of the flux vector at the boundary face no artificial dissipation is being used. This helps imposing the boundary condition in a more effective way, since the dissipation tends to excessively diffuse the important gradients that are found near the body boundary, reducing the accuracy of the solution. The wall boundary condition given by (3.47) can also be thought as a central approximation between a “mirror” cell located at the other side of the wall and the boundary cell 𝑖𝑏. This means that this formulation is also valid as a boundary condition for planes of symmetry in symmetric flow fields. Although this boundary condition is theoretically valid for any inviscid flow, it has to be taken into account that at the beginning of the transient simulation, the impulsive start from the initial conditions can lead to very large gradients in the wall, causing numerical instabilities. This is especially true for the case of high-speed flows using a uniform flow field as initial condition. In order to prevent this from occurring, the wall boundary condition can be imposed in a relaxed way, which can be accomplished by defining the following corrected velocity (Lyra & Morgan, 2002; Flores, Ortega, & Oñate, 2011; Ortega, 2014) 𝐕𝑐𝑜=𝐕𝑖𝑏−𝜅(𝐕∙𝐧𝑤)𝐧𝑤 (3.48) where 𝜅 is a parameter that has a value of zero at the start of the simulation and is progressively ramped up to one after a given number of time steps. Then, the solution is allowed to penetrate the wall at the start but, as time evolves, the normal velocity at the wall goes to zero. With this velocity correction, the body surface boundary flux eventually becomes Iván Padilla Montero 27 𝐅𝑛𝑤=𝑛𝑥𝑤 [ 𝜌𝑖𝑏𝑢𝑐𝑜 𝜌𝑖𝑏𝑢𝑐𝑜 2+𝑝𝑖𝑏 𝜌𝑖𝑏𝑢𝑐𝑜𝑣𝑐𝑜 𝜌𝑖𝑏𝑢𝑐𝑜ℎ0𝑖𝑏 ] +𝑛𝑦𝑤 [ 𝜌𝑖𝑏𝑣𝑐𝑜 𝜌𝑖𝑏𝑢𝑐𝑜𝑣𝑐𝑜 𝜌𝑖𝑏𝑣𝑐𝑜 2+𝑝𝑖𝑏 𝜌𝑖𝑏𝑣𝑐𝑜ℎ0𝑖𝑏 ] (3.49) where the total specific enthalpy ℎ0𝑖𝑏 is also computed based on the corrected velocity. 3.3.2 Far-field boundary condition As explained before, the numerical treatment of the far-field boundary condition has to account for the propagation of wave components across the boundary of the computational domain. Only those components which travel towards the interior of the domain can be prescribed, whereas the ones moving outwards have to be determined from the interior solution. Furthermore, as seen previously, the direction of propagation of one of the acoustic components changes depending on whether the flow is locally subsonic or supersonic. One way to achieve the correct setup of the far-field boundary condition is to work with characteristic variables, i.e., changing from conservative to characteristic variables, prescribing the necessary components as a function of the local direction of propagation, and transforming back to conservative variables. However, this is not desirable in practice since it is a computationally expensive procedure. Fortunately, the same result can be achieved by means of the positive flux Jacobian, following the principles of upwind schemes discussed before. Denoting by 𝐅𝑛𝑓 the flux vector across a far-field boundary face, the farfield boundary condition can be weakly enforced as (Lyra & Morgan, 2002; Flores, Ortega, & Oñate, 2011; Ortega, 2014) 𝐅𝑛𝑓=12(𝐅𝑛𝑖𝑏+𝐅𝑛∞)−12|𝐀𝑛 𝑅𝑜𝑒|(𝐔∞−𝐔𝑖𝑏) (3.50) where the superscripts 𝑖𝑏 and ∞ respectively denote the flow conditions at the far-field boundary cell and at the freestream. The upwind correction |𝐀𝑛 𝑅𝑜𝑒|(𝐔∞−𝐔𝑖𝑏) is evaluated at the Roe averages between the states at 𝑖𝑏 and ∞ by means of the algorithm described in section 3.2.1, but in this case no limitation is applied on the eigenvalues. This can be easily done by setting 𝜈𝑙=0 and 𝜈𝑛𝑙=0 at the corresponding far-field faces. To ensure that the condition given by (3.50) works properly, it is important for the far-field boundary to be located enough distance away from the body, so as to have near freestream conditions at the far-field boundary cells. With this formulation, the correct propagation of information is accounted for in the solution, without the need of directly prescribing any characteristic variables. This should guarantee that no wave reflections are produced at the far-field boundary, thus minimizing undesirable Iván Padilla Montero 28 perturbations in the flow field. Moreover, using this method avoids the need of differentiating between the subsonic/supersonic or inlet/outlet portions of the boundary. 3.4 Time integration Once the spatial discretization presented in (3.4) is complete, a semi-discrete scheme is obtained in the form of a system of nonlinear ordinary partial differential equations. The residual 𝐑𝑖 of a given cell 𝑖 can now be defined as 𝐑𝑖=1 𝐴𝑖∑𝐅𝑛𝑒 𝑁𝑒 𝑒=1 𝑙𝑒=−𝜕𝐔𝑖 𝜕𝑡 (3.51) which expresses the flux balance over the cell faces per unit area. There are many different options available for the time integration of a system defined by (3.51). Since our purpose is to obtain steady state solutions, the drivers for choosing a time integration scheme here are robustness and a fast convergence rate. Attending to these considerations, an explicit multistage Runge-Kutta scheme has been selected, given by 𝐔𝑠𝑖=𝐔𝑛𝑖−𝛼𝑠Δ𝑡𝑖 𝐑𝑖(𝐔𝑠−1 𝑖) 𝐔𝑛+1 𝑖=𝐔𝑛𝑠 𝑖 (3.52) for 𝑠=1,…,𝑛𝑠, with 𝑛𝑠 being the number of stages of the scheme. The subscripts 𝑛 and 𝑛+1 respectively denote the current and the next time levels, Δ𝑡𝑖 stands for the local time step associated to cell 𝑖, and 𝛼𝑠 are coefficients that depend on the number of stages employed. Note that at each stage, the residual has to be evaluated, as expressed by 𝐑𝑖(𝐔𝑠−1 𝑖). This scheme was originally introduced in the calculation of inviscid flows by (Jameson, Schmidt, & Turkel, 1981) for central schemes with artificial dissipation, and is mainly designed to allow the use of relatively high Courant numbers, in some cases higher than one. The most commonly used option, which is also the one adopted here, corresponds to the choice (Hirsch, 1990; Löhner, 2008; Lyra & Morgan, 2002) 𝛼1=14 , 𝛼2=13 , 𝛼3=12 , 𝛼4=1 (3.53) that usually offers a good trade-off between the allowable time step and the computational cost per time iteration (Flores, Ortega, & Oñate, 2011). This method is one of the so-called minimal storage Runge-Kutta schemes, which only require one extra copy of the conservative variables at each time step, that is, 𝐔𝑛𝑖 and 𝐔𝑠𝑖 for each cell in the mesh. Iván Padilla Montero 29 In order to minimize the computational cost of the scheme, a common practice is not to update the artificial dissipation terms at each stage of the scheme, but only at specific ones. A usual strategy is to calculate the dissipation terms only at the first stage, although this can cause some instabilities at the beginning of the transient when the flow properties change significantly between stages, especially if a low order dissipation scheme is being used. To increase the robustness of the solution, in this work the artificial dissipation is updated at the first and third stages of the time integration scheme. In practical implementations, the numerical solution is advanced in time until a specific convergence criterion is satisfied. A suitable criterion for the convergence of inviscid compressible flows is based on the norm of the density residual, see for instance (Hirsch, 2007). When the norm of the residual decreases by a given number of orders of magnitude, the numerical solution is stopped and a steady state is assumed. The 𝐿2 norm, or Euclidean norm, is often considered for that purpose, defined as ‖𝐑𝜌‖=√∑(𝑅𝜌𝑖)2 𝑁𝑐 𝑖=1 (3.54) where 𝑅𝜌𝑖 is the density residual of cell 𝑖, as given by the first component of the residual vector 𝐑𝑖, and 𝑁𝑐 is the number of cells in the mesh. Another relevant topic related with time integration is the choice of the initial conditions. In the case of high-speed inviscid flows below the hypersonic regime, the way to proceed is almost always to assume the initial flow-field to be the same as the freestream conditions. This creates strong gradients near the walls at the beginning of the time marching, which can in turn produce severe numerical instability. However, as discussed in the previous section, and as will be seen further below, some strategies have been devised in order to control such impulsive start. On the other side, it is interesting to comment that in the case of hypersonic flows, the availability of local surface inclination methods such as the Newtonian theory gives the ability of setting up much more appropriate initial conditions on common geometries. More details can be found in (Anderson, 2006). 3.4.1 Calculation of the time step Due to the fact that the chosen time integration scheme is explicit, the maximum allowable time step is constrained by the Courant-Friedrichs-Levy (CFL) stability criterion, which for the Euler equations can be formulated as (Hirsch, 1990) Δ𝑡𝑖=𝐶𝐹𝐿 ℎ𝑖 (𝑉𝑖+𝑐𝑖) (3.55) Iván Padilla Montero 30 where 𝐶𝐹𝐿 is the allowable Courant number, 𝑉𝑖 and 𝑐𝑖 are the velocity magnitude and the speed of sound at cell 𝑖, and ℎ𝑖 is a characteristic size of the cell. A reasonable characteristic cell size is be taken to be the minimum cell height (Flores, Ortega, & Oñate, 2011), which for triangular cells is defined as ℎ𝑖=2𝐴𝑖 max(𝑙𝑒) (3.56) where 𝑙𝑒 is the length of cell side 𝑒, with 𝑒=1,…,3. Since we are not interested in the transient solution, a local time stepping is chosen in order to accelerate the convergence to steady state, so that each cell progresses at its maximum possible time step. Attending to the limitation given by the CFL condition in (3.55), it is reasonable to think that an implicit scheme will pay-off in the case of steady state solutions. Actually, it may do. However, its implementation is much more complex than the explicit case, and their flexibility to adapt to changes in the numerical scheme is lower. Besides being simpler to implement, explicit schemes also have the advantage of being straightforward to parallelize (Löhner, 2008). This is very attractive when considering future improvements in the numerical solution, and is one of the main reasons why the explicit option has been adopted here. 3.4.2 A relaxed update procedure to promote the positivity of thermodynamic variables Due to the impulsive start from freestream conditions, during the first time integration steps negative values of the thermodynamic variables can appear in the flow field, with the subsequent failure of the numerical solution. To prevent the appearance of these local, spurious negative values during the convergence process and increase the robustness of the solution, the density and pressure can be updated using a relaxation process to help them remain positive. One of such relaxed update procedures is given by (Lyra & Morgan, 2002) 𝑝𝑛+1 𝑖=𝑝𝑛𝑖+Δ𝑝𝑖[1+𝜂(𝜃−Δ𝑝𝑖 𝑝𝑛𝑖)]−1 (3.57) whenever Δ𝑝𝑖𝑝𝑛𝑖 ⁄≤𝜃, with the recommended values being 𝜃=−0.2 and 𝜂=2. The increment Δ𝑝𝑖 simply denotes Δ𝑝𝑖=𝑝𝑛+1 𝑖−𝑝𝑛𝑖. The same procedure is then applied for the density 𝜌. The use of this special update method (and similar ones) is found to play a significant role in the solution of high-speed inviscid flows; see (Yee, Klopfer, & Montagné, 1988). Iván Padilla Montero 31 4 NUMERICAL IMPLEMENTATION AND TEST CASES 4.1 Description of the solver developed With the objective of putting in practice all the theoretical concepts presented before, a numerical solver has been developed using Fortran. This constitutes a direct implementation of the unstructured finite volume method for the solution of the Euler equations. Fortran has been chosen for its high performance in scientific computing and its simplicity of implementation. The code is designed to be able to solve two-dimensional steady inviscid flows around closed geometries (external flows), at speeds ranging from subsonic to moderate supersonic regimes, up to about Mach 3. The solution of higher Mach number flows has not been considered due to the complex physical phenomena that takes place when high-temperature effects become important, which would require including real gas effects in the solution (Anderson, 2006). In order to minimize the memory requirements and reduce the number of calculations needed, a side-based strategy has been taken, see for instance (Löhner, 2008). This approach consists on calculating the residual (3.51) by iterating over faces instead of cells. Any mesh considered has a given number of global faces, which, excluding the boundaries, are shared by two triangular cells each. Similarly, each triangular cell has three local faces. Then, for the numerical scheme chosen, the flux approximation at a global face shared between two cells contributes in the same amount to the residual of the each of the two cells, but with different sign. As a consequence, the use of a side-based solution almost halves the calculations otherwise required in the case of calculating the residual cell by cell. A diagram of the code structure is shown in Figure 4.1. As can be observed, the program follows a modular approach, which has been adopted to favor the possible implementation of future improvements in the solution. For the generation of unstructured triangular meshes, the pre and postprocessing tool GiD has been used, for which a problem type was created in order to automatically generate the input files for the solver. On the other end, MATLAB has been chosen for the analysis of results and the generation of different types of plots. The program starts by reading the input data files of the solution. These contain the node coordinates and the cell connectivities of the mesh, as well as the simulation parameters and the freestream flow conditions. Next, the necessary data structures are built, which include the definition of the global faces of the mesh and the calculation of necessary data for the reconstruction of the fourth difference stencils, to be used later in the evaluation of Iván Padilla Montero 32 the JST artificial dissipations terms. Once the face connectivities have been defined, a subroutine that calculates the geometrical data is executed, which produces the lists of cell areas, side lengths and normals, and the characteristic cell sizes required for evaluating the time step. Upon successfully reaching this point, the fluid variables are initialized and the Initialize fluid variables Read input data: Node coordinates Cell connectivities Freestream conditions Simulation parameters Build data structures: Calculate face connectivities Reconstruct fourth difference stencils Calculate geometry data: Cell area Face lengths Face unit normals Characteristic cell size Time loop: Calculate time step Integrate in time: oEvaluate residual oUpdate fluid variables Evaluate convergence criterion Write output results file and end program GiD: generate mesh input file MATLAB: postprocess results Figure 4.1 Diagram showing the different structural blocks and flow of the developed code. Iván Padilla Montero 33 time marching process starts. At each time iteration, the multistage scheme is executed, computing the residual and updating the variables at each stage. Then, the convergence criterion is evaluated. When convergence is achieved, the time loop stops and a function writes an output file with the results prepared to be postprocessed with MATLAB. Many different solution parameters can be controlled by the user, which mainly include: the stages of the time integration method, the Courant number, the target density residual for the convergence criterion and the different coefficients and form of the artificial dissipation terms. 4.2 Test cases A set of five different test cases have been solved in order to validate the numerical results obtained with the developed solver. The first three of them are based on the transonic and supersonic flow around a NACA0012 airfoil, which correspond to a series of rigorous airfoil benchmark tests performed by (Pulliam & Barton, 1985). The fourth case deals with the supersonic flow past a double wedge airfoil, which is an interesting benchmark problem since it has an analytical solution, and the fifth and last case tackles the difficult problem of the inviscid supersonic flow past a circular cylinder. A summary of the different cases, specifying the flow configuration for each one, can be found in Table 4.1. Case ID Geometry 𝑀∞ 𝛼 (degrees) 1 NACA0012 0.85 1 2 NACA0012 0.95 0 3 NACA0012 1.2 7 4 Double wedge airfoil 2 0 5 Circular cylinder 3 - Table 4.1 Summary of the different test cases considered. The reference length 𝐿 for these cases is taken to be the chord of the body, which is defined as the linear distance between the leading and trailing edges. For convenience, a unit value of the chord is adopted in all five cases, so 𝐿=1. There are some other parameters that are also common to all the examples. On one side, the ratio of specific heats is always chosen to be the standard value 𝛾=1.4. On the other side, the convergence criterion is always set to be a reduction of 8 orders of magnitude in the density residual with respect to its initial value. Moreover, as commented before, a four- Iván Padilla Montero 34 stage integration scheme is selected in all the tests, with the coefficients given by (3.53), in which the artificial dissipation terms are only calculated at the first and third stages. Note that thanks to the use of the nondimensional equations, for a fixed value of 𝛾 only the freestream Mach number and the angle of attack are needed to completely define the flow in each case. Regarding the body surface boundary condition (refer to section 3.3.1), the parameter 𝜅 is implemented so as to increase with the number of time iterations 𝑖𝑡 following 𝜅=1−0.9𝑖𝑡 (4.1) which can be assumed to reach a value of one at about 100 iterations. 4.2.1 Calculation of nondimensional quantities for the analysis of results In order to properly validate the numerical results obtained, there are different important nondimensional quantities that should be verified. The analysis of nondimensional results is more general and offers a simpler and unambiguous framework to perform comparisons against reference results. In this study, the quantities considered are the Mach number, the pressure coefficient, the nondimensional change in entropy and the aerodynamic force and moment coefficients. Recalling the nondimensional variables defined in section 2.4, the Mach number can be calculated from the results of the simulation as 𝑀=√𝑢2+𝑣2 𝑇 (4.2) Similarly, the following expression can be derived for the pressure coefficient 𝐶𝑝=2(𝑝−1𝛾) 𝑀∞ 2 (4.3) In a similar fashion, the change in entropy for a calorically perfect gas can be given by Δ𝑠=𝑠−𝑠∞=𝑐𝑝ln𝑇−𝑅ln(𝛾𝑝) (4.4) which can be made dimensionless using for example the specific heat at constant pressure, that is Δ𝑠=Δ𝑠 𝑐𝑝=ln𝑇−𝛾−1 𝛾ln(𝛾𝑝) (4.5) Iván Padilla Montero 41 Figure 4.4 Pressure coefficient and Mach number results for case 2. Iván Padilla Montero 42 Figure 4.5 Pressure coefficient and Mach number results for case 3. Iván Padilla Montero 43 Figure 4.6 Convergence history and entropy generation for cases 1 (top), 2 (middle) and 3 (bottom). Iván Padilla Montero 44 4.2.5 Supersonic flow past a double wedge airfoil This fourth test case aims to the solution of the supersonic flow past a double wedge airfoil, also known as a diamond-wedge airfoil. The wedge angle is 15 degrees, which forces the flow to turn a value of 30 degrees at the expansion points located at half of the chord, and the freestream Mach number and angle of attack are respectively 𝑀∞=2 and 𝛼=0º. This case is attractive as a benchmark problem because it has an analytical solution, and contains both shock and expansion waves. The analytical solution is given by the RankineHugoniot relationships and the Prandtl-Meyer function, see for instance (Anderson, 2011) for the specific details. The mesh considered for this problem is sketched in Figure 4.9, which contains 23778 cells and 12076 nodes. Only the cells near the body are shown because the rest of the domain is equal the one for the NACA0012 airfoil discussed before. The outer boundary is once again 50 chords away from the airfoil, and some refinement is applied in the leading edge and the expansion points. This time, the matrix form of the JST artificial dissipation terms have been used, which improve the accuracy of the solution. The parameters chosen are 𝐶𝐹𝐿=0.6, 𝑘2=1, 𝑘4=1/16, 𝜈𝑙=0.05 and 𝜈𝑛𝑙=0.1. A summary of the numerical results is provided in Figure 4.7 and Figure 4.8. As can be observed, the surface quantities show a good agreement with the analytical solution. The pressure coefficient is perfectly predicted, but the Mach number after the expansion wave presents a small deviation. It has been found in numerical experiments that at the expansion point, the numerical scheme in use tends to introduce an excessive amount of artificial dissipation at the surface, affecting only the surface Mach number predicted by the solution (the values computed far from the surface are correct). Actually, the use of a scalar artificial dissipation for this problem leads to significant errors in the Mach number behind the expansion wave. On the other hand, some small oscillations are also present at the shock and expansion waves. Although the JST scheme mimics a first-order upwind scheme in the vicinity of discontinuities, it does not lead to the same exact formulation as the first-order artificial dissipation scheme given in section 3.2.1, and use of a more elaborated pressure sensor may lead to better results. It has been checked that with the use of the first order scheme, these oscillations disappear. The aerodynamic force coefficients also have been calculated for this problem, and can be found in Table 4.5. The values are close to the analytical solution, demonstrating that the overall quality of the solution is satisfactory. 𝑐𝑙 𝑐𝑑 Numerical 0.0003 0.1704 Analytical 0 0.1715 Table 4.5 Comparison of aerodynamic force coefficients for case 4. Iván Padilla Montero 45 Figure 4.7 Pressure coefficient and Mach number results for case 4. Iván Padilla Montero 46 Figure 4.8 Convergence history and entropy change for case 4. Figure 4.9 Details of the mesh used for the numerical solution of case 4. Iván Padilla Montero 47 4.2.6 Supersonic flow past a circular cylinder The last test case that has been attempted is the solution of the complete supersonic flow field around a circular cylinder. The configuration of this case is based on a solution that was performed by (Lyra & Morgan, 2002), with a freestream Mach number of 𝑀∞=3. This problem is challenging in terms of stability behavior due to the presence of a strong bow shock close to the cylinder surface and a complex, rarefaction zone at back. The strong perturbations created by a very blunt object such as the circular cylinder, combined with the relatively high Mach number considered creates very strong gradients in the solution which are usually difficult to manage. As a result, this constitutes a firm benchmark to test the robustness of the numerical solution. It is important, however, to be aware that the solution of this case through the Euler equations is by no means realistic. In the real case, the flow field presents a large detachment region at the back of the cylinder, so the consideration of viscous phenomena is essential for reproducing the correct physical behavior. As usual, a mesh with a radius of 50 chords has been constructed for this case, as represented in Figure 4.12. The cell size is on the same order as the previous cases, resulting in 9999 nodes and 19779 cells. The reference solution is based on a finite element high-order upwind scheme, which is quite different from the approach considered in this analysis. Different numerical trials have been performed with the code developed. However, due to the difficulty of this case, only the solution using the matrix low-order artificial dissipation scheme has been obtained. The application of the JST dissipation terms has always resulted in bad convergence behavior and computational instability. The solution obtained here is presented in Figure 4.10 and Figure 4.11. It was computed with 𝐶𝐹𝐿=0.5, 𝜈𝑙=0.2, 𝜈𝑛𝑙=0.2 and the use of the special update procedure described in section 3.4.2. This, along with the relaxed imposition of the body boundary condition, have been found to be of significant importance in order to be able to advance the solution in time at the first steps. Observing the results, it can be seen that despite the use of a low-order scheme, the important flow field features are captured with an acceptable resolution. Actually, although the solution is over-diffusive due to the low-order upwinding, the obtained pressure coefficient distribution agrees very well with the values found in the reference calculation. Both the strength and location of the detached shock are correct, as well as the predicted pressure values at the back of the cylinder. Regarding the Mach number, the jump across the bow shock is well captured. However, in a similar way as in the previous case, the large expansion that takes place at the second half of the cylinder is not accurately resolved. The strong shock wave that is encountered at the rear part is in the correct location, but the surface Mach number just before the shock is about 4.5, whereas in the reference solution the value reached is only 4. In their studies, (Yee, Klopfer, & Montagné, 1988) and (Lyra & Morgan, 2002) suggest that the use a more elaborated entropy fix than that given by equation (3.36) is required when the freestream Mach number starts to become relatively high. This may be one of the possible explanations why an excessive amount of artificial dissipation is added to the solution in the strong expansion region. Iván Padilla Montero 48 Figure 4.10 Pressure coefficient and Mach number results for case 5. Iván Padilla Montero 49 Figure 4.11 Convergence history and entropy change for case 5. Figure 4.12 Detail views of the mesh for case 5. Iván Padilla Montero 50 5 CONCLUSIONS AND FUTURE WORK In summary, the main objective of this work has been achieved upon the obtention of satisfactory results with the numerical solver developed. Due to the hyperbolic mathematical nature of the Euler equations, the fact of taking into account the propagation of information in the fluid is found to be of paramount importance in order to obtain meaningful numerical results. The wave components that travel through the flow field have been described and analyzed, emphasizing the relevance of carefully imposing the boundary conditions and justifying the need of using artificial dissipation when central numerical schemes are considered. The unstructured finite volume discretization described has been successfully applied to all the test cases considered, demonstrating that the decisions taken regarding the piecewise constant approximation of values inside the cells, and the implementation of boundary conditions are perfectly valid for the solution of inviscid compressible flows. Besides, the use of the Jameson-Schmidt-Turkel artificial dissipation terms has provided satisfactory results in most of the test cases considered, which also confirms the feasibility of the stencil reconstruction technique adopted and the choices made regarding the calculation of the positive flux Jacobian. Focusing on the results obtained, it has been found that the overall numerical solution has a very good performance in the three NACA0012 airfoil examples calculated, showing an excellent agreement with the reference results. A good solution has also been achieved for the supersonic flow past a double wedge airfoil, matching the analytical results with very reasonable accuracy, although some small deviation is found in the Mach number behind the expansion waves. Finally, for the complete supersonic flow around a circular cylinder, only the solution with a low-order artificial dissipation scheme could be obtained. The numerical difficulties associated to the large perturbations created by the cylinder and the presence of a strong bow shock and rarefaction zones make this problem quite challenging. In order to obtain acceptable results for this case, the use of a more robust numerical scheme should be considered, such as a high-order upwind method. Significant discrepancies are also observed in this case with respect to the Mach number predicted in the large expansion region at the second half of the cylinder. This, similarly to the expansion in the double wedge airfoil case, appears to be caused by the introduction of excessive artificial dissipation in the surface. A possible explanation for this behavior may be the need for a better entropy fix, and is an aspect that should be further investigated. Other improvements of the numerical solver can be considered for future implementation, which mainly include mesh adaptation, parallelization, the extension to three dimensional problems and the introduction of a thermally perfect gas behavior, to allow modeling the variation of specific heats with temperature in the case of higher Mach number flows.