A Case-Study in Open-Source CFD Code Verification, Part I: Convergence Rate Loss Diagnosis

Available online at www.sciencedirect.
com
ScienceDirect
Mathematics and Computers in Simulation ( ) –

www.elsevier.com/locate/matcom
Original articles
A case-study in open-source CFD code verification,

Part I: Convergence rate loss diagnosis
H. Noriega a , ∗, F. Guibault a , M. Reggio a , R. Magnan b
a 2900, boul. Édouard-Montpetit, Campus de l’Université de Montréal, 2500, chemin de Polytechnique, Montréal (Québec) H3T 1J4, Canada
b Institut de recherche d’Hydro-Québec, 1800 Boulevard Lionel-Boulet, Varennes, QC J3X 1S1, Québec, Canada
Received 18 January 2017; received in revised form 11 October 2017; accepted 4 December 2017
Available online xxxxx
Abstract
This study analyzes the influence of cell geometry on the numerical accuracy of convection–diffusion operators in OpenFOAM.
The large variety of solvers and boundary conditions in this tool, as well as the precision of the finite-volume method in terms
of mesh quality, call for a verification process performed in steps. The work is divided into two parts. In the first (the current
manuscript), we focus on the diffusion operator, which has been found to exhibit a loss in convergence rate. Although the
cell-centered finite volume approach underlying OpenFOAM should preserve a theoretical second order convergence rate, loss
of convergence order is observed when non-orthogonal meshes are used at the boundaries. To investigate the origins of this
problem, the method of manufactured solutions is applied to yield analytical solutions for the Poisson equation and compute
the numerical error. The root cause is identified and corrections to recover second-order convergence are proposed. In part two
of this investigation, we show how convergence can be improved, and present results for problems described by the Poisson and
Navier–Stokes equations.
⃝c 2017 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights
reserved.
Keywords: Verification; Manufactured solution; OpenFOAM; CFD; Poisson equation
1. Introduction
Numerical calculations in fluid mechanics have developed considerably in recent decades, in part due to advances
in computational power and in part because of advances in numerical methods. However, despite this progress, the
simulation of real-world problems with complex geometries remains a major challenge. In fact, analytical solutions
do not exist for these problems, and instead the equations representing a particular physical problem are solved via
numerical simulation, which involves the use of approximations from various sources.
∗ Corresponding author.
E-mail addresses: [email protected] (H. Noriega), [email protected] (F. Guibault), [email protected] (M.
Reggio), [email protected] (R. Magnan).
https://doi.org/10.1016/j.matcom.2017.12.002
0378-4754/⃝ c 2017 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights
reserved.
Please cite this article in press as: H. Noriega, et al., A case-study in open-source CFD code verification, Part I: Convergence rate loss diagnosis, Mathematics and Computers in
Simulation (2017), https://doi.org/10.1016/j.matcom.2017.12.002.
2 H. Noriega et al. / Mathematics and Computers in Simulation ( ) –
The finite volume method (FVM), based on the laws of conservation, is currently one of the methods most often
used to solve the Navier–Stokes equations for predicting industrial flows. Theoretical support for FVM on complex
geometries can be found in the Computational Fluid Dynamics (CFD) literature [25,2], and [9].
OpenFOAM⃝ R
, a trademark of The OpenFOAM Foundation, is a free, open source software package that is widely
used for fluid flow simulations in both industry and academia. OpenFOAM’s theoretical support, based on cell-
centered finite volume discretization, includes a large variety of solvers, schemes, and processing tools for solving
a wide range of problems in computational fluid mechanics. OpenFOAM support and theoretical background can be
found in [11–13,15,16,14,17,18,20].
Flow simulations are routinely performed in industry, and OpenFOAM accuracy must be verified on the meshes
that usually accompany the complex geometries in these simulations. Because of the complexity and extent of
OpenFOAM libraries and tools, the verification process has to be performed in steps. However, there is little in
the literature on the verification of OpenFOAM schemes, which is why we are addressing this important issue here.
A background discussion, along with definitions and descriptions for some terms related to confidence building in
CFD, was presented by Roache [23] and extended by Stern et al. [26]. There, validation is described as solving the
right equations, and verification as solving the equations in the right way. The various aspects discussed in Roache’s
paper include the distinction between code verification and code validation, between the verification of code and the
verification of calculations, grid convergence vs. iterative convergence, and numerical error vs. conceptual modeling
error.
Along these lines, Abanto et al. [1] presented a grid convergence study on some atypical CFD cases using a
number of commercial CFD packages. Their verification test cases are distinctive in that exact solutions are known.
Verification test cases include the classical Poiseuille flow, an incompressible recirculating flow, a manufactured
incompressible laminar boundary layer flow, and an incompressible annular flow. Different convergence rates are
determined using structured and unstructured meshes. Stern et al. [26] describe a set of verification, validation, and
certification methodologies and procedures for numerical simulations. Examples of the application of quantitative
certification of Reynolds-Averaged Navier–Stokes (RANS) codes are presented for ship hydrodynamics.
In Tremblay et al. [27], the method of manufactured solutions (MMS) for fluid-structure interaction code
verification is applied. These researchers observed that, when used with systematic grid refinement, MMS provides
strong code verification. Roy et al. [24] used MMS to verify the order of accuracy of two finite-volume Euler and
Navier–Stokes codes. These exact solutions were used to accurately evaluate the discretization error in the numerical
solutions. Through global discretization error analysis, the spatial order of accuracy was observed to be second order
for a node-centered approach using unstructured meshes. More recently, Iannelli [10] compared an exact solution
and CFD solutions of the Navier–Stokes equations. They determined the convergence rates and orders of accuracy of
these solutions and illustrated the utility of the exact solution developed for verification purposes. There is a series of
published papers on MMS for 2D incompressible Navier–Stokes and turbulent models, which includes those by Eca
and Hoekstra [6,7] and Eca et al. [8].
A recent study regarding verification for an unstructured finite volume CFD code has been published by Veluri
and Roy [29], in which MMS is used to generate exact solutions to both the Euler and Navier–Stokes equations to
verify the order of accuracy of the code on different grid types in 2D and 3D. The various options for code verification
include the baseline steady-state governing equations, transport models, turbulence models, boundary conditions, and
unsteady flows. Diskin et al. [5] studied the accuracy and complexity in finite volume discretization schemes for
viscous fluxes on general grids using a node-centered scheme and three cell-centered schemes (a node-averaging
scheme and two schemes using least-squares face-gradient reconstruction). Among several interesting results, they
found that the node-averaging scheme has the highest complexity and can fail to converge to the exact solution when
the node-averaged values are clipped. On highly anisotropic grids, the least-squares schemes, the node-averaging
scheme without clipping, and the node-centered scheme demonstrated similar second-order accuracies per degree
of freedom. Overall, the accuracies of the node-centered and the best cell-centered schemes are comparable at an
equivalent number of degrees of freedom on isotropic and curved anisotropic grids. A similar work for inviscid
fluxes is presented by Diskin and Thomas [3], the second-order cell-centered and node-centered approaches for finite
volumes were compared for unstructured grid discretizations in two dimensions. Some weaknesses were observed in
all these schemes, such as instability, accuracy degradation, and/or poor convergence of defect-correction iterations.
In a more recent work, Diskin and Thomas [4] studied the effects of mesh regularity on the accuracy of unstructured
node-centered finite-volume discretizations. In their paper, they focused on an edge-based approach which uses
unweighted least-squares gradient reconstruction with a quadratic fit.
H. Noriega et al. / Mathematics and Computers in Simulation ( ) – 3
Along these lines, the goal of the present two-part study is the verification of OpenFOAM libraries used by solvers,
and detection of the influence of the cell geometry on the numerical accuracy of OpenFOAM Convection–Diffusion
operators. Here, in the first part of our study, we focus on the diffusion operator. To ensure the correctness of
the code libraries, the theoretical second order convergence rate for the cell-centered approach must be preserved.
We have considered the verification of the OpenFOAM convergence rate under three types of boundary conditions
(Dirichlet, Neumann, and Periodic). Three criteria of mesh quality (non-orthogonality, non-uniformity, and skewness)
are taken into consideration in the verification process. For this purpose, we investigate how the error induced by mesh
distortions affects the convergence rate of numerical simulations under systematic mesh refinement.
The paper is organized as follows: first, we present the formulation of the mathematical models, followed by the
meshes used, the manufactured solutions applied, and the errors introduced. The following section presents numerical
results: first on orthogonal meshes, and then on non-orthogonal ones. For the latter, a reduction in convergence rate
compared to what is predicted by theory is observed. This loss of accuracy is further explored and related to the
treatment of non-orthogonality at the boundary. The root cause of this loss of accuracy is identified and a solution is
proposed. The final section presents our concluding remarks and perspectives.
2. Mathematical models and discretization

