A Truncation Error Estimation Scheme for the Finite Volume Method on Unstructured Meshes

TL;DR

This paper introduces a new approximate scheme for estimating volume-based truncation errors in finite volume analysis on unstructured meshes, aiming to improve accuracy for mesh adaptation in laminar flow simulations.

Contribution

It develops a novel truncation error estimation method based on neglected terms, offering better accuracy on unstructured meshes compared to previous approaches.

Findings

The scheme provides more accurate truncation error estimates.

It is suitable for use in mesh adaptation processes.

The method is applicable to laminar flow simulations on unstructured meshes.

Abstract

This work is an attempt to develop an approximate scheme for estimating the volume-based truncation errors in the finite volume analysis of laminar flows. The volume-based truncation error is the net flow error across the faces of a control volume. Unfortunately, truncation error is not a natural outcome of the finite volume solution and needs to be estimated separately. Previous works in the literature estimate truncation error using either higher order interpolation schemes, higher order discretization schemes, or neglected terms in the discretization scheme. The first two approaches become complicated on general unstructured meshes and the third approach provides inaccurate results. This work proposes a truncation error estimation scheme, which is based on the third approach, but provides more accurate results compared to the existing results in the literature. The potential application of such a truncation error estimation scheme is in mesh adaptation.

Equations

$$\rho \nabla \cdot \mathbf{v}=0$$

$$\displaystyle\nabla(\rho{\bf v}\otimes{\bf v})=-\nabla p+\nabla\cdot\Big{[}\mu\left(\nabla{\bf v}+\nabla{\bf v}^{\rm T}\right)\Big{]}$$

$${\oint }_{\partial {\Omega }_i}\rho \mathbf{v}\cdot \hat{\mathbf{n}}dA=0$$

$${\oint }_{\partial {\Omega }_i}(\rho \mathbf{v}\otimes \mathbf{v})\cdot \hat{\mathbf{n}}dA=-{\oint }_{\partial {\Omega }_i}p\hat{\mathbf{n}}dA+{\oint }_{\partial {\Omega }_i}\mu (\nabla \mathbf{v}+\nabla {\mathbf{v}}^T)\cdot \hat{\mathbf{n}}dA$$

$$\sum _{\mathrm{faces}}{J}_{\mathrm{mass}}=0$$

$$\sum _{\mathrm{faces}}{\mathbf{J}}_{\mathrm{adv}}=-\sum _{\mathrm{faces}}{\mathbf{F}}_{\mathrm{pres}}+\sum _{\mathrm{faces}}{\mathbf{F}}_{\mathrm{visc}}$$

$${\psi }_f\approx \frac{1}{2}({\psi }_1+{\psi }_2)+\frac{1}{4}(\nabla {\psi }_1+\nabla {\psi }_2)\cdot ({\mathbf{r}}_1+{\mathbf{r}}_2)+\mathcal{O}(h^2)$$

$$\nabla {\psi }_f\approx \frac{1}{2}(\nabla {\psi }_1+\nabla {\psi }_2)+\mathcal{O}(h)$$

$${\frac{\partial \psi }{\partial n}|}_f\approx (\nabla \psi \cdot \hat{\mathbf{n}}{)}_f\approx \alpha \left(\frac{{\psi }_2-{\psi }_1}{|\mathbf{s}|}\right)+\nabla {\psi }_f\cdot (\hat{\mathbf{n}}-\alpha \hat{\mathbf{s}})+\mathcal{O}(h)$$

$$J_{{\mathrm{mass}}_{12}}={\int }_{\mathrm{face}}\rho \mathbf{v}\cdot \hat{\mathbf{n}}dA\approx \rho \left[{\bar{V}}_n-d_f\left({\frac{\partial p}{\partial n}|}_{\mathrm{face}}^{\mathrm{active}}-{\frac{\partial p}{\partial n}|}_{\mathrm{face}}^{\mathrm{lagged}}\right)\right]$$

