A modified Taylor Galerkin scheme for turbulent compressible flow simulation.

==Abstract==

As already demonstrated by different authors, the ''Taylor-Galerkin (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 non-viscous 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 ''Modified-Taylor-Galerkin'' 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 Spalart-Allmaras (S-A) turbulence model using a conservative version of the transport equation, as proposed by the authors of the original S-A 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, Taylor-Galerkin, CBS, Spalart-Allmaras.

==1. Introduction==

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 [1-2]. 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 non-viscous 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 [7-10]. 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/project.html 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 Spalart-Allmaras (S-A) 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 S-A 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 <math display="inline">\tilde{\nu }</math> in conservative form.

''TG'', ''CBS'' and ''MTG'' schemes are presented here with some cases of laminar and turbulent flows, comparing their performances.

==2. Governing equations with Spalart-Allmaras turbulence model==

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:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \boldsymbol{U}}{\partial t}+\frac{\partial {\boldsymbol{F}}_{i}}{\partial {x}_{i}}+</math><math>\frac{\partial {\boldsymbol{G}}_{i}}{\partial {x}_{i}}+\boldsymbol{Q}=\mathit{\boldsymbol{0}}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(1)
|}

where

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\boldsymbol{U}=\left\{ \begin{matrix}\rho \\ \rho {v}_{1} \\ \rho {v}_{2} \\ \rho {v}_{3} \\ \rho e \\ \rho \tilde{\nu }\end{matrix}\right\} ;{\boldsymbol{F}}_{i=}\left\{ \begin{matrix}\rho {v}_{i} \\ \rho {v}_{1}{v}_{i}+p{\delta }_{i1} \\ \rho {v}_{2}{v}_{i}+p{\delta }_{i2} \\ \rho {v}_{3}{v}_{i}+p{\delta }_{i3} \\ {v}_{i}\left( \rho e+p\right) \\ \rho \tilde{\nu }{v}_{i}\end{matrix}\right\} ;{\boldsymbol{G}}_{i}=</math><math>\left\{ \begin{matrix} 0 \\ -{\tau }_{i1} \\ -{\tau }_{i2} \\ -{\tau }_{i3} \\ -{\tau }_{ij}{v}_{j}-{q}_{i}-{\psi }_{i} \\ {\varphi }_{i}\end{matrix}\right\} ;\boldsymbol{Q}=</math><math>\left\{ \begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\P+D+{D}_{T}\end{matrix}\right\}</math> 
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(2)
|}

The dimensionless form has been obtained using a reference length, the freestream density and the freestream speed of sound. In Eq. (2):

:* <math display="inline">{v}_{i}\,</math> are velocity vector components.

:* <math display="inline">\rho</math>  is the specific mass.

:* <math display="inline">p</math>'' '' is the thermodynamic pressure.

:* <math display="inline">{\tau }_{ij}\,</math> are the viscous components of the stress tensor.

:* <math display="inline">e</math>''  ''is the specific total energy.

:* <math display="inline">\tilde{\nu }</math> is the transport equation variable of the S-A conservative turbulence model.

:* <math display="inline">{q}_{i}</math> are the conduction heat flux vector components.

:* <math display="inline">{\delta }_{ij}\,</math> is the Kronecker delta.

:* <math display="inline">{\varphi }_{i}</math> is the diffusive term of the S-A transport equation.

:* <math display="inline">{\psi }_{i}</math> is the diffusive term associated with turbulent kinetic energy ( <math display="inline">k</math>) in the energy equation.

:* <math display="inline">P</math> is the production term and ''D'' is the destruction term of the S-A transport equation.

:* <math display="inline">{D}_{T}\,</math> contains additional terms of the S-A transport equation.

:* <math display="inline">{x}_{i}\,</math> and ''t'' are the spatial and temporal coordinates, respectively.

:* <math display="inline">i,j</math> =1,2,3.

The viscous stress components of a Newtonian fluid are given by:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\tau }_{ij}=\left( \mu +{\mu }_{t}\right) \left[ \frac{\partial {v}_{i}}{\partial {x}_{j}}+\right. </math><math>\left. \frac{\partial {v}_{j}}{\partial {x}_{i}}\right] -\frac{2}{3}\left( \mu +\right. </math><math>\left. {\mu }_{t}\right) \frac{\partial {v}_{k}}{\partial {x}_{k}}{\delta }_{ij}-</math><math>\frac{2}{3}\rho k{\delta }_{ij}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(3)
|}


where <math display="inline">\mu</math>  and <math display="inline">{\mu }_{t}</math> are the molecular viscosity and the eddy viscosity, respectively, <math display="inline">k</math> is the turbulent kinetic energy and the subscript <math display="inline">k</math> runs from 1 to 3.

The components of the conduction heat flux vector are:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{q}_{i}=\left( K+{K}_{t}\right) \frac{\partial u}{\partial {x}_{i}}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(4)
|}

where <math display="inline">K</math> and <math display="inline">{K}_{t}</math> are the thermal conductivity and the turbulent thermal conductivity respectively and ''u'' is the specific internal energy.

The terms <math display="inline">{\psi }_{i}</math> in the energy equation are given by the following expression:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\psi }_{i}=\left( \mu +\frac{{\mu }_{t}}{{\sigma }_{k}}\right) \frac{\partial k}{\partial {x}_{i}}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(5)
|}


where <math display="inline">{\sigma }_{k}</math> is a constant which depends on the turbulence model. Since the S-A turbulence model does not include the calculation of turbulent kinetic energy, in this work, <math display="inline">{\psi }_{i}</math> in Eq. (2) and the third term in Eq. (3) are neglected.

The terms <math display="inline">{\varphi }_{i}</math> in the S-A transport equation are given by the following expression:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\varphi }_{i}=\frac{1}{\sigma }\left( \mu +\rho \tilde{\nu }\right) \frac{\partial \tilde{\nu }}{\partial {x}_{i}}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(6)
|}


The transport equation also includes the additional term <math display="inline">{D}_{T}</math> , expressed by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{D}_{T}=-\frac{1}{\sigma }\, \rho \, {cb}_{2}\frac{\partial \tilde{\nu }}{\partial {x}_{i}}\frac{\partial \tilde{\nu }}{\partial {x}_{i}}+</math><math>\, \frac{1}{\sigma }(\nu +\tilde{\nu })\, \frac{\partial \rho }{\partial {x}_{i}}\frac{\partial \tilde{\nu }}{\partial {x}_{i}}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(7)
|}

where <math display="inline">\sigma</math> '' ''and'' '' <math display="inline">{cb}_{2}</math> are constants of the S-A turbulence model.

Finally, the production and destruction terms, in the conservative transport equation of S-A are given by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>P=-{cb}_{1}\, S\, \rho \tilde{\nu };\, \, D={cw}_{1}\, {f}_{w}\, \rho \, {\left( \frac{\tilde{\nu }}{d}\right) }^{2}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(8)
|}

where the first term <math display="inline">P</math>'' ''represents the production and the second term'' D'' the destruction, being <math display="inline">{cb}_{1}</math>''<sub>,</sub> '' <math display="inline">S</math>'', '' <math display="inline">{cw}_{1}</math>'' ''and'' '' <math display="inline">{f}_{w}</math> constants and variables of the S-A turbulence model and <math display="inline">d</math> 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 trip-term 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 S-A turbulence model; they are given by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\sigma =\frac{2}{3};{\, cb}_{1}=0.1355;\, {\, cb}_{2}=0.622</math>

