Study of Conservation on Implicit Techniques for Unstructured Finite Volume Navier-Stokes Solvers

ABSTRACT:

The work is a study of conservation on linearization techniques of time-marching schemes for the unstructured finite volume Reynolds-averaged Navier-Stokes formulation. The solver used in this work calculates the numerical flux applying an upwind discretization based on the flux vector splitting scheme. This numerical treatment results in a very large sparse linear system. The direct solution of this full implicit linear system is very expensive and, in most cases, impractical. There are several numerical approaches which are commonly used by the scientific community to treat sparse linear systems, and the point-implicit integration was selected in the present case. However, numerical approaches to solve implicit linear systems can be non-conservative in time, even for formulations which are conservative by construction, as the finite volume techniques. Moreover, there are physical problems which strongly demand conservative schemes in order to achieve the correct numerical solution. The work presents results of numerical simulations to evaluate the conservation of implicit and explicit time-marching methods and discusses numerical requirements that can help avoiding such non-conservation issues.

KEYWORDS:
Computational fluid dynamics; Time marching methods; Flux vector splitting scheme; Conservative discretization

INTRODUCTION

This work discusses issues associated with the coupling of implicit integration methods, for unstructured finite volume formulations, with the spatial discretization based on flux vector splitting schemes. The material discussed in the work is an extension of the work of Barth (1987). The original paper discussed the approximate local time linearizations of nonlinear terms for the finite difference formulation using Total Variation Diminishing (TVD) (Harten, 1983) and upwind algorithms. Here, the time conservation of the point-implicit integration for unstructured meshes is analyzed and discussed. Finite volume formulations have the tremendously important property of being conservative by construction. However, some time integration approaches may present numerical issues that can destroy this property, at least for unsteady applications or during the process of converging to a steady state. This shortcoming was observed when performing simulations of the flow inside closed systems in which there is no addition or extraction of fluid mass to/from the system. For such systems, the use of the point-implicit schemes discussed in this present paper has led to heat generation and non-conservation of mass in the interior of the computational domain. Results of simulations using explicit and implicit integration methods are presented in the present paper in order to better understand the non conservation issue of the time-marching methods. An analysis of the problem is also performed by a detailed study of the backward Euler method linearization.

THEORETICAL FORMULATION

The formulation used in the present work is based on the Reynolds-averaged Navier-Stokes set of equations, also known by the Computational Fluid Dynamics (CFD) community as the Reynolds Averaged Navier-Stokes (RANS) equations. The RANS equations are obtained by time filtering the Navier-Stokes set of equations. The compressible RANS equations are written in the algebraic vector form as

(1) $\frac{\partial \mathit{Q}}{\partial \mathit{t}}+\nabla ·(\mathit{Fe}-\mathit{Fv})=0.$

The conserved variables vector, Q, the inviscid flux vector, F_e, and the viscous flux vector, F_v, are given by

(2) $\mathit{Q}=\left[\begin{array}{c}\mathit{\rho }\\ \mathit{\rho }\mathit{u}\\ \mathit{\rho }\mathit{w}\\ \mathit{e}\end{array}\right],$

(3) $\mathit{Fe}=\left[\begin{array}{c}\mathit{\rho }\mathit{v}\\ \mathit{pu}\mathit{v}+\mathit{p}{\hat{\mathit{i}}}_{\mathit{x}}\\ \mathit{pv}\mathit{v}+\mathit{p}{\hat{\mathit{i}}}_{\mathit{y}}\\ \mathit{pw}\mathit{v}+\mathit{p}{\hat{\mathit{i}}}_{\mathit{z}}\\ \left(\mathit{e}+\mathit{p}\right)\mathit{v}\end{array}\right],$

(4) $\mathit{Fv}=\left[\begin{array}{c}0\\ \left({\mathit{\tau }}_{\mathit{xi}}\right){\hat{\mathit{i}}}_{\mathit{i}}\\ \left({\mathit{\tau }}_{\mathit{yi}}\right){\hat{\mathit{i}}}_{\mathit{i}}\\ \left({\mathit{\tau }}_{\mathit{zi}}\right){\hat{\mathit{i}}}_{\mathit{i}}\\ \left({\mathit{\beta }}_{\mathit{i}}\right){\hat{\mathit{i}}}_{\mathit{i}}\end{array}\right],$

in which ρ stands for density, v = {u,v,w} for the velocity vector in Cartesian coordinates, p for static pressure, τ for viscous stress tensor, q_H for heat flux vector, e for the total energy per unit volume and β_i is given by

(5) ${\beta }_{\mathit{i}}={\mathit{\tau }}_{\mathit{ij}}{\mathit{u}}_{\mathit{j}}-{\mathit{q}}_{\mathit{Hi}·}$

Theî_x, î_y and î_z terms are the Cartesian-coordinate orthonormal vector basis. It is very important to emphasize that field forces, such as gravity, are neglected here.

Other equations are necessary in order to close the system of equations given by Eq. (4). These additional equations are called constitutive relations. The first constitutive equation presented to close the Navier-Stokes set of equations is known as the equation of state. This equation considers the perfect gas law, and it is written as

(6) $\mathit{p}=(\gamma -1)\left[\mathit{e}-\frac{1}{2}\mathit{\rho }({\mathit{u}}^2+{\mathit{v}}^2+{\mathit{w}}^2)\right],$

in which the mean total energy per unit volume, e, is given by

(7) $\mathit{e}=p\left[{\mathit{e}}_{\mathit{i}}+\frac{1}{2}({\mathit{u}}^2+{\mathit{v}}^2+{\mathit{w}}^2)\right],$

and ei stands for internal energy, defined as

(8) ${\mathit{e}}_{\mathit{i}}={\mathit{C}}_{\mathit{v}}\mathit{T},$

in which T stands for the mean static temperature and C_v is the specific heat at constant volume. The heat flux from Eq. (5) is obtained from the Fourier law for heat conduction, and it is given by

(9) ${\mathit{q}}_{{\mathit{H}}_{\mathit{j}}}=-\frac{\gamma \mu }{\mathrm{Pr}}\frac{\partial ({\mathit{e}}_{\mathit{i}})}{\partial {\mathit{x}}_{\mathit{j}}},$

in which γ is the ratio of specific heats and Pr is the Prandtl number. Typically, for air, it is assumed that γ= 1.4 and Pr = 0.72. C_p is the gas specific heat at constant pressure and µ is the dynamic molecular viscosity coefficient, calculated as a function of the temperature by the Sutherland law equation (Anderson, 1991), written as

(10) $\mu ={\mu }_{\infty }{\left(\frac{\mathit{T}}{{\mathit{T}}_{\infty }}\right)}^{\frac{2}{3}}\frac{{\mathit{T}}_{\infty }+\mathbf{S}}{\mathit{T}+\mathbf{S}},$

where S = 110K , and µ_∞ is the dynamic molecular viscosity coefficient of the fluid at temperature T_∞.

The components of the viscous stress tensor, for a Newtonian fluid, are given by

(11) ${\mathit{\tau }}_{\mathit{ij}}=\mu \left[\left(\frac{\partial {\mathit{u}}_{\mathit{i}}}{\partial {\mathit{x}}_{\mathit{j}}}+\frac{\partial {\mathit{u}}_{\mathit{j}}}{{\partial }_{{\mathit{x}}_{\mathit{i}}}}\right)-\frac{2}{3}\frac{\partial {\mathit{u}}_{\mathit{m}}}{\partial {\mathit{x}}_{{\mathit{x}}_{\mathit{m}}}}{\mathit{\delta }}_{\mathit{ij}}\right],$

in which δ_ij stands for the Kronecker delta.

NUMERICAL FORMULATION

The numerical formulation applied in this work is briefly presented in this section. The study is performed using the RANS equations discretized in the context of a cell-centered finite-volume. The solver used in the present work is Le Michigan Aerothermodynamics Navier-Stokes Solver (LeMANS). The original framework was developed by Scalabrin (2007) at the University of Michigan to simulate hypersonic flows over reentry capsules. It is a three dimensional (3-D) numerical solver that uses the finite volume formulation for unstructured meshes to solve the RANS equations coupled with non-equilibrium chemical reaction equations and the Spalart-Allmaras (SA) turbulence model (Spalart and Allmaras, 1992; 1994). The numerical flux is calculated using an upwind scheme based on the Steger-Warming flux vector splitting (Steger and Warming, 1981). The time integration is performed by point implicit and Runge-Kutta methods. This section describes the numerical formulations and discretizations applied in the present work.

FINITE VOLUME FORMULATION

The finite volume formulation is a numerical method applied to represent and evaluate partial differential equations. It is applied by the CFD community to find the solution of the RANS, equations. The method is obtained integrating the flow equations for each control volume within a given mesh. The RANS equations, written in the context of a cell-centered finite-volume formulation, are given by

(12) ${\int }_{{\mathit{V}}_{\mathit{i}}}\frac{\partial \mathit{Q}}{\partial \mathit{t}}\mathit{dV}+{\int }_{{\mathit{V}}_{\mathit{i}}}▽·\left(\mathit{Fe}-\mathit{Fv}\right)\mathit{dV}=0.$

Considering a cell-centered formulation, V_i is a given cell of the given grid. After the integration it is possible to apply the Gauss theorem over the equation above, resulting in:

