Development and Validation of an Open Source CFD Code for Analysis of Aerospace Vehicles

In this work, a review of the theoretical aspects and an assessment to validate a Computational Fluid Dynamics (CFD) open-source code for applications in aerospace problems are discussed. Th e code uses a fi nite volume method, with cell-centered implementation, and it is suitable for simulations of inviscid, laminar, and turbulent fl ows. Th e code considers two-dimensional cases with unstructured meshes and employs the turbulence model known as Spalart-Allmaras. Th e implementation is detailed presenting the spatial discretization, including the upwind scheme, the linear reconstruction algorithm, and the calculation applying the method of gradients. Th e temporal discretization considers the application of a multistage explicit algorithm using a 5 stages Runge-Kutta method. Th e validation was done considering three cases of study: the inviscid shock tube, the laminar fl at plate, and the fl ow over a rocket fairing. Th ese cases are simulated using the so ft ware developed and the results are compared with analytical and experimental results. Th e rocket fairing case is related to the analysis of the Brazilian VLS launch during its transonic fl ight and it exempli fi es the e ff ect of the shock wave/boundary-layer interaction in its pressure distribution. Th e simulation results present a good agreement with the experimental results.


INTRODUCTION
e aerospace applications of aerodynamics are a challenging area.In some cases, it is possible to use simpli ed methodologies like Newton's method or tangent-wedge methods to model aerospace vehicles (Rolim et al. 2020).But it is not always possible and, therefore, the implementation of more complex methods based in Computational Fluid Dynamics -CFD -are necessary.ese simulations can be highly complex when, for example, an entire aircra is modeled in order to study the e ect of the fuselage in the vortex dynamic (Sustrino et al. 2020) or even when an entire rocket is modeled with the grid using chimera technique (Oliveira Neto et al. 2011).It is also possible to consider an intermediary situation in that CFD is applied in a more simpli ed and accessible way. is corresponds to the analysis of two-dimensional cases.
e 2-D approach reduces the computational costs while keeping several fundamental characteristics of the problem.ere are several studies that attest this fact.For example, Azevedo and Korzenowski (2009) have presented the modeling of hypersonic inlet where the sensibility of the solver to several characteristics were analyzed in order to evaluate the con gurations with best performances.McNamara et al. (2010) has presented an approximate modeling of hypersonic aeroelasticity using the ow over a double-wedge airfoil.
e results of several simpli ed approaches were compared with the results obtained from CFD simulations considering a bi-dimensional model.e work of Silva and Pimenta (2019) has presented an analysis of aerodynamic properties of a reentry capsule using a twodimensional model.In another study, Farrokhfal and Pishevar (2014) developed a CFD optimization method for airfoils based on the adjoint method.e objective was to reduce the compressibility drag or pitch moment of transonic airfoils without compromising its li coe cient.All these cases show the importance of bi-dimensional analyses in the research of aerospace problems.
e CSU (Code of Simulation for Unstructured meshes) is a code of CFD developed to be used in tests of new models and to verify the results from other CFD so ware, considering bi-dimensional problems.e CSU is based on the nite volume method, using the density-based formulation and unstructured meshes.e methods used are common methods of CFD and the objective of this study is not to contribute with originality, but to present a review of the methodology used to develop this kind of code and make a computational tool available for the scienti c community.
e so ware was developed to be an accessible open-source code that allows the fast assimilation and use.e implementation of the computationally expensive parts of the algorithm is done in C and the managing of the compilation and data visualization is done in Python.e so ware has been released in Linux and can be obtained at https://github.com/CarlosCHMS/CSU. e implementation of CSU makes use of the OpenMP package to parallelize the calculation processing.e use of parallelization reduces the time of processing, which is quite helpful for a CFD code.OpenMP allows it to subdivide the process into several threads, which is e cient if this number is smaller than the number of physical processor cores.Special attention must be given to the memory access since processing in di erent threads and trying to write in regions of a vector close together can result in errors and reduction of the performance.e visualization of the results corresponds to plotting the graphics of ow contours and data.e contours are plotted using specialized algorithms based on unstructured meshes.However, it is necessary to know the values at the cell nodes not at the cell-centers.is requires the cell-centered values obtained from the simulation to be converted to node-centered values and it is made using the inverse of distance interpolation method (Tasri and Susilawati 2021).
SPATIAL DISCRETIZATION e complete system of Navier-Stokes equations constitutes the governing equations of the ow and are presented in the integral form in Eq. 1 (Blazek 2005): (1) where F c is the vector of convective uxes, F v is the vector of viscous uxes, Q are the source terms, Ω is the control volume and the vector W is called "vector of conservative variables".In the bi dimensional approach with a one equation, the turbulence model has 5 components as presented in Eq. 2: (2) where ρ is the density, u is the velocity in the x direction, v is the velocity in the y direction, E is the total energy per unit of mass and η _ is the modi ed turbulent viscosity, that is part of the Spalart-Allmaras implementation, and it is used in the turbulent cases.
e code uses the nite volume method to solve the governing equations.e basic idea is to divide the ow domain in a set of small control volumes and, in each of them, the system of Eq. 1 is solved.For each of these volumes, the Eq. 3 is solved: where the index I refers to each of these volumes, Ω l is the volume of the control volume (area in 2-D formulation) and R l is the residuals that are calculated with the numerical integration of the convective and viscous uxes over the surface of the volume.
e source terms will be considered null for all variables, except for the modi ed turbulent viscosity.
In general, the discretization of the volume can be made using a structured or unstructured approach.e structured meshes have the advantage of being easier to be generated in simple domains and of having a more direct implementation in the code.
However, these advantages do not overcome the di culties in generating meshes for more complex domains.Indeed, it can be said that the maturation of unstructured mesh technologies has ushered CFD into a new era. is advancement enables an increased speed of mesh generation, along with the ability to automate speci c processes and employ adjoint-based mesh adaptation and shape optimization techniques (Mani and Dorgan 2023).us, the CSU was developed to use unstructured meshes.
e mesh generation is performed using the gmsh so ware (Geuzaine and Remacle 2000).e mesh format used as input is .su2. is mesh format was developed for the use of the SU2 so ware, which is an open-source CFD solver also currently in development and that has its characteristics detailed in Economon et al. (2016).is format was chosen because it was considered a format simple to read.
ere are two alternatives to de ne the control volumes for an unstructured mesh: cell-centered scheme and cell-vertex scheme.
In the cell-centered scheme, the control volumes are the same as the mesh cells and the values of the ow variables refer to the centroids of these cells.For the cell-vertex scheme, the variables are associated to the mesh vertices and the control volume is, in general, formed as some sort of combination among the cells that share the vertex.Both methods are similar in terms of accuracy and computational work but the cell-centered scheme has advantages in non-conformal grids (Blazek 2005).e cell-centered scheme was chosen in CSU so ware because its implementation was considered more direct and simpler.e calculation of the ow volumes and areas is made considering the alternatives of planar or axially symmetric solution.In the axisymmetric case, the symmetry axis is considered as the x axis.
For the calculation of the residual in the boundary of each mesh cell, the ux vectors must be presented.e convective ux vector can be written accordingly to Eq. 4: (4) where V is the component of the velocity normal to the cell face, n x and n y are the components of the vector normal to the cell face, p is the static pressure and H is the total enthalpy per unit of mass.In fact, there are Eq. 5 and 6: (5) (6) and for perfect gases, the Eq.7: (7 where γ is the ratio of speci c heats. It is also important to observe that velocities and the moment fluxes must be considered with components normal and tangential to the cell face.This requires a transformation that consists in applying rotations before and after the flux calculation.That rotations make the flux coherent to the flux splitting method that will be applied (Mulder 1989).Making the rotation in the way that the u velocity component and moment component stay perpendicular to the cell border, the Eq. 4 became the Eq.8: (8) Equation 8 together with a simple interpolation of variables in the cell face is not enough to calculate the ux correctly.It is well known that simple central schemes for convective terms are not stable.It is necessary to include an arti cial dissipation to the central scheme or use an upwind scheme.In this work, the upwind scheme was chosen since it delivers a much better resolution of shocks and boundary layers than the central schemes (Blazek 2005).
e upwind schemes are also divided in several options of ux-vector splitting or ux-di erence splitting schemes like Roe, HLLE, Osher, Leer, AUSM, among others.But each one of them has its limitations.Roe, for example, diverges at strong expansions even if entropy x is used and it su ers from the carbuncle phenomenon, that consists of a numerical instability in capturing shock waves in a multidimensional domain.e HLLE does not incorporate information about contact discontinuities which makes it too dissipative.Osher's scheme fails in near vacuum condition.Leer scheme is very dissipative in 1-D contact discontinuities and generates glitches in the pressure near the edge of the boundary layer.AUSM (Advection upstream splitting method) is a quite interesting scheme since it can calculate strong shock waves.However, it generates numerical overshooting behind the shocks (Wada and Liou 1994).In fact, a variant of AUSM seems to be the best choice in terms of performance.Wada and Liou (1994) proposed the AUSMDV that has some advantages like: high-resolution at contact discontinuities; conservation of enthalpy for steady ows and numerical e ciency. is scheme consists in a mixture of two other variants of AUSM that are the AUSMD and AUSMV.e AUSMD is a ux-Di erence-splitting-biased scheme at the time that AUSMV is a ux-Vector-splitting-biased scheme.e convective ux coe cients are calculated as presented in the Eq.9: (9) where the ψ assume the values 1, v, H and η-, for each of the sides le L or right R of the cell face.e value of (ρu) 1 2 is calculated using the Eq.10: (10) where the values of u + L and u - R are calculated using a series of equations that takes into account sound velocity, normal velocity to the cell, pressure, and density in both sides of the cell face.
e only exception to Eq. 9 is the component referring to the moment normal to the cell face which is calculated using Eq.11: (11) where (ρu) 1 2 AUSMV corresponds to the calculation using the AUSMV scheme and (ρu) 1 2 AUSMD the procedure related to the AUSMD scheme, S is a switching function calculated using a relation of pressures in the le and right side of the border and ρ 1 2 is calculated using a system of equations similar to the used for u + L and u - R .e de nition of (ρu 2 ) 1 2 AUSMV and (ρu 2 ) 1 2 AUSMD are presented in Eq. 12 and 13: (12) (13) e characteristic of the switch function was de ned from numerical study cases based on the shock-tube problem. is analysis showed that AUSMV has a higher shock-capturing capability, however this scheme can result in strong oscillations in some cases.
e switching function proposed is de ned by Eq. 14: (14 where p R and p L are the pressures in the right and le side of the cell face and K is a constant parameter de ned as 10.More details on the equations used in the implementation of the AUSMDV scheme can be obtained in the (Wada and Liou 1994).e values of ow variables to be considered on the le and right sides of the cell face must be determined, considering the ow variables constant inside the cell result in solutions that are only rst-order accurate.First-order-accurate solutions present a good convergence, however they are too much dissipative.us, it is necessary to achieve the second order accuracy.It can be done by using piecewise linear reconstruction of the solution (Barth and Jespersen 1989).
e basic idea of the method is to assume that the ow variables vary linearly inside the cell.So, the values at the cell faces are given by the Eq. 15 and 16: where U L and U R represent the values on the le and right side of the cell face.ese values are the ones used in the calculation of the ux.U l and U J are the ow variables at the center of the cells l and J. U l and U J are the gradients of the ow variables calculated for each cell-center.r L and r R are the vectors from the cell-centers to the center of the cells faces.Finally, ψ l and ψ l are the values of the limiters.e gradients calculation is not only necessary for linear reconstruction, but it is used in the calculation of viscous uxes, making it a key feature of applying the nite volume method.e most popular methods for calculating gradients are Green-Gauss and least-squares.
e Green-Gauss method is based in the theorem of Green in 2-D (or Gauss theorem in 3-D) and it requires the numerical integration of the variable over the cell's surface.e least-square method is based on linear regression using the cell and its neighbors.ese methods are popular mainly because they are not dependent on a particular cell's geometry and can be used in cells with arbitrary numbers of faces.It is quite convenient for the application in unstructured meshes as the ones used in CSU. e choice between these two methods is guided by the ndings of Syrakos et al. (2017), who analyzed their performance across several mesh cases.eir conclusion is that, for arbitrary grids, the least-squares method is rst-order accurate, while the Green-Gauss method is zero-order accurate.Higher accuracy orders can be achieved for speci c grids.Considering these results, the least-squares method was chosen for implementation in CSU.
e limiter function ψ is crucial to the implementation of the second order upwind algorithm.Limiters avoid the wiggles that occur in the solution near shocks or discontinuities.e basic idea is that the limiter is close to 1 in regions of smooth solution, which results in a second-order scheme.In regions of intense variation of the ow variables, the limiter is close to zero, which results in a rst-order solution.e Venkatakrishnan's limiter was chosen to be implemented because of its superior convergence properties (Blazek 2005).e limiter can be de ned accordingly to Eq. 17 (Venkatakrishnan 1993): (17) where the Eqs.18-21 de ne the additional terms: U max and U min are the values maximal and minimum of the ow variable U considering the cell I and all its neighbors.e min l in Eq. 16 means that the value of the limiter must be calculated considering all possible values of ( U l • r L ) for each cell face and choosing the smaller one.e parameter є avoid a possible division by 0 in Eq. 17 and can be used to improve the convergence of the solution.
Until this point the formulation was focused in the calculation of the convective uxes, however there is also the viscous uxes.e current version of the so ware has the implementation of the viscous uxes relative to the laminar and turbulent ow.
e viscous ux is given by the Eq.22: (22 where the Eqs.23-29 de ne its components: where μ is the e ective dynamic viscosity coe cient and k is the e ective thermal conductivity.e terms τ xx and τ yy are the normal stresses and τ xy is the shear stress.e terms Θ x and Θ y correspond to the work of viscous stresses plus the heat uxes.e terms τ x and τ y are associated to the components of the viscous uxes related to the turbulence model.It is possible to observe that it is necessary to calculate the gradient of the primitive variables u, v, T and η _ .e value η laminar is the laminar kinematic viscosity.e least-square method can be used for this, but the results are related to the cell-center.But, to calculate the uxes, it is necessary to obtain the values of the gradients in the cell faces. A rst approach to calculate the gradient in the cell face could be to use the average of the gradient values of the neighboring cell to estimate the gradient in the faces.However, Haselbacher and Blazek (2000) have shown that this simple approach results in the decoupling of the solution on quadrilateral meshes.e solution proposed to tackle this problem is to interpolate the gradient according to Eq. 30: where the Eq.31-33 de ne its components: e value of U refers to the velocities, the temperature are inputs of the equations and r IJ is the radii vector that goes from the center of the cell I to the center of the cell J.

TURBULENCE MODELING
e ows considered in the modeling of aerospace vehicles are, in most cases, turbulent.erefore, it is necessary to include an approach in the CSU code to deal with these problems.ere are several common ways to tackle this problem, such as the use of turbulence models associated with Reynolds-averaged Navier-Stokes equations (RANS), the use of Large Eddy Simulation (LES), or even the use of Direct Numerical Simulation (DNS).However, the last two methods have a higher computational cost and generally require three-dimensional domains, which are not available for CSU in its current state of development.In fact, the use of RANS and turbulence models is the preferred choice for an initial approach to modeling turbulence.e de nition of RANS equations can be obtained from Rodriguez (2019).It involves modeling the turbulent ow by assuming that the ow variables can be represented as a sum of mean and oscillatory components.e ow is then solved by considering only the mean terms and variables.However, it is still necessary to determine the values of the mean terms that appear in the RANS equations, and this is where turbulence models come into play.
ere is a wide variety of turbulent models in the literature.ey are classi ed as Zero, One, or Two-equation models, depending on the number of ow variables introduced in the RANS equations.According to Blazek (2005) where Pr laminar is the laminar Prandtl number assumed as 0.72 and Pr turb is the turbulent Prandtl number assumed as 0.9.e values of the constants present in the model are c b1 = 0.1355, σ = 2/3, c b2 = 0.622 , k = 0.41, c w2 = 0.3 and c w1 that is calculated by Eq. 45: (45)

TEMPORAL DISCRETIZATION
For the implementation of the temporal discretization there are two possible classes: the explicit and the implicit methods.e advantage of the implicit methods is the possibility of achieving faster and stable solutions by de ning a higher CFL number.
However, the explicit methods are easier to implement and parallelize.Explicit schemes generally require more calculation steps than the implicit ones, but each iteration is relatively cheap in computational terms (Leer et al. 1989).Because of these advantages, an explicit method was chosen to be implemented in the CSU.e explicit method chosen uses a multistage scheme.e multistage method corresponds to the application of the Runge-Kutta method in which the new solution su ers several updates in a sequence of stages.is implementation can be seen in the equations of 5 stages presented in Eq. 46, below: (46) Equation 46show that in each stage, a preliminary solution W J (i) is calculated and then used to calculate the residual R J (i) of the next stage.At the end of the ve steps, the nal solution W J n+1 is calculated.It is important to observe that this implementation of Runge-Kutta requires that only the initial solution and the newest preliminary solution need to be stored in the memory.
It represents an advantage in terms of computational e ort.
Using a multistage scheme is convenient since it improves the stability of the spatial discretization at the same time that allows an increase in the time-step used.e ve order scheme applied to the so ware allows an increase of 2.5 times in the time step.e coe cients α can be de ned considering rst-or second-order spatial discretization schemes.However, for strong shocks, rst-order schemes are preferred (Blazek 2005) since they prevent oscillations in the solution.e coe cients α i implemented in the so ware are obtained from Leer et al. (1989) and are given in Table 1.e time step is de ned using the approach proposed by Vijayan and Kallinderis (1994) and modi ed by Blazek (2005).
e time step is given by Eq. 47: (47) where Δt l is the local time step and is given by Eq. 48: (48) e Λ are the spectral radii.e convective spectral radii are given by Eqs.49 and 50: where c is the velocity of sound in the cell, dŜ x and dŜ y are the projections of the cell volume in the direction x and y respectively.
ey can be calculated by Eqs.51 and 52: (51) (52) e summation corresponds to the sum over all the faces of the cell.e expressions for the viscous spectral radii Λ x vl and Λ y vl are more complex and can be found in Blazek (2005).

BOUNDARY CONDITIONS
e so ware supports four types of boundary conditions: wall, inlet, outlet and symmetry.e wall boundary condition is also divided in two cases, the inviscid and the viscous cases.In the inviscid case, all terms of the convective ux (Eq.8) are null unless there is pressure.e pressure on the wall is calculated using the formulation presented by Whit eld and Janus (1984) that uses the characteristic boundary conditions.ey use concepts of method of characteristics to estimate the pressure on the wall from the available ow variables in the closest cell-center.e expression of pressure is given by Eq. 53: e convective ux used for the viscous wall is similar to the one used for the inviscid wall.e main di erence is that now there are the viscous uxes to be calculated.In order to do that, the gradient of the velocities must be calculated considering the velocity in the wall boundary null.at corresponds to the no slip condition on the wall.e heat ux is another condition to be de ned in the viscous wall.ere are some possibilities as: constant temperature wall, radiative wall or adiabatic wall.e adiabatic wall corresponds to the case in which no exchange of heat occurs through the wall, which implies in considering the heat ux or the gradient of temperature in the Eq.23 and 24 to be zero.e adiabatic wall is the condition currently implemented in CSU. e boundary conditions for subsonic inlets and outlets have also been implemented.e implementation follows the same approach proposed by Whit eld and Janus (1984) and is also based on the characteristic method approach.For the supersonic cases, the implementation di ers.e input values are directly inserted as boundary values for the supersonic inlet boundary.For the supersonic outlet, the conditions at the boundary correspond to the ow variables of the cells adjacent to the boundary.is approach ensures that information ows upwind, ensuring the stability of the code.
e nal boundary condition is symmetry, which assumes that the boundary divides the ow in two symmetrical and opposite regions.e implementation is done by assuming the existence of dummy cells.is approach supposes the existence of auxiliary cells outside the boundary with ow variables that are used in the calculation of the uxes.For the symmetry condition, the ow variables at the adjacent wall are considered accordingly, Eq. 54: (54) where the velocity u is considered normal and v tangent to the cell face.Using the ow variables of the cell and its dummy counterpart, it is possible to calculate the ux in the face using the same approach used for the interior cells.
When considering the turbulent model, the condition at the wall is η _ = 0 and for the inlet the conditions is η _ = 5η laminar .ese conditions agree with what is proposed in Allmaras et al. (2012) for a fully turbulent Spalart-Allmaras model.

NUMERICAL SIMULATIONS AND VALIDATION
e cases chosen for code validation have the objective of evaluating the behavior of the so ware considering di erent aspects.
e shock tube problem can be used to verify the performance of the inviscid solver and the accuracy of the rst and the second order algorithm.e at plate problem has the objective of evaluating the performance of the laminar viscous implementation.
Finally, the rocket fairing cases have as objective to evaluate the so ware globally, considering inviscid and turbulent solutions and comparing them to experimental results.

Sod shock tube
A typical shock tube is a tube generally made of metal with a circular or rectangular cross-section, in which two regions with gas at low pressure and high pressure are separated by a diaphragm.Under certain conditions, the diaphragm bursts and produces a shock wave that propagates in the low-pressure region.A schematic representation of the problem is presented in Fig. 1.
is problem is closely related to the Riemann problem and is commonly used to test CFD algorithms.e problem can be solved using the Euler equations in one spatial dimension; however, it can also be solved using a rectangular tube in two dimensions.(55) e mesh used covers a rectangular domain measuring 50m x 1m and is generated using 3022 triangular elements.e division of the initial regions occurs at the position of 25m.e simulation results are compared with the analytical solution obtained by applying the method of characteristics.e transient simulation is conducted until a physical time of 0.02 seconds. is time was chosen to be su ciently long to observe the discontinuities in the solution well separated, yet short enough to avoid wave re ections with the walls at the tube's extremities, which could complicate the solution using the characteristics method.
e impact of using a rst or second-order method can be observed in Fig. 2. While the rst-order method presents some discrepancies regarding the characteristic solution in the shock regions, the second-order solution exhibits good agreement across all regions.ese results align with expectations for this type of simulation and demonstrate the so ware code's e ectiveness.
Additionally, it should be noted that some preliminary results of this problem exhibited overshoots in the second-order solutions near the discontinuities.is issue was resolved by reducing the value of the parameter є in Venkatakrishnan's limiter in Eq.
18.A value of є = 0.001 resulted in solutions without overshoots.is observation is consistent with Venkatakrishnan (1993), which states that the increasing of є can lead to discrepancies in the solutions.Flat plate e laminar ow on a at plate is a convenient problem to verify the implementation of the viscous ux since the Blasius analytical solution can be used for comparisons.
e problem is set as a ow with velocity u = 34.75m/s or Mach number of 0.1.e pressure is set as p = 1 • 10 5 Pa and the temperature is T = 300K.e domain is rectangular of 0.5 m x 0.006 m and the mesh has 2,786 elements.e results are presented in Figs. 3 and 4.
Fig. 4 on the left, shows a good agreement for the velocity profile at x = 0.5m obtained from the CFD code and the Blasius solution.The Fig. 4 on the right, also shows a good agreement of the shear stress in the region far from the plate leading edge, however the curves do not fit well for values lower than 0.07 m.This disagreement is related to the fact that, close to the leading edge, the boundary layer has a thickness of the same order of the height of the cell adjacent to the wall.Because of this, the solver is not able to resolve the boundary layer close to the edge and incorrect results are obtained.
The validation will be made by using wind tunnel results related to the VLS (Veículo Lançador de Satélites), a former Brazilian project of vehicle launchers.The experiments were conducted at the ONERA (Office National d'Etudes et de Recherches Aérospatiales), France, using a 1/15 scale model (Gauthier, 1989).The CFD simulations considered the same geometry of the scale model used in the wind tunnel experiments.
e experiments generated pressure coe cient c p distributions for several cases of Mach and Reynolds numbers.However, a detailed analysis of these results permitted the veri cation of two cases that have a remarkable e ect of interaction between shock and boundary layer.ese cases are: rst one Mach = 0.9 and R є = 1.68 • 10 6 and the second one Mach = 1.05 and R e = 1.57• 10 6 .It is desirable to use these cases in the validation process because they allow the observation of the e ect of the boundary layer that is linked to the turbulence model implemented.e Reynolds numbers were calculated using the vehicle diameter as reference length.e explanation of the e ects observed in these two cases is as follows: the case of Mach = 0.9 is transonic, with two near normal shocks being formed a er the expansion corners in the fairing.ese two shocks are strong enough to generate observable e ects in the distribution of C p when compared with the inviscid case.e case of Mach = 1.05 has a unique shock interacting with the boundary layer at the end of the fairing.Curiously, greater Mach numbers do not generate remarkable interactions.e reason is that, despite the increase in the external ow Mach number, the pressure ratio in the shock at the end of the fairing is reduced, resulting in weak shocks and a less pronounced interaction with the boundary layer.
e mesh created has 87,210 triangular and quadrilateral elements.e quadrilateral elements are used close to the surface of the vehicle that has as objective to improve the solution of the boundary layer.e quadrilateral layer has a growth rate of 1.2 and was carefully made in order to result in y+~1 in the wall.e meshes can be seen in Fig. 5. e domain is rectangular with the following boundary conditions: the le side being the input, the right side being the output, the upper side being wall and the lower side being the symmetry and wall.e walls are considered adiabatic.
Figure 6 shows the Mach number contours obtained from the solution for Mach 0.9, considering the turbulence model, alongside a Schlieren photograph of the ow under the same conditions.It is possible to observe the formation of a shock wave caused by the ow acceleration, which becomes supersonic a er the rst corner.A er that, the shock happens and makes the ow subsonic again.It is also remarkable that a relatively small region of higher Mach closes the second corner.is is associated with a second expansion through the corner and a second shock.In the Mach contours, it is also possible to observe a detail of the boundary layer in the shock region.
e comparison of the contours and the photography allows the conclusion that the ow patterns were captured by the CFD simulation.A er the shocks, there is an increase in pressure that propagates upstream through the subsonic part of the boundary layer.
is alteration in the ow around the shock region leads to a smoother pressure transition compared to what would be expected in an inviscid shock.is clari es the di erences obtained by the inviscid solutions and shows the importance of incorporating a turbulence model in these simulations.In fact, the agreement observed is not exclusive to the Spalart-Allmaras model and could be obtained by any other turbulence model that correctly represents the turbulent boundary layer.However, the results permit the conclusion that the Spalart-Allmaras turbulence model, as implemented in CSU, was su ciently accurate to represent the ow and match the experimental results.Additionally, Fig. 8 illustrates the temperature distribution on the fairing obtained from CSU using the Spalart-Allmaras model.
It reveals a peak near the initial position, corresponding to the stagnation point on the nose, followed by a valley representing the circular part of the nose.Around x = 0.1 m, there is another valley attributed to the rst fairing corner, which corresponds to an expansion fan.e curves remain relatively close until approximately x = 0.12 m, a er which they begin to diverge due to di erences in shock con gurations for Mach 0.9 and 1.05 cases.Even with these di erences, the curves approximate each other again a er the end of the fairing at x = 0.23 m.It is also remarkable that the total range of temperature variation was around 16 K.

CONCLUSIONS
e study presented the theory involved in the implementation of the code CSU and its validation through some test cases.
e spatial discretization was detailed using the AUSMVD scheme combined with the application of least-squares for gradient calculation and Venkatakrishnan's limiter.e turbulence modeling based on Spalart-Allmaras is presented.e temporal discretization using multistage scheme based in the application of the Runge-Kutta method is also detailed.
In terms of validation of the so ware, CSU was applied to three problems.e rst was the shock tube problem and the results were compared with the analytical results obtaining an excellent agreement considering the second order algorithm.e simulation of the at plate was made considering the laminar solver, which resulted in a good agreement with the Blasius solution.Finally, the solver using inviscid and Spalart-Allmaras approach is used in the simulation of the ow over a rocket fairing.e CFD results agreed with the experimental data and permitted the validation of the implementation of the turbulent model.
, some of the most popular models include the Spalart-Allmaras one-equation model, the kє two-equation model, and the Menter SST (shear stress transport) two-equation model.In their analysis, Zingg and Godin (2009) evaluated di erent turbulence models and their applicability, concluding that the Spalart-Allmaras model performs best for aerodynamic ows.It provides accurate results for attached and mildly separated ows.However, ows with large-scale separation may require the use of the Menter SST model.Given this, the Spalart-Allmaras model is chosen for implementation in CSU, with the possibility of a future implementation of the Menter SST model.e turbulence model is based on the one proposed by Spalart and Allmaras (1992) with the alterations proposed for a compressible implementation presented in Allmaras et al. (2012).e compressible version of the governing equation is given by Eq. 34: (34)where the production and destruction terms are de ned by Eqs.35 and 36, respectively: laminar suppression terms are not considered.e variable S is the modi ed vorticity and it is calculated by using the modi cation proposed byAllmaras et al. (2012) in order to prevent negative values.e function f w is de ned by Eqs.37-39: viscosity and thermal conductivity are calculated by using the Eqs.43 and 44: is calculated using the values of pressure p d , velocity u d and v d related to the cell adjacent to the border.e values of ρ 0 and c 0 are reference values of density and sound velocity, respectively, and are assumed from the cell adjacent to the border.

Figure 1 .
Figure 1.Schematic representation of the shock tube problem. Figure is not in scale.e initial conditions used in the problem were obtained from Sod (1978) and correspond to what is called Sod shock tube problem.ey correspond to assume in the region 1 ρ 1 = 1.0 and p 1 = 1.0 of the reference values, and for the region 2 the values are ρ 2 = 0.125 and p 2 = 0.1 of the reference values.e initial velocity is equal to zero for both regions.e CSU code works with atmospheric air and requires a di erent set of variables in the input data.e inputs are the initial values of pressure, temperature and both velocity components.Considering the Sod inputs and the reference values of pressure and temperature of 1.0 • 10 5 Pa and 300K, the initial conditions for the two regions of the 2-D case are given by Eq. 55:

Figure 2 .
Figure 2. Results from the simulation of the shock tube for the instant of 0.02s.(a) Mach number along the shock tube; (b) density along the shock tube.
by the authors.

Figure 3 .
Figure 3. Velocity contours for the at plate laminar ow.Out of scale gure.

Figure 4 .
Figure 4. Results of the at plate simulation.(a) e velocity pro le on the at plate at position x=0.5 m; (b) e shear stress along the at plate.
by the authors.

Figure 5 .
Figure 5. Mesh used in the simulations: (a) the entire domain; (b) a detail of the quadrilateral layer on the fairing tip.

Figure 6 .
Figure 6.Comparison between computational and experimental results: (a) contours of Mach number for the case of Mach 0.9 for the fairing considering the turbulence model; (b) the corresponding Schlieren photography.

Figure 7
Figure7presents a comparison between the simulation and experimental results for the pressure coe cient, cp.Up to x = 0.15 m for Mach 0.9 and x = 0.2 m for Mach 1.1, both inviscid and viscous solutions match the experimental data.However, beyond these positions, disparities emerge between the inviscid solutions and experimental results, while the results obtained with the Spalart-Allmaras turbulence model continue to exhibit good agreement.e cause of these discrepancies can be attributed to the occurrence of shocks in those regions and their interaction with the boundary layer.is phenomenon is well documented inHoughton and Carpenter (2003) and has been observed in a rocket fairing byMata et al. (2017).
by the authors using data fromGauthier (1989).

Figure 7 .
Figure 7. Distribution of c p on the fairing surface.(a) case of Mach number 0.9.(b) case of Mach number 1.05.

Figure 8 .
Figure 8. Temperature distribution on the fairing surface.Results from the application of CSU with the Spalart-Allmaras turbulence model.

Table 1 .
Coe cients α i for rst order implementation.
Source: Elaborated by the authors.