<math display="inline">S=\Omega +\, {f}_{v2}\, \frac{\tilde{\nu }}{{{\Kappa }^{2}d}^{2}}\, ;\, \Kappa =</math><math>0.41</math> <math display="inline"></math> <math>\boldsymbol{\Omega }=\, \, \frac{1}{2}\, rot\mathit{\boldsymbol{V;\, }}\Omega =\sqrt{\boldsymbol{\Omega }\cdot \boldsymbol{\Omega }}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(9)
|}

where <math display="inline">\Omega</math>  is the magnitude of the angular velocity vector, and

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{f}_{v2}=1-\, \, \frac{\chi }{1+\chi \, {f}_{v1}}\, ;\, \chi =\frac{\rho \tilde{\nu }}{\mu }\, ;\, {f}_{v1}=</math><math>\frac{{\chi }^{3}}{{\chi }^{3}+\, {{c}_{v1}}^{3}};\, {c}_{v1}=7.1</math>

<math>{f}_{w}=g{\left( \frac{1+{{c}_{w3}}^{6}}{{g}^{6}+{{c}_{w3}}^{6}}\right) }^{\frac{1}{6}}\, ;\, {c}_{w3}=2.0\, ;\, g=r+</math><math>{c}_{w2}\left( {r}^{6}-r\right) ;\, {c}_{w2}=0.3</math>

<math>r=\frac{\tilde{\nu }}{{{S\, \Kappa }^{2}d}^{2}}\, ;\, {c}_{w1}=\frac{{\, cb}_{1}}{{\Kappa }^{2}}+</math><math>\, \frac{1+{\, cb}_{2}}{\sigma }</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(10)
|}

The turbulent viscosity and turbulent thermal conductivity are defined as follows:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mu }_{t}=\rho \, \tilde{\nu }\, {f}_{v1};\, {K}_{t}={\mu }_{t}\frac{\gamma }{{Pr}_{t}};\, \, \gamma =</math><math>\frac{{c}_{p}}{{c}_{v}}=1.4</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(11)
|}


where <math display="inline">{Pr}_{t}</math> is the turbulent Prandtl number, which is adopted as <math display="inline">{Pr}_{t}=</math><math>0.9</math> and <math display="inline">{c}_{v}\,</math>  and <math display="inline">{c}_{p}\,</math>  are the specific heat coefficients at constant volume and at constant pressure, respectively.

The state equation for a perfect gas (air) is given by:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>p=\left( \gamma -1\right) \rho u;</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(12)
|}


The internal energy ''u'' and the temperature are related to the independent field variables by the following expression:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>u={c}_{v}\, T=e-\frac{1}{2}{v}_{i}{v}_{i}-k</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(13)
|}

where <math display="inline">T</math> is the temperature, being <math display="inline">k</math> neglected in the S-A 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:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\mu ={\mu }_{\infty }\frac{S+{u}_{\infty }}{S+u}{\left( \frac{u}{{u}_{\infty }}\right) }^{\frac{3}{2}};\quad \quad K=</math><math>{K}_{\infty }\, \frac{{S}_{k}+{u}_{\infty }}{{S}_{k}+u}{\left( \frac{u}{{u}_{\infty }}\right) }^{\frac{3}{2}}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(14)
|}


Values of <math display="inline">S</math> and <math display="inline">{S}_{k}</math>'' ''may be found in reference [14] for the Sutherland law based on temperature; they must be multiplied by <math display="inline">{c}_{v}</math> 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:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\begin{matrix}\begin{matrix}{v}_{i}=\overline{{v}_{i}}\quad on\quad {\Gamma }_{v}\\\rho =\overline{\rho }\quad \, \, on\quad {\Gamma }_{\rho }\\u=\overline{u}\quad \, \, on\quad {\Gamma }_{u}\end{matrix}\\\tilde{\nu }=\overline{\tilde{\nu }}\quad \quad on\quad {\Gamma }_{\tilde{\nu }}\end{matrix}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(15)
|}

where <math display="inline">\overline{{v}_{i}}</math> , <math display="inline">\overline{\rho }\, ,\, \overline{u}</math> and <math display="inline">\overline{\tilde{\nu }}</math> are prescribed values of the velocity vector components, the density (specific mass), the specific internal energy and the variable <math display="inline">\tilde{\nu }</math>, at the boundary surfaces <math display="inline">{\Gamma }_{v}</math> , <math display="inline">{\Gamma }_{\rho }</math> , <math display="inline">{\Gamma }_{u}</math>  and  ''Г'' <math display="inline">\tilde{v}\quad</math> respectively.

The natural boundary conditions (or Neumann boundary conditions) are defined by the following expressions:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{{\sigma }_{ij}\, n}_{j}=\hat{{t}_{i}}\quad on\quad {\Gamma }_{\sigma }</math>

<math display="inline">{{q}_{i}\, n}_{i}\, =\hat{q}\quad on\quad {\Gamma }_{q}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(16)
|}


In Eq. (16)

:* <math display="inline">{n}_{i}\,</math> are the cosines of the angles formed by the normal vector to the surface <math display="inline">{\Gamma }_{\sigma }</math> or <math display="inline">{\Gamma }_{q}</math> and the global reference axes  <math display="inline">{x}_{i}</math>.

:* <math display="inline">\hat{{t}_{i}}\,</math> are the tractions acting on <math display="inline">{\Gamma }_{\sigma }</math>.

:* <math display="inline">{\sigma }_{ij}\, =\, -p{\delta }_{ij}+\, {\tau }_{ij}</math>  are the components of the Cauchy stress tensor.

:* <math display="inline">\hat{q}\,</math>   is the heat flow acting on the perpendicular direction to <math display="inline">\, {\Gamma }_{q}</math>.

:* Effects of thermal radiation on boundary surfaces have not been considered.

For non-viscous flows the transport equation of the variable <math display="inline">\tilde{\nu }</math> and the diffusive terms <math display="inline">{\boldsymbol{G}}_{i}</math> 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.

==3. Finite element formulation==

In this section the well-known ''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.

===3.1 Taylor-Galerkin (TG) scheme===

To discretize Eq. (1) in time, the unknown vector <math display="inline">\boldsymbol{U}</math> is expanded in Taylor series omitting higher order terms. Then the space discretization is done using the classic Bubnov-Galerkin technique in the context of the finite element method. The Green-Gauss 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 [1-3].

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 <math display="inline">\, n+</math><math>1\,</math>  with the ''TG'' scheme depend of flow variables at the previous time step <math display="inline">n\,</math>  and at the current time step <math display="inline">\, n+</math><math>1</math>.  So, an iterative scheme to solve terms involving flow variables at the new time step is required [1,3].

===3.2 Characteristic Based Split (CBS) method===

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 S-A 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 Bubnov-Galerkin method. The Green-Gauss 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.

<span id='_Hlk7438137'></span>This method may be applied with implicit, semi-implicit 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].

===3.3 Single-Step CBS method===

An alternative of the ''CBS'' method consists in the elimination of the ''split'' in the momentum equation in order to obtain a one-step algorithm. Some authors [5,6] suggest that, although the single-step ''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 single-step ''CBS'' was implemented in this work. The common matrix expressions obtained with single-step ''CBS'' scheme can be found in reference [5].