(13) ${\int }_{{\mathit{V}}_{\mathit{i}}}\frac{\partial \mathit{Q}}{\partial \mathit{t}}\mathit{dV}+{\int }_{{\mathit{S}}_{\mathit{i}}}\left(\mathit{Fe}-\mathit{Fv}\right).\mathit{n}\mathit{dS}=0,$

in which S_i is the outward-oriented area vector and it is defined as

(14) ${\mathit{S}}_{\mathit{i}}=\{{\mathit{S}}_{\mathit{x}},{\mathit{S}}_{\mathit{y}},{\mathit{S}}_{\mathit{z}}\mathrm{\}.}$

Considering the mean value of the conserved variables within the i-th control volume, one can write the first term of Eq. (13) as

(15) ${\mathit{Q}}_{\mathit{i}}=\frac{1}{{\mathit{V}}_{\mathit{i}}}{\int }_{{\mathit{V}}_{\mathit{i}}}{\mathit{QdV}}_{\mathit{i}}.$

The second term of Eq. (13) can be written as the sum of all faces of a cell

(16) ${\int }_{{\mathit{S}}_{\mathit{i}}}(\mathit{Fe}+\mathit{Fv}\mathrm{).}\mathit{n}\mathit{dS}=\sum _{\mathit{k}=1}^{\mathit{nf}}\left(\overrightarrow{\mathbf{F}{\mathit{e}}_{\mathit{k}}}-\overrightarrow{\mathbf{F}{\mathit{v}}_{\mathit{k}}}\right).\overrightarrow{{\mathit{n}}_{\mathit{k}}}{\mathit{S}}_{\mathit{k}},$

in which the k subscript is the index of the cell face, and nf indicates the number of faces of the i-th volume. Finally, the RANS equations discretized with a finite volume approximation are given by

(17) $\frac{\partial {\mathit{Q}}_{\mathit{i}}}{\partial \mathit{t}}=-\frac{1}{{\mathit{V}}_{\mathit{i}}}\sum _{\mathit{k}=1}^{\mathit{nf}}(\overrightarrow{{\mathbf{Fe}}_{\mathit{k}}}-\overrightarrow{{\mathbf{Fv}}_{\mathit{k}}})·\overrightarrow{{\mathit{n}}_{\mathit{k}}}{\mathit{S}}_{\mathit{k}}.$

For this formulation, the fluxes are computed at the faces of the control volume, and the conserved variables are computed in the cell.

INVISCID FLUX CALCULATION

The inviscid fluxes are calculated using a method based on a classical flux vector splitting formulation, the Steger-Warming Scheme (SW) (Steger and Warming, 1981). This method is an upwind scheme which uses the homogeneous property of the inviscid flux vectors to write

(18) $\overrightarrow{{\mathbf{Fe}}_{\mathit{k}}}.\overrightarrow{{\mathit{n}}_{\mathit{k}}}={\mathit{Fe}}_{\mathit{n}}=\frac{\partial {\mathit{Fe}}_{\mathit{n}}}{\partial \mathit{Q}}\mathit{Q}=\mathit{AQ},$

where Fe_n is the normal flux at the k-th face, and A is the Jacobian matrix of the inviscid flux which can be diagonalized by the matrices of its eigen vectors from the left and from the right, L and R

(19) $\mathit{A}=\mathit{R}\Lambda \mathit{L},$

and Λ is the diagonal matrix of the eigenvalues of the Jacobian matrix. The A matrix can be split into positive and negative parts as

(20) ${\mathit{A}}^{+}=\mathit{R}{\Lambda }^{+}\mathit{L}\mathit{and}{\mathit{A}}^{-}=\mathit{R}{\Lambda }^{-}\mathit{L}.$

The splitting separates the flux into two parts, the downstream and the upstream flux, in relation to the face orientation as:

(21) $\overrightarrow{\mathit{Fe}}·\overrightarrow{\mathit{n}}={\mathit{Fe}}_{\mathit{n}}^{+}+{\mathit{Fe}}_{\mathit{n}}^{-}=({\mathit{A}}_{\mathit{cl}}^{+}{\mathit{Q}}_{\mathit{cl}}+{\mathit{A}}_{\mathit{cr}}^{-}{\mathit{Q}}_{\mathit{cr}}\mathrm{),}$

where the cl and cr subscripts are the cells on the left and right sides of the face. The split eigen values of the Jacobian matrix are given by

(22) ${\lambda }^{\pm }=\frac{1}{2}(\lambda \pm |\lambda \mathrm{|),}$

In order to avoid sudden sign transition, and then discontinous derivative, the split eigen values receive a small number, ε, turning Eq. (22) into

(23) ${\lambda }^{\pm }=\frac{1}{2}\left(\mathit{\lambda }\pm \sqrt{{\mathit{\lambda }}^2+{\varepsilon }^2}\right).$

The soften sign transition turns the derivative continous at the transition point.

Numerical studies performed in the present work indicated that this flux vector splitting is too dissipative and it can deteriorate the boundary layer profiles (Junqueira-Junior, 2012; Junqueira-Junior et al., 2013; MacCormack and Candler, 1989). Therefore, a pressure switch is implemented to smoothly shift the Steger-Warming scheme into a centered one. Then, the artificial dissipation is controlled and the numerical stability is maintained, as presented in the following formulation:

(24) $\overrightarrow{{\mathit{Fe}}_{\mathit{k}}}·\overrightarrow{{\mathit{n}}_{\mathit{k}}}={\mathit{F}}_{\mathit{k}}^{+}+{\mathit{F}}_{\mathit{k}}^{-}=({\mathit{A}}_{\mathit{k}+}^{+}{\mathit{Q}}_{\mathit{k}+}+{\mathit{A}}_{\mathit{k}-}^{-}{\mathit{Q}}_{\mathit{k}-}\mathrm{).}$

in which

(25) $\begin{array}{c}{\mathit{Q}}_{\mathit{k}+}=(1-\mathit{w}){\mathit{Q}}_{\mathit{cl}}+{\mathit{wQ}}_{\mathit{cr}},\mathit{and},\hfill \\ {\mathit{Q}}_{{\mathit{k}}^{-}}=(1-\mathit{w}){\mathit{Q}}_{\mathit{cr}}+{\mathit{wQ}}_{\mathit{cl}}·\hfill \end{array}$

The switch, w, is given by

(26) $\mathit{w}=\frac{1}{2}\frac{1}{{\left(\mathit{\alpha }▽\mathit{p}\right)}^2+1}\mathit{and}▽\mathit{p}=\frac{\left|{\mathit{p}}_{\mathit{cl}}-{\mathit{p}}_{\mathit{cr}}\right|}{\min \left({\mathit{p}}_{\mathit{cl}},{\mathit{p}}_{\mathit{cr}}\right)},$

where ∇p is a scalar number, a numerical approximation of the pressure gradient. For small ∇p, w = (1 − w) = 0.5, the code runs with a centered scheme, and for larger values of ∇p, w = 0 and (1 − w) = 1, the code runs with the original Steger-Warming scheme. For Eq. (26) one suggests α = 6, but some problems may require larger values (Scalabrin, 2007).

The applied formulation was originally created with interest on studying flows over reentry capsules. For such particular cases, with very strong shock waves, it is very common to find solutions with numerical and non-physical structures such as carbuncles (Ramalho et al., 2011). To avoid such numerical problems, artificial dissipation has to be added to the method. The dissipation was included into the split eigenvalues, Eq. (23), using an ε factor, which is given by:

