As already demonstrated by different authors, the TaylorGalerkin (TG) scheme, in the context of the Finite Element Method (FEM), is particularly suitable for the solution of supersonic flows. The TG scheme, using hexahedral finite elements with analytical evaluation of element matrices, is applied in this work. Tools to avoid locking and a shock capturing technique for the solution of supersonic viscous and nonviscous compressible flows are also employed. However, TG scheme usually presents instabilities in subsonic flows. Even in cases in which the free stream Mach number corresponds to supersonic flows, there are always flow regions, specifically near the walls of the immersed obstacles, where the speed is lower, and the local Mach number corresponds to a subsonic flow. The Characteristic Based Split (CBS) scheme was developed in order to obtain a single method to improve the behavior with respect to TG method in subsonic and supersonic regimes. In the last two decades some works have shown advantages in convergence rates of the CBS method when compared to the TG algorithm. However, simulation time increases in the CBS method since split operations, typical of this algorithm, imply in additional element loops. In this paper a hybrid algorithm called ModifiedTaylorGalerkin scheme (MTG) is proposed. This algorithm presents advantages with respect to TG and CBS schemes in terms of convergence properties and computational processing time. In order to get an efficient algorithm, the element matrices are analytically integrated. This is performed with two different approaches. In the first approach the inverse matrix and the determinant of the Jacobian matrix at element level are evaluated with a reduced integration form, using the point located in the center of the element for mass, convective, diffusive and stabilization element matrices; all these matrices are integrated analytically. In the second approach, mass and convective matrices are calculated by a complete integration scheme (including the inverse matrix and the determinant of the Jacobian matrix at element level in the analytical expression to be integrated) and the diffusive and stabilization matrices are calculated with a reduced integration form, using the point located at the center of the element to calculate the inverse matrix and the determinant of the Jacobian matrix at element level. Finally, this work incorporates the SpalartAllmaras (SA) turbulence model using a conservative version of the transport equation, as proposed by the authors of the original SA model in a later paper. Algorithms are tested to determine convergence rate improvements in both laminar and turbulent cases and for different Mach numbers (supersonic, transonic and subsonic flows).
Keywords: Compressible turbulent, FEM, TaylorGalerkin, CBS, SpalartAllmaras.
Different authors have shown that, in the context of FEM, the TG scheme is particularly suitable for the solution of compressible flows involving supersonic shock waves at high Mach numbers [12]. Burbridge and Awruch [3] applied this scheme using hexahedral finite elements with analytical evaluation of the element matrices but evaluating the inverse matrix and the determinant of the Jacobian matrix at the center of the element, together with a shock capturing technique for the solution of viscous and nonviscous compressible laminar flows [4]. However, for subsonic regimes (with M <0.9) or for incompressible flows, the TG scheme presents certain instabilities, as reported by Zienkiewicz and Taylor [5].
Even in supersonic flows there are always some regions of subsonic flow near the walls of immersed obstacles. It is therefore necessary to improve the convergence of the algorithms for low Mach numbers, even when they are applied to flows where Mach number at the freestream is supersonic.
To address this problem the CBS scheme was developed by some researchers, such as Zienkiewicz and Codina [6], to obtain a single algorithm with a good behavior, both in subsonic and supersonic regimes. During the last two and a half decades this algorithm was established as an acknowledged tool for the computation of a wide spectrum of incompressible or compressible flows at different Mach numbers [710]. The work presented by Zienkiewicz et al. in Reference [11] shows the explicit form of the algorithm and its performance for subsonic, transonic and supersonic flows.
However, in spite of convergence rates improvements, it is observed that computational processing time increases due to split operations involved in the CBS scheme. For this reason, an alternative hybrid algorithm is presented here, looking for improvements in convergence properties with respect to the TG algorithm without an increase of the computational processing time.
It is important to mention that, as in previous works [3,12], an analytical integration of the element matrices was used here, evaluating the inverse matrices and the determinants of the Jacobian matrices at the center of the element. But in this work, an alternative approach was also applied, integrating mass and convective matrices in an exact analytical form (including the inverse matrices and the determinants of the Jacobian matrices in the analytical expressions) using Maxima symbolic resolution program (http://maxima.sourceforge.net/). Stabilization and convective matrices were integrated in reduced form (evaluating the inverse matrices and the determinants of the Jacobian matrices at the center of the elements), as performed in [3,12].
Finally, the turbulence model of SpalartAllmaras (SA) was implemented here in order to test the behavior of these schemes, and their convergence properties, in the context of the FEM for turbulent flows. The form of the SA model implemented here was proposed for the first time by Allmaras et al. [13] in 2012, using the transport equation for the turbulent modified kinematic viscosity variable in conservative form.
TG, CBS and MTG schemes are presented here with some cases of laminar and turbulent flows, comparing their performances.
In an Eulerian description, the dimensionless system of partial differential equations governing the fluid dynamics problem in conservative form, and filtered with a Favre filter (FANS), may be written as follows:
(1) 
where
(2) 
The dimensionless form has been obtained using a reference length, the freestream density and the freestream speed of sound. In Eq. (2):
The viscous stress components of a Newtonian fluid are given by:
(3) 
where and are the molecular viscosity and the eddy viscosity, respectively, is the turbulent kinetic energy and the subscript runs from 1 to 3.
The components of the conduction heat flux vector are:
(4) 
where and are the thermal conductivity and the turbulent thermal conductivity respectively and u is the specific internal energy.
The terms in the energy equation are given by the following expression:
(5) 
where is a constant which depends on the turbulence model. Since the SA turbulence model does not include the calculation of turbulent kinetic energy, in this work, in Eq. (2) and the third term in Eq. (3) are neglected.
The terms in the SA transport equation are given by the following expression:
(6) 
The transport equation also includes the additional term , expressed by
(7) 
where and are constants of the SA turbulence model.
Finally, the production and destruction terms, in the conservative transport equation of SA are given by
(8) 
where the first term represents the production and the second term D the destruction, being _{,} , and constants and variables of the SA turbulence model and is the smallest distance from a field point (or node) to the nearest solid boundary (wall). The transport equation used here is expressed in its conservative form (Allmaras et al. [13]). Here the tripterm and the term of laminar suppression, both included in the original expressions of the model, are not considered.
To close the system of equations it is still necessary to define all parameters associated with the SA turbulence model; they are given by

(9) 
where is the magnitude of the angular velocity vector, and

(10) 
The turbulent viscosity and turbulent thermal conductivity are defined as follows:
(11) 
where is the turbulent Prandtl number, which is adopted as and and are the specific heat coefficients at constant volume and at constant pressure, respectively.
The state equation for a perfect gas (air) is given by:
(12) 
The internal energy u and the temperature are related to the independent field variables by the following expression:
(13) 
where is the temperature, being neglected in the SA model.
The Sutherland law is used in this work to establish the dependence of viscosity and thermal conductivity with respect to temperature (dimensionless internal energy). This law can be expressed as follows:
(14) 
Values of and may be found in reference [14] for the Sutherland law based on temperature; they must be multiplied by and divided by the square of the freestream sound speed to obtain the dimensionless form, which was used in this work.
Initial and boundary conditions must be applied to Eq. (1) in order to obtain a mathematical model correctly defined. The forced boundary conditions (or Dirichlet boundary conditions) are given by:
(15) 
where , and are prescribed values of the velocity vector components, the density (specific mass), the specific internal energy and the variable , at the boundary surfaces , , and Г respectively.
The natural boundary conditions (or Neumann boundary conditions) are defined by the following expressions:
(16) 
In Eq. (16)
For nonviscous flows the transport equation of the variable and the diffusive terms in Eq. (1) should be omitted. In this case, only normal velocity components to the solid boundaries are prescribed with null values and only the normal components of the surface tractions on the outlet boundaries are considered.
In this section the wellknown TG, CBS and Single Step CBS algorithms are mentioned but the respective equations are not shown for the sake of brevity. The reader can look at the references [1,3,5,6,8] to get more details. Instead, the full matrix equations for the proposed MTG scheme are presented in the section 4.
To discretize Eq. (1) in time, the unknown vector is expanded in Taylor series omitting higher order terms. Then the space discretization is done using the classic BubnovGalerkin technique in the context of the finite element method. The GreenGauss theorem is applied to terms containing second order derivatives (terms with third or higher order derivatives are neglected) to weaken the demands with respect to continuity of the element interpolation functions and their derivatives, obtaining the matrix equations [13].
A conditionally stable explicit time integration version of the TG scheme [4,5,15], may be obtained working with a lumped mass matrix (where only elements of the main diagonal are different from zero).
Also, it is important to point out that the increments of the flux variables obtained for the current time step with the TG scheme depend of flow variables at the previous time step and at the current time step . So, an iterative scheme to solve terms involving flow variables at the new time step is required [1,3].
The CBS method consists of four steps [6]. In the first step, intermediate values of the conservative flow variables of the momentum equations are calculated, omitting the terms containing pressure gradients. In the second step, the continuity equation is solved to obtain the density increment. In the third step, the flow variables of the momentum equations are updated. In the fourth step, the energy equation and the SA transport equation in its conservative form are solved. Finally, pressure is calculated using the state equation.
Once the temporal discretization is performed, spatial domain discretization is carried out using the classic BubnovGalerkin method. The GreenGauss theorem is applied to terms containing second order derivatives to weaken the demands with respect to continuity of the element interpolation functions and their derivatives.
This method may be applied with implicit, semiimplicit and explicit time integration schemes [11]. An explicit scheme was implemented in this work in order to compare its behavior with respect to other explicit algorithms presented here.
Matrix expressions obtained with the CBS method may be found in references [5,6,8].
An alternative of the CBS method consists in the elimination of the split in the momentum equation in order to obtain a onestep algorithm. Some authors [5,6] suggest that, although the singlestep CBS scheme saves computational processing time with respect to the classical CBS method, it is not suitable to simulate flows with high Mach numbers. The explicit version of the singlestep CBS was implemented in this work. The common matrix expressions obtained with singlestep CBS scheme can be found in reference [5].
It may be observed that stabilization problems arise when TG scheme is used, and they are linked with density and pressure values, being these variables basically controlled by the continuity equation. The TG method is an efficient technique in terms of computational processing time, but schemes such as CBS have much better convergence properties for problems involving a wide range of Mach numbers. In order to link both features (computational efficiency and good convergence properties), a modification of the TG scheme is proposed. The new scheme is called Modified TaylorGalerkin (MTG) scheme, which consists to change the stabilization term on the continuity equation by a similar expression, but defined in terms of pressure, as used in the CBS method. Since the TG scheme has good convergence rates for supersonic flows, the modification is only implemented in regions where the local Mach number is less than one. For this reason, a blending function is proposed such that, the continuity equation discretized in time results:

(17) 
where

(18) 
The blending function depends of the Mach number, which is evaluated locally. Momentum and energy equations as well as their iterative terms are not modified, remaining as in the original TG scheme.
Spatial discretization is performed using the classical BubnovGalerkin method. The GreenGauss theorem is applied to terms containing second order derivatives to weaken the demands with respect to continuity of the element interpolation functions and their derivatives.
The following matrix expressions are obtained for the element vectors of independent unknowns’ increments:

(19a) 

(19b) 

(19c) 

(19d) 
An iteration procedure is required for those righthand terms evaluated on the current time step , where the subscript indicates the previous iteration and corresponds to the current one. Element matrices are indicated below (with ):

(20a) 
The element vectors are:

(20b) 
where , are the velocity components, pressure and internal energy element vectors, containing their corresponding eight nodal values. The quantities between braces are element vectors built with the nodal values of the variables inside the braces. The boundary vectors are:

(20c) 
In these expressions and are the element volume and its boundary surface respectively; expressions with the subscript are element averages; = [Ф_{1}, Ф_{2}_{, …} Ф_{8}] is the row matrix containing the interpolation functions for each local node and is the row matrix containing the interpolation functions evaluated at the boundary surface, is the consistent mass matrix, while is the lumped mass matrix and is a vector where all its components are equal to one. Terms that involve turbulent kinetic energy have been neglected.
The variables with superscript are evaluated on previous time step and the variables with superscript are evaluated at current time step. Boundary vectors appearing in the SA transport equation and in the iterative terms of all the equations were neglected.
Once these equations are assembled for all the elements of the mesh and the corresponding boundary conditions are applied, the nodal values of , and can be calculated at each time step using an iterative scheme. Thus, the nodal values of the components of velocity, the specific total energy and the variable can be obtained immediately. Then the nodal values of the specific internal energy are calculated and, using the state equation, the nodal pressure values are also obtained.
When viscous terms are not considered, terms containing matrices , and as well as vectors and and the transport equation corresponding to the turbulence model are omitted.
Stabilization matrices for TG, CBS and MTG schemes are evaluated using a local (internal) time step . For the TG scheme there are no appreciable differences in using a local internal time step or a single time step [16].
Considering that these algorithms were implemented in their explicit form, they are conditionally stable and the stability condition of CourantFriedrichsLewy (CFL) for each element, in dimensionless form, is used to ensure stability. A safety coefficient is adopted with values [3].
To capture the strong discontinuities and eliminate high frequency oscillations near shock waves, a wellknown method is used to add artificial viscosity selectively in regions where highpressure gradients are observed [17]. An artificial damping coefficient is used to control de amount of artificial viscosity.
The convergence of the iterative process is obtained when the following conditions are simultaneously satisfied:

(21) 
where is an index of nodes ( ), is the total number of nodes of the finite element mesh and is the tolerance (in this paper was adopted). The solution process ends when the variable (time) reaches a value previously established by the user, or when the steady state is reached. This last condition is considered fulfilled when the following expression is satisfied:

(22) 
In this work was adopted.
Eight nodes isoparametric hexahedral elements were used in this work. Figures 1 and 2 show the element in the physical and computational spaces, respectively. The interpolation functions of this element are given by the following expressions:

(23) 
where and are the natural coordinates of the node . The derivatives of the interpolation functions with respect to the global coordinates are given by:

(24) 
being the components of the inverse of the Jacobian matrix , which can be expressed as follows:

(25) 
Figure 1. Physical space 
Figure 2. Computational space 
Here is the vector containing global coordinates of the eight nodes of each element, is the determinant of the Jacobian matrix and is the differential of the element volume.
By substituting Eqs. (23), (24) and (25) on the element matrices and vectors defined by the TG, CBS and MTG schemes, integrals in the computational domain are obtained. These integrals are commonly calculated by means of a numerical integration scheme to obtain components of the element matrices that will be later assembled to get the global system of equations. However, in this work, matrices were calculated by means of analytical integration using the Maxima symbolic resolution software. This implies an improvement in computational processing time efficiency.
In a first approach the inverse matrix elements and determinants of the Jacobian matrix were evaluated at the center of the element, where and these values were used to integrate matrices corresponding to variables time derivatives, convective, diffusive and stabilization terms [3,12]. Special attention must be given to elements distortion and an hourglass control technique to avoid spurious modes on diffusive terms are also necessary [12].
Alternatively, in a second approach, the mass matrix (representing variables time derivatives) and matrix (representing convective terms) were integrated analytically, including the inverse matrix elements and determinant of the Jacobian matrix in the analytical expressions (without evaluating these terms at the center of the element). However, to integrate matrices corresponding to stabilization and diffusive terms, the procedure used in the first approach was followed (inverse matrix elements and determinants of the Jacobian matrix were evaluated at the center of the element). Both integration approaches were used for the implementation of TG, CBS and MTG schemes.
The laminar transonic flow around a NACA 0012 airfoil with and is analyzed. Results obtained with version 17.1 of ANSYS/Fluent are compared with those obtained with the MTG scheme for the same problem. All simulations performed for this case were obtained with TG, MTG and CBS schemes with the second integration approach, that is, a complete analytical integration (including inverse matrix and determinant of the Jacobian matrix in the expressions to be evaluated analytically) for the mass and convective matrices, while the inverse matrix and determinant of the Jacobian matrix were evaluated at the center of the elements for the stabilization and diffusive matrices.
This problem is described with dimensionless quantities. The mesh used to solve this problem with 386 560 hexahedral elements and 775 684 nodes is shown in Figure 3. The coordinates axis are and . Figures 4 and 5 show details of the finite element mesh. It is considered as being threedimensional mesh with a single element in the direction. The same mesh is used to run this problem with ANSYS/Fluent.
Considering that the mesh is threedimensional with a single element in the direction, the condition is imposed on the lateral boundaries to avoid mass flow through those boundaries, ensuring the twodimensional characteristic of the flow.
Figure 3. Side view of whole domain and the finite element mesh 
Figure 4. Side view of the mesh 
Figure 5. Details of the airfoil leading edge 
The dimensionless inlet boundary conditions are:

On the solid boundaries the nonslip condition is applied. Then, on the airfoil wall:

together with the stagnation temperature, which is specified using the specific internal energy, according to the following expression:

where is the recovery factor (an empirical factor introduced because in practice the energy recovery is not perfect [14]). The conditions at the output boundaries are specified using the vectors and .
Initial conditions are given by:

These initial conditions are defined for all nodes of the finite element mesh, excepting those where Dirichlet boundary conditions are prescribed.
Using the safety coefficient , the effective dimensionless time step is . Artificial viscosity is not applied for this case .
Figure 6 shows pressure contours obtained with the MTG scheme. Pressure contours obtained with CBS scheme are practically coincident with results obtained with MTG scheme and are not shown here. Figure 7 shows pressure contours obtained with ANSYS/Fluent 17.1 available at the GFC (Grupo de Fluidodinámica Computacional) of the UNLP (Universidad Nacional de La Plata). An excellent agreement was obtained between MTG and ANSYS/Fluent simulations. Figure 8 shows details of pressure contours at the leading edge of the airfoil obtained with the MTG scheme (CBS shows no appreciable difference). Figure 9 shows the same image, but obtained with the TG scheme, where pressure oscillations in zones of low air speeds can be observed. CBS and MTG schemes do not present these oscillations.
Figure 6. MTG pressure contours 
Figure 7. ANSYS / Fluent pressure contours 
Figure 8. MTG pressure contours 
Figure 9. TG pressure contours 
Figure 10 shows an image of the skin friction coefficient comparing the MTG scheme with values obtained with ANSYS/Fluent 17.1 showing a good agreement. Figure 11 shows the comparison of pressure coefficient obtained with ANSYS/Fluent 17.1 and the MTG scheme on the NACA 0012 airfoil extrados.
Here, the skin friction coefficient and the pressure coefficient were calculated using the following expressions:

where is the wall tangent stress and is the pressure, both acting on each point of the solid boundary surface.
Results obtained for these coefficients with MTG scheme applying an internal time step are very close to values obtained with ANSYS/Fluent 17.1.
Figure 10. Skin friction coefficient over the airfoil extrados surface 
Figure 11. Pressure coefficient over the airfoil extrados surface 
Figure 12 shows the convergence of the different schemes. The CBS method presents a clear improvement in density and pressure convergence when used with a local interior time step [16]. The use of an internal time step with the original TG method does not show comparative advantage with respect to the same method but taking a single time step. Notably, the performance of the MTG scheme with a local internal time step, is equivalent to that obtained with the CBS algorithm with a local internal time step.
Figure 12. Residues showing convergence of the different algorithms 
Additional simulations were carried out with the singlestep CBS method. Results are not shown here, but they do not show significant differences compared to those obtained with the CBS or MTG schemes.
Relative times for the different simulations are indicated in Table 1, assuming a reference time with a value equal to one for the fastest scheme (the singlestep CBS).
Algorithm  Relative wall clock time 
Singlestep CBS  1.0 
TG  1.01 
MTG  1.04 
CBS  1.37 
It should be noted that the cases shown in Table 1 correspond to simulations of the same problem, with the same mesh and boundary conditions, on the same computer (Dell Precision 7610, 16 Cores, 64 Gb ram), and using the same number of threads (16) in parallel. The times that are compared are wall clock times between two recordings of results with the same number of time steps. This problem was also solved using the first integration approach, that is, analytical integration with evaluation of the inverse matrices and determinants of the Jacobian matrices at the center of the elements for all the mass, convective, stabilization and diffusive matrices. There are not significant differences between both integration approaches.
It should also be pointed out that the codes for the different schemes were developed starting from the same original TG MPI parallel code, including only those modifications required by CBS, singlestep CBS and MTG schemes. This ensures that the code structure is basically the same, with smallest differences required by the different algorithms.
The turbulent compressible flow on a flat plate in supersonic regime with and is presented.
This problem is described at the NASA Langley Research Center Turbulence Modeling Resource on the web (https://turbmodels.larc.nasa.gov/), where results are shown for several turbulence models and they are compared with the semiempirical theoretical correlations of Van Driest.
The geometry of this problem is shown in Figure 13. Considering that the mesh is threedimensional with a single element in the direction, the condition is imposed on the lateral boundaries to avoid the mass flow through them, keeping the twodimensional characteristic of the flow.
Figure 13. 3D view and geometric characteristics 
Figure 14 shows the side view of the finite element mesh, which has 13 056 hexahedral elements and 26 578 nodes. The mesh is the same found in the NASA web repository for turbulence models (3D with nodes and points on the flat plate).
Figure 14. Side view of the mesh 
The mesh was refined near the wall such that the estimate of is close to 0.3.
The boundary conditions at the inlet are:

On the solid boundaries (over the plate) the nonslip condition is applied:

together with the stagnation temperature, which is specified using the specific internal energy, according to the following expression:

with . The conditions in the output boundary are specified using the vectors and . The initial conditions are given by:

These initial values are defined in all nodes of the mesh, except for the nodes where Dirichlet boundary conditions were prescribed.
Using the safety factor , the dimensionless time step is . The artificial viscosity coefficient used here is .
The pressure contours obtained with the TG, MTG and CBS schemes shown no appreciable differences. However, the researchers who developed the CBS scheme [6] point out that the singlestep CBS scheme is not suitable for high Mach numbers. The simulations carried out in this work with this scheme for high Mach numbers indicate that the convergence of the singlestep CBS presents some instabilities. Therefore, observations given in the reference [6] is confirmed here: singlestep CBS scheme is not recommended for supersonic flows.
In Figure 15 an image of the skin friction coefficient, comparing with values provided by the semiempirical correlation of Van Driest II (see the NASA Langley Research Center Turbulence Modeling Resource) is shown. Good agreement is observed. There are no appreciable differences between results obtained with different schemes (TG, CBS and MTG). All these schemes seem to be adequate in terms of quality of the results.
Figure 15. Skin friction coefficient over the plate 
Finally, Figure 16 shows the dimensionless velocity profile, corresponding to the position of the plate where . As can be observed, results obtained in this work are very close to values shown at the NASA repository for the velocity profile, corresponding to the semiempirical correlation of Van Driest I.
Figure 16. Velocity profile at 
The convergence between the TG, MTG and CBS schemes does not differ significantly when the residuals are plotted (although they are not shown here); however, there are differences in the relative computational times, as shown in Table 2.
Algorithm  Relative wall clock time 
Singlestep CBS  Not recommended 
TG  1.0 
MTG  1.03 
CBS  1.30 
In this way, it is observed that, for supersonic Mach numbers, the three algorithms are suitable from the point of view of convergence rates, but there are appreciable differences in terms of wall clock computational time, being the CBS scheme the most expensive.
The following case corresponds to the turbulent flow over an ogiveshaped obstacle in a cold hypersonic regime: , . This problem is described at the NASA Langley Research Center Turbulence Modeling Resource where results for several turbulence models are compared with experimental data. The purpose is to provide a validation case for turbulence models. Unlike verification, which seeks to establish whether a numerical model has been implemented correctly by comparison with other software or code, the validation of CFD results implies to compare with those obtained by experimental data, establishing the ability of a numerical model to reproduce the physics of the problem.
The experimental study conducted by Kussoy and Horstman [18] involved a cone cylinder with a taper angle of degrees. The leading edge is a 10degree cone that makes the transition through a circular arc to a constantarea cylinder with a radius of . The cone cylinder with a taper angle of degrees, designed to produce an oblique shock and an interaction zone between the boundary layer and the shock wave, is located downstream from the leading edge, as can be seen in Figure 17. As described in reference [19], the leading edge can be excluded from the CFD calculations, provided that the conditions of the input limits are slightly adjusted (the Mach number is changed from to . The wall (cylinder and cone) has a constant temperature of . The fluid domain and the boundary conditions are shown in the Figure 18.
Figure 17. Solid body (from [18]) 
Figure 18. Domain and boundaries 
Although this is an axisymmetric problem, its resolution is made here with a 3D mesh, taking a quarter of the geometry by symmetry considerations. In addition, this problem was simulated by the dimensionless codes developed for this work. In the two lateral boundaries is imposed that the normal velocity component is zero. Inlet boundary conditions are:

The nonslip condition is applied on solid boundaries. Then

The temperature at the wall, which has a known value , is specified using the dimensionless specific internal energy:

and the outlet boundary conditions are specified using and .
The initial conditions are given by:

These initial values are defined in all nodes of the finite element mesh, except for the nodes where Dirichlet boundary conditions were prescribed.
Figure 19 shows the side view of the finite element mesh, which has 1 600 000 hexahedral elements. The mesh is the same as found in the NASA web repository for turbulence models (2D with 2zones of 81x101 and 81x101 nodes) but rotated to form a threedimensional domain mesh corresponding to a quarter of the space surrounding the obstacle with 100 subdivisions in the quarter circumferential length. The total number of elements is . Figure 20 shows a detailed view of the mesh where the flare starts.
Figure 19. Finite element mesh 
Figure 20. Mask in the zone where the cone begins 
The mesh was refined near the wall such that the estimate of is close to 0.5.
Using the safety factor , the dimensionless time step is . The artificial viscosity coefficient used here is .
Figure 21 compares the experimental velocity profiles with those obtained in this work and by the NASA WindUS program. The position in which these data were obtained corresponds to the plane (the origin of coordinates is located where the flare starts). Since there are no appreciable differences between the results obtained with TG, CBS and MTG, only the curve corresponding to the MTG scheme is shown here.
Figure 21. Comparison of velocity profiles 
Figure 22 compares the temperature profiles for the same location. Both figures show proximity between the values obtained with the WindUS program [19] and the MTG code implemented in this work. In turn, results of both codes show a good agreement with the experimental data obtained in reference [18].
Figure 22. Comparison of temperature profiles 
Figure 23 shows a comparison of the pressure distribution along the wall of the ogive. The coordinate is measured along the surface of the solid with its origin located when the flare starts.
Figure 23. Pressure distribution on the wall 
Finally, Figure 24 shows a comparison of the heat transfer through the wall of the solid obstacle. The definition of wall heat transfer is given by: , where is the thermal conductivity of the fluid (but here it is calculated from the corresponding dimensionless quantities).
Figure 24. Heat transfer through the wall 
A good agreement between results obtained with the MTG scheme and "WindUs" codes. As observed from Figure 23, the SA model tends to slightly overpredict the pressure immediately downstream of the shock. Then the model starts to underpredict this quantity at approximately. The SA model also underpredicts the heat transfer though the solid wall after position . The differences between the results of both codes and the experimental data for positions above can be attributed to the adopted RANS/FANS turbulence model (SA), provided that all references show similar results [19,20].
Like the previous case, the convergence between the TG, MTG and CBS schemes does not differ significantly when the residuals are considered. However, there are differences in the relative wall clock computational times, as shown in Table 3.
Algorithm  Relative wall clock time 
Singlestep CBS  Not recommended 
TG  1.0 
TGM  1.03 
CBS  1.30 
It is observed that, for hypersonic flows, the three schemes (i.e. TG, MTG and CBS) are suitable from the point of view of convergence rates, but there are appreciable differences in terms of wall clock computational time, being the CBS scheme the most expensive.
The singlestep CBS scheme was used for this case showing some numerical oscillations in the flare starting zone, showing again that this algorithm is not recommended for high Mach numbers.
Comparing density residuals, CBS and MTG methods show better convergence rate than the TG scheme. But the best convergence rate is obtained using CBS and MTG with local Δt_{in}. The TG method exhibits an oscillatory behavior in transonic and subsonic flows. This behavior is not observed when CBS or MTG schemes are used.
Using a local internal time step with the original TG method, density convergence is not improved, and pressure oscillation remains. However, by making a change in the stabilization term of the continuity equation with a blending function and using a local internal time step, the MTG scheme presents a convergence rate comparable to that obtained with the CBS method.
Although the CBS method does not require an iterative process, as required by the TG and MTG schemes, split operations are performed to calculate first an auxiliary momentum variable, which is used to obtain the density to finally update the momentum values. Then the energy and the SA variable are calculated. These split operations imply in additional computational time cost. In addition, iterative terms of TG and MTG schemes converge rapidly with very few iterative steps and split operations are not necessary. These features make TG and MTG schemes more efficient than CBS method.
In this way, the MTG scheme is more efficient in terms of computational processing time than the CBS method and with a comparable convergence rate for all ranges of Mach numbers.
The alternative of the singlestep CBS method, although efficient and suitable for transonic regimes, is not recommended for high Mach numbers, such as supersonic or hypersonic regimes.
With the modifications proposed here for the TG scheme, a modified method was obtained (the MTG scheme), which is the best option when convergence properties and computational processing times are considered for a wide range of Mach numbers.
[5] Zienkiewicz O.C., Taylor R.L. The finite element method. Vol. 3: Fluid Dynamics. Fifth edition, Butterworth Heinemann, Oxford, 2000.
[14] White F. Viscous fluid flow. McGraw Hill Book Co, First edition, New York, 1974.
[15] Huebner K.H., Thornton E.A., Byrom T.G. The finite element method for engineers. Third Edition, John Wiley and Sons, Inc. New York, 1995.
[17] Argyris J., Doltsinis I.S., Friz H. Hermes space shuttle: exploration of reentry aerodynamics. Computer Methods in Applied Mechanics and Engineering, 73:151, 1989.
Published on 26/09/19
Accepted on 19/09/19
Submitted on 25/05/19
Volume 35, Issue 3, 2019
DOI: 10.23967/j.rimni.2019.09.006
Licence: CC BYNCSA license
Are you one of the authors of this document?