Accuracy analysis of unstructured finite volume discretization schemes for diffusive fluxes

Computers & Fluids 101 (2014) 220–232

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Accuracy analysis of unstructured finite volume discretization schemes for diffusive fluxes Alireza Jalali ⇑, Mahkame Sharbatdar, Carl Ollivier-Gooch Department of Mechanical Engineering, The University of British Columbia, 2054-6250 Applied Science Lane, Vancouver, BC V6T 1Z4, Canada

a r t i c l e

i n f o

Article history: Received 11 February 2014 Received in revised form 29 May 2014 Accepted 5 June 2014 Available online 18 June 2014 Keywords: Finite volume methods Diffusive flux Unstructured meshes Truncation error Discretization error

a b s t r a c t Numerical errors form a major source of uncertainty in CFD simulations. This source of error, which is more significant on unstructured meshes, can be alleviated by a careful choice of discretization scheme. In the context of finite volume methods, there are many suggestions for the discretization of diffusive fluxes appearing in viscous flow simulations. In this paper, a wide range of discretization schemes commonly used for diffusive flux approximation in cell-centered unstructured finite volume solvers are compared in terms of discretization and truncation errors. In addition, a novel eigenanalysis tool is proposed to relate these two forms of numerical error to each other and to interpret the error behaviors obtained by each scheme. The error comparisons are performed on unstructured meshes consisting of both isotropic and anisotropic triangles. Our results suggest that adding a solution jump term to the baseline face gradient (determined as the average of two adjacent cell gradients) provides the most accurate approximation for diffusive fluxes. Also, this term produces sufficient damping to stabilize the discretization even on highly skewed meshes. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The design of aircraft, from commercial and military transport to combat airplanes, depends increasingly on the use of computational fluid dynamics. In particular, computational aerodynamics, which is concerned with the accurate prediction of airplanes’ force and moment coefficients, is progressively being used as a practical tool in the design and optimization of aerodynamic objects. The numerical simulation of fluid flows is highly affected by the presence of three error sources: geometry modeling, physical modeling and numerical errors. Geometry modeling error is unavoidable in those cases where complex geometry can not be accurately modeled. In those circumstances, some minor geometrical features of objects are ignored to facilitate the numerical procedure. As an example, the wing root fillet is often neglected in the CFD simulation of fluid flow around airplanes. Recent studies have proved that numerical errors are at least as large as physical modeling errors for computations of the flow around transport aircraft. Mavriplis [1] carried out a grid convergence and sensitivity study on a wing-body configuration ⇑ Corresponding author. Tel.: +1 6046002324. E-mail addresses: [email protected] (A. Jalali), mahkame@interchange. ubc.ca (M. Sharbatdar), [email protected] (C. Ollivier-Gooch). http://dx.doi.org/10.1016/j.compfluid.2014.06.008 0045-7930/Ó 2014 Elsevier Ltd. All rights reserved.

which was the subject of the second AIAA Drag Prediction Workshop (DPW). Physical modeling errors were quantified based on the sensitivity of computed drag coefficient to the formulation of the viscous terms and to the turbulence model. The results revealed that using the thin-layer approximation near the wall rather than solving the full Navier–Stokes equations affects the drag coefficient by less than 2%. Likewise, the sensitivity to the turbulence modeling was examined by varying the way by which the distance from wall is obtained. The variations in finding this quantity, which is important in Spalart–Allmaras turbulence model [2], rarely influenced the drag coefficient. On the other hand, using topologically different grids with equal resolutions near the trailing edge substantially changed the solution. These observations reinforced the notion that numerical errors are still the dominant source of error in most aerodynamic simulations especially for unstructured flow solvers. In fact, mesh features – including cell size, anisotropy, shape, connectivity, and variations between adjacent cells – can have an adverse interaction with discretization schemes and this interaction may affect the solution accuracy. One way to control the accuracy of spatial discretization is improving the local quality of cells in an unstructured mesh by node displacement [3]. Although this approach is effective in reducing the local error in the solution, it negates the advantages of automatic unstructured mesh generation. However, the need for precise vertex movement can be alleviated provided that

221

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

careful attention is given to the discretization schemes used to solve the governing equations. In the context of finite volume methods, discretization schemes are divided into two major categories: convective and diffusive fluxes. Convective fluxes are those that depend only on the solution at the cell interfaces while diffusive fluxes are dependent upon the interface solution gradient. While there is a clear consensus over the use of flux difference splittings for convective fluxes [4], the choice of diffusive flux discretization and its effect on the accuracy of the solution is less understood. Diskin et al. [5] compared the accuracy of a limited number of the node-centered and cell-centered discretization schemes for viscous (diffusive) fluxes with the aim of improving turbulent simulations. This comparison was done on a range of grids from regular to irregular composed of arbitrary mixtures of triangles and quadrilaterals. The comparison showed that there is little difference in accuracy between node-centered schemes and the best cell-centered schemes, but poorly-designed cell-centered schemes behave much worse. They also compared the nodecentered and cell-centered unstructured finite volume discretization schemes for inviscid (convective) fluxes [6]. The schemes, being different in the least-squares gradient reconstruction, were compared for a constant coefficient convection equation, linear advection. The authors concluded that carefully-designed cell-centered schemes offer the best options for second-order discretization. Nishikawa [7] introduced a general principle for constructing diffusion schemes which is applicable to various discretization methods, including finite volume, residual distribution, discontinuous Galerkin, and spectral-volume methods. This principle is derived based on a hyperbolic relaxation-system model for diffusion problems and results in a damping term, which is essential for effective high-frequency error damping and, in some cases, for consistency. He also demonstrated for the first time that the diffusion schemes with a finite difference term are special cases of the new damping scheme in the context of finite volume discretizations [8]. Moreover, he extended the damping scheme to the Navier–Stokes equations [9] with an optimal value of the damping coefficient which has already been obtained by Fourier analysis and truncation error analysis for scalar problems. In this paper, we compare the accuracy of common secondorder-accurate discretization schemes for calculating diffusive fluxes on unstructured meshes. In particular, we consider the Poisson equation as a model of viscous discretization and we use variations on the cell-centered finite volume approach. We describe in Section 2 a wide range of discretization schemes which differ in how cell gradients are combined to compute a face gradient; in the presence and form of finite difference correction terms to the gradient; and in treatment of the discontinuity at faces. The initial source of numerical errors is the truncation error, defined as the difference between the continuous PDEs and the finite discretized equations. Truncation error plays an important role in error quantification since the truncation error can be used to estimate the discretization error that occurs during the approximate numerical solution of differential equations. Discretization error is defined as the difference between the exact solution of governing equations and the numerically approximated solution. Roy [10] showed that the discretization error can be directly related to the truncation error by the error transport equation where the truncation error serves as a source term to the linearized system of equation whose solution is the discretization error. For a linear problem such as Poisson’s equation, it is easy to show that the linear operator applied to the discretization error in the error equation is equivalent to the flux Jacobian. In Section 3, an eigendecomposition paradigm will be introduced where the truncation error is expanded as a combination of linearly independent right eigenvectors of the flux Jacobian matrix. The weights