(27) ${\varepsilon }_{\mathit{k}}=\left\{\begin{array}{cc}0.3({\mathit{a}}_{\mathit{k}}+|\overrightarrow{{\mathit{u}}_{\mathit{k}}}\mathrm{|)}\hfill & {\mathit{d}}_{\mathit{k}}>{\mathit{d}}_0,\hfill \\ 0.3(1-|\overrightarrow{{\mathit{n}}_{\mathit{k}}}·\overrightarrow{{\mathit{m}}_{\mathit{k}}}\mathrm{|)(}{\mathit{a}}_{\mathit{k}}+|\overrightarrow{{\mathit{u}}_{\mathit{k}}}\mathrm{|)}\hfill & {\mathit{d}}_{\mathit{k}}<{\mathit{d}}_0,\hfill \end{array}\right.$

where d_k is the distance of the k-th face, to the nearest wall boundary. d₀ is set by the user and must be smaller than the boundary layer thickness and larger than the shock stand-off distance.⃗m_k is the normal vector of the nearest wall, and ⃗n_k is the normal vector to the k-th face. Equation (27) applies the term (1 -|⃗n_k . ⃗m_k|) to decrease the value at the faces parallel to the wall inside the boundary layer (Scalabrin, 2007). This artificial dissipation model has shown an important role in the prediction of boundary layer profiles (Junqueira-Junior, 2012; Junqueira-Junior et al., 2013).

VISCOUS FLUX CALCULATION

The viscous terms are based on derivative of properties on the faces. To build the derivative terms, two volumes are created over the face where the derivative is being calculated. At the center of each new volume, the derivative is calculated using the Green-Gauss theorem. This computation is used to find the derivative at the desired face.

A two dimensional (2-D) example is used in this section to better explain the derivative calculation. Consider the two cells, S1 and S2, in Fig. 1. Two new cells, S3 and S4, are created using node points, P1 and P3, and cell centered points, P2 and P4, to calculate the derivative on 1-3 face. The properties at the faces are calculated using simple averages. For example, in Figs. 1 and 2, the properties are given by

Figure 1
2-D example of new volume creation.

Figure 2
2-D example of derivative calculation.

(28) $\begin{array}{c}{\mathit{Q}}_{12}=\frac{1}{2}({\mathit{Q}}_1+{\mathit{Q}}_2)\hfill \\ {\mathit{Q}}_{23}=\frac{1}{2}({\mathit{Q}}_2+{\mathit{Q}}_3)\hfill \\ {\mathit{Q}}_{13}={\mathit{Q}}_{31}=\frac{1}{2}({\mathit{Q}}_1+{\mathit{Q}}_3)\hfill \\ {\mathit{Q}}_{34}=\frac{1}{2}({\mathit{Q}}_3+{\mathit{Q}}_4)\hfill \\ {\mathit{Q}}_{41}=\frac{1}{2}({\mathit{Q}}_4+{\mathit{Q}}_1)\hfill \end{array}$

(29) ${\int }_{\mathit{V}}▽\mathit{Q}\mathit{dV}={\oint }_{\mathit{S}}\overrightarrow{Q}·\overrightarrow{\mathit{n}}\mathit{dS}.$

Considering ∇Q as a constant over the cell, the equation above yields

(30) $\nabla \stackrel{─}{\mathit{Q}}=\frac{1}{\mathit{V}}{\oint }_{\mathit{S}}\overrightarrow{\mathbf{Q}}·\overrightarrow{\mathit{n}}\mathit{dS},$

in which ∇Q is the constant cell-centered gradient. Using the derivatives in the S3 and S4 cells, the derivative at faces 1-3 is computed using

(31) $\nabla {\mathit{Q}}_{13}=\frac{{\mathit{V}}_3\nabla {\mathit{Q}}_7+{\mathit{V}}_4\nabla {\mathit{Q}}_8}{{\mathit{V}}_3+{\mathit{V}}_4}.$

The derivative computation for other types of element, 2-D or 3-D, is straightforward.

TIME INTEGRATION

The explicit second-order Runge-Kutta (Lomax et al., 2001) and the point implicit scheme are the two time-marching methods applied in the present work.

The second-order Runge-Kutta integration method used in this work is given by

(32) $\begin{array}{c}\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}=-\frac{\Delta \mathit{t}}{{\mathit{V}}_{\mathit{cl}}}{\left[\sum _{\mathit{k}=1}^{\mathit{nf}}(\overrightarrow{\mathbf{F}{\mathit{e}}_{\mathit{k}}}-\overrightarrow{\mathbf{F}{\mathit{v}}_{\mathit{k}}}\mathrm{).}\overrightarrow{{\mathit{n}}_{\mathit{k}}}{\mathit{S}}_{\mathit{k}}\right]}^{\mathit{n}}\hfill \\ {\mathit{Q}}_{\mathit{cl}}^{\stackrel{─}{\mathit{n}+1}}={\mathit{Q}}_{\mathit{cl}}^n+\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}.\hfill \\ \Delta {\mathit{Q}}_{\mathit{cl}}^{\stackrel{─}{\mathit{n}+\perp }}=-\frac{\Delta \mathit{t}}{{\mathit{V}}_{\mathit{cl}}}{\left[\sum _{\mathit{k}=1}^{\mathit{nf}}(\overrightarrow{\mathbf{F}{\mathit{e}}_{\mathit{k}}}-\overrightarrow{\mathbf{F}{\mathit{v}}_{\mathit{k}}}\mathrm{).}\overrightarrow{{\mathit{n}}_{\mathit{k}}}{\mathit{S}}_{\mathit{k}}\right]}^{\stackrel{─}{\mathit{n}+1}}\hfill \\ {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}+\perp }=\frac{1}{2}({\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}+{\mathit{Q}}_{\mathit{cl}}^{\stackrel{─}{\mathit{n}+1}}+\Delta {\mathit{Q}}_{\mathit{cl}}^{\stackrel{─}{\mathit{n}+1}}\mathrm{).}\hfill \end{array}$

In order to simplify the forthcoming equations, the right hand side of Eq. (32) is written as

(33) $\begin{array}{c}{\mathit{R}}_{\mathit{cl}}^{\mathit{n}}=-\frac{\Delta \mathit{t}}{{\mathit{V}}_{\mathit{cl}}}{\left[\sum _{\mathit{k}=1}^{\mathit{nf}}(\overrightarrow{\mathbf{F}{\mathit{e}}_{\mathit{k}}}-\overrightarrow{\mathbf{F}{\mathit{v}}_{\mathit{k}}}\mathrm{).}\overrightarrow{{\mathit{n}}_{\mathit{k}}}{\mathit{S}}_{\mathit{k}}\right]}^{\mathit{n}}\hfill \\ {\mathit{R}}_{\mathit{cl}}^{\stackrel{─}{\mathit{n}+1}}=-\frac{\Delta \mathit{t},}{{\mathit{V}}_{\mathit{cl}}}\left[\sum _{\mathit{k}=1}^{\mathit{nf}}\left(\overrightarrow{\stackrel{─}{{\mathit{Fe}}_{\mathit{k}}}}-\overrightarrow{\stackrel{─}{{\mathit{Fv}}_{\mathit{k}}}}\right).\overrightarrow{{\mathit{n}}_{\mathit{k}}}{\mathit{S}}_{\mathit{k}}\right]\underset{}{\overset{\stackrel{─}{\mathit{n}+1}}{,}}\hfill \end{array}$

in which R_cl is the residue of the i-th cell. The implicit integration applied in the present work is based on the backward Euler method, which is given by.

(34) $\frac{\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}}{\Delta \mathit{t}}{\mathit{V}}_{\mathit{cl}}={\left[-\sum _{\mathit{k}=1}^{\mathit{nf}}(\overrightarrow{\mathbf{F}{\mathit{e}}_{\mathit{k}}}-\overrightarrow{\mathbf{F}{\mathit{v}}_{\mathit{k}}})·\overrightarrow{{\mathit{n}}_{\mathit{k}}}{\mathit{S}}_{\mathit{k}}\right]}^{\mathit{n}+1}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}+1}$

One can linearize the residue at time n + 1 as a function of properties at time n.

(35) $\begin{array}{c}\frac{\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}}{\Delta \mathit{t}}{\mathit{V}}_{\mathit{cl}}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}}-{\sum }_{\mathit{k}=1}^{\mathit{nf}}\left\{{\left[\frac{\partial (\overrightarrow{\mathit{F}\mathit{e}}.\overrightarrow{\mathit{n}})}{\partial \mathit{Q}}\right]}_{\mathit{k}}^{\mathit{n}}\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}\right.\\ \left.-{\left[\frac{\partial \left(\overrightarrow{\mathit{F}\mathit{v}}.\overrightarrow{\mathit{n}}\right)}{\partial \mathit{Q}}\right]}_{\mathit{k}}^{\mathit{n}}\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}\right\}.{\mathit{S}}_{\mathit{k}}.\end{array},$

From the spatial discretization, the inviscid terms can be written as

(36) $\begin{array}{c}\hfill \left[\frac{\partial (\overrightarrow{\mathit{F}\mathit{e}}·\overrightarrow{\mathit{n}})}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}=\left[\frac{\partial \left({\overrightarrow{\mathit{F}\mathit{e}}}^{+}.\overrightarrow{\mathit{n}}\right)}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}\\ \hfill +\left[\frac{\partial \left({\overrightarrow{\mathit{F}\mathit{e}}}^{-}.\overrightarrow{\mathit{n}}\right)}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}.\end{array}$

It is a common practice to assume

(37) $\frac{\partial \left({\overrightarrow{\mathit{F}\mathit{e}}}^{\pm }·\overrightarrow{\mathit{n}}\right)}{\partial \mathit{Q}}={\mathit{A}}^{\pm },$

and then write

(38) $\begin{array}{c}\left[\frac{\partial ({\overrightarrow{\mathit{F}\mathit{e}}}^{+}·\overrightarrow{\mathit{n}})}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}={\mathit{A}}_{\mathit{k}+}^{+}\Delta {\mathit{Q}}_{\mathit{cl}}\hfill \\ \left[\frac{\partial ({\overrightarrow{\mathit{F}\mathit{e}}}^{-}·\overrightarrow{\mathit{n}})}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}={\mathit{A}}_{\mathit{k}}▵{\mathit{Q}}_{\mathit{cr},\mathit{k},}\hfill \end{array}$

which is not true and different from the real definition of the Jacobian matrix written in Eq. (20). Using the approximate form, Eq. (38) to calculate the inviscid terms may decrease the numerical stability of the method. Hence, in this work, the true Jacobian matrices of the split fluxes, given by Eq. (20), are implemented in order to calculate the implicit operator. The approximate Jacobian matrices, which are calculated using Eqs. (37) and (38), are implemented at the right hand side of the linear system. Issues involving the true Jacobians matrices have a major importance in the context of numerical stability for computational methods (Anderson et al., 1986; Hirsch, 1990; Steger and Warming, 1981).