==4. The modified Taylor-Galerkin scheme==

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 Taylor-Galerkin'' (''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'' <math display="inline">f</math> is proposed such that, the continuity equation discretized in time results:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\begin{matrix}\Delta {\left( \rho \right) }_{I+1}^{n+1}=\Delta t\left\{ -\frac{{\partial \left( {\rho v}_{i}\right) }^{n}}{{\partial x}_{i}}+\frac{\Delta t}{2}\left[ f\frac{\partial }{{\partial x}_{i}}\left( \frac{{\partial p}^{n}}{{\partial x}_{i}}\right) +\left( 1-f\right) \frac{\partial }{\partial {x}_{k}}\left( {{v}_{k}}^{n}\, \frac{{\partial \left( {\rho v}_{i}\right) }^{n}}{\partial {x}_{i}}\right) \right] \right\} +\\\quad \quad \quad \quad \quad \quad \quad \quad \quad +\frac{\Delta t}{2}\left\{ -\frac{\partial \Delta {\left( {\rho v}_{i}\right) }_{I}^{n+1}}{{\partial x}_{i}}+\frac{\Delta t}{2}\left[ f\frac{\partial }{{\partial x}_{i}}\left( \frac{{\partial \Delta p}_{I}^{n+1}}{{\partial x}_{i}}\right) +\left( 1-f\right) \frac{\partial }{\partial {x}_{k}}\left( {{v}_{k}}^{n}\, \frac{{\partial \left( \Delta {\rho v}_{i}\right) }_{I}^{n+1}}{\partial {x}_{i}}\right) \right] \right\} \end{matrix}\,</math> 
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(17)
|}

where

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>f=\left\{ \begin{matrix}\, \, 1\quad \quad \quad \quad \, if\quad \quad \, \quad M\leq 1\, \\-M+2\quad \quad \, \, if\quad \, 1<M<2\quad \, \, \\\, 0\quad \quad \quad \, \, \, \, \, if\quad \quad \quad \, M\geq 2\end{matrix}\, \right.</math> 
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(18)
|}


The ''blending function '' <math display="inline">f</math>'' ''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 Bubnov-Galerkin method. The Green-Gauss 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:

* Mass conservation equation:  

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
|<math>\begin{matrix}{\left. \boldsymbol{\Delta \rho }\right| }_{I+1}^{n+1}=f\left[ {\Delta t}_{ext}\, {{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -{\boldsymbol{B}}_{i}\, {\left. {\boldsymbol{f}}_{i}^{\rho }\right| }^{n}-\boldsymbol{P}{\boldsymbol{\, p}}^{n}-{\boldsymbol{r}}^{n}\right) +\frac{{\Delta t}_{ext}}{2}{{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -\, {\boldsymbol{B}}_{i}\, {\left. {\boldsymbol{\Delta f}}_{i}^{\rho }\right| }_{I}^{n+1}-\boldsymbol{P}\, {\boldsymbol{\Delta p}}_{I}^{n+1}\right) \right] +\\\quad \quad \quad +\left( 1-f\right) \left[ {\Delta t}_{ext}{{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -{\boldsymbol{H}}_{i}\, {\left. {\boldsymbol{f}}_{i}^{\rho }\right| }^{n}\right) +\frac{{\Delta t}_{ext}}{2}{{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -\, {\boldsymbol{H}}_{i}\, {\left. {\boldsymbol{\Delta f}}_{i}^{\rho }\right| }_{I}^{n+1}\right) \right] \end{matrix}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(19a)
|}

* Momentum conservation equations: 

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
|<math>\quad \begin{matrix}{\left. \boldsymbol{\Delta \rho }{\boldsymbol{v}}_{j}\right| }_{I+1}^{n+1}={\Delta t}_{ext}\, {{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -\, {\boldsymbol{H}}_{i}\, {\left. {\boldsymbol{f}}_{ij}^{\rho v}\right| }^{n}-{\boldsymbol{D}}_{ij}{{\boldsymbol{\, v}}_{i}}^{n}+{{\boldsymbol{t}}_{j}}^{n}\right) +\\\quad \quad \quad \quad \, +\frac{{\Delta t}_{ext}}{2}{{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -{\boldsymbol{H}}_{i}\, \boldsymbol{\Delta }{\left. {\boldsymbol{f}}_{ij}^{\rho v}\right| }_{I}^{n+1}-{\boldsymbol{D}}_{ij}{\left. {\boldsymbol{\, \Delta v}}_{i}\right| }_{I}^{n+1}\right) \end{matrix}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(19b)
|}

* Energy conservation equation:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
|<math>\, \begin{matrix}{\left. \boldsymbol{\Delta \rho e}\right| }_{I+1}^{n+1}={\Delta t}_{ext}{{\boldsymbol{\, M}}_{\boldsymbol{D}}}^{-1}\left( -{\boldsymbol{H}}_{i}\, {\left. {\boldsymbol{f}}_{i}^{\rho e}\right| }^{n}-{\boldsymbol{E}}_{i}\, {{\boldsymbol{v}}_{i}}^{n}-\boldsymbol{K}\, {\boldsymbol{u}}^{n}+{\boldsymbol{h}}^{n}\right) +\\\quad \quad \quad \quad \quad \quad +\frac{{\Delta t}_{ext}}{2}{{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -{\boldsymbol{H}}_{i}\, {\left. {\boldsymbol{\Delta f}}_{i}^{\rho e}\right| }_{I}^{n+1}-{\boldsymbol{E}}_{i}\, {\left. {\boldsymbol{\Delta v}}_{i}\right| }_{I}^{n+1}-\boldsymbol{K\, }{\left. \boldsymbol{\Delta u}\right| }_{I}^{n+1}\right) \end{matrix}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(19c)
|}

* Transport equation of  <math display="inline">\tilde{\nu }</math> (being <math display="inline">\tilde{\nu }</math> a variable used in the S-A model):

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
|<math>\begin{matrix}{\left. \boldsymbol{\Delta \rho }\tilde{\boldsymbol{\nu }}\right| }_{I+1}^{n+1}=\, {\Delta t}_{ext}{{\boldsymbol{\, M}}_{\boldsymbol{D}}}^{-1}\left( -{\boldsymbol{H}}_{i}\, {\left. {\boldsymbol{f}}_{i}^{\rho \tilde{\nu }}\right| }^{n}\boldsymbol{-M}\, {\boldsymbol{q}}^{n}-{\boldsymbol{A}}_{i}\, {\left. {\boldsymbol{q}}_{{\mathit{\boldsymbol{v}}}_{\mathit{\boldsymbol{i}}}}\right| }^{n}-{\boldsymbol{D}}_{\nu }\, {\tilde{\boldsymbol{\nu }}}^{n}\right) +\\\quad \quad \quad \quad \quad +\, \frac{{\Delta t}_{ext}}{2}{{\boldsymbol{M}}_{\boldsymbol{D}}}^{-1}\left( -{\boldsymbol{H}}_{i}\, {\left. {\boldsymbol{\Delta f}}_{i}^{\rho \tilde{\nu }}\right| }_{I}^{n+1}-\boldsymbol{M}\, {\boldsymbol{\Delta q}}_{I}^{n+1}-\, {\boldsymbol{D}}_{\nu }\, \boldsymbol{\Delta }{\tilde{\boldsymbol{\nu }}}_{I}^{n+1}\right) \end{matrix}\quad \, \,</math> 
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(19d)
|}


An iteration procedure is required for those right-hand terms evaluated on the current time step <math display="inline">n+</math><math>1</math>, where the subscript <math display="inline">I</math>'' ''indicates'' ''the previous iteration and <math display="inline">I+</math><math>1</math> corresponds to the current one. Element matrices are indicated below (with <math display="inline">i,j,k=</math><math>1,2,3</math>):

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
|style="text-align: center;"|<math>{\boldsymbol{B}}_{i}=\int_{{\Omega }_{E}}^{}{\boldsymbol{\phi }}^{T}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ;\quad {\boldsymbol{C}}_{i}=\frac{{\Delta t}_{in}}{2}\int_{{\Omega }_{E}}^{}\left( \boldsymbol{\phi }{{\boldsymbol{v}}_{k}}^{n}\right) \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{k}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ;\, \, \boldsymbol{P}=</math><math>\frac{{\Delta t}_{in}}{2}\int_{{\Omega }_{E}}^{}\frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{j}}\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}\, d\Omega \,</math>

<math>\, {\boldsymbol{H}}_{i}={\boldsymbol{B}}_{i}+{\boldsymbol{C}}_{i};\quad {\boldsymbol{A}}_{i}=</math><math>\frac{{\Delta t}_{in}}{2}{\boldsymbol{B}}_{i}\,</math>

<math>{\boldsymbol{D}}_{ij}=\left\{ \begin{matrix}\int_{{\Omega }_{E}}^{}\eta \left( 2+\frac{\lambda }{\eta }\right) \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\frac{\partial \boldsymbol{\phi }}{\partial {x}_{(i)}}d\Omega +\int_{{\Omega }_{E}}^{}\eta \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{k}}\frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}d\Omega ,\, \, if\, i=j\, and\left\{ \begin{matrix}i=1\rightarrow k=2,3\\i=2\rightarrow k=1,3\\i=3\rightarrow k=1,2\end{matrix}\right. \, \, \\\int_{{\Omega }_{E}}^{}\eta \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}\, d\Omega +\int_{{\Omega }_{E}}^{}\lambda \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{j}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ,\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad if\, \, i\, \not =j\end{matrix}\right.</math>