corresponding to the right eigenvectors differ by the choice of discretization scheme and influence the error obtained in the numerical solution. The distribution of these weights are helpful in explaining how the change of discretization scheme results in a smaller or larger discretization error. In Section 4, the method of manufactured solution is used to compare the discretization and truncation errors produced by different schemes. The eigenanalysis tool is utilized to explain how the choice of schemes affects the eigenvalue spectra and weights associated with the eigendecomposition of error forms. It turns out that each class of discretization method alters these measures in a unique pattern that can be interpreted as the method signature. The accuracy analysis tests are performed on both isotropic and anisotropic unstructured triangular grids with known properties. 2. Discretization scheme To discretize the flow equations using the finite volume method, the governing equations should be recast in fully conservative form as

~ @U F¼S þ r ~ @t

ð1Þ

in which U denotes the conservative solution vector, F is the flux vector and S represents the source term vector. Assuming the discretized physical domain does not change in time and integrating Eq. (1) over an arbitrary control volume, the 2D finite volume formulation of the governing equations can be written in the form of a surface and an area integral:

dU i 1 ¼ ACV i dt

I

1 ~ ^ ds þ Fn ACV i CSi

Z

~  S dA ¼ R U

! ð2Þ

Ai

The right hand side is called the residual which is comprised of the flux integral representing the spatial discretization corresponding to each control volume and the control volume average of the source term. The flux integral is calculated by the Gauss quadrature along the faces of the control volumes. In the present work, we consider the cell-centered approach in which the control volumes are formed as the cells of the mesh. To evaluate the numerical fluxes at the interfaces, primitive variables must be reconstructed to give a series expansion of the solution about the cell’s reference point ðxi ; yi Þ to the desired order of accuracy:

/Ri ðx; yÞ

   @/ @/ @ 2 / ðx  xi Þ2 ¼ /ji þ  ðx  xi Þ þ  ðy  yi Þ þ 2  @x i @y i @x  2 i   @ 2 /  @ 2 / ðy  yi Þ2 þ ... þ  ðx  xi Þðy  yi Þ þ 2  @x@y @y  2 i

ð3Þ

i

kþl

@ /i In Eq. (3), /i is the value of the reconstructed solution and @x k @yl are its derivatives at the reference point of control volume i. A secondorder accurate solution can be achieved by knowing the gradient of the solution at the cell’s reference point and reconstructing a linear polynomial in the control volume. Conservation of the mean within a control volume requires that

1 Ai

Z

Ai

i /Ri dA ¼ /

ð4Þ

By expanding the left-hand side of Eq. (4) term by term, one can easily show that

1 Ai

Z Ai

    @/ @ 2 / x2 i @ 2 /  @ 2 / y2 i xi þ  y i þ 2  þ þ ...  xyi þ 2  @x i @y i @x  2 @x@y @y  2  

@/ /Ri dA ¼ /ji þ 

i

i

i

ð5Þ

222

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

where

xn y m

i

1 ¼ Ai

r/F ¼

Z

n

Ai

m

ðx  xi Þ ðy  yi Þ dA

ð6Þ

 i ¼ / for a second-order scheme if Considering Eq. (5) reveals that / i the cell centroid is chosen as the reference point. Also, the cell gradient can be computed solving a least-squares system which minimizes the error in predicting the mean value of the reconstructed solution for control volumes in the stencil Aj i . The mean value for a single control volume Aj of the reconstructed solution /Ri is

1 Aj

Z Aj

    @/    @/  xj þ xj  xi þ  yj þ yj  yi  @x i @y i     2 ! @ 2 / x2 j þ 2xj xj  xi þ xj  xi þ 2 @x  2 i         @ 2 /   þ  xyj þ xj yj  yi þ xj  xi yj þ xj  xi yj  yi @x@y i     2 ! @ 2 / y2 j þ 2yj yj  yi þ yj  yi þ 2 þ . .. @y  2

/Ri dA ¼ /i þ

i

ð7Þ This equation is written for every control volume within the stencil of control volume i. In this paper, control volumes within the stencil are chosen to be those three triangular cells which are the face neighbors of the target cell. Therefore, the least-squares system for linear reconstruction becomes:

2

3

wi1 ðx1  xi Þ wi1 ðy1  yi Þ 6 7 4 wi2 ðx2  xi Þ wi2 ðy2  yi Þ 5 wi3 ðx3  xi Þ wi3 ðy3  yi Þ

! @/ @x @/ @y

i

  1  / i 1 wi1 /   B 2  / i C ¼ @ wi2 / A   3  / i wi3 / 0

ð8Þ

In Eq. (8) the weights are set to emphasize geometrically nearby data:

1 wij ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2 xi  xj þ yi  yj

ð9Þ

The evaluation of diffusive fluxes requires estimates of the solution gradient at the face. Even though for a structured scheme, a central differencing based on the cell averages of neighboring control volumes is usually used [11], this approach is problematical for unstructured grids, because the cell centroids are often far from lying on the perpendicular bisector of the face between them. This drawback decreases the order of accuracy of computed gradient and therefore affects the total accuracy of solvers. Recasting the Poisson !equation, r2 / ¼ S, in the form of Eq. (1) yields the flux vector as F ¼ r/. One way to construct the face gradient is taking the average of the two adjacent cell gradients obtained by least squares (Eq. (8)). The most obvious candidate is arithmetic averaging in which equal weights are set for the two cell gradients:

r/F ¼

 1 r/i þ r/j 2

   1   ~ r jF r/i þ j~ r iF jr/j   ~ ~ jriF j þ rjF

Nevertheless, it is well understood that the schemes that only use the averaged gradients are susceptible to high frequency errors. As a result, a finite difference term is typically incorporated with the averaged component to enhance the robustness [12]. This idea, which is borrowed from the vertex-centered discretization schemes, first appeared in Crumpton et al. [13] and has been followed by other researchers [14,15] for cell-centered discretization. This correction has also been studied in detail by Puigt et al. [16] for cell vertex codes. In such schemes, the face gradient is constructed such that the projection of the face gradient along the line between the cell centers equals the finite difference term, /j /i j~rij j