The viscous terms can be written in the same form as

(39) $\begin{array}{c}\hfill \left[\frac{\partial (\overrightarrow{\mathit{F}\mathit{v}}·\overrightarrow{\mathit{n}})}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta \mathit{Qcl}=\left[\frac{\partial ({\overrightarrow{\mathit{F}\mathit{v}}}^{-}·\overrightarrow{\mathit{n}})}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}\\ \hfill -\left[\frac{\partial ({\overrightarrow{\mathit{F}\mathit{v}}}^{+}·\overrightarrow{\mathit{n}})}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}.\end{array}$

In this work, the viscous Jacobian matrices are represented by B. Hence, the Jacobian matrix splitting is written as

(40) $\left[\frac{\partial (\overrightarrow{\mathbf{F}\mathit{v}}·\to \mathit{n})}{\partial \mathit{Q}}\right]\underset{\mathit{k}}{\overset{}{}}\Delta {\mathit{Q}}_{\mathit{cl}}={\mathit{B}}_{{\mathit{k}}^{-}}^{-}\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k}}-{\mathit{B}}_{\mathit{k}+}^{+}\Delta {\mathit{Q}}_{\mathit{cl}}.$

One can write the system as

(41) $\begin{array}{c}\left[\frac{{\mathit{V}}_{\mathit{cl}}}{\Delta \mathit{t}}+\sum _{\mathit{k}=1}^{\mathit{nf}}({\mathit{A}}_{\mathit{k}+}^{+}+{\mathit{B}}_{\mathit{k}+}^{+}){\mathit{S}}_{\mathit{k}}\right]\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}+\hfill \\ \left[\sum _{\mathit{k}=1}^{\mathit{nf}}({\mathit{A}}_{{\mathit{k}}^{-}}^{-}-{\mathit{B}}_{{\mathit{k}}^{-}}^{-}){\mathit{S}}_{\mathit{k}}\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k}}^{\mathit{n}}\right]={\mathit{R}}_{\mathit{cl}}^{\mathit{n}}\hfill \end{array}$

It is, then, possible to write

(42) ${\mathit{M}}_{\mathit{cl}}\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}+\sum _{\mathit{k}=1}^{\mathit{nf}}{\mathit{N}}_{\mathit{k}}^{-}\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k}}^{\mathit{n}}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}},$

with

(43) ${\mathit{M}}_{\mathit{cl}}=\frac{{\mathit{V}}_{\mathit{cl}}}{\Delta \mathit{t}}+\sum _{\mathit{k}=1}^{\mathit{nf}}{\mathit{N}}_{\mathit{k}}^{+},$

(44) ${\mathit{N}}_{\mathit{k}}^{+}=({\mathit{A}}_{{\mathit{k}}^{+}}^{+}+{\mathit{B}}_{{\mathit{k}}^{+}}^{+}){\mathit{S}}_{\mathit{k}}.$

and

(45) ${\mathit{N}}_{\mathit{k}}^{-}=({\mathit{A}}_{{\mathit{k}}^{-}}^{-}-{\mathit{B}}_{{\mathit{k}}^{+}}^{+}){\mathit{S}}_{\mathit{k}}.$

As the code is an unstructured solver, this system of equations results in a sparse block matrix, where each block is a square matrix of size equal to the number of equations to be solved in each control volume. The solution of such system is typically very expensive and, depending on the size of the mesh, it is not even practical. A less expensive implicit method is applied in the present work, the implicit point integration (Gnoffo, 2003; Venkatakrishnan, 1995; Wright, 1997).

The main idea of the implicit point integration is to move all the off-diagonal terms to the right hand side and to solve the resulting system iteratively, i.e.,

(46) ${\mathit{M}}_{\mathit{cl}}\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}+1,\mathit{p}}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}}-\sum _{\mathit{k}=1}^{\mathit{nf}}{\mathit{N}}_{\mathit{k}}^{-}\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k}}^{\mathit{n}+1,\mathit{p}-1}.$

It is assumed that ΔQ^n+1.0 = 0 and four iterations are taken in the process as suggested in the literature (Wright, 1997). The given sparse linear system illustrates the point

$\begin{array}{c}\left[\begin{array}{ccccc}⬜& & ⬜& & \\ & ⬜& & ⬜& ⬜\\ ⬜& & ⬜& & ⬜\\ & ⬜& & ⬜& \\ & ⬜& ⬜& & ⬜\end{array}\right]\left[\begin{array}{c}\\ \\ \Delta {\mathit{Q}}^{(\mathit{n}+1\mathrm{),}\mathit{p}}\\ \\ \end{array}\right]=\left[\begin{array}{c}\hfill \\ \hfill \\ {\mathit{R}}^{\mathit{n}}\hfill \\ \hfill \\ \hfill \end{array}\right]\hfill \\ \left[\begin{array}{ccccc}⬜& & & & \\ & ⬜& & & \\ & & ⬜& & \\ & & & ⬜& \\ & & & & ⬜\end{array}\right]\left[\begin{array}{c}\\ \\ \Delta {\mathit{Q}}^{(\mathit{n}+1\mathrm{),}\mathit{p}}\\ \\ \end{array}\right]=\hfill \\ \left[\begin{array}{c}\hfill \\ \hfill \\ {\mathit{R}}^{\mathit{n}}\hfill \\ \hfill \\ \hfill \end{array}\right]-\left[\begin{array}{ccccc}& & ⬜& & \\ & & & ⬜& ⬜\\ ⬜& & & & ⬜\\ & ⬜& & & \\ & ⬜& ⬜& & \end{array}\right]\left[\begin{array}{c}\\ \\ \Delta {\mathit{Q}}^{(\mathit{n}+1\mathrm{),}\mathit{p}-1}\\ \\ \end{array}\right]\hfill \end{array}$

where each □ is a 5 × 5 block matrix.

The time step is computed by

(47) $\Delta \mathit{t}=\mathit{CFL}\frac{\ell }{‖\overrightarrow{\mathit{v}}‖+\mathit{a}},$

in which CFL is a parameter set to ensure the stability of the time integration method, l is the size of the cell and ||⃗υ|| + α is the largest wave speed in the cell.

BOUNDARY CONDITIONS

The boundary conditions are implemented using ghost cells. The solver creates the ghost cells to hold properties that satisfy the correct flux calculation at the boundaries. The implementation assigns properties that satisfy the Euler boundary conditions to calculate the inviscid fluxes, and properties that satisfy the Navier-Stokes boundary conditions to calculate the viscous fluxes. Therefore, the ghost volumes store two different types of fluxes for the correct computation of the RANS equations.

WALL INVISCID BOUNDARY CONDITIONS

The ghost cells hold the properties in the same manner to calculate the inviscid fluxes at the wall and at the symmetry boundaries. Mass and energy fluxes should yield zero, and the momentum flux is equal to the pressure flux. This is accomplished by setting the normal velocity component to the boundary face zero. In the present work, the interior domain is defined as the properties at the left side of the boundary face and the ghost cells are defined as the properties at the right side of the boundary.

In order to simplify the implementation, the properties at the left side of the boundary face are rotated to the face coordinates using

(48) ${\mathit{Q}}_{\mathit{cl}}^{\mathit{root}}=ℝ{\mathit{Q}}_{\mathit{cl}}.$

The vector of conservative properties, Q, is written as

(49) $\mathit{Q}={\left[\rho \rho \mathit{u}\rho \mathit{v}\rho \mathit{w}\mathit{e}\right]}^{\mathit{T}}$

In Eq. (48), ℝ is the rotation matrix given by,

(50) $ℝ=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& {\mathit{n}}_{\mathit{x}}& {\mathit{n}}_{\mathit{y}}& {\mathit{n}}_{\mathit{z}}& 0\\ 0& {\mathit{t}}_{\mathit{x}}& {\mathit{t}}_{\mathit{y}}& {\mathit{t}}_{\mathit{z}}& 0\\ 0& {\mathit{r}}_{\mathit{x}}& {\mathit{r}}_{\mathit{y}}& {\mathit{r}}_{\mathit{z}}& 0\\ 0& 0& 0& 0& 1\end{array}\right]$

and the ⃗n, ⃗t, ⃗r vectors define the face-based reference frame. The properties at the ghost cells are set to

(51) $\begin{array}{ccc}{\rho }_{\mathit{cr}}& =& {\rho }_{\mathit{cl}}\\ {\rho }_{\mathit{cr}}{\mathit{u}}_{\mathit{cr},\mathit{n}}^{\mathit{rot}}& =& {\rho }_{\mathit{cr}}{\mathit{u}}_{\mathit{cr},\mathit{n}}^{\mathit{rot}}\\ {\rho }_{\mathit{cr}}{\mathit{u}}_{\mathit{cr},\mathit{t}}^{\mathit{rot}}& =& {\rho }_{\mathit{cr}}{\mathit{u}}_{\mathit{cr},\mathit{t}}^{\mathit{rot}}\\ {\rho }_{\mathit{cr}}{\mathit{u}}_{\mathit{cr},\mathit{r}}^{\mathit{rot}}& =& {\rho }_{\mathit{cr}}{\mathit{u}}_{\mathit{cr},\mathit{r}}^{\mathit{rot}}\\ {\left({\mathit{e}}_{\mathit{i}}\right)}_{\mathit{cr}}& =& {\left({\mathit{e}}_{\mathit{i}}\right)}_{\mathit{cl}}\end{array}$