$$\begin{array}{rl}{J}_{{\mathrm{mass}}_{12}}=& \frac{\rho A_f}{2}({\mathbf{v}}_1+{\mathbf{v}}_2)\cdot \hat{\mathbf{n}}+\frac{\rho \alpha d_f{A}_f}{|\mathbf{s}|}(p_1-p_2)+\\ & \frac{\rho A_f}{4}[\nabla ({\mathbf{v}}_1\cdot \hat{\mathbf{n}})+\nabla ({\mathbf{v}}_2\cdot \hat{\mathbf{n}})]\cdot ({\mathbf{r}}_1+{\mathbf{r}}_2)+\frac{\rho \alpha d_f{A}_f}{2}(\nabla p_1+\nabla p_2)\cdot \hat{\mathbf{s}}\end{array}$$

$${\mathbf{J}}_{{\mathrm{adv}}_{12}}={\int }_{\mathrm{face}}(\rho \mathbf{v}\otimes \mathbf{v})\cdot \hat{\mathbf{n}}dA\approx \rho J_{{\mathrm{mass}}_{12}}({\mathbf{v}}_{\mathrm{up}}+\nabla {\mathbf{v}}_{\mathrm{up}}\cdot {\mathbf{r}}_{\mathrm{up}})$$

$${\mathbf{F}}_{{\mathrm{pres}}_{12}}$$

$${\mathbf{F}}_{{\mathrm{visc}}_{12}}$$

$$\approx \mu \left[\alpha \left(\frac{{\mathbf{v}}_2-{\mathbf{v}}_1}{|\mathbf{s}|}\right)+\frac{1}{2}(\nabla {\mathbf{v}}_1+\nabla {\mathbf{v}}_2)\cdot (\hat{\mathbf{n}}-\alpha \hat{\mathbf{s}})+\frac{1}{2}(\nabla {\mathbf{v}}_1^T+\nabla {\mathbf{v}}_2^T)\cdot \hat{\mathbf{n}}\right]A_f$$

$${\varphi }_i+\sum _{j\in \mathrm{nb}}{c}_{ij}{\varphi }_j=b_i$$

$${\mathcal{L}}_h(\Phi )=0$$

$${\mathcal{L}}_h^2({\varphi }_h^2)=0$$

$${\epsilon }_h^2={\varphi }_h^2-\Phi$$

$${\mathcal{L}}_h({\varphi }_h^2)={\delta }_h^2$$

$${\psi }_f\approx \frac{1}{2}({\psi }_1+{\psi }_2)+\frac{1}{4}(\nabla {\psi }_1+\nabla {\psi }_2)\cdot ({\mathbf{r}}_1+{\mathbf{r}}_2)+$$

$$\frac{1}{4}{\mathbf{r}}_1^T\cdot (\nabla \nabla {\psi }_1+\nabla \nabla {\psi }_2)\cdot {\mathbf{r}}_2+\mathcal{O}(h^3)$$

$$\nabla {\psi }_f\approx \frac{1}{2}(\nabla {\psi }_1+\nabla {\psi }_2)+\frac{1}{4}(\nabla \nabla {\psi }_1+\nabla \nabla {\psi }_2)\cdot ({\mathbf{r}}_1+{\mathbf{r}}_2)+\mathcal{O}(h^2)$$

$${\frac{\partial \psi }{\partial n}|}_f=(\nabla \psi \cdot \hat{\mathbf{n}}{)}_f\approx \alpha \left(\frac{{\psi }_2-{\psi }_1}{|\mathbf{s}|}\right)+\nabla {\psi }_f\cdot (\hat{\mathbf{n}}-\alpha \hat{\mathbf{s}})+$$

$$\frac{1}{4}({\mathbf{r}}_1+{\mathbf{r}}_2{)}^T\cdot (\nabla \nabla {\psi }_1+\nabla \nabla {\psi }_2)\cdot \hat{\mathbf{n}}+\mathcal{O}(h^2)$$

$$\Delta J_{\mathrm{mass}}\approx \frac{\rho A_f}{4}\left[\nabla \nabla ({\mathbf{v}}_1\cdot \hat{\mathbf{n}})+\nabla \nabla ({\mathbf{v}}_2\cdot \hat{\mathbf{n}})\right]:({\mathbf{r}}_1\otimes {\mathbf{r}}_2+\frac{1}{12}\mathbf{t}\otimes \mathbf{t})$$

$$\Delta {\mathbf{J}}_{\mathrm{adv}}\approx$$