(see Fig. 1). Special care should be taken to combine the finite

difference term with the averaged component not to destroy the consistency of discretization. There are two well-known methods for this construction: edge-normal and face-tangent. In the edge-normal   construction, the face gradient in the edge direction, ^ eij ¼ ~ r ij ~ r ij , is replaced with the finite difference term while the edge normal derivative (along ^e? ij ) is specified by the averaged cell gradient [5]:

/ /





r/F ¼ j  i ^eij þ r/  ^e?ij ^e?ij ~ r ij

 1  Ai r/i þ Aj r/j Ai þ Aj

ð13Þ

Alternatively, the finite difference term can be incorporated in the face normal direction for the face-tangent construction:

!  F /j  /i  1    ^t F  ^eij @ / n ^F r/F ¼ ~ ^ F  ^eij n r ij 

ð14Þ

where @ F / represents the face derivative which can be obtained by the projection of the averaged gradient onto the face tangent direction, ^t F :

@ F / ¼ r/  ^t F

ð15Þ

It should be noted that in Eqs. (13 and 15), r/ is the averaged component obtained by one of the three methods mentioned earlier. Although the use of finite difference term alone is common on quadrilateral structured elements, Cary et al. [17] showed that dropping the averaged component results in inconsistent schemes for unstructured triangular grids. Furthermore, Cary et al. [17] introduced a new construction of the face gradient with a finite difference term in a so-called corrected normal direction. In this scheme, the face gradient is directly obtained by a finite difference approximation based on the reconstructed solutions along the face normal direction:

ð10Þ

However, the weight of averaging can be tuned by some geometrical considerations within the cell. Volume (area in 2D) weighted averaging (Eq. (11)) or linear interpolation based on the distance of the opposite cell centroid to the face midpoint (Eq. (12)) can be used as alternatives:

r/F ¼

ð12Þ

ð11Þ Fig. 1. Mesh geometric properties.

223

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

   /cj  /ci ^ F ; lc ¼ n ^ F  ^eij ~ n rij  lc     /ci ¼/i þ r/i  ~ tF ~ tF riF  ~    c ~ ~ ~ /j ¼/j þ r/j  tF rjF  t F

r/F ¼

ð16Þ

ð17Þ

Even though this construction is independent of the averaged component, it is possible to show that this scheme is a variant of the face-tangent construction (Eq. (14)), where the face derivative is obtained by a more complicated averaging strategy, and thus is consistent (Note that ~ r Fj ¼ ~ r jF and ~ rij ¼ ~ riF þ ~ r Fj ):

@F / ¼











r/i  ~t F ~ t F þ r/j  ~ tF ~ tF riF  ~ rFj  ~ ~ rij  ~ tF

 ð18Þ

In addition to the finite difference term, a jump term which comes from the discontinuous solution at the face center can be added to the face gradient:

r/F; jump ¼ ð/F þ  /F  ÞG n^ F

ð19Þ

In the jump term formulation, /F þ is the value at the face center consistent with the representation in the neighbor cell and /F  is the value at the face center obtained from the original cell representation:

/F þ ¼ /j þ ~ r jF  r/j ;

/F  ¼ /i þ ~ riF  r/i

ð20Þ

Therefore, the jump term has an order property and can be simply added to the other components of the face gradient without violating the consistency of a discretization scheme. Also, G is a geometric term needed to be Oð1=hÞ to provide the correct dimension for gradient. Cary et al. [17] suggested this   terms as lF = Ai þ Aj where lF is the length of the face for which the gradient is estimated and Ai and Aj are the areas of the two adjacent cells. Doing some geometric work proves that this   ^ F  . So we will write the geometric term is equal to 2= 3~ r ij  n jump term as Eq. (21) to reproduce Nishikawa’s diffusion scheme [8] derived based on a hyperbolic model for diffusion problems as well. In the context of finite volume, this scheme happens to be in the form of arithmetic averaging and jump term while it is different for other numerical methods. The jump term embedded in this scheme has the role of stabilizing by damping high-frequency errors.

a

 ð/ þ  /F  Þn ^F r/F; jump ¼  ~ ^F  F r ij :n

ð21Þ

where a is an arbitrary value which is called the jump term coefficient. Moreover, Nishikawa [8] showed that the edge-normal and face-tangent constructions of the face gradient can be re-cast in the form of Eq. (21) with special choices of a. In the edge normal construction, the jump term coefficient is

   a ¼ n^ F  ^eij n^F  ^eij 

ð22Þ

and thus the face gradient becomes:



r/F ¼ n^ F  ^eij

!  /F þ  /F    ^F n ~ rij 

  n ^  ^e 

ð24Þ

that yields:

!

r/F ¼

1 /F þ  /F    ^F n ~ ^F  ^eij n r ij 

3. Error eigenanalysis Recent work by Jalali and Ollivier-Gooch [19] revealed that the truncation error for unstructured finite volume discretizations has a noisy appearance due to the discontinuous behavior of coefficients in the Taylor series expansion of the truncation error. This is in contrast to structured meshes where the truncation error is smooth. To address how this rough pattern influences the asymptotic order of accuracy of the truncation error, Ollivier-Gooch and Roy [3] Table 1 Face gradient component. Key

Description

0 A V L

Cell gradient averaging No cell averaging Arithmetic averaging Volume-weighted Linear interpolation

0 J

Face discontinuity correction No jump term Jump term

0 E N C

Finite difference correction No finite difference Cell-center ‘‘Edge’’ direction Face normal direction Corrected normal direction

ð23Þ

Likewise, the face-tangent construction can be written as a jump term with

a ¼ ^ F ^ij nF  eij

The face tangent construction is known to be more robust. This was explained by Thomas et al. [18] by showing that in such a construction the contribution of the finite difference term, which has the damping role, is much larger than with the edge-normal formulation on skewed grids. Expressing the finite difference schemes in the form of jump term makes this more clear as explained by ^ F  ^eij , which is present in the denominator Nishikawa [8] since n of the damping term (Eq. (25)), becomes smaller on skewed grids and thus amplifies the damping. On the other hand, the damping term vanishes on highly stretched triangles in the edge-normal construction (Eq. (23)). This also simplifies the construction of finite difference term as the jump term formulation can be easily added to the averaged component. Meanwhile, the definition of the face tangent vector, which is ambiguous particularly in three dimensions, is removed in the face-tangent construction [8]. In this paper, an extensive family of second-order-accurate discretization schemes that are used for diffusive flux calculation on unstructured meshes are compared in terms of the truncation and discretization errors produced. The schemes differ in the way by which the average of cell gradients is taken, the presence or absence of the face jump term and the type of finite difference correction term. Methods are described by an alphanumeric triplet according to the keys listed in Table 1 that indicates the averaging, whether the jump condition is used and the type of finite difference correction used (if any). For instance, VJE is deciphered as a scheme in which a volume-weighted average is performed of the cell center gradients. Then the finite difference component in the direction of the line connecting the cell-centers is computed. The averaged gradient is projected onto the direction perpendicular to the line connecting the cell centers and added to the finite difference term. Finally, the jump term is added to yield the face gradient used in the flux integration. In addition, Table 2 gives a complete definition of all those schemes tested in this paper. For the sake of simplicity, the finite difference constructions have been described in the form of jump formulation. Also, the schemes that have the finite difference term in the corrected normal direction do not need an averaging scheme and thus the averaging key is replaced with 0.