Here, in this first part of our study, we reduce the verification procedure for the convection and diffusion schemes to
a convergence analysis on the diffusion scheme, for which a reduction in convergence rate has already been observed
on non-orthogonal meshes, as we will see later. We will begin by considering the Poisson equation, written as:
− ∇ · (ν∇φ) = q (1)
in which φ represents a scalar field, ν the diffusion coefficient, and q the source term.
To discretize Eq. (1) using the FVM, we integrate over the elementary control volume V p limited by the control
surface S p , as follows:
∫ ∮ ∫
∇ · (ν∇φ) dτ = d S · (ν∇φ) =
⃗ qφ (⃗
x )dτ. (2)
Vp Sp Vp
A brief introduction to the FVM discretization of Eq. (2) is provided below.
2.1. Numerical schemes
An overview of the theoretical basis of the OpenFOAM schemes analyzed in this work is presented below
(see [11,17]).
2.1.1. Spatial schemes

OpenFOAM is based on finite volume schemes that use the center x⃗P of the cell to define the variation in space of
a scalar function φ = φ(⃗x ). In order to yield second order accurate schemes, this variation needs to be linear inside a
control volume V p (see [11,17]). So,
φ(⃗
x ) = φ(⃗
x P ) + (⃗
x − x⃗P ) · (∇φ)(⃗ x − x⃗P |2 )
x P ) + O(|⃗ (3)
where · is the inner product.
In fact, considering P = x⃗P as the centroid of V P to perform a volume integral on Eq. (3), the integral of the linear
term on the right-hand side disappears and a second order integral approximation is obtained:
∫
φ(⃗
x )dτ = φ P V P + O(|⃗
x − x⃗P |2 ) (4)
Vp
where φ P = φ(⃗x P ).
Similarly, considering a scalar quantity φC f evaluated at the center of a face (Fig. 1a), we can obtain:
∫
φ(⃗
x )d S⃗ = φC f S⃗C f + O(|⃗
x − x⃗C f |2 ) (5)
Sf
(a) Orthogonal mesh. (b) Mesh non-orthogonality.
(c) Mesh skewness. (d) Mesh non-uniformity.
Fig. 1. Mesh properties.
in which S⃗C f indicates the surface area vector S⃗ evaluated at the center C f of the face f , and the magnitude of SC f is
the area of the face.
Along a face f , a scalar property φC f can be evaluated by linear interpolation of φ using information from adjacent
elements stored on nodes P = x⃗P and N = x⃗N (see Fig. 1):
φC f = αφ P + (1 − α)φ N (6)
where α is an interpolation factor that can be defined as the ratio of distances C f N over |d|
⃗ = P N , and so,
Cf N
α= . (7)
PN
Alternatively, we can use midpoint interpolation, that is, a symmetric weighting (α = 0.5) interpolation from
nodes P and N , or a cubic interpolation derived from the linear interpolation, with the addition of an explicit higher
correction. This correction is calculated using the gradient values (∇φ) P and (∇φ) N .
If the control volume V p is bounded by a control surface composed of flat faces f , the divergence terms can be
approximated using the Gauss–Ostrowski theorem (see [11,17]). The divergence term of the vector property v⃗ can be
approximated as follows:
∫ ∮ ∑∫ ∑
(∇ · v⃗)dτ = d S · v⃗ =
⃗ d S⃗ · v⃗ ≈ S⃗C f · v⃗C f (8)
Vp Sp f Sf f
where the components v⃗C f can be evaluated using the interpolation Eq. (6).
For a scalar, the gradient term can be obtained using either the divergence theorem or a least-squares fit
(see [11,17]). Discretization based on the divergence theorem can be written as follows:
∫ ∮ ∑∫ ∑
∇φdτ = ⃗ =
d Sφ ⃗ ≈
d Sφ S⃗C f φC f (9)
Vp Sp f Sf f
where φC f is evaluated using the interpolation Eq. (6).

Taking into account a linear variation of φ over the control volume P, the cell-centered approach for the gradient
at P is calculated as follows:
1 ∑⃗
(∇φ) P = S f φC f . (10)
Vp f
This discretization generates the OpenFOAM Gauss scheme for the gradient.
The least squares fit approach is based on minimizing the following functional with respect to the gradient (∇φ) P
(see [17]):
( )2
∑ φN − φP d⃗ · (∇φ) P
L(∇φ P ) = − . (11)
⃗
|d| ⃗
|d|
f
The numerical expression φ N|−φ⃗

d|
P
in Eq. (11) represents a linear approximation of the gradient (∇φ) P projected
on vector d⃗ (see Fig. 1b). This approximation generates a second order least squares scheme for the gradient. The
optimum gradient (∇φ) P , defined as the gradient that minimizes Eq. (11), can be found by solving the system of
linear equations generated by the minimizing condition:
dL
= 0. (12)
d(∇φ) P
The diffusion term is discretized using Eq. (8) (see [17]):
∫ ∮ ∑∫ ∑ ( )
∇ · (ν∇φ) dτ = d S⃗ · (ν∇φ) = d S⃗ · (ν∇φ) ≈ S⃗C f · νC f ∇φ C f (13)
Vp Sp f Sf f
where ν represents the diffusivity (considered constant in this work). This discretization defines the Gauss Laplacian
scheme, which is the only means available in OpenFOAM for discretizing the Laplacian operator.
In our case, the source term q = q(⃗ x ) on the right-hand side of Eq. (1) depends only on spatial coordinates, and is
obtained by implementing a manufactured solution for verification purposes. A second order integral approximation
for the control volume V P is obtained as follows:
∫
x )dτ = q P V P + O(|⃗
q(⃗ x − x⃗P |2 ) (14)
Vp
where q P = q(⃗ x P ).
The second order discretized expression for the Poisson equation (2) can now be written as:
∑
S⃗C f · (ν∇φ)C f = q P V P . (15)
f
The gradient (∇φ)C f in Eq. (15) is expressed through the values of φ P and the values of φ in the surrounding cells.
Details of the evaluation of the left hand side of Eq. (15) are given in Section 2.3.
Accordingly, Eq. (15) leads to an algebraic equation (see [11]) for cell P:
∑
aP φP + aN φN = R P (16)
N
where a P , a N and R P are constants produced by the discretization process.

Assembling Eq. (16) for all cells centers P, a system of algebraic equations is obtained:
[A] [φ] = [R] (17)

which is solved iteratively.
2.2. The meshes
Mesh shape and size are key factors determining the accuracy of a numerical approximation. Juretic [17] defines
three properties that affect mesh quality:
• Non-orthogonality. An orthogonal mesh can be defined as a mesh for which the vector d⃗ joining the cell center
P and the center N of an adjacent cell is parallel to the vector S⃗ normal to a cell face f (see Fig. 1a). A mesh
that does not satisfy the orthogonality condition is defined as non-orthogonal (see Fig. 1b).
• Mesh skewness. When the segment P N does not intercept the centroid point C f of the face f , the mesh is
defined as skewed (see Fig. 1c).
(a) Hexahedral orthogonal. (b) Equilateral orthogonal. (c) Stretched orthogonal.
(d) Regular left prism. (e) Regular alternate prism. (f) Delaunay prism.
(g) Hexahedral slanted. (h) Hexahedral slanted stretched.
Fig. 2. Orthogonal and non-orthogonal prism meshes in 2D representation.
• Non-uniformity. A mesh is uniform when d⃗ intersects the face midway between the nodes P and N . A mesh
that does not satisfy the uniformity condition is defined as non-uniform (see Fig. 1d).
We have chosen a variety of meshes for the verification procedure to enable us to test these three properties. After
performing some basic tests, we will focus on the non-orthogonality of the grid for which a loss of convergence order
is observed. We use orthogonal meshes (shown in Figs. 2a, 2b, and 2c), and non-orthogonal prism meshes (shown
in Figs. 2d, 2e, and 2f). Non-orthogonal hexahedral meshes are depicted in Figs. 2g and 2h. A physical domain
[0,π ][0,π ] is used for the square domains. For the slanted domains (parallelograms), the upper boundary is moved
to the right by approximately 16% of the characteristic length, and the rhombus domain is chosen in order to obtain
equilateral prisms.
2.3. Treatment of non-orthogonality
For non-orthogonal grids (see Fig. 3), the diffusive flux on a given face f has more than one component, as is the
case for orthogonal grids. Several multi-point face schemes have been presented by Schafer [25], however using these
schemes may increase the complexity of the computations and cause difficulties in converging the solvers.
To deal with non-orthogonality, a one-point centroid-based scheme has been proposed by Jasak [11]. This scheme
uses a first term that is calculated using the cell centers P and N neighboring the face f (see Fig. 1b). This term is
Fig. 3. Non-orthogonality treatment.
Fig. 4. Non-orthogonality treatment: “minimum correction” approach.
called the orthogonal contribution. A second term is then added to provide the gradient projection of the outward-
pointing face area vector S⃗ f with a “correction”. This is called the non-orthogonal contribution. The mathematical
expression of this idea (see Eq. (18)) is:
⃗ φ N − φ P + k⃗ · (∇φ) f
S⃗ · (∇φ) f = |∆| (18)
⃗
|d|
⃗ φ N −φ P the
where f denotes the face of a control volume, (∇φ) f the gradient of a physical property φ at face f , |∆| ⃗
|d|
orthogonal contribution and (k⃗ · (∇φ) f ) the non-orthogonal correction.
The non-orthogonal contribution term in Eq. (18) corresponds to a deferred correction approach. As explained
in [9], deferred correction approaches are used in FVM to, among other things, define higher-order schemes.
Muzaferija [21] suggested a second order deferred correction scheme to deal with non-orthogonality that uses a
correction vector.
Among a set of possible non-orthogonal decomposition contributions in Eq. (18), Jasak [11] examines three, which
are called the “minimum correction”, “orthogonal correction”, and “over-relaxed” approaches (see [12]). These three
ways to compute the non-orthogonal correction are illustrated in Figs. 4, 5, and 6 respectively. What these three
approaches have in common is that they decompose the surface normal vector S⃗ at the center of face f as follows
(see [11]):
⃗ + k.
S⃗f = ∆ ⃗ (19)
⃗ is specifically calculated for each of the three approaches:
The vector ∆
• For the “minimum correction” (see Fig.( 4)), the calculation is performed in such a way as to keep the non-
orthogonal correction as small as possible:
⃗ = cos(α N )|S|d⃗u
∆ (20)
P⃗N
where vector d⃗u = ⃗ |S| is the area of face f , and cos(α N ) is the cosine of the angle
is the unit vector d,
| P⃗N |
S⃗
between the unitary vector d⃗u and the normal face vector n⃗ = |S| .
• For the “orthogonal correction” (see Fig.( 5)), the contribution from φ P and φ N is kept the same as on the
orthogonal mesh:
⃗ = |S|d⃗u .
∆ (21)
Fig. 5. Non-orthogonality treatment: “orthogonal correction” approach.
Fig. 6. Non-orthogonality treatment: “over-relaxed” approach.
• In the “over-relaxed approach” (see Fig. 6), the importance of the terms φ P and φ N increases with the increase
in non-orthogonality:
|S| ⃗
⃗ =
∆ du . (22)
cos(α N )
The vector k⃗ is calculated from the equation:
⃗
k⃗ = S⃗f − ∆. (23)
Finally, note that the scheme presented in Eq. (18) is a non-linear equation owing to the non-orthogonal
contribution. The implementation of this scheme needs the value of the gradient in the non-orthogonal contribution.
This is achieved by means of a recursive process in which the non-orthogonal contribution is explicitly calculated
on the right-hand side of the matrix of the algebraic system of equations (see Eq. (17)). Another important aspect
is the limitation of the non-orthogonal contribution, which can slow down the convergence of the iterative method
solving the matrix system. All three methods shown by Jasak for computing k⃗ yield a constraint on such a vector.
An alternative that will ease this problem is to apply a limiter to the non-orthogonal contribution. In OpenFOAM, a
constraint is added to the over-relaxed approach in order to prevent the non-orthogonal part from becoming dominant,
as we will see later in this work (see Eq. (32)).
3. Methodology
In this section, we begin by recalling the method of manufactured solutions (MMS), which is a powerful verification
tool for determining the convergence rate through the calculation of exact numerical error using mesh refinement.
Next, we present the norms used to obtain the average and maximum errors. Finally, we present the mathematical
expression for calculating the order of accuracy in the logarithmic scale.
3.1. Method of manufactured solution
The general idea behind MMS is to use a custom built source term that allows Eq. (1) to have an analytic solution.
The steps of the MMS methodology can be described as follows (see [19]):
• Propose a manufactured solution for a system of equations.
• Obtain a source term for the initial equations.
• Implement that source term in the solver.
• Solve the new equation.
• Compute the error.
In the present work, a scalar manufactured solution is used to verify the Poisson equation when using Dirichlet and
Periodic (Cyclic) boundaries:
φ(x, y) = φ0 [cos(ωx) sin(ωy)] (24)

where the frequency ω and the amplitude φ0 of the sinusoid are set to 1 for Dirichlet boundaries. Periodic boundaries
are easily tested if ω = 2.0 is set to ensure the periodicity of the solution.
The advantage of this choice of solution is that the source term obtained is multiplied by ω2 , which is a product
that can be used to avoid convergence problems in the explicit part of the matrix system equation (see Eq. (17)).
Other manufactured solutions will also be used in this work to test the Neumann boundary condition. The source term
produced by the manufactured solution for Poisson’s Eq. (1) is:
q(x, y) = 2νω2 φ(x, y). (25)
3.2. Error estimation
The exact error ei at the cell center “i” is the difference between the analytical manufactured solution φ̃ and the
numerical solution φh :
ei = φ̃ − φh . (26)
To calculate the order of convergence, we use two types of norm. The first, for the average error, is calculated as
follows:
∫ NC
1 ∑
∥e∥ L 1 (Ω) = |e| dΩ ≈ Ωi |ei | (27)
Ω |Ω |
i=1
where Ω indicates the domain definition, Ωi is the control volume of cell i, and N C indicates the number of cells. The
second norm, for the maximum error, is calculated as follows:
∥e∥∞ = sup |e(x)| ≈ max |ei | (28)
x∈Ω i
where i = 1, . . . , N C.
Among these two norms, the infinite norm is generally the most sensitive, which justifies its use. Other norms for
calculating the average error, for instance the Euclidean L 2 norm, will not be used in the present verification process.
The choice of the L 1 norm can be justified by the inequality :
∥e∥ L 2 ≤ ∥e∥ L 1 (29)
where the convergence of the L 1 norm implies the convergence of L 2 norm. Indeed, to justify Eq. (29) we can consider
the inequality (in a discrete finite space):
[ n ] 1p [ n ]
∑ ∑
p
∥e∥ L p = |ei | ≤ |ei | = ∥e∥ L 1 (30)
i=1 i=1
where p = 2, 3, . . ..
3.3. Convergence rate
The correctness of a computer code is verified via systematic mesh refinement. In our case, we uniformly
divide a parallelogram domain, with Nk being the number of nodes along each boundary, taking the values
Nk = [10, 20, 40, 60, 80, 100]. This sequence refines the control cell at each step for all meshes except for the
Delaunay mesh (see Fig. 2f). For this mesh, the refinement is applied only on boundary faces, after which the cells are
automatically constructed using a Delaunay algorithm without a strict guarantee that the cells will be systematically
refined by a constant factor inside the domain. Strictly speaking, the simulations using this mesh only illustrate the
trend of convergence when the boundary is refined, and do not constitute a measure of convergence based on the
systematic internal refinement of cells.
If ek represents the error produced by partitioning the domain Nk times (with [10, 20, 40, 60, 80, 100]), conver-
gence is defined as:
log (ek /ek−1 )
Rk = . (31)
log (Nk /Nk−1 )
3.4. Boundary conditions (BCs)
In this work, the verification procedure is only performed for three basic types of BC: the Dirichlet, Neumann, and
Periodic boundaries. The convergence rate issues revealed with respect to the Dirichlet and Neumann BCs explain
why we initially restrict the analysis to these types of BC.
3.5. OpenFOAM spatial schemes for the diffusion operator
In this section, the OpenFOAM implementation of the Laplacian scheme for the diffusion operator is described.
This scheme is based on the expression on the right-hand side of Eq. (13) for which an interpolation, a gradient, and
a non-orthogonal scheme are required. A number of options for each of these schemes are presented below.
As emphasized in the OpenFOAM documentation (see [28,22]) the “Gauss scheme” is the only choice for the
Laplacian operator. The Gauss scheme syntax is composed of the name “Gauss”, an interpolation scheme for the
diffusion coefficient ν considered constant in this study, and a selection of a limiter to the surface normal gradient for
the non-orthogonal treatment:
Gauss <interpolationScheme> <snGradScheme>
The many options for the surface normal gradient scheme, snGradScheme, are the following (see [28,22]):
• corrected — unbounded, second order, conservative
• uncorrected — bounded, first order, non-conservative
• fourth — fourth order, conservative
• limited ψ — blend of corrected and uncorrected
where the parameter ψ limits the non-orthogonal part. Typical values are:
• ψ = 1.0 corresponding to a corrected, unbounded, second order, conservative scheme
• ψ = 0.5 corresponding to a non-orthogonal correction ≤ of the orthogonal part
• ψ = 0.333 corresponding to a non-orthogonal correction ≤ 12 of the orthogonal part
• ψ = 0.0 corresponding to an uncorrected, bounded, first order, non-conservative scheme.
The role of the parameter ψ is to serve as a constraint to the non-orthogonal part, as follows:
( )
λ k⃗ · ∇φ (32)
where λ is calculated as:
( )
(ψ)|orthogonalpart|
λ = min , 1 (33)
(1 − ψ)|non-orthogonalcorrection| + ϵ
where ϵ is a “small” value.
The many schemes of the gradient operator, which is part of the non-orthogonal treatment, are: Gauss, least squares,
and fourth, and include cell-limited and face-limited versions of each. The Gauss and least-squares schemes are based
on Eqs. (10) and (11) respectively. The fourth scheme is a fourth order least squares approach. Following is a brief
description of each of these (see [28,22]):
• Gauss second order, based on Gaussian integration

• Least Squares second order, based on the least-squares approach
• Fourth, fourth order least squares approach.
OpenFOAM offers three main schemes to interpolate the values from cell centers to face centers: “linear”, based
on Eq. (6); “midpoint” which is linear with symmetrical weights α = 0.5 and “cubic”, which is a cubic scheme
implemented as a correction to the linear scheme.
Before ending this section, we wish to point out that the verification tests presented in this work were performed
on OpenFOAM-2.3.x and Foam-Extend-3.1, which serve as references for other OpenFOAM versions.
4. Verification results
This section reports on the convergence tests that make it possible to detect and track the origin of the loss of
convergence accuracy. First, we comment on the verification procedure for orthogonal meshes. Then, we verify the
convergence rate on various non-orthogonal meshes for different BCs. Finally, we show the verification strategy for
detecting the origin of the convergence loss when non-orthogonal meshes are used.
In the figures presented in this section, the error is represented as a function of the refinement steps on a logarithmic
scale for which the slope of functions represents the convergence order (see Section 3.3). The labels Avg. and Max.
indicate that the average and the maximum error norms respectively are used to calculate the convergence rate for
the numerical solutions. To perform the tests on 2D meshes, we use a single layer of cells with empty BCs in the
front and back. Note that, unless explicitly stated otherwise, the manufactured solution illustrated by Eq. (24) is used
hereafter on a square domain [0; π ][0; π ], applying the Dirichlet BC along the boundaries. The diffusion coefficient
is supposed to be constant and its interpolation scheme is not tested.
4.1. Initial test verification
The strategy for verifying the convection–diffusion operators follows the verification methodology proposed by
Knupp and Salari [19] to test Burger’s equation using the Dirichlet and Neumann BCs. In this way, an initial set
of tests for Burger’s equation is successfully performed in terms of the theoretical convergence rate on orthogonal
meshes (see Figs. 2a–2c). In these tests, the Dirichlet BC is first applied on all the boundaries, and then simulations are
performed combining the Dirichlet and Neumann BCs on alternating boundaries. Theoretical second order accuracy
is obtained for all the BC combinations. These results are not illustrated in this work.
This first set of verification tests is also performed for Left non-orthogonal regular prism meshes (see Fig. 2d) using
second order schemes for Burger’s equation. A reduction to first-order accuracy is observed when only the Dirichlet
BC is used. For this reason, to simplify the problem, the initial goal of verifying the convection–diffusion operators
focuses on verifying the diffusion term via a Poisson equation.
4.2. Verification of the Poisson equation on orthogonal meshes for the Dirichlet BC
In this section, a set of verification tests for Poisson’s equation on orthogonal meshes using the Dirichlet BC is
presented. To perform these tests, the orthogonal meshes presented in Figs. 2a–2c) √ are used.
√ The first is composed of
regular hexahedra, the second of equilateral prisms (for this, the domain [−1; 1][− 3; 3] is used), and the third of
stretched hexahedra with a ratio of 4:1 (i.e. each edge of the bottom-left cell is one-quarter the length of each side of
the upper-right cell). The Gauss linear orthogonal scheme is used for the Laplacian.
In Fig. 7, the convergence rates for the Avg. and Max. errors on hexahedral (Hex.), hexahedral stretched (Stret.),
and equilateral (Equil.) prism meshes are presented. The six curves in this plot have a slope of 2, which corresponds
to the predicted theoretical second order accuracy. The results obtained for the Stret. meshes verify that, for these
simulations, the convergence is not impacted by non-uniformity. Following this step, tests on elongated orthogonal
Fig. 7. Curves on orthogonal meshes.
Fig. 8. Curves on orthogonal elongated meshes.
meshes are presented (on domain [−0.4 : 0.4]; [−0.3 : 0.4]), first with a refinement by a factor of 10 in the x-
direction using 40 × 4, 80 × 8, 160 × 15, 320 × 32, 640 × 64, and 1280 × 128 partitions, and then with refinement
by a factor of 10 in the y-direction. Although the Hex. meshes are strongly elongated in one direction, the theoretical
convergence rate of 2 is preserved for the Avg. and Max. errors, as observed in Fig. 8. In this graph, [10:1] means 10
times more refined in the x-direction. Second-order accuracy is also observed in simulations combining the Dirichlet
and Neumann BCs on alternating boundaries on square domains, however those results are not illustrated here.
Fig. 9. Curves on non-orthogonal meshes with varying limited parameter ψ.
4.3. Grid refinement results on non-orthogonal meshes
On non-orthogonal meshes, we test all combinations of the scheme options supported for the Laplacian scheme,
as explained in Section 3.5: the interpolation, gradient, and non-orthogonal schemes. For this reason, the verification
process is divided into several stages. Next, we present the most representative results. For these simulations, unless
otherwise specified, the Gauss linear scheme is used for the gradient appearing in the non-orthogonal contribution.
4.3.1. Testing the OpenFOAM Laplacian non-orthogonal constraint using the Dirichlet BC
In the first step, the parameter ψ = [1.0, 0.5, 0.333, 0.0], explained in Section 3.5, for the non-orthogonal treatment
is tested. The verification tests are performed using only a Left non-orthogonal regular prism mesh, as presented in
Fig. 2d.
Significant differences in the results are found for the various values of ψ, and the best performance in terms of
convergence rate is obtained for the value ψ = 1.0, but, unfortunately, this is only a first-order convergence. As can
be appreciated in Fig. 9, the Left Avg. ψ = 1.0 and Left Avg. ψ = 0.5 curves are very close in terms of convergence
and accuracy, as a convergence rate of 1.0 is observed. The Left Avg. ψ = 0.333 curve has an initial slope of 1.0,
but, when refined, its convergence rate decreases considerably. The curve with non-orthogonal correction, Left Avg.
ψ = 0.0, has a zero convergence rate. In these plots, first-order accuracy is also observed for the maximum error in
the Left Max 1.0 curve (using ψ = 1.0), but this latter has less precision relative to that of the Left Avg. 1.0 curve.
When using Gauss cubic and Gauss midpoint instead of Gauss linear for the gradient reconstruction implicitly
appearing in the non-orthogonal contribution, no improvement in convergence is observed. Neither is there any
improvement in convergence when other gradient schemes are used, such as Fourth and Least Squares, including
their CellLimited and FaceLimited versions. Some of these results are presented in Fig. 10. In this plot, the Gauss,
Least Squares, and Fourth labels indicate the scheme used for the gradient inside the non-orthogonal term. The labels
linear and cubic indicate the interpolation scheme used for the Gauss scheme. The value for ψ indicates the value in
the limited scheme. As can be appreciated in Fig. 10, these curves present a convergence rate of 1. In this graph, we
can also compare some curves when the Least Squares scheme is used. The Least Squares curves with ψ = 1.0 and
ψ = 0.5 present a convergence rate of 1.0. For the least squares curve with ψ = 0.333, an asymptotic convergence
rate of less than 1.0 is observed.
Unless otherwise explicitly stated, the verification tests are hereafter conducted using the Gauss linear limited
scheme 1.0 for the Laplacian.
Fig. 10. Curves on non-orthogonal meshes using a set of gradient schemes for non-orthogonal treatment.
4.3.2. Results on different types of non-orthogonal triangular mesh