$$\frac{\rho A_f}{4}({\mathbf{v}}_1\cdot \hat{\mathbf{n}}+{\mathbf{v}}_2\cdot \hat{\mathbf{n}})\left[\nabla \nabla {\mathbf{v}}_{\mathrm{up}}:\left({\mathbf{r}}_{\mathrm{up}}\otimes {\mathbf{r}}_{\mathrm{up}}+\frac{1}{12}\mathbf{t}\otimes \mathbf{t}\right)\right]+$$

$$\frac{\rho A_f}{24}\left[\nabla ({\mathbf{v}}_1\cdot \hat{\mathbf{n}})+\nabla ({\mathbf{v}}_2\cdot \hat{\mathbf{n}})\cdot \mathbf{t}\right]\cdot (\nabla {\mathbf{v}}_{\mathrm{up}}\cdot \mathbf{t})+\mathcal{O}(h^3)$$

$$\Delta {\mathbf{F}}_{\mathrm{pres}}\approx$$

$$\Delta {\mathbf{F}}_{\mathrm{visc}}\approx$$

$$\Delta {\mathbf{J}}_{\mathrm{mom}}=\Delta {\mathbf{J}}_{\mathrm{adv}}+\Delta {\mathbf{F}}_{\mathrm{pres}}-\Delta {\mathbf{F}}_{\mathrm{visc}}$$

Full text

A Truncation Error Estimation Scheme for the Finite Volume Method on Unstructured Meshes

A. R. Baserinia

Plethora Corporation, San Francisco, CA 94124, United States

1 Abstract

This work is an attempt to develop an approximate scheme for estimating the volume-based truncation errors in the finite volume analysis of laminar flows. The volume-based truncation error is the net flow error across the faces of a control volume. Unfortunately, truncation error is not a natural outcome of the finite volume solution and needs to be estimated separately. Previous works in the literature estimate truncation error using either higher order interpolation schemes, higher order discretization schemes, or neglected terms in the discretization scheme. The first two approaches become complicated on general unstructured meshes and the third approach provides inaccurate results. This work proposes a truncation error estimation scheme, which is based on the third approach, but provides more accurate results compared to the existing results in the literature. The potential application of such a truncation error estimation scheme is in mesh adaptation.

2 Introduction

In numerical simulation of fluid flows, solution error is inevitable. This error is primarily due to the discretization of the governing equations. In other words, the discrete solution obtained from a numerical simulation is slightly different from the exact solution of the governing equations. Therefore substituting the numerical solution into the exact governing equations results in an apparent source term, called the truncation error.

The main application of solution truncation error is in CFD problems in which the assessment and reduction of the solution error is crucial. In general, there are two approaches that adopt the concept of truncation error to reduce the discretization error. The first approach is the use of the solution truncation error as a source term in the governing equations to construct a higher order discretization scheme [4, 6, 11, 17]. The second approach is the use of the solution truncation error as an indicator for driving a mesh adaptation process [13, 15, 16, 18]. This work concentrates on a truncation error estimation scheme that could potentially be used for the second purpose.

In the application of the solution truncation error as an indicator for mesh adaptation, we are interested in having a reasonable estimate of the distribution of the solution truncation error. However, this estimate does not need to be highly accurate. The reason is that the truncation error estimate will eventually be used for determining the characteristic sizes of mesh elements. In practice, we do not have complete control over these characteristic size. Therefore even an approximate truncation error estimator that correctly predicts the overall trend of truncation error distribution is adequate. The objective of this paper is to present a simple and straightforward truncation error estimator that can be used in mesh adaptation applications.

3 Conservation Equations and FVM Discretization

In this work, we focus on the conservation equations for two-dimensional steady-state laminar incompressible isothermal flows of constant-property Newtonian fluids in the absence of gravity, which are governed by the mass and momentum equations:

[TABLE]

where ${\bf v}=(u,v)$ is the flow velocity vector, $p$ is the pressure, $\rho$ is the density, $\mu$ is the viscosity, and the superscript $[\;]^{\rm T}$ is the transpose operator. Note that the mass equation (1) is scalar while the momentum equation (2) is vectorial. Therefore the entire system consists of three equations.