ð25Þ

224

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

A similar expansion for the truncation error gives:

Table 2 Mathematical relations for different face gradient constructions.

s¼

Mathematical formulation for r/F   1 2 r/i þ r/j   1 Ai þAj Ai r/i þ Aj r/j

Scheme A00 V00 L00

j~rjF j r iF jþj~ rjF j j~

AJ0

 1

LJ0



r/i þ j~r 

A0E AJE 00C 0JC

F

þ

ai ¼ bi ki

riF j j~ r/j j~rjF j

which links the eigendecomposition of the two forms of error. Since the discrete Jacobian is not symmetric, the right xj and left eigenvectors yi are different. We choose to normalize the right  eigenvectors so that xj  ¼ 1. With proper normalization of the left eigenvectors, the two sets of eigenvectors are orthonormal:

yi  xj ¼ dij

proposed an eigenanalysis strategy for unstructured finite volume schemes. In this paper, we will use the same strategy to show how the parameters in different discretization schemes change the accuracy of the numerical solution. Considering that this analysis requires an estimate of both the truncation and discretization errors, we will use the method of manufactured solutions. In this way, we can easily compute the discretization error, e, by comparing the exact manufactured solution, U e , with the numerical solution, U h :

e ¼ Uh  Ue

ð26Þ

Also, the truncation error, s, is obtained by computing the residual based on the exact solution. It should be mentioned that the control volume average of the source term vector in the formulation of the residual, Eq. (1), will be unchanged using either the numerical or exact solutions for the Poisson problem where the source term is not a function of solution. Hence,

I

~  ~ F Ue

!

^ ds n

ð27Þ

CSi

e¼s

ð28Þ

@Rh in which the @U is the discrete Jacobian matrix. Our flow solver is h capable of calculating the global Jacobian explicitly, for schemes up to fourth-order accuracy based on the way the discretization schemes estimates the flux vector and thus different schemes discussed in Section 2 result in distinct flux Jacobians. Comprehensive details of explicit flux Jacobian calculation have been given by Michalak [20]. Both truncation and discretization errors can be expanded as a combination of right eigenvectors of the Jacobian. The corresponding eigenvalues, ki , and right eigenvectors, xi , are, by definition:

@Rh xi ¼ ki xi @U h

X X X ai xi ¼ ai yi  xi ¼ ai dij ¼ aj i

i

ð35Þ

i

and consequently the discretization error coefficients by Eq. (33). The flow solver in this work uses the PETSc numerical linear algebra library [21], so the SLEPc eigenvalue package [22], which uses the same data structures as PETSc, is an obvious choice. The left and right eigenvectors of the Jacobian matrix are computed by the eigensolver built in SLEPc. The eigenvalue spectrum for a second-order cell-centered discretization of the Poisson equation using the A00 scheme, which employs arithmetic averaging without the jump term or finite difference correction, is shown in Fig. 2 for a mesh with 1810 cells. Since the Jacobian matrix is asymmetric, the eigenvalues are complex, though the magnitude of the imaginary part never exceeds about 10% of the magnitude of the real part. Therefore, the magnitude of the real part is a good measure of the eigenvalues magnitude. Also, the real part of all eigenvalues are negative which assures the stability of the discretization considering that our flow solver uses implicit time-stepping with Newton-GMRES iteration for converging to the solution. 4. Results and discussion

The discretization error can be directly related to the truncation error by the error transport equation [10]:

@Rh @U h

ð34Þ

This property can be used to evaluate the coefficients in the eigendecomposition of the truncation error:

^F ð /F þ  /F  Þ n

yi  s ¼ yi 

1 s¼ ACV i

ð33Þ

iF jþ

ij

j~rij :^nF j

ð32Þ

and substitution into the error equation leads to

r/i þ r/j þ

a

ai xi

i

a ð/ þ  /  Þn F ^F j~rij :n^ F j F

j~rjF j r iF j j~ ^F r/i þ j~r jþ ~r r/j þ ~r a:^n ð/F þ  /F  Þn r iF jþj~ r jF j j~ j jF j j ij F j iF

    /F þ /F  1 ^ ^ ^F n 2 r/i þ r/j þ nF  eij j~rij j

    / þ /F 1 a ð/ þ  /  Þn F ^F  ^ ^F eij n F ^F þ n F 2 r/i þ r/j þ j~ ^F j rij :n j~rij j        1 ~ ~ ~ ~ ~ ~ ^ r  / r n / þ r /  t  t  r /  t  t F F F F F jF iF j j i i ^ F :~ n r ij        1 ^F r jF  ~ r iF  ~ /j þ r /j  ~ tF ~ t F  /i  r/i  ~ tF ~ tF n ^ :~ r n 2

X

In this section, a wide range of discretization schemes used for diffusive fluxes are studied in terms of accuracy on isotropic and anisotropic unstructured meshes. The schemes are named based on the triplet discussed in Section 2 and used to discretize the Poisson equation as a scalar model of diffusion problems. The truncation and discretization errors are obtained by the method

ð29Þ

By expanding the discretization error as a linear combination of these eigenvectors

e¼

X

bi xi

ð30Þ

i

and substituting into the left-hand side of the error equation, we get:

X @Rh X @Rh @R X e¼ h bi xi ¼ bi xi ¼ bi ki xi @U h @U h i @U h i i

ð31Þ Fig. 2. Eigenvalue spectra using A00 on an unstructured isotropic mesh.

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

