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

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


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
Th e 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.Th e RANS equations are obtained by time fi ltering the Navier-Stokes set of equations.Th e compressible RANS equations are written in the algebraic vector form as .
(1) Th e conserved variables vector, Q, the inviscid fl ux vector, F e , and the viscous fl ux vector, F v , are given by , ( 2) in which ρ stands for density, for the velocity vector in Cartesian coordinates, p for static pressure, τ for viscous stress tensor, q H for heat fl ux vector, e for the total energy per unit volume and β i is given by . (5) Th e î x , î y and î z terms are the Cartesian-coordinate orthonormal vector basis.It is very important to emphasize that fi eld forces, such as gravity, are neglected here.
Other equations are necessary in order to close the system of equations given by Eq. ( 4).Th ese additional equations are called constitutive relations.Th e fi rst constitutive equation presented to close the Navier-Stokes set of equations is known as the equation of state.Th is equation considers the perfect gas law, and it is written as in which the mean total energy per unit volume, e, is given by , ( 7) and e i stands for internal energy, defi ned as in which T stands for the mean static temperature and C v is the specifi c heat at constant volume.Th e heat fl ux from Eq. ( 5) is obtained from the Fourier law for heat conduction, and it is given by , in which γ is the ratio of specifi c 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 specifi c heat at constant pressure and μ is the dynamic molecular viscosity coeffi cient, calculated as a function of the temperature by the Sutherland law equation (Anderson, 1991), written as where S = 110K , and µ ∞ is the dynamic molecular viscosity coeffi cient of the fl uid at temperature T ∞ .Th e components of the viscous stress tensor, for a Newtonian fl uid, are given by , (11) in which δ ij stands for the Kronecker delta.