In this section, the influence of non-orthogonality is illustrated for a variety of non-orthogonal meshes that include
a Left regular prism (see Fig. 2d), an Alternate regular prism (see Fig. 2e); a Delaunay regular prism (see Fig. 2f); a
Slanted hexahedral (see Fig. 2g), and a Slanted stretched hexahedral (see Fig. 2h) mesh.
In Fig. 11, we see some convergence curves for the average error and one curve for the maximum error. As
presented in this figure, an asymptotic convergence rate of 1 is observed in the Left Avg. curve. The convergence rate
for the Alternate Avg. curve is the expected theoretical value of 2.0, but this value is only 1.0 for the Alternate Max.
curve. That is because the Alternate regular mesh keeps the orthogonality inside the domain, but this property is not
preserved for the cells neighboring the boundaries (see Fig. 2e). With this result, we are able to detect the origins of
non-orthogonality issues. It is also remarkable that the Alternate mesh verifies the skewness inside the domain. So,
it seems that neither uniformity nor skewness participates in the convergence issues, which are likely caused only by
non-orthogonal treatment on boundaries.
The Delaunay Avg. curve presents a critical loss in convergence rate for the more refined meshes. As explained in
Section 3.3, the simulations using this mesh only illustrate the convergence trend when the boundary is refined. Finally,
the Slanted Avg. curve and the Stretched Avg. curve, using hexahedral meshes, present an asymptotic convergence
accuracy of 1.0 for the average error, which is in agreement with the trend of the previous non-orthogonal meshes.
4.3.3. Testing non-orthogonality that systematically inclines the border