of manufactured solutions and compared on a sequence of triangular meshes to study the asymptotic order of accuracy for these two forms of error. The eigenanalysis tool introduced in Section 3 is also utilized to explain the error behaviors of different schemes. This tool is also helpful in identifying the schemes that are not stable for anisotropic meshes as will be shown later. 4.1. Isotropic meshes To compare the truncation and discretization errors produced by various discretization schemes, the L2 -norm of errors are evaluated. The Poisson equation is solved inside a unit square exposed to homogeneous boundary conditions for which the exact solution is assumed to be sin ðpxÞ sin ðpyÞ. Having the exact solution, it is possible to obtain the corresponding source term. The comparisons are performed for a sequence of meshes where isotropic triangles tessellate the unit square and their length scales are reduced by a factor of 2 at every level of refinement. All the cells have geometrically high quality obtained by Delauny refinement in which small angles are avoided and their sizes remain nearly constant across the domain which is ideal for an isotropic solution. Fig. 3 shows the exact solution on one of these meshes with 1810 triangles. The discretization schemes used for diffusive fluxes can be categorized into four different families: pure averaging; averaging and the jump term; averaging and finite difference correction; and averaging plus the jump term and finite difference correction.

The error results are presented in this section for each of these families in the described order to investigate how different variables change the solution accuracy. In the pure averaging family, the schemes are distinguished by the method of averaging. Based on the description of Table 1, the reconstructed gradients at the two sides of a face can be averaged by arithmetic method which assigns equal weights, A00; linear interpolation, L00; or volume-weighted method, V00. Fig. 4 shows the L2 -norm of truncation and discretization errors for these three schemes versus the average length scales of the triangles in the mesh. Fig. 4b shows that all three schemes produce second-order accurate solutions. While arithmetic averaging and linear interpolation are comparable in terms of discretization error norm, the error that volume weighted averaging produces is larger by about a factor of two. This behavior can be explained by considering the fact that volume-weighted averaging emphasizes the effect of larger control volume although the centroid of the larger triangle is farther than the centroid of the smaller triangle from the face midpoint. This defect has been modified in linear interpolation where the weight of averaging is determined by the distance of the opposite cell centroid to the face midpoint. Likewise, Fig. 4a reveals that the truncation error produced by volume-weighted averaging is the largest, though it is Oð1Þ for all schemes. This is in contrast to the structured meshes where both the discretization and truncation errors have the same asymptotic order of accuracy. For unstructured meshes, the discretization order typically matches the order of solution reconstruction, whereas the truncation error is asymptotically worse [5,6]. To explain this behavior, it is useful to consider the flux integral for the Poisson equation:

FIi ¼

Fig. 3. Manufactured solution for Poisson’s equation on isotropic meshes.

225

1 Ai

I

r/  n^ ds

ð36Þ

@Ai

For both structured and unstructured meshes, the length

of the con2 trol volume interface and its area are OðhÞ and O h , respectively. However, the gradient on a structured mesh is evaluated by a second-order central differencing. Furthermore, the error in the gradient varies smoothly which produces additional cancellation during integration and provides a second order flux integral. For an unstructured mesh, the gradient is obtained by the reconstructed data which is only first-order accurate and the gradient error is not smooth. Consequently, no cancellation occurs during integration and the flux integral is zero-order accurate. The second family of schemes include averaging and the jump term. The schemes using volume-weighted averaging are

Fig. 4. L2 -Norm of errors for pure averaging discretizations.

226

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

Fig. 5. L2 -Norm of discretization errors for schemes with averaging and jump term.

eliminated from further investigation due to their poor performance. Therefore, this family is limited to AJ0 and LJ0. A challenge in applying the jump term is the choice of the jump term coefficient, a, that influences the flux integral. To address the effect of a on the accuracy of these schemes, the L2 -norm of discretization error resulting from the AJ0 and LJ0 schemes have been plotted in terms of a for the meshes whose average length scales are shown in Fig. 5. As seen in this figure, the use of jump term does not change the asymptotic order of accuracy as the error norms are reduced consistently by a factor of about four after each level of refinement. Also, any value set for the jump term coefficient in the range of ð0; 2, improves the accuracy of the numerical solution. The optimal value of the jump term coefficient by which the smallest error is obtained is not the same for the two schemes. In AJ0; a ¼ 43 gives the most accurate solution which is in agreement with Nishikawa’s suggested value [7] for this scheme that was obtained by considering the coefficient that gives a fourth-order accurate flux integral in one-dimension and multi-dimensional Cartesian grids. On the other hand, the numerical experiment suggests a ¼ 1 for the LJ0 scheme. Even though these optimal values are presented based on the particular solution of Poisson’s equation assumed in this section, they have proven to be consistent for all isotropic solutions on unstructured meshes. The next class of schemes are those which utilize arithmetic averaging or linear interpolation of cell gradients and have a finite difference correction term in the edge, normal or corrected normal

direction. In this study, we only consider the finite difference term in the edge and corrected normal directions. Note that the finite difference construction in the corrected normal direction is not dependent on the type of averaging and there will be no significance in choosing an averaging scheme. The schemes remain under this category are A0E, L0E and 00C. Fig. 6 compares the L2 -norm of truncation and discretization errors for this family of schemes with methods from the previous tests that performed well. As seen in Fig. 6b, all these schemes are perfectly second-order accurate except 00C which shows some variation on coarse to medium meshes. For the schemes based on the edge direction finite differencing, the choice of averaging method is less significant such that the error plots for A0E and L0E can not be easily distinguished. While the use of finite difference correction increases the numerical accuracy of the solution, it is not as effective as the schemes with the optimal jump term and without finite difference correction. Furthermore, the asymptotic order of truncation error for all these schemes is zero as expected but in conflict with the discretization error behavior as the L0E produces the largest truncation error. This is also in contrast with the pure averaging methods where the schemes had the same trends in both forms of error. This behavior can be explained by the eigenanalysis paradigm described earlier. We performed a full eigendecomposition of both types of error on the 1810 triangle mesh. Fig. 7 shows the magnitudes of the weights, ai and bi , plotted against the real part

Fig. 6. L2 -Norm of errors for schemes with jump term or finite difference correction.

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

227

Fig. 7. Eigendecomposition of truncation (left) and discretization (right) errors.

of the eigenvalue for four different schemes. To visualize the error weights more effectively, the logarithmic scale for both weights and real part of the eigenvalue is used; however, ReðkÞ is negative

and hence the negative of the real part of the eigenvalue, ReðkÞ , is plotted instead. A careful investigation of this figure reveals that the modes corresponding to the large eigenvalues are dominant

228

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

Table 3 Leftmost eigenvalues for different mesh sizes and L00 discretization scheme. Number of Control Volumes (NCV)

Min EV (leftmost EV)

Min kEV k NCV

462 1810 7176 28,670

2497.36 12801.21 53466.68 209864.41

5.40 7.07 7.45 7.32

