# On Numerical Solution of the Incompressible Navier-Stokes Equations with Static or Total Pressure Specified on Boundaries

Authors: N. P. Moshkin; D. Yambangwai
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

The purpose of this article is to develop and validate a computational method for the solution of viscous incompressible flow in a domain with specified static or total pressure on the flow-through boundaries (inflow and outflow). The computational algorithm is based on the Finite Volume Method in nonstaggered boundary-fitted grid. The implementations of the boundary conditions on the flow-through parts of the boundary are discussed. Test examples illustrate the main features and validity of the proposed method to study viscous incompressible flow through a bounded domain with specified static pressure (or total pressure) on boundary as a part of well-posed boundary conditions.


---

## Body

## 1. Introduction

A flow of a viscous incompressible fluid through a given domain is rather interesting for its numerous engineering applications. Typically, these include tube and channel flows with a variety of geometries. The difficulties in mathematical modeling and numerical simulation of such flows arise in the flow-through boundaries (inflow and outflow). If the domain of interest is completely bounded by impermeable walls, there is no ambiguity in the boundary conditions for the incompressible Navier-Stokes equations. However, when flow-through (inflow and outflow) boundaries are present, there is no general agreement on which kind of boundary conditions is both mathematically correct and physically appropriate on these flow-through boundaries. Traditionally, such problems are treated with specified velocity on the domain boundaries. However, in many applications the boundary velocities are not known; instead the pressure variation is given at the boundaries, and the flow within the domain has to be determined. For example, in the central air-conditioning or air-heating system of a building, a main supply channel branches into many subchannels that finally open into the different rooms, which can be at a different constant pressure. The distribution of the flow into various branches depends on the flow resistances of these branches, and in general, it is even impossible to predict the direction of flow.The problem of solvability and uniqueness of an initial boundary value problem for the incompressible Navier-Stokes equations is one of the various problems considered, for example, in [1–6] and many others. Major part of research deals with proper formulation of boundary conditions for pressure which are needed in numerical simulation but absence in the mathematical statement of problem (see, e.g., more recent [7, 8] and therein references). The object of our study is a boundary value problem in which the pressure is known on boundary as a part of boundary conditions in the mathematical statement of problem.Antontsev et al. [1], Ragulin [4], and Ragulin and Smagulov [5] have studied initial boundary value problems in which the values of pressure or total pressure are specified on flow-through boundaries. Ragulin [4] and Ragulin and Smagulov [5] have considered problems for the homogeneous Navier-Stokes equations. Antontsev et al. [1] have studied well-posedness of the nonhomogeneous Navier-Stokes equations. As these results are not well known, we will shortly represent the well-posed statement of initial boundary value problems with specified pressure boundaries.To the best of the authors’ knowledge, the research on numerically treated pressure boundary conditions for the incompressible Navier-Stokes equations is limited. Some of the research conducted is discussed below. Kuznetsov et al. [9] and Moshkin [10–12] developed finite difference algorithms to treat incompressible viscous flow in a domain with given pressure on flow-through parts of the boundary. Finite-difference numerical algorithms were developed for primitive variables and for stream function vorticity formulation of 2D Navier-Stokes equations.In the finite-element study by Hayes et al. [13], a brief discussion of the specified pressure on the outflow region of the boundary is presented. Kobayashi et al. [14] have discussed the role of pressure specified on open boundaries in the context of the SIMPLE algorithm.The prescription of a pressure drop between the inlet and the outlet of the flow was also considered by Heywood et al. [15], where a variational approach with given mean values of the pressure across the inflow and outflow boundaries was used.The construction of the discretized equations for unknown velocities on specified pressure boundaries and the solution of the discretized equations using the SIMPLE algorithm are discussed in [16]. The computational treatment of specified pressure boundaries in complex geometries is presented within the framework of a nonstaggered technique based on curvilinear boundary-fitted grids. The proposed method is applied for predicting incompressible forced flows in branched ducts and in buoyancy-driven flows.A finite-difference method for solving the incompressible time-dependent three-dimensional Navier–Stokes equations in open flows where Dirichlet boundary conditions for the pressure are given on part of the boundary is presented in [17]. The equations in primitive variables (velocity and pressure) are solved using a projection method on a nonstaggered grid with second-order accuracy in space and time. On the inflow and outflow boundaries the pressure is obtained from its given value at the contour of these surfaces using a two-dimensional form of the pressure Poisson equation, which enforces the incompressibility constraint ∇·v=0. The pressure obtained on these surfaces is used as Dirichlet boundary conditions for the three-dimensional Poisson equation inside the domain. The solenoidal requirement imposes some restrictions on the choice of the open surfaces.Barth and Carey [18] discussed the choice of appropriate inflow and outflow boundary conditions for Newtonian and generalized Newtonian channel flows. They came to conclusion that “…For real-world problems that are fundamentally pressure driven and involve complex geometries, it is desirable to impose a pressure drop by means of specified pressures at the inflow and outflow boundaries…” At the inflow and outflow boundaries one of the conditions specifies the normal component of the surface traction force, and the other two imply that there is no tangential flow at these boundaries; that is, flow is normal to the inflow and outflow boundaries. But no mathematical justification was given.Let us call problems where fluid can enter or leave a domain through parts of the boundary, a “flowing-through problem” for viscous incompressible fluid flow. In [17] these problems are called problems with “open” boundaries. We think that the term flowing-through problem is more suitable. The purpose of our research is not to add new insight into the mathematical statement of the problem but to develop a finite volume method for solving a flowing-through problem for the incompressible Navier-Stokes equations for which questions of existence and uniqueness have been considered in [1, 4, 5].In the following sections of this paper, a brief overview of various kinds of well-posed flowing-through problems for the incompressible Navier-Stokes equations is presented. This is followed by a description of the finite volume numerical method with strength on implementation of boundary conditions on the flow-through parts. The numerical method is then validated by a comparison of analytical and numerical solutions for the laminar flow driven by pressure drop in the2D plane channel, in the 2D gap between two cylinders, in U-bend channel, and in a planar T-junction channel.


## 2. Mathematical Formulation of Flowing-Through Problems

We present here the various kinds of well-posed flowing-through boundary value problems for the incompressible Navier-Stokes equation. In our explanation, we follow Antontsev et al. [1], Ragulin [4], and Ragulin and Smagulov [5]. Let us consider the flow of viscous liquid through bounded domain Ω of Rd (d=2 or 3), t∈[0,T], where T>0 is a fixed time. Let Γk1,k=1,…,K denote parts of the boundary Γ=∂Ω where the fluid enter or leave the domain. Let Γl0,l=1,…,L be an impermeable parts of the boundary, D=Ω×(0,T), S=Γ×(0,T), Siα=Γiα×(0,T), α=0,1. Scheme of the domain is depicted in Figure 1.Figure 1
Sketch of the flowing-through domain.The flowing-through problem is to find a solution of the Navier-Stokes equations(2.1)∂u⃗∂t+(u⃗·∇)u⃗=-1ρ∇p+νΔu⃗,∇·u⃗=0,
in the domain D=Ω×(0,T) with appropriate initial and boundary conditions, where u⃗ is the velocity vector, p is the pressure, ρ is the density, and ν is the kinematic viscosity. The initial data are
(2.2)u⃗∣t=0=u⃗0(x⃗),∇·u⃗0=0,x⃗∈Ω.
On the solid walls Γl0, the no-slip condition holds
(2.3)u⃗=0,(x⃗,t)∈Sl0,l=1,…,L.
On the flow-through parts Γk1, k=1,…,K three types of boundary conditions can be set up to make problem wellposed. As shown in [1, 4], the conditions are the followings.(i)
On the flow-through partsΓj1,j=j1,…,jn, the tangent components of the velocity vector and the total pressure are prescribed as
(2.4)u⃗·τ⃗m=Gjm(x⃗,t),m=1,2,p+12ρ|u⃗·u⃗|=Hj(x⃗,t),(x⃗,t)∈Sj1,j=j1,…,jn.
Here τ⃗1,τ⃗2 are the linearly independent vectors tangent to Γj1. Functions Gjm(x⃗,t), and Hj(x⃗,t) are given on Sj1=Γj1×(0,T).(ii)
On the flow-through partsΓl1,l=l1,…,ln, the tangent components of the velocity vector and pressure are known as
(2.5)u⃗·τ⃗m=Glm(x⃗,t),m=1,2,p=Hl(x⃗,t),(x⃗,t)∈Sl1,l=l1,…,ln.
Here Glm(x⃗,t) and Hl(x⃗,t) are given on Sl1=Γl1×(0,T).(iii)
On the flow-through partsΓs1,s=s1,…,sn, the velocity vector (all three components) has to be prescribed as
(2.6)u⃗=u⃗s1(x⃗,t),(x⃗,t)∈Ss1,s=s1,…,sn.
Here u⃗s1(x⃗,t) is given on Ss1=Γs1×(0,T).It should be mentioned that various combinations of boundary conditions onSk1,k=1,…,K give well-posed problems. For example, on the portion of the flow-through parts Sj1,j=j1,…,jn one kind of boundary condition may hold, and other portions another kinds may hold.


## 3. Finite Volume Approximation of Flowing-Through Problems

Let us present the numerical algorithm for the flowing-through problem. Numerous variation of projection methods have been developed and have been successfully utilized in computing incompressible flow problems. To emphasize on pressure boundary conditions, we used here simple explicit projection method. Although some of the main aspects are well known in literature, for the sake of completeness details are given.


### 3.1. Time Discretization

