(60 intermediate revisions by 4 users not shown) | |||
Line 1: | Line 1: | ||
− | + | ==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 | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | '''Keywords: '''Compressible | + | |
− | + | ||
− | + | ||
− | + | ||
− | + | ==1. Introduction== | |
− | + | ||
− | Different authors | + | 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. | 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 | + | 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. | 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 | + | 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 | + | 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. | ''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 | + | ==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: | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <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> | | <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;width: 5px;text-align: right;white-space: nowrap;"|(1) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | + | where | |
− | where | + | |
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
|- | |- | ||
− | + | | <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;width: 5px;text-align: right;white-space: nowrap;"|(2) | |
− | + | ||
− | | <math>\boldsymbol{U}=\left\{ | + | |
|} | |} | ||
− | |||
− | |||
− | |||
The dimensionless form has been obtained using a reference length, the freestream density and the freestream speed of sound. In Eq. (2): | The dimensionless form has been obtained using a reference length, the freestream density and the freestream speed of sound. In Eq. (2): | ||
Line 75: | Line 45: | ||
:* <math display="inline">\rho</math> is the specific mass. | :* <math display="inline">\rho</math> is the specific mass. | ||
− | :* <math display="inline">p</math> | + | :* <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">{\tau }_{ij}\,</math> are the viscous components of the stress tensor. | ||
− | :* <math display="inline">e</math> | + | :* <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">\tilde{\nu }</math> is the transport equation variable of the S-A conservative turbulence model. | ||
Line 89: | Line 59: | ||
:* <math display="inline">{\varphi }_{i}</math> is the diffusive term of the S-A transport equation. | :* <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">{\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">P</math> is the production term and ''D'' is the destruction term of the S-A transport equation. | ||
Line 102: | Line 72: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <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> | | <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;width: 5px;text-align: right;white-space: nowrap;"|(3) | |
− | | style="text-align: right | + | |
|} | |} | ||
Line 117: | Line 83: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
− | |||
| <math>{q}_{i}=\left( K+{K}_{t}\right) \frac{\partial u}{\partial {x}_{i}}</math> | | <math>{q}_{i}=\left( K+{K}_{t}\right) \frac{\partial u}{\partial {x}_{i}}</math> | ||
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(4) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | |||
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. | 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. | ||
Line 132: | Line 92: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <math>{\psi }_{i}=\left( \mu +\frac{{\mu }_{t}}{{\sigma }_{k}}\right) \frac{\partial k}{\partial {x}_{i}}</math> | | <math>{\psi }_{i}=\left( \mu +\frac{{\mu }_{t}}{{\sigma }_{k}}\right) \frac{\partial k}{\partial {x}_{i}}</math> | ||
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(5) | |
− | | style="text-align: right | + | |
|} | |} | ||
Line 147: | Line 103: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <math>{\varphi }_{i}=\frac{1}{\sigma }\left( \mu +\rho \tilde{\nu }\right) \frac{\partial \tilde{\nu }}{\partial {x}_{i}}</math> | | <math>{\varphi }_{i}=\frac{1}{\sigma }\left( \mu +\rho \tilde{\nu }\right) \frac{\partial \tilde{\nu }}{\partial {x}_{i}}</math> | ||
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(6) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | The transport equation also includes the additional term <math display="inline">{D}_{T}</math> , expressed by | + | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <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> | | <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;width: 5px;text-align: right;white-space: nowrap;"|(7) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | + | where <math display="inline">\sigma</math> and <math display="inline">{cb}_{2}</math> are constants of the S-A turbulence model. | |
− | where <math display="inline">\sigma</math> | + | |
Finally, the production and destruction terms, in the conservative transport equation of S-A are given by | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <math>P=-{cb}_{1}\, S\, \rho \tilde{\nu };\, \, D={cw}_{1}\, {f}_{w}\, \rho \, {\left( \frac{\tilde{\nu }}{d}\right) }^{2}</math> | | <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;width: 5px;text-align: right;white-space: nowrap;"|(8) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | + | 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>, <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. | |
− | where the first term <math display="inline">P</math> | + | |
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 | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <math>\sigma =\frac{2}{3};{\, cb}_{1}=0.1355;\, {\, cb}_{2}=0.622</math> | | <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> | <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;width: 5px;text-align: right;white-space: nowrap;"|(9) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | |||
where <math display="inline">\Omega</math> is the magnitude of the angular velocity vector, and | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <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}_{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> | ||
Line 214: | Line 148: | ||
<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> | <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;width: 5px;text-align: right;white-space: nowrap;"|(10) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | |||
The turbulent viscosity and turbulent thermal conductivity are defined as follows: | The turbulent viscosity and turbulent thermal conductivity are defined as follows: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <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> | | <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;width: 5px;text-align: right;white-space: nowrap;"|(11) | |
− | | style="text-align: right | + | |
|} | |} | ||
Line 234: | Line 162: | ||
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. | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <math>p=\left( \gamma -1\right) \rho u;</math> | | <math>p=\left( \gamma -1\right) \rho u;</math> | ||
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(12) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | |||
The internal energy ''u'' and the temperature are related to the independent field variables by the following expression: | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <math>u={c}_{v}\, T=e-\frac{1}{2}{v}_{i}{v}_{i}-k</math> | | <math>u={c}_{v}\, T=e-\frac{1}{2}{v}_{i}{v}_{i}-k</math> | ||
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(13) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | |||
where <math display="inline">T</math> is the temperature, being <math display="inline">k</math> neglected in the S-A model. | where <math display="inline">T</math> is the temperature, being <math display="inline">k</math> neglected in the S-A model. | ||
Line 265: | Line 183: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
− | |||
− | |||
− | |||
|- | |- | ||
| <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> | | <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;width: 5px;text-align: right;white-space: nowrap;"|(14) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | Values | + | 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: | 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: | ||
Line 281: | Line 195: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
|- | |- | ||
− | + | | <math>\begin{array}{ll}{v}_{i}=\overline{{v}_{i}}\quad & \hbox{on} \quad {\Gamma }_{v}\\ | |
− | + | \rho =\overline{\rho }\quad & \hbox{on} \quad {\Gamma }_{\rho }\\ | |
− | + | u=\overline{u}\quad & \hbox{on} \quad {\Gamma }_{u}\\ | |
− | | <math>\begin{ | + | \tilde{\nu }=\overline{\tilde{\nu }}\quad & \hbox{on}\quad {\Gamma }_{\tilde{\nu }}\end{array}</math> |
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(15) | |
− | | style="text-align: right | + | |
|} | |} | ||
− | |||
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. | 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. | ||
Line 296: | Line 208: | ||
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
|- | |- | ||
− | + | | <math>\begin{array}{ll}{{\sigma }_{ij}\, n}_{j}=\hat{{t}_{i}}\quad & \hbox{on} \quad {\Gamma }_{\sigma }\\[.4cm] | |
− | + | {{q}_{i}\, n}_{i}\, =\hat{q}\quad & \hbox{on} \quad {\Gamma }_{q}\end{array}</math> | |
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(16) | |
− | | <math>{{\sigma }_{ij}\, n}_{j}=\hat{{t}_{i}}\quad on\quad {\Gamma }_{\sigma } | + | |
− | + | ||
− | + | ||
|} | |} | ||
− | |||
− | |||
− | |||
In Eq. (16) | In Eq. (16) | ||
Line 321: | Line 227: | ||
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. | 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 | + | ==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 | + | 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 | + | 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 | + | 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 | + | 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. | 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. | ||
Line 339: | Line 245: | ||
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. | 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. | ||
− | + | 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 | + | Matrix expressions obtained with the ''CBS'' method may be found in references [5,6,8]. |
− | 3.3 | + | ===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 | + | 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 | + | ==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: | 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: | ||
Line 356: | Line 262: | ||
{| style="text-align: center; margin:auto;width: 100%;" | {| 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> | + | | <math>\begin{matrix}\Delta {\left( \rho \right) }_{I+1}^{n+1}=\Delta t\left\{ - \displaystyle \frac{{\partial \left( {\rho v}_{i}\right) }^{n}}{{\partial x}_{i}}+\displaystyle\frac{\Delta t}{2}\left[ f\frac{\partial }{{\partial x}_{i}}\left( \displaystyle\frac{{\partial p}^{n}}{{\partial x}_{i}}\right) +\left( 1-f\right) \displaystyle\frac{\partial }{\partial {x}_{k}}\left( {{v}_{k}}^{n}\, \displaystyle\frac{{\partial \left( {\rho v}_{i}\right) }^{n}}{\partial {x}_{i}}\right) \right] \right\} +\\\quad \quad \quad \quad \quad \quad \quad \quad \quad +\displaystyle\frac{\Delta t}{2}\left\{ -\displaystyle\frac{\partial \Delta {\left( {\rho v}_{i}\right) }_{I}^{n+1}}{{\partial x}_{i}}+\displaystyle\frac{\Delta t}{2}\left[ f\displaystyle\frac{\partial }{{\partial x}_{i}}\left( \displaystyle\frac{{\partial \Delta p}_{I}^{n+1}}{{\partial x}_{i}}\right) +\left( 1-f\right) \displaystyle\frac{\partial }{\partial {x}_{k}}\left( {{v}_{k}}^{n}\, \displaystyle\frac{{\partial \left( \Delta {\rho v}_{i}\right) }_{I}^{n+1}}{\partial {x}_{i}}\right) \right] \right\} \end{matrix}\,</math> |
|} | |} | ||
− | | style="text-align: right | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(17) |
|} | |} | ||
− | |||
where | where | ||
Line 369: | Line 274: | ||
{| style="text-align: center; margin:auto;width: 100%;" | {| 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> | + | | <math>f=\left\{ \begin{matrix}\, \, 1\quad \quad \quad \quad \, \hbox{if}\quad \quad \, \quad M\leq 1\, \\-M+2\quad \quad \, \, \hbox{if}\quad \, 1<M<2\quad \, \, \\\, 0\quad \quad \quad \, \, \, \, \, \hbox{if}\quad \quad \quad \, M\geq 2\end{matrix}\, \right.</math> |
|} | |} | ||
− | | style="text-align: right | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(18) |
|} | |} | ||
Line 381: | Line 286: | ||
The following matrix expressions are obtained for the element vectors of independent unknowns’ increments: | The following matrix expressions are obtained for the element vectors of independent unknowns’ increments: | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | * 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; | + | |<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) +\displaystyle\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) +\displaystyle\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;width: 5px;text-align: right;white-space: nowrap;"|(19a) | ||
|} | |} | ||
+ | * Momentum conservation equations: | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | {| 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; | + | |<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 \, +\displaystyle\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;width: 5px;text-align: right;white-space: nowrap;"|(19b) | ||
|} | |} | ||
+ | * Energy conservation equation: | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | {| 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; | + | |<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 +\displaystyle\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;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): | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | {| 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; | + | |<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 +\, \displaystyle\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;width: 5px;text-align: right;white-space: nowrap;"|(19d) | ||
|} | |} | ||
Line 419: | Line 337: | ||
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>): | 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>): | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | {| 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}}\, \displaystyle\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{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\ | + | <math>{\boldsymbol{D}}_{ij}=\left\{ \begin{matrix}\int_{{\Omega }_{E}}^{}\eta \left( 2+\displaystyle\frac{\lambda }{\eta }\right) \displaystyle\frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {x}_{(i)}}d\Omega +\int_{{\Omega }_{E}}^{}\eta \displaystyle\frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{k}}\displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}d\Omega ,\, \, if\, i=j\qquad \quad \hbox{and }\quad \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 \displaystyle\frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}\, d\Omega +\int_{{\Omega }_{E}}^{}\lambda \displaystyle\frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{j}}\, \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ,\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \hbox{if}\quad \, i\, \not =j\end{matrix}\right.</math> |
<math>\eta ={\mu +\mu }_{t};\, \lambda =-\frac{2}{3}\eta</math> | <math>\eta ={\mu +\mu }_{t};\, \lambda =-\frac{2}{3}\eta</math> | ||
Line 433: | Line 354: | ||
<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{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{ | + | <math>\boldsymbol{M}=\int_{{\Omega }_{E}}^{}{\boldsymbol{\phi }}^{T}\boldsymbol{\phi }d\Omega ;\, \, {\boldsymbol{M}}_{\boldsymbol{D}}=</math><math>\left\{ \begin{array}{cl}\frac{{\Omega }_{E}}{8}\quad \quad & \hbox{for diagonal terms.}\, \\\\0\quad \quad & \hbox{off-diagonal terms.}\end{array}\right.</math> |
− | + | |} | |
− | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(20a) | |
− | | style="text-align: right; | + | |
|} | |} | ||
Line 442: | Line 362: | ||
The element vectors are: | The element vectors are: | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | {| 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}{\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\} +\\+\, \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}\displaystyle{\boldsymbol{q}}^{n}=-{cb}_{1}\, {{S}_{E}}^{n}\, {\left\{ \rho \tilde{\nu }\right\} }^{n}+{cw}_{1}\, {{{f}_{w}}_{E}}^{n}\, {\left\{ \rho \, {\left( \displaystyle\frac{\tilde{\nu }}{d}\right) }^{2}\right\} }^{n}-\displaystyle\frac{{cb}_{2}}{\sigma }\, {\left( {\rho }_{E}\, {\left. \displaystyle\frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}{\left. \displaystyle\frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right) }^{n}\left\{ 1\right\} +\\ \displaystyle +\, \displaystyle\frac{1}{\sigma }{\left[ \left( {\nu }_{E}+{\tilde{\nu }}_{E}\right) {\left. \displaystyle\frac{\partial \rho }{\partial {x}_{i}}\right| }_{E}{\left. \displaystyle\frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right] }^{n}\left\{ 1\right\} \end{matrix}</math> |
− | <math>\begin{matrix}{\left. {\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}+\\+\, \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> | + | <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( \displaystyle\frac{\tilde{\nu }}{d}\right) }^{2}{v}_{i}\right\} }^{n}-\displaystyle\frac{{cb}_{2}}{\sigma }\, {\left( {\rho }_{E}\, {\left. \displaystyle\frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}{\left. \displaystyle\frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right) }^{n}{{\boldsymbol{v}}_{i}}^{n}+\\ \displaystyle+\, \displaystyle\frac{1}{\sigma }{\left[ \left( {\nu }_{E}+{\tilde{\nu }}_{E}\right) {\left. \displaystyle\frac{\partial \rho }{\partial {x}_{i}}\right| }_{E}{\left. \displaystyle\frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right] }^{n}{{\boldsymbol{v}}_{i}}^{n}\end{matrix}</math> |
− | | style="text-align: right; | + | |} |
+ | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(20b) | ||
|} | |} | ||
Line 455: | Line 379: | ||
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: | 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: | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" |
|- | |- | ||
− | | style="vertical-align: top;"|<math>{\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> | + | | |
+ | {| 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>{{\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>\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}{\boldsymbol{h}}^{n}=\ | + | <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( \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, {{\boldsymbol{v}}_{j}}^{n}+\displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}{{\boldsymbol{v}}_{i}}^{n}\right) +\lambda \left( \displaystyle\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( \displaystyle\frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}{\boldsymbol{u}}^{n}\right) \, {n}_{i}\, \, d\Gamma \end{matrix}</math> |
− | | style="text-align: right; | + | |} |
+ | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(20c) | ||
|} | |} | ||
Line 470: | Line 398: | ||
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. | 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 | + | 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. | 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. | ||
Line 476: | Line 404: | ||
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]. | 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 | + | ==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> | + | 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. | 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. | ||
Line 495: | Line 423: | ||
<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> | <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 | + | | style="text-align: right;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: | 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: | ||
Line 508: | Line 435: | ||
| <math>{R}^{n+1}={\left[ \sum _{N}^{}{\left| {\rho }^{n+1}-{\rho }^{n}\right| }^{2}\right] }^{\frac{1}{2}}\leq TO{L}_{t}</math> | | <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 | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(22) |
|} | |} | ||
− | |||
In this work <math display="inline">TO{L}_{t}={10}^{-5}</math>was adopted. | In this work <math display="inline">TO{L}_{t}={10}^{-5}</math>was adopted. | ||
− | =6. Explicit integration of element matrices | + | ==6. Explicit integration of element matrices== |
− | Eight nodes isoparametric hexahedral elements were used in this work. Figures 1 and | + | Eight nodes isoparametric hexahedral elements were used in this work. [[#img-1|Figures 1]] and [[#img-2|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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
Line 525: | Line 451: | ||
| <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> | | <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 | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(23) |
|} | |} | ||
− | where | + | 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;" | {| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" | ||
Line 538: | Line 464: | ||
| <math>\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}={I}_{ij}^{-1}\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{i}}</math> | | <math>\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}={I}_{ij}^{-1}\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{i}}</math> | ||
|} | |} | ||
− | | style="text-align: right | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(24) |
|} | |} | ||
Line 549: | Line 475: | ||
{| style="text-align: center; margin:auto;width: 100%;" | {| 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}\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{1}}{\boldsymbol{x}}_{1}\, \, \frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{1}}{\boldsymbol{x}}_{2}\, \, \frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{1}}{\boldsymbol{x}}_{3}\\\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{2}}{\boldsymbol{x}}_{1}\, \, \frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{2}}{\boldsymbol{x}}_{2}\, \, \frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{2}}{\boldsymbol{x}}_{3}\\\frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{3}}{\boldsymbol{x}}_{1}\, \, \frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{3}}{\boldsymbol{x}}_{2}\, \, \frac{\partial \boldsymbol{\phi }}{\partial {\xi }_{3}}{\boldsymbol{x}}_{3}\end{matrix}\right]</math> | + | | <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 | + | | style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(25) |
+ | |} | ||
+ | <div id='img-1'></div> | ||
+ | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;" | ||
+ | |- | ||
+ | |style="padding-top:10px;"| [[Image:Review_556075373152-image1.jpg|222px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" style="padding-bottom:10px;"| '''Figure 1'''. Physical space | ||
+ | |} | ||
+ | <div id='img-2'></div> | ||
+ | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;" | ||
+ | |- | ||
+ | |style="padding-top:10px;"| [[Image:Review_556075373152-image2.jpg|228px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" style="padding-bottom:10px;"| '''Figure 2'''. Computational space | ||
|} | |} | ||
Line 557: | Line 499: | ||
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. | 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 | + | 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. | 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. 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. | 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 | + | 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 [[#img-3|Figure 3]]. The coordinates axis are <math display="inline">x=</math><math>{x}_{1}</math> and <math display="inline">y={x}_{2}</math>. [[#img-4|Figures 4]] and [[#img-5|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. | 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. | ||
− | + | <div id='img-3'></div> | |
− | {| style=" | + | {| 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%;" | |
− | Figure 3 | + | | colspan="1" style="padding-bottom:10px;"| '''Figure 3'''. Side view of whole domain and the finite element mesh |
|} | |} | ||
− | + | <div id='img-4'></div> | |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;" | |
− | {| style=" | + | |
|- | |- | ||
− | | | + | |style="padding:10px;"| [[Image:Review_556075373152-image4.jpeg|288px]] |
− | + | |- style="text-align: center; font-size: 75%;" | |
− | Figure 4 | + | | colspan="1" style="padding-bottom:10px;"| '''Figure 4'''. Side view of the mesh |
− | | | + | |} |
− | + | <div id='img-5'></div> | |
− | Figure 5 | + | {| 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: | The dimensionless inlet boundary conditions are: | ||
− | {| style="width: 100%; | + | {| 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: | On the solid boundaries the non-slip condition is applied. Then, on the airfoil wall: | ||
− | {| style="width: 100%; | + | {| 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: | together with the stagnation temperature, which is specified using the specific internal energy, according to the following expression: | ||
− | {| style="width: 100%; | + | {| 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>. | 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>. | ||
Line 633: | Line 573: | ||
Initial conditions are given by: | Initial conditions are given by: | ||
− | {| style="width: 100%; | + | {| 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. | These initial conditions are defined for all nodes of the finite element mesh, excepting those where Dirichlet boundary conditions are prescribed. | ||
Line 644: | Line 586: | ||
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>. | 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. | + | [[#img-6|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. [[#img-7|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. [[#img-8|Figure 8]] shows details of pressure contours at the leading edge of the airfoil obtained with the ''MTG'' scheme (''CBS'' shows no appreciable difference). [[#img-9|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. |
− | + | <div id='img-6'></div> | |
− | {| style=" | + | {| 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%;" | |
− | Figure 6 | + | | colspan="1" style="padding-bottom:10px;"| '''Figure 6'''. ''MTG'' pressure contours |
− | | | + | |} |
− | + | <div id='img-7'></div> | |
− | + | {| 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%;" | |
− | Figure 8 | + | | colspan="1" style="padding-bottom:10px;"| '''Figure 7'''. ANSYS / Fluent pressure contours |
− | | | + | |} |
− | + | <div id='img-8'></div> | |
− | Figure 9 | + | {| 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 | ||
+ | |} | ||
+ | <div id='img-9'></div> | ||
+ | {| 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 | ||
|} | |} | ||
− | + | [[#img-10|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. [[#img-11|Figure 11]] shows the comparison of pressure coefficient obtained with ANSYS/Fluent 17.1 and the ''MTG'' scheme on the NACA 0012 airfoil extrados. | |
− | Figure 10 shows an image of the skin friction coefficient comparing the ''MTG ''scheme with values | + | |
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: | 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: | ||
− | {| style="width: 100%; | + | {| 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. | |
− | + | <div id='img-10'></div> | |
− | Results obtained for these coefficients with ''MTG'' scheme applying an internal time step are very close to values | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;" |
− | + | ||
− | {| style="width: | + | |
|- | |- | ||
− | | | + | |style="padding:10px;"| [[Image:Review_556075373152-image10-c.png|534px]] |
− | + | |- style="text-align: center; font-size: 75%;" | |
− | Figure 10 | + | | colspan="1" style="padding-bottom:10px;"| '''Figure 10'''. Skin friction coefficient over the airfoil extrados surface |
+ | |} | ||
+ | <div id='img-11'></div> | ||
+ | {| 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|534px]] |
− | + | |- style="text-align: center; font-size: 75%;" | |
− | Figure 11 | + | | colspan="1" style="padding-bottom:10px;"| '''Figure 11'''. Pressure coefficient over the airfoil extrados surface |
|} | |} | ||
− | + | [[#img-12|Figure 12]] shows the convergence of the different schemes. The ''CBS'' method presents a clear improvement in density and pressure convergence when used with a local interior time step [16]. The use of an internal time step with the original ''TG'' method does not show comparative advantage with respect to the same method but taking a single time step. Notably, the performance of the ''MTG'' scheme with a local internal time step, is equivalent to that obtained with the ''CBS'' algorithm with a local internal time step. | |
− | Figure 12 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. | + | <div id='img-12'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;" | |
− | {| style=" | + | |
|- | |- | ||
− | | | + | |style="padding:10px;"| [[Image:Review_556075373152-image12-c.png|414px]] |
− | + | |- style="text-align: center; font-size: 75%;" | |
− | Figure 12 | + | | 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. | 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''). | + | Relative times for the different simulations are indicated in [[#tab-1|Table 1]], assuming a reference time with a value equal to one for the fastest scheme (the single-step ''CBS''). |
− | + | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> | |
− | + | <span style="text-align: center; font-size: 75%;">'''Table 1'''. Relative wall clock times</span></div> | |
− | + | ||
− | | | + | <div id='tab-1'></div> |
− | |- | + | {| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" |
− | | style=" | + | |-style="text-align:center" |
− | | style=" | + | ! Algorithm !! Relative wall clock time |
− | |- | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|Single-step CBS |
− | | style=" | + | | style="text-align: center;vertical-align: top;"|1.0 |
− | |- | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|TG |
− | | style=" | + | | style="text-align: center;vertical-align: top;"|1.01 |
− | |- | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|MTG |
− | | style=" | + | | style="text-align: center;vertical-align: top;"|1.04 |
+ | |-style="text-align:center" | ||
+ | | style="vertical-align: top;"|CBS | ||
+ | | style="text-align: center;vertical-align: top;"|1.37 | ||
|} | |} | ||
− | + | It should be noted that the cases shown in [[#tab-1|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 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. | 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 | + | ===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. | 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. | ||
Line 737: | Line 692: | ||
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. | 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. | + | The geometry of this problem is shown in [[#img-13|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=" | + | {| 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%;" | |
− | Figure 13 | + | | colspan="1" style="padding-bottom:10px;"| '''Figure 13'''. 3D view and geometric characteristics |
|} | |} | ||
− | + | [[#img-14|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). | |
− | 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 | + | <div id='img-14'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;" | |
− | {| style=" | + | |
|- | |- | ||
− | | | + | |style="padding:10px;"| [[Image:Review_556075373152-image14-c.png|72px]] [[Image:Review_556075373152-image15.png|414px]] |
− | + | |- style="text-align: center; font-size: 75%;" | |
− | Figure 14 | + | | 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 mesh was refined near the wall such that the estimate of <math display="inline">{y}^{+}</math> is close to 0.3. | ||
Line 761: | Line 714: | ||
The boundary conditions at the inlet are: | The boundary conditions at the inlet are: | ||
− | {| style="width: 100%; | + | {| 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: | On the solid boundaries (over the plate) the non-slip condition is applied: | ||
− | {| style="width: 100%; | + | {| 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: | together with the stagnation temperature, which is specified using the specific internal energy, according to the following expression: | ||
− | {| style="width: 100%; | + | {| 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="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. | |
− | These initial values | + | |
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>. | 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 [6] is confirmed here: single-step ''CBS'' scheme is not recommended for supersonic flows. | + | 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. |
− | + | In [[#img-15|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. | |
+ | <div id='img-15'></div> | ||
+ | {| 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 | ||
+ | |} | ||
− | {| style=" | + | Finally, [[#img-16|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. |
+ | <div id='img-16'></div> | ||
+ | {| 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> | |
− | | | + | |
− | + | ||
− | Figure 16 | + | |
|} | |} | ||
+ | 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 [[#tab-2|Table 2]]. | ||
− | + | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> | |
+ | <span style="text-align: center; font-size: 75%;">'''Table 2'''. Relative wall clock times</span></div> | ||
− | + | <div id='tab-2'></div> | |
− | + | {| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" | |
− | {| style=" | + | |-style="text-align:center" |
− | + | ! Algorithm!! Relative wall clock time | |
− | + | |-style="text-align:center" | |
− | | | + | | style="vertical-align: top;"|Single-step CBS |
− | |- | + | | style="text-align: center;vertical-align: top;"|Not recommended |
− | | style=" | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|TG |
− | |- | + | | style="text-align: center;vertical-align: top;"|1.0 |
− | | style=" | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|MTG |
− | |- | + | | style="text-align: center;vertical-align: top;"|1.03 |
− | | style=" | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|CBS |
− | |- | + | | style="text-align: center;vertical-align: top;"|1.30 |
− | | style=" | + | |
− | | style=" | + | |
|} | |} | ||
− | |||
− | |||
− | |||
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. | 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 | + | ===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 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 | + | 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 [[#img-17|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 [[#img-18|Figure 18]]. |
+ | <div id='img-17'></div> | ||
+ | {| 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=" | + | <div id='img-18'></div> |
+ | {| 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 | |
− | | | + | |
− | + | ||
− | Figure 18 | + | |
|} | |} | ||
Line 863: | Line 830: | ||
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: | 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: | ||
− | {| style="width: 100%; | + | {| 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 | The non-slip condition is applied on solid boundaries. Then | ||
− | {| style="width: 100%; | + | {| 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: | The temperature at the wall, which has a known value <math display="inline">(311\, K)</math>, is specified using the dimensionless specific internal energy: | ||
− | {| style="width: 100%; | + | {| 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>. | and the outlet boundary conditions are specified using <math display="inline">{\boldsymbol{t}}_{j}</math> and <math display="inline">\boldsymbol{q}</math>. | ||
Line 892: | Line 865: | ||
The initial conditions are given by: | The initial conditions are given by: | ||
− | {| style="width: 100%; | + | {| 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 | + | 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 | + | [[#img-19|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 81<math>\times</math>101 and 81<math>\times</math>101 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>. [[#img-20|Figure 20]] shows a detailed view of the mesh where the flare starts. |
− | + | <div id='img-19'></div> | |
− | {| style=" | + | {| 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 | ||
+ | |} | ||
− | + | <div id='img-20'></div> | |
− | | | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;" |
− | + | |- | |
− | Figure 20 | + | |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 y+ is close to 0.5. | + | 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>. | 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. | + | [[#img-21|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. |
− | + | <div id='img-21'></div> | |
− | + | {| 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 | |
+ | |} | ||
− | |||
− | {| style=" | + | [[#img-22|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]. |
+ | <div id='img-22'></div> | ||
+ | {| 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 | ||
+ | |} | ||
− | |||
− | |||
− | + | [[#img-23|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. | |
+ | <div id='img-23'></div> | ||
+ | {| 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 | ||
|} | |} | ||
− | {| style=" | + | Finally, [[#img-24|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). |
+ | <div id='img-24'></div> | ||
+ | {| 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 | ||
+ | |} | ||
− | |||
− | |||
− | Figure | + | A good agreement between results obtained with the ''MTG'' scheme and "WindUs" codes. As observed from [[#img-23|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 [[#tab-3|Table 3]]. | ||
− | + | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;"> | |
+ | <span style="text-align: center; font-size: 75%;">'''Table 3'''. Relative computation times</span></div> | ||
− | {| style=" | + | <div id='tab-3'></div> |
− | + | {| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" | |
− | + | |-style="text-align:center" | |
− | | | + | ! Algorithm !! Relative wall clock time |
− | |- | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|Single-step CBS |
− | | style=" | + | | style="text-align: center;vertical-align: top;"|Not recommended |
− | |- | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|TG |
− | | style=" | + | | style="text-align: center;vertical-align: top;"|1.0 |
− | |- | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|TGM |
− | | style=" | + | | style="text-align: center;vertical-align: top;"|1.03 |
− | |- | + | |-style="text-align:center" |
− | | style=" | + | | style="vertical-align: top;"|CBS |
− | | style=" | + | | style="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. | 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. | ||
Line 978: | Line 969: | ||
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. | 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. | + | ==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. | 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. | ||
Line 992: | Line 983: | ||
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. | 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.] | [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 | + | [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 | + | [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 | + | [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 | + | [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 | + | [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 | + | [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 | + | [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 | + | [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. | + | [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 | + | [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.] | [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 | + | [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 | + | [14] White F. Viscous fluid flow. McGraw Hill Book Co, First edition, New York, 1974. |
− | [15] Huebner | + | [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 | + | [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 | + | [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 | + | [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 | + | [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. | + | [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> |
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
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/). 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 in conservative form.
TG, CBS and MTG schemes are presented here with some cases of laminar and turbulent flows, comparing their performances.
In an Eulerian description, the dimensionless system of partial differential equations governing the fluid dynamics problem in conservative form, and filtered with a Favre filter (FANS), may be written as follows:
(1) |
where
(2) |
The dimensionless form has been obtained using a reference length, the freestream density and the freestream speed of sound. In Eq. (2):
The viscous stress components of a Newtonian fluid are given by:
(3) |
where and are the molecular viscosity and the eddy viscosity, respectively, is the turbulent kinetic energy and the subscript runs from 1 to 3.
The components of the conduction heat flux vector are:
(4) |
where and are the thermal conductivity and the turbulent thermal conductivity respectively and u is the specific internal energy.
The terms in the energy equation are given by the following expression:
(5) |
where is a constant which depends on the turbulence model. Since the S-A turbulence model does not include the calculation of turbulent kinetic energy, in this work, in Eq. (2) and the third term in Eq. (3) are neglected.
The terms in the S-A transport equation are given by the following expression:
(6) |
The transport equation also includes the additional term , expressed by
(7) |
where and are constants of the S-A turbulence model.
Finally, the production and destruction terms, in the conservative transport equation of S-A are given by
(8) |
where the first term represents the production and the second term D the destruction, being , , and constants and variables of the S-A turbulence model and is the smallest distance from a field point (or node) to the nearest solid boundary (wall). The transport equation used here is expressed in its conservative form (Allmaras et al. [13]). Here the 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
|
(9) |
where is the magnitude of the angular velocity vector, and
|
(10) |
The turbulent viscosity and turbulent thermal conductivity are defined as follows:
(11) |
where is the turbulent Prandtl number, which is adopted as and and are the specific heat coefficients at constant volume and at constant pressure, respectively.
The state equation for a perfect gas (air) is given by:
(12) |
The internal energy u and the temperature are related to the independent field variables by the following expression:
(13) |
where is the temperature, being neglected in the 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:
(14) |
Values of and may be found in reference [14] for the Sutherland law based on temperature; they must be multiplied by and divided by the square of the freestream sound speed to obtain the dimensionless form, which was used in this work.
Initial and boundary conditions must be applied to Eq. (1) in order to obtain a mathematical model correctly defined. The forced boundary conditions (or Dirichlet boundary conditions) are given by:
(15) |
where , and are prescribed values of the velocity vector components, the density (specific mass), the specific internal energy and the variable , at the boundary surfaces , , and Г respectively.
The natural boundary conditions (or Neumann boundary conditions) are defined by the following expressions:
(16) |
In Eq. (16)
For non-viscous flows the transport equation of the variable and the diffusive terms in Eq. (1) should be omitted. In this case, only normal velocity components to the solid boundaries are prescribed with null values and only the normal components of the surface tractions on the outlet boundaries are considered.
In this section the 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.
To discretize Eq. (1) in time, the unknown vector is expanded in Taylor series omitting higher order terms. Then the space discretization is done using the classic 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 with the TG scheme depend of flow variables at the previous time step and at the current time step . So, an iterative scheme to solve terms involving flow variables at the new time step is required [1,3].
The CBS method consists of four steps [6]. In the first step, intermediate values of the conservative flow variables of the momentum equations are calculated, omitting the terms containing pressure gradients. In the second step, the continuity equation is solved to obtain the density increment. In the third step, the flow variables of the momentum equations are updated. In the fourth step, the energy equation and the 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.
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].
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].
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 is proposed such that, the continuity equation discretized in time results:
|
(17) |
where
|
(18) |
The blending function depends of the Mach number, which is evaluated locally. Momentum and energy equations as well as their iterative terms are not modified, remaining as in the original TG scheme.
Spatial discretization is performed using the classical 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:
|
(19a) |
|
(19b) |
|
(19c) |
|
(19d) |
An iteration procedure is required for those right-hand terms evaluated on the current time step , where the subscript indicates the previous iteration and corresponds to the current one. Element matrices are indicated below (with ):
|
(20a) |
The element vectors are:
|
(20b) |
where , are the velocity components, pressure and internal energy element vectors, containing their corresponding eight nodal values. The quantities between braces are element vectors built with the nodal values of the variables inside the braces. The boundary vectors are:
|
(20c) |
In these expressions and are the element volume and its boundary surface respectively; expressions with the subscript are element averages; = [Ф1, Ф2, … Ф8] is the row matrix containing the interpolation functions for each local node and is the row matrix containing the interpolation functions evaluated at the boundary surface, is the consistent mass matrix, while is the lumped mass matrix and is a vector where all its components are equal to one. Terms that involve turbulent kinetic energy have been neglected.
The variables with superscript are evaluated on previous time step and the variables with superscript are evaluated at current time step. Boundary vectors appearing in the 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 , and can be calculated at each time step using an iterative scheme. Thus, the nodal values of the components of velocity, the specific total energy and the variable can be obtained immediately. Then the nodal values of the specific internal energy are calculated and, using the state equation, the nodal pressure values are also obtained.
When viscous terms are not considered, terms containing matrices , and as well as vectors and and the transport equation corresponding to the turbulence model are omitted.
Stabilization matrices for TG, CBS and MTG schemes are evaluated using a local (internal) time step . For the TG scheme there are no appreciable differences in using a local internal time step or a single time step [16].
Considering that these algorithms were implemented in their explicit form, they are conditionally stable and the stability condition of Courant-Friedrichs-Lewy (CFL) for each element, in dimensionless form, is used to ensure stability. A safety coefficient is adopted with values [3].
To capture the strong discontinuities and eliminate high frequency oscillations near shock waves, a well-known method is used to add artificial viscosity selectively in regions where high-pressure gradients are observed [17]. An artificial damping coefficient is used to control de amount of artificial viscosity.
The convergence of the iterative process is obtained when the following conditions are simultaneously satisfied:
|
(21) |
where is an index of nodes ( ), is the total number of nodes of the finite element mesh and is the tolerance (in this paper was adopted). The solution process ends when the variable (time) reaches a value previously established by the user, or when the steady state is reached. This last condition is considered fulfilled when the following expression is satisfied:
|
(22) |
In this work was adopted.
Eight nodes isoparametric hexahedral elements were used in this work. Figures 1 and 2 show the element in the physical and computational spaces, respectively. The interpolation functions of this element are given by the following expressions:
|
(23) |
where and are the natural coordinates of the node . The derivatives of the interpolation functions with respect to the global coordinates are given by:
|
(24) |
being the components of the inverse of the Jacobian matrix , which can be expressed as follows:
|
(25) |
![]() |
Figure 1. Physical space |
![]() |
Figure 2. Computational space |
Here is the vector containing global coordinates of the eight nodes of each element, is the determinant of the Jacobian matrix and is the differential of the element volume.
By substituting Eqs. (23), (24) and (25) on the element matrices and vectors defined by the TG, CBS and MTG schemes, integrals in the computational domain are obtained. These integrals are commonly calculated by means of a numerical integration scheme to obtain components of the element matrices that will be later assembled to get the global system of equations. However, in this work, matrices were calculated by means of analytical integration using the Maxima symbolic resolution software. This implies an improvement in computational processing time efficiency.
In a first approach the inverse matrix elements and determinants of the Jacobian matrix were evaluated at the center of the element, where and these values were used to integrate matrices corresponding to variables time derivatives, convective, diffusive and stabilization terms [3,12]. Special attention must be given to elements distortion and an hourglass control technique to avoid spurious modes on diffusive terms are also necessary [12].
Alternatively, in a second approach, the mass matrix (representing variables time derivatives) and matrix (representing convective terms) were integrated analytically, including the inverse matrix elements and determinant of the Jacobian matrix in the analytical expressions (without evaluating these terms at the center of the element). However, to integrate matrices corresponding to stabilization and diffusive terms, the procedure used in the first approach was followed (inverse matrix elements and determinants of the Jacobian matrix were evaluated at the center of the element). Both integration approaches were used for the implementation of TG, CBS and MTG schemes.
The laminar transonic flow around a NACA 0012 airfoil with and is analyzed. Results obtained with version 17.1 of ANSYS/Fluent are compared with those obtained with the MTG scheme for the same problem. All simulations performed for this case were obtained with TG, MTG and CBS schemes with the second integration approach, that is, a complete analytical integration (including inverse matrix and determinant of the Jacobian matrix in the expressions to be evaluated analytically) for the mass and convective matrices, while the inverse matrix and determinant of the Jacobian matrix were evaluated at the center of the elements for the stabilization and diffusive matrices.
This problem is described with dimensionless quantities. The mesh used to solve this problem with 386 560 hexahedral elements and 775 684 nodes is shown in Figure 3. The coordinates axis are and . Figures 4 and 5 show details of the finite element mesh. It is considered as being three-dimensional mesh with a single element in the direction. The same mesh is used to run this problem with ANSYS/Fluent.
Considering that the mesh is three-dimensional with a single element in the direction, the condition is imposed on the lateral boundaries to avoid mass flow through those boundaries, ensuring the two-dimensional characteristic of the flow.
![]() |
Figure 3. Side view of whole domain and the finite element mesh |
![]() |
Figure 4. Side view of the mesh |
![]() |
Figure 5. Details of the airfoil leading edge |
The dimensionless inlet boundary conditions are:
|
On the solid boundaries the non-slip condition is applied. Then, on the airfoil wall:
|
together with the stagnation temperature, which is specified using the specific internal energy, according to the following expression:
|
where is the recovery factor (an empirical factor introduced because in practice the energy recovery is not perfect [14]). The conditions at the output boundaries are specified using the vectors and .
Initial conditions are given by:
|
These initial conditions are defined for all nodes of the finite element mesh, excepting those where Dirichlet boundary conditions are prescribed.
Using the safety coefficient , the effective dimensionless time step is . Artificial viscosity is not applied for this case .
Figure 6 shows pressure contours obtained with the MTG scheme. Pressure contours obtained with CBS scheme are practically coincident with results obtained with MTG scheme and are not shown here. Figure 7 shows pressure contours obtained with ANSYS/Fluent 17.1 available at the GFC (Grupo de Fluidodinámica Computacional) of the UNLP (Universidad Nacional de La Plata). An excellent agreement was obtained between MTG and ANSYS/Fluent simulations. Figure 8 shows details of pressure contours at the leading edge of the airfoil obtained with the MTG scheme (CBS shows no appreciable difference). Figure 9 shows the same image, but obtained with the TG scheme, where pressure oscillations in zones of low air speeds can be observed. CBS and MTG schemes do not present these oscillations.
![]() |
Figure 6. MTG pressure contours |
![]() |
Figure 7. ANSYS / Fluent pressure contours |
![]() |
Figure 8. MTG pressure contours |
![]() |
Figure 9. TG pressure contours |
Figure 10 shows an image of the skin friction coefficient comparing the MTG scheme with values obtained with ANSYS/Fluent 17.1 showing a good agreement. Figure 11 shows the comparison of pressure coefficient obtained with ANSYS/Fluent 17.1 and the MTG scheme on the NACA 0012 airfoil extrados.
Here, the skin friction coefficient and the pressure coefficient were calculated using the following expressions:
|
where is the wall tangent stress and is the pressure, both acting on each point of the solid boundary surface.
Results obtained for these coefficients with MTG scheme applying an internal time step are very close to values obtained with ANSYS/Fluent 17.1.
![]() |
Figure 10. Skin friction coefficient over the airfoil extrados surface |
![]() |
Figure 11. Pressure coefficient over the airfoil extrados surface |
Figure 12 shows the convergence of the different schemes. The CBS method presents a clear improvement in density and pressure convergence when used with a local interior time step [16]. The use of an internal time step with the original TG method does not show comparative advantage with respect to the same method but taking a single time step. Notably, the performance of the MTG scheme with a local internal time step, is equivalent to that obtained with the CBS algorithm with a local internal time step.
![]() |
Figure 12. Residues showing convergence of the different algorithms |
Additional simulations were carried out with the 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).
Algorithm | Relative wall clock time |
---|---|
Single-step CBS | 1.0 |
TG | 1.01 |
MTG | 1.04 |
CBS | 1.37 |
It should be noted that the cases shown in Table 1 correspond to simulations of the same problem, with the same mesh and boundary conditions, on the same computer (Dell Precision 7610, 16 Cores, 64 Gb ram), and using the same number of threads (16) in parallel. The times that are compared are wall clock times between two recordings of results with the same number of time steps. This problem was also solved using the first integration approach, that is, analytical integration with evaluation of the inverse matrices and determinants of the Jacobian matrices at the center of the elements for all the mass, convective, stabilization and diffusive matrices. There are not significant differences between both integration approaches.
It should also be pointed out that the codes for the different schemes were developed starting from the same original TG MPI parallel code, including only those modifications required by CBS, single-step CBS and MTG schemes. This ensures that the code structure is basically the same, with smallest differences required by the different algorithms.
The turbulent compressible flow on a flat plate in supersonic regime with and is presented.
This problem is described at the NASA Langley Research Center Turbulence Modeling Resource on the web (https://turbmodels.larc.nasa.gov/), where results are shown for several turbulence models and they are compared with the 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 direction, the condition is imposed on the lateral boundaries to avoid the mass flow through them, keeping the two-dimensional characteristic of the flow.
![]() |
Figure 13. 3D view and geometric characteristics |
Figure 14 shows the side view of the finite element mesh, which has 13 056 hexahedral elements and 26 578 nodes. The mesh is the same found in the NASA web repository for turbulence models (3-D with nodes and points on the flat plate).
![]() ![]() |
Figure 14. Side view of the mesh |
The mesh was refined near the wall such that the estimate of is close to 0.3.
The boundary conditions at the inlet are:
|
On the solid boundaries (over the plate) the non-slip condition is applied:
|
together with the stagnation temperature, which is specified using the specific internal energy, according to the following expression:
|
with . The conditions in the output boundary are specified using the vectors and . The initial conditions are given by:
|
These initial values are defined in all nodes of the mesh, except for the nodes where Dirichlet boundary conditions were prescribed.
Using the safety factor , the dimensionless time step is . The artificial viscosity coefficient used here is .
The pressure contours obtained with the TG, MTG and CBS schemes shown no appreciable differences. However, the researchers who developed the CBS scheme [6] point out that the 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.
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.
![]() |
Figure 15. Skin friction coefficient over the plate |
Finally, Figure 16 shows the dimensionless velocity profile, corresponding to the position of the plate where . As can be observed, results obtained in this work are very close to values shown at the NASA repository for the velocity profile, corresponding to the semi-empirical correlation of Van Driest I.
![]() |
Figure 16. Velocity profile at |
The convergence between the TG, MTG and CBS schemes does not differ significantly when the residuals are plotted (although they are not shown here); however, there are differences in the relative computational times, as shown in Table 2.
Algorithm | Relative wall clock time |
---|---|
Single-step CBS | Not recommended |
TG | 1.0 |
MTG | 1.03 |
CBS | 1.30 |
In this way, it is observed that, for supersonic Mach numbers, the three algorithms are suitable from the point of view of convergence rates, but there are appreciable differences in terms of wall clock computational time, being the CBS scheme the most expensive.
The following case corresponds to the turbulent flow over an ogive-shaped obstacle in a cold hypersonic regime: , . This problem is described at the NASA Langley Research Center Turbulence Modeling Resource where results for several turbulence models are compared with experimental data. The purpose is to provide a validation case for turbulence models. Unlike verification, which seeks to establish whether a numerical model has been implemented correctly by comparison with other software or code, the validation of CFD results implies to compare with those obtained by experimental data, establishing the ability of a numerical model to reproduce the physics of the problem.
The experimental study conducted by Kussoy and Horstman [18] involved a cone cylinder with a taper angle of degrees. The leading edge is a 10-degree cone that makes the transition through a circular arc to a constant-area cylinder with a radius of . The cone cylinder with a taper angle of degrees, designed to produce an oblique shock and an interaction zone between the boundary layer and the shock wave, is located downstream from the leading edge, as can be seen in Figure 17. As described in reference [19], the leading edge can be excluded from the CFD calculations, provided that the conditions of the input limits are slightly adjusted (the Mach number is changed from to . The wall (cylinder and cone) has a constant temperature of . The fluid domain and the boundary conditions are shown in the Figure 18.
![]() |
Figure 17. Solid body (from [18]) |
![]() |
Figure 18. Domain and boundaries |
Although this is an axisymmetric problem, its resolution is made here with a 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:
|
The non-slip condition is applied on solid boundaries. Then
|
The temperature at the wall, which has a known value , is specified using the dimensionless specific internal energy:
|
and the outlet boundary conditions are specified using and .
The initial conditions are given by:
|
These initial values are defined in all nodes of the finite element mesh, except for the nodes where Dirichlet boundary conditions were prescribed.
Figure 19 shows the side view of the finite element mesh, which has 1 600 000 hexahedral elements. The mesh is the same as found in the NASA web repository for turbulence models (2-D with 2-zones of 81101 and 81101 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 . Figure 20 shows a detailed view of the mesh where the flare starts.
![]() |
Figure 19. Finite element mesh |
![]() |
Figure 20. Mask in the zone where the cone begins |
The mesh was refined near the wall such that the estimate of is close to 0.5.
Using the safety factor , the dimensionless time step is . The artificial viscosity coefficient used here is .
Figure 21 compares the experimental velocity profiles with those obtained in this work and by the NASA WindUS program. The position in which these data were obtained corresponds to the plane (the origin of coordinates is located where the flare starts). Since there are no appreciable differences between the results obtained with TG, CBS and MTG, only the curve corresponding to the MTG scheme is shown here.
![]() |
Figure 21. Comparison of velocity profiles |
Figure 22 compares the temperature profiles for the same location. Both figures show proximity between the values obtained with the WindUS program [19] and the MTG code implemented in this work. In turn, results of both codes show a good agreement with the experimental data obtained in reference [18].
![]() |
Figure 22. Comparison of temperature profiles |
Figure 23 shows a comparison of the pressure distribution along the wall of the ogive. The coordinate is measured along the surface of the solid with its origin located when the flare starts.
![]() |
Figure 23. Pressure distribution on the wall |
Finally, Figure 24 shows a comparison of the heat transfer through the wall of the solid obstacle. The definition of wall heat transfer is given by: , where is the thermal conductivity of the fluid (but here it is calculated from the corresponding dimensionless quantities).
![]() |
Figure 24. Heat transfer through the wall |
A good agreement between results obtained with the MTG scheme and "WindUs" codes. As observed from Figure 23, the 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 approximately. The S-A model also under-predicts the heat transfer though the solid wall after position . The differences between the results of both codes and the experimental data for positions above can be attributed to the adopted RANS/FANS turbulence model (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.
Algorithm | Relative wall clock time |
---|---|
Single-step CBS | Not recommended |
TG | 1.0 |
TGM | 1.03 |
CBS | 1.30 |
It is observed that, for hypersonic flows, the three schemes (i.e. TG, MTG and CBS) are suitable from the point of view of convergence rates, but there are appreciable differences in terms of wall clock computational time, being the CBS scheme the most expensive.
The 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.
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 Δtin. 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.
[5] Zienkiewicz O.C., Taylor R.L. The finite element method. Vol. 3: Fluid Dynamics. Fifth edition, Butterworth Heinemann, Oxford, 2000.
[14] White F. Viscous fluid flow. McGraw Hill Book Co, First edition, New York, 1974.
[15] Huebner K.H., Thornton E.A., Byrom T.G. The finite element method for engineers. Third Edition, John Wiley and Sons, Inc. New York, 1995.
[17] Argyris J., Doltsinis I.S., Friz H. Hermes space shuttle: exploration of reentry aerodynamics. Computer Methods in Applied Mechanics and Engineering, 73:1-51, 1989.
Published on 26/09/19
Accepted on 19/09/19
Submitted on 25/05/19
Volume 35, Issue 3, 2019
DOI: 10.23967/j.rimni.2019.09.006
Licence: CC BY-NC-SA license
Are you one of the authors of this document?