in the eigendecomposition of truncation error and are associated with lager weights. Table 3 shows the largest eigenvalue magnitude for different mesh sizes and L00 discretization scheme. The last column of the table illustrates that the leftmost eigenvalue scales nearly linearly with the number of control volumes, and 2 therefore is proportional to h . This implies that the high frequency modes in the truncation which are dominant, will error,

2 be de-amplified by a factor of O h as we solve the linear system, restoring the design order of accuracy for the discretization error. It is worth noting that discretization schemes alter the flux Jacobian matrix and thus its eigenvalues. Moreover, the choice of discretization method changes the truncation error at the first place and produces a distinct set of weights, ai . Discretization error, which is the consequence of inaccurate PDE approximation (truncation error), is the secondary effect as no error appears in the solution if the truncation error becomes zero (Eq. (28)). The weights in the eigendecomposition of the discretization error can be automatically obtained by knowing the truncation error weights and the flux Jacobian eigenvalues, bi ¼ ai =ki . The weights computed in this way were verified to identically match the weights obtained directly from the discretization error eigendecomposition. This brief discussion suggests that the eigendecomposition of the truncation error can be used to explain the accuracy results. These plots can be interpreted as the fingerprint of each discretization family as they remain qualitatively the same for different meshes. One can compare the eigendecomposition of truncation error for the two schemes of the pure averaging family, L00 and V00 in Fig. 7 (left-side). Since the nature of discretization in pure averaging schemes is the same, minor changes occur in the flux Jacobian matrix and its eigenvalues and eigenvectors. So the eigendecomposition plots have similar pattern as the weights corresponding to medium eigenvalues dominate the truncation error although they are slightly larger for the V00 scheme. Also, the range of eigenvalues for the two schemes are comparable and consequently the weights in the eigendecomposition of discretization error become larger for V00. This is in agreement with the accuracy results shown in Fig. 4 where the schemes with higher truncation error are less accurate. Likewise, the comparison between the eigendecomposition of truncation error for L00 and LJ0 reveals that adding a jump term pushes the eigenvalues towards larger magnitudes and enhances the stability of the solution process. Also, the jump term changes the eigenvectors in a way that the weights of truncation error decomposition become smaller by one order and this is the reason for which LJ0 exhibits the smallest truncation error in Fig. 6a. The significant distinction caused by the jump term in the pattern of eigendecomposition is that the weights corresponding to the larger eigenvalues dominate the truncation error. All these effects together yields a stronger de-amplification of truncation error weights such that the weights of discretization error are smaller by one order of magnitude compared to the pure averaging schemes. These smaller weights result in a smaller discretization error as shown in Fig. 6b. Finally, a study of eigendecomposition plot for L0E clarifies the error properties obtained by this scheme. As seen in Fig. 7d, the

weights corresponding to the truncation error are comparable with pure averaging schemes, though a larger fraction of them have large values and thus the truncation error associated with L0E is larger than L00. Meanwhile, the finite difference term increases the magnitude of large eigenvalues and makes ai s corresponding to them dominant. In this way, the discretization error coefficients become smaller and yield a smaller discretization error. Also, a comparison between the eigendecomposition plots LJ0 and L0E reveals that the ranges of eigenvalues for the two schemes are fairly similar. However, the discretization error produced by LJ0 is smaller because the weights in the eigendecomposition of truncation error are smaller. Now, we turn our attention back to the discretization scheme families and present the error results for those where both the jump term and finite difference correction are applied. For the sake of brevity, we just consider LJE and 0JC. It is worth mentioning that adding the finite difference correction term does not preserve the optimal jump term coefficients previously obtained for schemes without it. It is possible to tune the jump term coefficient in a way to produce more accurate results compared to those with single jump term or finite difference correction. The optimal jump term corresponding to each scheme is reported in Fig. 8. As expected, the implementation of the two correction terms makes the schemes more accurate provided that the right jump term coefficient is applied. This can be explained by the stabilizing property of the jump term which moves the eigenvalues to the left and increasing their magnitudes and also the ability of finite difference correction in dominating larger eigenvalues. In other words, these schemes combine the positive features of finite difference correction and jump term to improve the overall accuracy of numerical solution. 4.2. Anisotropic meshes To assess the accuracy of different discretization schemes on anisotropic meshes, an irregular anisotropic grid is used as a test case. This kind of mesh is generated on a rectangular domain ðx; yÞ 2 ½0; 1  ½0; 0:05 using the procedure that has been described in detail by Diskin et al. [5]. Fig. 9 depicts the anisotropic mesh test case with 20  20 cells. For generating the test case, the first step is stretching a regular rectangular grid with ðN þ 1Þ  ðN þ 1Þ nodes toward the bottom line y ¼ 0 using a stretching factor s and the maximum aspect ratio ARmax . The y-coordinates of the horizontal grid lines in the top half of the domain are defined as:

yk ¼ yk1 þ

sk ARmax  N

ð37Þ

in which

y0 ¼ 0;

k ¼ 1; 2; . . . ; N

ð38Þ

Then, irregularities are introduced by random shifts on interior nodes in vertical and horizontal directions and finally each perturbed quadrilateral is randomly triangulated with one of the two choices for its diagonal. In Fig. 9 where N ¼ 20, the stretching factor and the maximum aspect ratio are set as s ¼ 1:14 and ARmax ¼ 100 to yield triangles whose aspect ratios vary approximately between 10 and 100. This mesh consists of anisotropic triangles which are long in the x-direction and thin in the y-direction. In this study, these meshes with N ¼ 10; 20; 40; 80 are used where the stretching factors are s ¼ 1:28, 1.14, 1.065, 1.0335, respectively. The s values are chosen so as to keep the length scale reduction in the y-direction the same as the x-direction. Since the anisotropic meshes are used for anisotropic solutions, the function being used for test purposes is manufactured to comply with the

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

229

Fig. 8. L2 -Norm of errors for schemes with jump term and/or finite difference correction.

Fig. 9. Anisotropic solution and derivative on a stretched grid with 20  20 cells.

cells’ aspect ratio. The anisotropic mesh of Fig. 9 is stretched toward the horizontal line y ¼ 0 using Eq. (37) which is strongly dependent on the y-coordinates. The variations in the x-direction are introduced by an exponential function to yield isotropic behavior in this direction and homogeneous boundary conditions:

/ðx; yÞ ¼

exp ðxðx  1Þyðy  0:05ÞÞ  1   s þ N 1  1s y ARmax

ð39Þ

In this part, Eq. (39) is used as the test function with geometric values corresponding to the coarsest mesh (N ¼ 10 and s ¼ 1:28). Fig. 9 depicts the manufactured solution on one of the constructed

