# Fluid Flow Simulations Based on an Equation-Solving Solution Gradient Strategy

Authors: Ho-Shuenn Huang; Yao-Hsin Hwang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A compact and accurate discretization for fluid flow simulations is introduced in this paper. Contrary to the common wisdom in a convectional scheme, the solution gradient required for a high-resolution scheme is provided by solving its corresponding difference equation rather than by interpolation from solution values at neighboring computational nodes. To achieve this goal, a supplementary equation and its associated control volume are proposed to retain a compact and accurate discretization. Scheme essentials are exposed by numerical analyses on simple one-dimensional modeled problems to reveal its formal accuracy. Several test problems are solved to illustrate the feasibility of present formulation. From the obtained numerical results, it is evident that the proposed scheme will be a useful tool to simulate fluid flow problems in arbitrary domains.


---

## Body

## 1. Introduction

To provide a feasible simulation tool for incompressible flow in a complicated domain, we have successfully developed associated solution procedures on arbitrary Lagrangian-Eulerian (ALE) [1, 2] and subsequently staggered triangular grid systems [3, 4]. Nevertheless, its formal accuracy is quite limited and can be proved to be only first-order with the truncation error analyses [5]. Since a high-resolution scheme is demanded in practical flow simulations, it will be beneficial to raise its discretization accuracy and simulation efficiency.In a finite volume method to simulate the convection-diffusion equation, solution gradient in the computational cell must be appropriately defined to yield a high-order accurate representation of convection term [6, 7]. Meanwhile, solution gradient must also be defined at the computational cell boundary to realize the diffusion term [8]. These solution gradients are generally obtained with the Taylor-series expansion based on the solution values at neighboring computational nodes, which will be an easy task on a regular grid. However, difficulties may be arisen if it is applied on an unstructured mesh to cope with a complicated domain. This process may bring on some computational disadvantages in the solution procedure.(1)
The resulting expression for solution gradient is indifferent to the peculiarity of governing equation.(2)
It will result in an extended difference stencil.(3)
Its accuracy heavily depends upon grid arrangement.To circumvent these inconveniences, we will employ an alternative approach to determine the solution gradient: it should be directly calculated by solving its corresponding governing equation just as the solution variable. In this way, the characteristics pertaining to the governing equation can be directly deliberated in the solution gradient and consequently increase the simulation accuracy. To achieve this goal, the present paper is designed to unveil the underlying basics to be used in fluid flow problems. It can be summarized as follows:(1)
introducing a supplementary equation for solution gradient,(2)
constructing a control volume for the supplementary equation,(3)
establishing a practical expression for diffusion flux at cell boundaries,(4)
conducting numerical analyses on simple modeled problems,(5)
detailing the solution procedure in fluid flow calculations,(6)
verifying the solution accuracy by solving test problems,It is noted that only the essential ingredients are delineated in the present paper. Detailed descriptions on associated formulation fundamentals, numerical analyses, solution procedure and more extensive validations can be consulted in our previous work [9–13].


## 2. Formulation

### 2.1. Governing Equations

