B 14873473
B 14873473
B 14873473
Channel
Roger R. Coroas
B.Sc. University of Havana, Cuba, 1988.
MASTEROF SCIENCE
in the Department
of
Mathematics & Statistics
APPROVAL
Name:
Roger R. Coroas
Degree:
Master of Science
Title of thesis:
I
8
Dr. T. Tang
Senior Supervisor
L/
Dr. R m . Russell
Dr. R. W. Lardner
Karageorghis
External Examiner
Date Approved:
~ u l y2 5 ,
-,
my t h e s i s ,
l i b r a r y o f any o t h e r u n i v e r s i t y , o r o t h e r e d u c a t i o n a l
i t s own b e h a l f o r f o r one o f i t s users.
i n s t i t u t i o n , on
I f u r t h e r agree t h a t permission
I t i s understood t h a t copying
o r p u b l i c a t i o n o f t h i s work f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d
w i t h o u t my w r i t t e n p e r m i s s i o n .
T i t l e o f Thesis/Project/Extended
Author:
(name
Essay
Abstract
The steady, two-dimensional, laminar flow of an incompressible fluid through a channel with an abrupt constriction of ratio 2:l is numerically calculated. The governing
equations are solved using a second-order central difference scheme and a fourth-order
compact scheme. The Interior Constraint (IC) Method is employed for the numerical treatment of the values at the boundary. The Point Successive Over-Relaxation
(SOR) Method is used to solve the discrete equations. Due to the combination of
the IC Method at the boundary and the compact scheme at the interior points overrelaxation is now possible. Faster covergence is achieved for higher Reynolds number.
An evaluation of different approximations for the boundary vorticity is given. Several
flow parameters are calculated and compared to previous works.
Dedication
Acknowledgements
I would like to thank all the people who contributed in the completion of this thesis,
among whom are:
My senior supervisor Dr. T. Tang, for his guidance and help during this period.
Dr. R. D. Russell and Dr. M. R. Trummer for their friendship and encouragement
since the beginning of this program.
Dr. H. Huang for the invaluable discussions and comments.
Dr. G. A. C. Graham and Mr. P. Dobud for making possible my arrival to Simon
Fraser University.
My family and friends for their moral support.
Jocelyn Laidlaw for her love and constant understanding.
Stephen and Victoria Trotter for their help.
Antonio and Tatiana for those good meals.
Kostas Korontinis for his unquestionable friendship and for introducing me in the
art of ouzo.
Mansilla, Rafa, Leo y Pablo for their support and the good times.
Mrs. S. Holmes for her assistance and the Department of Mathematics and
Statistics of Simon Fraser University for financial support.
Contents
Approval .
....................................
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iv
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
v
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
...
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V l l l
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
1
1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2
Governing Equations of Fluid Mechanics . . . . . . . . . . . . . . . .
4
2.1
Derivation of Governing Equations . . . . . . . . . . . . . . .
4
2.1.1
Preliminaries . . . . . . . . . . . . . . . . . . . . . .
6
2.1.2
Conservation of Mass . . . . . . . . . . . . . . . . . .
7
2.1.3
Conservation of Momentum . . . . . . . . . . . . . .
8
2.1.4
Conservation of Energy . . . . . . . . . . . . . . . .
...
2.1.5
2.1.6
....................
2.2
.......
..............
2.2.1
Discretization of Equations
2.2.2
3.2
......................
Review of Previous Work . . . . . . . . . . . . . . . . . . . . .
The Interior Constraint Method in a Constricted Channel . .
3.3
3.4
............
Numerical Boundary Conditions . . . . . . . . . . . . . . . . .
3.4.1
Boundary Vortici ty . . . . . . . . . . . . . . . . . . .
3.4.2
Boundary Stream Function Extrapolation . . . . . .
3.4.3
3.5
.............
........................
Numerical Results
..........
Length of Recirculation Regions . . . . . . . . . . . .
3.5.2
4
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1
vii
List of Tables
3.1 Separation point at upper corner (L1)
..
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
................
Reattachment point at the upper corner (L2)
..............
Separation point after the reentrant corner (L3)
............
Reattachment point after the reentrant corner (L4)
..........
Stress values on boundary W l
.......................
Separation values at the upper corner (L1)for Re=O, 250 and 500. . .
Stress values on boundary W 2. . . . . . . . . . . . . . . . . . . . . . .
Reattachment values at the upper corner (L2)
for Re=O, 250 and 500.
Stress values on boundary W3. . . . . . . . . . . . . . . . . . . . . . .
3.10 Separation values after the re-entrant corner (L3)for Re = 250 and 500.
3.11 Reattachment values after the reentrant corner (L4)for Re = 250 and
500. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
viii
List of Figures
...... . ...........................
Higher order stencils . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Stencils
2.2
............................
Grid and computational boundary. . . . . . . . . . . . . . . . . . . .
Outlier values for the VIC method. . . . . . . . . . . . . . . . . . . .
Characteristic lengths for the flow in a constricted channel. . . . . . .
Iterations number for schemes. . . . . . . . . . . . . . . . . . . . . . .
3.6 Stress on wall Wl for different schemes and vorticity boundary approximations.
3.7 Stress on wall W2for different schemes and vorticity boundary approximations.
..................................
3.8 Stress on wall W3for different schemes and vorticity boundary approximations.
3.9 Vorticity distribution on wall W3for CB scheme (Re = 250 and 500).
3.10 Vorticity distribution on wall W3 for C2 scheme (Re = 250 and 500).
3.11 Vorticity distribution on wall W3 for C3 scheme (Re = 250 and 500).
3.12 Vorticity distribution on wall W3 for T3 scheme (Re = 250 and 500).
3.13 Vorticity distribution on wall W3 for T4 scheme (Re = 250 and 500).
3.14 Vorticity distribution on wall W3 for TB scheme (Re = 250 and 500) .
42
43
........
45
46
49
..........................
Chapter 1
Introduction
Fluid dynamics, prior to the computer era, was mainly divided into two branches:
theoretical and experimental. The mathematical models resulting from fluid mechanics problems are highly non-linear. Therefore, theoretical solutions are difficult to
obtain and mainly found for linear problems. On the other hand to reproduce a physical phenomenon through experiments is a difficult matter. Wind tunnels have been
designed to model flows around objects, e.g. cars, airplanes and rockets. Their main
disadvantage is cost. Wind tunnels are very expensive to run, and a change in the
parameters flow might imply changes in the layout of the experiment. Even for small
scale experiments, e.g. flows through channels, the situation is not different. Very
sophisticated instruments are needed to measure the velocity of the particles, and the
accuracy of this measurements suffer near boundaries.
The introduction of computers resulted in the growth of a completely new field,
termed Computational Fluid Dynamics (CFD). This field has led to the development
of new mathematical theories for numerical solutions of non-linear differential equations. CFD is close to the experimental branch. Improvements in the computer speed
and the search for more efficient numerical solution procedures allow us to model
more complex phenomena. But it is still necessary to rely on rigorous mathematical
analysis of simpler linearized models, on physical intuition and trial-error procedures.
CHAPTER 1 . INTRODUCTION
In Computational Fluid Dynamics the complex non-linear equations are approximated by a system of discrete equations. Many procedures may be used to reduce the
continuous equations to a set of discrete ones. Finite differences is one of the most
used methods. A grid or mesh is placed over the domain of interest. At the mesh
intersections the finite-difference solutions will be defined. Taylor series are applied
at each grid point, and a set of discrete equations are obtained. The order of accuracy of the approximation may vary. Schemes of first and second-order of accuracy
are extensively used. A disadvantage is that, for complex flows, very fine meshes are
needed and, therefore, the number of discrete equations larger. Usually higher order
of accuracy means more grid points for the scheme, except for high order compact
schemes.
Constricted tubes are found in mechanical devices and in our body (due to atherosclerosis of blood vessel). In this thesis we shall study the flow of a tw+dimensional,
steady, viscous, incompressible fluid through a constricted channel. Numerical calculations will be made for a channel with constriction ratio 2 : l. A parabolic velocity profile is imposed upstream, while horizontal flow and constant pressure are
assured for the outflow downstream of the channel. Flows in contracting channels
have attracted considerable attention. Numerical experiments suggest the presence
of a recirculation region downstream. The flow singularity at the reentrant corner
makes the prediction of such a region difficult. The values of the Reynolds number
at which the downstream recirculation appears and the wake's dimensions are still
open problems. In an effort to solve these problems several flow parameters will be
analyzed. We will use a second order scheme and a fourth order compact scheme
and compare their performance. A careful examination of the effects of the boundary
conditions will be also presented.
In Chapter 2 we will present an overview of the continuous governing equations
and their discretization. A derivation of the governing equations of fluid flow, with
all the hypotheses involved, will be given in Section 2.1. In Section 2.2 we will give
a description of the numerical schemes employed to compute the flow. The solution
technique used will also be described.
CHAPTER 1. INTRODUCTION
Chapter 3 begins with a description of the model problem, and a review of the previous work is given in Section 3.1. The Interior Constraint (IC) method is formulated
in Section 3.2 for a constricted channel. A variation of the IC method, named the
Vorticity Interior Constraint (VIC) method, is proposed in Section 3.3. The boundary
conditions employed on this problem are discussed in Section 3.4 and a series of vorticity boundary approximations are also given. Section 3.4 is dedicated to the results.
The effects of different vorticity boundary conditions are studied. Characteristic flow
parameters are computed, and comparisons with published ones are also.given.
Chapter 2
Governing Equations of Fluid
Mechanics
2.1
2.1.1
Preliminaries
A fluid consists of many particles in interaction. The behavior of the particles ultimately characterizes the motion of the fluid. One could concentrate on the description
of the interaction of such particles and, using the laws of mechanics and probability
theory, infer the macroscopic behavior of the fluid. Such an approach is known as
the statistical method, and is valid for light gases yet incomplete for dense gases and
liquids [5].
Recently a new microscopic approach has been developed. The lattice Boltzmann
equation (LBE) method [16, 171, a derivation of the lattice gas automaton method
[30], constructs a microworld in which space, time and velocities are all discrete. The
particles reside on the lattice's nodes and change their velocities directions according
to scattering rules. This method allows the use of fully parallel algorithms and is
easily implemented for complex boundary flows.
A different approach is to consider a small volume of continuous fluid, where
individual particles are ignored. The microscopic properties of the particles within
the volume are averaged to obtain the macroscopic values such as: velocity, pressure,
density, or other variables that might describe the fluid. The validity of this continuow
method is conditional upon the satisfaction of the continuum hypothesis. The physical
quantities of interest must be spread uniformly over the volume in consideration, and
the mean free path of the particles must be small in comparison with the length of
the flow field. Only under this continuum hypothesis can the microscopic structure
of the fluid be ignored, and the average considered valid.
In this work, the fluid under consideration satisfies the continuum hypothesis.
Therefore, the continuous method will be used to state the conservation laws.
The conservation laws could be formulated in either one of the two known coordinate systems: the Eulerian coordinate system or the Lagrangian coordinate system.
In the Eulerian system a control volume independent of time is fixed in space. The
conservation laws are applied to the fluid which passes through the control volume.
The dependent variables are functions of the spatial coordinates and time. On the
other hand, the Lagrangian system follows a control volume through its motion in the
flow. The particles of fluid do not vary, but the shape of the volume might change
over time. In this formulation, the spatial variables are not independent of time. The
choice of coordinate system varies from author to author, and mainly depends on the
problem under consideration. Here the Eulerian coordinate system is used.
Recalling that the independent variables are the spatial coordinates and time,
it is very useful to define the total derivative, also known as convective or material
derivative. The total derivative gives the rate of change of any physical magnitude 3
when moving with the flow, and it is defined as
2.1.2
Conservation of Mass
Consider a fixed control volume V entirely occupied by fluid, and S its surface. The
volume V is fixed in space. The particles will enter and leave the volume through the
control surface S. These conditions will also apply in the subsequent sections of this
work.
The law of conservation of mass states that the mass within any volume varies
only as a consequence of the flux across the surface of the volume. The integral form
of the conservation of mass applied to the fluid confined in V yields the equation
where p is the fluid density, and n the outward unit normal of the surface S.
The previous equation holds for any control volume V. A more restrictive form of
the conservation of mass can be obtained applying the divergence theorem
+ v .(pu) = 0
which is known as the continuity equation in the conservative form. Equation (2.3) is
not valid when the quaetities involved are not continuous.
A non-conservative form of equation (2.3) can be obtained by expanding the divergence term V - ( p u ) in (2.3) and using the expression for the total derivative
In most liquid flows the total variation of the density p may be considered irrelevant
and, therefore, assumed to be uniform over the entire fluid. In such cases
holds and the fluid is said to be incompressible. The continuity equation (2.4) then
reduces to
v.u =
2.1.3
0.
(2.6)
Conservation of Momentum
The principle of conservation of momentum states that the rate of change of momentum within a control volume V is equal to the sum of forces acting in such a volume,
and it is expressed by
where u is the velocity, F the external force acting on the fluid contained in V, and
a ='( a i j )represents the surface stress tensor.
Equation (2.7) holds for any choice of volume V. If the integrands are continuous
functions of x on V the differential equation
where p is the kinematic pressure, 4,is the Kronecker delta, p the dynamic coefficient
i(e+z)
+ pV2u.
(2.10)
Substitution of (2.10) into (2.8) gives the well known expression of the Navier-Stokes
equations
In the majority of the cases the force F represents the gravitational force. The
term pF in the above equation can be absorbed by the pressure term Vp when the
pressure is written as
(2.12)
p = PO pF-x p
where po is a constant, pF.x the pressure to balance the force per unit of volume due
to F and p is the pressure arising from the effect of the motion of the fluid- it is known
as the modified pressure [I].
Substitution of (2.12) into (2.11) gives the formulation of the Navier-Stokes equations
in ~rimitivevariables
Some authors refer to equations (2.13)-(2.6) as the Navier-Stokes equations. The same
will be done throughout this work, unless otherwise indicated.
2.1.4
Conservation of Energy
the total work done on the system and any heat added. We are interested in fluids
in which the temperature is uniform, and therefore the energy equation will not be
necessary here.
2.1.5
We are interested in exploring the effects that changing the values of the parameters
p and C( have on the flow. In this case it is useful to write equations (2.13) and (2.6)
in terms of dimensionless variables. Taking L as the reference characteristic length
and U as the reference characteristic velocity, the new non-dimensional variables are
defined
Without any loss of generality the same notation for the new variables is kept.
With these new variables the governing equations (2.13)-(2.6) become
where Re =
is the Reynolds number.
The system of equations (2.15)-(2.16) is known as the dimensionless Navier-Stokes
equations in primitive form.
In many problems it is important to study the long time behavior of the flow.
The steady state equations are obtained by assuming that the physical magnitudes
that describe the flow are independent of time. That is, in equations (2.15)-(2.16) the
velocity u and pressure p are only functions of the position, u = u(x), p = p(x).
In this case the total derivative can be written as
Du- -
Dt -
10
2.1.6
In many cases it is enough to consider planar flows. Due to symmetries in the region
of the. problem, three-dimensional flows could be studied as two-dimensional flows.
This means that the velocity u is only described by two components u = (u, v, 0).
It is common practice in fluid dynamics to visualize the flow field by drawing the
streamlines. Therefore, it is convenient to introduce the scalar magnitude +, named
the stream function. The stream function is constant along the streamlines and is
related to the velocity through the equations
which satisfy the continuity equation (2.19). Vortex creation is also a distinctive
feature of some flows. It is useful to describe the flow in terms of the distributions of
the vorticity R, defined as
R = - - -
av
ax
au
ay
which is created on the boundary and diffused away by the viscosity. The NavierStokes equations in the vorticity-stream function formulation for steady, incompressible, planar flows take the form
where V2 is the Laplacian operator. Note that the pressure term disappears from the
governing equations. The elimination of the vorticity in equations (2.22)-(2.23) leads
to the formulation of the Navier-Stokes equations in terms of the stream function
where
V4 = & +2,-#
a4
+w
2.2
11
The equations of fluid mechanics are highly nonlinear. Analytical solutions are rarely
obtained. The appearance of computers opened a new era in fluid dynamics research.
A new field, termed Computational Fluid Dynamics (CFD), has emerged. This field
has provided new mathematical methods to solve numerous problems of fluid mechanics. Solutions are obtained by first replacing the derivatives appearing in the
governing equations by appropriate discrete approximations and then solving the discrete equations by a suitable numerical method.
2.2.1
Discretization of Equations
Numerous methods have been developed to discretize the governing equations. Among
them, the finite difference method is the most commonly used. Finite difference
approximations of different orders are found with the use of Taylor expansions of the
quantities involved.
XN =
b and h =
by evaluation at the discrete points xi, i.e. fi = f (xi). Different finite difference
approximations for the derivatives of f (x) at x = xi can be obtained from these. For
instance the first derivative could be approximated by any of the following expressions
centered second-order
backward first-order
12
10
6
13
(b)compact
fi+l - fi
h
+ O(h)
forward second-order
Approximations can also be found for higher order derivatives. The classical difference
approximation for the second-order derivative at x = xi is
which has second-order accuracy. Combinations of these finite difference approximations may be put together to obtain the finite difference approximation of the
governing equations.
To discretize equations (2.22)-(2.23) consider a uniform mesh in both the x and y
directions. Let h be the mesh step size.
The well known central difference or 5-point scheme for the Poisson equation can
be obtained using the central difference formula (2.29) in equation (2.22) to obtain
where the subscripts 0,1,2,3 and 4 represent the grid points (x, y), (x h, Y),(x, y
h), (x - h, y ) and (x,y - h) respectively (see Figure 2.la). Similarly the vorticity
(a) Wide
13
(b) Compact
where
Re (s1 s 2 ) -
14
+ s3
where
S1 =
s2 =
s3 =
are now needed, and they are numbered 9,10,11 and 12, respectively (see, e.g., [21]).
On this stencil the scheme for equation (2.24) is
where
and
2.2.2
15
Once the governing equations are discretized, it is necessary to employ a suitable numerical method to solve the system of discrete equations. Several methods are available. Direct methods are, generally, memory consuming and difficult to implement
for complex geometries. Iterative methods are more frequently used in computational
fluid dynamics.
The Point SOR method is an iterative method that can be used to solve discretized elliptic problems. This method is very easy to implement, requires very little
memory and is one of the most efficient point-iterative procedures for large systems
of equations. A sufficient condition for the convergence of this method is diagonal
dominance.
wept
has been obtained (see, e.g., [27]). For complex elliptic problems it is not
possible to obtain
wept
JS OF FLUID
where A is the relaxation parameter and 0; is calculated from any of the given approximations.
The iterative process (2.36)-(2.37) is stopped when, for a given value e, the condit ions
n;+l
where 0 < 6
puted.
sn; + (1 - qn;,
(2.40)
Chapter 3
Flow in a Constricted Channel
3.1
Flows in constricted tubes are frequently found in many technological devices and
even inside the human body, as in the case of blood vessels with stenoses. Such flows
exhibit complex patterns. Recirculating flow is usually found before and after the
constriction. This behavior is, in most cases, not desirable. A good understanding
of the flow structure is vital to prevent unwanted behavior and properly control the
process.
Flows in channels with a sudden contraction have been studied by several authors
[ l l , 10, 8, 20, 22, 15, 181 and present numerous interesting problems. A downstream
wake is observed for a certain value of the Reynolds number Re and higher values.
The fine structure of the flow is not well known yet. The singularity of the flow
at the reentrant corner makes it difficult to predict the length of the downstream
recirculation region and the value of Re at which first separation appears.
The first effort to numerically study separating flows in channels was produced
by Greenspan [l11 and Friedman [lo]. They investigated the flow in a channel with
a square constriction on one wall. Using a finite difference discretization for the governing equations in the stream function-vorticity formulation, Friedman [lo] obtained
<
<
ft,
< <
points are densely distributed in a vicinity of the corner and near the walls. They
reported that the downstream wake first appears at around Re = 175 and that the
separation occurs at the corner and not to the right as predicted by Hunt [20]. They
found a second vortex at the upper corner upstream of the step, which is in agreement
with Moffatt's theory [25].
Hawken et al. [15] studied this contraction flow using a finite element discretization
for the governing equations in primitive variables. A high concentration of elements
is placed close to the re-entrant corner. Even though they failed to visualize the
downstream vortex, the velocity gradient revealed a very thin recirculation region at
Re = 250 near the corner. They suggested that separation occurs to the right of the
corner. Their downstream wake reattachment values are slightly higher than Hunt's
values, but smaller than Karageorghis and Phillips'.
In a more recent paper Huang and Seymour [19] employed the Interior Constraint
(IC) method, proposed by Huang [IS], to solve this problem. A computational boundary is set up near the physical boundary and the use of vorticity values at the physical
boundary is avoided. Both uniform and non-uniform grids were employed in calcula-
A.
Re increases.
The treatment of the singularity of the vorticity function at the sharp corner has
also been a concern. Roache [27] reviews different methods of handling the corner
vorticity. A local coordinate transformation near the re-entrant corner was proposed
by Dennis and Smith [S]. The finite difference approximation was applied along lines
with 45" of inclination with respect to the main grid lines, to avoid the use of the
vorticity at the corner. The use of the vorticity at the boundary was avoided by Hunt
[20], who used the equations in terms of the stream function only. But as noted by
him the singularity of the vorticity leads to a singularity of the second derivative of
the stream function, that might affect the accuracy of the results. Hunt also proposed
the use of Moffat t 's expansion [25], successfully employed by Bramley and Dennis [2]
in a branching channel.
The vorticity does not appear in Hawken e t al. [15] equations, since the primitive
3.2
We are interested in the simulation of the flow of an incompressible fluid in a constricted channel, with constriction ratio 2 : 1. Considering that the flow is symmetric
with respect to the mid line of the channel, we can restrict our calculations to the
upper half of the channel (see Figure 3.1). The constricted wall of the channel is Wl=
{ y = 1, for x
21
and y = 0 is a symmetry wall. The upstream and downstream boundaries are set
at x = xin and x = xmt respectively. A parabolic Poiseuille flow is imposed on the
upstream boundary.
The Navier-Stokes equations formulated in terms of the stream function and vorticity are the governing equations
where J ( $ , Q) = +,R,
- $,R,
+ = o , R=O,
on y = 0 ,
+=1,%=0,
1
on y = l f o r x < O and y = - f o r x > O ,
(3.3)
(3.4)
on z=x,t
1
and O < y < 2'
av a9
0 and
outflow boundary.
A uniform mesh is set up with N and M grid points in the x and y directions,
respectively. The mesh size is taken to be
(3.10)
vi+ij= -Rij,
(3.11)
for i = 2 ,...,H N - 1 ,
and i = H N ,..., N - 1 ,
and i = 2 ,...,H N - 2 ,
j = 2 ,...,H M - 1 ,
j = 2 ,...,H M - 2 ,
j = H M,...,M M - 2 ,
for i = 2 ,..., H N - 1 , j = M - 1 ,
and i = H N l , j = H M ,..., M - 2 ,
and i = H N ,...,HN-1, j = H M - 1 .
The one sided discretizations of the non-slip conditions on the physical boundary are
used to interpolate the values of the stream function at the computational boundary
Di$ij = 0
(3.13)
...,H N - 1, j = M
i = HN, ...,H N - 1, j = HM
for i = 2 ,
and
where Dz$ij and Di$ij are the discretization at the grid point (xi, yj) of $, and $,
respectively. If, for instance, a second-order approximation for 111, = 0 is employed
at the physical boundary, the value of the stream function at the computational
boundary +i-lj can be isolated from the expression QtDij= 3tDij - 4+i-lj $i-2j = 0.
3.3
The purpose of the IC method is to avoid using the boundary vorticity values in the
calculations. The values of the stream function at the computational boundary are obtained by interpolation using the non-slip conditions at the wall. Here we will propose
a variation of the IC method, which will be named the Vorticity Interior Constraint
(VIC) method. The computational boundary is kept when vorticity values are calculated. For the stream function, the grid points adjacent to the contracting wall are
0 Grid Points
Outliers
Figure 3.3: Outlier values for the
VIC method.
now treated as interior points. The grid point 8 is now added to the computational
boundary.
The biharmonic formulation [27] of the Navier-Stokes equations (2.24) is now imposed at these grid points and interpolation from the non-slip condition is no longer
necessary. The non-slip conditions at the wall are now used to eliminate the outlier
values (see Figure 3.3).
For the VIC method the discretized vorticity values at the computational boundary
are calculated as in (3.12), and the stream function values are calculated as follows:
vt$ij= Re
Q$ij
(3.15)
JBh($ij),
(3.16)
= 0,
for i = H N , j = H M , ...,M
-1
(3.17)
and
V;$ij = Re JBh($ij),
(3.18)
D:$ij = 0,
(3.19)
for i = 2, ..., H N
- 1, j
=M
and i = H N , ..., H N - 1 , j = H M
where V t and JBh are defined as before. Note that this will give second-order accuracy
for the stream function for points adjacent to the contracting wall.
3.4
3.4.1
Boundary Vorticity
As mentioned earlier, if the IC or the VIC methods are used, the vorticity values
on the wall are not needed to calculate the flow inside the domain. Nevertheless,
some characteristic parameters, such as reattachment points, are usually provided to
establish the validity of the solution. These parameters can be obtained from the
values of the vorticity at the boundary. Roache [27] mentioned the importance of
choosing adequate formulae for the vorticity boundary. A wide variety of boundary
approximation formulae were tested by Gupta and Manohar [12] for the square cavity
problem.
It is our concern that the comparison parameters provided might depend on the
boundary vorticity formula used, even though such formulae are not used throughout
the calculations. A survey of different formulae for the vorticity at the wall is pr*
duced. Careful comparisons are made to determine the nonreliable approximations.
Let us assume that the grid point (xi, y j ) lies on the boundary Wl= { y = 1 and
xi,
< HN.
One of the first boundary vorticity formulae was produced by Thom [28]
This first-order formula has been used extensively in flow calculations and often gives
results as accurate as higher order formulae [27].
The second-order formula, developed by Woods [31], has been widely used
It also involves a vorticity value adjacent to the boundary. Some unstabilities and
low convergence rate have been reported for this formula.
Here we will also consider a whole class of first-order and second-order formulae,
called the (p, 0 ) and (p, q ) formulae, respectively, proposed in reference [12].
- 9)".
In this work we will consider from the class (3.22) the first-order formulae
Gupta ant Manohar [12] pointed out that for larger values of p and q convergence
is achieved faster but the solution becomes less accurate. In our work the boundary
approximations of the vorticity will not affect the convergence process, since they
are not used throughout the calculations. However they are expected to affect the
accuracy of the flow parameters calculated afterward.
27
In previous works the sizes of the recirculation regions are taking as comparison
parameters. These values depend strongly on the boundary vorticity approximation
used, since they are usually calculated using the boundary coordinate where the vorticity is equal to zero. Isolated values of the vorticity are not reliable parameters
to measure accuracy [12]; therefore, one should provide some other parameters for
comparison purposes.
As an alternative parameter the total shear stress on the contracting wall is calculated. Here a non-dimensional value of the shear stress, introduced in [12], will be
used. On the contracting walls Wl, W2 and W3 the shear force is defined by
As shown in [12], the shear stress at the wall reflects the accuracy of the vorticity
boundary approximations employed.
3.4.2
As ,mentioned earlier when the IC method is used the value of the stream function at
the computational boundary are extrapolated from the non-slip boundary conditions.
The most interesting features of the flow in a channel with a contraction appear at
the contracting wall. In reference [19],one sided second-order approximation formulae
are used to approximate the non-slip conditions and second-order central differences
are used at the interior points. Higher order approximations for the non-slip conditions
are proposed and used here.
Assuming that our boundary point lies on the upstream contracting wall Wl,
second, third and fourth-order one sided formulae [6] are defined by
28
3.4.3
One might think that the outflow boundary condition will be the easiest to impose,
setting the end of the channel far away from the constriction. It has been found
numerically that instabilities might propagate upstream and invalidate the solution
if the outflow conditions are too restrictive [27]. On the other hand, the use of large
values for xOutincreases the cost of the calculations and, to reduce the cost, some
kind of domain transformation will be needed. Therefore, it is very important to use
economical outflow boundary conditions.
29
by Greenspan [ll] and Friedman [lo], are employed. Horizontal flow is guaranteed
downstream by equation (3.8) and equation (3.9) gives constant pressure along the
outflow boundary at steady state.
3.5
Numerical Results
The calculations were carried out for Reynolds number up to 500 on a uniform mesh
with size h =
= -2
LI,of
the upper corner vortex and its width L2 , the separation and reattachement values
of the downstream wake, L3 and L4, are all measured from the corners (see Figure
3.4). The stress on each wall section of the contracting boundary is also monitored.
Numerical experiments show that the parameters L3 and L4 are sensitive to the grid
size, and L1 and L2 behave in a more stable way.
For comparison purposes the outflow boundary condition was also set at
xout
= 3.
Calculations were made using the fourth-order compact method with extrapolation of
fourt h-order. For the Reynolds number employed significative differences in the flow
parameters were not found, a s suggested in [8] and [19].
The solution is calculated for Re with increments of 50. For Re = 0 the initial
approximation is taken to be zero at every interior point. For Re
30
calculated for previous Reynolds number is used as initial approximation. This means,
for example, that the converged numerical solution for Re = 400 is used as the initial
approximation to calculate the solution for Re = 450. Continuation on Re is necessary
because for large values of Re the iterative process is not robust enough to converge
from an initial approximation of zero everywhere.
The iterative process considered requires the use of relaxation parameters. For
both schemes, central differences and compact, the boundary vorticity damping parameter b was set at 0.6.
At the interior points, when using the central differences scheme, overrelaxation
is possible for the stream function equation and w was set to be 1.4 throughout the
calculations. Under-relaxation is necessary for the vorticity equation. Initially the
value of X is set at 0.8, but the iteration process becomes divergent for Re = 450. For
450
p = 0.6.
The iterative process is stopped when conditions (2.38)-(2.39) are satisfied at the
interior points. The value of c is set as 0.0005 and
Figure 3.5 shows the number of iterations required for the iterative process to
converge for the schemes at different Reynolds numbers. Recall that continuation on
the Reynolds number is used to calculate the solution. The solution for Re = 0 is
calculated from an initial interior approximation of zero everywhere. This is the reason
why the iteration number for Re = 0 is considerably higher than for the rest. Note
that the compact scheme requires less iterations overall, but the cost per iteration is
higher than for the central difference scheme. A good feature of the compact scheme
3.5.1
Gupta and Manohar [12] found that the numerical boundary vorticity might affect the
flow structure. It is our interest to find the most appropriate vorticity formula for every
interior scheme used in the calculations. A theoretical solution is not known for the
flow in a constricted channel. Reference papers provide, as comparison parameters,
the lengths of the upstream and downstream wakes. These parameters are calculated
from isolated values of the vorticity or the velocity. As mentioned before, isolated
values should not be used to test accuracy. They should be used in combination with
more reliable parameters. Gupta and Manohar [12] mentioned that, even though
isolated values of the vorticity on the boundary are not reliable, the integral over the
boundary is a reliable parameter. This integral gives a non-dimensional stress on the
wall. The stress on each of the three sections of the contracting wall is calculated. The
vorticity at the re-entrant corner is considered bivalued and therefore discontinuous.
Most figures and tables shown in this chapter make use of abbreviations for schemes
and boundary vorticity formulae. Refer to pages 32-33 for full description of abbreviations.
33
< 500.
sponds to a different numerical scheme. The stress calculated with different boundary
vorticity approximations is plotted.
Figure 3.6: Stress on wall Wl for different schemes and vorticity boundary approximations.
It is obvious that there are no significant differencesfor Re
all the boundary vorticity approximations should produce similar results for Reynolds
numbers up to 100.
34
Significant differences are not found for L1 and L2 among the different bound-
<
Figure 3.7: Stress on wall W2 for different schemes and vorticity boundary approximations.
The values for the boundary stress F3 on W3 are plotted in Figure 3.8. The
value of F3 decreases from the first-order to the second-order formulae. A fairly
good agreement among the boundary vorticity conditions for Re
5 50 is found. Now
Woods' boundary condition gives the smallest stress values overall. In reference [12],
Gupta and Manohar described that Woods' boundary condition had the tendency to
over-estimate the stress values on the moving boundary for the square cavity problem.
Figure 3.8: Stress on wall W3 for different schemes and vorticity boundary approximations.
We presume that the same behavior is present in this case. Woods' formula might
be under-estimating the value of the stress on W3. Since the upstream results are
more reliable, and the same pattern is not observed on Wl and Wz,we cannot assure
that Woods' condition is "better" in this case. Values of F3 for Re = 0,250 and 500
are given in Table 3.9. Since the downstream recirculation region was not found for
Re 5 250, the values of L3 and L4 for Re = 250 and 500 are given in Table 3.10 and
Table 3.11. Note that the plots for CB differ totally from the rest. This indicates that
37
scheme CB is not trustworthy on this boundary. A quick inspection of Table 3.10 and
Table 3.11 shows that the length values obtained for this scheme disagree with the
values obtained by other schemes.
In the case of the flow in a channel with a sudden contraction, as reported by
Gupta and Manohar [12] for the cavity flow, there are differences in the values of the
wall stress for different boundary vorticity conditions for a fixed scheme. The values of
Fl and F3 decrease as the accuracy of the boundary condition increases. The opposite
behavior is found for FZ.
Figure 3.9 gives the vorticity distribution on W3 for Re = 250 and Re = 500.
Figure 3.9: Vorticity distribution on wall W3 for CB scheme (Re = 250 and 500).
38
'aR=
wood..T)wm md (2.0)
-- - - (43)
csn
.-.-.-.
..........
(3.1)
( 2 1 L wood..nam ad (2.9
Figure 3.10: Vorticity distribution on wall W3 for C2 scheme (Re = 250 and 500).
Examining the plots we conclude that scheme CB does not reproduce the expected
behavior, which is in agreement with the results previously found for the stress on
that boundary.
The (4,3) boundary formula suits the best scheme C2. However, this scheme is
unable to capture the recirculation region for Re = 250 (see Figure 3.10).
Scheme C3 employs the same interior approximation as scheme C2. They only
differ in the order of the extrapolation formula for the stream function on the computational boundary. The use of a third-order formula changes dramatically the solution
(see Figure 3.11).
39
Figure 3.11: Vorticity distribution on wall W3for C3 scheme (Re = 250 and 500).
Now the recirculation region is visible for Re = 250. Also note that now for
Thom's formula the boundary vorticity comes close to zero at Re = 250. This means
that Thom's formula is not good enough to produce recirculation downstream, but
clearly is the best of all first-order formulae. The analysis of the boundary stress
suggests that the (3,2) boundary formula suits this scheme the best.
Schemes T3 and T4 produce very close results (see Figure 3.12 and Figure 3.13).
It seems that the order of the extrapolation formula used for the stream function on
the computational boundary does not have a significant effect on the compact scheme.
Figure 3.12: Vorticity distribution on wall W3for T3 scheme (Re = 250 and 500).
The only difference occurs for the (4,3) formula for both Reynolds numbers considered. Recall that the (4,3) formula is the less accurate of all second-order formulae.
For both schemes the (3,2) vorticity boundary formula seems to work better. Furthermore, all the second-order formulae seem to work fine with this scheme. Based
upon the results for the stress on the boundaries the (3,2) formula has been chosen
as the most reliable.
There is a significant difference in the minimum value of the boundary vorticity
for both values of Reynolds numbers Re = 250 and Re = 500. The more accurate the
scheme, the smaller the minimum value become.
Figure 3.13: Vorticity distribution on wall W3 for T4 scheme ( R e = 250 and 500).
Note that none of the first-order formulae reproduced the wake after the contraction for Re = 250. Also Thom's formula seems to be the best of these.
The TB scheme produces a very small wake for Thom's formula at Re = 250 (see
Figure 3.14). This scheme gives the smallest vorticity values of all at the boundary
Figure 3.14: Vorticity distribution on wall W3for T B scheme (Re = 250 and 500).
So far we have determined which vorticity boundary approximation is best to use
with each scheme. All the boundary formulae chosen are of second-order of accuracy.
An interesting feature occurs near the upper corner for every scheme. An oscillation of the vorticity values is observed for most of the second-order formulae. Figure
3.15 shows the vorticity values for scheme T4 for two second-order boundary formulae
at Re = 500. Formula (4,3) behaves in the same way formula (3,2) does on this plot.
The rest of the second-order formulae, and the first-order Thom7s formula, give a
similar behavior to the shown in Figure 3.15 for formula (2,l).
Figure 3.16 shows the boundary vorticity values for Re = 0,150,250,350 and 500
on Wl and on W3for schemes C3, T3 and TB. Formula (2,O) is used for the plots on
Wl and the (3,2) formula for the plots on W3. TB uses the (2,l) formula on W3.
Figure 3.6 and Figure 3.8 might lead us to think that Fl has a maximum for
0 < Re < 50 and that F3 has a maximum for 100 < Re < 200. Further calculations
were made in that range of Reynolds number using the T4 scheme. It was found that
Fl attains the maximum value of 4.487 at Re = 15 and that F3 attains a maximum
of 29.189 at Re = 140. The values of the Reynolds number at which the wall stress
attains a maximum on W3 might be related to the value of the Reynolds number at
which the recirculation first appears.
Figure 3.16: Vorticity distribution on walls Wl and W3 for schemes C3, T4 and TB.
3.5.2
Figure 3.17 shows the streamline contours for different values of Re. The calculations
were made using the compact scheme and fourth-order extrapolation on the boundary.
A magnification of the streamlines and vorticity contours near the re-entrant corner
is also given (see Figure 3.18).
Let us first consider the solution upstream of the step. The vortex in the upper
corner is present for all values of Re. Fine meshes are not necessary to obtain reliable
results. In reference [lo] Friedman obtained this vortex using h = &. He stated that
the vortex reduces its size when 0
monotonically for Re
> Remin.
Figure 3.18: Magnification of streamlines and vorticity contours for Re = 500 near
the re-entrant corner for scheme T4-(3,2).
Dennis and Smith [8] used mesh size h = & and found that L1 attains a minimum
at about Re = 50, and it increases monotonically for 50
produced the vortex length and width for different values of Re. In reference [22]
Karageorghis and Phillips set the minimum value of L1 at Re = 45. Our calculations
using the scheme T4 and (3,2) formula show that Re,;, = 45 and the wake length L1
is 0.1286.
In Table 3.1 and Table 3.2 the values for L1 and L2 as well as those previously
obtained by other authors are shown. In general L1 and L2 show good agreement.
Re
Dennis & Smith [8]
Hunt (upwind) [20]
Hunt (central) [20]
Karageorghis & Phillips [8]
Hawken et al. [15]
Huang & Seymour [19]
Burggraf [19]
Dennis & Smith [19]
C3-(3.2)
47
Re
Dennis & Smith [8]
Hunt (upwind) (201
Hunt (central) [20]
Karageorghis & Phillips [8]
Hawken et al. 1151
Huang & Seymour [19]
Burmraf 1191
'Dennis & Smith 1191
(L2).
Karageorghis and Phillips [22], observed a second vortex near the upper corner
for Reynolds number as low as Re = 10. This second vortex was not observed in our
calculations.
The most interesting and controversial feature for the flow channel appears at
the re-entrant corner. A downstream recirculation region is formed after the tip of
the corner. Dennis and Smith in their early work [8], obtained a downstream eddy
formation for Re = 1000 and Re = 2000 with h = $. They predicted a recirculating
flow for Re = 500 with separation at a point on the interval 0
reattachement point to the right of x = 1.2. They were certain that no downstream
separation occurs for 0
5 Re 5 100 and
100 < Re < 500. Hunt [20] found a downstream vortex for Re= 250, with separation
to the right of the corner. Karageorghis and Phillips [22] placed the apparition of
the recirculating flow at about Re = 175, with separation at the corner. Their
reattachement values 0.500 and 0.995 are larger than the values obtained by other
authors.
In reference [15] Hawken employed a finite element method for primitive variables
formulation. He predicted a vortex, with separation to the right of the corner, starting
at Re = 250. His values are slightly larger than predicted by Hunt.
- J
1
Dennis & Smith 1191
I
I
0.119 1 - 1 - 1 0.414 1
0.500 ' 0.700 ' 0.924 0.995 '
0.155 0.301 0.423 0.477
0.126 1 - 1 - 1 0.444
I
0.-
0.0
0.-
49
3 om0.01s.
0.01
Tables
Table 3.6: Separation values at the upper corner ( L 1 )for Re=O, 250 and 500.
Table 3.8: Reattachment values at the upper corner ( L z ) for Re=O, 250 and 500.
Table 3.10: Separation values after the re-entrant corner ( L 3 )for Re = 250 and 500.
(L4)for Re
= 250 and
Chapter 4
Conclusions
In this thesis we have studied the flow through a constricted channel for moderate-high
Reynolds number. The Interior Constraint (IC) Method and the Vorticity Interior
Constraint (VIC) Method were employed to solve the steady Navier-Stokes equations.
Both methods were applied using a second-order central difference scheme and a
fourt h-order compact scheme. The results show that, for the second-order met hod,
an increase of the order for the stream function extrapolation formula, from second
to third-order, produces a qualitative change in the solution. The solution using the
compact scheme is not visibly affected by an increase of order for the extrapolation
formula from third to fourth-order. The compact scheme allows the use of overrelaxation and the number of SOR iterations is reduced considerably. It also allows
us to compute solutions for higher values of the Reynolds number. Iteration-wise the
VIC method is similar to the IC method for both schemes.
Several vorticity boundary approximations are tested for each method. The analysis of the boundary stress allows us to choose the most adequate boundary formula.
Several flow parameters, e.g. separation and reattachement points, are also calculated, and comparisons with published results are given. Our results show that the
combination of the VIC method and the second order central differences scheme is
not suitable to solve this problem. The compact scheme seems to work well for both
IC and VIC methods. The results obtained, for each case, agree qualitatively with
CHAPTER 4. CONCLUSIONS
To use the most reliable scheme and boundary vorticity to investigate the solutions for higher values of Reynolds numbers.
Bibliography
[I] G. K. Batchelor. An Introduction to Fluid Dynamics. Cambridge Univ. Press,
1967.
[2] J. S. Bramley and S. C. R. Dennis. The numerical solution of two-dimensional
flow in a branching channel. Computers and Fluids, l2(4) :339-355, 1984.
[3] J. S. Bramley and D. Sloan. A comparison of an upwind difference scheme with
a central difference scheme for moderate reynolds number. Technical Report 3,
Department of Mathematics, The University of Strathclyde, Glasgow, Scotland,
1988.
[4] L. Collatz. The Numerical Treatment of Differential Equations. Springer-Verlag,
Berlin/New York, 3rd edition, 1966.
[5] I. G. Currie. Fundamental Mechanics of Fluids. McGraw-Hill Inc., 1974.
161 B. P. Demidovich and I. A. Maron. Computational Mathematics. MIR Publishers,
Moscow, Russia, 1973.
[7] S. C. R. Dennis and J. D. Hudson. A difference method for solving the Navier
Stokes equations. In Proc. 1st Conf. on Numerical Methods in Laminar and
Turbulent Flow. Pentech Press, London, 1978.
[8] S. C. R. Dennis and F. T. Smith. Steady flow through a channel with a symmetrical constriction in the form of a step. Proceedings Royal Society London A,
372:393-414, 1980.
BIBLIOGRAPHY
61
[9] F. Durst and T. Loy. Investigations of laminar flow in a pipe with sudden contraction of cross sectional area. Computers and Fluids, l3(l):15-36, 1985.
[lo] M. Friedman. Laminar flow in a channel with step. Journal of Engineering
BIBLIOGRAPHY
62
[20] R. Hunt. The numerical solution of the laminar flow in a constricted channel at
moderately high Reynolds number using Newton iteration. International Journal
for Numerical Methods in Fluids, 11:247-259, 1990.
[21] L. V. Kantorovich and V. I. Krylov. Approximate Methods for Higher Analysis.
BIBLIOGRAPHY
63