# Symplectic Exact Solution for Stokes Flow in the Thin Film Coating Applications

Authors: Yan Wang; Zi-Chen Deng; Wei-Peng Hu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The symplectic analytical method is introduced to solve the problem of the stokes flow in the thin film coating applications. Based on the variational principle, the Lagrangian function of the stokes flow is established. By using the Legendre transformation, the dual variables of velocities and the Hamiltonian function are derived. Considering velocities and stresses as the basic variables, the equations of stokes flow problems are transformed into Hamiltonian system. The method of separation of variables and expansion of eigenfunctions are developed to solve the governing equations in Hamiltonian system, and the analytical solutions of the stokes flow are obtained. Several numerical simulations are carried out to verify the analytical solutions in the present study and discuss the effects of the driven lids of the square cavity on the dynamic behavior of the flow structure.


---

## Body

## 1. Introduction

The problems of the stokes flow in rectangular-shaped cavities have received considerable attention from engineers and researchers. For example, in coating systems [1], polymer melts [2], ceramic tape casting [3], and so on. Thus, the analysis of mechanical properties of the stokes flow in rectangular-shaped cavities has become a subject of primary interest in recent research.For the solutions of the problems of the stokes flow, there are three formulations, in terms of primitive variables of velocity-pressure and velocity, velocity-vorticity, and vorticity-potential forms. By using the primitive variables of velocity-pressure formulations, Zeb et al. [4], Alves and Silvestre [5], Chen et al. [6], and Tsai et al. [7] studied the two- and three-dimensional stokes flow problems. By adopting vorticity-velocity approach, a system of Laplace and Poisson type equations for the components of vorticity and velocity fields [8, 9], respectively can be transferred from the governing equations of primitive variable Stokes flow. The third formulations of the vorticity-potential formulations are used to deal with two-dimensional stokes flow problem, and they can be applied to derive the vorticity-stream functions, which are governed by the Laplace and Poisson equations or even only the biharmonic streamfunction equation. Gaskell solved the biharmonic equation of streamfunction to analyze Stokes flow in a double-lid-driven cavity with free surface side walls [10, 11], and he investigated the flow structure further by considering the behavior of the streamfunction close to stagnation points for two control parameters: the cavity aspect ratio and the speed ratio [12, 13]. Smyrlis and Karageorghis [14] and Onyejekwe [15] used the numerical methods for the biharmonic equation.However, each formulation will suffer different degrees of difficulty [16]. The velocity-pressure formulation requires complicated integrals to be evaluated and the use of the tensorial fundamental solution for the full stokes system [17], whilst the vorticity-velocity direct computation will have no immediate solutions for the pressure and needs explicitly the vorticity boundary conditions. The streamfunction-vorticity formulation is restricted to 2D and it is not straightforward to extend 2D streamfunction formulation to 3D flows, whilst the solution is based on the semi-inverse method which requires trial functions to satisfy the biharmonic equation of streamfunction and the boundary conditions.The present paper develops a symplectic analytical method to study the stokes flow problem in a double-lid-driven cavity with free surface side walls. Zhong first introduced the Hamiltonian formulation into the theory of elasticity and put forward a direct method [18]. Since the early 1990s, the researchers developed some problems in solid mechanics and elasticity [19–21]. Xu et al. used Hamiltonian systematic method to solve stokes flow problem in the rectangular domain with no-slip boundaries [22, 23] and axisymmetric problems of stokes flow in a tube [24]. Although Xu et al. solved for Stokes flow in a lid-driven rectangular cavity, they did not study the Stokes flow problem in a rectangular cavity with free surface side walls. In the present study we consider this problem through symplectic analytical method. Symplectic analytical method is based on Hamiltonian principle with Legengre’s transformation where by analytical solutional solutions could be obtained by expansion of eigenfunctions. The methodology is rational and systematic with clear step-by-step derivation procedure. The special properties under Hamiltonian system of this problem is that it’s eigenvalue-equations give the real repeated double roots hence existing first-order Jordan form.Two-dimensional variational principle for the Stokes flow in a rectangular domain is established. The Lagrange function is derived from the variational principle. Taking original variables (velocities) and their dual variables as the basic variables, Hamiltonian canonical equation of stokes flow using state variable is obtained and the symplectic solutions under free-surface conditions are also determined with the coordinatey simulated as time. The eigensolutions corresponding to zero-eigenvalue and nonzero-eigenvalues can be solved analytically, and the unknown expansion coefficients can be determined by using the end conditions based on the symplectic orthonormalization relationships.


## 2. Problem Specification and Hamiltonian System

### 2.1. Problem Specification

The flows [1] in the coating bead are modeled as double lid-driven cavity stokes flows which are bounded at the top and below by moving rigid roll surfaces and at the two sides by two free surfaces. The geometry properties of the rectangular are width b, depth a, the velocities of the free surfaces, and moving lids S
u and S
d as in Figure 1. There are some assumptions: (a) flow rate is negligibly small (based on experimental observation [25]); (b) the roll curvature can be neglected; (c) liquid-gas interfaces are planar; (d) gas adjacent to liquid remains at uniform pressure, its viscosity is neglected; (e) the fluid is Newtonian and incompressible with density ρ and dynamic viscosity G; (f) the flow is steady in the (x
,
y) plane with velocity q
=
(
u
,
v
); (g) the body forces and inertia effects are neglected.Figure 1
(a) Schematic of a roll coating system; (b) idealized model for the creeping flow, located between the rolls.The shear stress and velocityu on the two free surfaces are expressed as
(1)
τ
x
y
=
0
,
u
=
0
(
x
=
0
,
b
)
.For stokes flow in steady-state, the equilibrium equations can be written as(2)
∂
σ
x
∂
x
+
∂
τ
x
y
∂
y
=
0
,
∂
τ
x
y
∂
x
+
∂
σ
y
∂
y
=
0
,
where σ
x, σ
y and τ
x
y are the stress tensor components.The constitutive equations can be given by(3)
σ
x
=
-
p
f
+
2
G
ε
x
,
σ
y
=
-
p
f
+
2
G
ε
y
,
τ
x
y
=
τ
y
x
=
2
G
ε
x
y
,
where p
f is pressure, and ε
x, ε
y and ε
x
y are the strain tensor components.The strain-velocity relationships are denoted by(4)
ε
x
=
∂
u
∂
x
,
ε
y
=
∂
v
∂
y
,
ε
x
y
=
1
2
(
∂
u
∂
y
+
∂
v
∂
x
)
.The continuity equation is given by(5)
∂
u
∂
x
+
∂
v
∂
y
=
0
.The velocity components are expressed in terms of a stream functionΓ:
(6)
u
=
∂
Γ
∂
y
,
v
=
-
∂
Γ
∂
x
,
where the stream function is a constant (taken to be zero) on the boundaries.The force boundary conditions are given on the part of the surfaceS
f:
(7)
F
n
x
=
F
-
n
x
,
F
n
y
=
F
-
n
y
on
⁡
S
f
,
where F
-
n
x and F
-
n
y are the extra forces on the part of the surface S
f.On the other part of the boundaryS
u  (S
=
S
f
∪
S
u), velocities are defined as
(8)
u
=
u
-
,
v
=
v
-
on
⁡
S
u
,
where u
- and v
- are the velocities on the part of the surface S
u.