The first step in a finite volume solution is to discretize the governing equations on a mesh. In this work, we use the cell-centered finite volume method on unstructured meshes with arbitrarily shaped cells, shown in Figure 1. The discretization starts with integrating Equations (1) and (2) over control volumes to obtain the integral form of the governing equations. For this purpose we use a cell-centered finite volume approach, shown in Figure 1. Then we use the Gauss theorem to transform the volume integrals into surface integrals.

[TABLE]

where ${\Omega_{i}}$ is the $i$-th control volume and $\partial{\Omega_{i}}$ is its boundary, or rather its faces.

The next step in the discretization process is to approximate the surface integrals of Equations (5) and (6) in terms of the average flow across the control volumes faces.

[TABLE]

where $J_{\rm mass}$ is the mass flow rate, ${\bf J}_{\rm mom}$ is the vector of momentum flow rate, and ${\bf F}_{\rm pres}$ and ${\bf F}_{\rm visc}$ are the vectors of pressure and viscous forces, respectively. Note that the values of these flows and forces at each face are unknown and need to be determined.

To determine the values of mass and momentum flows across a face and those of forces acting on it, we need to interpolate the pressure and velocities from the nodes of the neighboring control volumes. For example if we use a second-order central discretization scheme, the values of a typical variable, $\psi$, and its gradient, $\nabla\psi$, at the interface of control volumes 1 and 2, shown in Figure 2, become [3]:

[TABLE]

where $h$ is the local mesh characteristic size. Following Ferziger and Peric [7] and Zwart et al. [19], we may discretize the face-normal component of the gradient vector differently to improve the numerical stability of the scheme:

[TABLE]

where $\alpha=\hat{\bf n}\cdot\hat{\bf s}$ represents the mesh nonorthogonality [19]. Note that the approximations in Equations (8) and (9) are formally first-order accurate on a nonuniform mesh. However both approximations exhibit a second-order convergence as the mesh is refined [7]. Therefore the overall order of accuracy of the discretization method would still be two.

Using Equations (7–9), we can approximate the values of velocity and pressure at an interior face. For a boundary face we may use a similar approach but the central discretization scheme must be replaced by a one-sided extrapolation scheme [3]. Once the velocity and pressure values are determined at a face, we can calculate the mass and momentum flows and pressure and viscous forces at the face.

To calculate the mass flow rate across a face, we use the Rhie and Chow velocity pressure interpolation scheme [14]. This scheme prevents the infamous velocity-pressure decoupling from happening. The result is:

[TABLE]

where $\bar{V}_{n}$ is the advected velocity and $d_{f}$ is the dissipation coefficient [12, 19]. Using Equations (7), (8), and (9) for calculating $\bar{V}_{n}$, ${\partial p}/{\partial n}|^{\rm lagged}$, and ${\partial p}/{\partial n}|^{\rm active}$, respectively, results in the discretized form of the mass flow across an interior face:

[TABLE]

where $A_{f}$ is the face area in Fiqure 2. To calculate the momentum flow vector across a face, ${\bf J}_{\rm adv}$, we need to use the upwind approximation:

[TABLE]

where ${\bf v}_{\rm up}={\bf v}_{1}$ if $J_{{\rm mass}_{12}}>0$ and ${\bf v}_{\rm up}={\bf v}_{2}$ otherwise. We can also calculate the pressure and viscous forces acting on a face using Equations (7–9):

[TABLE]

Substituting Equations (11–14) into Equations (5) and (6) results in the discretized form of the mass and momentum equations. Although the final result is tedious, we can represent it in the following symbolic form:

[TABLE]

where $\phi_{i}=(p_{i},u_{i},v_{i})^{T}$ is the vector of primitive variables at the $i$-th control volume and $j$ is the index of its neighboring control volumes, denoted by ‘nb’ in Figure 1. Note that in Equation (15), the coefficients $c_{ij}$ and $b_{i}$ are generally functions of $\phi$. Therefore it is a nonlinear equation which should be solved using iterative methods.

4 Truncation Error in the Finite Volume Method

Let us describe the concept of truncation error in the finite volume context. Suppose we represent the integral form of the governing equations, Equations (3) and (4), in the following symbolic form:

[TABLE]

where ${\mathscr{L}}_{h}$ is the integral conservation operator on a mesh with the characterisitic size $h$, acting on the exact solution vector, $\Phi=(p,u,v)^{\rm T}$. Nevertheless in the finite volume method we solve the discretized form of the above equation. In the case of a second-order discretization, which is the most common choice in CFD applications, the symbolic form of the discretized equation is:

[TABLE]

where ${\cal L}_{h}^{2}$ is the second-order discrete conservation operator on a mesh with the characteristic size $h$ and $\phi_{h}^{2}$ is numerical solution. The difference between the numerical solution, $\phi_{h}^{2}$, and the exact solution, $\Phi$, is the discretization error.

[TABLE]

Note that the numerical solution, $\phi_{h}^{2}$, may not satisfy the exact conservation operator, ${\mathscr{L}}_{h}$. Therefore substituting the numerical solution, $\phi_{h}^{2}$, into Equation (16) results in an apparent source term, called the truncation error and shown by $\delta_{h}^{2}$:

[TABLE]

The latter shows that the truncation error of a numerical solution is the apparent source term required in the exact equations to satisfy the numerical solution. On the other hand, using the divergence theorem on a control volume, we can prove that the apparent source term is equal to the net flow imbalance across the boundaries of the control volume. Therefore the truncation error is simply the net flow error across the boundaries of a control volume. This definition of truncation error is on a volume-based basis, contrary to the face-based definitions in Roe and Nishikawa [15] and Reuss and Stubley [13].

5 Calculating the Volume-Based Truncation Error

In the finite volume context, the truncation error is equal to the net imbalance of mass and momentum flows across the faces of each control volume as mentioned in the previous section. Although this definition of truncation error is based on Equation (19), we cannot simply use it to calculate the volume-based truncation error of the numerical solution, $\phi_{h}^{2}$. Note that $\phi_{h}^{2}$ is only peace-wise continuous within control volumes but discontinuous across control volume faces. Therefore it is not feasible to evaluate ${\mathscr{L}}_{h}(\phi_{h}^{2})$ since it involves evaluating integrals on the faces.

The alternative approach to truncation error calculation is to replace the exact operator ${\mathscr{L}}_{h}$ by an accurate discrete operator. Due to the consistency condition [5], the most obvious choice for an accurate discrete operator is ${\cal L}_{h^{\prime}}^{2}$ where $h^{\prime}\ll h$; in other words the second-order discrete operator on a very fine mesh. To understand how to estimate the truncation error using this operator, imagine each face of a two-dimensional mesh is divided into $n$ smaller faces. Therefore each control volume is split into $n^{2}$ smaller control volumes, shown in Figures 3(b) and 3(c). If we interpolate the second-order solution $\phi_{h}^{2}$ from the original mesh of Figure 3(b) to the fine mesh of Figure 3(c) and apply the conservation operator ${\cal L}_{h^{\prime}}^{2}$ on it, we can obtain the truncation errors of individual control volumes on the fine mesh. To obtain the truncation error on the original mesh of Figure 3(b), it suffices to find the summation of truncation errors of the smaller control volumes in Figure 3(c). The reason is that the face flow errors of the smaller control volumes cancel out each other except at the faces of the original control volume of Figure 3(b), which provides the truncation error. Unfortunately this approach is both hard to implement and intense on computer memory requirement. Therefore it is rarely used in practical applications.

In the literature, there are three techniques for estimating the face flow errors [9]:

• •

Estimating the neglected terms in the discretization scheme;

• •

Recovering a higher order accurate solution based on the discrete solution; and

• •

Use of a higher order accurate discretization scheme.

The first technique for estimating the face flow errors is based on analyzing the neglected terms in the Taylor series expansion in the discretization scheme. The neglected terms are estimated and used for establishing a higher order accurate face flow calculation. Therefore these neglected terms provide an estimate for the mass and momentum face flow errors. Examples of this technique in the literature are the works of Ilinca et al. [10], Reuss and Stubley [13], and Zhang et al. [18].