To illustrate the influence of non-orthogonality on the convergence order loss, the initial square domain [0; π ][0; π ]
is systematically inclined moving the upper boundary to the right. Initially, the upper boundary is moved one percent
of the characteristic distance L (i.e., ∆x = 0.01L) to the right, then this process is repeated taking ∆x = n(0.01L),
with n = 2, 3, 4, 5.
The results for average error are presented in Fig. 12. In this representation, the index 0,1,2,3,4,5 corresponds
to the number of steps moved, i.e. Avg-0 corresponds to the initial orthogonal mesh, Avg-1 is obtained from the
initial orthogonal moving the top boundary one step, etc. As seen in Fig. 12, the convergence rate, which is initially
2 (as we have seen for the Avg-0 curve) decreases when the number of steps is increased, and this convergence
loss is strengthened as non-orthogonality increases. These results clearly indicate that convergence is influenced by
non-orthogonality.
4.3.4. Non-orthogonal mesh with a periodic (cyclic) BC

In this section, simulations with a periodic BC are performed. Issues with these types of BC have not been found,
either for orthogonal or non-orthogonal meshes. Two pairs of neighboring periodic BCs are used (left and right, and
Fig. 11. Curves on non-orthogonal meshes.
Fig. 12. Incidence on non-orthogonal meshes inclining an initial square domain.
up and down boundaries). An Orthogonal hexahedral mesh, and two non-orthogonal Left and Alternate prism meshes,
which are presented in Figs. 2a, 2d, and 2e respectively, are used. The manufactured solution specified by Eq. (24) is
implemented, using the space frequency ω = 2.0 to guarantee continuity of the solution.
Some results are shown in Fig. 13. In this representation, the Orthogonal Avg. and Orthogonal Max. curves
correspond to the average and maximum error convergence on orthogonal meshes. Both curves have a slope of 2.
The same trend is observed for some curves on a Left regular and an Alternate regular prism mesh for the average and
maximum errors, as we have seen in the Left Avg. ψ = 1.0, Left Max. ψ = 1.0, and Alternate Avg. ψ = 1.0 curves
(note that the Orthogonal Avg. and Alternate Avg. ψ = 1.0 curves are superposed. The value ψ refers to the parameter
Fig. 13. Curves using periodic (cyclic) BCs.
used for the Laplacian scheme). Good convergence is also observed for ψ = 0.5 on Left regular meshes. However,
as we observed in other verification tests, a poor convergence rate is observed on Left meshes for the average and
maximum error (this latter result is not presented here) when ψ = 0.333, as can be seen in the Left Avg. ψ = 0.333
curve. These results on a Periodic BC confirm that the convergence issues originated at the boundaries.
4.4. Convergence rate for the gradient
Generally speaking, the trend seen up to now for the convergence rate of the basic variable φ in Poisson’s equation
is preserved for the gradient (∇φ): using theoretically second order schemes, the convergence rate of 2.0 is observed
for both gradient components [∂φ/∂ x, ∂φ/∂ y] on orthogonal meshes, whereas that on non-orthogonal meshes is 1.0
for both components. The numerical gradient is calculated with the help of the OpenFOAM libraries fvc::grad on an
Orthogonal regular hexahedral mesh and on a Left regular prism mesh. The gradient reconstruction is performed using
the Gauss Linear, Least Squares, and Fourth schemes, but the results are presented only for the Gauss Linear scheme.
The Fourth and Least Squares schemes do not show any significant improvement in terms of convergence behavior.
In Fig. 14, we can see the curves with labels Avg. ∂φ/∂ x Ortho, Avg. ∂φ/∂ y Ortho, Max. ∂φ/∂ x Ortho, and Max.
∂φ/∂ y Ortho for the two gradient components for the average and maximum errors on orthogonal meshes. For these,
we observe a theoretical convergence rate of 2.0. For non-orthogonal meshes, the curves with labels Avg. ∂φ/∂ x
nonOrtho and Avg. ∂φ/∂ y nonOrtho for the average error have an asymptotic convergence rate of 1. By contrast, the
Max. ∂φ/∂ x nonOrtho and Max. ∂φ/∂ y nonOrtho curves have an asymptotic convergence rate of 0, which originates
from the boundaries.
4.5. Implementing a non-orthogonal treatment at the BC
A few factors guide us in identifying the source of convergence reduction on non-orthogonal meshes: theoretical
convergence rates are observed in simulations using a Periodic BC and in simulations using Alternate meshes.
Convergence can be improved by using an external treatment, i.e. by not modifying OpenFOAM FVM libraries,
and this approach is based on the use of the manufactured solution. We use the fact that the non-orthogonal term is
implemented explicitly on the right-hand side of matrix equations (Eq. (17)) as a source term, and so this correction
can be implemented as a source term in the Poisson equation using the manufactured gradient. For the non-orthogonal
term k⃗ · (∇φ) f , both the non-orthogonal correction vector k⃗ and the gradient (∇φ) f need to be reconstructed. For
the latter, the gradient produced by the MMS can be used directly. We only perform a reconstruction on boundary
Fig. 14. Curves on non-orthogonal meshes for the gradient of φ.
faces for these two fields, so this implementation does not interfere in any way with the non-orthogonal scheme for
the Laplacian inside the domain. For the Surface Field vector k, ⃗ an over-relaxed correction approach at the BC is
computed. By implementing k⃗ · (∇φ) f on the laplacianFoam solver, an immediate improvement in the convergence
rate is obtained.
The results using the Left regular prism, the Alternate regular prism, and the Delaunay prism meshes (see Figs. 2d,
2e, and 2f respectively) are presented. In Fig. 15, some curves for the average error Avg. are shown. In the Left Avg.,
Alternate Avg., and Delaunay Avg. curves, the surface gradient (∇φ) f is replaced by the analytical gradient evaluated
at the faces. For these curves, we obtain the expected theoretical convergence rate of 2.0.
4.6. Non-orthogonal mesh with Neumann BC
In this section, we consider the Dirichlet BC combined with the Neumann BC (in the upper and in the lower
boundaries), alternating opposite sides in a square domain [0; π ][0; π ]. We verify two ways of handling the Neumann
BC: zeroGradient and fixedGradient, which are separate classes in OpenFOAM. For zeroGradient, the projection of
the gradient along the normal to a boundary face is 0.
First, we test zeroGradient using the following manufactured solution. The projection along the normal on the
upper and lower boundaries is nil:
[( )]
1 − x 2 2y 3 − 3π y 2 + π 3
)(
φ(x, y) = . (34)
π3
Some results using a Left non-orthogonal mesh type (see Fig. 2d) are shown in Fig. 16. This representation shows
the curves labeled Left Avg. and Left Max., which show a convergence in the order of 1.0 (using the original code).
The plots Left Avg. Corr. and Left Max. Corr. were obtained by applying the non-orthogonal analytical correction to
the Dirichlet BC on the left and right boundaries. An improvement in accuracy and convergence rate can clearly be
seen. However, both curves display an asymptotic convergence in the order of 1.0. This makes it possible to conclude
that the zeroGradient condition is the cause of a lower order of convergence.
Second, we applied a fixedGradient BC on the upper and lower boundaries using the following manufactured
solution:
)( 3 )⎤
1 − x 2 y3 − π2 y
⎡( 2
φ(x, y) = ⎣ ⎦. (35)
π2
Fig. 15. Curves using an analytical gradient on the non-orthogonal correction.
Fig. 16. Curves using Dirichlet and adiabatic BCs.
The results are shown in Fig. 17. Once again, the Left Avg. and Left Max. curves obtained with non-orthogonal
grids of the Left type, show a convergence rate of the first order only (using the original code). The Left Avg. Corr.
and Left Max. Corr. plots with non-orthogonal analytical correction at the Dirichlet BC show better accuracy and
convergence; however, the asymptotic convergence rate is only 1.0.
Fig. 17. Curves using Dirichlet and Neumann BCs.
5. Conclusions
A manufactured solution procedure has been applied to verify the OpenFOAM diffusion operator on different types
of mesh. To achieve this, a manufactured Poisson’s equation was solved using OpenFOAM schemes. Convergence
analysis was performed by applying systematic grid refinement. Predicted second order accuracy was found when
orthogonal grids and various alternatives for handling BCs (Dirichlet, Neumann, and periodic) are used.
On non-orthogonal grids, convergence deterioration was observed for several OpenFOAM schemes. It was found
that even minor non-orthogonality significantly reduces the convergence order. It was also discovered that, in the
domain’s core, second-order accuracy is maintained, and that convergence reduction essentially derives from boundary
treatment.
The verification procedure brought to light two factors which enable us to detect the origin of convergence
reduction. First, it was found that when periodic BCs are used, second-order accuracy is obtained for all types of
grid, orthogonal and non orthogonal. Second, a convergence loss was found for grids of the Alternate type, but only in
the neighborhood of the boundary (In fact, the theoretical second order convergence is preserved in the domain). Thus,
convergence reduction derives from non-orthogonal treatment at boundaries using the Dirichlet BC. A strategy for the
treatment of non-orthogonality has been proposed by taking into account the analytical gradient (i.e. the gradient
calculated from the manufactured solution).
Besides, on the one hand, when applying a combination of the Neumann BC and improved treatment of the
Dirichlet BC, it was found that the BC zeroGradient and fixedGradient implementations also lead to a convergence
degradation for non-orthogonal boundaries. On the other hand, it was found that the presence of skewness and non-
uniformity seems to be well handled.
Another factor that has a significant impact on convergence reduction, and which does not originate from the
boundaries, involves the choice of the constraint of the non-orthogonal correction scheme. In particular, the Gauss
limited ψ scheme shows the best results for the value ψ = 1.0 for simulations using a Laplacian solver.
Alternatives for reconstructing numerical schemes on boundaries for second-order convergence will be presented
in the second part of this work. The well-founded nature of the proposed options will be verified by performing
convergence tests on problems described by the Laplace and Navier–Stokes equations.
This work has revealed the need for a verification of the numerical algorithms, as these are strongly influenced by
the kind of mesh involved. It has been shown that verification is a complex task, even though it is only applied to
the diffusion operator when a simple solver like laplacianFOAM is used, and which has been applied on prisms and
regular hexahedrons with the Dirichlet, Neumann, and Periodic boundary conditions.
References
[1] J. Abanto, D. Pelletier, A. Garon, J. Trepanier, M. Reggio, Verification of some commercial CFD codes on atypical CFD problems, in: 43rd
AIAA Aerosp. Sci. Meet. and Exhib., Vol. 1, 2005, pp. 15693–15725.
[2] J. Blazek, Computational Fluid Dynamics: Principles and Applications, second ed., Elsevier Sci., 2006.
[3] B. Diskin, J. Thomas, Comparison of node-centered and cell-centered unstructured finite-volume discretizations inviscid fluxes, AIAA J.
49 (4) (2011) 836–854.
[4] B. Diskin, J. Thomas, Effects of mesh regularity on accuracy of finite-volume schemes, in: 50th AIAA Aerosp. Sci. Meet., 2012.
[5] B. Diskin, J. Thomas, E. Nielsen, H. Nishikawa, J. White, Comparison of node-centered and cell-centered unstructured finite-volume
discretizations: viscous fluxes, AIAA J. 48 (7) (2010) 1326–1338.
[6] L. Eca, M. Hoekstra, An introduction to CFD code verification including eddy-viscosity Models, in: Eur. Conf. on Comput. Fluid Dyn.,
ECCOMAS CFD, 2006.
[7] L. Eca, M. Hoekstra, Code verification of unsteady flow solvers with the method of the manufactured solutions, in: Int. offShore and Polar
Eng., ISOPE 2007, 2007, pp. 2012–2019.
[8] L. Eca, M. Hoekstra, A. Hay, D. Pelletier, Verification of RANS solvers with manufactured solutions’, Eng. Comput. 23 (2007) 253–270.
[9] J. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, third ed., Springer, 2002.
[10] J. Iannelli, An exact non-linear NavierStokes compressible-flow solution for CFD code verification, Internat. J. Numer. Methods Fluids 72 (2)
(2013) 157–176.
[11] H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows (Ph.D. thesis), Imperial College,
1996.
[12] H. Jasak, A. Gosman, Automatic resolution control for the finite-volume method, part 1: a-posteriori error estimates, Numer. Heat Transf.
Part B: Fundam. 38 (3) (2000) 237–256.
[13] H. Jasak, A. Gosman, Automatic resolution control for the finite-volume method, part 2: adaptive mesh refinement and coarsening, Numer.
Heat Transf. Part B Fundam. 38 (3) (2000) 257–271.
[14] H. Jasak, A. Gosman, Automatic resolution control for the finite-volume method, Part 3: turbulent flow applications, Numer. Heat Transf. Part
B Fundam. 38 (3) (2000) 273–290.
[15] H. Jasak, A. Gosman, Rsidual error estimate for the finite volume method, Numer. Heat Transf. Part B: Fundam. 39 (1) (2001) 1–19.
[16] H. Jasak, A. Gosman, Element residual error estimate for the finite volume method, Comput. & Fluids (2003) 223–248.
[17] F. Juretic, Error Analysis in Finite Volume CFD (Ph.D. thesis), Department of Mechanical Engineering, Imperial College, 2004.
[18] F. Juretic, A. Gosman, Error analysis of the finite-volume method with respect to mesh type, Numer. Heat Transf. Int. J. Comput. Methodol.
57 (6) (2010) 414–439.
[19] P. Knupp, K. Salari, Verification of Computer Codes in Computational Science and Engineering, Chapman and Hall/CRC, 2002.
[20] F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamic, Springer, 2015.
[21] S. Muzaferija, Adaptive Finite Volume Method for Flow Predictions using Unstructured Meshes and Multigrid Approach (Ph.D. thesis),
University of London, 1994.
[22] Programmer-Guide, OpenFOAM. The Open Source CFD ToolBox. Programmers Guide. Version 2.2. September 2013. OpenFOAM-
Foundation, 2013.
[23] P. Roache, Verification of codes and calculations, AIAA J. 36, 5 (5) (1998) 696–702.
[24] C.J. Roy, C.C. Nelson, T.M. Smith, C.C. Ober, Verification of Euler/NavierStokes codes using the method of manufactured solutions, Internat.
J. Numer. Methods Fluids 44 (6) (2004) 599–620.
[25] M. Schafer, Computational Engineering. Introduction to Numerical Methods, Springer, 2006.
[26] F. Stern, R. Wilson, J. Shao, Quantitative V-V of CFD simulations and certification of CFD codes, Internat. J. Numer. Methods Fluids 50 (11)
(2006) 1335–1355.
[27] D. Tremblay, S. Etienne, D. Pelletier, Code verification and the method of manufactured solutions for fluid-structure interaction problems, in:
36th AIAA Fluid Dyn. Conf. and Exhib., Vol. 2, 2006, pp. 882–892.
[28] User-Guide, OpenFOAM. The Open Source CFD ToolBox. User Guide. Version 2.2. September 2013. OpenFOAM Foundation, 2013.
[29] S. Veluri, C. Roy, Comprehensive code verification for an unstructured finite volume CFD code, in: 48th AIAA Aerosp. Sci. Meet., 2010.

A Case-Study in Open-Source CFD Code Verification, Part I: Convergence Rate Loss Diagnosis

Uploaded by

Copyright:

Available Formats

A Case-Study in Open-Source CFD Code Verification, Part I: Convergence Rate Loss Diagnosis

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

A Case-Study in Open-Source CFD Code Verification, Part I: Convergence Rate Loss Diagnosis

Uploaded by

Copyright:

Available Formats

Available online at www.sciencedirect.

Mathematics and Computers in Simulation ( ) –

A case-study in open-source CFD code verification,

Keywords: Verification; Manufactured solution; OpenFOAM; CFD; Poisson equation

2. Mathematical models and discretization

A brief introduction to the FVM discretization of Eq. (2) is provided below.

2.1. Numerical schemes

2.1.1. Spatial schemes

(a) Orthogonal mesh. (b) Mesh non-orthogonality.

(c) Mesh skewness. (d) Mesh non-uniformity.

Fig. 1. Mesh properties.

where φC f is evaluated using the interpolation Eq. (6).

The numerical expression φ N|−φ⃗

where a P , a N and R P are constants produced by the discretization process.

[A] [φ] = [R] (17)

2.2. The meshes

(a) Hexahedral orthogonal. (b) Equilateral orthogonal. (c) Stretched orthogonal.

(g) Hexahedral slanted. (h) Hexahedral slanted stretched.

Fig. 2. Orthogonal and non-orthogonal prism meshes in 2D representation.

2.3. Treatment of non-orthogonality

Fig. 3. Non-orthogonality treatment.

Fig. 4. Non-orthogonality treatment: “minimum correction” approach.

Fig. 5. Non-orthogonality treatment: “orthogonal correction” approach.

Fig. 6. Non-orthogonality treatment: “over-relaxed” approach.

3.1. Method of manufactured solution

φ(x, y) = φ0 [cos(ωx) sin(ωy)] (24)

q(x, y) = 2νω2 φ(x, y). (25)

3.2. Error estimation

3.3. Convergence rate

3.4. Boundary conditions (BCs)

3.5. OpenFOAM spatial schemes for the diffusion operator

• Gauss second order, based on Gaussian integration

4.1. Initial test verification

Fig. 7. Curves on orthogonal meshes.

Fig. 8. Curves on orthogonal elongated meshes.

Fig. 9. Curves on non-orthogonal meshes with varying limited parameter ψ.

4.3. Grid refinement results on non-orthogonal meshes

4.3.2. Results on different types of non-orthogonal triangular mesh

4.3.3. Testing non-orthogonality that systematically inclines the border

4.3.4. Non-orthogonal mesh with a periodic (cyclic) BC

Fig. 11. Curves on non-orthogonal meshes.

Fig. 12. Incidence on non-orthogonal meshes inclining an initial square domain.

Fig. 13. Curves using periodic (cyclic) BCs.

4.4. Convergence rate for the gradient

4.5. Implementing a non-orthogonal treatment at the BC

Fig. 14. Curves on non-orthogonal meshes for the gradient of φ.

4.6. Non-orthogonal mesh with Neumann BC

Fig. 15. Curves using an analytical gradient on the non-orthogonal correction.

Fig. 16. Curves using Dirichlet and adiabatic BCs.

Fig. 17. Curves using Dirichlet and Neumann BCs.

You might also like