The time discretization used here is based upon the simplest projection scheme originally proposed by Chorin and Temam (see, e.g., [19, 20]). This scheme has an irreducible splitting error of order O(Δt). Hence, using a higher-order time stepping scheme for the operator ∂/∂t-νΔ does not improve overall accuracy. Using the explicit Euler time stepping, the marching steps in time are the following.Setu⃗|t=0=u⃗0, then for n≥0 compute u⃗*,u⃗n+1, and pn+1 by solving first substep:
(3.1)u⃗*-u⃗nΔt+(u⃗n·∇)u⃗n=ν△u⃗n,
and second substep:
(3.2)u⃗n+1-u⃗*Δt=-∇pn+1,(3.3)∇·u⃗n+1=0,(u⃗n+1)Γ0=0,
where Δt=T/N is the time step, N is the integer, u⃗n≈u⃗(x⃗,nΔt), and pn+1≈p(x⃗,(n+1)Δt). Without loss of generality, density is equal to one, ρ=1. The explicit approximation of convective and viscous terms in (3.1) introduces restriction on the time step for stability. This is analyzed by many (see, e.g., [20, 21] and therein references).


### 3.2. Space Discretization

For the sake of simplicity and without loosing generality, the formulation of the numerical algorithm is illustrated for a two-dimensional domain. Letu⃗=(ux,uy) be the velocity vector, where ux and uy are the Cartesian components in x and y direction, respectively. The finite volume discretization is represented for nonorthogonal quadrilaterals grid. The collocated variable arrangement is utilized. Each discrete unknown is associated with the center of control volume Ω. First, we discretize the convection and diffusion parts of the Navier-Stokes equation. One can recast (3.1) in the form
(3.4)ϕ*-ϕnΔt+∇·(ϕnu⃗n)=νΔϕn,
where the variable ϕ can be either ux or uy, and u⃗n is such that ∇·u⃗n=0.The discrete form of (3.4) is obtained by integrating on each control volume Ω, followed by the application of the Gauss theorem:
(3.5)∫Ωϕ*-ϕnΔtdΩ+∮Sϕn(u⃗n·n⃗)dS=ν∮S∇ϕn·n⃗dS,
where S is the boundary of control volume Ω (e.g., in the case shown in Figure 2, S is the union of the control volume faces s,e,n, and w), and n⃗ is the unit outward normal vector to S. Using the midpoint rule to approximation, the surface and volume integrals yield
(3.6)∫Ωϕ*-ϕnΔtdΩ≈(ϕ*-ϕnΔt)PΔΩ,(3.7)∮Sϕn(u⃗n·n⃗)dS≈∑c=e,s,n,wϕcn(u⃗n·n⃗)cSc,(3.8)∮S∇ϕn·n⃗dS=∮SDnϕndS≈∑c=e,s,n,w(Dnϕn)cSc,
where ΔΩ is the volume of control volume Ω, Sc is the area of the “c” control volume face, and (Dnϕ)c is the derivative of Cartesian velocity components in the normal direction at the center of the “c” control volume face. To estimate the right-hand side in (3.7) and (3.8), we need to know the value of Cartesian velocity components and its normal derivative on the faces of each control volume. The implementation of Cartesian velocity components on nonorthogonal grids requires special attention because the boundary of the control volume is usually not aligned with the Cartesian velocity components. The 2D interpolation of irregularly-spaced data (see, e.g., [22]) is used to interpolate Cartesian velocity components on the boundary of each control volume in (3.7). Only the east side of a 2D control volume shown in Figure 2(a) will be considered. The same approach applies to other faces, only the indices need to be changed. For example, let ϕk be the value of Cartesian velocity components at point k where k=N,P,S,SE,E,NE, and let L(e,k) be the Cartesian distance between e and k. Using 2D interpolation yields
(3.9)ϕe=(∑kL(e,k)-2ϕk)(∑kL(e,k)-2),k=N,P,S,SE,E,NE,
where L(e,k)-2=1/[(xe-xk)2+(ye-yk)2]. The derivative of Cartesian velocity components in the normal direction at the center of the control volume face in (3.8) can be calculated by using the central difference approximation (see Figure 2(a)):
(3.10)(Dnϕ)e≈ϕE′-ϕP′L(P′,E′).
The auxiliary nodes P′ and E′ lie at the intersection of the line passing through the point “e" in the direction of normal vector n⃗ and the straight lines which connect nodes P and N or E and NE, respectively, and L(P′,E′) stands for the distance between P′ and E′. The values of ϕE′ and ϕP′ can be calculated by using the gradient at control volume center:
(3.11)ϕE′=ϕE+∇ϕE·(x⃗E′-x⃗E),ϕP′=ϕP+∇ϕP·(x⃗P′-x⃗P),
where x⃗P, x⃗E, x⃗P′, and x⃗E′ are the radius vectors of P, E, P′, and E′, respectively. The kth Cartesian components of ∇ϕP are approximated using Gauss’s theorem:
(3.12)∇ϕP·i⃗k=(∂ϕ∂xk)P=1ΔΩ∑c=e,s,n,wϕcn+1Sck,Sck=Sc(n⃗·ik⃗),
where Sc is the area of “c” control volume face, n⃗ is the unit outward normal vector to Sc, and i⃗k is the unit basis vector of Cartesian coordinate system (x1,x2)=(x,y). Using (3.6)–(3.12) to approximate (3.5), one can determine velocity field u⃗* (which is not solenoidal) at each grid node, even on the boundary.A typical 2D control volume and the notation used. A way of calculating cell face values and gradients.
(a)(b)In the first substep the continuity (3.3) is not used so that the intermediate velocity field is, in general, nondivergence free. The details of the setting and discretization of the second substep developed on nonuniform, collocated grid are discussed below. Equation (3.2) applies both in continuous and discrete sense. Taking the divergence of both sides of (3.2) and integrating over a control volume Ω, after applying the Gauss theorem and setting the update velocity filed, u⃗n+1, to be divergence free, one gets the equation
(3.13)0=1ΔΩ∮Su⃗n+1·n⃗dS=1ΔΩ∮Su⃗*·n⃗dS-Δt1ΔΩ∮S∇pn+1·n⃗dS,
that has to be discretized while collocating the variables in the control volume centers. Here n⃗ is outward normal to the boundary S of control volume Ω. At this stage of the projection procedure, the discrete values of ux* and uy* are already known and represent the source term in (3.13). A second-order discretization of the surface integrals can be obtained by utilizing the mean value formula. This means that the surface integrals in (3.13) can be approximated as
(3.14)1ΔΩ∮Su⃗n+1·n⃗dS≅1ΔΩ∑c=e,s,w,n(u⃗n+1·n⃗)cSc,1ΔΩ∮S∇pn+1·n⃗dS≅1ΔΩ∑c=e,s,w,n(∇pn+1·n⃗)cSc.
It follows that by substituting (3.14) into (3.13), one gets the discrete pressure equation
(3.15)1ΔΩ∑c=e,s,w,n(u⃗*·n⃗)cSc-ΔtΔΩ∑c=e,s,w,n(Dnpn+1)cSc=0.
The iterative method is utilized to approximate (Dnpn+1)c and solve (3.15). The normal-to-face intermediate velocities (u⃗*·n⃗)c,c=e,s,w,n are not directly available. They are found using interpolation. The derivative of pressure with respect to the direction of the outward normal n⃗ through the cell face “c”, (Dnp)cn+1 is approximated by on iterative technique (see, e.g., [23]) to reach a higher order of approximation and preserved compact stencil in the discrete equation (3.15). Only the east face of a 2D control volume shown in Figure 2(a) will be considered. The same approach applies to other faces. Using second upper index “s” to denote the number of iteration, one writes
(3.16)(Dnp)en+1,s+1=(Dξp)en+1,s+1+[(Dnp)e-(Dξp)e]n+1,s,s=0,…,Ŝ,(Dnp)n+1,0=(Dnp)n,
where ξ is the direction along the line connecting nodes P and E (see Figure 2(a)). The terms in the square brackets are approximated with high order and are evaluated by using values known from the previous iteration. Once the iterations converge, the low-order approximation term (Dξp)en+1,s+1 drops out, and the solution obtained corresponds to the higher order of approximation. The derivatives of pressure in the square brackets are written as
(3.17)(Dnp)en+1,s=(∇p·n⃗)en+1,s,(Dξp)en+1,s=(∇p·ξ⃗)en+1,s,
where n⃗ is the unit outward normal vector to cell face “e”, and ξ⃗ is the unit vector in ξ direction from point P to E. The term (∇p)en+1,s is approximated similar to (3.9) as
(3.18)(∇p)en+1,s=(∑lL(e,l)-2∇pln+1,s)(∑lL(e,l)-2),l=N,P,S,SE,E,NE,
where ∇pln+1,s is the gradient of the pressure at grid node l and L(e,l)-2=1/[(xe-xl)2+(ye-yl)2]. The kth components of ∇pln+1,s are discretized by using Gauss theorem (e.g., at grid node P):
(3.19)∇pPn+1,s·i⃗k=(∂pn+1,s∂xk)P≅1ΔΩ∑c=e,s,n,wpcn+1,sSck,Sck=Sc(n⃗·ik⃗).
The first term in the right-hand side of (3.16) is treated implicitly, and a simple approximation is used (that gives a compact stencil):
(3.20)(Dξp)en+1,s+1≈pEn+1,s+1-pPn+1,s+1L(P,E),
where L(P,E) is the distance between nodes P and E. The final expression for the approximation of the derivative of pressure with respect to n⃗ through the cell face “e” can now be written as
(3.21)(Dnp)en+1,s+1=pEn+1,s+1-pPn+1,s+1L(P,E)+∇pn+1,s·(n⃗-ξ⃗)e.
The terms labeled “n+1,s” become zero when ξ⃗=n⃗ is required. Repeating steps similar to (3.16)–(3.21) for other faces of control volume and substitute result into (3.15), one generates the equation for finding the pressure at next iteration (n+1,s+1) as
(3.22)1ΔΩ∑c=e,s,w,n(u⃗*·n⃗)cSc-ΔtΔΩ∑c=e,s,w,n(∇pn+1,s)c(n⃗-ξ⃗)c=ΔtΔΩ{(pE-pPL(P,E))n+1,s+1-(pP-pWL(P,W))n+1,s+1+(pNn+1,s-pPn+1,s+1L(P,N))-(pP-pSL(P,S))n+1,s+1}.
We use pNn+1,s instead of pNn+1,s+1 to make matrix of algebraic system to be tridiagonal.