The second technique for estimating the face flow errors is based on recovering a higher order accurate solution from the available discrete solution in order to obtain a more accurate face flow estimation. In this technique, the discrete solution on an initial mesh is used along with a higher order accurate interpolation scheme. Since the discrete solution is only given at the control volume nodes, the higher order interpolation provides a more accurate estimate of the solution variables at the faces. Therefore an estimate for the face flow error can be obtained. Examples of this method can be seen in recent works of Hay and Visonneau [9] and Hay et al. [8]

The third technique for estimating the face flow errors is based on using a higher order discretization scheme. In this technique, one uses the numerical solution obtained from an $n$-th order accurate discretization scheme along with another discretization scheme which is $m$-th order accurate where $m>n$. Substitution of the $n$-th order solution into the $m$-th order discretization scheme results in an estimate of the truncation error. This technique, which is referred to as the defect correction method, is elaborated on by Ervin and Layton [6] and Pierce and Giles [11].

Note that all the above techniques are conceptually equivalent. Figure 3(d) shows the basic concept of these techniques. As seen, these techniques try to estimate the face flow errors of Figure 3(b) by using higher order terms in the discretization scheme. The advantage of the first two techniques is the ease of implementation while the third technique is more cumbersome to develop and implement. In addition, the first technique does not involve complicated interpolations, which is an advantage in terms of implementation especially on unstructured meshes. Therefore we use the first technique to estimate the face flow errors.

6 The Proposed Truncation Error Estimation Scheme

This section explains a truncation error estimation scheme that is based on evaluating the first neglected Taylor series terms in the calculation of the mass and momentum flow errors. A second-order accurate discretization scheme uses Equations (7–9) to calculate the face flows. Therefore one needs to account for the first neglected terms in these equations:

[TABLE]

[TABLE]

In these equations $\nabla\nabla\psi$ is the tensor of second-order derivatives, called the Hessian tensor. Using the equations above and following the procedure of Section 3, one can derive new equations for the mass and momentum flows across a face. These relations would be similar to Equations (10–14) with some extra terms, which account for the face flow errors. The final result for the mass face flow error across an interior face, shown in Figure 2, is:

[TABLE]

where ${\bf t}=A_{f}\hat{\bf t}$ in Figure 2, and the symbols $(\otimes)$ and $(:)$ represent the tensor product of two vectors and the Frobenius product of two tensors, respectively. Similarly one can obtain the face flow and force errors in the momentum equation:

[TABLE]

And the total face error for the momentum flow across a face based on the equation above and Equation (6) becomes:

[TABLE]

Note that Equations (23) and (27) only provide the mass and momentum flow errors across a face. To calculate the total truncation error for a control volume, one has to find the summation of flow errors across all the control volume faces.

Equations (23) and (27) provide a straightforward scheme for estimating mass and momentum truncation errors in the context of finite volume methods on general unstructured meshes. The only problem with these equations is the requirement to calculate of the Hessian tensor of the primitive variables, which do not appear in a second-order discretization. However, this is not a major issue since one can estimate the solution Hessians by applying a gradient reconstruction procedure to the already known solution gradients [1, 13]. In the following section, we apply the proposed truncation error scheme to two academic test cases and discuss it performance.

7 Results and Discussion

To show the performance of the truncation error estimation scheme, discussed in the previous section, we apply it to two academic test cases: incompressible laminar flows in a lid-driven square cavity at ${\rm Re}=1000$ and over a backward facing step at ${\rm Re}=400$. Then we discuss the strengths and shortcomings of the proposed scheme.

7.1 Flow in a Lid-Driven Cavity

The lid-driven cavity flow is an important test case since it exhibits many important phenomena of fluid flows in a simple geometry. Figure 4 shows the schematic of the lid-driven square cavity flow. The Reynolds number of the flow $\mbox{Re}=UL/\nu$, which is defined based on the lid velocity and the cavity size, is 1000. The boundary conditions of the flow are no-slip walls. To set the pressure level, the pressure at the lower left corner is set to zero.