<math>\eta ={\mu +\mu }_{t};\, \lambda =-\frac{2}{3}\eta</math>

<math>{\boldsymbol{E}}_{i}=\int_{{\Omega }_{E}}^{}\left[ \eta \left( \boldsymbol{\phi }{{\boldsymbol{\, v}}_{i}}^{n}\right) \, \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{k}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}+\eta \left( \boldsymbol{\phi \, }{{\boldsymbol{v}}_{k}}^{n}\right) \, \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}\right] d\Omega +</math><math>\int_{{\Omega }_{E}}^{}\lambda \left( \boldsymbol{\phi }{{\boldsymbol{\, v}}_{k}}^{n}\right) \, \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{k}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega</math>

<math>\boldsymbol{K}=\int_{{\Omega }_{E}}^{}\left( {K+K}_{t}\right) \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ;\quad {\boldsymbol{D}}_{\nu }=</math><math>\int_{{\Omega }_{E}}^{}\frac{1}{\sigma }\left[ \mu +\left( \boldsymbol{\phi }{\left\{ \rho \nu \right\} }^{n}\right) \right] \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega</math>

<math>\boldsymbol{M}=\int_{{\Omega }_{E}}^{}{\boldsymbol{\phi }}^{T}\boldsymbol{\phi }d\Omega ;\, \, {\boldsymbol{M}}_{\boldsymbol{D}}=</math><math>\left\{ \begin{matrix}\frac{{\Omega }_{E}}{8}\quad \quad \quad \, \, \hbox{for diagonal terms.}\, \\\\0\quad \quad \quad \, \hbox{off-diagonal terms.}\end{matrix}\right.</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(20a)
|}


The element vectors are:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
|style="text-align: center;"|<math>{\boldsymbol{f}}_{i}^{\rho }=\left\{ \rho {v}_{i}\right\} ;\, {\boldsymbol{f}}_{ij}^{\rho v}=</math><math>\left\{ \rho {v}_{i}{v}_{j}+p{\delta }_{ij}\right\} ;\, {\boldsymbol{f}}_{i}^{\rho e}=</math><math>\left\{ \rho {ev}_{i}+p{v}_{i}\right\} ;\, {\boldsymbol{f}}_{i}^{\rho \tilde{\nu }}=</math><math>\left\{ \rho {v}_{i}\tilde{\nu }\right\}</math>

<math>\begin{matrix}\displaystyle{\boldsymbol{q}}^{n}=-{cb}_{1}\, {{S}_{E}}^{n}\, {\left\{ \rho \tilde{\nu }\right\} }^{n}+{cw}_{1}\, {{{f}_{w}}_{E}}^{n}\, {\left\{ \rho \, {\left( \frac{\tilde{\nu }}{d}\right) }^{2}\right\} }^{n}-\frac{{cb}_{2}}{\sigma }\, {\left( {\rho }_{E}\, {\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right) }^{n}\left\{ 1\right\} +\\ \displaystyle +\, \frac{1}{\sigma }{\left[ \left( {\nu }_{E}+{\tilde{\nu }}_{E}\right) {\left. \frac{\partial \rho }{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right] }^{n}\left\{ 1\right\} \end{matrix}</math>

<math>\begin{matrix}{\left. \displaystyle {\boldsymbol{q}}_{{\mathit{\boldsymbol{v}}}_{\mathit{\boldsymbol{i}}}}\right| }^{n}=-{cb}_{1}\, {{S}_{E}}^{n}\, {\left\{ \rho \tilde{\nu }{v}_{i}\right\} }^{n}+{cw}_{1}\, {{{f}_{w}}_{E}}^{n}\, {\left\{ \rho \, {\left( \frac{\tilde{\nu }}{d}\right) }^{2}{v}_{i}\right\} }^{n}-\frac{{cb}_{2}}{\sigma }\, {\left( {\rho }_{E}\, {\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right) }^{n}{{\boldsymbol{v}}_{i}}^{n}+\\ \displaystyle+\, \frac{1}{\sigma }{\left[ \left( {\nu }_{E}+{\tilde{\nu }}_{E}\right) {\left. \frac{\partial \rho }{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right] }^{n}{{\boldsymbol{v}}_{i}}^{n}\end{matrix}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(20b)
|}


where <math display="inline">\, {\boldsymbol{v}}_{i}=\left\{ {v}_{i}\right\}</math> , <math display="inline">\boldsymbol{p}=</math><math>\left\{ p\right\} \boldsymbol{,\, }\boldsymbol{u}=\left\{ u\right\}</math> ''' '''are <math display="inline">\,</math> 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:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
|  style="text-align: center;vertical-align: top;"|<math>\displaystyle {\boldsymbol{r}}^{n}=\frac{{\Delta t}_{in}}{2}\int_{{\Omega }_{E}}^{}{{\boldsymbol{\phi }}^{\ast }}^{T}\left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}{\boldsymbol{p}}^{n}\right) {n}_{i}\, d\Gamma</math>

<math>\displaystyle{{\boldsymbol{t}}_{j}}^{n}=\int_{{\Gamma }_{E}}^{}{{\boldsymbol{\phi }}^{\ast }}^{T}\, \left[ \eta \, \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, {{\boldsymbol{v}}_{j}}^{n}+\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}{{\boldsymbol{v}}_{i}}^{n}\right) +\, \lambda \, \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}\, {{\boldsymbol{v}}_{k}}^{n}\right) {\delta }_{ij}\right] {n}_{i}\, \, d\Gamma</math>

