Burbridge et al 2019a - Revision history

Rimni at 12:25, 21 May 2021

2021-05-21T12:25:52Z

Rimni at 12:25, 21 May 2021

2021-05-21T12:25:19Z

Rimni at 11:36, 12 March 2021

2021-03-12T11:36:43Z

Rimni: /* 7. Numerical simulations and comparative results */

2021-03-12T11:33:01Z

‎7. Numerical simulations and comparative results

Rimni: /* 7. Numerical simulations and comparative results */

2021-03-12T11:31:46Z

‎7. Numerical simulations and comparative results

Show changes

Rimni: /* 6. Explicit integration of element matrices */

2021-03-12T11:12:30Z

‎6. Explicit integration of element matrices

Rimni at 12:51, 18 March 2020

2020-03-18T12:51:59Z

Rimni at 09:42, 26 September 2019

2019-09-26T09:42:51Z

Rimni at 09:40, 26 September 2019

2019-09-26T09:40:12Z

Rimni at 09:37, 26 September 2019

2019-09-26T09:37:04Z

@@ Line 59: / Line 59: @@
 :* <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.

@@ Line 49: / Line 49: @@
 :* <math display="inline">{\tau }_{ij}\,</math> are the viscous components of the stress tensor.
-:* <math display="inline">e</math>''  ''is the specific total energy.
+:* <math display="inline">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.
@@ Line 109: / Line 109: @@
-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;"
@@ Line 117: / Line 117: @@
 |}
-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> and <math display="inline">{cb}_{2}</math> are constants of the S-A turbulence model.
 Finally, the production and destruction terms, in the conservative transport equation of S-A are given by

@@ Line 45: / Line 45: @@
 :* <math display="inline">\rho</math>  is the specific mass.
-:* <math display="inline">p</math>'' '' is the thermodynamic pressure.
+:* <math display="inline">p</math> is the thermodynamic pressure.
 :* <math display="inline">{\tau }_{ij}\,</math> are the viscous components of the stress tensor.
@@ Line 127: / Line 127: @@
 |}
-where the first term <math display="inline">P</math>'' ''represents the production and the second term'' D'' the destruction, being <math display="inline">{cb}_{1}</math>''<sub>,</sub> '' <math display="inline">S</math>'', '' <math display="inline">{cw}_{1}</math>'' ''and'' '' <math display="inline">{f}_{w}</math> constants and variables of the S-A turbulence model and <math display="inline">d</math> is the smallest distance from a field point (or node) to the nearest solid boundary (wall). The transport equation used here is expressed in its conservative form (Allmaras et al. [13]). Here the trip-term and the term of laminar suppression, both included in the original expressions of the model, are not considered.
+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.
 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
@@ Line 245: / 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.
-<span id='_Hlk7438137'></span>This method may be applied with implicit, semi-implicit and explicit time integration schemes [11]. An explicit scheme was implemented in this work in order to compare its behavior with respect to other explicit algorithms presented here.
+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].
@@ Line 398: / 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.
 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.
@@ Line 499: / 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.
 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].

@@ Line 667: / Line 667: @@
 |-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;"|1.0
-|-
+|-style="text-align:center"
 |  style="vertical-align: top;"|TG
 |  style="text-align: center;vertical-align: top;"|1.01
-|-
+|-style="text-align:center"
 |  style="vertical-align: top;"|MTG
 |  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
@@ Line 789: / Line 789: @@
 |-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="text-align:center"
 |  style="vertical-align: top;"|TG
 |  style="text-align: center;vertical-align: top;"|1.0
-|-
+|-style="text-align:center"
 |  style="vertical-align: top;"|MTG
 |  style="text-align: center;vertical-align: top;"|1.03
-|-
+|-style="text-align:center"
 |  style="vertical-align: top;"|CBS
 |  style="text-align: center;vertical-align: top;"|1.30
@@ Line 950: / Line 950: @@
 |-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="text-align:center"
 |  style="vertical-align: top;"|TG
 |  style="text-align: center;vertical-align: top;"|1.0
-|-
+|-style="text-align:center"
 |  style="vertical-align: top;"|TGM
 |  style="text-align: center;vertical-align: top;"|1.03
-|-
+|-style="text-align:center"
 |  style="vertical-align: top;"|CBS
 |  style="text-align: center;vertical-align: top;"|1.30

@@ Line 442: / Line 442: @@
 ==6. Explicit integration of element matrices==
-Eight nodes isoparametric hexahedral elements were used in this work. Figures 1 and 2 show the element in the physical and computational spaces, respectively. The interpolation functions of this element are given by the following expressions:
+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;"
@@ Line 455: / Line 455: @@
-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:
+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;"
@@ Line 481: / Line 481: @@
 |  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;"
 |-
@@ Line 488: / Line 488: @@
 | 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;"
 |-

@@ Line 3: / Line 3: @@
 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 turbulent, FEM, Taylor-Galerkin, CBS, Spalart-Allmaras
 ==1. Introduction==

@@ Line 369: / Line 369: @@
 |style="text-align: center;"|<math>{\boldsymbol{f}}_{i}^{\rho }=\left\{ \rho {v}_{i}\right\} ;\, {\boldsymbol{f}}_{ij}^{\rho v}=</math><math>\left\{ \rho {v}_{i}{v}_{j}+p{\delta }_{ij}\right\} ;\, {\boldsymbol{f}}_{i}^{\rho e}=</math><math>\left\{ \rho {ev}_{i}+p{v}_{i}\right\} ;\, {\boldsymbol{f}}_{i}^{\rho \tilde{\nu }}=</math><math>\left\{ \rho {v}_{i}\tilde{\nu }\right\}</math>
