# Numerical Performance of Higher-Order Semicompact Scheme for Arbitrary Triangular Cavity Flow

Authors: Xiaofeng Wang; Dongyang Shi
Journal: International Journal of Engineering Mathematics
 (2014)

---

## Abstract

An efficient fourth-order semicompact finite difference scheme has been
developed to solve steady incompressible Navier-Stokes (N-S) equations in stream function and vorticity formulation in a triangular cavity of arbitrary geometry. The governing equations are transformed into curvilinear coordinates by a simple linear transformation to handle the nonregular geometry of the problem. The main feature of the new higher-order semicompact scheme is that it can calculate a triangle flow with arbitrary shape for high Reynolds numbers. It is found that the solutions obtained with the present scheme are in good agreement with the analytical results or with the existing results depending on the availability.


---

## Body

## 1. Introduction

Numerical simulations for the solution of steady incompressible viscous flow within a driven cavity are a complex and significant topic, and a lot of researchers contributed to this subject [1–4]. Although there are still some minor discrepancies in the results, the square cavity flow has been essentially characterized [5]. As pointed out in McQuain et al. [6], the results for the square cavity may not be applied to other important geometries such as a triangular cavity. Also the latter shapes are more common in practice.The problem under consideration is the steady motion of an incompressible viscous flow in a triangular cavity of arbitrary geometry. This flow was studied numerically by McQuain et al. [6], Jyotsna and Vanka [7], Li and Tang [8], and Gaskell et al. [9]. McQuain et al. [6] applied the Batchelor’s mean square law to triangular cavity flow and analytically obtained the inviscid core vorticity for infinite Reynolds number. Recent calculations of the steady problem in an equilateral triangular cavity have been given in Ribbens et al. [10] for Re
≤
500. Erturk and Gokcol [11] have presented high accurate, fine grid solutions of 2D steady incompressible flow in triangle cavities. The governing equations are solved up to very low residuals at various Reynolds numbers. The results of these studies, however, show some discrepancies. Furthermore, the scheme in [11] is only second-order spatial accuracy. This constitutes the main motivation of this study.The main object of this study is the development of an accurate and efficient scheme for solving the N-S problems in triangular geometries. In the present paper, we have attempted an efficient fourth-order semicompact scheme (compactness of the scheme is relaxed for few terms of the governing equations) for stream function vorticity form of incompressible Navier-Stokes equations in curvilinear coordinates inside a triangular cavity. A geometric transformation handling triangles of arbitrary shape is also presented. The developed equations are solved numerically subject to appropriate boundary conditions by a fourth-order accurate finite difference method. Overall, besides opening up new possibilities, the method may be considered an efficient one for computation of flow for this physical configuration and the results thus obtained are not only useful for engineers in the process design; they also supplement the existing literature.


## 2. Governing Equations