<math>\begin{matrix}\displaystyle{\boldsymbol{h}}^{n}={\displaystyle\int}_{{\Gamma }_{E}}^{}{{\boldsymbol{\phi }}^{\ast }}^{T}\, \left( \boldsymbol{\phi }{{\boldsymbol{\, v}}_{j}}^{n}\right) \left[ \eta \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, {{\boldsymbol{v}}_{j}}^{n}+\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}{{\boldsymbol{v}}_{i}}^{n}\right) +\lambda \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}\, {{\boldsymbol{v}}_{k}}^{n}\right) {\delta }_{ij}\right] {n}_{i}\, \, d\Gamma +\\ \displaystyle +\, {\displaystyle\int}_{{\Gamma }_{E}}^{}{{\boldsymbol{\phi }}^{\ast }}^{T}\, \left( {K+K}_{t}\right) \, \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}{\boldsymbol{u}}^{n}\right) \, {n}_{i}\, \, d\Gamma \end{matrix}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(20c)
|}


In these expressions <math display="inline">{\Omega }_{E}</math> and <math display="inline">{\Gamma }_{E}</math> are the element volume and its boundary surface respectively; expressions with the subscript <math display="inline">E</math> are element averages; <math display="inline">\boldsymbol{\phi }</math> = [''Ф<sub>1</sub>'', ''Ф<sub>2</sub>''<sub>, …</sub> ''Ф<sub>8</sub>''] is the row matrix containing the interpolation functions for each local node and <math display="inline">{\boldsymbol{\phi }}^{\ast }</math> is the row matrix containing the interpolation functions evaluated at the boundary surface, <math display="inline">\boldsymbol{M}</math> is the consistent mass matrix, while <math display="inline">{\boldsymbol{M}}_{\boldsymbol{D}}</math> is the lumped mass matrix and <math display="inline">\left\{ 1\right\}</math>  is a vector where all its components are equal to one. Terms that involve turbulent kinetic energy have been neglected.

The variables with superscript <math display="inline">n</math> are evaluated on previous time step and the variables with superscript <math display="inline">n+</math><math>1</math> are evaluated at current time step. Boundary vectors appearing in the S-A 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 <math display="inline">\rho</math> '', '' <math display="inline">\rho {v}_{j}</math>'' '' <math display="inline">\rho e</math>'' ''and'' '' <math display="inline">\rho \tilde{\nu }</math> 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 <math display="inline">\tilde{\nu }</math> 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 <math display="inline">{\boldsymbol{D}}_{ij}</math>, <math display="inline">{\boldsymbol{E}}_{i}</math> and <math display="inline">\boldsymbol{K}</math> as well as vectors <math display="inline">{{\boldsymbol{t}}_{j}}^{n}</math> and <math display="inline">{\boldsymbol{h}}^{n}</math> 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 <math display="inline">{\Delta t}_{in}</math>. For the ''TG'' scheme there are no appreciable differences in using a local internal time step or a single time step <math display="inline">{\Delta t}_{ext}</math> [16].

==5. Stability condition, shock capture and convergence==

Considering that these algorithms were implemented in their explicit form, they are conditionally stable and the stability condition of Courant-Friedrichs-Lewy (CFL) for each element, in dimensionless form, is used to ensure stability. A safety coefficient <math display="inline">\beta</math>  is adopted with values <math display="inline">0.2\leq \beta \leq \, 0.5</math> [3].

To capture the strong discontinuities and eliminate high frequency oscillations near shock waves, a well-known method is used to add artificial viscosity selectively in regions where high-pressure gradients are observed [17]. An artificial damping coefficient <math display="inline">{C}_{AD}</math> is used to control de amount of artificial viscosity.

The convergence of the iterative process is obtained when the following conditions are simultaneously satisfied:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{r}_{I+1}^{\rho }={\left[ \frac{\sum _{N}^{}{\left| {\rho }_{I+1}-{\rho }_{I}\right| }^{2}}{\sum _{N}^{}{{\rho }_{I}}^{2}}\right] }^{\frac{1}{2}}\leq To</math>

<math>{r}_{I+1}^{\rho v}={\left[ \frac{\sum _{N}^{}{\left| {\rho v}_{iI+1}-{\rho v}_{iI}\right| }^{2}}{\sum _{N}^{}{{\rho v}_{iI}}^{2}}\right] }^{\frac{1}{2}}\leq To</math>

<math>{r}_{I+1}^{\rho e}={\left[ \frac{\sum _{N}^{}{\left| {\rho e}_{I+1}-{\rho e}_{I}\right| }^{2}}{\sum _{N}^{}{{\rho e}_{I}}^{2}}\right] }^{\frac{1}{2}}\leq To</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(21)
|}

where <math display="inline">N</math> is an index of nodes ( <math display="inline">N=</math><math>1,2,3,\ldots ,NNODES</math>), <math display="inline">NNODES</math> is the total number of nodes of the finite element mesh and'' '' <math display="inline">To</math>'' ''is the tolerance (in this paper <math display="inline">To=</math><math>{10}^{-3}</math> was adopted). The solution process ends when the variable <math display="inline">t</math> (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:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{R}^{n+1}={\left[ \sum _{N}^{}{\left| {\rho }^{n+1}-{\rho }^{n}\right| }^{2}\right] }^{\frac{1}{2}}\leq TO{L}_{t}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(22)
|}

In this work <math display="inline">TO{L}_{t}={10}^{-5}</math>was adopted.

==6. Explicit integration of element matrices==

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:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\, \Phi }_{N}=\frac{1}{8}\left( 1+{\xi }_{1N}{\xi }_{1}\right) \left( 1+{\xi }_{2N}{\xi }_{2}\right) \left( 1+\right. </math><math>\left. {\xi }_{3N}{\xi }_{3}\right)</math> 
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(23)
|}


where'' '' <math display="inline">N=1,2,3,\ldots ,8</math> and <math display="inline">{\xi }_{jN}</math> are the natural coordinates of the node <math display="inline">N</math>. The derivatives of the interpolation functions with respect to the global coordinates are given by:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}={I}_{ij}^{-1}\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{i}}</math>
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(24)
|}