anisotropic mesh along with its y-derivative which is too large near the bottom line with highly stretched cells and is reduced smoothly in the vertical direction compliant with cells’ aspect ratio. For meshes with highly anisotropic cells, the pure averaging family of discretization schemes fails due to stability issues. For these scheme (such as A00 and L00), there are some eigenvalues in the right-half plane that make the solution procedure unstable. The eigenanalysis of the flux Jacobian matrix shows that adding a jump term is effective in pushing all the eigenvalues to the left-half plane. This is in agreement with the stabilizing property of the damping term appears as a jump formulation in finite volume discretization. However, the damping that occurs by adding a finite difference correction in the edge direction is not sufficient to make the scheme stable. This has already been verified by Nishikawa [8], who showed that the damping term in Avg-LSQ scheme (corresponding to A0E scheme here) will be vanishing on highlyskewed triangular grids. Fig. 10 compares a portion of eigenvalue spectra between L00 and LJ0 as well as L00 and L0E. As is clear in these figures, some eigenvalues are positive using L00 and the damping produced by L0E is not sufficient to push the eigenvalues to the left-half plane and thus they are both unstable. Therefore, the number of schemes can be used for anisotropic meshes is less compared to the isotropic ones. It should be noted that applying the finite difference correction in the corrected normal direction adds sufficient damping to the scheme and stabilizes 00C and 0JC. Fig. 11 shows the discretization error produced by several discretization schemes on a sequence of anisotropic meshes where h ¼ 1=N. The error plots prove that all these scheme lead to nearly second-order solutions except 00C that does not converge perfectly linearly. The results obtained by AJ0 and LJ0 are very similar as their difference can not be distinguished on this plot. Our numerical tests showed that a ¼ 43 exhibits the optimal result for AJ0 previously verified by Nishikawa [7] for anisotropic meshes and since LJ0 behaves the same as AJ0 for anisotropic cells, the optimal jump term coefficient does not change. This is in contrast to the isotropic cases where we obtained different optimal values for the two schemes. The use of 00C causes a less accurate solution and shows more sensitivity to mesh perturbations and refinements; however, adding the jump term to this scheme improves the accuracy of the solution as well as convergence behavior. 0JC converges almost linearly and is more accurate than any other scheme on the coarse meshes although its performance degrades on the finer grids. For the schemes with the jump and finite difference correction, we just consider AJE that can be seen as adding a finite difference correction to AJ0 with the same coefficient. Nevertheless, the

230

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

Fig. 10. Comparison of eigenvalue spectra for some discretizations on anisotropic meshes.

Fig. 11. Discretization error on anisotropic meshes.

difference between the accuracy of AJE and AJ0 is quite small. This can be explained by considering the idealized anisotropic triangles depicted in Fig. 12. This ideal mesh consists of regular right triangles whose aspect ratios are AR. Recognizing that the contribution of flux vector determined at edge ab (a very small edge) is negligible in the flux integral at high aspect ratio, two possible geometrical configurations for edges bc and ca have been illustrated in Fig. 12a and b. In Fig. 12a, it is quite straightforward to show that the cosines of h and c, the angles between the finite difference correction vectors and corresponding edges, are

Fig. 12. Anisotropic mesh configurations.

cos h ¼

AR2  1 AR2 þ 1

AR cos c ¼

12 2 AR þ 4

ð40Þ

For large values of AR, these relations can be approximated as

cos hu cos cu1 



1 þ O AR2 AR4 2

ð41Þ

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

Hence, these two angles are very close to zero for anisotropic mesh with high aspect ratio. This implies that the finite difference vectors are nearly tangential to edges bc and ca and thus the finite difference correction does not play a significant role in flux integration. Likewise, in Fig. 12b, cos h is evaluated through Eq. (40) and the finite difference term has an insignificant effect on the face gradient of edge ca provided that high aspect ratio triangles are used. For this configuration, c ¼ p2 , so the finite difference term is the only term in the gradient of face bc. Under this condition, the averaged component coming from the least-squares reconstruction is projected on the perpendicular direction to ~ r bc ; this component, which is along edge bc, does not influence the flux integral. Therefore, we expect a different flux integral by adding a finite difference term in the edge direction. The key point is that in computing the cell gradient using the least-squares reconstruction, nearby control volume averages are emphasized by introducing the geometrical weights in Eq. (8). For the mesh configuration of Fig. 12b, j~ r ca j, r bc j  j~ r ab j; j~ so the y-component of the computed cell gradient for CV i is mainly dependent on CV 1 and vice versa. This suggests that the average of the cell gradients for face bc is approximately equal to the finite difference approximation for that face. Therefore, replacing the face normal part of the cell gradient with the finite difference term in the edge direction does not change the truncation error for these canonical configurations of anisotropic mesh. Although this argument was described for an idealized mesh, it is easily generalizable to irregular anisotropic grids. 5. Conclusions Different discretization schemes that are widely used for the computation of diffusive fluxes on unstructured meshes were compared in terms of truncation and discretization errors. The schemes differ in the way by which the cell gradients are combined to approximate the face gradient, the presence or absence of a face jump term and the type of finite difference correction term. The method of manufactured solutions was employed to compute the two forms of error obtained by each discretization scheme in a second-order cell-center finite volume solver and Poisson’s equation was used as a model of viscous discretization. The comparisons were carried out for both isotropic and anisotropic unstructured meshes for which the solution anisotropy was consistent with the cells’ aspect ratio. In addition, an eigendecomposition paradigm was introduced by which both the truncation and discretization errors can be written as a linear combination of flux Jacobian matrix eigenvectors. The weights corresponding to each eigenvector in the linear expansions of the two forms of error are related by the eigenvalues of the flux Jacobian matrix. The choice of discretization schemes alters the Jacobian matrix, its eigenvalues and eigenvectors and thus the weights corresponding to each eigenvector. The eigendecomposition tool is able to interpret how each family of schemes influences accuracy behavior by showing the distribution of weights and eigenvalues. Moreover, the eigendecomposition is helpful in identifying those scheme that have eigenvalues in the right-half plane and consequently are not stable for anisotropic meshes. Note that the eigendecomposition is a posteriori process and requires the Jacobian matrix to provide the eigenvalue spectrum. Even though one can use eigendecomposition to obtain the Jacobian matrix for a desired eigenvalue spectrum, the construction of the diffusion scheme that yields this Jacobian is not trivial. Our accuracy analysis suggests that volume-weighted averaging for combining the two adjacent cells’ reconstructed gradients results in larger errors compared to arithmetic averaging and linear interpolation. This is in agreement with our eigenanalysis results where the coefficients in this case are larger. Thus, we eliminated from further considerations those discretization schemes that