Governing equations for incompressible flow can be generalized by the following convection-diffusion-source form:(1)∇·f^=S,
with flux function f^=ρvΦ-μ∇Φ. To make the differential equation (1) numerically tractable, the interested domain is divided into nonoverlapping polygons in the finite volume formulation. Figure 1(a) shows a typical polygon delegated by a set of irregularly distributed grid nodes. The computational node for solution variable is assigned as the centroid of its associated polygon. Meanwhile, solution variable is mimicked by a combination of piecewise linear distributions:(2a)Φ(r)≈∑P=1nCϕ(r;P),
where ϕ(r;P) is the locally supported approximate linear function in association with Pth computational cell,
(2b)ϕ(r;P)={ϕP+∇ϕP·(r-rP)ifr∈VP0otherwise,Grid arrangement.
(a)
Control volumes
(b)
Local coordinate
(c)
Integration regionwith∇ϕP=pPi+qPj. Physically, parameters ϕP, pP, and qP can be related to the solution variable and its spatial derivatives at the computational node.In the present proposition, these spatial derivatives are directly computed based on their respective governing equations. An additional equation set derived from the original equation (1) is employed as the transporting equations for solution gradients:
(3)∇·(∇f^)=∇S.
These are the supplementary equations whose discretizations are also of vital importance in our formulation.


### 2.2. Control Volumes and Flux-Balanced Forms

In the finite volume formulation, difference equation is derived by integrating the governing equations over their corresponding control volumes. The control volume for original equation (1) is designated as the computational cell, CVo,P=VP, which will ensure global conservation for solution variable. Nevertheless, the control volume for the supplementary equation (3), CVs,P, is chosen as a subset of VP. These control volumes are displayed in Figure 1(a), where corner nodes of CVs,P are located at the midpoints of computational and grid nodes. Therefore, we denote CVs,P as the supplementary control volume. With standard procedure in finite volume methodology, the difference counterpart for original equation (1) can be expressed by the following flux balanced form:
(4)∑m=1nVfm·Am=s-PVP,
where fm·Am resembles the volume flux across the mth boundary surface of  CVo,P and s-PVP is the resulting integration of source term over the control volume. Similarly, the flux balanced form for supplementary equation (3) reads
(5)∑m′=1nV(∇f)m′·Am′=∑m′=1nVsm′·Am′,
where m′ denotes the boundary face of the supplementary control volume, CVs,P.To express the flux derivative term,(∇f)m′·Am′, in (5) with flux functions, we introduce a local coordinate (ξ,η) defined by the unit base vectors:
(6)eξ=rm-rP|rm-rP|,eη=rk-rkp|rk-rkp|.
These unit vectors are illustrated in Figure 1(b) for reference. By manipulating the coordinate transformation between the local and global systems, (∇f)m′ can be rearranged with the following expression:
(7)(∇f)m′=fm-fPΔξmP∇ξ+12(∇fP+∇fm)·∂r∂η∇η.
Substituting these terms into (5) will yield
(8)∑m=1nV1Δξmp[(fm·Am)-(fP·Am)]∇ξ+12{[(∇f)P·Am+(∇f)m·Am]·∂r∂η}∇η=∑m=1nVsm′Am.
In this way, (4) and (8) form a system of difference equations to imitate the original and supplementary equations (1) and (3), respectively.


### 2.3. Numerical Flux and Difference Equations

The numerical flux function at the computational node,fP·Am, can be realized without difficulty, because the approximate solution (2a), (2b) is continuous at this point:
(9)fP·Am=(ρv)P·AmϕP-μP∇ϕP·Am.
However, at the computational cell boundary, it should be handled with care since the approximate solution (2a), (2b)  becomes discontinuous. First, the upwind treatment is applied to derive the convection part:
(10a)(ρvϕ·A)m=(ρv·A)m+ϕ(rm;P)+(ρv·A)m-ϕ(rm;C).Solution gradient in the diffusion term is resorted to the divergence theorem:(11)∇Φm≈∇ϕm=1CVd∫CVd∇ϕdV=1CVd∫A(CVd)ϕdA.
The integration volume, CVd, consists of a similar region as that enclosed by the lines connecting the grid and computational nodes surrounding the boundary surface, which is illustrated in Figure 1(c). This integration volume can be divided into two parts pertaining to neighboring computational cells:
(12)CVd=VP,m+VC,m.
The resulting solution gradient reads
(10b)∇ϕm=VP,mVP,m+VC,m∇ϕP+VC,mVP,m+VC,m∇ϕC+ωAmVP,m+VC,m[ϕ(rm;C)-ϕ(rm;P)].
In the previous equation, scheme parameter ω indicates the integration factor to yield  CVd. Further interpolations of this scheme factor will be given in a subsequent section.Collecting the discretized convection and diffusion terms will lead to the resulting numerical flux term at cell boundary:(13)fm·Am=[(ρv·A)m++μmω|Am|2VP,m+VC,m]ϕP+[(ρv·A)m--μmω|Am|2VP,m+VC,m]ϕC+{[(ρv·A)m++μmω|Am|2VP,m+VC,m](rm-rP)-μmωVP,mVP,m+VC,mAm}·∇ϕP+{[(ρv·A)m--μmω|Am|2VP,m+VC,m](rm-rC)-μmωVC,mVP,m+VC,mAm}·∇ϕC.
Only variables at two adjacent cells are employed to determine the numerical flux function. With flux functions (9) and (13), their derivatives can then be obtained:
(14)∇fP·Am=[(ρv)P·Am]∇ϕP,∇fm·Am=[(ρv·A)m++μmω|Am|2VP,m+VC,m]∇ϕP+[(ρv·A)m--μmω|Am|2VP,m+VC,m]∇ϕC.
By substituting these numerical functions into flux balanced equations (4) and (8), one can obtain the resulting difference equations for the original and supplementary governing equations, which can be expressed with the following form:
(15)aPϕP+bPpP+cPqP+∑m=1nVaCϕC+bCpC+cCqC=d.


## 2.1. Governing Equations

Governing equations for incompressible flow can be generalized by the following convection-diffusion-source form:(1)∇·f^=S,
with flux function f^=ρvΦ-μ∇Φ. To make the differential equation (1) numerically tractable, the interested domain is divided into nonoverlapping polygons in the finite volume formulation. Figure 1(a) shows a typical polygon delegated by a set of irregularly distributed grid nodes. The computational node for solution variable is assigned as the centroid of its associated polygon. Meanwhile, solution variable is mimicked by a combination of piecewise linear distributions:(2a)Φ(r)≈∑P=1nCϕ(r;P),
where ϕ(r;P) is the locally supported approximate linear function in association with Pth computational cell,
(2b)ϕ(r;P)={ϕP+∇ϕP·(r-rP)ifr∈VP0otherwise,Grid arrangement.
(a)
Control volumes
(b)
Local coordinate
(c)
Integration regionwith∇ϕP=pPi+qPj. Physically, parameters ϕP, pP, and qP can be related to the solution variable and its spatial derivatives at the computational node.In the present proposition, these spatial derivatives are directly computed based on their respective governing equations. An additional equation set derived from the original equation (1) is employed as the transporting equations for solution gradients:
(3)∇·(∇f^)=∇S.
These are the supplementary equations whose discretizations are also of vital importance in our formulation.


## 2.2. Control Volumes and Flux-Balanced Forms

In the finite volume formulation, difference equation is derived by integrating the governing equations over their corresponding control volumes. The control volume for original equation (1) is designated as the computational cell, CVo,P=VP, which will ensure global conservation for solution variable. Nevertheless, the control volume for the supplementary equation (3), CVs,P, is chosen as a subset of VP. These control volumes are displayed in Figure 1(a), where corner nodes of CVs,P are located at the midpoints of computational and grid nodes. Therefore, we denote CVs,P as the supplementary control volume. With standard procedure in finite volume methodology, the difference counterpart for original equation (1) can be expressed by the following flux balanced form:
(4)∑m=1nVfm·Am=s-PVP,
where fm·Am resembles the volume flux across the mth boundary surface of  CVo,P and s-PVP is the resulting integration of source term over the control volume. Similarly, the flux balanced form for supplementary equation (3) reads
(5)∑m′=1nV(∇f)m′·Am′=∑m′=1nVsm′·Am′,
where m′ denotes the boundary face of the supplementary control volume, CVs,P.To express the flux derivative term,(∇f)m′·Am′, in (5) with flux functions, we introduce a local coordinate (ξ,η) defined by the unit base vectors:
(6)eξ=rm-rP|rm-rP|,eη=rk-rkp|rk-rkp|.
These unit vectors are illustrated in Figure 1(b) for reference. By manipulating the coordinate transformation between the local and global systems, (∇f)m′ can be rearranged with the following expression:
(7)(∇f)m′=fm-fPΔξmP∇ξ+12(∇fP+∇fm)·∂r∂η∇η.
Substituting these terms into (5) will yield
(8)∑m=1nV1Δξmp[(fm·Am)-(fP·Am)]∇ξ+12{[(∇f)P·Am+(∇f)m·Am]·∂r∂η}∇η=∑m=1nVsm′Am.
In this way, (4) and (8) form a system of difference equations to imitate the original and supplementary equations (1) and (3), respectively.


## 2.3. Numerical Flux and Difference Equations

The numerical flux function at the computational node,fP·Am, can be realized without difficulty, because the approximate solution (2a), (2b) is continuous at this point:
(9)fP·Am=(ρv)P·AmϕP-μP∇ϕP·Am.
However, at the computational cell boundary, it should be handled with care since the approximate solution (2a), (2b)  becomes discontinuous. First, the upwind treatment is applied to derive the convection part:
(10a)(ρvϕ·A)m=(ρv·A)m+ϕ(rm;P)+(ρv·A)m-ϕ(rm;C).Solution gradient in the diffusion term is resorted to the divergence theorem:(11)∇Φm≈∇ϕm=1CVd∫CVd∇ϕdV=1CVd∫A(CVd)ϕdA.
The integration volume, CVd, consists of a similar region as that enclosed by the lines connecting the grid and computational nodes surrounding the boundary surface, which is illustrated in Figure 1(c). This integration volume can be divided into two parts pertaining to neighboring computational cells:
(12)CVd=VP,m+VC,m.
The resulting solution gradient reads
(10b)∇ϕm=VP,mVP,m+VC,m∇ϕP+VC,mVP,m+VC,m∇ϕC+ωAmVP,m+VC,m[ϕ(rm;C)-ϕ(rm;P)].
In the previous equation, scheme parameter ω indicates the integration factor to yield  CVd. Further interpolations of this scheme factor will be given in a subsequent section.Collecting the discretized convection and diffusion terms will lead to the resulting numerical flux term at cell boundary:(13)fm·Am=[(ρv·A)m++μmω|Am|2VP,m+VC,m]ϕP+[(ρv·A)m--μmω|Am|2VP,m+VC,m]ϕC+{[(ρv·A)m++μmω|Am|2VP,m+VC,m](rm-rP)-μmωVP,mVP,m+VC,mAm}·∇ϕP+{[(ρv·A)m--μmω|Am|2VP,m+VC,m](rm-rC)-μmωVC,mVP,m+VC,mAm}·∇ϕC.
Only variables at two adjacent cells are employed to determine the numerical flux function. With flux functions (9) and (13), their derivatives can then be obtained:
(14)∇fP·Am=[(ρv)P·Am]∇ϕP,∇fm·Am=[(ρv·A)m++μmω|Am|2VP,m+VC,m]∇ϕP+[(ρv·A)m--μmω|Am|2VP,m+VC,m]∇ϕC.
By substituting these numerical functions into flux balanced equations (4) and (8), one can obtain the resulting difference equations for the original and supplementary governing equations, which can be expressed with the following form:
(15)aPϕP+bPpP+cPqP+∑m=1nVaCϕC+bCpC+cCqC=d.


## 3. Numerical Analysis

From the previous derivations, the essential characteristics in our proposition can be summarized as follows:(1)
adoption of a supplementary equation to derive difference equation for solution gradient,(2)
selection of a particular control volume for the supplementary equation,(3)
usage of a variable integration range to obtain solution gradient at cell boundary.With these special features, numerical analyses should be conducted to ensure its feasibility before it can be applied in more complicated circumstances.


### 3.1. Interpretation of Solution Gradient at Cell Boundary

In the one-dimensional linear problem, solution gradient at cell boundary(10b) can be simplified as (16a)ϕx,i+1/2=hihi+hi+1pi+hi+1hi+hi+1pi+1+2ωhi+hi+1[ϕ(xi+1/2;i+1)-ϕ(xi+1/2;i)].
It can be interpreted as the volume averaged between two adjacent cells added by a penalty term arising from solution discontinuity. Therefore, besides the integration factor to yield integration volume, the scheme parameter ω can be regarded as a penalty factor in this expression. Meanwhile, it can also be rewritten as
(16b)ϕx,i+1/2=(1-ω)(hihi+hi+1pi+hi+1hi+hi+1pi+1)+ω(ϕi+1-ϕi)hi/2+hi+1/2.
This is a weighting average between those interpolated from solution variable and its gradient at computational nodes and ω stands for the weighting factor. Furthermore, this solution gradient can be further rearranged as:
(16c)ϕx,i+1/2=ϕ(xi+1/2+hi+1/2ω;i+1)-ϕ(xi+1/2-hi/2ω;i)hi/2ω+hi+1/2ω.This is a simple two-point representation andω can then be designated as distance factor.


### 3.2. Increasing Factor

Consider the case with constant positive convective velocity and viscosity in a uniform grid system (hi=h). Difference equations can then be reduced as
(17)BΘi-1+AΘi+CΘi+1=di,
with the solution vector, Θ=(ϕ,ph/2)T, and coefficient matrixes,
(18)B=(-Rh-ω-Rh+(1-ω)00),A=(Rh+2ωRhωRh+1+ω),C=(-ω-(1-ω)-ω-(1-ω)),di=(h2μs-ih24μ(si+1/4+2s-i-si-1/4)),
where Rh=ρuh/μ is the cell Reynolds number. In the case of ω=1, the solution gradient can be explicitly derived as
(19)pi,ω=1=2Rh+2[ϕi+1-ϕih+h4μ(si+1/4+2s-i-si-1/4)].
This expression clearly depicts that the solution gradient will adapt itself to flow condition (cell Reynolds number and source term). The resulting characteristic (nontrivial) increasing factor of our proposed scheme can be derived as [14](20a)r=Rh2+Rh+ω(2+Rh)Rh+ω(2-Rh).
Obviously, the increasing factor will be dependent on the scheme factor, ω. Meanwhile, the critical cell Reynolds number, Rh,cr, for occurrence of possible oscillatory solution (r=0) is
(20b)Rh,cr=2ωω-1.From the viewpoint of increasing factor, it is interesting to note that UD and CD schemes can be recovered by the proposed scheme withω=0 and ω=∞, respectively. That is, if ω is assigned with a small value, the scheme may become too diffusive as the UD scheme, whereas, if it is assigned with a large value, the resulting solution may be prone to being contaminated by oscillatory distribution as the CD scheme.To determine the possible range of the scheme factor, we turn to the exact increasing factor,rex=eRh,  which determines the exact scheme factor
(21)ωex(Rh)=eRhRh-Rh2-RheRh(Rh-2)+Rh+2.
This is a monotone distribution with ωex(0)=3 and ωex(∞)=1. Therefore, distribution range of this exact scheme factor is
(22)1≤ωex≤3.
This constraint is selected as the bound for scheme factor in our study. Figure 2 illustrates the distribution of increasing factors for ω=1, 2 and 3 as well as other conventional schemes. It is clearly shown that the present formulation can provide more accurate increasing factor than the conventional schemes. In fact, the increasing factor (20a) can be approximated as
(23)r=1+Rh+12Rh2+14(1-1ω)Rh3+18(1+1ω2-2ω)Rh4+⋯,
which is second-order accurate representation as compared with the exact one, rex. Meanwhile, third-order accurate approximation can be obtained with ω=3.Figure 2
Increasing factors.


### 3.3. Solution of a Modeled Problem

Consider the modeled problem with the following source term:(24a)S(x)={ρuhxM-h2≤x≤xM+h20otherwise,
with boundary conditions:
(24b)Φ(0)=Φ(1)=0.That is, the source term exists within the computational cell located atx=xM and its strength increases as the grid spacing is refined. The total amount of source terms remain unchanged. This equation will demonstrate effects of existing source on solution characteristics. Exact and numerical solutions of this problem can be easily secured with simple arithmetic manipulations. Figure 3 depicts the resulting solution at the source cell ϕM with respect to cell Reynolds number with various grid sizes. Generally speaking, the present formulation can provide more accurate results than UD and 2UD, especially when Rh is large. In fact, as Rh→∞, it is easy to prove that UD and 2UD will approach to 1 and 2/3, respectively, while both the exact solution and present formulation will approach to 1/2.Solution of modeled problem.
(a)
(b)
(c)


## 3.1. Interpretation of Solution Gradient at Cell Boundary

In the one-dimensional linear problem, solution gradient at cell boundary(10b) can be simplified as (16a)ϕx,i+1/2=hihi+hi+1pi+hi+1hi+hi+1pi+1+2ωhi+hi+1[ϕ(xi+1/2;i+1)-ϕ(xi+1/2;i)].
It can be interpreted as the volume averaged between two adjacent cells added by a penalty term arising from solution discontinuity. Therefore, besides the integration factor to yield integration volume, the scheme parameter ω can be regarded as a penalty factor in this expression. Meanwhile, it can also be rewritten as
(16b)ϕx,i+1/2=(1-ω)(hihi+hi+1pi+hi+1hi+hi+1pi+1)+ω(ϕi+1-ϕi)hi/2+hi+1/2.
This is a weighting average between those interpolated from solution variable and its gradient at computational nodes and ω stands for the weighting factor. Furthermore, this solution gradient can be further rearranged as:
(16c)ϕx,i+1/2=ϕ(xi+1/2+hi+1/2ω;i+1)-ϕ(xi+1/2-hi/2ω;i)hi/2ω+hi+1/2ω.This is a simple two-point representation andω can then be designated as distance factor.


## 3.2. Increasing Factor

Consider the case with constant positive convective velocity and viscosity in a uniform grid system (hi=h). Difference equations can then be reduced as
(17)BΘi-1+AΘi+CΘi+1=di,
with the solution vector, Θ=(ϕ,ph/2)T, and coefficient matrixes,
(18)B=(-Rh-ω-Rh+(1-ω)00),A=(Rh+2ωRhωRh+1+ω),C=(-ω-(1-ω)-ω-(1-ω)),di=(h2μs-ih24μ(si+1/4+2s-i-si-1/4)),
where Rh=ρuh/μ is the cell Reynolds number. In the case of ω=1, the solution gradient can be explicitly derived as
(19)pi,ω=1=2Rh+2[ϕi+1-ϕih+h4μ(si+1/4+2s-i-si-1/4)].
This expression clearly depicts that the solution gradient will adapt itself to flow condition (cell Reynolds number and source term). The resulting characteristic (nontrivial) increasing factor of our proposed scheme can be derived as [14](20a)r=Rh2+Rh+ω(2+Rh)Rh+ω(2-Rh).
Obviously, the increasing factor will be dependent on the scheme factor, ω. Meanwhile, the critical cell Reynolds number, Rh,cr, for occurrence of possible oscillatory solution (r=0) is
(20b)Rh,cr=2ωω-1.From the viewpoint of increasing factor, it is interesting to note that UD and CD schemes can be recovered by the proposed scheme withω=0 and ω=∞, respectively. That is, if ω is assigned with a small value, the scheme may become too diffusive as the UD scheme, whereas, if it is assigned with a large value, the resulting solution may be prone to being contaminated by oscillatory distribution as the CD scheme.To determine the possible range of the scheme factor, we turn to the exact increasing factor,rex=eRh,  which determines the exact scheme factor
(21)ωex(Rh)=eRhRh-Rh2-RheRh(Rh-2)+Rh+2.
This is a monotone distribution with ωex(0)=3 and ωex(∞)=1. Therefore, distribution range of this exact scheme factor is
(22)1≤ωex≤3.
This constraint is selected as the bound for scheme factor in our study. Figure 2 illustrates the distribution of increasing factors for ω=1, 2 and 3 as well as other conventional schemes. It is clearly shown that the present formulation can provide more accurate increasing factor than the conventional schemes. In fact, the increasing factor (20a) can be approximated as
(23)r=1+Rh+12Rh2+14(1-1ω)Rh3+18(1+1ω2-2ω)Rh4+⋯,
which is second-order accurate representation as compared with the exact one, rex. Meanwhile, third-order accurate approximation can be obtained with ω=3.Figure 2
Increasing factors.


## 3.3. Solution of a Modeled Problem

Consider the modeled problem with the following source term:(24a)S(x)={ρuhxM-h2≤x≤xM+h20otherwise,
with boundary conditions:
(24b)Φ(0)=Φ(1)=0.That is, the source term exists within the computational cell located atx=xM and its strength increases as the grid spacing is refined. The total amount of source terms remain unchanged. This equation will demonstrate effects of existing source on solution characteristics. Exact and numerical solutions of this problem can be easily secured with simple arithmetic manipulations. Figure 3 depicts the resulting solution at the source cell ϕM with respect to cell Reynolds number with various grid sizes. Generally speaking, the present formulation can provide more accurate results than UD and 2UD, especially when Rh is large. In fact, as Rh→∞, it is easy to prove that UD and 2UD will approach to 1 and 2/3, respectively, while both the exact solution and present formulation will approach to 1/2.Solution of modeled problem.
(a)
(b)
(c)


## 4. Fluid Flow Calculations

### 4.1. Governing Equations

For convenience, the steady incompressible Navier-Stokes equations are expressed by the following conservation law form:(25a)∇·f^c=0,(25b)∇·f^u=-∂p∂x,(25c)∇·f^v=-∂p∂y,where mass,x- and y-momentum flux functions read
(26)f^c=ρv,f^u=ρvu-μ∇u,f^v=ρvv-μ∇v.
The supplementary equations for solution gradients are derived from the original ones:(27a)∇·∇f^c=0,(27b)∇·∇f^u=-∂∇p∂x,(27c)∇·∇f^v=-∂∇p∂y.To derive the difference counterparts for these equations, the finite volume method is also employed by integrating the supplementary equations over their control volumes. However, apart from those for original equations, the control volumes for these equations are specified in a different way.


### 4.2. Domain Discretization and Approximate Solutions

The interested domain is divided into nonoverlapping polygons in the finite volume formulation. Figure4(a) shows some typical polygons delegated by a set of irregular grid nodes. These polygons are regarded as scalar cells (VC) and their centroids the scalar points. In our formulation, these centroids are assigned to be the computational locations for pressure, which are also designated as pressure nodes. Therefore, the control volume for continuity equation is coincident with the scalar cell. To prevent the occurrence of checkerboard pressure filed, we adopt a staggered polygonal grid system as shown in Figure 4(b). The computational locations for velocity vector are staggered to those for pressure (pressure nodes). They are defined at the boundary face centers of scalar cells. For brevity, these locations are then designated as vector points or velocity nodes which are also demonstrated in Figures 4(a) and 4(b). The control volume associated with a velocity node is defined by connecting neighboring grid and pressure nodes. Figure 4(b) depicts a typical control volume for a velocity node, which is the region enclosed by four dashed lines connecting  C1-N2-C2-N1-C1 surrounding the velocity node, E. This control volume is regarded as the vector cell (VE) for momentum equation. From this particular arrangement, there are four boundary faces and two neighboring pressure nodes for a vector cell irrespective of the initial grid distribution.Staggered polygonal grid.
(a)
Control volumes for scalar variable
(b)
Control volumes for vector variableTo have a second-order accurate discretization, the approximate solutions for pressure and velocity vector are mimicked by a combination of piecewise linear distributions:(28a)p(r)=∑C=1nCp(r;C),(28b)v(r)=∑E=1nEv(r;E),
with
(28c)p(r;C)={pC+∇pC·(r-rC)ifr∈VC0otherwise,(28d)v(r;E)={vE+∇vE·(r-rE)ifr∈VE0otherwise.Following the derivation procedure given for a convection-diffusion-source equation, the resulting difference equations for fluid flow calculations read(29a)∑E(C)(ρAxU+ρAyV)E(C)=0,(29b)CEUE+∑Ee(E)CEe(E)UEe(E)=∑C(E)DC(E)UPC(E),(29c)CEVE+∑Ee(E)CEe(E)VEe(E)=∑C(E)DC(E)VPC(E),with unknown variable vectors:(30)P=(p,px,py)T,U=(u,ux,uy)T,V=(v,vx,vy)T.
In the previous difference equations, Ax and Ay are the x- and y-direction components of surface vector A=Axi+Ayj, respectively. In addition, C(E) and Ee(E) are the collection set of neighboring scalar and vector nodes for vector node, E, respectively, while E(C) is the set of neighboring vector nodes for scalar node, C. As a reference shown in Figure 4(a), E(C) = {E1,E2,E3,E4,E5}  and in Figure 4(b), C(E) = {C1,C2} and Ee(E) = {E1,E2,E3,E4}.


### 4.3. Pressure-Velocity Coupling and Solution Procedure

In the incompressible flow calculations, a pressure-velocity coupling relation must be established to proceed with the solution procedure [15]. In this study, we employ the SIMPLE algorithm to circumvent this difficulty. However, the conventional SIMPLE algorithm is used for pressure-velocity coupling; it must be extended to include their derivatives as well in our formulation. First, all the solution variables and their derivatives are expressed with a tentative quantity added by a correction term:(31a)UE=U~E+UE′,VE=V~E+VE′,PC=P~C+PC′.
In general, this tentative field (U~,V~,P~) does not satisfy the continuity constraint and needed to be corrected by the correction term (U′,V′,P′). The pressure-velocity coupling is built by the difference momentum equations (29b) and (29c) via their correction terms:
(31b)UE′=CE-1∑C(E)DC(E)UPC(E)′,VE′=CE-1∑C(E)DC(E)VPC(E)′.Substituting these pressure-velocity correction relations into difference continuity equation (29a) will lead to the pressure correction equation:
(32)∑E(C)∑C(E)(ρAE(C)xCE(C)-1DC(E)V+ρAE(C)yCE(C)-1DC(E)U)PC(E)′=-∑E(C)(ρAxU~+ρAyV~)E(C).
With the given difference momentum and pressure correction equations, simulation can be carried out with the standard SIMPLE procedure.


## 4.1. Governing Equations

For convenience, the steady incompressible Navier-Stokes equations are expressed by the following conservation law form:(25a)∇·f^c=0,(25b)∇·f^u=-∂p∂x,(25c)∇·f^v=-∂p∂y,where mass,x- and y-momentum flux functions read
(26)f^c=ρv,f^u=ρvu-μ∇u,f^v=ρvv-μ∇v.
The supplementary equations for solution gradients are derived from the original ones:(27a)∇·∇f^c=0,(27b)∇·∇f^u=-∂∇p∂x,(27c)∇·∇f^v=-∂∇p∂y.To derive the difference counterparts for these equations, the finite volume method is also employed by integrating the supplementary equations over their control volumes. However, apart from those for original equations, the control volumes for these equations are specified in a different way.


## 4.2. Domain Discretization and Approximate Solutions

The interested domain is divided into nonoverlapping polygons in the finite volume formulation. Figure4(a) shows some typical polygons delegated by a set of irregular grid nodes. These polygons are regarded as scalar cells (VC) and their centroids the scalar points. In our formulation, these centroids are assigned to be the computational locations for pressure, which are also designated as pressure nodes. Therefore, the control volume for continuity equation is coincident with the scalar cell. To prevent the occurrence of checkerboard pressure filed, we adopt a staggered polygonal grid system as shown in Figure 4(b). The computational locations for velocity vector are staggered to those for pressure (pressure nodes). They are defined at the boundary face centers of scalar cells. For brevity, these locations are then designated as vector points or velocity nodes which are also demonstrated in Figures 4(a) and 4(b). The control volume associated with a velocity node is defined by connecting neighboring grid and pressure nodes. Figure 4(b) depicts a typical control volume for a velocity node, which is the region enclosed by four dashed lines connecting  C1-N2-C2-N1-C1 surrounding the velocity node, E. This control volume is regarded as the vector cell (VE) for momentum equation. From this particular arrangement, there are four boundary faces and two neighboring pressure nodes for a vector cell irrespective of the initial grid distribution.Staggered polygonal grid.
(a)
Control volumes for scalar variable
(b)
Control volumes for vector variableTo have a second-order accurate discretization, the approximate solutions for pressure and velocity vector are mimicked by a combination of piecewise linear distributions:(28a)p(r)=∑C=1nCp(r;C),(28b)v(r)=∑E=1nEv(r;E),
with
(28c)p(r;C)={pC+∇pC·(r-rC)ifr∈VC0otherwise,(28d)v(r;E)={vE+∇vE·(r-rE)ifr∈VE0otherwise.Following the derivation procedure given for a convection-diffusion-source equation, the resulting difference equations for fluid flow calculations read(29a)∑E(C)(ρAxU+ρAyV)E(C)=0,(29b)CEUE+∑Ee(E)CEe(E)UEe(E)=∑C(E)DC(E)UPC(E),(29c)CEVE+∑Ee(E)CEe(E)VEe(E)=∑C(E)DC(E)VPC(E),with unknown variable vectors:(30)P=(p,px,py)T,U=(u,ux,uy)T,V=(v,vx,vy)T.
In the previous difference equations, Ax and Ay are the x- and y-direction components of surface vector A=Axi+Ayj, respectively. In addition, C(E) and Ee(E) are the collection set of neighboring scalar and vector nodes for vector node, E, respectively, while E(C) is the set of neighboring vector nodes for scalar node, C. As a reference shown in Figure 4(a), E(C) = {E1,E2,E3,E4,E5}  and in Figure 4(b), C(E) = {C1,C2} and Ee(E) = {E1,E2,E3,E4}.


## 4.3. Pressure-Velocity Coupling and Solution Procedure

In the incompressible flow calculations, a pressure-velocity coupling relation must be established to proceed with the solution procedure [15]. In this study, we employ the SIMPLE algorithm to circumvent this difficulty. However, the conventional SIMPLE algorithm is used for pressure-velocity coupling; it must be extended to include their derivatives as well in our formulation. First, all the solution variables and their derivatives are expressed with a tentative quantity added by a correction term:(31a)UE=U~E+UE′,VE=V~E+VE′,PC=P~C+PC′.
In general, this tentative field (U~,V~,P~) does not satisfy the continuity constraint and needed to be corrected by the correction term (U′,V′,P′). The pressure-velocity coupling is built by the difference momentum equations (29b) and (29c) via their correction terms:
(31b)UE′=CE-1∑C(E)DC(E)UPC(E)′,VE′=CE-1∑C(E)DC(E)VPC(E)′.Substituting these pressure-velocity correction relations into difference continuity equation (29a) will lead to the pressure correction equation:
(32)∑E(C)∑C(E)(ρAE(C)xCE(C)-1DC(E)V+ρAE(C)yCE(C)-1DC(E)U)PC(E)′=-∑E(C)(ρAxU~+ρAyV~)E(C).
With the given difference momentum and pressure correction equations, simulation can be carried out with the standard SIMPLE procedure.


## 5. Numerical Verification

### 5.1. One-Dimensional Problems

The first test problem to be solved is the linear convectiondiffusion without source term with boundary conditions(33)Φ(0)=0,Φ(1)=1.
Figure 5 illustrates the maximum computational error with ρuℓ/μ=500 and characteristic length ℓ=1. As is shown in this figure, the computational results confirm the numerical analyses. Implementation of a constant scheme factor of  1≤ω<3 can yield a second-order accurate solution. Solution with ω=3 becomes third-order accurate. In the range of large cell Reynolds number, calculations with ω=2 and ω=3 depict larger errors due to the occurrence of oscillatory solution. However, these errors are still less than those for the conventional CD and QUICK schemes. Therefore, we can summarize that the present formulation can provide much more accurate numerical solution than the conventional schemes.Figure 5
Solution error in linear convection-diffusion problem.The second test problem considered in this category is a steady Burger’s equation:(34a)
d
d
x
(
1
2
Φ
2
-
d
Φ
d
x
)
=
0
.
The corresponding exact solution reads
(34b)Φ(x)=-ϕctanh(ϕcx2).The calculated profiles of solution variable and its spatial gradient along with the exact solution are shown in Figure6. The calculations are preformed with 40 computational cells. In this figure, the normalizing value, ϕo,  is defined as the solution at the left boundary, x=-1, or ϕo=Φ(-1). Formulation with 1≤ω≤3 can provide almost identically well-predicted solution profiles as well as its spatial gradient. The computational errors of solution variable with respect to characteristic cell Reynolds number, Rc=ϕch, are demonstrated in Figure 6(b) for ϕc=500. It clearly shows that the proposed scheme can achieve second-order accurate simulation.Solution and error distributions in one-dimensional second test problem.
(a)
Solution
(b)
ErrorThe last test problem is a modeled equation to simulate water transport in membrane [16], which reads(35a)ddx(Φ-ΦdΦdx)=0.
In this model equation, the nonlinear effect is caused by the diffusion coefficient. The interested domain is defined in the range of ℓn2≤x≤1 and the boundary conditions are assigned as
(35b)Φ(ℓn2)=0,Φ(1)=1.
The exact solution for this equation can be obtained as
(35c)x=Φ-ℓn(1+Φ)+ℓn2.Difficulty in simulation arises due to the solution singular behavior at the inlet boundary,x=ℓn2, where the value of solution variable vanishes and its spatial gradient becomes infinite. Figure 7 demonstrates the computational profiles of solution variable and its spatial gradient. These calculations are performed with 20 computational cells. It is shown that all the scheme factors can yield very accurate results.Solution distribution in one-dimensional third test problem.
(a)
(b)


### 5.2. Two-Dimensional Scalar Problems

To investigate the effects of grid distribution, especially in unstructured grid arrangements, we adopted four types of grid arrangement for comparison. These are rectangular (mesh A), regularly triangular grid (mesh B), unstructured triangular grid (mesh C), and randomly quadrilateral grid (mesh D). The regularly triangular grid is derived by dividing a uniform rectangular cell into two triangles and the unstructured triangular grid is generated by the advancing front method [17]. The randomly distributed quadrilateral grid is also derived from a rectangular grid by randomly varying its grid locations:
(36)x=xrec+α(ψx-0.5)Δx,y=yrec+α(ψy-0.5)Δy,
where xrec and yrec are the original grid location coordinates in rectangular grid, ψx and ψy are the random numbers in the range of  0<ψx,ψy<1, and α is a grid distortion parameter. Figure 8(a) shows a typical 40×40 rectangular grid, while Figure 8(b) a regularly triangular grid derived from this rectangular mesh. Figure 8(c) illustrates an unstructured triangular grid system. This grid system is generated based on uniformly distributed grid density same as in Figures 8(a) and 8(b). The resulting randomly quadrilateral grid with α=0.8 is displayed in Figure 8(d).Four types of grid arrangement.
(a)
Mesh A
(b)
Mesh B
(c)
Mesh C
(d)
Mesh DThe first test problem is a nonlinear convection-diffusion problem:(37a)Φcosθ∂Φ∂x+Φsinθ∂Φ∂y-μ(∂2Φ∂x2+∂2Φ∂y2)=0,
with the exact solution
(37b)Φ(x,y)=-tanh(xcosθ+ysinθ2μ).Fluid viscosity ofμ=0.02 is employed. Figure 9(a) depicts the solution contours imbedded with the adopted grid arrangement (mesh D), which can be compared with the exact solution in Figure 9(b). As evident in these figures, the computational results are in good agreement with the exact solution. This observation can be reinforced by the comparisons of solution profiles along y=0, which are shown in Figures 10(a) and 10(b). In these figures, both the solution variable and its spatial gradients obtained with various grid arrangements are compared with the exact solution. It is clearly shown that well-predicted results can be obtained regardless of grid arrangement.Solution contours in two-dimensional first test problem.
(a)
Mesh D
(b)
ExactSolution profile alongy=0 in two-dimensional first test problem.
(a)
Solution
(b)
GradientThe second problem is a linear convection-diffusion problem:(38a)cosθ∂Φ∂x+sinθ∂Φ∂y-μ(∂2Φ∂x2+∂2Φ∂y2)=0,
which may simulate cross-stream diffusion effects. If the streamwise diffusion is neglected, we can obtain the approximate exact solution
(38b)Φ(x,y)=erf[12μ-xsinθ+ycosθ(x+1)cosθ+(y+1)sinθ+εd],where erf[] is the error function. In our computations, Fluid viscosity of μ=0.01 is employed. Meanwhile, εd=10-3 is assigned to avoid the possible solution singularity. Figure 11(a) demonstrates the computational solution contours accompanied by adopted grid arrangement (mesh D) and Figure 11(b) depicts the approximate exact solution for θ=45°. Detailed solution profiles along x=0 for the computational and approximate exact solutions are displayed in Figures 12(a) and 12(b). As shown in these figures, very agreed results can be obtained for all grid arrangements.Solution contours in two-dimensional second test problem.
(a)
Mesh D
(b)
ExactSolution profile alongx=0 in two-dimensional second test problem.
(a)
Solution
(b)
GradientThe last test problem considered in this study is a convection-diffusion problem in a rotating flow field:(39a)u∂Φ∂x+v∂Φ∂y-μ(∂2Φ∂x2+∂2Φ∂y2)=S,
with the velocity field
(39b)(u,v)=(y,-x).
Source term in the governing equation (39a) is assigned to yield the exact solution
(39c)Φ(x,y)=ysin(πx).Computations are conducted with dimensionless fluid viscosity ofμ=0.01. Figure 13(a) shows the solution contours with the adopted grid system (mesh D), while Figure 13(b) depicts the exact solution contours for comparison. Comparisons of the computational solution profiles along the line of x-y=0 on various types of grid arrangements and exact solution are shown in Figure 14. From these figures, it is clearly shown that both the solution variable and its gradients can yield reasonable results as compared with the exact solution.Solution contours in two-dimensional second test problem.
(a)
Mesh D
(b)
ExactSolution profile alongx=y in the two-dimensional third test problem.
(a)
(b)
(c)


### 5.3. Lid-Driven Cavity Flow

The flow Reynolds number based on driving velocity, fluid viscosity, and cavity width isReld=400 in present computations. To investigate the effects of grid distribution, four types of grid arrangement are employed. They are rectangular (mesh R), randomly quadrilateral grid (mesh R′), regularly triangular grid (mesh T), and randomly triangular grid (mesh T′). Figure 15 illustrates the convergence history for a typical computation with 40×40 uniform rectangular grid (mesh R). Equation residuals of the mass, x-momentum (u) and y-momentum (v) equations are depicted in this figure. From the monotonely decreasing equation residues depicted in this figure, it is clearly shown that the present formulation can provide a very stable computation. The random meshes (R′ and T′) are obtained with the grid distortion parameter α=0.5. With this mesh density, Figures 16(a), 16(b), 16(c), and 16(d) demonstrate the resulting streamlines as well as the adopted grid arrangements for mesh R, R′, T, and T′, respectively. From these figures, quite agreed streamlines can be observed irrespective  of grid arrangement. This observation can be reinforced by the velocity profile distributions across the cavity geometrical center in Figures 17(a) and 17(b). All computational results are in very good agreement with the benchmark solution [18]. Figure 18 illustrates the vortex strength in terms of adopted grid distortion. Only slight variations in vortex strength are observed for 0≤α≤0.5.Figure 15
Convergence history of a typical lip-driven cavity flow simulation.Grid arrangement and streamlines in lip-driven cavity flow computations.
(a)
Mesh R
(b)
Mesh R′
(c)
Mesh T
(d)
Mesh T′Velocity profile across cavity geometrical center.
(a)
x-component
(b)
y-componentFigure 18
Effects of grid distortion on vortex strength.


### 5.4. Backward Facing Step Flow

The second test problem considered in the present study is a backward facing step flow with expansion ratio,(hin+hstep)/hin=1.942, where hin and hstep are the inlet and step heights, respectively. In our computations, the grid distribution in the x-direction is defined by the following rule:
(40)0≤x≤15hstepuniformΔx,15hstep<x≤30hstepeexpansionΔx,30hstep<x≤60hstepuniformΔx.
The initial grid spacing is set as Δx=0.2hstep and expansion ratio is 1.1. Uniform grid is set in the y-direction with Δy=0.1hstep. Figure 19 shows the convergence history of a typical flow simulation, which demonstrates a very stable computation. The resulting streamlines of computations with Reynolds numbers of  Restep=100 and 389 are displayed in Figures 20(a) and 20(b), respectively. This Reynolds number is defined by the inlet mean velocity, fluid viscosity, and inlet hydraulic diameter. Comparisons with available experimental data [19] on axial velocity are shown in Figures 21(a) and 21(b) for Restep=100 and 389, respectively. In these figures, the axial velocity has been normalized by the mean velocity in the duct, um. Very good agreements are observed between computational and experimental results. Figure 22 illustrates the resulting reattachment length in our simulation along with experimental and other numerical results for comparison. The numerical results provided by Erturk [20] were conducted with a very fine grid, Δy=0.01(hin+hstep). As shown in this figure, quite agreed results can be observed with the present formulation.Figure 19
Convergence history of a typical backward facing step flow simulation.Streamlines in backward facing step flow.
(a)
R
e
step
=
100
(b)
R
e
step
=
389Comparison of axial velocity profiles.
(a)
R
e
step
=
100
(b)
R
e
step
=
389Figure 22
Comparison of reattachment length.


## 5.1. One-Dimensional Problems

The first test problem to be solved is the linear convectiondiffusion without source term with boundary conditions(33)Φ(0)=0,Φ(1)=1.
Figure 5 illustrates the maximum computational error with ρuℓ/μ=500 and characteristic length ℓ=1. As is shown in this figure, the computational results confirm the numerical analyses. Implementation of a constant scheme factor of  1≤ω<3 can yield a second-order accurate solution. Solution with ω=3 becomes third-order accurate. In the range of large cell Reynolds number, calculations with ω=2 and ω=3 depict larger errors due to the occurrence of oscillatory solution. However, these errors are still less than those for the conventional CD and QUICK schemes. Therefore, we can summarize that the present formulation can provide much more accurate numerical solution than the conventional schemes.Figure 5
Solution error in linear convection-diffusion problem.The second test problem considered in this category is a steady Burger’s equation:(34a)
d
d
x
(
1
2
Φ
2
-
d
Φ
d
x
)
=
0
.
The corresponding exact solution reads
(34b)Φ(x)=-ϕctanh(ϕcx2).The calculated profiles of solution variable and its spatial gradient along with the exact solution are shown in Figure6. The calculations are preformed with 40 computational cells. In this figure, the normalizing value, ϕo,  is defined as the solution at the left boundary, x=-1, or ϕo=Φ(-1). Formulation with 1≤ω≤3 can provide almost identically well-predicted solution profiles as well as its spatial gradient. The computational errors of solution variable with respect to characteristic cell Reynolds number, Rc=ϕch, are demonstrated in Figure 6(b) for ϕc=500. It clearly shows that the proposed scheme can achieve second-order accurate simulation.Solution and error distributions in one-dimensional second test problem.
(a)
Solution
(b)
ErrorThe last test problem is a modeled equation to simulate water transport in membrane [16], which reads(35a)ddx(Φ-ΦdΦdx)=0.
In this model equation, the nonlinear effect is caused by the diffusion coefficient. The interested domain is defined in the range of ℓn2≤x≤1 and the boundary conditions are assigned as
(35b)Φ(ℓn2)=0,Φ(1)=1.
The exact solution for this equation can be obtained as
(35c)x=Φ-ℓn(1+Φ)+ℓn2.Difficulty in simulation arises due to the solution singular behavior at the inlet boundary,x=ℓn2, where the value of solution variable vanishes and its spatial gradient becomes infinite. Figure 7 demonstrates the computational profiles of solution variable and its spatial gradient. These calculations are performed with 20 computational cells. It is shown that all the scheme factors can yield very accurate results.Solution distribution in one-dimensional third test problem.
(a)
(b)


## 5.2. Two-Dimensional Scalar Problems

To investigate the effects of grid distribution, especially in unstructured grid arrangements, we adopted four types of grid arrangement for comparison. These are rectangular (mesh A), regularly triangular grid (mesh B), unstructured triangular grid (mesh C), and randomly quadrilateral grid (mesh D). The regularly triangular grid is derived by dividing a uniform rectangular cell into two triangles and the unstructured triangular grid is generated by the advancing front method [17]. The randomly distributed quadrilateral grid is also derived from a rectangular grid by randomly varying its grid locations:
(36)x=xrec+α(ψx-0.5)Δx,y=yrec+α(ψy-0.5)Δy,
where xrec and yrec are the original grid location coordinates in rectangular grid, ψx and ψy are the random numbers in the range of  0<ψx,ψy<1, and α is a grid distortion parameter. Figure 8(a) shows a typical 40×40 rectangular grid, while Figure 8(b) a regularly triangular grid derived from this rectangular mesh. Figure 8(c) illustrates an unstructured triangular grid system. This grid system is generated based on uniformly distributed grid density same as in Figures 8(a) and 8(b). The resulting randomly quadrilateral grid with α=0.8 is displayed in Figure 8(d).Four types of grid arrangement.
(a)
Mesh A
(b)
Mesh B
(c)
Mesh C
(d)
Mesh DThe first test problem is a nonlinear convection-diffusion problem:(37a)Φcosθ∂Φ∂x+Φsinθ∂Φ∂y-μ(∂2Φ∂x2+∂2Φ∂y2)=0,
with the exact solution
(37b)Φ(x,y)=-tanh(xcosθ+ysinθ2μ).Fluid viscosity ofμ=0.02 is employed. Figure 9(a) depicts the solution contours imbedded with the adopted grid arrangement (mesh D), which can be compared with the exact solution in Figure 9(b). As evident in these figures, the computational results are in good agreement with the exact solution. This observation can be reinforced by the comparisons of solution profiles along y=0, which are shown in Figures 10(a) and 10(b). In these figures, both the solution variable and its spatial gradients obtained with various grid arrangements are compared with the exact solution. It is clearly shown that well-predicted results can be obtained regardless of grid arrangement.Solution contours in two-dimensional first test problem.
(a)
Mesh D
(b)
ExactSolution profile alongy=0 in two-dimensional first test problem.
(a)
Solution
(b)
GradientThe second problem is a linear convection-diffusion problem:(38a)cosθ∂Φ∂x+sinθ∂Φ∂y-μ(∂2Φ∂x2+∂2Φ∂y2)=0,
which may simulate cross-stream diffusion effects. If the streamwise diffusion is neglected, we can obtain the approximate exact solution
(38b)Φ(x,y)=erf[12μ-xsinθ+ycosθ(x+1)cosθ+(y+1)sinθ+εd],where erf[] is the error function. In our computations, Fluid viscosity of μ=0.01 is employed. Meanwhile, εd=10-3 is assigned to avoid the possible solution singularity. Figure 11(a) demonstrates the computational solution contours accompanied by adopted grid arrangement (mesh D) and Figure 11(b) depicts the approximate exact solution for θ=45°. Detailed solution profiles along x=0 for the computational and approximate exact solutions are displayed in Figures 12(a) and 12(b). As shown in these figures, very agreed results can be obtained for all grid arrangements.Solution contours in two-dimensional second test problem.
(a)
Mesh D
(b)
ExactSolution profile alongx=0 in two-dimensional second test problem.
(a)
Solution
(b)
GradientThe last test problem considered in this study is a convection-diffusion problem in a rotating flow field:(39a)u∂Φ∂x+v∂Φ∂y-μ(∂2Φ∂x2+∂2Φ∂y2)=S,
with the velocity field
(39b)(u,v)=(y,-x).
Source term in the governing equation (39a) is assigned to yield the exact solution
(39c)Φ(x,y)=ysin(πx).Computations are conducted with dimensionless fluid viscosity ofμ=0.01. Figure 13(a) shows the solution contours with the adopted grid system (mesh D), while Figure 13(b) depicts the exact solution contours for comparison. Comparisons of the computational solution profiles along the line of x-y=0 on various types of grid arrangements and exact solution are shown in Figure 14. From these figures, it is clearly shown that both the solution variable and its gradients can yield reasonable results as compared with the exact solution.Solution contours in two-dimensional second test problem.
(a)
Mesh D
(b)
ExactSolution profile alongx=y in the two-dimensional third test problem.
(a)
(b)
(c)


## 5.3. Lid-Driven Cavity Flow

The flow Reynolds number based on driving velocity, fluid viscosity, and cavity width isReld=400 in present computations. To investigate the effects of grid distribution, four types of grid arrangement are employed. They are rectangular (mesh R), randomly quadrilateral grid (mesh R′), regularly triangular grid (mesh T), and randomly triangular grid (mesh T′). Figure 15 illustrates the convergence history for a typical computation with 40×40 uniform rectangular grid (mesh R). Equation residuals of the mass, x-momentum (u) and y-momentum (v) equations are depicted in this figure. From the monotonely decreasing equation residues depicted in this figure, it is clearly shown that the present formulation can provide a very stable computation. The random meshes (R′ and T′) are obtained with the grid distortion parameter α=0.5. With this mesh density, Figures 16(a), 16(b), 16(c), and 16(d) demonstrate the resulting streamlines as well as the adopted grid arrangements for mesh R, R′, T, and T′, respectively. From these figures, quite agreed streamlines can be observed irrespective  of grid arrangement. This observation can be reinforced by the velocity profile distributions across the cavity geometrical center in Figures 17(a) and 17(b). All computational results are in very good agreement with the benchmark solution [18]. Figure 18 illustrates the vortex strength in terms of adopted grid distortion. Only slight variations in vortex strength are observed for 0≤α≤0.5.Figure 15
Convergence history of a typical lip-driven cavity flow simulation.Grid arrangement and streamlines in lip-driven cavity flow computations.
(a)
Mesh R
(b)
Mesh R′
(c)
Mesh T
(d)
Mesh T′Velocity profile across cavity geometrical center.
(a)
x-component
(b)
y-componentFigure 18
Effects of grid distortion on vortex strength.


## 5.4. Backward Facing Step Flow

The second test problem considered in the present study is a backward facing step flow with expansion ratio,(hin+hstep)/hin=1.942, where hin and hstep are the inlet and step heights, respectively. In our computations, the grid distribution in the x-direction is defined by the following rule:
(40)0≤x≤15hstepuniformΔx,15hstep<x≤30hstepeexpansionΔx,30hstep<x≤60hstepuniformΔx.
The initial grid spacing is set as Δx=0.2hstep and expansion ratio is 1.1. Uniform grid is set in the y-direction with Δy=0.1hstep. Figure 19 shows the convergence history of a typical flow simulation, which demonstrates a very stable computation. The resulting streamlines of computations with Reynolds numbers of  Restep=100 and 389 are displayed in Figures 20(a) and 20(b), respectively. This Reynolds number is defined by the inlet mean velocity, fluid viscosity, and inlet hydraulic diameter. Comparisons with available experimental data [19] on axial velocity are shown in Figures 21(a) and 21(b) for Restep=100 and 389, respectively. In these figures, the axial velocity has been normalized by the mean velocity in the duct, um. Very good agreements are observed between computational and experimental results. Figure 22 illustrates the resulting reattachment length in our simulation along with experimental and other numerical results for comparison. The numerical results provided by Erturk [20] were conducted with a very fine grid, Δy=0.01(hin+hstep). As shown in this figure, quite agreed results can be observed with the present formulation.Figure 19
Convergence history of a typical backward facing step flow simulation.Streamlines in backward facing step flow.
(a)
R
e
step
=
100
(b)
R
e
step
=
389Comparison of axial velocity profiles.
(a)
R
e
step
=
100
(b)
R
e
step
=
389Figure 22
Comparison of reattachment length.


## 6. Conclusion

In the present study, we have successfully developed a feasible solution procedure to simulate incompressible fluid flow in complicated domains. It is put into effect by discretizing the incompressible Navier-Stokes equations with primitive variables on staggered polygonal grids. The discretization procedure of its kernel scheme has been unveiled for convection-diffusion equation. Its main feature is to directly calculate variable gradients with their own governing equation without further interpolation from solution values at neighboring computational nodes. This process will yield a high-resolution discretization with compact stencil in difference equations. Meanwhile, the adoption of staggered grid will eliminate the occurrence of checkerboard pressure field in the resulting solution for its own sake. Pressure-velocity coupling problem in incompressible flow calculation is resolved by the SIMPLE algorithm. Finally, the proposed solution procedure is employed to solve several classic benchmark problems. From the computational results, it can be concluded that the present formulation will be one of the feasible tools to simulate incompressible flow problems.



---

*Source: 409846-2013-06-02.xml*