being <math display="inline">{I}_{ij}^{-1}</math> the components of the inverse of the Jacobian matrix <math display="inline">\mathit{\boldsymbol{\, }}\mathit{\boldsymbol{I}}</math>, which can be expressed as follows:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\mathit{\boldsymbol{I}}\left( {\xi }_{1},{\xi }_{2},{\xi }_{3}\right) =\left[ \begin{matrix}{I}_{11}\, {I}_{12}\, \, {I}_{13}\\{I}_{21}\, \, {I}_{22}\, \, {I}_{23}\\{I}_{31}\, \, {I}_{32}\, \, {I}_{33}\end{matrix}\right] =</math><math>\left[ \begin{matrix}\displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{1}}{\boldsymbol{x}}_{1}\, \, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{1}}{\boldsymbol{x}}_{2}\, \, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{1}}{\boldsymbol{x}}_{3}\\
\displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{2}}{\boldsymbol{x}}_{1}\, \, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{2}}{\boldsymbol{x}}_{2}\, \, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{2}}{\boldsymbol{x}}_{3}\\
\displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{3}}{\boldsymbol{x}}_{1}\, \, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{3}}{\boldsymbol{x}}_{2}\, \, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{3}}{\boldsymbol{x}}_{3}\end{matrix}\right]</math> 
|}
|  style="text-align: right;vertical-align: center;width: 5px;text-align: right;white-space: nowrap;"|(25)
|}

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding-bottom:10px;"| [[Image:Review_556075373152-image1.jpg|222px]] 
|style="padding-bottom:10px;"| [[Image:Review_556075373152-image2.jpg|228px]] 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 1'''. Physical space
|colspan="1" style="padding-bottom:10px;"| '''Figure 2'''. Computational space
|}

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image2.jpg|228px]] 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 2'''. Computational space
|}


Here <math display="inline">{\boldsymbol{x}}_{i}</math> is the vector containing global coordinates of the eight nodes of each element, <math display="inline">\, \left| \mathit{\boldsymbol{I}}\right|</math>  is the determinant of the Jacobian matrix and <math display="inline">d\Omega =</math><math>\left| \mathit{\boldsymbol{I}}\right| \, d{\xi }_{1}\, d{\xi }_{2\, }\, d{\xi }_{3\, }</math> 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 <math display="inline">{\xi }_{1}=</math><math>\, {\xi }_{2\, }=\, {\xi }_{3\, }=0,</math> 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 <math display="inline">\boldsymbol{\, }\boldsymbol{M}</math> (representing variables time derivatives) and matrix <math display="inline">{\boldsymbol{B}}_{i}</math> (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.

==7. Numerical simulations and comparative results==

===7.1 Transonic flow around a NACA 0012 airfoil===

The laminar transonic flow around a NACA 0012 airfoil with <math display="inline">M=</math><math>0.85\,</math>  and <math display="inline">Re=2000\,</math> 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 <math display="inline">x=</math><math>{x}_{1}</math> and <math display="inline">y={x}_{2}</math>. Figures 4 and  5 show details of the finite element mesh. It is considered as being three-dimensional mesh with a single element in the <math display="inline">z=</math><math>{x}_{3}</math> direction. The same mesh is used to run this problem with ANSYS/Fluent.

Considering that the mesh is three-dimensional with a single element in the <math display="inline">{x}_{3}</math> direction, the <math display="inline">\, {v}_{3}=</math><math>0.0</math> condition is imposed on the lateral boundaries to avoid mass flow through those boundaries, ensuring the two-dimensional characteristic of the flow.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image3.png|312px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 3'''. Side view of whole domain and the finite element mesh
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image4.jpeg|288px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 4'''. Side view of the mesh
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image5.jpeg|288px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 5'''. Details of the airfoil leading edge
|}


The dimensionless inlet boundary conditions are:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{v}_{1\infty }=0.85</math>; <math display="inline">{v}_{2\infty }=</math><math>{v}_{3\infty }=0.00</math>;  <math display="inline">{u}_{\infty }=1.7857</math>; <math display="inline">{\rho }_{\infty }=1.0</math>
|}
|}

On the solid boundaries the non-slip condition is applied. Then, on the airfoil wall:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>{v}_{1}={v}_{2}={v}_{3}=0.0</math>
|}
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>{u}_{stg}={u}_{\infty }\left( 1+r\frac{\gamma -1}{2}{{M}_{\infty }}^{2}\right) =</math><math>2.01535</math>
|}
|}

where <math display="inline">r=0.89</math> 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 <math display="inline">{\boldsymbol{t}}_{j}</math> and <math display="inline">\boldsymbol{q}</math>.

Initial conditions are given by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{{v}_{1}}^{0}=0.85</math>; <math display="inline">{{v}_{2}}^{0}=</math><math>{{v}_{3}}^{0}=0.00</math>; <math display="inline">{u}^{0}=1.7857</math>; <math display="inline">{\rho }^{0}=</math><math>1.0</math> y <math display="inline">{p}^{0}=0.71428</math>
|}
|}


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 <math display="inline">\beta =0.2</math>, the effective dimensionless time step is <math display="inline">\Delta t=</math><math>{0.5}^{\ast }{10}^{-4}</math>. Artificial viscosity is not applied for this case <math display="inline">({C}_{AD}=</math><math>0.0)</math>.

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.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image6-c.png|282px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 6'''. ''MTG'' pressure contours
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image7.jpeg|282px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 7'''. ANSYS / Fluent pressure contours
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image8-c.png|276px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 8'''. ''MTG'' pressure contours
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image9-c.png|282px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''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'' '' <math display="inline">{c}_{f}</math> and the pressure coefficient <math display="inline">{c}_{p}</math> were calculated using the following expressions:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>{c}_{f}=\frac{\tau }{\frac{1}{2}{\rho }_{\infty }{\left| {v}_{\infty }\right| }^{2}}\, ;\, {c}_{p}=</math><math>\frac{p-{p}_{\infty }}{\frac{1}{2}{\rho }_{\infty }{\left| {v}_{\infty }\right| }^{2}}</math>
|}
|}

where <math display="inline">\tau</math>  is the wall tangent stress and <math display="inline">p</math> 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.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Review_556075373152-image10-c.png|534px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 10'''. Skin friction coefficient over the airfoil extrados surface
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image11-c.png|600px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''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.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image12-c.png|414px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 12'''. Residues showing convergence of the different algorithms
|}


Additional simulations were carried out with the single-step ''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 single-step ''CBS'').

<div class="center" style="font-size: 75%;">'''Table 1'''.  Relative wall clock times
</div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;text-align:center;width:45%;" 
|-
|  style="border: 1pt solid black;vertical-align: top;"|'''Algorithm'''
|  style="border: 1pt solid black;vertical-align: top;"|'''Relative wall clock time'''
|-
|  style="border: 1pt solid black;vertical-align: top;"|Single-step CBS
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1.0
|-
|  style="border: 1pt solid black;vertical-align: top;"|TG
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1.01
|-
|  style="border: 1pt solid black;vertical-align: top;"|MTG
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1.04
|-
|  style="border: 1pt solid black;vertical-align: top;"|CBS
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|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'', single-step ''CBS'' and ''MTG ''schemes. This ensures that the code structure is basically the same, with smallest differences required by the different algorithms.

===7.2 Turbulent supersonic flow on flat plate===

The turbulent compressible flow on a flat plate in supersonic regime with <math display="inline">M=</math><math>2.0\,</math>  and <math display="inline">Re=15\ast {10}^{6}\,</math> is presented.