One can write in the matrix form as

(52) ${\mathit{Q}}_{\mathit{cr}}^{\mathit{rot}}=\mathit{W}{\mathit{Q}}_{\mathit{cl}}^{\mathit{rot}},$

in which W is the inviscid wall matrix given by

(53) $\mathit{W}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& -1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right].$

Therefore, the boundary condition can be written as

(54) ${\mathit{Q}}_{\mathit{cr}}={ℝ}^{-1}\mathit{W}ℝ{\mathit{Q}}_{\mathit{cl}},$

in which the R⁻¹ matrix is given by

(55) ${ℝ}^{-1}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& {\mathit{n}}_{\mathit{x}}& {\mathit{t}}_{\mathit{x}}& {\mathit{r}}_{\mathit{x}}& 0\\ 0& {\mathit{n}}_{\mathit{y}}& {\mathit{t}}_{\mathit{y}}& {\mathit{r}}_{\mathit{y}}& 0\\ 0& {\mathit{n}}_{\mathit{z}}& {\mathit{t}}_{\mathit{z}}& {\mathit{r}}_{\mathit{z}}& 0\\ 0& 0& 0& 0& 1\end{array}\right]$

It returns the properties to the Cartesian coordinate frame.

VISCOUS WALL BOUNDARY CONDITION WITH SPECIFIED TEMPERATURE

Viscous wall with specified temperature does not have necessarily zero heat conduction at the boundary face. For this boundary condition the user provides the wall temperature, T_wall , and the wall pressure is extrapolated from the interior in order to satisfy the zero normal wall pressure gradient condition Schlichting (1978). Hence,

(56) $\begin{array}{c}{\mathit{T}}_{\mathit{wall}}={\mathit{T}}_{\mathit{cl}}\hfill \\ {\mathit{p}}_{\mathit{wall}}={\mathit{p}}_{\mathit{cl}}.\hfill \end{array}$

One can computate the conservative properties at the wall boundary using the Cartesian components of the wall velocity, u_wall, v_wall and, w_wall which are set by the user.

(57) $\begin{array}{ccc}{\mathit{\rho }}_{\mathit{wall}}& =& \Re \frac{{\mathit{p}}_{\mathit{wall}}}{{\mathit{T}}_{\mathit{wall}}}\\ {\left(\mathit{\rho }\mathit{u}\right)}_{\mathit{wall}}& =& {\rho }_{\mathit{wall}}{\mathit{u}}_{\mathit{wall}}\\ {\left(\mathit{\rho }\mathit{v}\right)}_{\mathit{wall}}& =& {\rho }_{\mathit{wall}}{\mathit{v}}_{\mathit{wall}}\\ {\left(\mathit{\rho }\mathit{w}\right)}_{\mathit{wall}}& =& {\rho }_{\mathit{wall}}{\mathit{w}}_{\mathit{wall}}\\ {\mathit{e}}_{\mathit{wall}}& =& \begin{array}{c}{\rho }_{\mathit{wall}}\left[{\mathit{C}}_{\mathit{v}}{\mathit{T}}_{\mathit{wall}}\right.\\ \left.+\frac{1}{2}\left({\mathit{u}}_{\mathit{wall}}^2+{\mathit{v}}_{\mathit{wall}}^2+{\mathit{w}}_{\mathit{wall}}^2\right)\right].\end{array}\end{array}$

Then, the ghost cell conservative properties are calculated using an average between the left and right cells

(58) $\begin{array}{ccc}{\mathit{\rho }}_{\mathit{cr}}& =& 2{\mathit{\rho }}_{\mathit{wall}}-{\mathit{\rho }}_{\mathit{cl}}\\ \mathit{\rho }{\mathit{u}}_{\mathit{cr}}& =& 2\mathit{\rho }{\mathit{u}}_{\mathit{wall}}-\mathit{\rho }{\mathit{u}}_{\mathit{cl}}\\ \mathit{\rho }{\mathit{v}}_{\mathit{cr}}& =& 2\mathit{\rho }{\mathit{v}}_{\mathit{wall}}-\mathit{\rho }{\mathit{v}}_{\mathit{cl}}\\ \mathit{\rho }{\mathit{w}}_{\mathit{cr}}& =& 2\mathit{\rho }{\mathit{w}}_{\mathit{wall}}-\mathit{\rho }{\mathit{w}}_{\mathit{cl}}\\ {\mathit{e}}_{\mathit{cr}}& =& 2{\mathit{e}}_{\mathit{wall}}-{\mathit{e}}_{\mathit{cl}}.\end{array}$

IMPLICIT BOUNDARY CONDITIONS

Implicit boundary conditions are necessary in order to obtain a truly implicit time-marching method. The use of explicit boundary conditions can limitate substantially the stability of the numerical method in the marching procedure for the solution convergence.

BASIC IMPLICIT FORMULATION

A simplified form of the implicit equation, Eq. (42), is written in this section in order to detail the implementation of implicit boundary conditions for flux vector splitting schemes,

(59) $\begin{array}{c}\left[\frac{{\mathit{V}}_{\mathit{cl}}}{\Delta \mathit{t}}+\left({\mathit{A}}_{{\mathit{k}}^{+}}^{+}+{\mathit{B}}_{{\mathit{k}}^{+}}^{+}\right){\mathit{S}}_{\mathit{k}}\right]\Delta {\mathit{Q}}_{\mathit{cl}}\hfill \\ +[\left({\mathit{A}}_{{\mathit{k}}^{-}}^{-}-{\mathit{B}}_{{\mathit{k}}^{-}}^{-}\right){\mathit{S}}_{\mathit{k}}]\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k}}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}}\hfill \end{array}$

where the repeated k index, in the second term in left hand side of the equation, indicates summation over all the k faces of the control volume. The equation above is written only to present the relation between an internal cell, cl, and a boundary cell, cr, k. In the original formulation, Eq. (42), the cl-th cell has contributions from other faces, which may or may not be boundaries.

As discussed in the present work, the ghost cells hold different values for inviscid and viscous calculations. Hence, using the splitting definition, presented in section Inviscid Fluxes Calculation, Eq. (59) can be written as

(60) $\begin{array}{c}\left[\frac{{\mathit{V}}_{\mathit{cl}}}{▵\mathit{t}}+({\mathit{A}}_{\mathit{k}+}^{+}+{\mathit{B}}_{\mathit{k}+}^{+}){\mathit{S}}_{\mathit{k}}\right]▵{\mathit{Q}}_{\mathit{cl}}+\hfill \\ {\mathit{A}}_{{\mathit{k}}^{-}}^{-}{\mathit{S}}_{\mathit{k}}\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k},\mathit{inv}}-{\mathit{B}}_{{\mathit{k}}^{-}}^{-}{\mathit{S}}_{\mathit{k}}\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k},\mathit{visc}}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}}.\hfill \end{array}$

The contributions of the boundary face can be expressed in terms of the internal cell corrections as

(61) $\begin{array}{c}\Delta {\mathit{Q}}_{\mathit{cr},\mathit{k},\mathit{inv}}={\mathit{F}}_{k,\mathit{inv}}\Delta {\mathit{Q}}_{\mathrm{cl}}\hfill \\ \Delta {\mathit{Q}}_{\mathit{cr},\mathit{k},\mathit{visc}}={\mathit{F}}_{\mathit{k},\mathit{visc}}\Delta {\mathit{Q}}_{\mathit{cl}}·\hfill \end{array}$

Hence, Eq. (60), can be rewritten as

(62) $\begin{array}{c}\left[\frac{{\mathit{V}}_{\mathit{cl}}}{\Delta \mathit{t}}+({\mathit{A}}_{\mathit{k}+}^{+}+{\mathit{F}}_{\mathit{k},\mathit{inv}}{\mathit{A}}_{{\mathit{k}}^{-}}^{-}\right.\hfill \\ \left.-{\mathit{F}}_{\mathit{k},\mathit{visc}}{\mathit{B}}_{{\mathit{k}}^{-}}^{-}+{\mathit{B}}_{\mathit{k}+}^{+}){\mathit{S}}_{\mathit{k}}\right]\Delta {\mathit{Q}}_{\mathit{cl}}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}}·\hfill \end{array}$

The viscous Jacobians are calculated using primitive variables, then the corrections are set for the primitive variables and applied directly at the calculation of the viscous Jacobians.The matrix B⁺_k+ already includes the contribution from the boundary. The A⁺_k+, A⁻_k-, B⁺_k+ and B⁻_k- are presented in the work of Scalabrin (2007).

INVISCID WALL MATRICES FOR IMPLICIT BOUNDARY CONDITIONS

The matrix for an inviscid wall or for a symmetry boundary is the same matrix presented in Eq. (54),

(63) ${\mathit{F}}_{\mathit{k},\mathit{inv},\mathit{wall}}={ℝ}^{-1}\mathbf{W}ℝ.$

The matrices are applied to ΔQ for the implicit boundary condition according to Eq. (61).

VISCOUS WALL WITH SPECIFIED TEMPERATURE MATRIX FOR IMPLICIT BOUNDARY CONDITIONS

