Numerical Resolution of Supersonic Flow using the MacCormack Method

0 Aerospace Engineering and Fluid Mechanics Department Superior Technique School of Engineering University of Seville Seville, 2018 Author: Antonio González Carvajal Tutor: Miguel Pérez-Saborid Sánchez-Pastor Numerical Resolution of Supersonic Flow using the MacCormack Method Senior Research Project Degree in Aerospace Engineering 1 2 Resolución numérica de flujo supersónico usando el método de MacCormack Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Trabajo de Fin de Grado Grado en Ingeniería Aeroespacial Autor: Antonio González Carvajal Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor Titular 3 4 Trabajo de Fin de Grado: Numerical Resolution of Supersonic Flow using the MacCormack’s Technique Autor: Antonio González Carvajal Tutor: Miguel Pérez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2018 El Secretario del Tribunal 5 6 A mi familia, a mis amigos y especialmente a mis abuelas y a Lorena. 7 8 Agradecimientos Han sido varios meses los que le he dedicado a este trabajo y tengo que decir que a pesar de todo mi esfuerzo, no habría podido conseguirlo sin muchas personas. Porque a veces he estado al borde de un ataque de nervios, al borde de desistir, tirar la toalla, no verme capaz, y ellos me han sado fuerzas para seguir intentándolo. Este trabajo no habría sido terminado sin su apoyo. En primer lugar mi familia, que siempre ha estado ahí dándome ánimos no solo en el tiempo que he estado trabajando en el proyecto, sino durante toda la carrera. Han sido mis padres los que han conseguido calmarme cuando perdía los nervios, sentía que no avanzaba y casi quería llorar. Ha sido mi hermana la que me distraía y me hacía reír aun cuando tenía miles de preocupaciones en la cabeza. Y han sido mis abuelas las que siempre se preocupaban por mí y me han dado tanto, desde conversaciones tranquilizadoras con la perspectiva de la edad hasta platos de cuchareo para comer en familia. No solo a ellos tengo que agradecerles por supuesto. Tampoco habría conseguido terminar si no fuese por mis amigos. Mis amigos en Sevilla me recordaban en Julio que aún tenía tiempo y que vería las cosas mejor después de una cerveza en la Alameda. Amigos en Chipiona que consiguen hacer de los veranos una auténtica pasada. Gente tan maravillosa que a veces sientes que no las mereces de lo bien que lo pasas con ellos. A todos ellos, gracias. Os lo debo todo y más. 15 On the other hand, the detachment of the boundary layer happens when the pressure gradient in the boundary layer is against the direction of the flow, and is big enough to be stronger than the momentum transmitted to the inner layers from the outer flow by viscosity effects. When that happens, the inner layers stop moving and a recirculating area is created around which the rest of the flow slips. 1.3 Structure of the project This project will apply the numerical technique developed by McCormack into a program to solve supersonic problems with different geometries and conditions. It will begin after this introduction by explaining in the second chapter the theoretical knowledge that is needed to understand the project. First, the equations of the flow that have been implemented, followed by an explanation of the method used for uniform grids. After that, it will be explained the changes in the method necessary to develop the program for non-uniform grids. The third chapter will be focused on explaining the program that has been written section by section, first for the uniform grids, and after that for the non-uniform grids. At the end of the explanation of each section the program used in each case will be explicitly showed. Figure 1.4 . In the smoke of a candle it can be appreciated the jump from laminar to turbulent behaviour Figure 1.5 . Boundary layer detachment around a wing profile. 16 In the fourth chapter, the results of the different applications will be showed, beginning with the flat plate. After the plate, the forward and rearward steps will be explained, followed by the rearward slope, where it will be analysed the Prandtl-Meyer Expansion Waves. The last application will be the study of a leading edge and the shock wave created around it. Finally, the project will end with a brief discussion about conclusions and further developments in the fifth chapter, the references in the chapter six, and the summary in spanish language demanded by the University of Seville normative in Appendix 1. 17 2 Theoretical basis From the multiple methods proposed by CFD, in this project it has been chosen a time marching calculation using the MacCormack’s technique, for several reasons. It was first introduced by Stanford University professor Robert W. MacCormack in 1969. Among its advantages there is that it is one of the easiest and more understandable methods that allows the solving of the complete Navier-Stokes equations, which is the goal of this project, even though it is still as complex as expected for a full-on implementation of this equations. The current chapter is going to explain the method step by step as it has been used in this project, and after that, the following chapter will explain the coding implementation by going through the program used to obtain the results. 2.1 CDF equations If the resolution of the complete Navier-Stokes equations by numeric means is required, first an adaptation is recommended to make the equations more manageable. To begin with, the equations for a bidimensional flow are remembered here.  Conservation of mass: 𝜕𝜌 𝜕𝑡+𝜌𝛻⋅𝑉 󰇍  =0 Equation 2.1  Conservation of momentum in the x and y directions: 𝜕(𝜌𝑢) 𝜕𝑡 +𝛻⋅(𝜌𝑢⋅𝑉 󰇍  ) = −𝜕𝑝 𝜕𝑥+𝜕𝜏 𝜕𝑥 +𝜕𝜏 𝜕𝑦 +𝜌𝑓 Equation 2.2 𝜕(𝜌𝑣) 𝜕𝑡 +𝛻⋅(𝜌𝑣⋅𝑉 󰇍  ) = −𝜕𝑝 𝜕𝑦+𝜕𝜏 𝜕𝑥 +𝜕𝜏 𝜕𝑦 +𝜌𝑓 Equation 2.3  Conservation of energy: 𝜕 𝜕𝑡󰇩𝜌󰇧𝑒+𝑉 2󰇨󰇪+𝛻⋅󰇩𝜌󰇧𝑒+𝑉 2󰇨⋅𝑉 󰇍  󰇪 = 𝜌𝑞󰇗+𝜕 𝜕𝑥𝑘𝜕𝑇 𝜕𝑥+𝜕 𝜕𝑦𝑘𝜕𝑇 𝜕𝑦−𝜕(𝑢𝑝) 𝜕𝑥 −𝜕(𝑣𝑝) 𝜕𝑦 +𝜕(𝑢𝜏) 𝜕𝑥 +𝜕(𝑢𝜏) 𝜕𝑦 +𝜕(𝑣𝜏) 𝜕𝑥 +𝜕(𝑣𝜏) 𝜕𝑦 +𝜌𝑓⋅𝑉 󰇍  Equation 2.4 18 Just with a first glance of the equations, it can be recognised a simpler, more compacted way of expressing them by grouping together the derivative terms with respect to the same variable. The result would be an expression with the following form: 𝜕𝑈 𝜕𝑡+𝜕𝐹 𝜕𝑥+𝜕𝐺 𝜕𝑦=𝑆 Equation 2.5 Where U, F, G and S are vectors with one component for each of the equations that govern the behaviour of the flow. For the bidimensional flow that is being discussed, their expressions are shown here: 𝑈= ⎩ ⎪ ⎨ ⎪ ⎧ 𝜌 𝜌𝑢 𝜌𝑣 𝜌󰇧𝑒+𝑉 2󰇨 ⎭ ⎪ ⎬ ⎪ ⎫ Equation 2.6 𝐹= ⎩ ⎪ ⎨ ⎪ ⎧ 𝜌𝑢 𝜌𝑢+𝑝−𝜏 𝜌𝑣𝑢−𝜏 𝜌󰇧𝑒+𝑉 2󰇨𝑢+𝑝𝑢−𝑘𝜕𝑇 𝜕𝑥−𝑢𝜏−𝑣𝜏 ⎭ ⎪ ⎬ ⎪ ⎫ Equation 2.7 𝐺= ⎩ ⎪ ⎨ ⎪ ⎧ 𝜌𝑣 𝜌𝑢𝑣−𝜏 𝜌𝑣+𝑝−𝜏 𝜌󰇧𝑒+𝑉 2󰇨𝑣+𝑝𝑣−𝑘𝜕𝑇 𝜕𝑦−𝑢𝜏−𝑣𝜏 ⎭ ⎪ ⎬ ⎪ ⎫ Equation 2.8 𝑆= ⎩ ⎪ ⎨ ⎪ ⎧ 0 𝜌𝑓 𝜌𝑓 𝜌𝑢𝑓+𝑣𝑓+𝜌𝑞󰇗 ⎭ ⎪ ⎬ ⎪ ⎫ Equation 2.9 This way it can be obtained a particularly suitable expression to manage with CFD methods. This expression is often called the “Conservation Form”, and is extensively used throughout several CFD methods, not only in the case of the MacCormack’s technique. In the equations Equation 2.7 and Equation 2.8, an expression for 𝜏 , 𝜏 and 𝜏 is needed to successfully evaluate the variable vectors F and G. Assuming a newtonian fluid, these are given by the expressions obtained by Stokes for this kind of fluids: 19 𝜏=𝜆∇⋅𝑉 󰇍  +2𝜇𝜕𝑢 𝜕𝑥 Equation 2.10 𝜏=𝜆∇⋅𝑉 󰇍  +2𝜇𝜕𝑣 𝜕𝑦 Equation 2.11 𝜏=𝜏=𝜇𝜕𝑣 𝜕𝑥+𝜕𝑢 𝜕𝑦 Equation 2.12 In addition, it is going to be used the hypothesis made also by Stokes about the value of the volumetric viscosity coefficient, 𝜆=−𝜇. So, the expressions Equation 2.10 and Equation 2.11 change into: 𝜏=23𝜇2𝜕𝑢 𝜕𝑥−𝜕𝑣 𝜕𝑦 Equation 2.13 𝜏=23𝜇2𝜕𝑣 𝜕𝑦−𝜕𝑢 𝜕𝑥 Equation 2.14 With these expressions, the system consists of eight variables: 𝜌,𝑢,𝑣,𝑝,𝑇,𝑒,𝜇 and 𝑘; but only has four equations. Thus, we must provide four more equations to be able to solve the system, counting on the fact that we consider 𝑉 not an independent variable, as it can be easily obtained from the velocity components 𝑢 and 𝑣. The first of them is the assumption of a perfect gas, from which the perfect gas equation is derived: 𝑝=𝜌𝑅𝑇 Equation 2.15 Where R is the perfect gas constant of value 287 J/ (kg K). The second is an assumption built on top of the previous one, the calorically perfect gas. The relation that describes this assumption is the following: 𝑒=𝑐𝑇 Equation 2.16 Where 𝑐 is the specific heat at constant volume and is, by the calorically perfect definition, considered constant. Given the value of the ratio of specific heats, 𝛾, equal to 1.4, it is obtained with the following relation: 𝑐=𝑅 𝛾−1 20 The next one is a model for the viscosity provided by the Sutherland’s law: 𝜇=𝜇𝑇 𝑇/𝑇+110 𝑇+110 Equation 2.17 In which 𝜇 and 𝑇 are reference values taken at standard sea level. The last equation is given by assuming the Prandtl number constant and approximately equal to 0.71 in the case of air. From the Prandtl number definition it is obtained: 𝑃𝑟=𝜇𝑐 𝑘=0.71 Equation 2.18 Where 𝑐 is the specific heat at constant pressure and is easily obtained as: 𝑐=𝛾𝑐 21 2.2 MacCormack’s Method for Uniform Grids The MacCormack’s technique relies on transforming the governing equations into the conservative form and then performs a two-phase discretization and combines them to assure second order accuracy, as it is going to be demonstrated in this subchapter. From Equation 2.5, by isolating the time derivative in one side, it is obtained: 𝜕𝑈 𝜕𝑡=𝑆−𝜕𝐹 𝜕𝑥−𝜕𝐺 𝜕𝑦 Equation 2.19 Now, the flow domain is discretised into points following an equispaced rectangular grid, with each point being designated by its coordinates i and j. The number of points in the X and Y axis are Nx and Ny respectively, as shown in Figure 2.1. For each time step, the properties of the flow are calculated by finite difference approximation only for the inner points of the grid. The properties of the flow in the boundary must be treated differently and obtained by other means. That is, given by the boundary conditions or extrapolated from the flow properties inside the domain. The boundary conditions that have previously been mentioned have to be explicitly enforced in the method. However, the Forward-Rearward steps that are going to be explained later offer as their result the vector U. Moreover, the physical variables, that are going to be referred as “primitive” from this point on, also have to be obtained to assess the final results. The boundary conditions could be enforced by combining the enforced primitive variables and by combining U values when an extrapolation is needed. Despite that, in this project it has Figure 2.1 . Computational domain and Uniform Grid 22 been chosen the more straightforward and comprehensible approach of decoding the primitive variables in each step, then enforcing the boundary conditions on them, and after that, recalculating the vectors U, F, G and S. So, after each step, the primitive variables are decoded from U as follows: 𝜌=𝑈 Equation 2.20 𝑢=𝜌𝑢 𝜌=𝑈 𝑈 Equation 2.21 𝑣=𝜌𝑣 𝜌=𝑈 𝑈 Equation 2.22 𝑒=𝑈 𝑈−𝑢+𝑣 2 Equation 2.23 Using the Equation 2.15 the pressure is obtained, and from the Equation 2.16: 𝑇=𝑒 𝑐 Equation 2.24 After that, the viscosity is calculated from Sutherland’s law, Equation 2.17, and using Equation 2.18: 𝑘=𝜇𝑐 𝑃𝑟 Equation 2.25 For the first step, starting with Equation 2.19, the derivatives are approximated using finite differences, with the values Δx, Δy and Δt substituting 𝜕𝑥, 𝜕𝑦 and 𝜕𝑡 respectively. Their appropriate values are going to be discussed in the next chapter. The magnitudes are projected in the forward direction for both space directions, and the result is the predicted value of the vector U at the time 𝑡+ 𝛥𝑡, which is denoted by 𝑈. The result expression, already solved for 𝑈, is: 𝑈, =𝑈, −∆𝑡 ∆𝑥𝐹, −𝐹, −∆𝑡 ∆𝑦𝐺, −𝐺, +∆𝑡⋅𝑆 Equation 2.26 Now, there is an important remark that has to be made. The vectors F and G are here finitely differenced, but they are themselves formed by derivative terms. These terms also have to be calculated using a finite difference approximation, but, to ensure the second order accuracy of 23 the method, these finite differences must be in the opposite direction that of the vector they form part of. If the differentiation is not in the same direction, then central difference is used. For example, in the forward step, the viscosity 𝜏, that forms part of the vector F, must be calculated using rearward difference in the X direction, and using central difference in the Y direction. Thus, the final expression for 𝜏 is the following: 𝜏,= 23𝜇, 󰇧2𝑢, −𝑢,  𝛥𝑥 −𝑣, −𝑣,  2 𝛥𝑦 󰇨 Equation 2.27 The same reasoning with 𝜏 conducts to: 𝜏,= 23𝜇, 󰇧2𝑣, −𝑣,  𝛥𝑦 −𝑢, −𝑢,  2 𝛥𝑥 󰇨 Equation 2.28 The case of 𝜏 is different. As it takes part in both F and G calculations, two different calculations of 𝜏 have to be made, one for each vector, with distinct differentiation in each case. Thus, the expressions reached for 𝜏 used in F and G are: 𝜏,󰇻=𝜇, 󰇧𝑣, −𝑣,  𝛥𝑥 +𝑢, −𝑢,  2 𝛥𝑦 󰇨 Equation 2.29 𝜏,󰇻=𝜇, 󰇧𝑣, −𝑣,  2 𝛥𝑥 +𝑢, −𝑢,  𝛥𝑦 󰇨 Equation 2.30 Also, the terms referring heat conduction with the Fourier’s law are approximated with finite differences: 𝜕𝑇 𝜕𝑥,≈𝑇,−𝑇, 𝛥𝑥 Equation 2.31 𝜕𝑇 𝜕𝑦,≈𝑇,−𝑇, 𝛥𝑦 Equation 2.32 It must be taken into account that these expressions cannot be used in the limits of the computational domain, as some of them would require points that are outside of the boundaries. So, for example, in the lower boundary the central and rearward differences done in the Y direction are substituted by forward differentiation. In the upper border, the forward and central are substituted with rearward ones, and the same logic applies to the left and the right boundaries for the X direction, respectively. As stated by John D. Anderson in his book 24 “Computational Fluid Dynamics” [Ref 1], doing this only compromises slightly the second order accuracy of the method. After calculating the vectors F and G, the Equation 2.26 is applied, the primitive variables are extracted as shown in equations Equation 2.20 to Equation 2.25 and the boundary conditions are enforced onto these variables. The second step proceeds in a similar manner. However, there is a clear difference. Now the vectors F and G used are calculated using the predicted values of the variables. The equivalent to Equation 2.26 now is the following expression: 𝑈, =𝑈, −∆𝑡 ∆𝑥𝐹, −𝐹, −∆𝑡 ∆𝑦𝐺, −𝐺, +∆𝑡⋅𝑆 Equation 2.33 Where 𝑈 is called the corrected value of U. As it can be seen in the Equation 2.33, the vectors 𝐹 and 𝐺 are obtained using the predicted values calculated from the previous step. Also, the finite differences are now projected rearward. Similarly to what happened in the forward step, the derivative terms that form part of 𝐹 and 𝐺 must be approximated using the opposite direction in their finite differences. So, the expressions for 𝜏 , 𝜏 , 𝜏 ,   and   used in the corrector step are similar to the expressions Equation 2.27 to Equation 2.32, but with different differentiation. To avoid excessive repetition in this theoretical chapter, these expressions are not explicitly stated. However, if the reader so chooses, their implementation can be found in the section “Encode variables for corrector step” of the script that appears at the end of subchapter 3.1. Exactly as it happened in the predictor step, two calculations of 𝜏 must be made, one for the vector 𝐹 and other for 𝐺. Also, the same considerations regarding the substitution of some types of differences with others in the boundaries of the computational domain apply. Once the value of 𝑈 is calculated, the primitive variables are decoded to apply the boundary conditions on them, and then the vector 𝑈 is recalculated with the boundary conditions already applied. Having the values of 𝑈 and 𝑈 from the predictor and corrector steps respectively, the value of 𝑈 is obtained as their arithmetic mean. 𝑈, =12𝑈, +𝑈,  Equation 2.34 Finally, the physical variables are obtained as the final result of the method at time 𝑡+ 𝛥𝑡. The cycle is then restarted, advancing in time and capturing the behaviour of the flow. 31 𝐽𝜕𝐺 𝜕𝜉𝜕𝜉 𝜕𝑦=𝜕 𝜕𝜉𝐽𝐺𝜕𝜉 𝜕𝑦−𝐺𝜕 𝜕𝜉𝐽𝜕𝜉 𝜕𝑦 Equation 2.63 𝐽𝜕𝐺 𝜕𝜂𝜕𝜂 𝜕𝑦=𝜕 𝜕𝜂𝐽𝐺𝜕𝜂 𝜕𝑦−𝐺𝜕 𝜕𝜂𝐽𝜕𝜂 𝜕𝑦 Equation 2.64 Substituting them into the Equation 2.60, the result is this other equation: 𝐽𝜕𝑈 𝜕𝑡+𝜕 𝜕𝜉𝐽𝐹𝜕𝜉 𝜕𝑥+𝐽𝐺𝜕𝜉 𝜕𝑦+𝜕 𝜕𝜂𝐽𝐹𝜕𝜂 𝜕𝑥+𝐽𝐺𝜕𝜂 𝜕𝑦−𝐹𝜕 𝜕𝜉𝐽𝜕𝜉 𝜕𝑥+𝜕 𝜕𝜂𝐽𝜕𝜂 𝜕𝑥 −𝐺𝜕 𝜕𝜉𝐽𝜕𝜉 𝜕𝑦+𝜕 𝜕𝜂𝐽𝜕𝜂 𝜕𝑦=𝐽⋅𝑆 The terms multiplied by F and G can be rewritten by the means of the Equation 2.57, and be this way proven equal to zero. 𝜕 𝜕𝜉𝐽𝜕𝜉 𝜕𝑥+𝜕 𝜕𝜂𝐽𝜕𝜂 𝜕𝑥 = 𝜕 𝜕𝜉𝜕𝑦 𝜕𝜂−𝜕 𝜕𝜂𝜕𝑦 𝜕𝜉 = 𝜕𝑦 𝜕𝜉𝜕𝜂−𝜕𝑦 𝜕𝜂𝜕𝜉 =0 𝜕 𝜕𝜉𝐽𝜕𝜉 𝜕𝑦+𝜕 𝜕𝜂𝐽𝜕𝜂 𝜕𝑦 = 𝜕 𝜕𝜉−𝜕𝑥 𝜕𝜂+𝜕 𝜕𝜂𝜕𝑥 𝜕𝜉 = 𝜕𝑥 𝜕𝜂𝜕𝜉−𝜕𝑥 𝜕𝜉𝜕𝜂 =0 So, the initial goal is reached, and the governing equation is: 𝜕𝑈 𝜕𝑡+𝜕𝐹 𝜕𝜉 +𝜕𝐺 𝜕𝜂=𝑆 Equation 2.65 Where the new vectors 𝑈,𝐹,𝐺and 𝑆are related to the original vectors U, F, G and S as these expressions, using the direct metrics: 𝑈=𝐽⋅𝑈 Equation 2.66 𝐹=𝐽⋅𝐹𝜕𝜉 𝜕𝑥+𝐽⋅𝐺𝜕𝜉 𝜕𝑦 Equation 2.67 𝐺=𝐽⋅𝐹𝜕𝜂 𝜕𝑥+𝐽⋅𝐺𝜕𝜂 𝜕𝑦 Equation 2.68 𝑆=𝐽⋅𝑆 Equation 2.69 Alternatively, the vectors 𝐹 and 𝐺 can be calculated with the inverse metrics: 32 𝐹=𝐹𝜕𝑦 𝜕𝜂−𝐺𝜕𝑥 𝜕𝜂 Equation 2.70 𝐺=−𝐹𝜕𝑦 𝜕𝜉+𝐺𝜕𝑥 𝜕𝜉 Equation 2.71 The possibility of using both sets of expressions has been very important in the success of the project, for reasons that will be explained in the next chapter. Exactly as it happened in subchapter 2.2, the variables that form part of the new vectors 𝐹 and 𝐺 must be differentiated in the opposite direction than the vector the form part of. As now the original vectors F and G both are into the calculation of the new ones, keeping track of what variable must be differentiated in what direction has been in some situations very messy. However, it is not difficult to follow if the concepts explained previously are clear and the procedure is done orderly, step by step. If the viscosity terms stated in equations Equation 2.12, Equation 2.13 and Equation 2.14 and the temperature transfer terms   and   are transformed by the derivative relations Equation 2.50 and Equation 2.51, the results are the following: 𝜏=23𝜇2𝜕𝑢 𝜕𝜉𝜕𝜉 𝜕𝑥+𝜕𝑢 𝜕𝜂𝜕𝜂 𝜕𝑥−𝜕𝑣 𝜕𝜉𝜕𝜉 𝜕𝑦+𝜕𝑣 𝜕𝜂𝜕𝜂 𝜕𝑦 Equation 2.72 𝜏=23𝜇2𝜕𝑣 𝜕𝜉𝜕𝜉 𝜕𝑦+𝜕𝑣 𝜕𝜂𝜕𝜂 𝜕𝑦−𝜕𝑢 𝜕𝜉𝜕𝜉 𝜕𝑥+𝜕𝑢 𝜕𝜂𝜕𝜂 𝜕𝑥 Equation 2.73 𝜏=𝜏=𝜇𝜕𝑣 𝜕𝜉𝜕𝜉 𝜕𝑥+𝜕𝑣 𝜕𝜂𝜕𝜂 𝜕𝑥+𝜕𝑢 𝜕𝜉𝜕𝜉 𝜕𝑦+𝜕𝑢 𝜕𝜂𝜕𝜂 𝜕𝑦 Equation 2.74 𝜕𝑇 𝜕𝑥=𝜕𝑇 𝜕𝜉𝜕𝜉 𝜕𝑥+𝜕𝑇 𝜕𝜂𝜕𝜂 𝜕𝑥 Equation 2.75 𝜕𝑇 𝜕𝑦=𝜕𝑇 𝜕𝜉𝜕𝜉 𝜕𝑦+𝜕𝑇 𝜕𝜂𝜕𝜂 𝜕𝑦 Equation 2.76 The terms  ,  ,  ,  ,   and   now must be calculated with forward, rearward, or central differentiation depending on the vector they are going to be used into. What happens is the following. The terms of equations Equation 2.72 to Equation 2.76 form part of the vectors F and G, and at the same time they form part of 𝐹 and 𝐺, which are differentiated differently. 33 With this reasoning, it becomes clear that two variants of the vectors F and G must be calculated, which will be called 𝐹, 𝐹, 𝐺 and 𝐺. Each variant is designated to be used in the formation of one of the vectors 𝐹 or 𝐺, and is formed with variables differentiated accordingly to where they are going to be used. The best way to visualize this concept is by looking at an example, in this case the viscosity term 𝜏. These equations show explicitly the differentiation used for this term when calculating 𝐹 and 𝐺 in the predictor step. 𝜏,󰇻= 23𝜇, 󰇯2󰇭𝜕𝜉 𝜕𝑥󰇧𝑢, −𝑢,  𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑥󰇧𝑢, −𝑢,  2 𝛥𝜂 󰇨󰇮 −󰇭𝜕𝜉 𝜕𝑦󰇧𝑣, −𝑣,  𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑦󰇧𝑣, −𝑣,  2 𝛥𝜂 󰇨󰇮󰇰 Equation 2.77 𝜏,󰇻= 23𝜇, 󰇯2󰇭𝜕𝜉 𝜕𝑥󰇧𝑢, −𝑢,  2 𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑥󰇧𝑢, −𝑢,  𝛥𝜂 󰇨󰇮 −󰇭𝜕𝜉 𝜕𝑦󰇧𝑣, −𝑣,  2 𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑦󰇧𝑣, −𝑣,  𝛥𝜂 󰇨󰇮󰇰 Equation 2.78 As said, the expressions can get confusing, but the compliance with the differentiation requisite is important for the second order accuracy of the method. In order to alleviate the section of extremely large expressions, not all the equivalents to Equation 2.77 and Equation 2.78 with the other variables are going to be shown. Just for comparison, appreciate how the differentiation used changes in the corrector step. 𝜏,󰇻= 23𝜇, 󰇯2󰇭𝜕𝜉 𝜕𝑥󰇧𝑢, −𝑢,  𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑥󰇧𝑢, −𝑢,  2 𝛥𝜂 󰇨󰇮 −󰇭𝜕𝜉 𝜕𝑦󰇧𝑣, −𝑣,  𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑦󰇧𝑣, −𝑣,  2 𝛥𝜂 󰇨󰇮󰇰 Equation 2.79 𝜏,󰇻= 23𝜇, 󰇯2󰇭𝜕𝜉 𝜕𝑥󰇧𝑢, −𝑢,  2 𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑥󰇧𝑢, −𝑢,  𝛥𝜂 󰇨󰇮 −󰇭𝜕𝜉 𝜕𝑦󰇧𝑣, −𝑣,  2 𝛥𝜉 󰇨+𝜕𝜂 𝜕𝑦󰇧𝑣, −𝑣,  𝛥𝜂 󰇨󰇮󰇰 Equation 2.80 34 The correct implementation of these expressions has been one of the most difficult parts of the project, and its correctness is of vital importance. Despite not all the expressions appearing in this theoretical chapter, the reader can check their implementation by taking a look at the sections “Encode variables for predictor step” and “Encode variables for corrector step“ of the script that appears at the end of subchapter 3.2. 35 3 Numerical resolution To codify the technique and perform the computation required to obtain the results, a programming language must be chosen to develop the code and then run it. In this project the program of choice has been MATLAB®, due to several reasons. Among them, that it is a program in which a considerable amount of focus is put into during the Aerospace Degree, and so it didn’t require learning a new language. Also, the University of Seville has educational licenses and it can be easily accessed. There have been a few scripts written to achieve the numerical resolution. However, they are structured as is going to be explained. The base script is the one called “MainBasic.m”. It represents the bedrock of the coding and is the one who calls the other scripts, that are defined as functions. Those functions are the responsible of setting the starting conditions, the boundary conditions, and assess the results, for instance to check the validity of the results and represent the plots. “MainBasic.m”, as well as almost all the other scripts, have had multiple versions, denoted by the adding of, for example, “_v1” or “_v3” to the tittle. The script is structured in sections, and the following subchapters are going to explain the high-level functioning of the script section by section. The different subscripts, called “auxiliary functions”, vary depending on the problem that is going to be solved, thus will be treated in the next chapter for each application. Finally, the chapter will end by explaining the modifications of the code that are needed to implement the non-uniform grid transformations by the method described in the subchapter 2.3. This chapter is mostly written to provide a comprehensive guide about how to understand and read the code if any person wants to take this project as base for a future project. 3.1 The Program for Uniform Grids Before starting the MacCormack’s technique to advance the simulation in time, the whole problem and a lot of variables must be set. The purpose of the sections before the MacCormack’s technique begins is just that. After deleting all current plots any pre-existing variables in the work space, the script “MainBasic.m” begins with the first section called “Properties of the air”, in which the basic parameters of the air are set, such as the Prandtl number, the gas constant R, the ratio of specific heats 𝛾, and the reference values at sea level 𝑇 and 𝜇 needed for the Sutherland’s law. Then 𝑐 and 𝑐 are calculated from the values of R and 𝛾. The second section is called “Entry air data”, which specifies the parameters of the incoming air like Mach number and temperature. Following it begins the section “Geometry and grid definition”. It determines the number of points in the grid, the security factor for the time step and also the dimensions of the computational domain. 36 The most important of the set-up sections is the “Initial Conditions” section. It calls the auxiliary function “Initial_Conditions_Application.m”, that must return the value of the density, vertical and horizontal velocity, temperature and pressure at all grid points at the simulation starting time. The values of the internal energy 𝑒, viscosity 𝜇 and the heat transfer coefficient 𝑘 are calculated from this data. Due to the evolution of the program, the values of the predicted and the corrected variables, which are represented with the addition of “_p” and “_c” to the variables respectively, have to be set in the first iteration. So, they are simply made equal to the initial conditions. After the initial conditions are finished, the preparations end by allocating variables to speed up the code and specifying the maximum time and iterations for the simulation. The advance in time is achieved by a loop “while”, that checks in each iteration if the number of iterations or the time are greater than the maximums specified, if convergence has been achieved, or if a custom condition expressed by a “flag” variable has occurred. Before the loop, the time, iterations, convergence and flag are set to zero. All the stopping conditions are optional, and they can be disabled at will. Inside the loop, the first section is dedicated to the calculation of the time step. A mathematically proven criterion has not been demonstrated jet, however, following the recommendation made by John Anderson in chapter 10 of [Ref 1], the Courant-Friedrichs-Lewy (CFL) criterion has been used to determine the time step which has been empirically proven enough to preserve a convergent behaviour. (∆𝑡),=󰇯𝑢, ∆𝑥 +𝑣, ∆𝑦 +𝑎,1 ∆𝑥+1 ∆𝑦+2𝑣, 󰆒1 ∆𝑥+1 ∆𝑦󰇰 Equation 3.1. Courant-Friedrichs-Lewy criterion. Where 𝑣, 󰆒=𝑚𝑎𝑥󰇯43𝜇,(𝛾𝜇,/𝑃𝑟) 𝜌, 󰇰 Equation 3.2 The time step is the minimum of the values of ∆𝑡, with a security factor that ensures the numerical stability that was set in the set-up phase. ∆𝑡=𝐾⋅𝑚𝑖𝑛(∆𝑡), Equation 3.3 As to apply the boundary conditions, the variables have to be decoded after each step. The implementation of the predictor and corrector steps has multiple sections, and the structure of these sections is the following: 37 1. Calculation of the variables for the predictor step with the equations Equation 2.27 to Equation 2.32. 2. Formation of the vectors U, F, G and S as stated in the expressions Equation 2.6 to Equation 2.9. 3. Calculation of the predicted vector U_p with the Equation 2.26. The predicted variables have as common distinction the adding of “_p” to their names. 4. Decoding of the predicted primitive variables with the procedure explained at the beginning of section 2.2. 5. Enforcement of the boundary conditions on the predicted primitive variables. This enforcement is not done directly in the “MainBasic.m” script. It is done by calling an auxiliary function called “Boundary_Conditions_Application.m”, which changes with each application. This, together with the “Initial_Conditions_Application.m”, allows great flexibility with the coding, as the main program is general and for each application only the specific initial and boundary condition scripts have to be written. 6. Calculation of the variables for the corrector step with equations equivalent to Equation 2.27 to Equation 2.32. Now, also new predicted values of viscosity and heat transfer coefficient must be calculated to use them in the formation of the “conservative form” vectors. 7. Formation of the predicted vectors U, F, G and S with the expressions Equation 2.6 to Equation 2.9, only that now the predicted values of the variables are used. 8. Calculation of the corrected vector U_c with the Equation 2.33. The corrected variables have as common distinction the adding of “_c” to their names. 9. Decoding of the corrected primitive variables with the procedure explained at the beginning of section 2.2. 10. Enforcement of the boundary conditions on the corrected primitive variables with the same “Boundary_Conditions_Application.m” function used in the predictor step. 11. Reformation of the vector U_c with the boundary conditions already applied. 12. Combination of the steps to achieve the final results with Equation 2.34. The last part encompasses several sections with different purposes. The first section checks the convergence of the code by comparing the new obtained density at time 𝑡+∆𝑡 with the density at time t from the previous iteration in each computational point. If the maximum difference is less than 10 the convergence variable is activated and, if the option is enabled, the program exits the loop at the end of the current iteration. As a comment, this condition 38 could potentially be improved by substituting the static number with an expression that takes into consideration the time step. The next section decodes the obtained vector 𝑈∆, erasing the old variables from the previous iteration in the process and after that the time and iterations are updated. In the end, the current time and time step are displayed to allow a visual check of the correct execution of the program. This slightly slows the program, but not much, and is very helpful when dealing with a new script that may or may not work at the first time. Optionally, it can be enabled the plotting of a figure to see the evolution of the variables. However, this greatly slows the execution, and should only be used in exceptional situations. Lastly, the variable “flag” can be used to stop the program if a custom condition has been met and the option is enabled. In the code for this project, it checks if the viscosity has any complex component, as this is signal that the simulation has somehow broken. When the program ends the loop because one of the stopping conditions has been met, a check of the validity of the results has been implemented by integrating the mass flow on the boundaries of the computational domain. If the difference between the mass that enters the domain and the mass than exits it is less than one percent, the solution is considered valid. If the value is more than that, a warning is displayed and the amount of mass differential in percentage with respect to the mass inflow is also displayed. 39 Main program for rectangular and uniform grids “MainBasic.m” clear all; close all; clc; tic %%%%%%%%%%% MAIN PROGRAM OF THE MACCORMACK'S METHOD %%%%%%%%%%% % FOR UNIFORM RECTANGULAR GRIDS %---------------------------------------------------------- %% Air properties % Prandlt number Pr = 0.71; % Specific heats ratio gam = 1.4; % Ideal gas constant R = 287; % Heat coefficients c_v = R/(gam-1); c_p = gam*c_v; % Reference values T0 = 298; nu0 = 1.849e-5; %---------------------------------------------------------- %% Air flow stream data Minf = 4; % Choose the desired Mach number Tinf = 298; rhoinf = 1.225; pinf = rhoinf*R*Tinf; nuinf = nu0*(Tinf/T0).^(3/2).*((T0+110)./(Tinf+110)); uinf = Minf*sqrt(gam*R*Tinf); %---------------------------------------------------------- %% Geometry and grid fineness Nx = 70; Ny = 70; % Choose the adequate grid fineness m = Nx-1; n = Ny-1; K = 0.3; % Choose the adequate safe factor LHOR = 0.00001; % Modify to fit the desire dimensions Reinf = rhoinf*uinf*LHOR/nuinf; delta = 5*LHOR/sqrt(Reinf); LVER = 5*delta; % or LVER = fix value; % Choose at will dx = LHOR /(Nx-1); dy = LVER /(Ny-1); %---------------------------------------------------------- %% Initial conditions % Real variables % Choose the correct auxiliary function depending on the application [rho,u,v,T,p] = Initial_Conditions(rhoinf,uinf,Tinf,R,Nx,Ny); e = c_v*T; nu = nu0*T.*sqrt(T)/T0^(3/2).*((T0+110)./(T+110)); 40 k = nu*c_p/Pr; % Predictor variables rho_p = rho; u_p = u; v_p = v; T_p = T; p_p = p; e_p = e; nu_p = nu; k_p = k; % Corrector variables rho_c = rho; u_c = u; v_c = v; T_c = T; p_c = p; e_c = e; nu_c = nu; k_c = k; %---------------------------------------------------------- %% Maximum time and iterations tmax = 4e-7; itermax = 5000; %---------------------------------------------------------- %% Allocate Variables U = zeros(Nx,Ny,4); F = zeros(Nx,Ny,4); G = zeros(Nx,Ny,4); S = zeros(Nx,Ny,4); U_p = zeros(size(U)); U_c = zeros(size(U)); dudx = zeros(size(u)); dudy = zeros(size(u)); dvdx = zeros(size(u)); dvdy = zeros(size(u)); dTdx = zeros(size(u)); dTdy = zeros(size(u)); %---------------------------------------------------------- %% MacCormack's Method advance t = 0; iter = 0; conver = 0; flag = 0; % Enable and disable at will the conditions to stop the loop while t<tmax && flag==0 && conver==0 && iter<itermax %------------------------------------------------------ %% Time step calculation M = u./sqrt(gam*R*T); 47 3.2 Non-Uniform Grid Modifications This subchapter will explain the modifications that has been made to the code in order to implement the non-uniform grid transformations as they were explained in the subchapter 2.3. The new main program is now called “MainGridTrans_DirectMetrics.m”. The first difference is that the set-up section “Geometry and grid definition” is now much more complex. It makes use of an auxiliary, application specific function, which must return the coordinates x and y of the grid points, the Jacobian at each point, and the direct and inverse metrics. In principle, with the relations given by Equation 2.57, only one set of metrics could be used. However, as it will be seen later, the inverse metrics are necessary at one point of the coding. Using the x and y coordinates, the distance between points ∆𝑥 and ∆𝑦 must be calculated too to use them in the Courant-Friedrichs-Lewy criterion. As it has been said before, the fact that both vectors F and G form part of the new 𝐹 and 𝐺 forces the formation of two variants of F and G, to comply with the opposite differentiation needed for the second order accuracy. Thus, the vectors 𝐹, 𝐹, 𝐺 and 𝐺 are formed, each using the adequate variables. Then, the vectors 𝐹 and 𝐺 are used in the Equation 2.70, together with the inverse metrics, to form 𝐹. The same happens with 𝐹 and 𝐺 that are used to form 𝐺 with Equation 2.71. An important commentary has to be made here. In the formation of 𝑈 the Jacobian “J” is used, and it is calculated using the inverse metrics. When exact analytical expressions are used, the mixing of both direct and inverse metrics seems to be irrelevant. However, if the metrics are approximately calculated with central differentiation, as it is the case in most of the applications in this project, it exists a small error between the two sets. This very small error adds-up iteration after iteration and ends up making the results inexact. To avoid that, the expressions used to calculate 𝐹 and 𝐺 have been Equation 2.70 and Equation 2.71, which use inverse metrics, instead of Equation 2.68 and Equation 2.69, which because of using direct metrics produced this very subtle error. 48 Main program for non-uniform grids “MainGridTrans_DirectMetrics.m” clear all; close all; clc; tic %%%%%%%%%%% MAIN PROGRAM OF THE MACCORMACK'S METHOD %%%%%%%%%%% % FOR NON-UNIFORM GRID TRANSFORMATIONS %---------------------------------------------------------- %% Air properties % Prandlt number Pr = 0.71; % Specific heats ratio gam = 1.4; % Ideal gas constant R = 287; % Heat coefficients c_v = R/(gam-1); c_p = gam*c_v; % Reference values T0 = 298; nu0 = 1.849e-5; %---------------------------------------------------------- %% Air flow stream data Minf = 4; % Choose the desired Mach number Tinf = 298; rhoinf = 1.225; pinf = rhoinf*R*Tinf; nuinf = nu0*(Tinf/T0).^(3/2).*((T0+110)./(Tinf+110)); uinf = Minf*sqrt(gam*R*Tinf); %---------------------------------------------------------- %% Geometry and grid fineness Nxi = 100; Neta = 150; % Choose the adequate grid fineness m = Nxi-1; n = Neta-1; K = 0.1; % Choose the adequate safe factor LHOR = 0.00001; % Modify to fit the desire dimensions Reinf = rhoinf*uinf*LHOR/nuinf; delta = 5*LHOR/sqrt(Reinf); LVER = 1.6345e-05; % or LVER = value*delta; % Choose at will [x,y,dxi,deta,dxidx,dxidy,detadx,detady,J,dxdxi,dxdeta,dydxi,dydeta] = Geometry_Function(Nxi,Neta,LHOR,LVER); delta_x(1,:) = x(2,:) - x(1,:); delta_x(2:m,:) = ( x(3:m+1,:) - x(1:m-1,:) )/2; delta_x(Nxi,:) = x(Nxi,:) - x(Nxi-1,:); delta_y(:,1) = y(:,2) - y(:,1); delta_y(:,2:n) = ( y(:,3:n+1) - y(:,1:n-1) )/2; delta_y(:,Neta) = y(:,Neta) - y(:,Neta-1); 49 %---------------------------------------------------------- %% Initial conditions % Real variables % Choose the correct auxiliary function depending on the application [rho,u,v,T,p] = Initial_Conditions(rhoinf,uinf,Tinf,R,Nxi,Neta); e = c_v*T; nu = nu0*T.*sqrt(T)/T0^(3/2).*((T0+110)./(T+110)); k = nu*c_p/Pr; % Predictor variables rho_p = rho; u_p = u; v_p = v; T_p = T; p_p = p; e_p = e; nu_p = nu; k_p = k; % Corrector variables rho_c = rho; u_c = u; v_c = v; T_c = T; p_c = p; e_c = e; nu_c = nu; k_c = k; %---------------------------------------------------------- %% Maximum time and iterations tmax = 2.5e-07; %1.5e-7; itermax = 5000; %---------------------------------------------------------- %% Allocate Variables U_p = zeros(Nxi,Neta,4); U_c = zeros(Nxi,Neta,4); S = zeros(Nxi,Neta,4); %---------------------------------------------------------- %% MacCormack's Method advance t = 0; iter = 0; conver = 0; flag = 0; % Enable and disable at will the conditions to stop the loop while t<tmax && flag==0 && conver==0 && iter<itermax %------------------------------------------------------- %% Time step calculation M = u./sqrt(gam*R*T); 50 v_prim = max( max( (4*gam*nu.*nu)./(3.*rho*Pr) ) ); dt = 1./( abs(u)./delta_x + abs(v)./delta_y + sqrt(gam*R*T).*sqrt(1./delta_x.^2+1./delta_y.^2) + 2*v_prim.*(1./delta_x.^2+1./delta_y.^2) ); dt = K*min( min (dt) ); %------------------------------------------------------- %% Encode variables for predictor step V = sqrt(u.*u + v.*v); % DR means Rearward Difference % DC means Central Difference % DP means Forward Difference % xi:DR dudxi_R(1,:) = ( u(2,:) - u(1,:) )/dxi; dudxi_R(2:Nxi,:) = ( u(2:Nxi,:) - u(1:Nxi-1,:) )/dxi; dvdxi_R(1,:) = ( v(2,:) - v(1,:) )/dxi; dvdxi_R(2:Nxi,:) = ( v(2:Nxi,:) - v(1:Nxi-1,:) )/dxi; dTdxi_R(1,:) = ( T(2,:) - T(1,:) )/dxi; dTdxi_R(2:Nxi,:) = ( T(2:Nxi,:) - T(1:Nxi-1,:) )/dxi; % xi:DC dudxi_C(1,:) = ( u(2,:) - u(1,:) )/dxi; dudxi_C(2:m,:) = ( u(3:m+1,:) - u(1:m-1,:) )/(2*dxi); dudxi_C(Nxi,:) = ( u(Nxi,:) - u(Nxi-1,:) )/dxi; dvdxi_C(1,:) = ( v(2,:) - v(1,:) )/dxi; dvdxi_C(2:m,:) = ( v(3:m+1,:) - v(1:m-1,:) )/(2*dxi); dvdxi_C(Nxi,:) = ( v(Nxi,:) - v(Nxi-1,:) )/dxi; dTdxi_C(1,:) = ( T(2,:) - T(1,:) )/dxi; dTdxi_C(2:m,:) = ( T(3:m+1,:) - T(1:m-1,:) )/(2*dxi); dTdxi_C(Nxi,:) = ( T(Nxi,:) - T(Nxi-1,:) )/dxi; % eta:DR dudeta_R(:,1) = ( u(:,2) - u(:,1) )/deta; dudeta_R(:,2:Neta) = ( u(:,2:Neta) - u(:,1:Neta-1) )/deta; dvdeta_R(:,1) = ( v(:,2) - v(:,1) )/deta; dvdeta_R(:,2:Neta) = ( v(:,2:Neta) - v(:,1:Neta-1) )/deta; dTdeta_R(:,1) = ( T(:,2) - T(:,1) )/deta; dTdeta_R(:,2:Neta) = ( T(:,2:Neta) - T(:,1:Neta-1) )/deta; % eta:DC dudeta_C(:,1) = ( u(:,2) - u(:,1) )/deta; dudeta_C(:,2:n) = ( u(:,3:n+1) - u(:,1:n-1) )/(2*deta); dudeta_C(:,Neta) = ( u(:,Neta) - u(:,Neta-1) )/deta; dvdeta_C(:,1) = ( v(:,2) - v(:,1) )/deta; dvdeta_C(:,2:n) = ( v(:,3:n+1) - v(:,1:n-1) )/(2*deta); dvdeta_C(:,Neta) = ( v(:,Neta) - v(:,Neta-1) )/deta; dTdeta_C(:,1) = ( T(:,2) - T(:,1) )/deta; 51 dTdeta_C(:,2:n) = ( T(:,3:n+1) - T(:,1:n-1) )/(2*deta); dTdeta_C(:,Neta) = ( T(:,Neta) - T(:,Neta-1) )/deta; % tauxx tauxx_F1 = 2/3*nu.*( 2*( dxidx.*dudxi_R + detadx.*dudeta_C ) - ( dxidy.*dvdxi_R + detady.*dvdeta_C ) ); tauxx_G1 = 2/3*nu.*( 2*( dxidx.*dudxi_C + detadx.*dudeta_R ) - ( dxidy.*dvdxi_C + detady.*dvdeta_R ) ); % tauyy tauyy_F1 = 2/3*nu.*( 2*( dxidy.*dvdxi_R + detady.*dvdeta_C ) - ( dxidx.*dudxi_R + detadx.*dudeta_C ) ); tauyy_G1 = 2/3*nu.*( 2*( dxidy.*dvdxi_C + detady.*dvdeta_R ) - ( dxidx.*dudxi_C + detadx.*dudeta_R ) ); % tauxy tauxy_F1 = nu.*( ( dxidy.*dudxi_R + detady.*dudeta_C ) + ( dxidx.*dvdxi_R + detadx.*dvdeta_C ) ); tauxy_G1 = nu.*( ( dxidy.*dudxi_C + detady.*dudeta_R ) + ( dxidx.*dvdxi_C + detadx.*dvdeta_R ) ); % dTdx dTdx_F1 = dxidx.*dTdxi_R + detadx.*dTdeta_C; dTdx_G1 = dxidx.*dTdxi_C + detadx.*dTdeta_R; % dTdy dTdy_F1 = dxidy.*dTdxi_R + detady.*dTdeta_C; dTdy_G1 = dxidy.*dTdxi_C + detady.*dTdeta_R; %------------------------------------------------------- %% Predictor step: xi:DP eta:DP U(:,:,1) = rho; U(:,:,2) = rho.*u; U(:,:,3) = rho.*v; U(:,:,4) = rho.*(e+V.*V/2); F_F1(:,:,1) = rho.*u; F_F1(:,:,2) = rho.*u.*u + p - tauxx_F1; F_F1(:,:,3) = rho.*v.*u - tauxy_F1; F_F1(:,:,4) = rho.*(e+V.*V/2).*u + p.*u - k.*dTdx_F1 - u.*tauxx_F1 - v.*tauxy_F1; G_F1(:,:,1) = rho.*v; G_F1(:,:,2) = rho.*u.*v - tauxy_F1; G_F1(:,:,3) = rho.*v.*v + p - tauyy_F1; G_F1(:,:,4) = rho.*(e+V.*V/2).*v + p.*v - k.*dTdy_F1 - u.*tauxy_F1 - v.*tauyy_F1; F_G1(:,:,1) = rho.*u; F_G1(:,:,2) = rho.*u.*u + p - tauxx_G1; F_G1(:,:,3) = rho.*v.*u - tauxy_G1; F_G1(:,:,4) = rho.*(e+V.*V/2).*u + p.*u - k.*dTdx_G1 - u.*tauxx_G1 - v.*tauxy_G1; G_G1(:,:,1) = rho.*v; G_G1(:,:,2) = rho.*u.*v - tauxy_G1; G_G1(:,:,3) = rho.*v.*v + p - tauyy_G1; 52 G_G1(:,:,4) = rho.*(e+V.*V/2).*v + p.*v - k.*dTdy_G1 - u.*tauxy_G1 - v.*tauyy_G1; for i=1:4 U1(:,:,i) = J.*U(:,:,i); F1(:,:,i) = F_F1(:,:,i).*dydeta - G_F1(:,:,i).*dxdeta; G1(:,:,i) = - F_G1(:,:,i).*dydxi + G_G1(:,:,i).*dxdxi; end U1_p(2:m,2:n,:) = U1(2:m,2:n,:) - dt*(F1(3:m+1,2:n,:)- F1(2:m,2:n,:))/dxi - dt*(G1(2:m,3:n+1,:)-G1(2:m,2:n,:))/deta + dt*S(2:m,2:n,:); %------------------------------------------------------- %% Decode variables of predictor step for i=1:4 U_p(2:m,2:n,i) = U1_p(2:m,2:n,i)./J(2:m,2:n); end rho_p(2:m,2:n) = U_p(2:m,2:n,1); u_p(2:m,2:n) = U_p(2:m,2:n,2)./U_p(2:m,2:n,1); v_p(2:m,2:n) = U_p(2:m,2:n,3)./U_p(2:m,2:n,1); e_p(2:m,2:n) = U_p(2:m,2:n,4)./U_p(2:m,2:n,1) - (u_p(2:m,2:n).*u_p(2:m,2:n) + v_p(2:m,2:n).*v_p(2:m,2:n))/2; T_p(2:m,2:n) = e_p(2:m,2:n)/c_v; p_p(2:m,2:n) = R*rho_p(2:m,2:n).*T_p(2:m,2:n); %------------------------------------------------------- %% Boundary conditions of predictor step % Choose the correct auxiliary function depending on the application [rho_p,u_p,v_p,T_p,p_p] = Boundary_Conditions(rho_p,u_p,v_p,T_p,p_p,Nxi,Neta); e_p(Nxi,2:Neta-1) = c_v*T_p(Nxi,2:Neta-1); %------------------------------------------------------- %% Encode variables for corrector step V_p = sqrt(u_p.*u_p + v_p.*v_p); nu_p = nu0*(T_p/T0).^(3/2).*((T0+110)./(T_p+110)); k_p = nu_p*c_p/Pr; % xi:DP dudxi_P_p(1:m,:) = ( u_p(2:m+1,:) - u_p(1:m,:) )/dxi; dudxi_P_p(Nxi,:) = ( u_p(Nxi,:) - u_p(Nxi-1,:) )/dxi; dvdxi_P_p(1:m,:) = ( v_p(2:m+1,:) - v_p(1:m,:) )/dxi; dvdxi_P_p(Nxi,:) = ( v_p(Nxi,:) - v_p(Nxi-1,:) )/dxi; dTdxi_P_p(1:m,:) = ( T_p(2:m+1,:) - T_p(1:m,:) )/dxi; dTdxi_P_p(Nxi,:) = ( T_p(Nxi,:) - T_p(Nxi-1,:) )/dxi; % xi:DC dudxi_C_p(1,:) = ( u_p(2,:) - u_p(1,:) )/dxi; dudxi_C_p(2:m,:) = ( u_p(3:m+1,:) - u_p(1:m-1,:) )/(2*dxi); dudxi_C_p(Nxi,:) = ( u_p(Nxi,:) - u_p(Nxi-1,:) )/dxi; 53 dvdxi_C_p(1,:) = ( v_p(2,:) - v_p(1,:) )/dxi; dvdxi_C_p(2:m,:) = ( v_p(3:m+1,:) - v_p(1:m-1,:) )/(2*dxi); dvdxi_C_p(Nxi,:) = ( v_p(Nxi,:) - v_p(Nxi-1,:) )/dxi; dTdxi_C_p(1,:) = ( T_p(2,:) - T_p(1,:) )/dxi; dTdxi_C_p(2:m,:) = ( T_p(3:m+1,:) - T_p(1:m-1,:) )/(2*dxi); dTdxi_C_p(Nxi,:) = ( T_p(Nxi,:) - T_p(Nxi-1,:) )/dxi; % eta:DP dudeta_P_p(:,1:n) = ( u_p(:,2:n+1) - u_p(:,1:n) )/deta; dudeta_P_p(:,Neta) = ( u_p(:,Neta) - u_p(:,Neta-1) )/deta; dvdeta_P_p(:,1:n) = ( v_p(:,2:n+1) - v_p(:,1:n) )/deta; dvdeta_P_p(:,Neta) = ( v_p(:,Neta) - v_p(:,Neta-1) )/deta; dTdeta_P_p(:,1:n) = ( T_p(:,2:n+1) - T_p(:,1:n) )/deta; dTdeta_P_p(:,Neta) = ( T_p(:,Neta) - T_p(:,Neta-1) )/deta; % eta:DC dudeta_C_p(:,1) = ( u_p(:,2) - u_p(:,1) )/deta; dudeta_C_p(:,2:n) = ( u_p(:,3:n+1) - u_p(:,1:n-1) )/(2*deta); dudeta_C_p(:,Neta) = ( u_p(:,Neta) - u_p(:,Neta-1) )/deta; dvdeta_C_p(:,1) = ( v_p(:,2) - v_p(:,1) )/deta; dvdeta_C_p(:,2:n) = ( v_p(:,3:n+1) - v_p(:,1:n-1) )/(2*deta); dvdeta_C_p(:,Neta) = ( v_p(:,Neta) - v_p(:,Neta-1) )/deta; dTdeta_C_p(:,1) = ( T_p(:,2) - T_p(:,1) )/deta; dTdeta_C_p(:,2:n) = ( T_p(:,3:n+1) - T_p(:,1:n-1) )/(2*deta); dTdeta_C_p(:,Neta) = ( T_p(:,Neta) - T_p(:,Neta-1) )/deta; % tauxx tauxx_F1_p = 2/3*nu.*( 2*( dxidx.*dudxi_P_p + detadx.*dudeta_C_p ) - ( dxidy.*dvdxi_P_p + detady.*dvdeta_C_p ) ); tauxx_G1_p = 2/3*nu.*( 2*( dxidx.*dudxi_C_p + detadx.*dudeta_P_p ) - ( dxidy.*dvdxi_C_p + detady.*dvdeta_P_p ) ); % tauyy tauyy_F1_p = 2/3*nu.*( 2*( dxidy.*dvdxi_P_p + detady.*dvdeta_C_p ) - ( dxidx.*dudxi_P_p + detadx.*dudeta_C_p ) ); tauyy_G1_p = 2/3*nu.*( 2*( dxidy.*dvdxi_C_p + detady.*dvdeta_P_p ) - ( dxidx.*dudxi_C_p + detadx.*dudeta_P_p ) ); % tauxy tauxy_F1_p = nu.*( ( dxidy.*dudxi_P_p + detady.*dudeta_C_p ) + ( dxidx.*dvdxi_P_p + detadx.*dvdeta_C_p ) ); tauxy_G1_p = nu.*( ( dxidy.*dudxi_C_p + detady.*dudeta_P_p ) + ( dxidx.*dvdxi_C_p + detadx.*dvdeta_P_p ) ); % dTdx dTdx_F1_p = dxidx.*dTdxi_P_p + detadx.*dTdeta_C_p; dTdx_G1_p = dxidx.*dTdxi_C_p + detadx.*dTdeta_P_p; % dTdy dTdy_F1_p = dxidy.*dTdxi_P_p + detady.*dTdeta_C_p; dTdy_G1_p = dxidy.*dTdxi_C_p + detady.*dTdeta_P_p; %------------------------------------------------------- 54 %% Corrector step: xi:DR eta:DR U_p(:,:,1) = rho_p; U_p(:,:,2) = rho_p.*u_p; U_p(:,:,3) = rho_p.*v_p; U_p(:,:,4) = rho_p.*(e_p+V_p.*V_p/2); F_F1_p(:,:,1) = rho_p.*u_p; F_F1_p(:,:,2) = rho_p.*u_p.*u_p + p_p - tauxx_F1_p; F_F1_p(:,:,3) = rho_p.*v_p.*u_p - tauxy_F1_p; F_F1_p(:,:,4) = rho_p.*(e_p+V_p.*V_p/2).*u_p + p_p.*u_p - k_p.*dTdx_F1_p - u_p.*tauxx_F1_p - v_p.*tauxy_F1_p; G_F1_p(:,:,1) = rho_p.*v_p; G_F1_p(:,:,2) = rho_p.*u_p.*v_p - tauxy_F1_p; G_F1_p(:,:,3) = rho_p.*v_p.*v_p + p_p - tauyy_F1_p; G_F1_p(:,:,4) = rho_p.*(e_p+V_p.*V_p/2).*v_p + p_p.*v_p - k_p.*dTdy_F1_p - u_p.*tauxy_F1_p - v_p.*tauyy_F1_p; F_G1_p(:,:,1) = rho_p.*u_p; F_G1_p(:,:,2) = rho_p.*u_p.*u_p + p_p - tauxx_G1_p; F_G1_p(:,:,3) = rho_p.*v_p.*u_p - tauxy_G1_p; F_G1_p(:,:,4) = rho_p.*(e_p+V_p.*V_p/2).*u_p + p_p.*u_p - k_p.*dTdx_G1_p - u_p.*tauxx_G1_p - v_p.*tauxy_G1_p; G_G1_p(:,:,1) = rho_p.*v_p; G_G1_p(:,:,2) = rho_p.*u_p.*v_p - tauxy_G1_p; G_G1_p(:,:,3) = rho_p.*v_p.*v_p + p_p - tauyy_G1_p; G_G1_p(:,:,4) = rho_p.*(e_p+V_p.*V_p/2).*v_p + p_p.*v_p - k_p.*dTdy_G1_p - u_p.*tauxy_G1_p - v_p.*tauyy_G1_p; for i=1:4 F1_p(:,:,i) = F_F1_p(:,:,i).*dydeta - G_F1_p(:,:,i).*dxdeta; G1_p(:,:,i) = - F_G1_p(:,:,i).*dydxi + G_G1_p(:,:,i).*dxdxi; end U1_c(2:m,2:n,:) = U1(2:m,2:n,:) - dt*(F1_p(2:m,2:n,:)-F1_p(1:m- 1,2:n,:))/dxi - dt*(G1_p(2:m,2:n,:)-G1_p(2:m,1:n-1,:))/deta + dt*S(2:m,2:n,:); %------------------------------------------------------- %% Decode variables of corrector step for i=1:4 U_c(2:m,2:n,i) = U1_c(2:m,2:n,i)./J(2:m,2:n); end rho_c(2:m,2:n) = U_c(2:m,2:n,1); u_c(2:m,2:n) = U_c(2:m,2:n,2)./U_c(2:m,2:n,1); v_c(2:m,2:n) = U_c(2:m,2:n,3)./U_c(2:m,2:n,1); e_c(2:m,2:n) = U_c(2:m,2:n,4)./U_c(2:m,2:n,1) - (u_c(2:m,2:n).*u_c(2:m,2:n) + v_c(2:m,2:n).*v_c(2:m,2:n))/2; T_c(2:m,2:n) = e_c(2:m,2:n)/c_v; % p_c(2:m,2:n) = R*rho_c(2:m,2:n).*T_c(2:m,2:n); % not necessary decoding %------------------------------------------------------- %% Boundary conditions of corrector step 55 [rho_c,u_c,v_c,T_c,p_c] = Condiciones_Contorno_Rampa_ascendente(rho_c,u_c,v_c,T_c,p_c,Nxi,Neta); e_c(Nxi,2:Neta-1) = c_v*T_c(Nxi,2:Neta-1); %------------------------------------------------------- %% Encode variables for result and calculation of result U_c(:,:,1) = rho_c; U_c(:,:,2) = rho_c.*u_c; U_c(:,:,3) = rho_c.*v_c; U_c(:,:,4) = rho_c.*(e_c + (u_c.*u_c + v_c.*v_c)/2); U = (U_p + U_c)/2; %------------------------------------------------------- %% Check for convergence if max( max( abs( U(:,:,1) - rho ) ) )>1e-8 conver = 0; else conver = 1; end %------------------------------------------------------- %% Decode final variables rho = U(:,:,1); u = U(:,:,2)./U(:,:,1); v = U(:,:,3)./U(:,:,1); e = U(:,:,4)./U(:,:,1) - (u.*u + v.*v)/2; T = e./c_v; p = R*rho.*T; p(:,1) = 2*p(:,2) - p(:,3); [rho,u,v,T,p] = Boundary_Conditions(rho,u,v,T,p,Nxi,Neta); nu = nu0*(T/T0).^(3/2).*((T0+110)./(T+110)); k = nu*c_p/Pr; % Update time and iteration counts iter = iter + 1; t = t + dt; % Print at screen to check for correct execution display([t,dt]) if max(max(abs(imag(T))))>0 && flag==0 tbreak = t; flag = 1; end % mesh(y,x,p) % pause(0.001) end %---------------------------------------------------------- 56 %% Check the validity of the results by calculating the mass inflow and outflow balance Mass_in = trapz( y(1,:) , u(1,:).*rho(1,:) ); Mass_out = trapz( y(Nxi,:) , u(Nxi,:).*rho(Nxi,:) ); if abs( ( Mass_in - Mass_out ) / Mass_in*100 ) <= 1 display('The solution seems valid checking the mass balance') else display('The solution seems to violate the mass balance, it is not valid') end %---------------------------------------------------------- %% Viscosity force at the wall f_pared = nu(:,1).*( u(:,2) - u(:,1) )/dy; F_pared = trapz(0:dx:LHOR,f_pared); C_d = F_pared/0.5/rhoinf/uinf^2/LHOR; F_press = trapz(0:dx:LHOR,p(:,1)); %---------------------------------------------------------- %% Plots and ending plot( x(:,1) , p(:,1) ) figure(2) mesh(0:dy:LVER,0:dx:LHOR,real(p)) disp = [( Mass_in - Mass_out ) / Mass_in*100 conver] time = toc % Additional plots and calculations using the results may be done at % different scripts 63 As a final conclusion it can be said that the code programmed obtains valid results and is ready to analyse more complex situations. 4.1.1.2 Dependency with the grid points One of the tests that can be made now is for example how the results vary with the fineness of the grid. So, several simulations where made of the flat plate with Mach 2, only varying the number of grid points, for the results to be compared. Figure 4.11. Trailing edge temperature profile, Mach 4, Grid 70x70. Figure 4.12 . John D. Anderson’s “Computational Fluid Dynamics” [Ref 1] results. Mach number profile at trailing edge with grid 70x70 and Mach 4. 64 What has been found is that at the trailing edge the results are almost the same. For instance, the pressure profiles are shown in Figure 4.14. Simply the shock wave moves up slightly. But at the leading edge the results change significantly. The pressure spike becomes more intense, starts earlier and is thinner with the increase in number points. The resistance coefficient rises a little bit, and the total vertical force, in Newtons, is almost unchanged. Figure 4.14 . Trailing edge pressure profiles, Mach 2, dependency with the grid points. Figure 4.13 . Surface pressure, Mach 2, dependency with the grid points. 65 40 x 40 50 x 50 60 x 60 70 x 70 80 x 80 100x100 120x120 150x150 𝐂 𝐝 0.0569 0.0589 0.0604 0.0616 0.0626 0.0643 0.0655 0.0670 V ertical force (N) -1.383 -1.386 -1.388 -1.389 -1.389 -1.390 -1.390 -1.389 4.1.1.3 Dependency with the Mach number It is also interesting to make the analysis of how the Mach number affects the results. For instance, looking at Figure 4.15 and comparing it with Figure 4.1, it is clear that the shock wave has been made narrower. Also, the pressure spike is much greater, and other calculations like the forces change too. The maximum value of the pressure in the shock wave at the trailing edge rises with the Mach number according to the Figure 4.18. Also, as mentioned before, the point of this maximum gets closer to the flat plate, in what seems to be an asymptotic behaviour according to Figure 4.17. Figure 4.15. Pressure (Pa) field, Grid 70 x 70 points, Mach number 10. 66 The pressure at the plate is considerably altered by the Mach number, shown in Figure 4.16. Figure 4.18. Trailing edge maximum shock wave pressure value. Figure 4.17 . Trailing edge maximum shock wave pressure position. Figure 4.16. Pressure at the surface of the plate 67 The resistance coefficient decreases, even though the actual friction force increases, due to the increase in the external velocity that is used in the normalization. The downward vertical force increases due to the increase in the overall pressure on the plate surface. 4.1.2 Adiabatic Plate The second variant of the boundary conditions is the case of the adiabatic plate. This time, the temperature of the plate is not constant, instead there is no heat transfer crossing through it. The results obtained by John Anderson for the adiabatic plate have already been shown in figures Figure 4.5, Figure 4.8, Figure 4.9 and Figure 4.12. Figure 4.20 . Resistance coefficient. Figure 4.19 . Downward vertical force (N). 68 For this section not to be too large only the case of Mach number 4 and 70 x 70 grid points is going to be discussed to compare it with the results obtained in [Ref 1]. To begin with, the surface pressure obtained by John Anderson in Figure 4.5 and the one obtained in this project are identical. The pressure peak is slightly higher than five time the outside value and the behaviour, with an increase in pressure around the grid point 10 and a value of about 2.8 times the outside pressure. Figure 4.21. Surface pressure for adiabatic plate, Mach 4, Grid 70 x 70 points. Figure 4.22. Trailing edge pressure profile for adiabatic plate, Mach 4, Grid 70 x 70. 69 The trailing edge pressure profile is also the same in Figure 4.22 than in Figure 4.8, the temperature profile is the same in Figure 4.23 than in Figure 4.9, and finally the Mach number profile is the same in Figure 4.24 than in Figure 4.12. To close the adiabatic plate case the temperature field is displayed in Figure 4.25. The air is firstly heated by the shock wave and then heats more and more when approaching the surface. Figure 4.24. Trailing edge Mach number profile for adiabatic plate, Mach 4, Grid 70 x 70. Figure 4.23 . Trailing edge temperature profile for adiabatic plate, Mach 4, Grid 70 x 70. 70 4.1.3 Flat plate with compressed grid Not only the program used for rectangular and uniform grids needs testing, also the used for non-uniform grids does. To achieve this what has been done is solving the flat plate but now using a compressed grid. In particular, the grid transformation is given by the following relations: 𝑥=𝑒−1 𝑓 Equation 4.2 𝑦=𝑒−1 𝑓 Equation 4.3 Theses relations create a compressed grid around the point (0,0), where f is a corrector factor that determines the severity of this compression. Obtaining the analytical direct and inverse metrics from Equation 4.2 and Equation 4.3 is easily achieved: 𝜕𝑥 𝜕𝜉=𝑒 𝑓 ; 𝜕𝑦 𝜕𝜂=𝑒 𝑓 ; 𝜕𝜉 𝜕𝑥=𝑓 𝑓𝑥+1 ; 𝜕𝜂 𝜕𝑦=𝑓 𝑓𝑦+1 Figure 4.25. Temperature (K) field for adiabatic plate. 71 With the rest been equal to zero. As a demonstration, the Figure 4.26 shows the grid generated with this method for a factor of 10 and 15 x 15 points. The only thing that is left to say is that with this grid the results obtained by the “MainGridTrans_DirectMetrics.m” program have been exactly the same as the obtained by the uniform grid, so they are not going to be repeated here. By checking that the results coincide, the correct functioning of the program developed for non-uniform grid transformations has been ensured. 4.1.4 Auxiliary functions for the flat plate Initial Conditions Function (both for uniform and non-uniform grid) function [rho,u,v,T,p] = Initial_Conditions_FlatPlate(rhoinf,uinf,Tinf,R,Nxi,Neta) rho(1:Nxi,1:Neta) = rhoinf; u(1:Nxi,1) = 0; u(1:Nxi,2:Neta) = uinf; v(1:Nxi,1:Neta) = 0; T(1:Nxi,1:Neta) = Tinf; p(1:Nxi,1:Neta) = R*rho.*T; Boundary Conditions Function (both for uniform and non-uniform grid) function [rho,u,v,T,p] = Boundary_Conditions_FlatPlate(rho,u,v,T,p,Nxi,Neta) % CC at the right side rho(Nxi,2:Neta-1) = 2*rho(Nxi-1,2:Neta-1) - rho(Nxi-2,2:Neta-1); u(Nxi,2:Neta-1) = 2*u(Nxi-1,2:Neta-1) - u(Nxi-2,2:Neta-1); v(Nxi,2:Neta-1) = 2*v(Nxi-1,2:Neta-1) - v(Nxi-2,2:Neta-1); T(Nxi,2:Neta-1) = 2*T(Nxi-1,2:Neta-1) - T(Nxi-2,2:Neta-1); p(Nxi,2:Neta-1) = 2*p(Nxi-1,2:Neta-1) - p(Nxi-2,2:Neta-1); % CC at the plate surface Figure 4.26. Compressed grid of 15x15 points and f = 10  72 rho(2:Nxi,1) = 2*rho(2:Nxi,2) - rho(2:Nxi,3); p(2:Nxi,1) = 2*p(2:Nxi,2) - p(2:Nxi,3); T(2:Nxi,1) = T(2:Nxi,2); Geometry Function (only for compressed grid) function [x,y,dxi,deta,dxidx,dxidy,detadx,detady,J,dxdxi,dxdeta,dydxi,dydeta] = Geometry_CompressedGrid(Nxi,Neta,LHOR,LVER) f = 10^6; xi_LHOR = log(f*LHOR+1); eta_LVER = log(f*LVER+1); dxi = xi_LHOR/(Nxi-1); deta = eta_LVER/(Neta-1); for j=1:Neta xi(:,j) = 0:dxi:xi_LHOR; end for i=1:Nxi eta(i,:) = 0:deta:eta_LVER; end for i=1:Nxi for j=1:Neta x(i,j) = (exp(xi(i,j)) - 1)/f; y(i,j) = (exp(eta(i,j)) - 1)/f; end end dxidx = f./(f*x+1); dxidy = 0; detadx = 0; detady = f./(f*y+1); dxdxi = exp(xi)/f; dxdeta = 0; dydxi = 0; dydeta = exp(eta)/f; J = dxdxi.*dydeta - dydxi.*dxdeta; 79 4.3 Rearward Step The complementary problem is to see what happens in the opposite situation, when the step is now facing rearwards to the incoming flow. In this geometry, the flat plate above the solid square creates a shock wave just as the flat plate without the step would do. However, after that, the pressure falls significantly behind the step. The geometry has been delimited to 10 in the horizontal direction and to 1.5⋅10 in the vertical direction. This problem is less interesting than the previous one as the only meaningful thing that is to be measured is the pressure behind the step. The more interesting behaviour, that discusses the formation of the Prandtl-Meyer Expansion waves, is when the step becomes a slope and constitutes the next application. The pressure field for a Mach 2 flow with a step height of 30% of the total vertical dimension, 4.5⋅10, is shown in Figure 4.35. The pressure increase in the shock wave created by the flat plate above the solid square is identical to the normal flat plate. However, after that the pressure falls in an expansion wave to a value lower than the external pressure. Figure 4.35. Pressure (Pa) in rearward step, Mach 2, step (3e-6 x 4.5e-6). Twodimensional view. 80 The peak of low pressure happens just at the corner of the step, where it exists a point that almost reaches zero (and it has been a nightmare for the numerical stability of the program). After the peak, the pressure starts to recover both when advancing down close to the surface or back following the flow direction. It is seen clearly in Figure 4.36, that represents the pressure for a Mach 4 flow with a rearward step of height 3⋅10. The evolution of the mean pressure in the vertical wall is displayed in Figure 4.37. As it is showed by the two lines, now the mean pressure at the wall, in contrast with the forward step (Figure 4.33), highly depends on the step height. When the step is higher the overall mean pressure is also higher. That is probably because for the lower step the peak of low pressure represents more proportionally to the total vertical surface. Figure 4.36. Pressure (Pa) in rearward step, Mach 4, step (3e-6 x 3e-6). Twodimensional view. Figure 4.37. Mean pressure in (Pa) for rearward step. Dependency with step height and Mach number. 81 Also, a dependency is appreciated regarding the Mach number. With higher flow velocities, the vacuum effect grows behind the wall. This effect is stronger for higher steps and seems to have diminishing behaviour with the increase in Mach. The last piece of information that is needed to be discussed is the temperature field. Now it is not similar to the pressure field as it was for the forward step. Instead, the air is heated in two different regions, as it is seen in Figure 4.38 for a Mach 2 flow with a 3.5⋅10 step. The first region of hot air is caused by the shock wave, but after it, it appears another region in which the air is heated by friction. The proof of that can be obtained by taking a look at the velocity field. In Figure 4.39 it is seen that the heated region coincides with great velocity gradients caused by friction. Figure 4.38. Temperature (K) in rearward step, Mach 2, step (3e-6 x 3.5e-6). Twodimensional view. Figure 4.39 Velocity (m/s) in rearward step, Mach 2, step (3e-6 x 3.5e-6). Twodimensional view. 82 Also, it is appreciated that a recirculating air region is created at the proximity of the bottom corner of the vertical wall, as isolating and analysing only the horizontal velocity field, showed in Figure 4.40, reveals that this component reaches negative values and has a circular behaviour. The rearward step is a good application to test a demonstrated fact; that the stagnation enthalpy is not affected in a shock wave. The stagnation enthalpy has the following expression: ℎ=𝑐𝑇+𝑉 2 As it is shown in the figure below, the presence of the shock wave seen in Figure 4.35 does not affect the stagnation enthalpy. Only the viscosity effects in the high velocity gradients do. Stagnation enthalpy in rearward step, Mach 2, step (3e-6 x 4.5e-6). 4.3.1 Auxiliary functions for the rearward step Initial Conditions Function Figure 4.40 . Zoom at the bottom corner, horizontal velocity (m/s) in rearward step, Mach 2, step (3e-6 x 3.5e-6). Twodimensional view. 83 function [rho,u,v,T,p] = Initial_Conditions_RearwardStep(rhoinf,uinf,Tinf,R,Nx,Ny) % Step X direction beginning Nxs = floor(0.3*Nx); % Step height Nya = floor(0.2*Ny); u(1:Nx,1:Ny) = uinf; u(1:Nx,1:Nya) = 0; v(1:Nx,1:Ny) = 0; rho(1:Nx,1:Ny) = rhoinf; T(1:Nx,1:Ny) = Tinf; p(1:Nx,1:Ny) = R*rho.*T; Boundary Conditions Function function [rho,u,v,T,p] = Boundary_Conditions_RearwardStep(rho,u,v,T,p,Nx,Ny,Tinf,pinf) % Step X direction beginning Nxs = floor(0.3*Nx); % Step height Nya = floor(0.2*Ny); % Force variables in the solid square u(1:Nxs,1:Nya) = 0; v(1:Nxs,1:Nya) = 0; T(1:Nxs,1:Nya) = Tinf; p(1:Nxs-1,1:Nya-1) = pinf; % CC at the right side rho(Nx,1:Ny-1) = 2*rho(Nx-1,1:Ny-1) - rho(Nx-2,1:Ny-1); u(Nx,1:Ny-1) = 2*u(Nx-1,1:Ny-1) - u(Nx-2,1:Ny-1); v(Nx,1:Ny-1) = 2*v(Nx-1,1:Ny-1) - v(Nx-2,1:Ny-1); T(Nx,1:Ny-1) = 2*T(Nx-1,1:Ny-1) - T(Nx-2,1:Ny-1); p(Nx,1:Ny-1) = 2*p(Nx-1,1:Ny-1) - p(Nx-2,1:Ny-1); % CC at the solid surface rho(Nxs:Nx,1) = 2*rho(Nxs:Nx,2) - rho(Nxs:Nx,3); %rho_p(2:Nx,1) = rho(2:Nx,1) - dt/2/dy*( 3*rho(2:Nx,1).*u(2:Nx,1) - 4*rho(2:Nx,2).*u(2:Nx,2) + rho(2:Nx,3).*u(2:Nx,3)); p(Nxs:Nx,1) = 2*p(Nxs:Nx,2) - p(Nxs:Nx,3); % CC at the Nxs: vertical wall rho(Nxs,1:Nya) = 2*rho(Nxs+1,1:Nya) - rho(Nxs+2,1:Nya); p(Nxs,1:Nya) = 2*p(Nxs+1,1:Nya) - p(Nxs+2,1:Nya); % CC at the Nya: horizontal wall rho(1:Nxs,Nya) = 2*rho(1:Nxs,Nya+1) - rho(1:Nxs,Nya+2); p(1:Nxs,Nya) = 2*p(1:Nxs,Nya+1) - p(1:Nxs,Nya+2); 84 4.4 Rearward Slope and The Prandtl-Meyer Expansion Waves It is time to use the non-uniform grid transformation for the first time, turning the rearward vertical wall into a slope. In this problem the horizontal and vertical dimensions of the computational domain have been set to 10 and 1.63⋅10 respectively. The beginning and height of the slope can be modified at will, and the slope falls from the indicated point till the end of the computational domain. Firstly, it is going to be studied the case in which still exist a portion of flat plate above the solid surface. In this case the flat plate creates the typical shock wave that has been seen in previous applications. This pressure field showed in Figure 4.41 for a 35º slope with Mach 3 shows precisely that. The temperature field on the other hand, Figure 4.43, reveals that the second region of heated air, at least in this configuration, remains attached to the slope. That may be an indication that the recirculating air region has not been formed. In addition, the temperature field reveals that there exists a region in which the temperature falls below the external value. This is probably because there the air is accelerated due to the expansion wave to velocities higher than the incoming flow. Figure 4.41. Pressure (Pa) in rearward slope, Mach 3, Slope angle 35º. Two-dimensional view. 85 Looking at the velocity field, Figure 4.42, it is confirmed that, as it was hinted by the temperature, there is no recirculating area. Also, it is seen that at the area in which the temperature falls, the velocity of the air is accelerated to even a greater value than the external flow. Figure 4.43. Temperature (K) in rearward slope, Mach 3, Slope angle 35º. Two-dimensional view. Figure 4.42. Velocity (m/s) in rearward slope, Mach 3, Slope angle 35º. Threedimensional view. 86 It is therefore interesting to see the surface pressure distribution and its dependency with the slope angle and the Mach number. The Figure 4.44 shows the pressure surface profiles for different angles of the slope. The first curve corresponds to the flat plate and coincides with the results of section 4.1. It is seen that the presence of a minimum slope already affects greatly to the results after the slope corner. Then, the increase in the slope angle decreases the pressure in a diminishing fashion with the angle. The behaviour of the pressure at the surface after the step is asymptotical and is not expected to differ much more after the end of the computational domain. Figure 4.44. Pressure at surface for rearward slope, Mach 3, dependency with slope angle. Figure 4.45. Pressure at surface for rearward slope, slope angle 35º, dependency with Mach number. 87 When varying the Mach number, what is observed is that the pressure peak at the beginning is accentuated, as was expected due to Figure 4.16, and what happens after that is that, at least for an angle as big as 35 degrees, the pressure after the corner does not depend much on the Mach number. 4.4.1 Prandtl-Meyer expansion waves The most interesting part of this application is when the flat plate surface is eliminated, letting only the slope. In this case, if the slope angle is big enough, there is no shock wave created. It is interesting to investigate the transition between the behaviour of a flat plate and this other. This Figure 4.46 represents the pressure for a Mach three flow with an angle of 26º. Now the higher values of pressure correspond to the unperturbed flow, and the plate only creates an expansion wave. It is known that it is an expansion wave because it decreases the pressure of the air and accelerates it. Figure 4.46 . Pressure (Pa) in rearward slope, Mach 3, Slope angle 26º, no flat plate. Figure 4.47. Velocity (m/s) in rearward slope, Mach 3, Slope angle 26º, no flat plate. 88 Looking at Figure 4.47 it is seen that the air accelerates over 100 m/s with respect to the unperturbed air stream in the expansion wave. The elimination of the shock wave also removes the heated air area created by it. Instead, the air cools in the expansion wave. Close to the surface the air heats again due to friction effects, as it corresponds to high velocity gradients observed in Figure 4.47. Validating our results is very difficult as no exactly similar numerical experiment has been found, but it can be tried by comparing them with the results obtained by John Anderson, again in “Computational Fluid Dynamics” [Ref 1], in chapter 8. There, the analytical results for the expansion wave problem is demonstrated. The set-up and the problem in general that is solved in this chapter is quite different, so differences in the results are to be expected. Specifically, the problem solved in the reference has a free slip condition at the surface, and is treated with a step marching method, so its dimensions are unreachable by our code, of over 10 meters. However, recreating the 5,352 degree angle and the Mach 2 flow in this project’s code, we obtain the following results. The density at the surface at end of the computational domain is 1.071 kg/𝑚, that compared with the 0.992 kg/𝑚 of the analytical solution represents a 7.3% difference. Counting on the differences in problems mentioned, it is reasonable to present the results obtained as not far from reality. Making the same comparison it is obtained an 8.7% difference in pressure value, 0.81⋅10 Pa compared to the analytical 0.734⋅10 Pa, again reasonably far away given the differences in set up. Figure 4.48. Temperature (K) in rearward slope, Mach 3, Slope angle 26º, no flat plate. 95 4.5 Leading Edge: Study of Shock Wave Angle It is well known for the study of aerodynamic applications, and is comprehensibly explained by A. B. Ripoll and Miguel Pérez-Saborid in “Fundamentos y aplicaciones de la mecánica de fluidos” [Ref 3], that a leading edge that has a certain angle, frontally receiving a supersonic flow, creates an oblique shock wave to adequate the flow to its solid geometry. In this section this problem is going to be treated and the results obtained compared with the analytical solutions for oblique shock waves. Before anything else, a description of the geometry is given. It consists of an edge that has a certain angle facing the stream of air. The horizontal dimension is set to 10 𝑚 and the edge is situated at half this distance. The vertical dimension for the computational domain is 1.63⋅10 𝑚 in the positive direction and the same value in the negative direction. The flow created is symmetrical, as can be seen in Figure 4.56 for a Mach 2 flow with an edge semi-angle of 40 degrees, as the problem geometry and initial and boundary conditions also are. Figure 4.56. Pressure (Pa) for leading edge, Mach 2, edge semi-angle 40º. 96 The surface is set at constant temperature. That creates an effect seen in Figure 4.57, in which the air is not heated if is arrives exactly from the symmetrical line, while the rest of the flow is affected by the shock wave and heated substantially. However, not always happens that the oblique shock wave begins before the leading edge. For a Mach three flow with a 25 degrees semi-angle, the shock wave begins exactly at the leading edge. When that happens, it is said that the shock wave is anchored to the edge. It is interesting to know the pressure that the leading edge creates at the surface and its dependency with the Mach number and the edge semi-angle. For a Mach two flow the Figure 4.57. Temperature (K) for leading edge, Mach 2, edge semiangle 40º Figure 4.58. Pressure (Pa) for leading edge, Mach 3, edge semiangle 25º. 97 pressure profiles at the surface of the edge, plotted with respect to the distance in the X direction, are shown in Figure 4.59. There is a change of behaviour between 20 and 25 degrees. With low angles the pressure rises at the beginning of the surface, signalling that the shock wave is anchored. After that the pressure at the surface starts high and only falls, a sign that the shock wave is no longer anchored. Figure 4.59. Pressure at the surface for a leading edge, Mach 2, dependency with edge semi-angle. Figure 4.60. Pressure at the surface for a leading edge, Mach 3, dependency with edge semi-angle. 98 Looking at figures Figure 4.60 and Figure 4.61 that show the same pressure profiles for Mach three and four respectively, it is clear that this behaviour also depends on the Mach number. For Mach three and four, at a semi-angle of 30 degrees, the shock wave is still anchored. Also, it is appreciated that a peak of pressure followed by a peak of suction is accentuated by both the angle and the Mach number. At the beginning of the edge for Mach four, with the higher values of semi-angle, the flow seems to become turbulent. From [Ref 3] it is known that the analytical solution for the angle of an oblique shock wave depends on the angle of the flow after the wave and the Mach number with the following equation. 𝑡𝑎𝑛 𝛼 = 2 𝑐𝑜𝑡 𝛽 󰇩 𝑀 𝑠𝑖𝑛 𝛽 − 1 𝑀 ( 𝛾 + 𝑐𝑜𝑠 2𝛽 ) + 2󰇪 Equation 4.4 Where 𝛾 is the specific heat ratio and 𝑀 is the Mach number of the incoming flow before the shock wave. 𝛼 is the angle of the air stream after the shock wave, which in this case is the same as the edge semi-angle, and 𝛽 the angle that the shock wave forms with respect to the horizontal direction. Figure 4.61. Pressure at the surface for a leading edge, Mach 4, dependency with edge semiangle. 99 Looking at the pressure fields to locate the shock wave, the angle 𝛽 can be measured in each configuration that has been tried in this project. When plotting the lines corresponding to Equation 4.4, and plotting the pairs of values of edge semi-angle and shock wave angle created measured, the result is Figure 4.62. There, it can be seen that the results obtained match reasonably well with the theoretical curves. Specially for Mach three and four at high angles the results fit almost perfectly. However, when decreasing the edge semi-angle the accuracy decreases significantly, most notably when going under ten degrees, where the angle of the oblique shock wave is overestimated by the code. This also happens with Mach two, where overall the shock wave angle is overestimated, although the accuracy also grows with the leading edge semi-angle. The Figure 4.63 shows that the value of the mean pressure at the surface, and thus the resistant force applied to the leading edge, is highly dependent on the edge semi-angle. Moreover, the dependency is accentuated when increasing the Mach number. Figure 4.62. Oblique shock wave angle with respect to leading edge semiangle and Mach number. Comparison between analytical and numerical results. 100 When comparing the Figure 4.63 to Figure 4.33, it can be seen that the leading edge creates much lower pressures, about one order of magnitude lower, than the forward step. However, it is expected that their behaviours become similar, if not identical, when the edge semi-angle approaches 90 degrees. However, to calculate the pressure resistance that the leading edge suffers it is necessary to take into account the surface length and the surface angle. It is plotted at Figure 4.64, where it can be seen that the resistance approaches zero when the angle also does, an expected behaviour. Figure 4.63. Mean pressure at the surface for leading edge, dependency with Mach number and leading edge semiangle. Figure 4.64. Backward pressure resistance force for leading edge, dependency with Mach number and leading edge semi-angle. 101 Also, it is clear that the higher the Mach number the higher the resistance that the leading edge suffers, as well as happens when increasing the angle. However, when calculating the resistance coefficient, the behaviour with α remains, but the one with the Mach number is reversed, because of the increased velocity that takes part in the adimensionalization. 4.5.1 Auxiliary functions for the leading edge Geometry Function function [x,y,dxi,deta,dxidx,dxidy,detadx,detady,J,dxdxi,dxdeta,dydxi,dydeta] = Geometry_LeadingEdge(Nxi,Neta,LHOR,LVER) xi_LHOR = LHOR; eta_LVER = LVER; dxi = xi_LHOR/(Nxi-1); deta = eta_LVER/(Neta-1); for j=1:Neta xi(:,j) = 0:dxi:xi_LHOR; end for i=1:Nxi eta(i,:) = 0:deta:eta_LVER; end % Edge X direction beginning s = floor(0.5*Nxi); xis = LHOR*(s-1)/(Nxi-1); % Edge semi-angle (\alpha) a = 30*pi/180; Figure 4.65. Backward pressure resistance force for leading edge, dependency with Mach number and leading edge semiangle. 102 y0(1:s) = 0; y0(s+1:Nxi) = ( xi(s+1:Nxi) - xis )*tan(a); for i=1:Nxi for j=1:Neta x(i,j) = xi(i,j); end y(i,:) = linspace(y0(i),LVER,Neta); end % Direct metrics DC dxidx = Dx(xi,x); dxidy = Dy(xi,y); detadx = Dx(eta,x); detady = Dy(eta,y); % Inverse metrics DC dxdxi = Dx(x,xi); dxdeta = Dy(x,eta); dydeta = Dy(y,eta); dydxi = Dx(y,xi); J = dxdxi.*dydeta - dydxi.*dxdeta; Initial Conditions Function function [rho,u,v,T,p] = Initial_Conditions_LeadingEdge(rhoinf,uinf,Tinf,R,Nxi,Neta) % Beginning of the edge in the X direction s = floor(0.5*Nxi); rho(1:Nxi,1:Neta) = rhoinf; u(1:Nxi,1:Neta) = uinf; u(s:Nxi,1) = 0; v(1:Nxi,1:Neta) = 0; T(1:Nxi,1:Neta) = Tinf; p(1:Nxi,1:Neta) = R*rho.*T; Boundary Conditions Function function [rho,u,v,T,p] = Boundary_Conditions_LeadingEdge(rho,u,v,T,p,Nxi,Neta) % Beginning of the edge in the X direction s = floor(0.5*Nxi); % CC at the right side rho(Nxi,2:Neta-1) = 2*rho(Nxi-1,2:Neta-1) - rho(Nxi-2,2:Neta-1); u(Nxi,2:Neta-1) = 2*u(Nxi-1,2:Neta-1) - u(Nxi-2,2:Neta-1); v(Nxi,2:Neta-1) = 2*v(Nxi-1,2:Neta-1) - v(Nxi-2,2:Neta-1); T(Nxi,2:Neta-1) = 2*T(Nxi-1,2:Neta-1) - T(Nxi-2,2:Neta-1); p(Nxi,2:Neta-1) = 2*p(Nxi-1,2:Neta-1) - p(Nxi-2,2:Neta-1); % CC at the surface 103 rho(2:Nxi,1) = 2*rho(2:Nxi,2) - rho(2:Nxi,3); p(2:Nxi,1) = 2*p(2:Nxi,2) - p(2:Nxi,3); u(2:s-1,1) = 2*u(2:s-1,2) - u(2:s-1,3); 104 5 Conclusions and Further Developments Reaching the end of the project, it is safe to say that the overall objectives of the project, which were developing a custom program that could numerically solve supersonic flows for both uniform and non-uniform grids and obtaining those results for a few interesting problems, have been met. One of the conclusions that can be extracted from the hours of work put into this project is that the numerical stability of the program is a serious problem for the expansion of it. In general, it can be said that the Courant-Friedrichs-Lewy (CFL) criterion (Equation 3.1) is not enough to preserve a convergent behaviour of the numerical solution, specially when sharp edges like the ones that appear in the steps are involved. In those cases, a strong correction made with the safety factor had to be made in order to obtain any solutions. In this project it was considered the addition of numerical artificial viscosity to alleviate this problem, as John D. Anderson discusses in chapter 6.6 of [Ref 1], but finally it was discarded as it would have unnecessarily complicated the project, and in the end was not necessary to obtain the results. Other conclusion reached is that the computation times can be efficiently optimised using the program structure. One may wonder why the programs showed in chapter 3 are written the way they are, with almost no auxiliary functions that perform small tasks like the differentiation, which could have avoided some repetitions, also minimising the probability of human error and making easier the creation of the program. Well the truth is that at the beginning they were made this way. However, when any function was called it had to make a copy of the variables it worked with, and then copy them back once it finished. This made the original programs painfully slow, and with successive testing and improvements the computation time was reduced in some cases to even a half of the original time. Talking about the results I am very happy to see that they correspond with reasonable accuracy with the external results obtained by other authors or with the analytical solutions they have been compared to. It creates a profound satisfaction when something that has required so much time and effort finally works, and the solutions obtained are exactly what they were supposed to be, like what happened when comparing the results of the flat plate. Because it has not been an easy task. I would say that this project is the single work that I have invested more time into of all my life. When developing the programs and the project overall, a lot of focus has been put into explaining how the things work piece by piece, and making everything, specially the programs, as easy to understand and manage as possible. The objective of this is to make it easier for a future student to take this project as basis and expand the work with other ideas or analysing other problems using the programs developed here. That was one of the key motivations of creating the chapter 3 to explain the programs, and to make them modular, with the information of the specific problems restricted to the auxiliary functions, which are interchangeable. In the future I would like to see another student using the programs to generate and solve problems that due to time constrains have not been solved here, like the case of a rectangle or