This problem is described at the NASA Langley Research Center Turbulence Modeling Resource on the web ([https://turbmodels.larc.nasa.gov/ https://turbmodels.larc.nasa.gov/]), where results are shown for several turbulence models and they are compared with the semi-empirical theoretical correlations of Van Driest.

The geometry of this problem is shown in Figure 13. Considering that the mesh is three-dimensional with a single element in the <math display="inline">{x}_{2}</math> direction, the condition <math display="inline">{v}_{2}=</math><math>0.0</math> is imposed on the lateral boundaries to avoid the mass flow through them, keeping the two-dimensional characteristic of the flow.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image13-c.png|372px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''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 (3-D with <math>2\times 137 \times 97</math> nodes and <math>2\times 113</math> points on the flat plate).

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image14-c.png|72px]]  [[Image:Review_556075373152-image15.png|414px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 14'''. Side view of the mesh
|}


The mesh was refined near the wall such that the estimate of <math display="inline">{y}^{+}</math> is close to 0.3.

The boundary conditions at the inlet are:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{v}_{1\infty }=2.0</math>; <math display="inline">{v}_{2\infty }=</math><math>{v}_{3\infty }=0.00</math>;  <math display="inline">{u}_{\infty }=1.7795</math>; <math display="inline">{\rho }_{\infty }=</math><math>1.0</math>; <math display="inline">{\tilde{\nu }}_{\infty }=6.667\ast {10}^{-7}</math>
|}
|}


On the solid boundaries (over the plate) the non-slip condition is applied:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{v}_{1}={v}_{2}={v}_{3}=</math><math>0.0</math>   and   <math display="inline">{\tilde{\nu }}_{\infty }=0.0</math>
|}
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>{u}_{stg}={u}_{\infty }\left( 1+r\frac{\gamma -1}{2}{{M}_{\infty }}^{2}\right) =</math><math>3.0465</math>
|}
|}

with <math display="inline">r=0.89</math>. The conditions in the output boundary are specified using the vectors <math display="inline">{\boldsymbol{t}}_{j}</math> and <math display="inline">\boldsymbol{q}</math>.

The initial conditions are given by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{{v}_{1}}^{0}=2.0</math>; <math display="inline">{{v}_{2}}^{0}=</math><math>{{v}_{3}}^{0}=0.00</math>; <math display="inline">{u}^{0}=1.7795</math>; <math display="inline">{\rho }^{0}=</math><math>1.0</math>; <math display="inline">{p}^{0}=0.7118</math> y <math display="inline">{\tilde{\nu }}^{0}=</math><math>6.667\ast {10}^{-7}</math>
|}
|}


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 <math display="inline">\beta =0.2</math>, the dimensionless time step is <math display="inline">\Delta t=</math><math>{1.0}^{\ast }{10}^{-7}</math>. The artificial viscosity coefficient used here is <math display="inline">{C}_{AD}=</math><math>0.3</math>.

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 single-step ''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 single-step ''CBS'' presents some instabilities. Therefore, observations given in the reference [6] is confirmed here: single-step ''CBS'' scheme is not recommended for supersonic flows.

<span id='_Hlk7529192'></span>In Figure 15 an image of the skin friction coefficient, comparing with values provided by the semi-empirical 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.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image16.png|294px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 15'''. Skin friction coefficient over the plate
|}


Finally, Figure 16 shows the dimensionless velocity profile, corresponding to the position of the plate where <math display="inline">{Re}_{\theta }=</math><math>10\, 000</math>. 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 semi-empirical correlation of Van Driest I.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image17.png|294px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 16'''. Velocity profile at <math display="inline">{Re}_{\theta }=</math><math>10\, 000</math>
|}


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.

<div class="center" style="font-size: 75%;">'''Table 2'''. Relative wall clock times
</div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;text-align:center;width:45%;" 
|-
|  style="border: 1pt solid black;vertical-align: top;"|'''Algorithm'''
|  style="border: 1pt solid black;vertical-align: top;"|'''Relative wall clock time'''
|-
|  style="border: 1pt solid black;vertical-align: top;"|Single-step CBS
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Not recommended
|-
|  style="border: 1pt solid black;vertical-align: top;"|TG
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1.0
|-
|  style="border: 1pt solid black;vertical-align: top;"|MTG
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1.03
|-
|  style="border: 1pt solid black;vertical-align: top;"|CBS
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|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.

===7.3 Axisymmetric Shock Wave Boundary Layer Interaction (ASWBLI). Turbulent cold hypersonic flow===

The following case corresponds to the turbulent flow over an ogive-shaped obstacle in a cold hypersonic regime: <math display="inline">M=</math><math>7.11</math>, <math display="inline">\frac{Re}{L}=57\, 060\, \frac{1}{cm}</math>. 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 <math display="inline">20</math> degrees. The leading edge is a 10-degree cone that makes the transition through a circular arc to a constant-area cylinder with a radius of <math display="inline">100\, mm</math>. The cone cylinder with a taper angle of <math display="inline">20</math> degrees, designed to produce an oblique shock and an interaction zone between the boundary layer and the shock wave, is located <math display="inline">1390\, mm</math> 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 <math display="inline">7.05</math> to <math display="inline">7.11</math>. The wall (cylinder and cone) has a constant temperature of <math display="inline">311\, K</math>. The fluid domain and the boundary conditions are shown in the Figure 18.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image18-c.png|264px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 17'''. Solid body (from [18])
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image19-c.png|330px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 18'''. Domain and boundaries
|}


Although this is an axisymmetric problem, its resolution is made here with a 3-D 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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{v}_{1\infty }=7.11</math>; <math display="inline">{v}_{2\infty }=</math><math>{v}_{3\infty }=0.00</math>;  <math display="inline">{u}_{\infty }=1.7795</math>; <math display="inline">{\rho }_{\infty }=</math><math>1.0</math>; <math display="inline">{\tilde{\nu }}_{\infty }=6.230284\ast {10}^{-4}</math>
|}
|}

The non-slip condition is applied on solid boundaries. Then

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{v}_{1}={v}_{2}={v}_{3}=</math><math>0.0</math>   and   <math display="inline">{\tilde{\nu }}_{\infty }=0.0</math>
|}
|}

The temperature at the wall, which has a known value <math display="inline">(311\, K)</math>, is specified using the dimensionless specific internal energy:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math>{u}_{w}=6.918</math>
|}
|}

and the outlet boundary conditions are specified using <math display="inline">{\boldsymbol{t}}_{j}</math> and <math display="inline">\boldsymbol{q}</math>.

The initial conditions are given by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |<math display="inline">{{v}_{1}}^{0}=7.11</math>; <math display="inline">{{v}_{2}}^{0}=</math><math>{{v}_{3}}^{0}=0.00</math>; <math display="inline">{u}^{0}=1.7795</math>; <math display="inline">{\rho }^{0}=</math><math>1.0</math>; <math display="inline">\, {p}^{0}=0.7118</math>  and  <math display="inline">{\tilde{\nu }}^{0}=</math><math>6.230284\ast {10}^{-4}</math>
|}
|}


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 (2-D with 2-zones of 81x101 and 81x101 nodes) but rotated to form a three-dimensional 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 <math>80 \times (2\times 100) \times 100</math>. Figure 20 shows a detailed view of the mesh where the flare starts.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image20.png|288px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 19'''. Finite element mesh
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Review_556075373152-image21.png|288px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 20'''. Mask in the zone where the cone begins
|}


The mesh was refined near the wall such that the estimate of <math>y+</math> is close to 0.5.

Using the safety factor <math display="inline">\beta =0.2</math>, the dimensionless time step is <math display="inline">\Delta t=</math><math>{1.0}^{\ast }{10}^{-5}</math>. The artificial viscosity coefficient used here is <math display="inline">{C}_{AD}=</math><math>0.3</math>.

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 <math display="inline">x=</math><math>-6</math> (the origin of coordinates <math display="inline">x</math> 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.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image22.png|294px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''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].

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Review_556075373152-image23.png|294px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''Figure 22'''. Comparison of temperature profiles
|}