### 2.2. Symplectic Dual Equations

Multiplying (2) and (7) with virtual velocities field (δ
u
,
δ
v) and integrating the resultant equations over Ω and S
f, respectively, the following expression can be obtained as
(9)
-
∬
Ω
[
(
∂
σ
x
∂
x
+
∂
τ
y
x
∂
y
)
δ
u
+
(
∂
τ
y
x
∂
x
+
∂
σ
y
∂
y
)
δ
v
]
d
Ω
+
∫
S
f
[
(
F
n
x
-
F
-
n
x
)
δ
u
+
(
F
n
y
-
F
-
n
y
)
δ
v
]
d
S
=
0
,
where d
Ω
=
d
x
d
y (d
Ω is the area element of domain Ω) and d
S indicates the length element.The virtual velocities satisfy the geometric boundary conditions (u
=
u
-, v
=
v
- on S
u),
(10)
δ
u
=
0
,
δ
v
=
0
on
⁡
S
u
.We use the Gauss integral theorem to receive(11)
∬
Ω
(
σ
x
δ
ε
x
+
2
τ
x
y
δ
ε
x
y
+
σ
y
δ
ε
y
)
d
x
d
y
-
∫
S
f
(
F
-
n
x
δ
u
+
F
-
n
y
δ
v
)
d
S
=
∬
Ω
δ
(
G
ε
x
2
+
G
ε
y
2
+
2
G
ε
x
y
2
-
p
f
ε
x
-
p
f
ε
y
)
d
x
d
y
-
∫
S
f
(
F
-
n
x
δ
u
+
F
-
n
y
δ
v
)
d
S
=
δ
∬
Ω
(
G
ε
x
2
+
G
ε
y
2
+
2
G
ε
x
y
2
-
p
f
ε
x
-
p
f
ε
y
)
d
x
d
y
-
∫
S
f
(
F
-
n
x
δ
u
+
F
-
n
y
δ
v
)
d
S
=
0
,
where δ
ε
x
=
∂
δ
u
/
∂
x, δ
ε
y
=
∂
δ
v
/
∂
y, δ
ε
x
y
=
(
1
/
2
)
(
∂
δ
u
/
∂
y
+
∂
δ
v
/
∂
x
).Hamilton’s principle is(12)
δ
∬
Ω
L
d
x
d
y
=
0
.The Lagrange function of the stokes flow can be calculated by(13)
L
=
G
ε
x
2
+
G
ε
y
2
+
2
G
ε
x
y
2
-
p
f
ε
x
-
p
f
ε
y
=
G
[
(
∂
u
∂
x
)
2
+
(
∂
v
∂
y
)
2
]
+
1
2
G
(
∂
u
∂
y
+
∂
v
∂
x
)
2
-
p
f
(
∂
u
∂
x
+
∂
v
∂
y
)
.Introducing the dual variables(14)
p
1
=
∂
l
∂
u
˙
=
G
(
u
˙
+
∂
v
∂
x
)
,
p
2
=
∂
l
∂
v
˙
=
2
G
v
˙
-
p
f
,
the dot represents differential with respect to y; namely, (
·
)
=
∂
/
∂
y. The y coordinate is analogous to the time coordinate in the dynamic problem. The dual variables are stresses τ
x
y and σ
y. Write the vector forms of velocity variables and their dual variables as q
=
{
u
,
v
}
T and p
=
{
τ
x
y
,
σ
y
}
T
=
{
τ
,
σ
}
T, respectively. The Hamiltonian function H
(
q
,
p
) can be calculated as
(15)
Η
(
q
,
p
)
=
p
T
q
˙
-
L
(
q
,
q
˙
)
=
τ
u
˙
+
σ
v
˙
-
L
(
u
,
v
,
u
˙
,
v
˙
)
=
-
σ
∂
u
∂
x
-
τ
∂
v
∂
x
+
1
2
G
τ
2
-
2
G
(
∂
u
∂
x
)
2
.
Substituting (15) into (12) gives [23]:
(16)
q
˙
=
∂
H
∂
p
p
˙
=
-
∂
H
∂
q
.
Equation (16) is the Hamiltonian dual equations. Rewriting (16) in matrix-vector form, one has
(17)
w
˙
=
H
w
,
where w
=
(
u
,
v
,
τ
,
σ
)
T is the state vector for variables; the Hamiltonian operator matrix H is defined as
(18)
H
=
[
0
-
∂
∂
x
1
G
0
-
∂
∂
x
0
0
0
-
4
G
∂
2
∂
x
2
0
0
-
∂
∂
x
0
0
-
∂
∂
x
0
]
.


### 2.3. Eigenvalue Solutions