The viscous Jacobians are created using primitive variables. Hence, the implementation of implicit viscous matrices is performed using primitive variables. They are applied directly at the calculation of the Jacobians matrices.

The code was originally created for flow simulation with chemical reactions, hence the mass fraction, Y , is present in the primitive variable vector, which is given by

(64) $V=[\mathit{Y}\mathit{u}\mathit{v}\mathit{w}\mathit{T}{]}^{\mathit{T}}.$

In the present work, it is always considered that there is only one species in the flow, then, Y = 1.

The implicit _F matrix for wall boundary with specified temperature is derived from the average between the left and right cells,

(65) ${\mathit{z}}_{\mathit{cr}}=2{\mathit{z}}_{\mathit{wall}}-{\mathit{z}}_{\mathit{cl}},$

in which z is given property. The velocity components and temperature at the wall are considered constants in time, hence,

(66) $\begin{array}{c}\Delta {\mathit{u}}_{i_{\mathrm{wall}}}^n=0\\ \Delta {\mathit{t}}_{\mathrm{wall}}^n=0\end{array}.$

The mass fraction, Y, at the wall, is given by the left cell value, Y_wall = Y_cl Y. Therefore, the implicit matrix for wall boundary

^_Ywallcl

with specified temperature is given by.

(67) ${\mathit{F}}_{\mathit{k},\mathit{visc}}=\left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& -1& 0& 0& 0\\ 0& 0& -1& 0& 0\\ 0& 0& 0& -1& 0\\ 0& 0& 0& 0& -1\end{array}\right].$

CONSERVATION ANALYSIS OF IMPLICIT METHODS

GENERAL CONSERVATION ANALYSIS

After this brief review, the backward Euler time marching method, Eq. (42), is re-written using a simplified notation.

(68) $\left[\mathit{M}\right]\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}+\frac{{\mathit{V}}_{\mathit{cl}}}{\Delta \mathit{t}}[\mathit{I}]\Delta {\mathit{Q}}_{\mathit{cl}}^{\mathit{n}}={\mathit{R}}_{\mathit{cl}}^{\mathit{n}}$

The integration is applied on a generic grid, which is a representative 2-D finite volume mesh using nine quad-cells, as illustrated in Fig. 3, in order to present the conservative issues of such linearization. The [M] matrix, for the mesh used here, is represented by

Figure 3
Representative 2-D finite volume mesh.

(69) $\begin{array}{c}[\mathit{M}]=\left[\begin{array}{cccc}\sum _{{\mathit{A}}_{2,1}^{-}}{\mathit{A}}_1^{+}\hfill & \hfill & \hfill & {\mathit{A}}_{1,4}^{-}\hfill \\ {\mathit{A}}_{4,1}^{-}\hfill & \sum _{{\mathit{A}}_{3,2}^{-}}^{{\mathit{A}}_{1,2}^{-}}{\mathit{A}}_{2,3}^{+}\hfill & \hfill & \hfill \\ \hfill & {\mathit{A}}_{5,2}^{-}\hfill & \sum _{}^{{\mathit{A}}_{2,3}^{-}}{\mathit{A}}_3^{+}\hfill & \hfill \\ \hfill & \hfill & {\mathit{A}}_{6,3}^{-}\hfill & \sum _{{\mathit{A}}_{5,4}^{-}}{\mathit{A}}_4^{+}\hfill \\ \hfill & \hfill & \hfill & {\mathit{A}}_{7,4}^{-}\hfill \end{array}\right.\hfill \\ \left.\begin{array}{ccccc}{\mathit{A}}_{2,5}^{-}\hfill & \hfill & \hfill & \hfill & \hfill \\ \sum _{{\mathit{A}}_{6,5}^{-}}^{{\mathit{A}}_{4,5}^{-}}{\mathit{A}}_5^{+}\hfill & {\mathit{A}}_{3,6}^{-}\hfill & \hfill & \hfill & \hfill \\ \hfill & \sum _{}^{{\mathit{A}}_{5,6}^{-}}{\mathit{A}}_6^{+}\hfill & {\mathit{A}}_{4,7}^{-}\hfill & \hfill & \hfill \\ {\mathit{A}}_{8,5}^{-}\hfill & \hfill & \sum _{{\mathit{A}}_{8,7}^{-}}^{}{\mathit{A}}_7^{+}\hfill & {\mathit{A}}_{5,8}^{-}\hfill & \hfill \\ \hfill & {\mathit{A}}_{9,6}^{-}\hfill & \hfill & \sum _{{\mathit{A}}_{9,8}^{-}}^{{\mathit{A}}_{7,8}^{-}}{\mathit{A}}_8^{+}\hfill & {\mathit{A}}_{6,9}^{-}\hfill \\ \hfill & \hfill & \hfill & \hfill & \sum _{}^{{\mathit{A}}_{8,9}^{-}}{\mathit{A}}_9^{+}\hfill \end{array}\right]\hfill \end{array}$

The formulation can be considered conservative when the sum of V_i ΔQ_iⁿ terms, for all cells, yields zero, i.e.,

(70) $\sum [{\mathit{V}}_{\mathit{i}}\Delta {\mathit{Q}}_{\mathit{i}}^{\mathit{n}}]=0,$

or yet

(71) $\sum [{\mathit{V}}_{\mathit{i}}{\mathit{Q}}_{\mathit{i}}^{\mathit{n}}]=\sum [{\mathit{V}}_{\mathit{i}}{\mathit{Q}}_{\mathit{i}}^{\mathit{n}+1}].$

Solving the linear system for the converged solution, which means that R_iⁿ = 0, one obtains

(72) $\begin{array}{cc}\hfill \left[\frac{{\mathit{V}}_1}{\Delta \mathit{t}}+\sum {\mathit{A}}_1^{+}+{\mathit{A}}_{2,1}^{-}+{\mathit{A}}_{4.1}^{-}\right]& \hfill \Delta \mathit{Q}1+\\ \hfill \left[\frac{{\mathit{V}}_2}{\Delta \mathit{t}}+{\mathit{A}}_{1,2}^{-}+\sum {\mathit{A}}_2^{+}+{\mathit{A}}_{3.2}^{-}+{\mathit{A}}_{5,2}^{-}\right]& \hfill \Delta {\mathit{Q}}_2+\\ \hfill \left[\frac{{\mathit{V}}_3}{\Delta \mathit{t}}+{\mathit{A}}_{2,3}^{-}+\sum {\mathit{A}}_3^{+}+{\mathit{A}}_{6,3}^{-}\right]& \hfill \Delta {\mathit{Q}}_3+\\ \hfill \left[\frac{{\mathit{V}}_4}{\Delta \mathit{t}}+{\mathit{A}}_{1,4}^{-}+\sum {\mathit{A}}_4^{+}+{\mathit{A}}_{5,4}^{-}+{\mathit{A}}_{7,4}^{-}]\right]& \hfill \Delta {\mathit{Q}}_4+\\ \hfill \left[\frac{{\mathit{V}}_5}{\Delta \mathit{t}}+{\mathit{A}}_{2,5}^{-}+{\mathit{A}}_{4,5}^{-}\sum {\mathit{A}}_5^{+}+{\mathit{A}}_{6,5}^{-}+{\mathit{A}}_{8,5}^{-}\right]& \hfill \Delta {\mathit{Q}}_5+\\ \hfill \left[\frac{{\mathit{V}}_6}{\Delta \mathit{t}}+{\mathit{A}}_{3,6}^{-}+{\mathit{A}}_{5,6}^{-}\sum {\mathit{A}}_6^{+}+{\mathit{A}}_{9,6}^{-}\right]& \hfill \Delta {\mathit{Q}}_6+\\ \hfill \left[\frac{{\mathit{V}}_7}{\Delta \mathit{t}}+{\mathit{A}}_{4,7}^{-}+\sum {\mathit{A}}_7^{+}+{\mathit{A}}_{8,7}^{-}\right]& \hfill \Delta {\mathit{Q}}_7+\\ \hfill \left[\frac{{\mathit{V}}_8}{\Delta \mathit{t}}+{\mathit{A}}_{5,8}^{-}+{\mathit{A}}_{7,8}^{-}+\sum {\mathit{A}}_8^{+}+{\mathit{A}}_{9,8}^{-}\right]& \hfill \Delta {\mathit{Q}}_8+\\ \hfill \left[\frac{{\mathit{V}}_9}{\Delta \mathit{t}}+{\mathit{A}}_{6,9}^{-}+{\mathit{A}}_{8,9}^{-}+\sum {\mathit{A}}_9^{+}\right]& \hfill \Delta {\mathit{Q}}_9=0.\end{array}$

It is possible to simplify this equation in order to write