### 3.3. Implementation of Boundary Conditions

The Finite Volume Method requires the boundary fluxes for each control volume to be either known or expressed through known quantities and interior nodal values. If the variables values are known at some boundary point, then there is no need to solve problem for it. A difficulty arises when approximations of normal derivatives are needed. Usually (see, e.g., [23]) these derivatives are approximated with lower order than the approximations used for interior point and may be one-sided differences. The accuracy of the results depended not only on the approximation near boundary but also on the accuracy of approximations at interior points. If higher accuracy is required, one has to use higher-order one-sided finite differences of derivatives at boundary and higher-order approximations at interior point. We used first-order one-sided finite differences near boundary.Impermeable Wall
The following condition is prescribed on the impermeable wall:(3.23)u⃗=u⃗wall.
This condition follows from the fact that a viscous fluid sticks to a solid wall. Since there is no flow through the wall, mass fluxes and convective fluxes of all quantities are zero. Diffusive fluxes in the momentum equation are approximated using known boundary values of the unknown and one-sided finite difference approximation for the gradients.Flow-Through Part
The implementations of three kinds of boundary conditions on the flow-through parts are addressed here. Only the case where the east face of the control volume aligns with flow-through boundaryΓ1 will be considered. A sketch of the grid and the notations used are shown in Figure 2(b). Other faces are treated similar.(a)
The velocity is set up (see (2.6) as
(3.24)u⃗Γ1=u⃗s1(x⃗,t).
Since the velocity vector is given, the mass flow rate and the convective fluxes can be calculated directly. The diffusive fluxes are not known, but they are approximated using known boundary values of the unknowns and one-sided finite difference approximation for the gradient. It is important to note how boundary condition (3.24) is involved in the derivation of the discrete pressure equation. Because (u⃗n+1)e is given by (3.24), the approximation of (3.13) becomes
(3.25)1ΔΩ[(u⃗n+1·n⃗)e+∑c=s,w,n(u⃗*·n⃗)cSc]-ΔtΔΩ∑c=s,w,n(Dnpn+1)cSc=0.
One does not need to approximate (Dnpn+1)e at face “e”. However, if pressure at the boundary Γ1 is needed at some stage, it can be obtained by extrapolation within the domain.(b)
The tangential velocity and pressure are prescribed (see, (2.5) as
(3.26)(u⃗·τ⃗)Γ1=G(x,y,t),pΓ1=H(x,y,t).
When the tangential velocity and pressure are specified on the flow-through part of boundary, the mass and convective fluxes are not known. One has to find them during the solution process. The solenoidal constraint ∇·u⃗=0 has to be applied at the boundary where the pressure is specified. Because the flow-through boundaries may not be aligned with the Cartesian coordinates, we will refer to the local coordinate system (n,τ) which is a rotated Cartesian frame with n in the direction of normal vector to the flow-through boundary and τ in the direction of the tangential vector to the flow-through boundary. The velocity vector u⃗=(ux,uy) can be expressed in terms of velocity components in local orthogonal coordinates u⃗=(Un,Uτ), where Un=u⃗·n⃗ is the normal velocity component to the flow-through boundary, and Uτ=u⃗·τ⃗ is the tangential velocity component to the flow-through boundary which is known at Γ1 from boundary condition (3.26). The continuity equation in terms of local orthogonal coordinates (n,τ) reads
(3.27)∂Un∂n+∂Uτ∂τ=0.
Using (3.26) and (3.27) yields
(3.28)(∂Un∂n)Γ1=-∂G∂τ.
To find the flux on the flow-through part, one needs to calculate the normal velocity (Un)e at the east cell face “e” (See Figure 2(b)). The normal derivative of Un at the east cell face is approximated by one-side difference:
(3.29)(∂Un∂n)e′=(Un)e′-(Un)PL(e′,P),
where e′ is the point of intersection of the line passing through node P parallel to normal vector to Γ1 at point “e” and the line coincide with boundary Γ1 (see Figure 2(b)). Following (3.28) and (3.29), the normal velocity component at point e′ is approximated as
(3.30)(u⃗n+1·n⃗)e′=(Unn+1)e′=(Unn+1)P-L(e′,P)(∂G∂τ)e′.
The discrete pressure equation for control volume Ω near flow-through boundary has the following form:
(3.31)1ΔΩ[(u⃗n+1·n⃗)e′Se+∑c=s,w,n(u⃗*·n⃗)cSc]-ΔtΔΩ∑c=s,w,n(Dnpn+1)cSc=0.
Here, the point “e′” is used instead of “e” to approximate the flux through the east face. In this case the order of approximation is reduced to first order. Moreover, in many cases, the grid is arranged such that “e′” coincides with the center of the east face. Substituting (3.30) into (3.31) and utilizing (3.13) at node P yields(3.32)1ΔΩ[(u⃗*·n⃗)P-Δt(∇pn+1,s+1·n⃗)P-L(e′,P)∂G∂τe′]Se+1ΔΩ[∑c=s,w,n(u⃗*·n⃗)cSc]-ΔtΔΩ∑c=s,w,n(Dnpn+1,s+1)cSc=0.
The derivative of pressure with respect to outward normal direction n at node P approximated by one-side difference is
(3.33)(Dnp)Pn+1,s+1=pe′n+1-pPn+1,s+1L(P,e′),
where L(P,e′) is the distance between nodes P and e′ on the boundary Γ1.(c)
The tangential velocity and total pressure are prescribed (see, (2.4) by
(3.34)(u⃗·τ⃗)Γ1=G(x,y,t),p+12|u⃗·u⃗|=H(x,y,t),(x,y)∈Γ1.
When the tangential velocity and total pressure are specified on the flow-through part, the situation arises where mass flux, convective flux, and pressure are not known. Let us use a local coordinates system (n,τ) as in the previous case. The flux (Un)e′=(u⃗·n⃗)e′ is approximated by (3.30). Since the pressure term on the flow-through boundary Γ1 (see Figure 2(b)) is unknown, one needs to approximate the pressure on the flow-through part by using the total pressure boundary condition, and one needs to calculate the pressure at point e′. The total pressure on flow-through part can be expressed in terms of local orthogonal coordinates (n,τ) in 2D at point e′ as
(3.35)pe′n+1,s+1+12|U⃗e′n+1|2=pe′n+1,s+1+12((Unn+1)e′2+(Uτ)e′2)=H.
Using boundary condition (3.34) the last equation recasts as
(3.36)pe′n+1,s+1+12(Unn+1)e′2=H-12Ge′2.
Substituting (Unn+1)e′ given by (3.30) yields
(3.37)pe′n+1,s+1+12((Unn+1)P-L(e′,P)2∂G∂τ)2=H-12Ge′2.
Using (Unn+1)P=(u⃗n+1·n⃗)P=(u⃗*·n⃗)P-Δt(∇pn+1,s+1·n⃗)P yields
(3.38)pe′n+1,s+1+12[(u⃗*·n⃗)P-Δt(∇pn+1,s+1·n⃗)P-L(e′,P)∂G∂τ]2=H-12Ge′.
Dropping terms of order O(Δt), one gets
(3.39)pe′n+1,s+1=H-12Ge′-12(u⃗*·n⃗)P2-12L(e′,P)2(∂G∂τ)e′2+(u⃗*·n⃗)PL(e′,P)(∂G∂τ)e′.
We have the previous case where pressure is given on the flow-through parts. When on the flow-through boundary n⃗=ξ⃗ and G=0, the expression for pe (3.38) reads
(3.40)pen+1,s+1=H-(u⃗*·n⃗)P2.


## 3.1. Time Discretization

The time discretization used here is based upon the simplest projection scheme originally proposed by Chorin and Temam (see, e.g., [19, 20]). This scheme has an irreducible splitting error of order O(Δt). Hence, using a higher-order time stepping scheme for the operator ∂/∂t-νΔ does not improve overall accuracy. Using the explicit Euler time stepping, the marching steps in time are the following.Setu⃗|t=0=u⃗0, then for n≥0 compute u⃗*,u⃗n+1, and pn+1 by solving first substep:
(3.1)u⃗*-u⃗nΔt+(u⃗n·∇)u⃗n=ν△u⃗n,
and second substep:
(3.2)u⃗n+1-u⃗*Δt=-∇pn+1,(3.3)∇·u⃗n+1=0,(u⃗n+1)Γ0=0,
where Δt=T/N is the time step, N is the integer, u⃗n≈u⃗(x⃗,nΔt), and pn+1≈p(x⃗,(n+1)Δt). Without loss of generality, density is equal to one, ρ=1. The explicit approximation of convective and viscous terms in (3.1) introduces restriction on the time step for stability. This is analyzed by many (see, e.g., [20, 21] and therein references).


## 3.2. Space Discretization

For the sake of simplicity and without loosing generality, the formulation of the numerical algorithm is illustrated for a two-dimensional domain. Letu⃗=(ux,uy) be the velocity vector, where ux and uy are the Cartesian components in x and y direction, respectively. The finite volume discretization is represented for nonorthogonal quadrilaterals grid. The collocated variable arrangement is utilized. Each discrete unknown is associated with the center of control volume Ω. First, we discretize the convection and diffusion parts of the Navier-Stokes equation. One can recast (3.1) in the form
(3.4)ϕ*-ϕnΔt+∇·(ϕnu⃗n)=νΔϕn,
where the variable ϕ can be either ux or uy, and u⃗n is such that ∇·u⃗n=0.The discrete form of (3.4) is obtained by integrating on each control volume Ω, followed by the application of the Gauss theorem:
(3.5)∫Ωϕ*-ϕnΔtdΩ+∮Sϕn(u⃗n·n⃗)dS=ν∮S∇ϕn·n⃗dS,
where S is the boundary of control volume Ω (e.g., in the case shown in Figure 2, S is the union of the control volume faces s,e,n, and w), and n⃗ is the unit outward normal vector to S. Using the midpoint rule to approximation, the surface and volume integrals yield
(3.6)∫Ωϕ*-ϕnΔtdΩ≈(ϕ*-ϕnΔt)PΔΩ,(3.7)∮Sϕn(u⃗n·n⃗)dS≈∑c=e,s,n,wϕcn(u⃗n·n⃗)cSc,(3.8)∮S∇ϕn·n⃗dS=∮SDnϕndS≈∑c=e,s,n,w(Dnϕn)cSc,
where ΔΩ is the volume of control volume Ω, Sc is the area of the “c” control volume face, and (Dnϕ)c is the derivative of Cartesian velocity components in the normal direction at the center of the “c” control volume face. To estimate the right-hand side in (3.7) and (3.8), we need to know the value of Cartesian velocity components and its normal derivative on the faces of each control volume. The implementation of Cartesian velocity components on nonorthogonal grids requires special attention because the boundary of the control volume is usually not aligned with the Cartesian velocity components. The 2D interpolation of irregularly-spaced data (see, e.g., [22]) is used to interpolate Cartesian velocity components on the boundary of each control volume in (3.7). Only the east side of a 2D control volume shown in Figure 2(a) will be considered. The same approach applies to other faces, only the indices need to be changed. For example, let ϕk be the value of Cartesian velocity components at point k where k=N,P,S,SE,E,NE, and let L(e,k) be the Cartesian distance between e and k. Using 2D interpolation yields
(3.9)ϕe=(∑kL(e,k)-2ϕk)(∑kL(e,k)-2),k=N,P,S,SE,E,NE,
where L(e,k)-2=1/[(xe-xk)2+(ye-yk)2]. The derivative of Cartesian velocity components in the normal direction at the center of the control volume face in (3.8) can be calculated by using the central difference approximation (see Figure 2(a)):
(3.10)(Dnϕ)e≈ϕE′-ϕP′L(P′,E′).
The auxiliary nodes P′ and E′ lie at the intersection of the line passing through the point “e" in the direction of normal vector n⃗ and the straight lines which connect nodes P and N or E and NE, respectively, and L(P′,E′) stands for the distance between P′ and E′. The values of ϕE′ and ϕP′ can be calculated by using the gradient at control volume center:
(3.11)ϕE′=ϕE+∇ϕE·(x⃗E′-x⃗E),ϕP′=ϕP+∇ϕP·(x⃗P′-x⃗P),
where x⃗P, x⃗E, x⃗P′, and x⃗E′ are the radius vectors of P, E, P′, and E′, respectively. The kth Cartesian components of ∇ϕP are approximated using Gauss’s theorem:
(3.12)∇ϕP·i⃗k=(∂ϕ∂xk)P=1ΔΩ∑c=e,s,n,wϕcn+1Sck,Sck=Sc(n⃗·ik⃗),
where Sc is the area of “c” control volume face, n⃗ is the unit outward normal vector to Sc, and i⃗k is the unit basis vector of Cartesian coordinate system (x1,x2)=(x,y). Using (3.6)–(3.12) to approximate (3.5), one can determine velocity field u⃗* (which is not solenoidal) at each grid node, even on the boundary.A typical 2D control volume and the notation used. A way of calculating cell face values and gradients.
(a)(b)In the first substep the continuity (3.3) is not used so that the intermediate velocity field is, in general, nondivergence free. The details of the setting and discretization of the second substep developed on nonuniform, collocated grid are discussed below. Equation (3.2) applies both in continuous and discrete sense. Taking the divergence of both sides of (3.2) and integrating over a control volume Ω, after applying the Gauss theorem and setting the update velocity filed, u⃗n+1, to be divergence free, one gets the equation
(3.13)0=1ΔΩ∮Su⃗n+1·n⃗dS=1ΔΩ∮Su⃗*·n⃗dS-Δt1ΔΩ∮S∇pn+1·n⃗dS,
that has to be discretized while collocating the variables in the control volume centers. Here n⃗ is outward normal to the boundary S of control volume Ω. At this stage of the projection procedure, the discrete values of ux* and uy* are already known and represent the source term in (3.13). A second-order discretization of the surface integrals can be obtained by utilizing the mean value formula. This means that the surface integrals in (3.13) can be approximated as
(3.14)1ΔΩ∮Su⃗n+1·n⃗dS≅1ΔΩ∑c=e,s,w,n(u⃗n+1·n⃗)cSc,1ΔΩ∮S∇pn+1·n⃗dS≅1ΔΩ∑c=e,s,w,n(∇pn+1·n⃗)cSc.
It follows that by substituting (3.14) into (3.13), one gets the discrete pressure equation
(3.15)1ΔΩ∑c=e,s,w,n(u⃗*·n⃗)cSc-ΔtΔΩ∑c=e,s,w,n(Dnpn+1)cSc=0.
The iterative method is utilized to approximate (Dnpn+1)c and solve (3.15). The normal-to-face intermediate velocities (u⃗*·n⃗)c,c=e,s,w,n are not directly available. They are found using interpolation. The derivative of pressure with respect to the direction of the outward normal n⃗ through the cell face “c”, (Dnp)cn+1 is approximated by on iterative technique (see, e.g., [23]) to reach a higher order of approximation and preserved compact stencil in the discrete equation (3.15). Only the east face of a 2D control volume shown in Figure 2(a) will be considered. The same approach applies to other faces. Using second upper index “s” to denote the number of iteration, one writes
(3.16)(Dnp)en+1,s+1=(Dξp)en+1,s+1+[(Dnp)e-(Dξp)e]n+1,s,s=0,…,Ŝ,(Dnp)n+1,0=(Dnp)n,
where ξ is the direction along the line connecting nodes P and E (see Figure 2(a)). The terms in the square brackets are approximated with high order and are evaluated by using values known from the previous iteration. Once the iterations converge, the low-order approximation term (Dξp)en+1,s+1 drops out, and the solution obtained corresponds to the higher order of approximation. The derivatives of pressure in the square brackets are written as
(3.17)(Dnp)en+1,s=(∇p·n⃗)en+1,s,(Dξp)en+1,s=(∇p·ξ⃗)en+1,s,
where n⃗ is the unit outward normal vector to cell face “e”, and ξ⃗ is the unit vector in ξ direction from point P to E. The term (∇p)en+1,s is approximated similar to (3.9) as
(3.18)(∇p)en+1,s=(∑lL(e,l)-2∇pln+1,s)(∑lL(e,l)-2),l=N,P,S,SE,E,NE,
where ∇pln+1,s is the gradient of the pressure at grid node l and L(e,l)-2=1/[(xe-xl)2+(ye-yl)2]. The kth components of ∇pln+1,s are discretized by using Gauss theorem (e.g., at grid node P):
(3.19)∇pPn+1,s·i⃗k=(∂pn+1,s∂xk)P≅1ΔΩ∑c=e,s,n,wpcn+1,sSck,Sck=Sc(n⃗·ik⃗).
The first term in the right-hand side of (3.16) is treated implicitly, and a simple approximation is used (that gives a compact stencil):
(3.20)(Dξp)en+1,s+1≈pEn+1,s+1-pPn+1,s+1L(P,E),
where L(P,E) is the distance between nodes P and E. The final expression for the approximation of the derivative of pressure with respect to n⃗ through the cell face “e” can now be written as
(3.21)(Dnp)en+1,s+1=pEn+1,s+1-pPn+1,s+1L(P,E)+∇pn+1,s·(n⃗-ξ⃗)e.
The terms labeled “n+1,s” become zero when ξ⃗=n⃗ is required. Repeating steps similar to (3.16)–(3.21) for other faces of control volume and substitute result into (3.15), one generates the equation for finding the pressure at next iteration (n+1,s+1) as
(3.22)1ΔΩ∑c=e,s,w,n(u⃗*·n⃗)cSc-ΔtΔΩ∑c=e,s,w,n(∇pn+1,s)c(n⃗-ξ⃗)c=ΔtΔΩ{(pE-pPL(P,E))n+1,s+1-(pP-pWL(P,W))n+1,s+1+(pNn+1,s-pPn+1,s+1L(P,N))-(pP-pSL(P,S))n+1,s+1}.
We use pNn+1,s instead of pNn+1,s+1 to make matrix of algebraic system to be tridiagonal.


## 3.3. Implementation of Boundary Conditions

The Finite Volume Method requires the boundary fluxes for each control volume to be either known or expressed through known quantities and interior nodal values. If the variables values are known at some boundary point, then there is no need to solve problem for it. A difficulty arises when approximations of normal derivatives are needed. Usually (see, e.g., [23]) these derivatives are approximated with lower order than the approximations used for interior point and may be one-sided differences. The accuracy of the results depended not only on the approximation near boundary but also on the accuracy of approximations at interior points. If higher accuracy is required, one has to use higher-order one-sided finite differences of derivatives at boundary and higher-order approximations at interior point. We used first-order one-sided finite differences near boundary.Impermeable Wall
The following condition is prescribed on the impermeable wall:(3.23)u⃗=u⃗wall.
This condition follows from the fact that a viscous fluid sticks to a solid wall. Since there is no flow through the wall, mass fluxes and convective fluxes of all quantities are zero. Diffusive fluxes in the momentum equation are approximated using known boundary values of the unknown and one-sided finite difference approximation for the gradients.Flow-Through Part
The implementations of three kinds of boundary conditions on the flow-through parts are addressed here. Only the case where the east face of the control volume aligns with flow-through boundaryΓ1 will be considered. A sketch of the grid and the notations used are shown in Figure 2(b). Other faces are treated similar.(a)
The velocity is set up (see (2.6) as
(3.24)u⃗Γ1=u⃗s1(x⃗,t).
Since the velocity vector is given, the mass flow rate and the convective fluxes can be calculated directly. The diffusive fluxes are not known, but they are approximated using known boundary values of the unknowns and one-sided finite difference approximation for the gradient. It is important to note how boundary condition (3.24) is involved in the derivation of the discrete pressure equation. Because (u⃗n+1)e is given by (3.24), the approximation of (3.13) becomes
(3.25)1ΔΩ[(u⃗n+1·n⃗)e+∑c=s,w,n(u⃗*·n⃗)cSc]-ΔtΔΩ∑c=s,w,n(Dnpn+1)cSc=0.
One does not need to approximate (Dnpn+1)e at face “e”. However, if pressure at the boundary Γ1 is needed at some stage, it can be obtained by extrapolation within the domain.(b)
The tangential velocity and pressure are prescribed (see, (2.5) as
(3.26)(u⃗·τ⃗)Γ1=G(x,y,t),pΓ1=H(x,y,t).
When the tangential velocity and pressure are specified on the flow-through part of boundary, the mass and convective fluxes are not known. One has to find them during the solution process. The solenoidal constraint ∇·u⃗=0 has to be applied at the boundary where the pressure is specified. Because the flow-through boundaries may not be aligned with the Cartesian coordinates, we will refer to the local coordinate system (n,τ) which is a rotated Cartesian frame with n in the direction of normal vector to the flow-through boundary and τ in the direction of the tangential vector to the flow-through boundary. The velocity vector u⃗=(ux,uy) can be expressed in terms of velocity components in local orthogonal coordinates u⃗=(Un,Uτ), where Un=u⃗·n⃗ is the normal velocity component to the flow-through boundary, and Uτ=u⃗·τ⃗ is the tangential velocity component to the flow-through boundary which is known at Γ1 from boundary condition (3.26). The continuity equation in terms of local orthogonal coordinates (n,τ) reads
(3.27)∂Un∂n+∂Uτ∂τ=0.
Using (3.26) and (3.27) yields
(3.28)(∂Un∂n)Γ1=-∂G∂τ.
To find the flux on the flow-through part, one needs to calculate the normal velocity (Un)e at the east cell face “e” (See Figure 2(b)). The normal derivative of Un at the east cell face is approximated by one-side difference:
(3.29)(∂Un∂n)e′=(Un)e′-(Un)PL(e′,P),
where e′ is the point of intersection of the line passing through node P parallel to normal vector to Γ1 at point “e” and the line coincide with boundary Γ1 (see Figure 2(b)). Following (3.28) and (3.29), the normal velocity component at point e′ is approximated as
(3.30)(u⃗n+1·n⃗)e′=(Unn+1)e′=(Unn+1)P-L(e′,P)(∂G∂τ)e′.
The discrete pressure equation for control volume Ω near flow-through boundary has the following form:
(3.31)1ΔΩ[(u⃗n+1·n⃗)e′Se+∑c=s,w,n(u⃗*·n⃗)cSc]-ΔtΔΩ∑c=s,w,n(Dnpn+1)cSc=0.
Here, the point “e′” is used instead of “e” to approximate the flux through the east face. In this case the order of approximation is reduced to first order. Moreover, in many cases, the grid is arranged such that “e′” coincides with the center of the east face. Substituting (3.30) into (3.31) and utilizing (3.13) at node P yields(3.32)1ΔΩ[(u⃗*·n⃗)P-Δt(∇pn+1,s+1·n⃗)P-L(e′,P)∂G∂τe′]Se+1ΔΩ[∑c=s,w,n(u⃗*·n⃗)cSc]-ΔtΔΩ∑c=s,w,n(Dnpn+1,s+1)cSc=0.
The derivative of pressure with respect to outward normal direction n at node P approximated by one-side difference is
(3.33)(Dnp)Pn+1,s+1=pe′n+1-pPn+1,s+1L(P,e′),
where L(P,e′) is the distance between nodes P and e′ on the boundary Γ1.(c)
The tangential velocity and total pressure are prescribed (see, (2.4) by
(3.34)(u⃗·τ⃗)Γ1=G(x,y,t),p+12|u⃗·u⃗|=H(x,y,t),(x,y)∈Γ1.
When the tangential velocity and total pressure are specified on the flow-through part, the situation arises where mass flux, convective flux, and pressure are not known. Let us use a local coordinates system (n,τ) as in the previous case. The flux (Un)e′=(u⃗·n⃗)e′ is approximated by (3.30). Since the pressure term on the flow-through boundary Γ1 (see Figure 2(b)) is unknown, one needs to approximate the pressure on the flow-through part by using the total pressure boundary condition, and one needs to calculate the pressure at point e′. The total pressure on flow-through part can be expressed in terms of local orthogonal coordinates (n,τ) in 2D at point e′ as
(3.35)pe′n+1,s+1+12|U⃗e′n+1|2=pe′n+1,s+1+12((Unn+1)e′2+(Uτ)e′2)=H.
Using boundary condition (3.34) the last equation recasts as
(3.36)pe′n+1,s+1+12(Unn+1)e′2=H-12Ge′2.
Substituting (Unn+1)e′ given by (3.30) yields
(3.37)pe′n+1,s+1+12((Unn+1)P-L(e′,P)2∂G∂τ)2=H-12Ge′2.
Using (Unn+1)P=(u⃗n+1·n⃗)P=(u⃗*·n⃗)P-Δt(∇pn+1,s+1·n⃗)P yields
(3.38)pe′n+1,s+1+12[(u⃗*·n⃗)P-Δt(∇pn+1,s+1·n⃗)P-L(e′,P)∂G∂τ]2=H-12Ge′.
Dropping terms of order O(Δt), one gets
(3.39)pe′n+1,s+1=H-12Ge′-12(u⃗*·n⃗)P2-12L(e′,P)2(∂G∂τ)e′2+(u⃗*·n⃗)PL(e′,P)(∂G∂τ)e′.
We have the previous case where pressure is given on the flow-through parts. When on the flow-through boundary n⃗=ξ⃗ and G=0, the expression for pe (3.38) reads
(3.40)pen+1,s+1=H-(u⃗*·n⃗)P2.


## 4. Results and Discussion

The proposed method is applied to test problems. The details of each of the problems and computed results are discussed in the following sections.


### 4.1. Flow between Two Parallel Plates

The purpose of this test is to estimate the potential and quality of the developed method in the case of unsteady flow. Considering the2D channel flow between two parallel plates, the Cartesian coordinate system (x,y,z) is chosen so that the x-axis is taken as the direction of flow, y is the coordinate normal to the plate, and z is the coordinate normal to x and y, respectively. The velocity field is assumed to be of the form u⃗=u(y,t)i⃗, where u is the velocity in the x-coordinate direction, and i⃗ is the unit vector in the x-coordinate direction. The Navier-Stokes equation implies that the pressure gradient is a function of time only, ∂p/∂x=f(t).Initial data att=0 is the fluid at the rest, u(y,0)=0. The flow is driven by pressure difference p2(L,t)-p1(0,t)=Δpcos(ωt) where L is the distance between the flow-through parts, ω is the frequency, and Δp is the characteristic pressure difference between two flow-through parts. The problem is dimensionalized with the height of the channel h as the length scale, Δp·h/L as the pressure scale, Δp·h/ρL as the velocity scale, and ρhL/Δp as the time scale. Nondimensional frequency is η=ωΔp/ρLh. Since the flow is driven by pressure difference and there is no velocity scale in the problem, we use ρU2=Δp·h/L in the traditional definition of the Reynolds number and call it the “Pressure Reynolds Number.”
(4.1)ReΔp=Uhν=hνΔphρL,
where ν is the kinematic viscosity. The analytical solution of the dimensionless problem obtained by separation of variables is
(4.2)u(y,t)=∑n=1∞[2(1-cos(nπ))nπsin(nπy)∫0te-λn2(t-τ)·cos(ητ)dτ],λn=nπReΔp1/2,0≤y≤1.
Computations are carried out with 1000 cells distributed in a uniform manner in the channel. A uniform grid having 20 lines across the channel and 50 lines in the direction of x was found to reproduce the flow parameters with good accuracy. In order to reduce computing cost, the distance between the flow-through parts was chosen to be one, L=1. The dependence between ReQ and ReΔp is plotted in Figure 3, for constant pressure drop p2(L,t)-p1(0,t)=Δp. The solid line represents the exact relation ReQ=ReΔp2/24, where ReQ=Q/2ν is the Reynolds number based on the flow rate, Q=∫01u(y)dy. Circle signs represent the results of our numerical simulations. The Reynolds number ReQ is not known a priori; it was computed at the end of the numerical simulation from the steady state flow rate obtained with the given ReΔp. As expected, the results are very close, and the velocity profile for all cases was the parabolic Poiseuille flow.Figure 3
The relation betweenReQ and ReΔp for η=0.From the analytical solution given by (4.2), it is obvious that the mass flow rate oscillation is a function of the oscillating frequency η and the pressure Reynolds number, ReΔp. In Figure 4, the variation of Q(t)=∫01u(y,t)dy with time is shown for given η=1 and 3, and ReΔp=150. Solid and dashed lines represent exact solutions for η=1 and 3, respectively. Circle and triangle signs correspond to the result of our numerical simulations for η=1and3, respectively. The numerical solution starts at t=0, and the time step is Δt=10-4. The above result corroborates that the proposed numerical method successfully predicts the volume rate for the constant and oscillated pressure drop.Figure 4
Volume rate versus time.


### 4.2. Flow with Circular Streamline

Another simple type of fluid motion through a bounded domain is one in which all the streamlines are circles centered on a common axis of symmetry. Steady motion can be generated by a circumferential pressure gradient in the domain between two concentric cylinders of radiir1 and r2. If the motion is to remain purely rotatory with the axial component of velocity to be zero, the axial pressure gradient must be zero, and the Navier-Stokes equations show that motion must be 2D. Using the equation of motion in polar coordinates (r,θ) and assuming that the velocity component in direction of the θ-coordinate line v=v(r) is a function of r only, and the radial velocity component is zero, one finds
(4.3)v2r=∂p∂r,∂2v∂r2+1r∂v∂r-vr2=1r∂p∂θ.
The variables in (4.3) are made nondimensional with h=r2-r1 as length scale, ν/h as velocity scale, and ρν2/h2 as pressure scale. Let d0=(r1+r2)/(2h)=R0/h be the nondimensional radius of centerline. Figure 5 represents a sketch of the problem geometry and main notations. It is easy to see from (4.3) that pressure has to be a linear function of θ:
(4.4)p(r,θ)=f(r)+K·θ,K=∂p∂θ-const,f(r)=∫d0-1/2rv2(ξ)ξdξ.Figure 5
The sketch of problem domain, flow with circular streamlines.With the boundary condition(4.5)v(d0±12)=0,
one obtains solution of (4.3)–(4.5) in the following form:
(4.6)v(r)=K8{C1r+C2r+4rlnr},d0-12≤r≤d0+12,(4.7)p(r,θ)=∫d0-1/2rv2(ξ)ξdξ+K·θ,0≤θ≤θ̂,d0-12≤r≤d0+12,(4.8)C1=(2d0-1)2ln(d0-1/2)-(2d0+1)2ln(d0+1/2)2d0,C2=(4d02-1)28d0ln(d0+1/2d0-1/2).
The nondimensional volume rate of flow becomes
(4.9)Q=∫d0-1/2d0+1/2vdr=K8E,E=2C1d0+C2ln((d0+1/2)(d0-1/2))-4d0+2[(d0+12)2ln(d0+12)-(d0-12)2ln(d0-12)].
Problem (4.3)–(4.5) can be considered as an example of the flowing-through problem where pressure and the tangential component of the velocity vector are given on flow-through parts AB and DC. It is worth to note here that the distribution of pressure is not constant at the flow-through parts and that the numerical solution uses the Navier-Stokes equation in terms of Cartesian coordinates and Cartesian velocity components u⃗=(ux,uy) where ux=-v(r)sin(θ), uy=v(r)cos(θ), r2=x2+y2, and θ=tan-1(y/x). Using exact solution (4.6)–(4.9), one can formulate the flowing-through problem where total pressure p+(1/2)[ux2(r,θ)+uy2(r,θ)] and tangent velocity are known on flow-through parts. It is also possible to consider problems where, in flow-through parts, different kinds of boundary conditions apply. The test cases of flowing-through problems computed in this section are summarized in Table 1. In all cases, we use 0≤θ≤π. Nonorthogonal logically rectangular boundary-fitted grids were constructed as follows. The impermeable boundaries AD and CB are partitioned equally into M subintervals. The flowing-though parts AB and CD are divided into an equal number of N subintervals. To reach steady flow, we used marching in time until the solution no longer changes. The grid independence study has been carried out for several values of circumferential pressure gradient, K, and for four cases of the flowing-through problems. The influence of the grid size on the difference between the exact velocity (4.6) and the approximate velocity in the maximum norm is shown in Table 2, for K=500. The convergence rates for the two finest grids are compared to the next coarser grid (see values in the brackets). Upper indices “ext” and “app” reference the exact and approximate solutions, respectively. It can be clearly seen from these results that the rate of convergence is near two. For Case  1, Figure 6 shows the variation of the dimensionless x-component of the velocity vector along the line θ=π/2 with circumferential pressure gradient ∂p/∂θ=K. The value of the circumferential pressure gradient varies from K=250 to K=1000.Table 1
Test cases0≤θ≤π,u→=(ux,uy), and H(r)=p+v(r)2/2.
CaseFlow-through boundarySolid wallAB(θ=0)CD(θ=π)AD,CB(0≤θ≤π)1ux=0,p=f(r)ux=0,p=f(r)+Kπux=uy=02ux=0,uy=v(r)ux=0,p=f(r)+Kπux=uy=03ux=0,uy=v(r)ux=0,p+uy2/2=H(r)ux=uy=04ux=0,p=f(r)ux=0,p+uy2/2=H(r)ux=uy=0Table 2
Rate of convergence ofuy for test cases.
GridM×N∥uapp-uext∥Case  1Case  2Case  3Case  420×101.508E-11.220E-11.218E-11.425E-140×203.979E-2(1.92)3.026E-2(2.01)3.018E-2(2.01)3.753E-2(1.93)80×409.955E-3(1.99)8.291E-3(1.87)7.590E-3(1.99)9.739E-3(1.95)Figure 6
Velocity profiles at the verticle lineθ=π/2(component ux), for Case  1 and different values of K.Figure7 shows pressure distribution for Case  1 along the line θ=π/2 and K=500. In both figures the solid lines represent the exact solutions (4.6) and (4.7), and the circle signs represent the numerical results. The calculated velocity profile and pressure along the line θ=const for Cases  2–4 are also in excellent agreement with the exact solution.Figure 7
Pressure along the lineθ=π/2 for Case  1.


### 4.3. Flowing-Through Problem for U-Bend Channel

For further validation, two-dimensional U-bend channel flow simulations are conducted. The flow configuration and main notations are shown in Figure8. The channel has a curvature ratio δ=R/d, where R is the radius of curvature, and d is the width of channel. The lengths of the channel before and after the bend L are taken sufficiently large to assume that pressure at sections A1A1′ and A2A2′ can be considered as constant, and fluid enters or leaves the channel legs with laminar, fully developed velocity profiles. The developed finite volume method has been utilized to simulate steady flow. For obtaining steady-state solution, the time is considered as pseudotime, and equations are iterated until the solution converges to steady state. Three kinds of the flowing-through problem have been considered. In all cases, no-slip boundary condition holds at the impermeable parts Γ10 and Γ20.Figure 8
Schematic diagram of  U-bend channel.The three flowing-through problems are formulated as follows.(P1)
On flow-through partsΓ11 and Γ21, the tangent components of velocity vector and pressure are specified (see (2.5) by
(4.10)u⃗·τ⃗=ux=0,p=p1,(x⃗)∈Γ11,u⃗·τ⃗=ux=0,p=p2,(x⃗)∈Γ21,
where τ⃗ is tangent unit vector to Γ11 and Γ21.(P2)
On flow-through partΓ11, the tangent and normal components of velocity vector are given (see (2.6) by
(4.11)u⃗=(ux,uy)=(0,us1(x)),x⃗∈Γ11,
where us1(x) is the parabolic Poiseuille velocity profile.On flow-through partΓ21, the tangent component of velocity and pressure are specified (see (2.5) by
(4.12)u⃗·τ⃗=ux=0,p=p2,(x⃗)∈Γ11.(P3)
On flow-through partΓ11, the tangent component of velocity vector and total pressure are prescribed (see (2.4) by
(4.13)u⃗·τ⃗=ux=0,p+12ρ|u⃗·u⃗|=H1(x⃗),(x⃗)∈Γ11,
where H1(x⃗) is a given function and is computed from the solution of P2. On the flow-through parts Γ21, the tangent component of the velocity vector and pressure are known (see (2.5) by
(4.14)u⃗·τ⃗=ux=0,p=p2,(x⃗)∈Γ21.The main characteristic of flow in curve channels is pressure loss. The pressure losses are presented in the form of friction factor versus Reynolds number:(4.15)fw=2ΔpρU2(L/d)=f(Re),
where U is the mean velocity, ρ is the density of the fluid, Re=Ud/ν is the Reynolds number, ν is the kinematic viscosity of fluid, and Δp is the pressure losses, Δp=p2-p1. Before the main computations were started, a test was executed with a straight channel. A very good agreement of the computed pressure losses with the theoretical solution based on the Poiseuille law fw≈36/Re was observed. Based on the preliminary experiments, the length of the channel legs l=L/d=5 was used in the main computations represented below. The impermeable boundaries A1A2 and A1′A2′ were equally partitioned into M subintervals. The flowing-through parts A1A1′ and A2A2′ were divided into an equal number of N subintervals. Three grid sequences of 100×10, 200×20, and 400×40 nodes were tested. Computations using these grid sequences are shown in Table 3. In the case of the flowing-through problem P1, the pressure losses are known a priori, and the Reynolds number was computed from the steady state flow rate. In problem P2 the Reynolds number is known a priori, and Δp was estimated from the steady state flow regime. In the problem P3 neither Δp nor Re is known a priori, and both of them were computed at the end of the numerical simulation from steady state.Table 3
Friction factor for three kinds of the flowing-through problem.Re=100, δ=3.
GridM×NFriction factor,fwP1P2P3100×100.431280.4331530.431782200×200.4296730.4316390.429832400×400.4296470.4315490.429587Total pressure losses of a U-bend channel flow are presented in the form of the friction factor versus Reynolds numberfw=f(Re) in Figure 9, where the effect of the dimensionless curvature ratio, δ=R/d, is shown. All three flowing-through problems P1, P2, and P3 give very close results. From Figure 9 it is seen that the effect of the channel curvature ratio on the friction factor is small for δ>3 for all tested flowing-through problems. The friction factor fw increases with decreasing δ. In Figure 10, streamline patterns are presented. Figure 10(a) is drawn for δ=1 and Re=200, Figure 10(b) shows the case of δ=1 and Re=300, Figure 10(c) depicts δ=0.6 and Re=200, and Figure 10(d) is drawn for δ=0.6 and Re=300. The sharp bend δ=0.6 and increasing Reynolds number cause separation which occurs on the right side of the bend. The size of the separation zone increases with increasing flow rate and decreasing δ.Figure 9
Friction factor as a function  of the Reynolds number.Streamline patterns of flow in the U-bend channel for variousδ and Re.
(a)
δ=1.0, Re=200(b)
δ=1.0, Re=300(c)
δ=0.6, Re=200(d)
δ=0.6, Re=300The velocity profile in the cross sectiony=1 of the right-hand side leg of U-bend is depicted in Figure 11 for Re=200 and 300 and δ=0.6.Figure 11
Velocity profile in sectiony=1 of right-hand side of U-bend channel. δ=0.6, Re=200 and 300.


### 4.4. Flow in Planar T-Junction Channel

The T-junction flow geometry is schematically represented in Figure12. The origin of the coordinate system is located in the lower horizontal boundary opposite the left corner of branch as demonstrated. The left-hand side branch, the upper branch, the right-hand side branch, and the junction area are denoted by Γ1, Γ2, Γ3, and Γ4, respectively. All branches have the same width w.Figure 12
Schematic geometry of T-junction bifurcation and coordinate system.The flow rate ratio is defined asβ≡Q3/Q1 where Q1 and Q3 are the inlet duct and branch duct flow rates per unit span, respectively. The following problem has been considered.(i)
On flow-through partΓ11a laminar, fully developed, parabolic velocity profile is prescribed by
(4.16)u⃗=(ux(y),0),(x,y)∈Γ11.(ii)
On flow-through partsΓ21 and Γ31 the tangent component of the velocity vector and the pressure are specified by
(4.17)u⃗·τ⃗2=uy=0,p=p2,(x,y)∈Γ21,u⃗·τ⃗3=ux=0,p=p3,(x,y)∈Γ31,
where τ⃗iis the unit tangent vector to Γi1,i=2,3.The calculations are compared with those of Hayes et al. [13], Kelkar and Choudhury [16], and Fluent [24]. A flowing-through problem with ux=4y-4y2 and equal static pressure p2=p3=0 is considered. The Navier-Stokes equation dimensionalized with the width, w, as characteristic length, the inlet centerline velocity Uc as the characteristic velocity, and ρUc2 as the scale of pressure. A range of Reynolds number Re=wUc/ν, where ν is the kinematic viscosity, is studied with Re∈[10,400]. The computational domain is set to have lengths of L1/w=2 and L2/w=L3/w=3 according to the results represented in Fluent Inc. [24]. The square meshes containing 20, 30, and 40  cells from wall to wall are used. The studied cases start from a motionless state. A steady flow is achieved if the following condition is held: ∥u⃗n+1-u⃗n∥≤ε=10-8. The maximum norm of grid function is used. Figure 13 shows the effect of increasing the Reynolds number on the flow split between the main and the side exit branches. The value of β increases from 0.5 for a small Reynolds number, Re<10, to about 0.9 at Re=400. Figure 14 shows the predicted streamline pattern and pressure contour plots for two Reynolds Numbers Re=100,400. Flow separation from the left wall of the upper branch occurs at all considered Reynolds numbers. These are very similar to those reported in Fluent Inc. (1998). The size and extent of flow separation zone are in a good agreement with results of Hayes et al. [13], Kelkar and Choudhury [16], and Fluent [24].Figure 13
The flow rate ratio,β, as a function of Reynolds number, Re.Streamline patterns and pressure contour of flow in the T-junction for various Re and equal static pressure at the existsΓ21and Γ31.
(a)(b)(c)(d)


## 4.1. Flow between Two Parallel Plates

The purpose of this test is to estimate the potential and quality of the developed method in the case of unsteady flow. Considering the2D channel flow between two parallel plates, the Cartesian coordinate system (x,y,z) is chosen so that the x-axis is taken as the direction of flow, y is the coordinate normal to the plate, and z is the coordinate normal to x and y, respectively. The velocity field is assumed to be of the form u⃗=u(y,t)i⃗, where u is the velocity in the x-coordinate direction, and i⃗ is the unit vector in the x-coordinate direction. The Navier-Stokes equation implies that the pressure gradient is a function of time only, ∂p/∂x=f(t).Initial data att=0 is the fluid at the rest, u(y,0)=0. The flow is driven by pressure difference p2(L,t)-p1(0,t)=Δpcos(ωt) where L is the distance between the flow-through parts, ω is the frequency, and Δp is the characteristic pressure difference between two flow-through parts. The problem is dimensionalized with the height of the channel h as the length scale, Δp·h/L as the pressure scale, Δp·h/ρL as the velocity scale, and ρhL/Δp as the time scale. Nondimensional frequency is η=ωΔp/ρLh. Since the flow is driven by pressure difference and there is no velocity scale in the problem, we use ρU2=Δp·h/L in the traditional definition of the Reynolds number and call it the “Pressure Reynolds Number.”
(4.1)ReΔp=Uhν=hνΔphρL,
where ν is the kinematic viscosity. The analytical solution of the dimensionless problem obtained by separation of variables is
(4.2)u(y,t)=∑n=1∞[2(1-cos(nπ))nπsin(nπy)∫0te-λn2(t-τ)·cos(ητ)dτ],λn=nπReΔp1/2,0≤y≤1.
Computations are carried out with 1000 cells distributed in a uniform manner in the channel. A uniform grid having 20 lines across the channel and 50 lines in the direction of x was found to reproduce the flow parameters with good accuracy. In order to reduce computing cost, the distance between the flow-through parts was chosen to be one, L=1. The dependence between ReQ and ReΔp is plotted in Figure 3, for constant pressure drop p2(L,t)-p1(0,t)=Δp. The solid line represents the exact relation ReQ=ReΔp2/24, where ReQ=Q/2ν is the Reynolds number based on the flow rate, Q=∫01u(y)dy. Circle signs represent the results of our numerical simulations. The Reynolds number ReQ is not known a priori; it was computed at the end of the numerical simulation from the steady state flow rate obtained with the given ReΔp. As expected, the results are very close, and the velocity profile for all cases was the parabolic Poiseuille flow.Figure 3
The relation betweenReQ and ReΔp for η=0.From the analytical solution given by (4.2), it is obvious that the mass flow rate oscillation is a function of the oscillating frequency η and the pressure Reynolds number, ReΔp. In Figure 4, the variation of Q(t)=∫01u(y,t)dy with time is shown for given η=1 and 3, and ReΔp=150. Solid and dashed lines represent exact solutions for η=1 and 3, respectively. Circle and triangle signs correspond to the result of our numerical simulations for η=1and3, respectively. The numerical solution starts at t=0, and the time step is Δt=10-4. The above result corroborates that the proposed numerical method successfully predicts the volume rate for the constant and oscillated pressure drop.Figure 4
Volume rate versus time.


## 4.2. Flow with Circular Streamline

Another simple type of fluid motion through a bounded domain is one in which all the streamlines are circles centered on a common axis of symmetry. Steady motion can be generated by a circumferential pressure gradient in the domain between two concentric cylinders of radiir1 and r2. If the motion is to remain purely rotatory with the axial component of velocity to be zero, the axial pressure gradient must be zero, and the Navier-Stokes equations show that motion must be 2D. Using the equation of motion in polar coordinates (r,θ) and assuming that the velocity component in direction of the θ-coordinate line v=v(r) is a function of r only, and the radial velocity component is zero, one finds
(4.3)v2r=∂p∂r,∂2v∂r2+1r∂v∂r-vr2=1r∂p∂θ.
The variables in (4.3) are made nondimensional with h=r2-r1 as length scale, ν/h as velocity scale, and ρν2/h2 as pressure scale. Let d0=(r1+r2)/(2h)=R0/h be the nondimensional radius of centerline. Figure 5 represents a sketch of the problem geometry and main notations. It is easy to see from (4.3) that pressure has to be a linear function of θ:
(4.4)p(r,θ)=f(r)+K·θ,K=∂p∂θ-const,f(r)=∫d0-1/2rv2(ξ)ξdξ.Figure 5
The sketch of problem domain, flow with circular streamlines.With the boundary condition(4.5)v(d0±12)=0,
one obtains solution of (4.3)–(4.5) in the following form:
(4.6)v(r)=K8{C1r+C2r+4rlnr},d0-12≤r≤d0+12,(4.7)p(r,θ)=∫d0-1/2rv2(ξ)ξdξ+K·θ,0≤θ≤θ̂,d0-12≤r≤d0+12,(4.8)C1=(2d0-1)2ln(d0-1/2)-(2d0+1)2ln(d0+1/2)2d0,C2=(4d02-1)28d0ln(d0+1/2d0-1/2).
The nondimensional volume rate of flow becomes
(4.9)Q=∫d0-1/2d0+1/2vdr=K8E,E=2C1d0+C2ln((d0+1/2)(d0-1/2))-4d0+2[(d0+12)2ln(d0+12)-(d0-12)2ln(d0-12)].
Problem (4.3)–(4.5) can be considered as an example of the flowing-through problem where pressure and the tangential component of the velocity vector are given on flow-through parts AB and DC. It is worth to note here that the distribution of pressure is not constant at the flow-through parts and that the numerical solution uses the Navier-Stokes equation in terms of Cartesian coordinates and Cartesian velocity components u⃗=(ux,uy) where ux=-v(r)sin(θ), uy=v(r)cos(θ), r2=x2+y2, and θ=tan-1(y/x). Using exact solution (4.6)–(4.9), one can formulate the flowing-through problem where total pressure p+(1/2)[ux2(r,θ)+uy2(r,θ)] and tangent velocity are known on flow-through parts. It is also possible to consider problems where, in flow-through parts, different kinds of boundary conditions apply. The test cases of flowing-through problems computed in this section are summarized in Table 1. In all cases, we use 0≤θ≤π. Nonorthogonal logically rectangular boundary-fitted grids were constructed as follows. The impermeable boundaries AD and CB are partitioned equally into M subintervals. The flowing-though parts AB and CD are divided into an equal number of N subintervals. To reach steady flow, we used marching in time until the solution no longer changes. The grid independence study has been carried out for several values of circumferential pressure gradient, K, and for four cases of the flowing-through problems. The influence of the grid size on the difference between the exact velocity (4.6) and the approximate velocity in the maximum norm is shown in Table 2, for K=500. The convergence rates for the two finest grids are compared to the next coarser grid (see values in the brackets). Upper indices “ext” and “app” reference the exact and approximate solutions, respectively. It can be clearly seen from these results that the rate of convergence is near two. For Case  1, Figure 6 shows the variation of the dimensionless x-component of the velocity vector along the line θ=π/2 with circumferential pressure gradient ∂p/∂θ=K. The value of the circumferential pressure gradient varies from K=250 to K=1000.Table 1
Test cases0≤θ≤π,u→=(ux,uy), and H(r)=p+v(r)2/2.
CaseFlow-through boundarySolid wallAB(θ=0)CD(θ=π)AD,CB(0≤θ≤π)1ux=0,p=f(r)ux=0,p=f(r)+Kπux=uy=02ux=0,uy=v(r)ux=0,p=f(r)+Kπux=uy=03ux=0,uy=v(r)ux=0,p+uy2/2=H(r)ux=uy=04ux=0,p=f(r)ux=0,p+uy2/2=H(r)ux=uy=0Table 2
Rate of convergence ofuy for test cases.
GridM×N∥uapp-uext∥Case  1Case  2Case  3Case  420×101.508E-11.220E-11.218E-11.425E-140×203.979E-2(1.92)3.026E-2(2.01)3.018E-2(2.01)3.753E-2(1.93)80×409.955E-3(1.99)8.291E-3(1.87)7.590E-3(1.99)9.739E-3(1.95)Figure 6
Velocity profiles at the verticle lineθ=π/2(component ux), for Case  1 and different values of K.Figure7 shows pressure distribution for Case  1 along the line θ=π/2 and K=500. In both figures the solid lines represent the exact solutions (4.6) and (4.7), and the circle signs represent the numerical results. The calculated velocity profile and pressure along the line θ=const for Cases  2–4 are also in excellent agreement with the exact solution.Figure 7
Pressure along the lineθ=π/2 for Case  1.


## 4.3. Flowing-Through Problem for U-Bend Channel

For further validation, two-dimensional U-bend channel flow simulations are conducted. The flow configuration and main notations are shown in Figure8. The channel has a curvature ratio δ=R/d, where R is the radius of curvature, and d is the width of channel. The lengths of the channel before and after the bend L are taken sufficiently large to assume that pressure at sections A1A1′ and A2A2′ can be considered as constant, and fluid enters or leaves the channel legs with laminar, fully developed velocity profiles. The developed finite volume method has been utilized to simulate steady flow. For obtaining steady-state solution, the time is considered as pseudotime, and equations are iterated until the solution converges to steady state. Three kinds of the flowing-through problem have been considered. In all cases, no-slip boundary condition holds at the impermeable parts Γ10 and Γ20.Figure 8
Schematic diagram of  U-bend channel.The three flowing-through problems are formulated as follows.(P1)
On flow-through partsΓ11 and Γ21, the tangent components of velocity vector and pressure are specified (see (2.5) by
(4.10)u⃗·τ⃗=ux=0,p=p1,(x⃗)∈Γ11,u⃗·τ⃗=ux=0,p=p2,(x⃗)∈Γ21,
where τ⃗ is tangent unit vector to Γ11 and Γ21.(P2)
On flow-through partΓ11, the tangent and normal components of velocity vector are given (see (2.6) by
(4.11)u⃗=(ux,uy)=(0,us1(x)),x⃗∈Γ11,
where us1(x) is the parabolic Poiseuille velocity profile.On flow-through partΓ21, the tangent component of velocity and pressure are specified (see (2.5) by
(4.12)u⃗·τ⃗=ux=0,p=p2,(x⃗)∈Γ11.(P3)
On flow-through partΓ11, the tangent component of velocity vector and total pressure are prescribed (see (2.4) by
(4.13)u⃗·τ⃗=ux=0,p+12ρ|u⃗·u⃗|=H1(x⃗),(x⃗)∈Γ11,
where H1(x⃗) is a given function and is computed from the solution of P2. On the flow-through parts Γ21, the tangent component of the velocity vector and pressure are known (see (2.5) by
(4.14)u⃗·τ⃗=ux=0,p=p2,(x⃗)∈Γ21.The main characteristic of flow in curve channels is pressure loss. The pressure losses are presented in the form of friction factor versus Reynolds number:(4.15)fw=2ΔpρU2(L/d)=f(Re),
where U is the mean velocity, ρ is the density of the fluid, Re=Ud/ν is the Reynolds number, ν is the kinematic viscosity of fluid, and Δp is the pressure losses, Δp=p2-p1. Before the main computations were started, a test was executed with a straight channel. A very good agreement of the computed pressure losses with the theoretical solution based on the Poiseuille law fw≈36/Re was observed. Based on the preliminary experiments, the length of the channel legs l=L/d=5 was used in the main computations represented below. The impermeable boundaries A1A2 and A1′A2′ were equally partitioned into M subintervals. The flowing-through parts A1A1′ and A2A2′ were divided into an equal number of N subintervals. Three grid sequences of 100×10, 200×20, and 400×40 nodes were tested. Computations using these grid sequences are shown in Table 3. In the case of the flowing-through problem P1, the pressure losses are known a priori, and the Reynolds number was computed from the steady state flow rate. In problem P2 the Reynolds number is known a priori, and Δp was estimated from the steady state flow regime. In the problem P3 neither Δp nor Re is known a priori, and both of them were computed at the end of the numerical simulation from steady state.Table 3
Friction factor for three kinds of the flowing-through problem.Re=100, δ=3.
GridM×NFriction factor,fwP1P2P3100×100.431280.4331530.431782200×200.4296730.4316390.429832400×400.4296470.4315490.429587Total pressure losses of a U-bend channel flow are presented in the form of the friction factor versus Reynolds numberfw=f(Re) in Figure 9, where the effect of the dimensionless curvature ratio, δ=R/d, is shown. All three flowing-through problems P1, P2, and P3 give very close results. From Figure 9 it is seen that the effect of the channel curvature ratio on the friction factor is small for δ>3 for all tested flowing-through problems. The friction factor fw increases with decreasing δ. In Figure 10, streamline patterns are presented. Figure 10(a) is drawn for δ=1 and Re=200, Figure 10(b) shows the case of δ=1 and Re=300, Figure 10(c) depicts δ=0.6 and Re=200, and Figure 10(d) is drawn for δ=0.6 and Re=300. The sharp bend δ=0.6 and increasing Reynolds number cause separation which occurs on the right side of the bend. The size of the separation zone increases with increasing flow rate and decreasing δ.Figure 9
Friction factor as a function  of the Reynolds number.Streamline patterns of flow in the U-bend channel for variousδ and Re.
(a)
δ=1.0, Re=200(b)
δ=1.0, Re=300(c)
δ=0.6, Re=200(d)
δ=0.6, Re=300The velocity profile in the cross sectiony=1 of the right-hand side leg of U-bend is depicted in Figure 11 for Re=200 and 300 and δ=0.6.Figure 11
Velocity profile in sectiony=1 of right-hand side of U-bend channel. δ=0.6, Re=200 and 300.


## 4.4. Flow in Planar T-Junction Channel

The T-junction flow geometry is schematically represented in Figure12. The origin of the coordinate system is located in the lower horizontal boundary opposite the left corner of branch as demonstrated. The left-hand side branch, the upper branch, the right-hand side branch, and the junction area are denoted by Γ1, Γ2, Γ3, and Γ4, respectively. All branches have the same width w.Figure 12
Schematic geometry of T-junction bifurcation and coordinate system.The flow rate ratio is defined asβ≡Q3/Q1 where Q1 and Q3 are the inlet duct and branch duct flow rates per unit span, respectively. The following problem has been considered.(i)
On flow-through partΓ11a laminar, fully developed, parabolic velocity profile is prescribed by
(4.16)u⃗=(ux(y),0),(x,y)∈Γ11.(ii)
On flow-through partsΓ21 and Γ31 the tangent component of the velocity vector and the pressure are specified by
(4.17)u⃗·τ⃗2=uy=0,p=p2,(x,y)∈Γ21,u⃗·τ⃗3=ux=0,p=p3,(x,y)∈Γ31,
where τ⃗iis the unit tangent vector to Γi1,i=2,3.The calculations are compared with those of Hayes et al. [13], Kelkar and Choudhury [16], and Fluent [24]. A flowing-through problem with ux=4y-4y2 and equal static pressure p2=p3=0 is considered. The Navier-Stokes equation dimensionalized with the width, w, as characteristic length, the inlet centerline velocity Uc as the characteristic velocity, and ρUc2 as the scale of pressure. A range of Reynolds number Re=wUc/ν, where ν is the kinematic viscosity, is studied with Re∈[10,400]. The computational domain is set to have lengths of L1/w=2 and L2/w=L3/w=3 according to the results represented in Fluent Inc. [24]. The square meshes containing 20, 30, and 40  cells from wall to wall are used. The studied cases start from a motionless state. A steady flow is achieved if the following condition is held: ∥u⃗n+1-u⃗n∥≤ε=10-8. The maximum norm of grid function is used. Figure 13 shows the effect of increasing the Reynolds number on the flow split between the main and the side exit branches. The value of β increases from 0.5 for a small Reynolds number, Re<10, to about 0.9 at Re=400. Figure 14 shows the predicted streamline pattern and pressure contour plots for two Reynolds Numbers Re=100,400. Flow separation from the left wall of the upper branch occurs at all considered Reynolds numbers. These are very similar to those reported in Fluent Inc. (1998). The size and extent of flow separation zone are in a good agreement with results of Hayes et al. [13], Kelkar and Choudhury [16], and Fluent [24].Figure 13
The flow rate ratio,β, as a function of Reynolds number, Re.Streamline patterns and pressure contour of flow in the T-junction for various Re and equal static pressure at the existsΓ21and Γ31.
(a)(b)(c)(d)


## 5. Conclusion

A mathematical formulation of well-posed initial boundary value problems for viscous incompressible fluid flow-through-bounded domain is described for the case where the values of static or total pressure and tangential components of the velocity vector on flow-through parts of the domain boundary are prescribed. A computational method for the approximate solution of these well-posed problems is developed within the framework of the finite volume approach. The robustness of the method is validated by its application for channel flows driven by pressure drop for which analytical solutions are available (2D Poiseuille flow, purely rotatory flow in the annular domain between cylinders). The effect of curvature ratio of planar U-bend channel is analyzed for various flowing-through problem formulations. The flow through planar T-junction channel is utilized as a benchmark test in the case of several flow-through parts of boundary. Results of all tests confirm the reliability and accuracy of developed method. The method is robust and accurate in simulating incompressible flows in domains with known boundary pressure (or total pressure) and with known velocity profiles in flow-through parts of boundary.



---

*Source: 372703-2009-06-15.xml*