A triangular driven cavity of general shape is given by locating its three vertices atA
^
(
x
a
,
h
), B
^
(
x
b
,
h
), and C
^
(
x
c
,
y
c
), with the upper side moving to the right via a constant velocity U
0; see Figure 1(a). The 2D steady incompressible flow inside a triangular cavity is governed by the N-S equations. We consider the N-S equations in streamfunction and vorticity formulation in Cartesian coordinates, such that [5]
(1)
∇
x
y
2
ψ
=
-
ω
,
∇
x
y
2
ω
=
Re
(
U
·
∇
x
y
ω
)
,
where
(2)
∇
x
y
2
=
∂
2
∂
x
2
+
∂
2
∂
y
2
,
∇
x
y
=
(
∂
∂
x
,
∂
∂
y
)
,
U
=
(
u
,
v
)
,
u
=
∂
ψ
∂
y
,
v
=
-
∂
ψ
∂
x
,
ω
=
∂
v
∂
x
-
∂
u
∂
y
,
where ψ is the streamfunction, ω is the vorticity, u and v are the components of the velocity in x- and y-directions, respectively, and Re is Reynolds number. We note that these equations are non dimensional and a length scale of (
h
-
c
y
)
/
3 and a velocity scale of U
0, that is, the velocity of the lid, are used to nondimensionalize the parameters and Reynolds number is defined accordingly.Geometric transformation of the triangular cavity: Physical domain (a) and Computational domain (b).
(a)
(b)By a simple linear transformation(3)
ξ
=
x
-
x
a
x
b
-
x
a
+
(
x
a
-
x
c
)
(
y
-
h
)
(
x
b
-
x
a
)
(
y
c
-
h
)
,
η
=
(
y
c
-
h
)
2
+
(
x
c
-
x
a
)
2
(
x
b
-
x
a
)
(
y
c
-
h
)
(
y
-
h
)
,
the triangle ∆
A
^
B
^
C
^ is transformed to a right-angled triangle with vertices A
(
0,0
), B
(
1,0
), and C
(
0
,
(
y
c
-
h
)
2
+
(
x
c
-
x
a
)
2
/
(
x
b
-
x
a
)
); see Figure 1(b). From these relations, we can calculate the transformation metrics as follows:
(4)
∂
ξ
∂
x
=
1
x
b
-
x
a
,
∂
ξ
∂
y
=
x
a
-
x
c
(
x
b
-
x
a
)
(
y
c
-
h
)
,
∂
η
∂
x
=
0
,
∂
η
∂
y
=
(
y
c
-
h
)
2
+
(
x
c
-
x
a
)
2
(
x
b
-
x
a
)
(
y
c
-
h
)
,
and also
(5)
∂
2
ξ
∂
x
2
=
∂
2
ξ
∂
y
2
=
∂
2
η
∂
x
2
=
∂
2
η
∂
y
2
=
0
.
Using the chain rule, the governing N-S equations (1) in general curvilinear coordinates are as follows:
(6)
(
(
∂
ξ
∂
x
)
2
+
(
∂
ξ
∂
y
)
2
)
∂
2
ψ
∂
ξ
2
+
(
(
∂
η
∂
x
)
2
+
(
∂
η
∂
y
)
2
)
∂
2
ψ
∂
η
2
+
(
∂
2
ξ
∂
x
2
+
∂
2
ξ
∂
y
2
)
∂
ψ
∂
ξ
+
(
∂
2
η
∂
x
2
+
∂
2
η
∂
y
2
)
∂
ψ
∂
η
+
2
(
∂
ξ
∂
x
∂
η
∂
x
+
∂
ξ
∂
y
∂
η
∂
y
)
∂
2
ψ
∂
ξ
∂
η
=
-
ω
,
(
(
∂
ξ
∂
x
)
2
+
(
∂
ξ
∂
y
)
2
)
∂
2
ω
∂
ξ
2
+
(
(
∂
η
∂
x
)
2
+
(
∂
η
∂
y
)
2
)
∂
2
ω
∂
η
2
+
(
∂
2
ξ
∂
x
2
+
∂
2
ξ
∂
y
2
)
∂
ω
∂
ξ
+
(
∂
2
η
∂
x
2
+
∂
2
η
∂
y
2
)
∂
ω
∂
η
+
2
(
∂
ξ
∂
x
∂
η
∂
x
+
∂
ξ
∂
y
∂
η
∂
y
)
∂
2
ω
∂
ξ
∂
η
=
Ω
,
where
(7)
Ω
=
Re
∂
ξ
∂
x
∂
η
∂
y
(
∂
ψ
∂
η
∂
ω
∂
ξ
-
∂
ψ
∂
ξ
∂
ω
∂
η
)
.
Substituting for the transformation metrics, we obtain the equations that govern the flow in a triangular cavity shown in Figure 1(b) as follows:
(8)
∇
ξ
η
2
ψ
=
-
ω
,
(9)
∇
ξ
η
2
ω
=
J
(
∂
ψ
∂
η
∂
ω
∂
ξ
-
∂
ψ
∂
ξ
∂
ω
∂
η
)
,
where
(10)
∇
ξ
η
2
=
α
∂
2
∂
ξ
2
+
β
∂
2
∂
η
2
+
γ
∂
2
∂
ξ
∂
η
,
α
=
(
∂
ξ
∂
x
)
2
+
(
∂
ξ
∂
y
)
2
,
β
=
(
∂
η
∂
y
)
2
,
γ
=
2
∂
ξ
∂
y
∂
η
∂
y
,
J
=
Re
∂
ξ
∂
x
∂
η
∂
y
,
and we obtain
(11)
α
=
β
,
u
=
∂
ψ
∂
ξ
∂
ξ
∂
y
+
∂
ψ
∂
η
∂
η
∂
y
,
v
=
-
∂
ψ
∂
ξ
∂
ξ
∂
x
.The boundary conditions inx
y-plane are follows [11]:
(12)
ψ
=
0
on
all
three
sides
,
(13)
(
u
,
v
)
·
τ
=
{
1
on
the
top
side
,
0
on
the
other
two
sides
,
(14)
(
u
,
v
)
·
n
=
0
on
all
three
sides
,
where n is the outside normal unit vector and τ is the tangential unit vector in clockwise direction. In the ξ
η-plane condition (11) remains the same form, while (13)-(14) need to be converted as follows. On the top side A
^
B
^, substituting (11) into (13)-(14), with τ
=
(
1,0
) and n
=
(
0,1
) yields
(15)
∂
ψ
∂
ξ
=
0
,
∂
ψ
∂
η
=
∂
y
∂
η
=
(
x
b
-
x
a
)
(
y
c
-
h
)
(
y
c
-
h
)
2
+
(
x
c
-
x
a
)
2
on
side
A
B
.
On side B
^
C
^ we have
(16)
τ
=
(
x
c
-
x
b
,
y
c
-
h
)
(
x
c
-
x
b
)
2
+
(
y
c
-
h
)
2
,
n
=
(
h
-
y
c
,
x
c
-
x
b
)
(
x
c
-
x
b
)
2
+
(
y
c
-
h
)
2
.
Combining this with (11) we obtain
(17)
∂
ψ
∂
ξ
=
0
,
∂
ψ
∂
η
=
0
on
side
B
C
.
Similarly on side C
^
A
^ by using
(18)
τ
=
(
x
a
-
x
c
,
h
-
y
c
)
(
x
a
-
x
c
)
2
+
(
y
c
-
h
)
2
,
n
=
(
y
c
-
h
,
x
a
-
x
c
)
(
x
a
-
x
c
)
2
+
(
y
c
-
h
)
2
.
we obtain
(19)
∂
ψ
∂
ξ
=
0
,
∂
ψ
∂
η
=
0
on
side
C
A
.


## 3. Higher-Order Semicompact Scheme

### 3.1. Discretization of the Stream Function Equation

It is known that the higher-order compact approximation for any Poisson type of operator acting onf is given by [12, 13]
(20)
D
x
2
(
1
+
Δ
h
2
12
D
y
2
)
f
+
D
y
2
(
1
+
Δ
h
2
12
D
x
2
)
f
=
(
1
+
Δ
h
2
12
(
D
x
2
+
D
y
2
)
)
g
,
where g is any known function, Δ
h is the constant step length in x and y directions, and
(21)
D
x
2
f
(
x
,
y
)
=
1
Δ
h
(
f
(
x
-
Δ
h
,
y
)
-
2
f
(
x
,
y
)
D
x
2
f
(
x
,
y
)
=
1
Δ
h
v
+
f
(
x
+
Δ
h
,
y
)
)
,
D
y
2
f
(
x
,
y
)
=
1
Δ
h
(
f
(
x
,
y
-
Δ
h
)
-
2
f
(
x
,
y
)
D
y
2
f
(
x
,
y
)
=
1
Δ
h
v
+
f
(
x
,
y
+
Δ
h
)
)
.
If we define
(22)
(
1
+
Δ
h
2
12
(
D
ξ
2
+
D
η
2
)
)
ω
=
ω
~
,
and then approximate the stream function (8) in the form of (25) we get
(23)
α
D
ξ
2
(
1
+
Δ
h
2
12
D
η
2
)
ψ
+
β
D
η
2
(
1
+
Δ
h
2
12
D
ξ
2
)
ψ
+
γ
~
ψ
=
-
ω
~
+
O
(
Δ
h
4
)
,
where
(24)
γ
~
=
γ
D
ξ
η
2
(
1
-
Δ
h
2
12
(
D
ξ
2
+
D
η
2
)
)
,
D
ξ
f
(
ξ
,
η
)
=
1
2
Δ
h
(
f
(
ξ
-
Δ
h
,
η
)
-
f
(
ξ
+
Δ
h
,
η
)
)
,
D
η
f
(
ξ
,
η
)
=
1
2
Δ
h
(
f
(
ξ
,
η
-
Δ
h
)
-
f
(
ξ
,
η
+
Δ
h
)
)
.
A 2D finite difference scheme is compact if it only involves at most the 8 nearest neighboring grid points (of the center point) in the approximation formula [14]. Then in (23), first and second terms on the left side of the equation can be discretized on the compact nine point stencil as in any other compact scheme. However, discretization of the third term γ
~
ψ needs some special attention. We refer this term as T
E
R
M
1 and the discretization of this term is considered after the discussion of the vorticity equation.


### 3.2. Discretization of the Vorticity Equation

Equation (9) with fourth-order approximation can be written as(25)
α
D
ξ
2
(
1
+
Δ
h
2
12
D
η
2
)
ω
+
β
D
η
2
(
1
+
Δ
h
2
12
D
ξ
2
)
ω
+
γ
~
ω
=
J
(
1
+
Δ
h
2
12
(
D
ξ
2
+
D
η
2
)
)
(
D
η
ψ
D
ξ
ω
-
D
ξ
ψ
D
η
ω
)
+
O
(
Δ
h
4
)
.
The first and second terms on the left hand side of (25) can be discretized within the compact stencil. We call the third term γ
~
ω as T
E
R
M
2 which needs a larger stencil to approximate to fourth-order. The term on the right hand side of (25) contains
(26)
(
D
ξ
2
+
D
η
2
)
(
D
η
ψ
D
ξ
ω
-
D
ξ
ψ
D
η
ω
)
=
2
D
ξ
η
2
ω
(
D
η
2
-
D
ξ
2
)
ψ
+
2
D
ξ
η
2
ψ
(
D
ξ
2
-
D
η
2
)
ω
-
D
ξ
(
D
ξ
2
+
D
η
2
)
ψ
D
η
ω
-
D
ξ
ψ
D
η
(
D
ξ
2
+
D
η
2
)
ω
+
D
η
(
D
ξ
2
+
D
η
2
)
ψ
D
ξ
ω
+
D
η
ψ
D
ξ
(
D
ξ
2
+
D
η
2
)
ω
.
The first two terms on the right hand side of (26) again can be discretized within the compact stencil to the fourth-order; however, the last four terms, call them as T
E
R
M
3, T
E
R
M
4, T
E
R
M
5, and T
E
R
M
6, need larger stencil.


### 3.3. Numerical Treatment of the TermsTERM1–TERM6

We obtain from (8) that
(27)
∂
2
ψ
∂
ξ
2
+
∂
2
ψ
∂
η
2
=
-
ω
α
-
γ
α
∂
2
ψ
∂
ξ
∂
η
.
Therefore,
(28)
D
ξ
η
2
(
D
ξ
2
+
D
η
2
)
ψ
=
-
1
α
D
ξ
η
2
ω
-
γ
α
D
ξ
ξ
η
η
4
ψ
.
Similarly, from (9) we obtain
(29)
D
ξ
η
2
(
D
ξ
2
+
D
η
2
)
ω
=
-
γ
α
D
ξ
ξ
η
η
4
ω
+
J
α
(
D
ξ
η
η
3
ψ
D
ξ
ω
+
D
η
ψ
D
ξ
ξ
η
3
ω
=
-
γ
α
D
ξ
ξ
η
η
4
ω
+
J
α
vvvvv
+
D
η
2
ψ
D
ξ
2
ω
)
-
J
α
(
D
ξ
ξ
η
3
ψ
D
η
ω
+
D
ξ
ψ
D
ξ
η
η
3
ω
+
D
ξ
2
ψ
D
η
2
ω
)
.
In (29), we have third derivatives such as D
ξ
(
D
ξ
2
+
D
η
2
)
ψ and D
η
(
D
ξ
2
+
D
η
2
)
ω which present some problems to the fourth-order compact formulation. In order to find an expression for these derivatives, we use (8)-(9) to obtain the following equations:
(30)
D
κ
(
D
ξ
2
+
D
η
2
)
ψ
=
-
1
α
D
κ
ω
-
γ
α
D
ξ
κ
η
3
ψ
,
κ
=
ξ
,
η
,
(31)
D
κ
(
D
ξ
2
+
D
η
2
)
ω
=
-
γ
α
D
ξ
κ
η
3
ω
+
J
α
(
D
ξ
ω
D
κ
η
2
ψ
+
D
ξ
κ
2
ω
D
η
ψ
)
-
J
α
(
D
ξ
κ
2
ψ
D
η
ω
+
D
ξ
ψ
D
κ
η
2
ω
)
,
D
ξ
κ
2
ψ
D
η
ω
+
D
ξ
i
κ
=
ξ
,
η
.


### 3.4. Numerical Boundary Conditions

The numerical implementation of the boundary conditions foru
,
v, and ψ are straightforward. On side A
B boundary, we have
(32)
ω
i
,
1
=
-
β
∂
2
ψ
i
,
1
∂
η
2
,
α
∂
2
ω
i
,
1
∂
ξ
2
+
β
∂
2
ω
i
,
1
∂
η
2
+
γ
∂
2
ω
i
,
1
∂
ξ
∂
η
=
J
K
∂
ω
i
,
1
∂
ξ
,
where i = 2,3
,
…
,
N and
(33)
K
=
∂
ψ
i
,
1
∂
η
=
(
x
b
-
x
a
)
(
y
c
-
h
)
(
y
c
-
h
)
2
+
(
x
c
-
x
a
)
2
.
Employing Taylor series expansion, we get
(34)
K
=
∂
ψ
i
,
1
∂
η
=
D
η
+
ψ
i
,
1
-
Δ
h
2
∂
2
ψ
i
,
1
∂
η
2
-
Δ
h
2
6
∂
3
ψ
i
,
1
∂
η
3
-
Δ
h
3
24
∂
4
ψ
i
,
1
∂
η
4
+
O
(
Δ
h
4
)
,
where D
η
+ is the first-order forward difference operator in the η-direction and Using (32) in (34), we get the fourth-order accuracy O
(
Δ
h
4
) expression as follows:
(35)
K
=
D
η
+
ψ
i
,
1
+
Δ
h
2
β
ω
i
,
1
+
Δ
h
2
6
β
D
η
+
ω
i
,
1
-
Δ
h
3
24
β
2
(
J
K
D
ξ
+
-
α
D
ξ
2
-
γ
δ
ξ
η
2
)
ω
i
,
1
,
where D
ξ
+ is the first-order forward difference operator in the ξ-direction, δ
ξ
η
2 is the first-order cross difference operator, that is
(36)
δ
ξ
η
2
ω
i
,
1
=
1
Δ
h
2
(
ω
i
,
2
-
ω
i
,
1
-
ω
i
-
1,2
+
ω
i
-
1,1
)
+
O
(
Δ
h
)
.
Substituting (36) into (35), we get
(37)
ω
i
,
1
=
Θ
-
(
J
K
Δ
h
-
α
)
ω
i
+
1,1
+
(
α
+
γ
)
ω
i
-
1,1
+
γ
ω
i
-
1,2
γ
-
J
K
Δ
h
+
2
α
-
8
β
,
where Θ
=
θ
(
ψ
i
,
2
-
K
Δ
h
)
+
(
4
β
+
γ
)
ω
i
,
2, θ
=
24
β
2
/
Δ
h
2.Similarly on sideA
C the vorticity is equal to
(38)
ω
1
,
j
=
θ
ψ
2
,
j
+
(
4
β
+
γ
)
ω
2
,
j
+
α
ω
1
,
j
+
1
+
(
α
+
γ
)
ω
1
,
j
-
1
+
γ
ω
2
,
j
-
1
γ
+
2
α
-
8
β
,
where j = 2,3
,
…
,
N. On side B
C, it is equal to
(39)
ω
i
,
k
=
θ
ψ
i
-
1
,
k
+
(
4
β
+
γ
)
ω
i
-
1
,
k
+
α
ω
i
,
k
-
1
+
(
α
+
γ
)
ω
i
,
k
+
1
+
γ
ω
i
-
1
,
k
+
1
γ
+
2
α
-
8
β
,
where k = N
+
2
-
i, i = 2,3
,
…
,
N. It is found that the point (
i
,
k
+
1
) is outside the 9-point domain being required and the value of ω
i
,
k
+
1 is not available; we replace ω
i
,
k
+
1 by using the following equation:
(40)
ω
i
,
k
+
1
=
2
ω
i
,
k
-
ω
i
,
k
-
1
+
ω
i
-
1
,
k
+
1
-
2
ω
i
-
1
,
k
+
ω
i
-
1
,
k
-
1
+
O
(
Δ
h
4
)
.
At the three singular corners, the values of vorticity need some special attention in the process of discretization. We use the average of vorticity on the two adjacent points with the singular corners.


## 3.1. Discretization of the Stream Function Equation

It is known that the higher-order compact approximation for any Poisson type of operator acting onf is given by [12, 13]
(20)
D
x
2
(
1
+
Δ
h
2
12
D
y
2
)
f
+
D
y
2
(
1
+
Δ
h
2
12
D
x
2
)
f
=
(
1
+
Δ
h
2
12
(
D
x
2
+
D
y
2
)
)
g
,
where g is any known function, Δ
h is the constant step length in x and y directions, and
(21)
D
x
2
f
(
x
,
y
)
=
1
Δ
h
(
f
(
x
-
Δ
h
,
y
)
-
2
f
(
x
,
y
)
D
x
2
f
(
x
,
y
)
=
1
Δ
h
v
+
f
(
x
+
Δ
h
,
y
)
)
,
D
y
2
f
(
x
,
y
)
=
1
Δ
h
(
f
(
x
,
y
-
Δ
h
)
-
2
f
(
x
,
y
)
D
y
2
f
(
x
,
y
)
=
1
Δ
h
v
+
f
(
x
,
y
+
Δ
h
)
)
.
If we define
(22)
(
1
+
Δ
h
2
12
(
D
ξ
2
+
D
η
2
)
)
ω
=
ω
~
,
and then approximate the stream function (8) in the form of (25) we get
(23)
α
D
ξ
2
(
1
+
Δ
h
2
12
D
η
2
)
ψ
+
β
D
η
2
(
1
+
Δ
h
2
12
D
ξ
2
)
ψ
+
γ
~
ψ
=
-
ω
~
+
O
(
Δ
h
4
)
,
where
(24)
γ
~
=
γ
D
ξ
η
2
(
1
-
Δ
h
2
12
(
D
ξ
2
+
D
η
2
)
)
,
D
ξ
f
(
ξ
,
η
)
=
1
2
Δ
h
(
f
(
ξ
-
Δ
h
,
η
)
-
f
(
ξ
+
Δ
h
,
η
)
)
,
D
η
f
(
ξ
,
η
)
=
1
2
Δ
h
(
f
(
ξ
,
η
-
Δ
h
)
-
f
(
ξ
,
η
+
Δ
h
)
)
.
A 2D finite difference scheme is compact if it only involves at most the 8 nearest neighboring grid points (of the center point) in the approximation formula [14]. Then in (23), first and second terms on the left side of the equation can be discretized on the compact nine point stencil as in any other compact scheme. However, discretization of the third term γ
~
ψ needs some special attention. We refer this term as T
E
R
M
1 and the discretization of this term is considered after the discussion of the vorticity equation.


## 3.2. Discretization of the Vorticity Equation

Equation (9) with fourth-order approximation can be written as(25)
α
D
ξ
2
(
1
+
Δ
h
2
12
D
η
2
)
ω
+
β
D
η
2
(
1
+
Δ
h
2
12
D
ξ
2
)
ω
+
γ
~
ω
=
J
(
1
+
Δ
h
2
12
(
D
ξ
2
+
D
η
2
)
)
(
D
η
ψ
D
ξ
ω
-
D
ξ
ψ
D
η
ω
)
+
O
(
Δ
h
4
)
.
The first and second terms on the left hand side of (25) can be discretized within the compact stencil. We call the third term γ
~
ω as T
E
R
M
2 which needs a larger stencil to approximate to fourth-order. The term on the right hand side of (25) contains
(26)
(
D
ξ
2
+
D
η
2
)
(
D
η
ψ
D
ξ
ω
-
D
ξ
ψ
D
η
ω
)
=
2
D
ξ
η
2
ω
(
D
η
2
-
D
ξ
2
)
ψ
+
2
D
ξ
η
2
ψ
(
D
ξ
2
-
D
η
2
)
ω
-
D
ξ
(
D
ξ
2
+
D
η
2
)
ψ
D
η
ω
-
D
ξ
ψ
D
η
(
D
ξ
2
+
D
η
2
)
ω
+
D
η
(
D
ξ
2
+
D
η
2
)
ψ
D
ξ
ω
+
D
η
ψ
D
ξ
(
D
ξ
2
+
D
η
2
)
ω
.
The first two terms on the right hand side of (26) again can be discretized within the compact stencil to the fourth-order; however, the last four terms, call them as T
E
R
M
3, T
E
R
M
4, T
E
R
M
5, and T
E
R
M
6, need larger stencil.


## 3.3. Numerical Treatment of the TermsTERM1–TERM6

We obtain from (8) that
(27)
∂
2
ψ
∂
ξ
2
+
∂
2
ψ
∂
η
2
=
-
ω
α
-
γ
α
∂
2
ψ
∂
ξ
∂
η
.
Therefore,
(28)
D
ξ
η
2
(
D
ξ
2
+
D
η
2
)
ψ
=
-
1
α
D
ξ
η
2
ω
-
γ
α
D
ξ
ξ
η
η
4
ψ
.
Similarly, from (9) we obtain
(29)
D
ξ
η
2
(
D
ξ
2
+
D
η
2
)
ω
=
-
γ
α
D
ξ
ξ
η
η
4
ω
+
J
α
(
D
ξ
η
η
3
ψ
D
ξ
ω
+
D
η
ψ
D
ξ
ξ
η
3
ω
=
-
γ
α
D
ξ
ξ
η
η
4
ω
+
J
α
vvvvv
+
D
η
2
ψ
D
ξ
2
ω
)
-
J
α
(
D
ξ
ξ
η
3
ψ
D
η
ω
+
D
ξ
ψ
D
ξ
η
η
3
ω
+
D
ξ
2
ψ
D
η
2
ω
)
.
In (29), we have third derivatives such as D
ξ
(
D
ξ
2
+
D
η
2
)
ψ and D
η
(
D
ξ
2
+
D
η
2
)
ω which present some problems to the fourth-order compact formulation. In order to find an expression for these derivatives, we use (8)-(9) to obtain the following equations:
(30)
D
κ
(
D
ξ
2
+
D
η
2
)
ψ
=
-
1
α
D
κ
ω
-
γ
α
D
ξ
κ
η
3
ψ
,
κ
=
ξ
,
η
,
(31)
D
κ
(
D
ξ
2
+
D
η
2
)
ω
=
-
γ
α
D
ξ
κ
η
3
ω
+
J
α
(
D
ξ
ω
D
κ
η
2
ψ
+
D
ξ
κ
2
ω
D
η
ψ
)
-
J
α
(
D
ξ
κ
2
ψ
D
η
ω
+
D
ξ
ψ
D
κ
η
2
ω
)
,
D
ξ
κ
2
ψ
D
η
ω
+
D
ξ
i
κ
=
ξ
,
η
.


## 3.4. Numerical Boundary Conditions

The numerical implementation of the boundary conditions foru
,
v, and ψ are straightforward. On side A
B boundary, we have
(32)
ω
i
,
1
=
-
β
∂
2
ψ
i
,
1
∂
η
2
,
α
∂
2
ω
i
,
1
∂
ξ
2
+
β
∂
2
ω
i
,
1
∂
η
2
+
γ
∂
2
ω
i
,
1
∂
ξ
∂
η
=
J
K
∂
ω
i
,
1
∂
ξ
,
where i = 2,3
,
…
,
N and
(33)
K
=
∂
ψ
i
,
1
∂
η
=
(
x
b
-
x
a
)
(
y
c
-
h
)
(
y
c
-
h
)
2
+
(
x
c
-
x
a
)
2
.
Employing Taylor series expansion, we get
(34)
K
=
∂
ψ
i
,
1
∂
η
=
D
η
+
ψ
i
,
1
-
Δ
h
2
∂
2
ψ
i
,
1
∂
η
2
-
Δ
h
2
6
∂
3
ψ
i
,
1
∂
η
3
-
Δ
h
3
24
∂
4
ψ
i
,
1
∂
η
4
+
O
(
Δ
h
4
)
,
where D
η
+ is the first-order forward difference operator in the η-direction and Using (32) in (34), we get the fourth-order accuracy O
(
Δ
h
4
) expression as follows:
(35)
K
=
D
η
+
ψ
i
,
1
+
Δ
h
2
β
ω
i
,
1
+
Δ
h
2
6
β
D
η
+
ω
i
,
1
-
Δ
h
3
24
β
2
(
J
K
D
ξ
+
-
α
D
ξ
2
-
γ
δ
ξ
η
2
)
ω
i
,
1
,
where D
ξ
+ is the first-order forward difference operator in the ξ-direction, δ
ξ
η
2 is the first-order cross difference operator, that is
(36)
δ
ξ
η
2
ω
i
,
1
=
1
Δ
h
2
(
ω
i
,
2
-
ω
i
,
1
-
ω
i
-
1,2
+
ω
i
-
1,1
)
+
O
(
Δ
h
)
.
Substituting (36) into (35), we get
(37)
ω
i
,
1
=
Θ
-
(
J
K
Δ
h
-
α
)
ω
i
+
1,1
+
(
α
+
γ
)
ω
i
-
1,1
+
γ
ω
i
-
1,2
γ
-
J
K
Δ
h
+
2
α
-
8
β
,
where Θ
=
θ
(
ψ
i
,
2
-
K
Δ
h
)
+
(
4
β
+
γ
)
ω
i
,
2, θ
=
24
β
2
/
Δ
h
2.Similarly on sideA
C the vorticity is equal to
(38)
ω
1
,
j
=
θ
ψ
2
,
j
+
(
4
β
+
γ
)
ω
2
,
j
+
α
ω
1
,
j
+
1
+
(
α
+
γ
)
ω
1
,
j
-
1
+
γ
ω
2
,
j
-
1
γ
+
2
α
-
8
β
,
where j = 2,3
,
…
,
N. On side B
C, it is equal to
(39)
ω
i
,
k
=
θ
ψ
i
-
1
,
k
+
(
4
β
+
γ
)
ω
i
-
1
,
k
+
α
ω
i
,
k
-
1
+
(
α
+
γ
)
ω
i
,
k
+
1
+
γ
ω
i
-
1
,
k
+
1
γ
+
2
α
-
8
β
,
where k = N
+
2
-
i, i = 2,3
,
…
,
N. It is found that the point (
i
,
k
+
1
) is outside the 9-point domain being required and the value of ω
i
,
k
+
1 is not available; we replace ω
i
,
k
+
1 by using the following equation:
(40)
ω
i
,
k
+
1
=
2
ω
i
,
k
-
ω
i
,
k
-
1
+
ω
i
-
1
,
k
+
1
-
2
ω
i
-
1
,
k
+
ω
i
-
1
,
k
-
1
+
O
(
Δ
h
4
)
.
At the three singular corners, the values of vorticity need some special attention in the process of discretization. We use the average of vorticity on the two adjacent points with the singular corners.


## 4. Results and Discussion

Following the numerical procedure described in the previous section, the new schemes (23) and (25) are fourth-order accurate; therefore, they need to be solved in an iterative manner such as SOR iteration with the relaxation parameter ϖ = 2/3. In all of the cases considered, we start the iterations from a homogeneous initial guess and continue until a certain condition of convergence is satisfied. As convergence criteria we decided to use the difference of the ψ-ω variables between two steps normalized by the previous value of the corresponding variable, such that
(41)
max
(
|
ψ
n
+
1
-
ψ
n
|
)
<
1
0
-
8
,
max
(
|
ω
n
+
1
-
ω
n
|
)
<
1
0
-
8
.
At this convergence level, this would indicate that the variables ψ and ω are changing less than 0.000001
% of their value between two iterations at every grid point in the mesh. These residuals indicate the degree to which the numerical solution has converged to steady state.We considered different triangle geometries with different Reynolds numbers. Numerical tests for a variety of triangular geometries have been investigated, withRe up to 2000 for an equilateral cavity and 1500 for scalene cavities. It should be pointed out that Reynolds number is based on h, which is consistent with the definition in the case of the equilateral cavity. We solved the flow in this triangle cavity at various Reynolds numbers ranging between 1 and 2000. We note that if the length of one side of the triangle 2
3 was used in nondimensionalization, as it was used by Erturk and Gokcol [11], then our Reynolds number of 2000 would 2
3-fold such that it would correspond to a Reynolds number of 6928.


### 4.1. Equilateral Triangular Cavity

Let us first consider a nondimensional equilateral triangle with coordinates of corner points(42)
x
a
=
-
3
,
x
b
=
3
,
h
=
1
,
x
c
=
0
,
y
c
=
-
2
,
which was also considered by McQuain et al. [6], Jyotsna and Vanka [7], Li and Tang [8], Gaskell et al. [9], and Ribbens et al. [10]. We solved the flow in this triangle cavity at various Reynolds numbers ranging between 1 and 2000.In order to verify the accuracy of the present numerical study, we consider arbitrary triangular cavity flow with the following analytical solutions [15]:
(43)
ψ
=
y
-
x
Re
-
e
x
+
y
,
ω
=
2
e
x
+
y
.
For this model problem, as the boundary conditions we decided to use the analytical solutions defined in (43) at the grid points on the boundaries. This way we would be able to avoid any effect of numerical boundary condition approximations on the numerical solution and concentrate on the accuracy of the solution of both formulations in the interior domain for given analytical values at the boundaries. We note that in (43), the vorticity ω is independent of Re and the solution of the streamfunction ψ at different Reynolds numbers looks almost the same in a contour plot. We solve this model problem using different grid meshes, 11
×
11
/
2, 21 × 21/2, 41 × 41/2, 81
×
81
/
2, respectively. In these higher-order solutions, the average absolute (E
A
A) error and the the maximum norm (E
MAX) error between the exact solution Ex
a
k given in (43) and the numerical solution So
l
k obtained from (23) and (25) are defined as
(44)
E
A
A
=
1
S
∑
k
‍
|
Ex
a
k
-
So
l
k
|
,
E
MAX
=
max
k
|
Ex
a
k
-
So
l
k
|
,
where S being the total number of grid points. Since the grid size is decreased by a factor of 2, we can calculate the convergence rate m using the following formula [16]:
(45)
m
=
log
(
E
Δ
h
/
E
Δ
h
/
2
)
log
2
.
As observed by Wan and Zhou [17], the accuracy of the present results is getting better as the grid number increases, and we can see that when the grid spacing is decreased progressively by half, the scheme maintains fourth-order of spatial accurate; the convergence rate is very close to m
≈
4. In Figure 2, the average absolute error E
A
A and the the maximum norm error E
MAX are plotted with respect to the grid spacing in a log-log scale [16]. From Figure 2, we can again clearly see that the fully HOC formulation indeed provide fourth-order accurate solutions; the slope between log
E and log
h is close to m
≈
4.The average absolute errorE
A
A and the maximum norm error E
MAX in streamfunction and vorticity variables with respect to grid spacing. (a) streamfunction and (b) vorticity.
(a)
(b)Figures3 and 4 show the streamline and vorticity contours of the triangular cavity flow for a variety of Reynolds numbers obtained with using (
128
×
128
)
/
2 grids points. These figures show the streams and vortices that as the Reynolds number increases. From these contour figures, we conclude that the fourth-order compact formulation provides very smooth solutions and it is seen that as Re becomes large, the location of the center of the primary eddy and its streamfunction value seem to have converged. In terms of quantitative analysis, Table 1 tabulates the center of the primary eddy and the streamfunction and vorticity values at the core, together with results found in the literature. Our results are in good agreement with the results in [6] and [9] up to the maximum Reynolds number (Re
=
500).Table 1
Properties of the center of the primary eddy, located at(
x
,
y
) with streamfunction value and vorticity for equilateral triangle.
Re
Source
ψ
ω
(
x
,
y
)
1
Present
−0.2348
−1.4335
(0.0271, 0.4844)
Li and Tang [8]
−0.2340
−1.3990
(0.0000, 0.4750)
Erturk and Gokcol  [11]
−0.2329
−1.3788
(0.0101, 0.4668)
Gaskell et al. [9]
−0.2320
−1.3750
(0.0100, 0.4670)
McQuain et al.   [6]
−0.2330
−1.3630
(0.0170, 0.4600)
50
Present
−0.2377
−1.4533
(0.3248, 0.4375)
Li and Tang [8]
−0.2350
−1.4380
(0.3680, 0.4380)
Erturk and Gokcol  [11]
−0.2369
−1.4689
(0.3484, 0.4434)
Gaskell et al. [9]
−0.2360
−1.4690
(0.3500, 0.4420)
McQuain et al.   [6]
−0.2370
−1.4640
(0.3460, 0.4450)
100
Present
−0.2489
−1.3795
(0.3291, 0.3700)
Li and Tang [8]
−0.2440
−1.2640
(0.3680, 0.3630)
Erturk and Gokcol  [11]
−0.2482
−1.3669
(0.3315, 0.3555)
Gaskell et al. [9]
−0.2490
−1.3650
(0.3330, 0.3560)
McQuain et al.   [6]
−0.2470
−1.3730
(0.3290, 0.3550)
350
Present
−0.2663
−1.2388
(0.1624, 0.2500)
Li and Tang [8]
−0.2740
−1.1530
(0.1520, 0.2130)
Erturk and Gokcol  [11]
−0.2724
−1.1985
(0.1556, 0.2383)
Gaskell et al. [9]
−0.2740
−1.2010
(0.1490, 0.2340)
McQuain et al.   [6]
−0.2680
−1.2320
(0.1730, 0.2650)
500
Present
−0.2748
−1.1992
(0.1299, 0.2350)
Li and Tang [8]
−0.2780
−1.1240
(0.1080, 0.2130)
Erturk and Gokcol  [11]
−0.2774
−1.1791
(0.1319, 0.2207)
Gaskell et al. [9]
−0.2810
−1.1870
(0.1280, 0.2170)
McQuain et al.   [6]
−0.2690
−1.2500
(0.1730, 0.2650)
1000
Present
−0.2902
−1.2067
(0.0866, 0.2200)
Li and Tang [8]
−0.2790
−1.0480
(0.1080, 0.1380)
Erturk and Gokcol  [11]
−0.2844
−1.1629
(0.1116, 0.1973)
1500
Present
−0.2962
−1.2282
(0.0664, 0.2250)
Li and Tang [8]
−0.2770
−0.9980
(0.1080, 0.1380)
Erturk and Gokcol  [11]
−0.2873
−1.1639
(0.1015, 0.1914)
2000
Present
−0.2951
−1.1827
(0.0693, 0.2029)
Li and Tang [8]
—
—
—
Erturk and Gokcol  [11]
—
—
—Streamline contours of triangular cavity flow at various Reynolds numbers.
(a)
(b)
(c)
(d)Vorticity contours of triangular cavity flow at various Reynolds numbers.
(a)
(b)
(c)
(d)In Figure5(a), we plot the vorticity values at the center of the primary vortex tabulated in Table 1, with respect to Reynolds number. The results of Li and Tang [8] start to behave differently starting from Re
=
100 from the rest of the results. Also in Figure 5(a), the results of McQuain et al. [6] start to deviate from present results and Gaskell et al.’s [9] results after Re
=
200. Having a different behavior, the vorticity value of McQuain et al. [6] shows an increase such that their vorticity value at Re
=
500 is greater than the vorticity value at Re
=
350, whereas the present results and Erturk et al.’s [11] results show a continuous decrease until Re
≤
1500. We believe that these behaviors are due to the coarse grids used in those studies and in order to resolve these behaviors we decided to solve the same flow problem using several coarse grid meshes. For a given grid mesh we have solved the flow for increasing Reynolds number until a particular Re, where the solution was not convergent but oscillating. For this particular Re when the number of grid points was increased, the convergence was recovered and we were able to obtain a solution.Comparison of the vorticity at the center of the primary eddy in the equilateral triangular cavity (a) and right-oriented triangular cavity (b).
(a)
(b)According to the mean square law [18], the value of vorticity is approximately constant at the primary vortex. For an equilateral cavity with length of side 2
3 this constant is found to be 1.054 by McQuain et al. [6]. This infinite Reynolds number value of the vorticity of 1.054 is shown with the dotted line in Figure 5(a). We note that especially the eddies at the bottom corner occupy some portion of the corner as Reynolds number increases. However, the increase in the size of the portion of the bottom corner eddies almost stops after Re
=
500 and the size of the primary eddy remains almost constant beyond and so is the value of the vorticity at the core of the primary eddy. It looks like, for an equilateral triangular cavity flow at high Reynolds numbers, the mean square law predicts the strength of the primary eddy within an error due to the effect of the secondary eddies. For circular or elliptic boundaries Ribbens et al. [10], for square cavity flow Erturk and Gokcol [11], and for rectangle cavities McQuain et al. [6] have shown that the mean square law is approximately valid and successful in predicting the strength of the primary eddy at high Reynolds numbers.


### 4.2. Scalene Triangle Cavity

Owing to the unsymmetric geometry, flow motion in a scalene triangular cavity is very difficult to simulate. It was conjectured in McQuain et al. [6] and Ribbens et al. [10] that it might be unavoidably ill-conditioned. Then we considered a scalene triangle with coordinates of corner points
(46)
x
a
=
-
3
,
x
b
=
3
,
h
=
1
,
x
c
=
3
2
,
y
c
=
-
2
.
Using this triangle geometry, we could carry out the calculations for Reynolds numbers as high as 800. In Figure 5(b), we plot the vorticity values at the center of the primary vortex obtained with present study and Li and Tang. [8]. We can see from Figure 5(b) that the vorticity values are getting more and more close to the theoretical value as the Reynolds number increases. In Figure 6, we plot the streamlines and vorticity contours for Re
=
350 and 500. Again it is seen that as Re becomes large, the location of the center of the primary eddy and its streamfunction value seem to have converged.Streamlines distributions for a right-oriented triangular cavity at different Reynolds numbers.
(a)
(b)


### 4.3. Isosceles Triangle Cavity

We first considered an isosceles triangle which was also considered by Jyotsna and Vanka [7], Gaskell et al. [9], and Erturk and Gokcol [11], where
(47)
x
a
=
-
1
,
x
b
=
1
,
h
=
0
,
x
c
=
0
,
y
c
=
-
4
.
We note that our definition of Reynolds number is equivalent to one fourth of the Reynolds number definition used by Jyotsna and Vanka [7]. Figure 7 shows the streamline contours of the flow in this triangle at various Reynolds numbers. Comparing the location of the primary eddy center with that of Jyotsna and Vanka [7], Gaskell et al. [9], and Erturk and Gokcol [11], we believe that our results are more accurate.Contour figures of isosceles triangle.
(a)
(b)Then we considered an isosceles triangle which was also considered by Gaskell et al. [9] and Erturk and Gokcol. [11], with
(48)
x
a
=
-
1
,
x
b
=
1
,
h
=
0
,
x
c
=
0
,
y
c
=
-
cot
θ
2
,
where θ
= 5°, 10°, 20°, 40°, 70°, 100°, 140°, 160°, and 170°. Figure 8 qualitatively shows the change in the flow with Re
=
200 in an isosceles triangle as the corner angle, θ, changes.Contour figures of isosceles triangle with various corner angle.
(a)
(b)
(c)
(d)
(e)
(f)We also considered an isosceles right triangle with the 90° corner being at the top right corner, such that with corner points [11]
(49)
x
a
=
0
,
x
b
=
1
,
h
=
1
,
x
c
=
1
,
y
c
=
0
.
Figure 9 shows the flow topology as a function of Reynolds numbers and at last, we also considered an isosceles right triangle with the 90° corner being at the top left corner, such that with corner points [11]
(50)
x
a
=
0
,
x
b
=
1
,
h
=
1
,
x
c
=
0
,
y
c
=
0
.
Using this triangle geometry, we were able to obtain the solutions and Figure 10 shows the flow topology as a function of Reynolds numbers. For the considered triangle geometry, we can see that the eddy is closest to the moving lid. As it is obvious from these figures, in both cases, the flow behaves very differently as Reynolds number increases, which shows that the flow structures in a triangle cavity are greatly affected by the triangle geometry.Contour figures of right hand side aligned right triangle at different Reynolds numbers.
(a)
(b)Contour figures of left hand side aligned right triangle at various Reynolds numbers.
(a)
(b)


## 4.1. Equilateral Triangular Cavity

Let us first consider a nondimensional equilateral triangle with coordinates of corner points(42)
x
a
=
-
3
,
x
b
=
3
,
h
=
1
,
x
c
=
0
,
y
c
=
-
2
,
which was also considered by McQuain et al. [6], Jyotsna and Vanka [7], Li and Tang [8], Gaskell et al. [9], and Ribbens et al. [10]. We solved the flow in this triangle cavity at various Reynolds numbers ranging between 1 and 2000.In order to verify the accuracy of the present numerical study, we consider arbitrary triangular cavity flow with the following analytical solutions [15]:
(43)
ψ
=
y
-
x
Re
-
e
x
+
y
,
ω
=
2
e
x
+
y
.
For this model problem, as the boundary conditions we decided to use the analytical solutions defined in (43) at the grid points on the boundaries. This way we would be able to avoid any effect of numerical boundary condition approximations on the numerical solution and concentrate on the accuracy of the solution of both formulations in the interior domain for given analytical values at the boundaries. We note that in (43), the vorticity ω is independent of Re and the solution of the streamfunction ψ at different Reynolds numbers looks almost the same in a contour plot. We solve this model problem using different grid meshes, 11
×
11
/
2, 21 × 21/2, 41 × 41/2, 81
×
81
/
2, respectively. In these higher-order solutions, the average absolute (E
A
A) error and the the maximum norm (E
MAX) error between the exact solution Ex
a
k given in (43) and the numerical solution So
l
k obtained from (23) and (25) are defined as
(44)
E
A
A
=
1
S
∑
k
‍
|
Ex
a
k
-
So
l
k
|
,
E
MAX
=
max
k
|
Ex
a
k
-
So
l
k
|
,
where S being the total number of grid points. Since the grid size is decreased by a factor of 2, we can calculate the convergence rate m using the following formula [16]:
(45)
m
=
log
(
E
Δ
h
/
E
Δ
h
/
2
)
log
2
.
As observed by Wan and Zhou [17], the accuracy of the present results is getting better as the grid number increases, and we can see that when the grid spacing is decreased progressively by half, the scheme maintains fourth-order of spatial accurate; the convergence rate is very close to m
≈
4. In Figure 2, the average absolute error E
A
A and the the maximum norm error E
MAX are plotted with respect to the grid spacing in a log-log scale [16]. From Figure 2, we can again clearly see that the fully HOC formulation indeed provide fourth-order accurate solutions; the slope between log
E and log
h is close to m
≈
4.The average absolute errorE
A
A and the maximum norm error E
MAX in streamfunction and vorticity variables with respect to grid spacing. (a) streamfunction and (b) vorticity.
(a)
(b)Figures3 and 4 show the streamline and vorticity contours of the triangular cavity flow for a variety of Reynolds numbers obtained with using (
128
×
128
)
/
2 grids points. These figures show the streams and vortices that as the Reynolds number increases. From these contour figures, we conclude that the fourth-order compact formulation provides very smooth solutions and it is seen that as Re becomes large, the location of the center of the primary eddy and its streamfunction value seem to have converged. In terms of quantitative analysis, Table 1 tabulates the center of the primary eddy and the streamfunction and vorticity values at the core, together with results found in the literature. Our results are in good agreement with the results in [6] and [9] up to the maximum Reynolds number (Re
=
500).Table 1
Properties of the center of the primary eddy, located at(
x
,
y
) with streamfunction value and vorticity for equilateral triangle.
Re
Source
ψ
ω
(
x
,
y
)
1
Present
−0.2348
−1.4335
(0.0271, 0.4844)
Li and Tang [8]
−0.2340
−1.3990
(0.0000, 0.4750)
Erturk and Gokcol  [11]
−0.2329
−1.3788
(0.0101, 0.4668)
Gaskell et al. [9]
−0.2320
−1.3750
(0.0100, 0.4670)
McQuain et al.   [6]
−0.2330
−1.3630
(0.0170, 0.4600)
50
Present
−0.2377
−1.4533
(0.3248, 0.4375)
Li and Tang [8]
−0.2350
−1.4380
(0.3680, 0.4380)
Erturk and Gokcol  [11]
−0.2369
−1.4689
(0.3484, 0.4434)
Gaskell et al. [9]
−0.2360
−1.4690
(0.3500, 0.4420)
McQuain et al.   [6]
−0.2370
−1.4640
(0.3460, 0.4450)
100
Present
−0.2489
−1.3795
(0.3291, 0.3700)
Li and Tang [8]
−0.2440
−1.2640
(0.3680, 0.3630)
Erturk and Gokcol  [11]
−0.2482
−1.3669
(0.3315, 0.3555)
Gaskell et al. [9]
−0.2490
−1.3650
(0.3330, 0.3560)
McQuain et al.   [6]
−0.2470
−1.3730
(0.3290, 0.3550)
350
Present
−0.2663
−1.2388
(0.1624, 0.2500)
Li and Tang [8]
−0.2740
−1.1530
(0.1520, 0.2130)
Erturk and Gokcol  [11]
−0.2724
−1.1985
(0.1556, 0.2383)
Gaskell et al. [9]
−0.2740
−1.2010
(0.1490, 0.2340)
McQuain et al.   [6]
−0.2680
−1.2320
(0.1730, 0.2650)
500
Present
−0.2748
−1.1992
(0.1299, 0.2350)
Li and Tang [8]
−0.2780
−1.1240
(0.1080, 0.2130)
Erturk and Gokcol  [11]
−0.2774
−1.1791
(0.1319, 0.2207)
Gaskell et al. [9]
−0.2810
−1.1870
(0.1280, 0.2170)
McQuain et al.   [6]
−0.2690
−1.2500
(0.1730, 0.2650)
1000
Present
−0.2902
−1.2067
(0.0866, 0.2200)
Li and Tang [8]
−0.2790
−1.0480
(0.1080, 0.1380)
Erturk and Gokcol  [11]
−0.2844
−1.1629
(0.1116, 0.1973)
1500
Present
−0.2962
−1.2282
(0.0664, 0.2250)
Li and Tang [8]
−0.2770
−0.9980
(0.1080, 0.1380)
Erturk and Gokcol  [11]
−0.2873
−1.1639
(0.1015, 0.1914)
2000
Present
−0.2951
−1.1827
(0.0693, 0.2029)
Li and Tang [8]
—
—
—
Erturk and Gokcol  [11]
—
—
—Streamline contours of triangular cavity flow at various Reynolds numbers.
(a)
(b)
(c)
(d)Vorticity contours of triangular cavity flow at various Reynolds numbers.
(a)
(b)
(c)
(d)In Figure5(a), we plot the vorticity values at the center of the primary vortex tabulated in Table 1, with respect to Reynolds number. The results of Li and Tang [8] start to behave differently starting from Re
=
100 from the rest of the results. Also in Figure 5(a), the results of McQuain et al. [6] start to deviate from present results and Gaskell et al.’s [9] results after Re
=
200. Having a different behavior, the vorticity value of McQuain et al. [6] shows an increase such that their vorticity value at Re
=
500 is greater than the vorticity value at Re
=
350, whereas the present results and Erturk et al.’s [11] results show a continuous decrease until Re
≤
1500. We believe that these behaviors are due to the coarse grids used in those studies and in order to resolve these behaviors we decided to solve the same flow problem using several coarse grid meshes. For a given grid mesh we have solved the flow for increasing Reynolds number until a particular Re, where the solution was not convergent but oscillating. For this particular Re when the number of grid points was increased, the convergence was recovered and we were able to obtain a solution.Comparison of the vorticity at the center of the primary eddy in the equilateral triangular cavity (a) and right-oriented triangular cavity (b).
(a)
(b)According to the mean square law [18], the value of vorticity is approximately constant at the primary vortex. For an equilateral cavity with length of side 2
3 this constant is found to be 1.054 by McQuain et al. [6]. This infinite Reynolds number value of the vorticity of 1.054 is shown with the dotted line in Figure 5(a). We note that especially the eddies at the bottom corner occupy some portion of the corner as Reynolds number increases. However, the increase in the size of the portion of the bottom corner eddies almost stops after Re
=
500 and the size of the primary eddy remains almost constant beyond and so is the value of the vorticity at the core of the primary eddy. It looks like, for an equilateral triangular cavity flow at high Reynolds numbers, the mean square law predicts the strength of the primary eddy within an error due to the effect of the secondary eddies. For circular or elliptic boundaries Ribbens et al. [10], for square cavity flow Erturk and Gokcol [11], and for rectangle cavities McQuain et al. [6] have shown that the mean square law is approximately valid and successful in predicting the strength of the primary eddy at high Reynolds numbers.


## 4.2. Scalene Triangle Cavity

Owing to the unsymmetric geometry, flow motion in a scalene triangular cavity is very difficult to simulate. It was conjectured in McQuain et al. [6] and Ribbens et al. [10] that it might be unavoidably ill-conditioned. Then we considered a scalene triangle with coordinates of corner points
(46)
x
a
=
-
3
,
x
b
=
3
,
h
=
1
,
x
c
=
3
2
,
y
c
=
-
2
.
Using this triangle geometry, we could carry out the calculations for Reynolds numbers as high as 800. In Figure 5(b), we plot the vorticity values at the center of the primary vortex obtained with present study and Li and Tang. [8]. We can see from Figure 5(b) that the vorticity values are getting more and more close to the theoretical value as the Reynolds number increases. In Figure 6, we plot the streamlines and vorticity contours for Re
=
350 and 500. Again it is seen that as Re becomes large, the location of the center of the primary eddy and its streamfunction value seem to have converged.Streamlines distributions for a right-oriented triangular cavity at different Reynolds numbers.
(a)
(b)


## 4.3. Isosceles Triangle Cavity

We first considered an isosceles triangle which was also considered by Jyotsna and Vanka [7], Gaskell et al. [9], and Erturk and Gokcol [11], where
(47)
x
a
=
-
1
,
x
b
=
1
,
h
=
0
,
x
c
=
0
,
y
c
=
-
4
.
We note that our definition of Reynolds number is equivalent to one fourth of the Reynolds number definition used by Jyotsna and Vanka [7]. Figure 7 shows the streamline contours of the flow in this triangle at various Reynolds numbers. Comparing the location of the primary eddy center with that of Jyotsna and Vanka [7], Gaskell et al. [9], and Erturk and Gokcol [11], we believe that our results are more accurate.Contour figures of isosceles triangle.
(a)
(b)Then we considered an isosceles triangle which was also considered by Gaskell et al. [9] and Erturk and Gokcol. [11], with
(48)
x
a
=
-
1
,
x
b
=
1
,
h
=
0
,
x
c
=
0
,
y
c
=
-
cot
θ
2
,
where θ
= 5°, 10°, 20°, 40°, 70°, 100°, 140°, 160°, and 170°. Figure 8 qualitatively shows the change in the flow with Re
=
200 in an isosceles triangle as the corner angle, θ, changes.Contour figures of isosceles triangle with various corner angle.
(a)
(b)
(c)
(d)
(e)
(f)We also considered an isosceles right triangle with the 90° corner being at the top right corner, such that with corner points [11]
(49)
x
a
=
0
,
x
b
=
1
,
h
=
1
,
x
c
=
1
,
y
c
=
0
.
Figure 9 shows the flow topology as a function of Reynolds numbers and at last, we also considered an isosceles right triangle with the 90° corner being at the top left corner, such that with corner points [11]
(50)
x
a
=
0
,
x
b
=
1
,
h
=
1
,
x
c
=
0
,
y
c
=
0
.
Using this triangle geometry, we were able to obtain the solutions and Figure 10 shows the flow topology as a function of Reynolds numbers. For the considered triangle geometry, we can see that the eddy is closest to the moving lid. As it is obvious from these figures, in both cases, the flow behaves very differently as Reynolds number increases, which shows that the flow structures in a triangle cavity are greatly affected by the triangle geometry.Contour figures of right hand side aligned right triangle at different Reynolds numbers.
(a)
(b)Contour figures of left hand side aligned right triangle at various Reynolds numbers.
(a)
(b)


## 5. Conclusion

In this work we have developed a new semicompact fourth-order scheme for the time-independentψ-ω form of the 2D, incompressible N-S equations governing the fluid flow in a triangular driven cavity. Our numerical scheme has been proved robust for a wide range of Reynolds numbers and applicable to triangular cavities with arbitrary shape. The key point with the present scheme is that it allows direct iteration for low-to-high Reynolds numbers. It is the success of the current method with the wide range of Re and mesh sizes that indicates the potential of this method as an accurate and stable numerical method applicable to a wide range of problems. We have tested the present method for both the equilateral triangular cavity problem and the scalene triangular cavity problem, and excellent agreement is found in all the cases, both qualitatively and quantitatively.



---

*Source: 349242-2014-02-04.xml*