To obtain the basic second-order solution, we solve the cavity flow on a uniform unstructured triangular mesh with 3602 control volumes. To validate the truncation error estimates we have to find the actual volume-based truncation error with a reasonable accuracy. For this purpose we interpolate the basic solution to a refined mesh and then apply the second-order conservation operator on the refined mesh on the basic solution, discussed in Section 3. Nevertheless we need to confirm that the truncation error obtained using this method is accurate enough. For this purpose we use a mesh refinement study by refining the mesh until the truncation error estimates converge. Based on this technique, if we subdivide each face of the original mesh into 16 smaller faces, the RMS error in the truncation error calculation would be smaller than one percent. Therefore we use the second-order conservation operator on a mesh with 922112 control volumes as an approximation to the exact conservation operator.

To obtain the truncation error estimates based on the proposed scheme we use Equations (23) and (27). Figures 5(a) through 5(f) show the distribution of the volume-based truncation error scaled by the average truncation error across the domain. Note that the contour levels in these figure are logarithmic and as a result these figures represent the truncation error order of magnitude.

A brief examination of Figure 5 shows that the distributions of the estimated truncation errors are very similar to those of the actual truncation errors. In all cases, the truncation errors near the driven lid and the right wall are orders of magnitude larger than those of the rest of the domain. This is due to the strong velocity and pressure gradients experienced in the near wall regions. In contrast, in the lower left and right corners where two low energy separation bubbles form, the truncation errors are small. These observations show that the proposed truncation error estimation scheme successfully predicts the overall trend of the truncation error distribution.

7.2 Flow over a Backward Facing Step

The second test case that we use to examine the performance of the proposed truncation error estimation scheme is the laminar flow over a confined backward step [2], shown in Figure 6. In this work, the Reynolds number based on the total height of the domain, $2H$, and the average inlet velocity, $U_{\rm ref}$ is 400 and $L=10H$. The inlet velocity profile is profile is parabolic and the outlet pressure is set to zero.

The main difference between the backward facing step flow and the cavity flow of the previous section is the role of the pressure gradient. In the backward facing step flow, the pressure gradient is mild everywhere in the domain except near the sharp corner of the step. Therefore it is easier to examine the effect of the viscous term on the performance of the proposed truncation error estimation method.

To obtain the basic second-order solution, we use a uniform unstructured triangular mesh with 5120 control volumes. To obtain the solution truncation error on the original mesh, we subdivide each face into 16 smaller faces. A mesh refinement analysis confirm that the truncation error calculated based on conservation operator on this mesh is adequately accurate with an error level smaller than one percent.

Figure 7 shows the final results for the volume-based truncation error of the mass and momentum equations. As seen, the overall trends of the estimated truncation errors are similar to those of actual truncation errors. However, there are certain discrepancies between the estimated and the actual truncation error. For example the results for the mass truncation error, Figures 7(a) and 7(b) show that the proposed scheme overestimates the truncation error in the separation bubble behind the step. In contrast, the proposed scheme underestimates the truncation error near the outlet boundary in the vicinity of walls. Similar discrepancies can also be found in the $x$ and $y$ momentum truncation error. The main reason for this discrepancies is the way that we estimated the viscous force error across a face.

Equation (26) shows the first neglected term in the discretization of the face viscous force, which is first order. A brief examination of this term reveals that on a uniform mesh, the vector ${\bf r}_{1}+{\bf r}_{2}=0$ (see Figure 2). Therefore the proposed truncation error estimation scheme fails to resolve face viscous force errors on uniform meshes. To fix this problem, higher order terms must be taken into account. Unfortunately this leads to the appearance of third order derivatives of the solution variables, which are not easy to estimate for a second-orders solution.

8 Conclusion

This work proposed a simple truncation error estimation scheme in the context of the finite volume method. In this context, the solution truncation error in each control volume is the sum of flow errors across the faces of the control volume. To calculate the face flow errors, the first neglected terms in the discretization of the mass and momentum flows across faces were used. Then these terms were assembled to obtain the solution truncation error. The numerical experiments with the lid-driven cavity and backward facing step flow showed that the proposed scheme successfully estimated the overall trend of truncation error distribution. However in the backward facing step flow, certain discrepancies were observed between the estimated and the actual truncation error due to the lack of accurate in resolving viscous features in the flow. The main advantage of the proposed scheme was its simplicity in the sense that only second-order derivatives need to be reconstructed. However in the cases that the viscous term in the momentum equation is sole source of error, the proposed truncation error estimation scheme may become less accurate.