NUMERICAL FORMULATION
Th e numerical formulation applied in this work is briefl y presented in this section.Th e 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).Th e original framework was developed by Scalabrin (2007) at the University of Michigan to simulate hypersonic fl ows over reentry capsules.It is a three dimensional (3-D) numerical solver that uses the fi nite 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).Th e numerical fl ux is calculated using an upwind scheme based on the Steger-Warming fl ux vector splitting (Steger and Warming, 1981).Th e time integration is performed by point implicit and Runge-Kutta methods.Th is 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 diff erential equations.It is applied by the CFD community to fi nd the solution of the RANS, equations.Th e method is obtained integrating the fl ow equations for each control volume within a given mesh.Th e RANS equations, written in the context of a cell-centered fi nite-volume formulation, are given by . ( Considering a cell-centered formulation, V i is a given cell of the given grid.Aft er the integration it is possible to apply the Gauss theorem over the equation above, resulting in: in which S i is the outward-oriented area vector and it is defi ned as Considering the mean value of the conserved variables within the i-th control volume, one can write the fi rst term of Eq. (13) as .( 15) Th e second term of Eq. ( 13) can be written as the sum of all faces of a cell , ( 16) 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 fi nite volume approximation are given by . ( For this formulation, the fl uxes are computed at the faces of the control volume, and the conserved variables are computed in the cell.

INVISCID FLUX CALCULATION
Th e inviscid fl uxes are calculated using a method based on a classical fl ux vector splitting formulation, the Steger-Warming Scheme (SW) (Steger and Warming, 1981).Th is method is an upwind scheme which uses the homogeneous property of the inviscid fl ux vectors to write where Fe n is the normal flux at the k-th face, and A is the Jacobian matrix of the inviscid fl ux which can be diagonalized by the matrices of its eigen vectors from the left and from the right, L and R and Λ is the diagonal matrix of the eigenvalues of the Jacobian matrix.Th e A matrix can be split into positive and negative parts as Th e splitting separates the fl ux into two parts, the downstream and the upstream fl ux, in relation to the face orientation as: where the cl and cr subscripts are the cells on the left and right sides of the face.Th e split eigen values of the Jacobian matrix are given by , In order to avoid sudden sign transition, and then discontinous derivative, the split eigen values receive a small number, ε , turning Eq. ( 22) into .( 23) Th e soft en 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: in which ( 25) Th e switch, w, is given by where ∇p is a scalar number, a numerical approximation of the pressure gradient.For small ∇p, = (1 − ) = 0.5, the code runs with a centered scheme, and for larger values of ∇p, = 0 and (1 − ) = 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).Th e applied formulation was originally created with interest on studying fl ows over reentry capsules.For such particular cases, with very strong shock waves, it is very common to fi nd solutions with numerical and non-physical structures such as carbuncles (Ramalho et al., 2011).To avoid such numerical problems, artifi cial dissipation has to be added to the method.Th e dissipation was included into the split eigenvalues, Eq. ( 23), using an ε factor, which is given by: ( 27) where d k is the distance of the k-th face, to the nearest wall boundary.d 0 is set by the user and must be smaller than the boundary layer thickness and larger than the shock stand-off distance.
is the normal vector of the nearest wall, and is the normal vector to the k-th face.Equation ( 27) applies the term 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
Th e 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.Th is computation is used to fi nd 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.Th e properties at the faces are calculated using simple averages.For example, in Figs. 1 and 2, the properties are given by ( 28) . (29) Considering ∇Q as a constant over the cell, the equation above yields 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 The derivative computation for other types of element, 2-D or 3-D, is straightforward.

TIME INTEGRATION
Th e 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.
Th e second-order Runge-Kutta integration method used in this work is given by (32) In order to simplify the forthcoming equations, the right hand side of Eq. ( 32) is written as (33) in which R cl is the residue of the i-th cell.
Th e implicit integration applied in the present work is based on the backward Euler method, which is given by .( 34  From the spatial discretization, the inviscid terms can be written as It is a common practice to assume and then write (38) 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).Th e viscous terms can be written in the same form as 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.Th e 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 (Gnoff o, 2003;Venkatakrishnan, 1995;Wright, 1997).
Th e 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., . ( It is assumed that ∆Q n+1.0 = 0 and four iterations are taken in the process as suggested in the literature (Wright, 1997).Th e given sparse linear system illustrates the point where each ☐ is a 5 × 5 block matrix.Th e time step is computed by 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
Th e boundary conditions are implemented using ghost cells.Th e solver creates the ghost cells to hold properties that satisfy the correct fl ux calculation at the boundaries.Th e implementation assigns properties that satisfy the Euler boundary conditions to calculate the inviscid fl uxes, and properties that satisfy the Navier-Stokes boundary conditions to calculate the viscous fl uxes.Th erefore, the ghost volumes store two diff erent types of fl uxes for the correct computation of the RANS equations.

WALL INVISCID BOUNDARY CONDITIONS
Th e ghost cells hold the properties in the same manner to calculate the inviscid fl uxes at the wall and at the symmetry boundaries.Mass and energy fl uxes 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 defi ned as the properties at the left side of the boundary face and the ghost cells are defi ned 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) Th e vector of conservative properties, Q, is written as (49) In Eq. ( 48), is the rotation matrix given by,

VISCOUS WALL BOUNDARY CONDITION WITH SPECIFIED TEMPERATURE
Viscous wall with specifi ed 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) One can computate the conservative properties at the wall boundary using the Cartesian components of the wall velocity, and , which are set by the user. (57) Th en, the ghost cell conservative properties are calculated using an average between the left and right cells (58)

IMPLICIT BOUNDARY CONDITIONS
Implicit boundary conditions are necessary in order to obtain a truly implicit time-marching method.Th e 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 fl ux vector splitting schemes, (59) 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.Th e 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 diff erent values for inviscid and viscous calculations.Hence, using the splitting defi nition, presented in section Inviscid Fluxes Calculation, Eq. ( 59) can be written as