(73) $\left[\frac{{\mathit{V}}_1}{\Delta \mathit{t}}+[{\mathit{C}}_1\right]\Delta {\mathit{Q}}_1+\left[\frac{{\mathit{V}}_2}{\Delta \mathit{t}}+[{\mathit{C}}_2]\right]\Delta {\mathit{Q}}_2+\mathrm{...}=0,$

where the [C_i] matrix is the sum of the inviscid and viscous Jacobian matrices in each column of the [M] matrix. Comparing Eqs. (71) and (73), one can easily state that, in order to preserve the conservative property of the finite volume formulation, the [C_i] matrices should be zero and Δ_t should be the same for all the cells in the domain, i.e.,

(74) $\left[\frac{{\mathit{V}}_1}{\Delta \mathit{t}}\right]\Delta {\mathit{Q}}_1+\left[\frac{{\mathit{V}}_2}{\Delta \mathit{t}}\right]\Delta {\mathit{Q}}_{\mathit{2}}+\mathrm{...}=\frac{1}{\Delta \mathit{t}}\sum [{\mathit{V}}_{\mathit{i}}\Delta {\mathit{Q}}_{\mathit{i}}^{\mathit{n}}]=0.$

POINT IMPLICIT CONSERVATION ANALYSIS

If one solves the linear system for the same mesh presented in Fig. 3 using the point-implicit integration, after a converged solution is obtained, i.e., R_i = 0, one could finally write

(75) $\begin{array}{cc}\hfill \left\{\left[\frac{{\mathit{V}}_1}{\Delta \mathit{t}}+\sum {\mathit{A}}_1^{+}\right]\Delta {\mathit{Q}}_1^{\mathit{p}}+[{\mathit{A}}_{2,1}^{-}+{\mathit{A}}_{4,1}^{-}]\Delta {\mathit{Q}}_1^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_2}{\Delta \mathit{t}}+\sum {\mathit{A}}_2^{+}\right]\Delta {\mathit{Q}}_2^{\mathit{p}}+[{\mathit{A}}_{1,2}^{-}+{\mathit{A}}_{3.2}^{-}+{\mathit{A}}_{5,2}^{-}]\Delta {\mathit{Q}}_2^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_3}{\Delta \mathit{t}}+\sum {\mathit{A}}_3^{+}\right]\Delta {\mathit{Q}}_3^{\mathit{p}}+[{\mathit{A}}_{2.3}^{-}+{\mathit{A}}_{6,3}^{-}]\Delta {\mathit{Q}}_3^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_4}{\Delta \mathit{t}}+\sum {\mathit{A}}_4^{+}\right]\Delta {\mathit{Q}}_4^{\mathit{p}}+[{\mathit{A}}_{1,4}^{-}+{\mathit{A}}_{5,4}^{-}+{\mathit{A}}_{7,4}^{-}]\Delta {\mathit{Q}}_4^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_5}{\Delta \mathit{t}}+\sum {\mathit{A}}_5^{+}\right]\Delta {\mathit{Q}}_5^{\mathit{p}}+[{\mathit{A}}_{2,5}^{-}+{\mathit{A}}_{4,3}^{-}+{\mathit{A}}_{6,5}^{-}+{\mathit{A}}_{8,5}^{-}]\Delta {\mathit{Q}}_5^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_6}{\Delta \mathit{t}}+\sum {\mathit{A}}_6^{+}\right]\Delta {\mathit{Q}}_6^{\mathit{p}}+[{\mathit{A}}_{3,6}^{-}+{\mathit{A}}_{5,6}^{-}+{\mathit{A}}_{9,6}^{-}]\Delta {\mathit{Q}}_6^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_7}{\Delta \mathit{t}}+\sum {\mathit{A}}_7^{+}\right]\Delta {\mathit{Q}}_7^{\mathit{p}}+[{\mathit{A}}_{4.7}^{-}+{\mathit{A}}_{8,7}^{-}]\Delta {\mathit{Q}}_7^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_8}{\Delta \mathit{t}}+\sum {\mathit{A}}_8^{+}\right]\Delta {\mathit{Q}}_8^{\mathit{p}}+[{\mathit{A}}_{5,8}^{-}+{\mathit{A}}_{7,8}^{-}+{\mathit{A}}_{9,8}^{-}]\Delta {\mathit{Q}}_8^{\mathit{p}-1}\right\}& \hfill +\\ \hfill \left\{\left[\frac{{\mathit{V}}_9}{\Delta \mathit{t}}+\sum {\mathit{A}}_9^{+}\right]\Delta {\mathit{Q}}_9^{\mathit{p}}+[{\mathit{A}}_{6,9}^{-}+{\mathit{A}}_{8,9}^{-}]\Delta {\mathit{Q}}_9^{\mathit{p}-1}\right\}=& \hfill 0.\end{array}$

One can state that, as long the sub-iterations of the point implicit method have not achieved the convergence for ΔQ_i^p and/or Δt is not the same for all the cells in the domain, it is not possible to find a general relation for any Δ Q_i^p and Δ Q_i^p−1 such that ∑ [V_i ΔQ_iⁿ] = 0. Therefore, the conservation property for the point-implicit integration can be achieved only if the [C_i] matrices are zero, the solution of ΔQ_i^p for the sub-iterations in p is converged, and Δt is constant over the entire mesh.

However, achieving convergence of the point-implicit subiterations can be as expensive as performing the fully implicit integration. Hence, all practical numerical solvers perform a limited number of sub-iterations. The authors, usually, perform up to 10 sub-iterations in p for the point-implicit integration and, then, move on to the next time step. Typically, this is not enough to achieve convergence for ΔQ_i^p, as discussed the forthcoming section of the paper.

RESULTS AND DISCUSSION

Results for the so-called "rigid body simulation" problem are presented in this section in order to expose the effects of the time-marching scheme on the mass conservation. The rigid body problem consists of the simulation of the flow contained between two concentric cylinders, in which both walls rotate at the same angular velocity. Therefore, after convergence, the fluid in the domain is rotating at the same angular velocity as if it were a rigid body. Here, the problem is addressed as a 2-D flow. The present work performed simulations using the point-implicit and the Runge-Kutta time marching methods. An analysis of the use of a constant CFL number on all mesh cells is also presented here, in comparison to the use of a constant Δt throughout the mesh.

The two dimensional geometry and the detailed 360,000 cell mesh, used in the simulations, are presented in Fig. 4. The external and internal cylinders are rotating walls with fixed angular velocity, ω, and fixed temperature, T₀. Air with zero velocity and T₀ temperature are considered as initial conditions. All simulations are conducted with the same initial and boundary conditions. Each simulation is performed using a different time marching method, as presented in Table 1.

Figure 4
Rigid body geometry and mesh detail.

Table 1
Defi nition of the test case confi gurations for the numerical simulations.

Case	Time marching scheme	Constant Δt or CFL
1	Point implicit	Constant CFL
2	Point implicit	Constant Δt
3	Runge-Kutta	Constant CFL
4	Runge-Kutta	Constant Δt