231

utilize volume-weighted average of the cell gradients. Furthermore, adding the jump term along with the optimal jump term coefficient reduces both the truncation and discretization errors considerably. Also, this term has a significant effect in stabilizing the scheme and pushing the eigenvalues to the left-half plane. The optimal jump term coefficient is dependent on the choice of averaging scheme but is independent of mesh size for isotropic triangles. Adding a finite difference term in the edge direction results in larger truncation error compared to pure averaging schemes with correspondingly larger expansion coefficients. However, this term has a damping effect and produces eigenvalues that are considerably larger in magnitude so as it produces smaller discretization errors. The schemes with finite difference correction are more accurate than the pure averaging family although they are not as accurate as those with jump term. So those schemes that combine the advantages of both terms present the best accuracy results for isotropic meshes. Based on our results, LJE with the optimal coefficient of 23 is recommended for isotropic meshes. For anisotropic meshes, the fewer schemes can be used for viscous discretizations due to stability issues. The pure averaging schemes fail to converge since some of the eigenvalues fall into the right-half plane. The damping effect of finite difference correction term in the edge direction is not sufficient to push all eigenvalues to the left. Therefore, only the schemes with the jump term and/or finite difference in the corrected normal direction are stable. Even though linear interpolation performs somewhat better than arithmetic averaging for schemes with the jump term on isotropic grids, the difference is less clear for high aspect ratio meshes. Adding a finite difference correction in the edge direction to this scheme does not help. Also, the scheme with finite difference correction without the jump term does not converge linearly and produces the least accurate results. Therefore, AJ0 with the optimal coefficient of 43 is recommended as is computationally cheaper and accurate enough for anisotropic meshes. This paper presents a new technique to analyze the accuracy of discretization schemes used for diffusive fluxes in cell-centered unstructured finite volume solvers. Future work includes the same type of analysis for various scheme choices used in vertex-centered solvers. In addition, the effect of cell gradient accuracy on the overall accuracy of each scheme should be investigated. This requires the comparison of different methods for cell gradient calculation such as unweighted/weighted least-squares and different Green-Gauss formulations.

References [1] Mavriplis DJ. Grid resolution study of a drag prediction workshop configuration using the NSU3D unstructured mesh solver. In: Proceedings of the twenty-third AIAA applied aerodynamics conference; June 2005 [AIAA paper 2005-4729]. [2] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. In: Proceedings of the thirtieth AIAA aerospace sciences meeting and exhibit; January 1992 [AIAA paper 92-0439]. [3] Ollivier-Gooch C, Roy C. Reducing truncation error on unstructured meshes by vertex movement. In: Proceedings of the twentieth AIAA computational fluid dynamics conference; 2012 [AIAA paper 2012-2712]. [4] Roe PL. Approximate Riemann solvers, parameter vectors, and difference schemes. J Comput Phys 1981;43:357–72. [5] Diskin B, Thomas JL, Nielsen EJ, Nishikawa H, White JA. Comparison of nodecentered and cell-centered unstructured finite-volume discretizations. Part I: Viscous fluxes. In: Proceedings of the forty-seventh AIAA aerospace sciences meeting; 2009 [AIAA paper 2009-597]. [6] Diskin B, Thomas JL. Comparison of node-centered and cell-centered unstructured finite-volume discretizations: inviscid fluxes. In: Proceedings of the forty-eighth AIAA aerospace sciences meeting; 2010 [AIAA paper 20101079]. [7] Nishikawa H. Robust and accurate viscous discretization via upwind scheme – I: Basic principle. Comput Fluids 2011;49:62–86. [8] Nishikawa H. Beyond interface gradient: a general principle for constructing diffusion scheme. In: Proceedings of the fortieth AIAA fluid dynamics conference and exhibit; 2010 [AIAA paper 2010-5093].

232

A. Jalali et al. / Computers & Fluids 101 (2014) 220–232

[9] Nishikawa H. Two ways to extend diffusion schemes to Navier–Stokes schemes: gradient formula or upwind flux. In: Proceedings of the twentieth AIAA computational fluid dynamics conference; 2011 [AIAA paper 2011-3043]. [10] Roy C. Review of discretization error estimators in scientific computing. In: Proceedings of the forty-eighth AIAA aerospace sciences meeting; 2010 [AIAA paper 2010-126]. [11] Hirsch C. Numerical computation of internal and external flows: computational methods for inviscid and viscous flows. Wiley; 1990. [12] Haselbacher A, McGuirk JJ, Page GJ. Finite volume discretization aspects for viscous flows on mixed unstructured grids. Am Inst Aeronaut Astronaut J 1999;37:177–84. [13] Crumpton PI, Moinier P, Giles M. An unstructured algorithm for high Reynolds number flows on highly stretched grids. In: Taylor C, Cross JT, editors. Numerical methods in laminar and turbulent flows. Pineridge Press; 1997. p. 561–72. [14] Haselbacher A, Blazek J. Accurate and efficient discretization of Navier Stokes equations on mixed grids. Am Inst Aeronaut Astronaut J 2000;38:2094–102. [15] Eliasson P. EDGE: a Navier–Stokes solver for unstructured grids. Technical report, FOI-R-0298-SE; 2001.

[16] Puigt G, Auffray V, Muller JD. Discretization of diffusive fluxes on hybrid grids. J Comput Phys 2010;229:1425–47. [17] Cary AW, Dorgan AJ, Mani M. Towards accurate flow predictions using unstructured meshes. In: Proceedings of the nineteenth AIAA computational fluid dynamics conference; 2009 [AIAA paper 2009-3650]. [18] Thomas JL, Diskin B, Nishikawa H. A critical study of agglomerated multigrid methods for diffusion on highly-stretched grids. Comput Fluids 2011;41: 82–93. [19] Jalali A, Ollivier-Gooch C. Accuracy assessment of finite volume discretizations of diffusive fluxes on unstructured meshes. In: Proceedings of the fiftieth AIAA aerospace sciences meeting; 2012 [AIAA paper 2012-0608]. [20] Michalak C, Ollivier-Gooch C. Globalized matrix-explicit Newton-GMRES for the high-order accurate solution of the Euler equations. Comput Fluids 2010;39:1156–67. [21] Smith B, Balay S, Brown J, Dalcin L, Karpeev D, Knepley M, et al. PETSc home page, 2011. . [22] Campos C, Roman JE, Romero E, Tomas A. SLEPc home page, 2011. .