INVISCID WALL MATRICES FOR IMPLICIT BOUNDARY CONDITIONS
Th e matrix for an inviscid wall or for a symmetry boundary is the same matrix presented in Eq. ( 54), k, inv, wall = -1 W . (63) 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
Th e viscous Jacobians are created using primitive variables.Hence, the implementation of implicit viscous matrices is performed using primitive variables.Th ey are applied directly at the calculation of the Jacobians matrices.
Th e code was originally created for fl ow simulation with chemical reactions, hence the mass fraction, Y , is present in the primitive variable vector, which is given by . (64) In the present work, it is always considered that there is only one species in the fl ow, then, Y = 1.Th e implicit matrix for wall boundary with specifi ed temperature is derived from the average between the left and right cells,

GENERAL CONSERVATION ANALYSIS
Aft er this brief review, the backward Euler time marching method, Eq. ( 42), is re-written using a simplifi ed notation . ( The integration is applied on a generic grid, which is a representative 2-D fi nite volume mesh using nine quad-cells, as illustrated in Fig. 3, in order to present the conservative issues of such linearization.Th e [M] matrix, for the mesh used here, is represented by ( 69) Th e formulation can be considered conservative when the sum of V i ∆Q i n terms, for all cells, yields zero, i.e., ,  71) and ( 73), one can easily state that, in order to preserve the conservative property of the fi nite volume formulation, the [C i ] matrices should be zero and ∆ t should be the same for all the cells in the domain, i.e., .

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 n = 0, one could finally write (75) 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 fi nd a general relation for any Th erefore, 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.Th e 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 eff ects of the time-marching scheme on the mass conservation.Th e rigid body problem consists of the simulation of the fl ow contained between two concentric cylinders, in which both walls rotate at the same angular velocity.Th erefore, aft er convergence, the fl uid 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 fl ow.Th e 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.
Th e two dimensional geometry and the detailed 360,000 cell mesh, used in the simulations, are presented in Fig. 4. Th e external and internal cylinders are rotating walls with fi xed angular velocity, ω, and fi xed temperature, T 0 .Air with zero velocity and T 0 temperature are considered as initial conditions.All simulations are conducted with the same initial and boundary conditions.Each simulation is performed using a diff erent time marching method, as presented in Table 1.
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 fi gures that the point implicit time marching method can signifi cantly deteriorate the important conservative property of the fi nite 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 infl uence 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 eff ects of using a constant CFL number are far less signifi cant than the eff ects of using the point implicit integration.Furthermore, it is correct to state that both simulations diverged aft er some time (not shown in Fig. 5 (a) and (b).
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 fi nite volume code.Th e 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.Th is 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.Th e 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.Th is is automatically enforced by the explicit integration schemes by construction.Th erefore, 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    Study of Conservation on Implicit Techniques for Unstructured Finite Volume Navier-Stokes Solvers 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.

Figure 1
Figure 1.2-D example of new volume creation.Figure 2. 2-D example of derivative calculation.
Figure 1.2-D example of new volume creation.Figure 2. 2-D example of derivative calculation.
this work, the viscous Jacobian matrices are represented by B. Hence, the Jacobian matrix splitting is written as .
properties to the Cartesian coordinate frame.
of the boundary face can be expressed in terms of the internal cell corrections as (61) Hence, Eq. (60), can be rewritten as (62) Th e 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.Th e matrix B + k+ already includes the contribution from the boundary.Th e A + k+ , A − k-, B + k+ and B − k-are presented in the work of Scalabrin (2007).
z is given property.Th e velocity components and temperature at the wall are considered constants in time, hence, .(66) Th e mass fraction, Y, at the wall, is given by the left cell value, Y wall = Y cl .Th erefore, the implicit matrix for wall boundary with specifi ed temperature is given by .(67) the linear system for the converged solution, which means that R i n = 0, one obtains (72) It is possible to simplify this equation in order to write , (73) where the [C i ] matrix is the sum of the inviscid and viscous Jacobian matrices in each column of the [M] matrix.Comparing Eqs. (

Figure 4 .
Figure 4. Rigid body geometry and mesh detail.

Figure 5 .
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 .
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.

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