Figure 23 shows a comparison of the pressure distribution along the wall of the ogive. The <math display="inline">s</math> coordinate is measured along the surface of the solid with its origin <math display="inline">s=</math><math>0\,</math> located when the flare starts.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|[[Image:Review_556075373152-image24.png|288px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''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 <math display="inline">{q}_{w}</math> is given by: <math display="inline">{q}_{w}=</math><math>-K\, {(dT/dy)}_{w}</math>, where <math display="inline">K</math> is the thermal conductivity of the fluid (but here it is calculated from the corresponding dimensionless quantities).

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|[[Image:Review_556075373152-image25.png|294px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;"| '''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 S-A model tends to slightly over-predict the pressure immediately downstream of the shock. Then the model starts to under-predict this quantity at <math display="inline">s=</math><math>7</math> approximately. The S-A model also under-predicts the heat transfer though the solid wall after position <math display="inline">s=</math><math>7</math>. The differences between the results of both codes and the experimental data for positions above <math display="inline">s=</math><math>7\,</math> can be attributed to the adopted RANS/FANS turbulence model (S-A), 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.

<div class="center" style="font-size: 75%;">'''Table 3'''. Relative computation times
</div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;text-align:center;width:45%;" 
|-
|  style="border: 1pt solid black;vertical-align: top;"|'''Algorithm'''
|  style="border: 1pt solid black;vertical-align: top;"|'''Relative wall clock time'''
|-
|  style="border: 1pt solid black;vertical-align: top;"|Single-step CBS
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Not recommended
|-
|  style="border: 1pt solid black;vertical-align: top;"|TG
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1.0
|-
|  style="border: 1pt solid black;vertical-align: top;"|TGM
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1.03
|-
|  style="border: 1pt solid black;vertical-align: top;"|CBS
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|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 single-step ''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.

==8. Conclusions==

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''<sub>in</sub>. 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 S-A 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 single-step ''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.

==References==

<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">


[1] [https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620200108 Donea J. A Taylor-Galerkin method for convective transport problems. International Journal for Numerical Methods in Engineering, 20:101-119, 1984.]

[2] [https://www.sciencedirect.com/science/article/pii/S004578259700114X?via%3Dihub Yoon K.T., Moon S.Y., Garcia S.A., Heard G.W., Chung T.J. Flowfield-dependent mixed explicit-implicit (FDMEI) methods for high and low speed and compressible and incompressible flows. Computer Methods in Applied Mechanics and Engineering, 151:75-104, 1998.]

[3] [http://www.dl.begellhouse.com/journals/1a37c0433d303e8b,7a0dd8b40eaafd93,529f21a1216620d2.html Burbridge H.P., Awruch A.M. A finite element Taylor-Galerkin scheme for 3-D numerical simulation of high compressible flow with analytical evaluation of element matrices. Hybrid Methods in Engineering, 2.I4:485-506, 2000.]

[4] [https://www.lume.ufrgs.br/bitstream/handle/10183/121789/000269871.pdf?sequence=1&isAllowed=y Burbridge H.P. O esquema explícito de Taylor-Galerkin na simulação numérica de escoamentos compressíveis tridimensionais utilizando elementos finitos hexaédricos de oito nós. Tese de mestrado. UFRGS, Porto Alegre, 1999.]

[5] Zienkiewicz O.C., Taylor R.L. The finite element method. Vol. 3: Fluid Dynamics. Fifth edition, Butterworth Heinemann, Oxford, 2000.

[6] [https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.1650200812 Zienkiewicz O.C., Codina R. A general algorithm for compressible and incompressible flow-Part I: The Split Charactestic-Based Scheme. International Journal for Numerical Methods in Fluids, 20:869-885, 1995.]

[7] [https://pdfs.semanticscholar.org/88bf/c1b51758a7e53cf7c38eb81603bb833f3c16.pdf Boonmarlert P., Phongthanapanich S., Dechaumpai P. Combined characteristic-based split algorithm and Mesh Adaptation Technique for high-speed compressible flow analysis. Indian Journal of Engineering & Material Sciences, 12:376-388, 2005.]

[8] [https://cronfa.swan.ac.uk/Record/cronfa42736/Download/0042736-02082018162517.pdf Liu C.B., The Characteristic Based Split (CBS) scheme for laminar and turbulent incompressible flow simulations. Doctoral Thesis, University of Wales, Swansea, 2005].

[9] [https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1700 Massarotti N., Arpino F., Lewis R.W., Nithiarasu P. Explicit and semi-implicit CBS procedures for incompressible viscous flows. International Journal for Numerical Methods in Engineering, 66:1618-1640, 2006.]

[10] [https://www.sciencedirect.com/science/article/pii/S1309104215302439 Scungio M., Arpino F, Cortellessa G., Buonanno G. Detached eddy simulation of turbulent flow in isolated street canyons of different aspect ratios. Atmospheric Pollution Research, 6.I2:351-364, 2015.]

[11] [https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.1650200813 Zienkiewicz O.C., Satya Sai B.V.K., Morgan K., Codina R., Vázquez M. A general algorithm for compressible and incompressible flow-Part II: Tests on the Explicit Form. International Journal for Numerical Methods in Fluids, 20:887-913, 1995.]

[12] [https://www.sciencedirect.com/science/article/pii/S0045782597000509 Christon M.A. A domain-decomposition message-passing approach to transient viscous incompressible flow using explicit time integration. Computer Methods in Applied Mechanics and Engineering, 148:329-352, 1997.]

[13] [http://iccfd.org/iccfd7/assets/pdf/papers/ICCFD7-1902_paper.pdf Allmaras S.R., Johnson F.T., Spalart P.R. Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model. Seventh International Conference on Computational Fluid Dynamics, ICCFD7-1902, 2012.]

[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.

[16] [https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291097-0207%2820000630%2948%3A6%3C875%3A%3AAID-NME907%3E3.0.CO%3B2-U Nithiarasu P.,  Zienwkiewicz O.C. On stabilization of the CBS algorithm: Internal and external time steps, International Journal for Numerical Methods in Engineering, 48:875-880, 2000.]

[17] Argyris J., Doltsinis I.S., Friz H. Hermes space shuttle: exploration of reentry aerodynamics. Computer Methods in Applied Mechanics and Engineering, 73:1-51, 1989.

[18] [https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19890010729.pdf Kussoy M.I., Horstman C.C. Documentation of two- and three-dimensional hypersonic shock wave boundary layer interaction flows. NASA TM 101075, January 1989]

[19] [https://arc.aiaa.org/doi/abs/10.2514/6.2015-0316 Georgiadis N.J., Rumsey Ch.L., Huang G.P. Revisiting turbulence model validation for high-mach number axisymmetric compression corner flows. 53rd AIAA Aerospace Sciences Meeting, AIAA 2015-0316, 2015.]

[20] [https://arc.aiaa.org/doi/abs/10.2514/6.2018-5195 Erb A., Hosder S. Investigation of turbulence model uncertainty for supersonic/hypersonic shock wave-boundary layer interaction predictions. 22nd AIAA International Space Planes and Hypersonics Systems and Technologies Conference. September 2018.]
</div>