Figures 5 and 6 present the total amount of mass (per unit of length in the cylinder axial direction) inside the computational domain, as a function of the iteration number, for the rigid body simulations. It is clear from the figures that the point implicit time marching method can significantly deteriorate the important conservative property of the finite volume formulation. Moreover, both simulations with the point-implicit timemarching scheme presented exactly the same non-conservative behavior. Moreover, for the point-implicit time integration test cases, i.e., test cases 1 and 2, or cases shown in Fig. 5 (a) and (b), there is almost no influence of the selection of constant CFL or of constant Δt in the time march. In other words, one could state that, for these test cases, the non-conservation effects of using a constant CFL number are far less significant than the effects of using the point implicit integration. Furthermore, it is correct to state that both simulations diverged after some time (not shown in Fig. 5 (a) and (b).

Figure 5
Effects of implicit time marching scheme on total mass conservation; (a) Constant CFL calculation for point-implicit scheme; (b) Constant Δt calculation for point-implicit scheme.

Figure 6
Effects of explicit time marching scheme on total mass conservation; (a) Constant CFL calculation for explicit RK scheme; (b) Constant Δt calculation for explicit RK scheme.

In contrast to that behavior, simulations performed using the explicit Runge-Kutta scheme and a constant time step throughout the domain perfectly conserve the mass in the computational domain, as one would expect from a finite volume code. The results for this test case (case 4) are shown in Fig. 6 (b). On the other hand, results in Fig.6 (a) indicate that, even with an explicit scheme, there is no mass conservation if a variable time step, or constant CFL number, is used in the time integration. This is a serious problem since most convergence acceleration procedures typically employed in aerospace CFD codes are based on the use of implicit integration or on variable time stepping, or both. The present results are clearly demonstrating that, for such cases, there is no mass conservation during the transient process of converging to a steady state solution.

Moreover, all results presented in Figs. 5 and 6 reinforce the previous analysis performed in this work. To assure the conservative property of a time marching method, the [C_i] matrix should be zero. This is automatically enforced by the explicit integration schemes by construction. Therefore, for the point implicit integration, the convergence of the sub-iterations has to be achieved in order to obtain a conservative scheme. Furthermore, all the mesh cells have to advance in time using the same Δt value in order to maintain the conservative property of the finite volume method. This is true even for the explicit time marching methods.

CONCLUSIONS

A discussion on the issues associated with the coupling of implicit integration methods, for unstructured finite volume formulations, with the spatial discretization based on flux vector splitting schemes, is presented in this work. The linearization of inviscid and viscous Jacobians may result in a non-conservative method during the transient phase of the flow simulation, even for the finite volume formulation, which is supposed to be conservative by construction. The present analysis of the numerical formulation and of the computational results obtained has indicated that the time integration can be considered conservative only if the sum of the Jacobian matrices in each column of the linear system matrix is zero and all mesh cells have the same Δt value. Moreover, very popular approximate methods used in many CFD codes, to solve sparse linear systems, such as the point implicit integration, need to achieve convergence of the sub-iterations in order to be conservative during the transient portion of the simulation. This is a very expensive proposition and it can make such approximate solvers as expensive as those which perform the direct solution of the full implicit linear system.

The important conclusion is that care must be exercised in the linearization of time-marching methods for simulations which demand the conservative property. It is important to be aware of the effects of such issue on the physical problem of interest. Physical problems with incoming and outcoming flow boundary conditions, very common in simulations for aerospace applications, do not necessarily need a conservative scheme during the transient portion of a steady state calculation. The amount of variation in the flow properties, during one time step, is negligible compared to the flux of the same properties that is crossing the open boundaries of the domain. Moverover, for many steady-state external aerodynamic applications, the initial conditions are strictly numerical, in the sense that they cannot be physically realized as implemented in the solver. In other words, they are typically associated to impulsively started flows. For such cases, the aspect of the lack of mass conservation in the entire computational domain, as the solution evolves from a non-physical initial condition to a physically relevant converged steady state condition, is not an issue. On the other hand, the conservative property is absolutely essential for closed systems, in which the total mass, and other properties, must always be conserved in order to achieve a physically relevant solution. In these critical cases, approximate numerical methods should not be used in order to solve the full implicit linear systems.

ABBREVIATIONS

• 2-D:   Two dimensional
• 3-D:   Three dimensional
• CFD:   Computational Fluid Dynamics
• CFL:   Courant-Friedrichs-Lewy number
• LeMANS:   Le Michigan Aerothermodynamics Navier- Stokes Solver
• RANS:   Reynolds Averaged Navier-Stokes
• SA:   Spalart-Allmaras Tubulence Model
• SW:   Sterger-Warming Scheme
• TVD:   Total Variation Diminishing

LIST OF SYMBOLS

English Characters

• a:  Speed of the sound
• A:   Jacobian matrix of the inviscid flux
• B:  Jacobian matrix of the viscous flux
• C:  Sum of Jacobian matrices
• C_p:  Specific heat at constant pressure
• C_v:  Specific heat at constant volume
• d_k:   Distance of k-th face to the nearest wall
• e:   Total energy per unit volume
• e_i:   Internal energy
• F_e:   Inviscid flux vector
• F_v:   Viscous flux vector
• î_x, î_y, î_z:   Cartesian-coordinate orthonormal vector basis
• I:   Identity matrix
• L:   Matrix of eigenvectors from the left
• m:   Normal vector to the nearest wall
• n:   Normal vector
• nf:   Number of faces of a given volume
• p:   Static pressure
• Pr:   Prandtl number
• qH:   Heat transfert vector
• Q:   Conserved variable vector
• R:   Matrix of eigenvectors from the right
• R:   Residue
• ℝ:  Rotation matrix
• ℜ:   Gas constant
• S:  Outward-oriented area vector
• S:   Constant for the Sutherland law equation
• t:  Time
• T:   Temperature
• v = {u, v, w}:   Velocity vector in the cartesian coordinate
• V:   Volume
• V:   Primitive variables vector
• W:   Inviscid wall matrix
• w:   Switch of the Steger and Warming scheme
• Y:   Mass fraction

Greek Characters

• α:   Switch factor
• β:   Viscous force work and heat transfer term
• δij:   Kronecker delta
• F:   Implicit matrix
• ∈:   Artificial dissipation
• γ:   Ratio of specific heats
• κ:   Thermal conductivity coefficient
• λ:   Eigenvalue of the Jacobian matrix
• ∧:   Diagonal matrix of the eigenvalues of the Jacobian
• µ:   Dynamic molecular viscosity coefficient
• ν:   Kinematic molecular viscosity coefficient
• ρ:   Density
• τ_ij:   Shea-stress tensor

Subscripts

• cl:   Cell on the left side of the fac
• cr:   Cell on the right side of the face
• ∞:   Free-stream property
• inv:   Inviscid property
• k:   Index of the cell face
• n:   Normal property at a given face
• wall:   Property at the wall
• visc:   Viscous property

Superscripts

• n:   Index of iteration in time
• p:   Index of iteration for the point-implicit time integration
• +:   Positive part of a matrix or vector
• -:   Negative part of a matrix or vector
• rot:   Rotated property

ACKNOWLEDGEMENTS

The authors would like to acknowledge Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq, which partially supported the project under the Research Grant No. 309985/2013-7. Partial support for the present research was also provided by Fundação de Amparo à Pesquisa do Estado de São Paulo, FAPESP, under Research Grants No. 2013/07375 and No. 2013/21535-0. Further partial support was provided by Fundação Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, CAPES, through a Ph.D. scholarship to the first author, which is also gratefully acknowledged.

REFERENCES

• Anderson, J.D.J., 1991, "Fundamentals of Aerodynamics", McGraw- Hill International Editions, New York, NY, USA.
• Anderson, W.K., Thomas, J.L. and Van Leer, B., 1986, "A Comparison of Finite Volume Flux Vector Splittings for the Euler Equations", AIAA Journal, Vol. 24, No. 6, pp. 1453-1460.
• Barth, T.J., 1987, "Analysis of Implicit Local Linearization Techniques for Upwind and TVD Algorithms", Proceedings of the 25th AIAA Aerospace Sciences Meeting, AIAA Paper No. 87-0595, Reno, NV, USA.
• Gnoffo, P.A., 2003, "Computational Aerothermodynamics in Aeroassist Applications", Journal of Spacecraft and Rockets, Vol. 40, No. 3, pp. 305-312. doi: 10.2514/2.3957.
  » https://doi.org/10.2514/2.3957
• Harten, A., 1983, "High Resolution Schemes for Hyperbolic Conservation Laws", Journal of Computational Physics, Vol. 49, No. 3, pp. 357-393. doi: 10.1016/0021-9991(83)90136-5.
  » https://doi.org/10.1016/0021-9991(83)90136-5
• Hirsch, C., 1990, "Numerical Computation of Internal and External Flows", Vol. 2, Computational Methods for Inviscid and Viscous Flows, John Wiley and Sons, New York, USA.
• Junqueira-Junior, C.A., 2012, "A Study on the Extension of an Upwind Parallel Solver for Turbulent Flow Applications" Masters Thesis, Instituto Tecnológico de Aeronáutica, São José dos Campos, SP, Brasil.
• Junqueira-Junior, C.A., Azevedo, J.L.F., Scalabrin, L.C. and Basso, E., 2013, "A Study of Physical and Numerical Effects of Dissipation on Turbulent Flows Simulations", Journal of Aerospace Technology and Management, Vol. 5, No. 2, pp. 145-168. doi: 10.5028/jatm.v5i2.179.
  » https://doi.org/10.5028/jatm.v5i2.179
• Lomax, H., Pulliam, T.H. and Zingg, D.W., 2001, "Fundamentals of Computation Fluid Dynamics", Springer, NY, USA.
• MacCormack, R.W. and Candler, G.V., 1989, "The Solution of the Navier-Stokes Equations Using Gauss-Seidel Line Relaxation", Computer and Fluids, Vol. 17, No. 1, pp. 135-150. doi: 10.1016/0045-7930(89)90012-1.
  » https://doi.org/10.1016/0045-7930(89)90012-1
• Ramalho, M.V.C., Azevedo, J.L.F. and Azevedo, J.H., 2011, "Further Investigation into the Origin of the Carbuncle Phenomenon in Aerodynamic Simulations", 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, AIAA Paper No. 2011-1184, Orlando, FL, USA.
• Scalabrin, L.C., 2007, "Numerical Simulation of Weakly Ionized Hypersonic Flow Over Reentry Capsules", PhD Thesis, Department of Aerospace Engineering, The University of Michigan, Michigan, USA.
• Schlichting, H., 1978, "Boundary-Layer Theory", 7th Edition, McGraw Hill, New York, USA.
• Spalart, P.R. and Allmaras, S.R., 1992, "A One-Equation Turbulence Model for Aerodynamic Flows", Proceedings of the 30th AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper No. 92-0439, Reno, NV, USA.
• Spalart, P.R. and Allmaras, S.R., 1994, "A One-Equation Turbulence Model for Aerodynamic Flows", La Recherche Aerospatiale, Vol. 1, pp. 5-21.
• Steger, J.L. and Warming, R.F., 1981, "Flux Vector Splitting of the Inviscid Gasdynamic Equations with Application to the Finite-Difference Method", Journal of Computational Physics, Vol. 40, No. 2, pp. 263-293. doi: 10.1016/0021-9991(81)90210-2.
  » https://doi.org/10.1016/0021-9991(81)90210-2
• Venkatakrishnan, V., 1995, "Implicit Schemes and Parallel Computing in Unstructured Grid CFD", Technical Report, ICASE Report No. 95-28, ICASE.
• Wright, M.J., 1997, "A Family of Data Parallel Relaxation Methods for the Navier-Stokes Equations", PhD Thesis, University of Minnesota, Minneapolis, MN, USA.

Publication Dates

• Publication in this collection
  Jul-Sep 2014

History

• Received
  24 Apr 2014
• Accepted
  19 Aug 2014