Applying the method of separation of variables tow yields
(19)
w
(
y
,
x
)
=
f
(
y
)
ψ
(
x
)
.
Substituting the above expression into (17) gives
(20)
f
(
y
)
=
e
μ
y
and the eigenvalue equation
(21)
H
ψ
(
x
)
=
μ
ψ
(
x
)
,
where μ is the eigenvalue and Ψ
(
x
) is the corresponding eigenvector.As to zero eigenvalueμ
=
0, (21) gives
(22)
H
ψ
(
x
)
=
0
.
Under lateral boundary condition (1), there does not exist the zero eigenvalue solution.The eigensolutions of nonzero eigenvalues in (21) may be obtained by expanding the eigenvalue equation. First, the eigenvalues α in the x-direction can be obtained by substituting
(23)
u
=
e
α
x
,
v
=
e
α
x
,
τ
=
e
α
x
,
σ
=
e
α
x
into (21) which yields a characteristic polynomial as follows:
(24)
det
⁡
(
-
μ
-
α
1
G
0
-
α
-
μ
0
0
-
4
G
α
2
0
-
μ
-
α
0
0
-
α
-
μ
)
=
0
with repeated roots α
=
±
μ
i as the eigenvalues. Hence, the general solutions of nonzero eigenvalues are
(25)
u
=
A
u
sin
(
μ
x
)
+
B
u
cos
⁡
(
μ
x
)
+
C
u
x
cos
⁡
(
μ
x
)
+
D
u
x
sin
(
μ
x
)
v
=
A
v
cos
⁡
(
μ
x
)
+
B
v
sin
(
μ
x
)
+
C
v
x
sin
(
μ
x
)
+
D
v
x
cos
⁡
(
μ
x
)
τ
=
A
τ
sin
(
μ
x
)
+
B
τ
cos
⁡
(
μ
x
)
+
C
τ
x
cos
⁡
(
μ
x
)
+
D
τ
x
sin
(
μ
x
)
σ
=
A
σ
cos
⁡
(
μ
x
)
+
B
σ
sin
(
μ
x
)
+
C
σ
x
sin
(
μ
x
)
+
D
σ
x
cos
⁡
(
μ
x
)
.The constants are not all independent. For convenience,A
u, B
u, C
u, and D
u may be chosen as the independent constants. Substituting (25) into (21) yields the relations between these constants. Further, substituting general solution (25) into the corresponding boundary conditions on both sides x
=
0 or b, then
(26)
si
n
2
(
μ
b
)
=
0
which gives real repeated double roots as
(27)
μ
n
=
n
π
b
(
n
=
±
1
,
±
2
,
…
)
.
The corresponding basic eigenvector is
(28)
ψ
n
(
0
)
=
{
u
,
v
,
τ
x
y
,
σ
y
}
T
=
{
sin
(
μ
n
x
)
,
-
cos
⁡
(
μ
n
x
)
,
2
G
μ
n
sin
(
μ
n
x
)
,
-
2
G
μ
n
cos
⁡
(
μ
n
x
)
}
T
.
Then, the solution to eigenvalue equation (17) is
(29)
w
n
(
0
)
=
e
μ
n
y
ψ
n
(
0
)
.
From velocity-streamfunction relation (6), the streamfunction Γ can be expressed as
(30)
Γ
n
(
0
)
=
e
μ
n
y
μ
n
sin
(
μ
n
x
)
.
Because the eigenvalue μ
n is a double root, the first-order Jordan form eigen-solution can be solved via
(31)
H
ψ
n
(
1
)
=
μ
n
ψ
n
(
1
)
+
ψ
n
(
0
)
(
n
=
±
1
,
±
2
,
…
)
.
Imposing boundary conditions (1) yields
(32)
ψ
n
(
1
)
=
{
1
μ
n
sin
(
μ
n
x
)
,
0,2
G
sin
(
μ
n
x
)
,
0
}
T
.
Hence, the solution to (17) is
(33)
w
n
(
1
)
=
e
μ
n
y
(
ψ
n
(
1
)
+
y
ψ
n
(
0
)
)
.
Again, from velocity-streamfunction relation (6), the streamfunction Γ can be expressed as
(34)
Γ
n
(
1
)
=
y
μ
n
e
μ
n
y
sin
(
μ
n
x
)
.From the eigenvalues and eigenvectors with adjoint symplectic orthogonality property, the general solution for stokes flow with free surface on lateral boundaries can be expressed as(35)
w
=
∑
n
=
1
∞
[
f
n
(
0
)
w
n
(
0
)
+
f
n
(
1
)
w
n
(
1
)
+
f
-
n
(
0
)
w
-
n
(
0
)
+
f
-
n
(
1
)
w
-
n
(
1
)
]
.
Equation (35) above strictly satisfies the homogeneous differential equation in the domain and homogeneous boundary conditions (1), while f
n
(
k
)
(
k
=
0,1
;
n
=
±
1
,
±
2
,
…
) are unknown constants which can be determined by imposing the remaining two boundary conditions at y
=
a
/
2 and y
=
-
a
/
2. According to the expansion theorem, the streamfunction Γ can be expressed as
(36)
Γ
=
∑
n
=
1
∞
[
f
n
(
0
)
Γ
n
(
0
)
+
f
n
(
1
)
Γ
n
(
1
)
+
f
-
n
(
0
)
Γ
-
n
(
0
)
+
f
-
n
(
1
)
Γ
-
n
(
1
)
]
.


### 2.4. End Conditions