-<math>\begin{matrix}\displaystyle{\boldsymbol{q}}^{n}=-{cb}_{1}\, {{S}_{E}}^{n}\, {\left\{ \rho \tilde{\nu }\right\} }^{n}+{cw}_{1}\, {{{f}_{w}}_{E}}^{n}\, {\left\{ \rho \, {\left( \frac{\tilde{\nu }}{d}\right) }^{2}\right\} }^{n}-\frac{{cb}_{2}}{\sigma }\, {\left( {\rho }_{E}\, {\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right) }^{n}\left\{ 1\right\} +\\ \displaystyle +\, \frac{1}{\sigma }{\left[ \left( {\nu }_{E}+{\tilde{\nu }}_{E}\right) {\left. \frac{\partial \rho }{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right] }^{n}\left\{ 1\right\} \end{matrix}</math>
+<math>\begin{matrix}\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. \displaystyle {\boldsymbol{q}}_{{\mathit{\boldsymbol{v}}}_{\mathit{\boldsymbol{i}}}}\right| }^{n}=-{cb}_{1}\, {{S}_{E}}^{n}\, {\left\{ \rho \tilde{\nu }{v}_{i}\right\} }^{n}+{cw}_{1}\, {{{f}_{w}}_{E}}^{n}\, {\left\{ \rho \, {\left( \frac{\tilde{\nu }}{d}\right) }^{2}{v}_{i}\right\} }^{n}-\frac{{cb}_{2}}{\sigma }\, {\left( {\rho }_{E}\, {\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right) }^{n}{{\boldsymbol{v}}_{i}}^{n}+\\ \displaystyle+\, \frac{1}{\sigma }{\left[ \left( {\nu }_{E}+{\tilde{\nu }}_{E}\right) {\left. \frac{\partial \rho }{\partial {x}_{i}}\right| }_{E}{\left. \frac{\partial \tilde{\nu }}{\partial {x}_{i}}\right| }_{E}\right] }^{n}{{\boldsymbol{v}}_{i}}^{n}\end{matrix}</math>
+<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;width: 5px;text-align: right;white-space: nowrap;"|(20b)
@@ Line 388: / Line 388: @@
 <math>\displaystyle{{\boldsymbol{t}}_{j}}^{n}=\int_{{\Gamma }_{E}}^{}{{\boldsymbol{\phi }}^{\ast }}^{T}\, \left[ \eta \, \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, {{\boldsymbol{v}}_{j}}^{n}+\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}{{\boldsymbol{v}}_{i}}^{n}\right) +\, \lambda \, \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}\, {{\boldsymbol{v}}_{k}}^{n}\right) {\delta }_{ij}\right] {n}_{i}\, \, d\Gamma</math>
-<math>\begin{matrix}\displaystyle{\boldsymbol{h}}^{n}={\displaystyle\int}_{{\Gamma }_{E}}^{}{{\boldsymbol{\phi }}^{\ast }}^{T}\, \left( \boldsymbol{\phi }{{\boldsymbol{\, v}}_{j}}^{n}\right) \left[ \eta \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, {{\boldsymbol{v}}_{j}}^{n}+\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}{{\boldsymbol{v}}_{i}}^{n}\right) +\lambda \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{k}}\, {{\boldsymbol{v}}_{k}}^{n}\right) {\delta }_{ij}\right] {n}_{i}\, \, d\Gamma +\\ \displaystyle +\, {\displaystyle\int}_{{\Gamma }_{E}}^{}{{\boldsymbol{\phi }}^{\ast }}^{T}\, \left( {K+K}_{t}\right) \, \left( \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}{\boldsymbol{u}}^{n}\right) \, {n}_{i}\, \, d\Gamma \end{matrix}</math>
+<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;width: 5px;text-align: right;white-space: nowrap;"|(20c)

@@ Line 342: / Line 342: @@
 {| style="text-align: center; margin:auto;width: 100%;"
 |-
-|style="text-align: center;"|<math>{\boldsymbol{B}}_{i}=\int_{{\Omega }_{E}}^{}{\boldsymbol{\phi }}^{T}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ;\quad {\boldsymbol{C}}_{i}=\frac{{\Delta t}_{in}}{2}\int_{{\Omega }_{E}}^{}\left( \boldsymbol{\phi }{{\boldsymbol{v}}_{k}}^{n}\right) \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{k}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ;\, \, \boldsymbol{P}=</math><math>\frac{{\Delta t}_{in}}{2}\int_{{\Omega }_{E}}^{}\frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{j}}\frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}\, d\Omega \,</math>
+|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{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\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 \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{i}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{j}}\, d\Omega +\int_{{\Omega }_{E}}^{}\lambda \frac{\partial {\boldsymbol{\phi }}^{T}}{\partial {x}_{j}}\, \frac{\partial \boldsymbol{\phi }}{\partial {x}_{i}}\, d\Omega ,\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \hbox{if}\quad \, i\, \not =j\end{matrix}\right.</math>
+<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>

@@ Line 293: / Line 293: @@
 {| 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>
+|<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)
@@ Line 305: / Line 305: @@
 {| 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>
+|<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)
@@ Line 317: / Line 317: @@
 {| 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>
+|<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)
@@ Line 329: / Line 329: @@
 {| 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>
+|<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)