The end boundary conditions are(37)
u
|
y
=
a
/
2
=
S
u
,
u
|
y
=
-
a
/
2
=
S
d
,
v
|
y
=
±
a
/
2
=
0
.The expression ofu above can be represented in Fourier series as [19]:
(38)
u
=
S
=
∑
n
=
1
∞
2
sin
⁡
(
n
π
x
/
b
)
b
∫
0
b
S
sin
n
π
x
b
d
x
=
-
4
S
π
∑
n
=
1,3
,
5
,
…
∞
1
n
sin
n
π
x
b
.
However, from (35), one has
(39)
u
=
∑
n
=
1
∞
{
f
n
(
0
)
e
μ
n
y
+
f
n
(
1
)
e
μ
n
y
(
1
μ
n
+
y
)
-
f
-
n
(
0
)
e
-
μ
n
y
-
f
-
n
(
1
)
e
-
μ
n
y
(
-
1
μ
n
+
y
)
}
sin
⁡
μ
n
x
v
=
∑
n
=
1
∞
[
-
f
n
(
0
)
e
μ
n
y
-
f
n
(
1
)
e
μ
n
y
y
-
f
-
n
(
0
)
e
-
μ
n
y
-
f
-
n
(
1
)
e
-
μ
n
y
y
]
sin
⁡
μ
n
x
.
Substituting y
=
±
a
/
2 into (39) for the left-hand-side of (37) and using the Fourier series representations of u in (38) on the right-hand-side, four set of equations can be derived. The constants f
n
(
0
), f
n
(
1
), f
-
n
(
0
), andf
n
(
0
) can be solved by comparing the coefficients of sin(μ
n
x) and cos(μ
n
x), which are
(40)
f
n
(
0
)
=
f
-
n
(
0
)
=
f
n
(
1
)
=
f
-
n
(
1
)
for
n
=
2,4
,
6
,
…
,
f
n
(
0
)
=
-
[
2
e
2
α
n
s
h
(
2
α
n
)
+
4
α
n
]
a
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
-
[
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
a
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
f
-
n
(
0
)
=
[
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
a
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
+
[
2
e
2
α
n
s
h
(
2
α
n
)
+
4
α
n
]
a
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
f
n
(
1
)
=
2
[
2
e
2
α
n
s
h
(
2
α
n
)
-
4
α
n
]
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
+
2
[
-
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
f
-
n
(
1
)
=
2
[
-
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
+
2
[
2
e
2
α
n
s
h
(
2
α
n
)
-
4
α
n
]
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
where
(41)
α
n
=
λ
n
b
2
(
n
=
1,3
,
5
,
…
)
.
Substituting (40)~(41) into (35)~(36), the general solutions for velocity components, stress components, and streamfunction can be derived. In practice, it is only possible to include a finite number (
N
) of terms in the series (35), it is necessary to first establish that the series converges and then determine the number of terms that need to be taken into account in order to ensure the convergence of the series was satisfactorily.


## 2.1. Problem Specification

The flows [1] in the coating bead are modeled as double lid-driven cavity stokes flows which are bounded at the top and below by moving rigid roll surfaces and at the two sides by two free surfaces. The geometry properties of the rectangular are width b, depth a, the velocities of the free surfaces, and moving lids S
u and S
d as in Figure 1. There are some assumptions: (a) flow rate is negligibly small (based on experimental observation [25]); (b) the roll curvature can be neglected; (c) liquid-gas interfaces are planar; (d) gas adjacent to liquid remains at uniform pressure, its viscosity is neglected; (e) the fluid is Newtonian and incompressible with density ρ and dynamic viscosity G; (f) the flow is steady in the (x
,
y) plane with velocity q
=
(
u
,
v
); (g) the body forces and inertia effects are neglected.Figure 1
(a) Schematic of a roll coating system; (b) idealized model for the creeping flow, located between the rolls.The shear stress and velocityu on the two free surfaces are expressed as
(1)
τ
x
y
=
0
,
u
=
0
(
x
=
0
,
b
)
.For stokes flow in steady-state, the equilibrium equations can be written as(2)
∂
σ
x
∂
x
+
∂
τ
x
y
∂
y
=
0
,
∂
τ
x
y
∂
x
+
∂
σ
y
∂
y
=
0
,
where σ
x, σ
y and τ
x
y are the stress tensor components.The constitutive equations can be given by(3)
σ
x
=
-
p
f
+
2
G
ε
x
,
σ
y
=
-
p
f
+
2
G
ε
y
,
τ
x
y
=
τ
y
x
=
2
G
ε
x
y
,
where p
f is pressure, and ε
x, ε
y and ε
x
y are the strain tensor components.The strain-velocity relationships are denoted by(4)
ε
x
=
∂
u
∂
x
,
ε
y
=
∂
v
∂
y
,
ε
x
y
=
1
2
(
∂
u
∂
y
+
∂
v
∂
x
)
.The continuity equation is given by(5)
∂
u
∂
x
+
∂
v
∂
y
=
0
.The velocity components are expressed in terms of a stream functionΓ:
(6)
u
=
∂
Γ
∂
y
,
v
=
-
∂
Γ
∂
x
,
where the stream function is a constant (taken to be zero) on the boundaries.The force boundary conditions are given on the part of the surfaceS
f:
(7)
F
n
x
=
F
-
n
x
,
F
n
y
=
F
-
n
y
on
⁡
S
f
,
where F
-
n
x and F
-
n
y are the extra forces on the part of the surface S
f.On the other part of the boundaryS
u  (S
=
S
f
∪
S
u), velocities are defined as
(8)
u
=
u
-
,
v
=
v
-
on
⁡
S
u
,
where u
- and v
- are the velocities on the part of the surface S
u.


## 2.2. Symplectic Dual Equations

Multiplying (2) and (7) with virtual velocities field (δ
u
,
δ
v) and integrating the resultant equations over Ω and S
f, respectively, the following expression can be obtained as
(9)
-
∬
Ω
[
(
∂
σ
x
∂
x
+
∂
τ
y
x
∂
y
)
δ
u
+
(
∂
τ
y
x
∂
x
+
∂
σ
y
∂
y
)
δ
v
]
d
Ω
+
∫
S
f
[
(
F
n
x
-
F
-
n
x
)
δ
u
+
(
F
n
y
-
F
-
n
y
)
δ
v
]
d
S
=
0
,
where d
Ω
=
d
x
d
y (d
Ω is the area element of domain Ω) and d
S indicates the length element.The virtual velocities satisfy the geometric boundary conditions (u
=
u
-, v
=
v
- on S
u),
(10)
δ
u
=
0
,
δ
v
=
0
on
⁡
S
u
.We use the Gauss integral theorem to receive(11)
∬
Ω
(
σ
x
δ
ε
x
+
2
τ
x
y
δ
ε
x
y
+
σ
y
δ
ε
y
)
d
x
d
y
-
∫
S
f
(
F
-
n
x
δ
u
+
F
-
n
y
δ
v
)
d
S
=
∬
Ω
δ
(
G
ε
x
2
+
G
ε
y
2
+
2
G
ε
x
y
2
-
p
f
ε
x
-
p
f
ε
y
)
d
x
d
y
-
∫
S
f
(
F
-
n
x
δ
u
+
F
-
n
y
δ
v
)
d
S
=
δ
∬
Ω
(
G
ε
x
2
+
G
ε
y
2
+
2
G
ε
x
y
2
-
p
f
ε
x
-
p
f
ε
y
)
d
x
d
y
-
∫
S
f
(
F
-
n
x
δ
u
+
F
-
n
y
δ
v
)
d
S
=
0
,
where δ
ε
x
=
∂
δ
u
/
∂
x, δ
ε
y
=
∂
δ
v
/
∂
y, δ
ε
x
y
=
(
1
/
2
)
(
∂
δ
u
/
∂
y
+
∂
δ
v
/
∂
x
).Hamilton’s principle is(12)
δ
∬
Ω
L
d
x
d
y
=
0
.The Lagrange function of the stokes flow can be calculated by(13)
L
=
G
ε
x
2
+
G
ε
y
2
+
2
G
ε
x
y
2
-
p
f
ε
x
-
p
f
ε
y
=
G
[
(
∂
u
∂
x
)
2
+
(
∂
v
∂
y
)
2
]
+
1
2
G
(
∂
u
∂
y
+
∂
v
∂
x
)
2
-
p
f
(
∂
u
∂
x
+
∂
v
∂
y
)
.Introducing the dual variables(14)
p
1
=
∂
l
∂
u
˙
=
G
(
u
˙
+
∂
v
∂
x
)
,
p
2
=
∂
l
∂
v
˙
=
2
G
v
˙
-
p
f
,
the dot represents differential with respect to y; namely, (
·
)
=
∂
/
∂
y. The y coordinate is analogous to the time coordinate in the dynamic problem. The dual variables are stresses τ
x
y and σ
y. Write the vector forms of velocity variables and their dual variables as q
=
{
u
,
v
}
T and p
=
{
τ
x
y
,
σ
y
}
T
=
{
τ
,
σ
}
T, respectively. The Hamiltonian function H
(
q
,
p
) can be calculated as
(15)
Η
(
q
,
p
)
=
p
T
q
˙
-
L
(
q
,
q
˙
)
=
τ
u
˙
+
σ
v
˙
-
L
(
u
,
v
,
u
˙
,
v
˙
)
=
-
σ
∂
u
∂
x
-
τ
∂
v
∂
x
+
1
2
G
τ
2
-
2
G
(
∂
u
∂
x
)
2
.
Substituting (15) into (12) gives [23]:
(16)
q
˙
=
∂
H
∂
p
p
˙
=
-
∂
H
∂
q
.
Equation (16) is the Hamiltonian dual equations. Rewriting (16) in matrix-vector form, one has
(17)
w
˙
=
H
w
,
where w
=
(
u
,
v
,
τ
,
σ
)
T is the state vector for variables; the Hamiltonian operator matrix H is defined as
(18)
H
=
[
0
-
∂
∂
x
1
G
0
-
∂
∂
x
0
0
0
-
4
G
∂
2
∂
x
2
0
0
-
∂
∂
x
0
0
-
∂
∂
x
0
]
.


## 2.3. Eigenvalue Solutions

Applying the method of separation of variables tow yields
(19)
w
(
y
,
x
)
=
f
(
y
)
ψ
(
x
)
.
Substituting the above expression into (17) gives
(20)
f
(
y
)
=
e
μ
y
and the eigenvalue equation
(21)
H
ψ
(
x
)
=
μ
ψ
(
x
)
,
where μ is the eigenvalue and Ψ
(
x
) is the corresponding eigenvector.As to zero eigenvalueμ
=
0, (21) gives
(22)
H
ψ
(
x
)
=
0
.
Under lateral boundary condition (1), there does not exist the zero eigenvalue solution.The eigensolutions of nonzero eigenvalues in (21) may be obtained by expanding the eigenvalue equation. First, the eigenvalues α in the x-direction can be obtained by substituting
(23)
u
=
e
α
x
,
v
=
e
α
x
,
τ
=
e
α
x
,
σ
=
e
α
x
into (21) which yields a characteristic polynomial as follows:
(24)
det
⁡
(
-
μ
-
α
1
G
0
-
α
-
μ
0
0
-
4
G
α
2
0
-
μ
-
α
0
0
-
α
-
μ
)
=
0
with repeated roots α
=
±
μ
i as the eigenvalues. Hence, the general solutions of nonzero eigenvalues are
(25)
u
=
A
u
sin
(
μ
x
)
+
B
u
cos
⁡
(
μ
x
)
+
C
u
x
cos
⁡
(
μ
x
)
+
D
u
x
sin
(
μ
x
)
v
=
A
v
cos
⁡
(
μ
x
)
+
B
v
sin
(
μ
x
)
+
C
v
x
sin
(
μ
x
)
+
D
v
x
cos
⁡
(
μ
x
)
τ
=
A
τ
sin
(
μ
x
)
+
B
τ
cos
⁡
(
μ
x
)
+
C
τ
x
cos
⁡
(
μ
x
)
+
D
τ
x
sin
(
μ
x
)
σ
=
A
σ
cos
⁡
(
μ
x
)
+
B
σ
sin
(
μ
x
)
+
C
σ
x
sin
(
μ
x
)
+
D
σ
x
cos
⁡
(
μ
x
)
.The constants are not all independent. For convenience,A
u, B
u, C
u, and D
u may be chosen as the independent constants. Substituting (25) into (21) yields the relations between these constants. Further, substituting general solution (25) into the corresponding boundary conditions on both sides x
=
0 or b, then
(26)
si
n
2
(
μ
b
)
=
0
which gives real repeated double roots as
(27)
μ
n
=
n
π
b
(
n
=
±
1
,
±
2
,
…
)
.
The corresponding basic eigenvector is
(28)
ψ
n
(
0
)
=
{
u
,
v
,
τ
x
y
,
σ
y
}
T
=
{
sin
(
μ
n
x
)
,
-
cos
⁡
(
μ
n
x
)
,
2
G
μ
n
sin
(
μ
n
x
)
,
-
2
G
μ
n
cos
⁡
(
μ
n
x
)
}
T
.
Then, the solution to eigenvalue equation (17) is
(29)
w
n
(
0
)
=
e
μ
n
y
ψ
n
(
0
)
.
From velocity-streamfunction relation (6), the streamfunction Γ can be expressed as
(30)
Γ
n
(
0
)
=
e
μ
n
y
μ
n
sin
(
μ
n
x
)
.
Because the eigenvalue μ
n is a double root, the first-order Jordan form eigen-solution can be solved via
(31)
H
ψ
n
(
1
)
=
μ
n
ψ
n
(
1
)
+
ψ
n
(
0
)
(
n
=
±
1
,
±
2
,
…
)
.
Imposing boundary conditions (1) yields
(32)
ψ
n
(
1
)
=
{
1
μ
n
sin
(
μ
n
x
)
,
0,2
G
sin
(
μ
n
x
)
,
0
}
T
.
Hence, the solution to (17) is
(33)
w
n
(
1
)
=
e
μ
n
y
(
ψ
n
(
1
)
+
y
ψ
n
(
0
)
)
.
Again, from velocity-streamfunction relation (6), the streamfunction Γ can be expressed as
(34)
Γ
n
(
1
)
=
y
μ
n
e
μ
n
y
sin
(
μ
n
x
)
.From the eigenvalues and eigenvectors with adjoint symplectic orthogonality property, the general solution for stokes flow with free surface on lateral boundaries can be expressed as(35)
w
=
∑
n
=
1
∞
[
f
n
(
0
)
w
n
(
0
)
+
f
n
(
1
)
w
n
(
1
)
+
f
-
n
(
0
)
w
-
n
(
0
)
+
f
-
n
(
1
)
w
-
n
(
1
)
]
.
Equation (35) above strictly satisfies the homogeneous differential equation in the domain and homogeneous boundary conditions (1), while f
n
(
k
)
(
k
=
0,1
;
n
=
±
1
,
±
2
,
…
) are unknown constants which can be determined by imposing the remaining two boundary conditions at y
=
a
/
2 and y
=
-
a
/
2. According to the expansion theorem, the streamfunction Γ can be expressed as
(36)
Γ
=
∑
n
=
1
∞
[
f
n
(
0
)
Γ
n
(
0
)
+
f
n
(
1
)
Γ
n
(
1
)
+
f
-
n
(
0
)
Γ
-
n
(
0
)
+
f
-
n
(
1
)
Γ
-
n
(
1
)
]
.


## 2.4. End Conditions

The end boundary conditions are(37)
u
|
y
=
a
/
2
=
S
u
,
u
|
y
=
-
a
/
2
=
S
d
,
v
|
y
=
±
a
/
2
=
0
.The expression ofu above can be represented in Fourier series as [19]:
(38)
u
=
S
=
∑
n
=
1
∞
2
sin
⁡
(
n
π
x
/
b
)
b
∫
0
b
S
sin
n
π
x
b
d
x
=
-
4
S
π
∑
n
=
1,3
,
5
,
…
∞
1
n
sin
n
π
x
b
.
However, from (35), one has
(39)
u
=
∑
n
=
1
∞
{
f
n
(
0
)
e
μ
n
y
+
f
n
(
1
)
e
μ
n
y
(
1
μ
n
+
y
)
-
f
-
n
(
0
)
e
-
μ
n
y
-
f
-
n
(
1
)
e
-
μ
n
y
(
-
1
μ
n
+
y
)
}
sin
⁡
μ
n
x
v
=
∑
n
=
1
∞
[
-
f
n
(
0
)
e
μ
n
y
-
f
n
(
1
)
e
μ
n
y
y
-
f
-
n
(
0
)
e
-
μ
n
y
-
f
-
n
(
1
)
e
-
μ
n
y
y
]
sin
⁡
μ
n
x
.
Substituting y
=
±
a
/
2 into (39) for the left-hand-side of (37) and using the Fourier series representations of u in (38) on the right-hand-side, four set of equations can be derived. The constants f
n
(
0
), f
n
(
1
), f
-
n
(
0
), andf
n
(
0
) can be solved by comparing the coefficients of sin(μ
n
x) and cos(μ
n
x), which are
(40)
f
n
(
0
)
=
f
-
n
(
0
)
=
f
n
(
1
)
=
f
-
n
(
1
)
for
n
=
2,4
,
6
,
…
,
f
n
(
0
)
=
-
[
2
e
2
α
n
s
h
(
2
α
n
)
+
4
α
n
]
a
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
-
[
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
a
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
f
-
n
(
0
)
=
[
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
a
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
+
[
2
e
2
α
n
s
h
(
2
α
n
)
+
4
α
n
]
a
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
f
n
(
1
)
=
2
[
2
e
2
α
n
s
h
(
2
α
n
)
-
4
α
n
]
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
+
2
[
-
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
f
-
n
(
1
)
=
2
[
-
2
s
h
(
2
α
n
)
+
4
α
n
e
2
α
n
]
S
u
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
+
2
[
2
e
2
α
n
s
h
(
2
α
n
)
-
4
α
n
]
S
d
b
e
α
n
[
c
h
(
4
α
n
)
-
1
-
8
α
n
2
]
for
n
=
1,3
,
5
,
…
,
where
(41)
α
n
=
λ
n
b
2
(
n
=
1,3
,
5
,
…
)
.
Substituting (40)~(41) into (35)~(36), the general solutions for velocity components, stress components, and streamfunction can be derived. In practice, it is only possible to include a finite number (
N
) of terms in the series (35), it is necessary to first establish that the series converges and then determine the number of terms that need to be taken into account in order to ensure the convergence of the series was satisfactorily.


## 3. Numerical Results

### 3.1. Verification of the Present Analytical Solutions

In order to validate the present formulation and the computing program developed by the author, the comparison is made to illustrate the convergence of the solution. The comparison results are shown in Table1, where the cavity aspect ratio A
=
a
/
b and the speed ratio S
=
S
u
/
S
d are taken 0.5 and 2. In Table 1, those data with superscripts * are calculated from [10]. By comparing results in Table 1, it can be obviously found that the analytical solutions in the present study are in good agreement with those in [10]. From Table 1, the average deviation from the lid velocities is approximately 3%, when the number of eigen-solution is equal to N
=
50. Meanwhile, it can be easily found that as N increases to 200, the average deviation decreases approximately to 0.3%, and the only places where we have any significant deviation are at the extremities where there is clearly a manifestation of the Gibbs phenomenon. The discontinuity in velocity will result in a Gibbs phenomenon at the point of contact between the stationary and moving boundaries but will not cause any significant numerical difficulties with the results we have in Section 3.2. The convergence of the present method has been proved in [26].Table 1
Comparison calculated results at the two moving end walls with that in Reference [10] (superscript*); A
=
0.5, S
=
2.
x
N
y
=
a
/
2
y
=
-
a
/
2
u
v
Γ
u
v
0.0
200
0.0
1.9968*
1.3323 × 10−15
0.0
0.0
0.9984*
−1.0333 × 10−16
50
0.0
1.9873*
1.3323 × 10−15
0.0
0.0
0.9936*
−1.0333 × 10−16
0.2
200
1.9892
1.9967*
−2.2204 × 10−16
1.9530 × 10−16
0.9946
0.9983*
−1.9723 × 10−15
50
1.9568
1.9866*
−2.2204 × 10−16
1.8865 × 10−16
0.9784
0.9933*
−1.9723 × 10−15
0.4
200
1.9933
1.9961*
2.2204 × 10−16
1.8596 × 10−16
0.9967
0.9980*
−7.0408 × 10−16
50
1.9732
1.9843*
2.2204 × 10−16
1.7981 × 10−16
0.9866
0.9921*
−7.0408 × 10−16
0.6
200
1.9933
1.9946*
−2.2204 × 10−16
1.8596 × 10−16
0.9967
0.9973*
7.0408 × 10−16
50
1.9732
1.9783*
−2.2204 × 10−16
1.7981 × 10−16
0.9866
0.9892*
7.0408 × 10−16
0.8
200
1.9892
1.9897*
2.2204 × 10−16
1.9530 × 10−16
0.9946
0.9949*
1.9723 × 10−15
50
1.9568
1.9589*
2.2204 × 10−16
1.8865 × 10−16
0.9784
0.9794*
1.9723 × 10−15
1.0
200
3.2826 × 10−14
3.2244 × 10−14*
−1.3323 × 10−15
0.0
1.6413 × 10−14
1.6122 × 10−14*
1.0333 × 10−16
50
9.0049 × 10−15
8.0239 × 10−15*
−1.3323 × 10−15
0.0
4.5024 × 10−15
4.0119 × 10−15*
1.0333 × 10−16


### 3.2. Stagnation Points and Flow Patterns

The stagnation points in the interior of the cavity, where the velocity is equal to zero, are important to study the flow structure. A full search of stagnation points in the cavity for any value of the top and bottom velocities can be done by analysing the contour levels of the streamfunction and by looking for its local extremum. Besides, by symmetry conditions, the position of a stagnation point (0.5
b
,
Yst) on x
=
0.5
b can be found by solving equation u
(
0.5
b
,
Y
st
)
=
0, because the component v is equal to zero on this line. The stagnation points can be divided into elliptic stagnation points, where Γ has an extremum value (maximum or minimum) and nearby streamlines surrounding it have the form of ellipses or hyperbolic (saddle) stagnation points, where Γ has a minimax value and nearby streamlines form two sets of hyperbolae which are locally not closed [27]. Gaskell et al. [11] identified the flow bifurcations at the stagnation points as control parameter A is continuously varied. Here, we find the position of stagnation points numerically as control parameter S is continuously varied from −40 to 40 with A
=
1.0.We start with the problem that the flow in a square cavity,A
=
1.0, is generated by moving the top and bottom wall in same directions with the same velocity. The flow structure simply consists of two eddies of equal size (Figure 2(a)). The steady streamline patterns for S
=
1.0 reveal closed lines around elliptical stagnation points with Cartesian coordinate (0.5
b, 0.2938
a) and (0.5
b, -
0.2938
a), respectively. These two stagnation points are equivalent to the zero points of velocity u at line x
=
0.5
b. It is clear that the streamfunction Γ, the velocity v, and the shear stress τ are antisymmetric about the horizontal centerline, whilst the velocity u and the normal stress σ are symmetric about this line. From Figure 2(b), we can see two zero-contour lines, which indicate that there are two zero-value points of velocity u for every section paralleling with line x
=
0.5
b. The velocity v in the domain is divided into four zones and there are also two zero-contour lines which correspond to the horizontal and vertical centerline. In the vicinity of the corners, there are stress concentration zones (Figures 3(a) and 3(b)) due to the discontinuity in velocity u. For S
=
2, two separate eddy structures (Figure 3(c)) still exist, the larger of which is that adjacent to the faster moving wall, and the stagnation points in Cartesian coordinate are (0.5
b, 0.2581
a) and (0.5
b, 0.3411
a). Figure 4(a) illustrates that, for S
>
0, there are two large eddies within the cavity with two stagnation points on the centre-line between the free surfaces. The effect on the position of the eddies for S
∈
[
0,40
] is shown in Figure 4(b). As S increases towards infinity, the stagnation point in the upper cavity moves towards the horizontal centerline and reaches to the point (0.5
b, 0.2069
a), while the stagnation point in the under cavity moves towards the bottom lid and reaches to (0.5
b, -
0.5
a) and the eddy in the under cavity reduces and disappears.Contour lines withS
=
1: (a) streamfunction Γ; (b) velocity u; (c) velocity v.
(a)
(b)
(c)Contour lines: (a) shear stressτ (S
=
0); (b) normal stress σ (S
=
0); (c) streamfunction Γ (S
=
2).
(a)
(b)
(c)(a) (u
,
S) diagram for S
>
0 at x
=
0.5
b; (b) (Y
st
,
S) diagram for S
>
0 with A
=
1.0.
(a)
(b)Consider, then, the case of flowS
<
0 with A
=
1, including the case S
=
0, which is generated by moving the top and bottom wall in opposite directions. As for flow patterns when S
=
0, the cavity consists only of a large recirculation. The streamline patterns for S
=
0 reveal closed lines around elliptical stagnation points with Cartesian coordinate (0.5
b, -
0.2069
a). We can note that the case S
=
0 is essentially the same as S
=
∞; the immobile and moving lids are merely interchanged. From Figure 5(b), we can see one zero-contour line, which indicates that there are one zero-value point of velocity u for every section parallel to line x
=
0.5
b. The velocity v in the domain is divided into two zones and there is also one zero-contour line which corresponds to the vertical centerline. In the vicinity of the bottom corners, there are stress concentration zones (Figure 6) due to the discontinuity in velocity u. For S
=
-
1, the streamlines in Figure 7(a) show a single eddy with its centre at the middle of the cavity. From Figure 7(b), we can see one zero-contour line which just corresponds to line y
=
0. The velocity v is antisymmetric about x
=
0.5
b and there is also one zero-contour line which corresponds to the vertical centerline. The shear stress τ is symmetric about both the horizontal and the vertical centerline and there are also two zero-contour lines, whilst the normal stress σ is symmetric about the center of the cavity and there are four zero-contour lines. Figure 8(c) shows the flow pattern for S
=
-
2, where it is observed that (compared with S
=
-
1) the eddy centre is displaced from (0.5
b, 0) along the line of x
=
0.5
b—towards the faster moving lid—and thus symmetry about y
=
0 is lost, and the stagnation points in Cartesian coordinate are (0.5
b, 0.1256
a). Figure 9(a) illustrates that, for S
≤
0, there is one eddy within the cavity with one stagnation point on the centre-line between the free surfaces. The effect on the position of the eddies for S
∈
[
-
40
,
0
] is shown in Figure 9(b). As S is reduced towards -
∞, the stagnation point moves towards the top lid and reaches the point (0.5
b, 0.2069
a). We note that the case S
=
-
∞ is essentially the same as the case S
=
∞, the direction of swirl is merely interchanged.Contour lines withS
=
0: (a) streamfunction Γ; (b) velocity u; (c) velocity v.
(a)
(b)
(c)Contour lines withS
=
0: (a) shear stress τ; (b) normal stress σ.
(a)
(b)Contour lines withS
=
-
1: (a) streamfunction Γ; (b) velocity u; (c) velocity v.
(a)
(b)
(c)Contour lines: (a) shear stressτ (S
=
-
1); (b) normal stress σ (S
=
-
1); (c) streamfunction Γ (S
=
-
2).
(a)
(b)
(c)(a) (u
,
S) diagram with S
≤
0 at x
=
0.5
b; (b) (Y
st
,
S) diagram with S
≤
0 for A
=
1.0.
(a)
(b)


## 3.1. Verification of the Present Analytical Solutions

In order to validate the present formulation and the computing program developed by the author, the comparison is made to illustrate the convergence of the solution. The comparison results are shown in Table1, where the cavity aspect ratio A
=
a
/
b and the speed ratio S
=
S
u
/
S
d are taken 0.5 and 2. In Table 1, those data with superscripts * are calculated from [10]. By comparing results in Table 1, it can be obviously found that the analytical solutions in the present study are in good agreement with those in [10]. From Table 1, the average deviation from the lid velocities is approximately 3%, when the number of eigen-solution is equal to N
=
50. Meanwhile, it can be easily found that as N increases to 200, the average deviation decreases approximately to 0.3%, and the only places where we have any significant deviation are at the extremities where there is clearly a manifestation of the Gibbs phenomenon. The discontinuity in velocity will result in a Gibbs phenomenon at the point of contact between the stationary and moving boundaries but will not cause any significant numerical difficulties with the results we have in Section 3.2. The convergence of the present method has been proved in [26].Table 1
Comparison calculated results at the two moving end walls with that in Reference [10] (superscript*); A
=
0.5, S
=
2.
x
N
y
=
a
/
2
y
=
-
a
/
2
u
v
Γ
u
v
0.0
200
0.0
1.9968*
1.3323 × 10−15
0.0
0.0
0.9984*
−1.0333 × 10−16
50
0.0
1.9873*
1.3323 × 10−15
0.0
0.0
0.9936*
−1.0333 × 10−16
0.2
200
1.9892
1.9967*
−2.2204 × 10−16
1.9530 × 10−16
0.9946
0.9983*
−1.9723 × 10−15
50
1.9568
1.9866*
−2.2204 × 10−16
1.8865 × 10−16
0.9784
0.9933*
−1.9723 × 10−15
0.4
200
1.9933
1.9961*
2.2204 × 10−16
1.8596 × 10−16
0.9967
0.9980*
−7.0408 × 10−16
50
1.9732
1.9843*
2.2204 × 10−16
1.7981 × 10−16
0.9866
0.9921*
−7.0408 × 10−16
0.6
200
1.9933
1.9946*
−2.2204 × 10−16
1.8596 × 10−16
0.9967
0.9973*
7.0408 × 10−16
50
1.9732
1.9783*
−2.2204 × 10−16
1.7981 × 10−16
0.9866
0.9892*
7.0408 × 10−16
0.8
200
1.9892
1.9897*
2.2204 × 10−16
1.9530 × 10−16
0.9946
0.9949*
1.9723 × 10−15
50
1.9568
1.9589*
2.2204 × 10−16
1.8865 × 10−16
0.9784
0.9794*
1.9723 × 10−15
1.0
200
3.2826 × 10−14
3.2244 × 10−14*
−1.3323 × 10−15
0.0
1.6413 × 10−14
1.6122 × 10−14*
1.0333 × 10−16
50
9.0049 × 10−15
8.0239 × 10−15*
−1.3323 × 10−15
0.0
4.5024 × 10−15
4.0119 × 10−15*
1.0333 × 10−16


## 3.2. Stagnation Points and Flow Patterns

The stagnation points in the interior of the cavity, where the velocity is equal to zero, are important to study the flow structure. A full search of stagnation points in the cavity for any value of the top and bottom velocities can be done by analysing the contour levels of the streamfunction and by looking for its local extremum. Besides, by symmetry conditions, the position of a stagnation point (0.5
b
,
Yst) on x
=
0.5
b can be found by solving equation u
(
0.5
b
,
Y
st
)
=
0, because the component v is equal to zero on this line. The stagnation points can be divided into elliptic stagnation points, where Γ has an extremum value (maximum or minimum) and nearby streamlines surrounding it have the form of ellipses or hyperbolic (saddle) stagnation points, where Γ has a minimax value and nearby streamlines form two sets of hyperbolae which are locally not closed [27]. Gaskell et al. [11] identified the flow bifurcations at the stagnation points as control parameter A is continuously varied. Here, we find the position of stagnation points numerically as control parameter S is continuously varied from −40 to 40 with A
=
1.0.We start with the problem that the flow in a square cavity,A
=
1.0, is generated by moving the top and bottom wall in same directions with the same velocity. The flow structure simply consists of two eddies of equal size (Figure 2(a)). The steady streamline patterns for S
=
1.0 reveal closed lines around elliptical stagnation points with Cartesian coordinate (0.5
b, 0.2938
a) and (0.5
b, -
0.2938
a), respectively. These two stagnation points are equivalent to the zero points of velocity u at line x
=
0.5
b. It is clear that the streamfunction Γ, the velocity v, and the shear stress τ are antisymmetric about the horizontal centerline, whilst the velocity u and the normal stress σ are symmetric about this line. From Figure 2(b), we can see two zero-contour lines, which indicate that there are two zero-value points of velocity u for every section paralleling with line x
=
0.5
b. The velocity v in the domain is divided into four zones and there are also two zero-contour lines which correspond to the horizontal and vertical centerline. In the vicinity of the corners, there are stress concentration zones (Figures 3(a) and 3(b)) due to the discontinuity in velocity u. For S
=
2, two separate eddy structures (Figure 3(c)) still exist, the larger of which is that adjacent to the faster moving wall, and the stagnation points in Cartesian coordinate are (0.5
b, 0.2581
a) and (0.5
b, 0.3411
a). Figure 4(a) illustrates that, for S
>
0, there are two large eddies within the cavity with two stagnation points on the centre-line between the free surfaces. The effect on the position of the eddies for S
∈
[
0,40
] is shown in Figure 4(b). As S increases towards infinity, the stagnation point in the upper cavity moves towards the horizontal centerline and reaches to the point (0.5
b, 0.2069
a), while the stagnation point in the under cavity moves towards the bottom lid and reaches to (0.5
b, -
0.5
a) and the eddy in the under cavity reduces and disappears.Contour lines withS
=
1: (a) streamfunction Γ; (b) velocity u; (c) velocity v.
(a)
(b)
(c)Contour lines: (a) shear stressτ (S
=
0); (b) normal stress σ (S
=
0); (c) streamfunction Γ (S
=
2).
(a)
(b)
(c)(a) (u
,
S) diagram for S
>
0 at x
=
0.5
b; (b) (Y
st
,
S) diagram for S
>
0 with A
=
1.0.
(a)
(b)Consider, then, the case of flowS
<
0 with A
=
1, including the case S
=
0, which is generated by moving the top and bottom wall in opposite directions. As for flow patterns when S
=
0, the cavity consists only of a large recirculation. The streamline patterns for S
=
0 reveal closed lines around elliptical stagnation points with Cartesian coordinate (0.5
b, -
0.2069
a). We can note that the case S
=
0 is essentially the same as S
=
∞; the immobile and moving lids are merely interchanged. From Figure 5(b), we can see one zero-contour line, which indicates that there are one zero-value point of velocity u for every section parallel to line x
=
0.5
b. The velocity v in the domain is divided into two zones and there is also one zero-contour line which corresponds to the vertical centerline. In the vicinity of the bottom corners, there are stress concentration zones (Figure 6) due to the discontinuity in velocity u. For S
=
-
1, the streamlines in Figure 7(a) show a single eddy with its centre at the middle of the cavity. From Figure 7(b), we can see one zero-contour line which just corresponds to line y
=
0. The velocity v is antisymmetric about x
=
0.5
b and there is also one zero-contour line which corresponds to the vertical centerline. The shear stress τ is symmetric about both the horizontal and the vertical centerline and there are also two zero-contour lines, whilst the normal stress σ is symmetric about the center of the cavity and there are four zero-contour lines. Figure 8(c) shows the flow pattern for S
=
-
2, where it is observed that (compared with S
=
-
1) the eddy centre is displaced from (0.5
b, 0) along the line of x
=
0.5
b—towards the faster moving lid—and thus symmetry about y
=
0 is lost, and the stagnation points in Cartesian coordinate are (0.5
b, 0.1256
a). Figure 9(a) illustrates that, for S
≤
0, there is one eddy within the cavity with one stagnation point on the centre-line between the free surfaces. The effect on the position of the eddies for S
∈
[
-
40
,
0
] is shown in Figure 9(b). As S is reduced towards -
∞, the stagnation point moves towards the top lid and reaches the point (0.5
b, 0.2069
a). We note that the case S
=
-
∞ is essentially the same as the case S
=
∞, the direction of swirl is merely interchanged.Contour lines withS
=
0: (a) streamfunction Γ; (b) velocity u; (c) velocity v.
(a)
(b)
(c)Contour lines withS
=
0: (a) shear stress τ; (b) normal stress σ.
(a)
(b)Contour lines withS
=
-
1: (a) streamfunction Γ; (b) velocity u; (c) velocity v.
(a)
(b)
(c)Contour lines: (a) shear stressτ (S
=
-
1); (b) normal stress σ (S
=
-
1); (c) streamfunction Γ (S
=
-
2).
(a)
(b)
(c)(a) (u
,
S) diagram with S
≤
0 at x
=
0.5
b; (b) (Y
st
,
S) diagram with S
≤
0 for A
=
1.0.
(a)
(b)


## 4. Conclusions

This paper has presented a new solution for stokes flow in a double lid-driven cavity with free side surfaces. It is based on a symplectic approach which has been used previously in elasticity mechanics. The variational principle is derived in a geometrical symplectic space. Using the essential Hamiltonian principle with Legendre’s transformation, an eigenvalue is obtained and thus solved in the symplectic space under the Hamiltonian system formulation but not in the Euclidian space under the traditional Lagrange system formulation. Stokes flows analysis requires solving rational first order ordinary differential equations. The examples show the applicability and validity of the symplectic approach for solving stokes flow problems.



---

*Source: 151470-2014-04-24.xml*