Citation
Finite element modeling for energy dissipation in a pressure flow

Material Information

Title:
Finite element modeling for energy dissipation in a pressure flow
Creator:
Guo, Chwen-geng
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
126 leaves : illustrations ; 28 cm

Subjects

Subjects / Keywords:
Streamflow -- Experiments ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 102-104).
Thesis:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Civil Engineering
Statement of Responsibility:
by Chwen-geng Guo.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
19783233 ( OCLC )
ocm19783233
Classification:
LD1190.E53 1988m .G86 ( lcc )

Downloads

This item has the following downloads:


Full Text
FINITE ELEMENT MODELING FOR ENERGY
DISSIPATION IN A PRESSURE FLOW
by
Chwen-geng Guo
.S., National Cheng-Kung University, Taiwan, 1978
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Civil Engineering
1988


This thesis for the Master of Science Degree by
Chwen-geng Guo
has been approved for the
Department of
civil Engineering
by
Date
(Zb
/fW
Ronald C. Scherer


Guo, Chwen-geng (M.S. Civil Engineering)
Finite Element Modeling for Energy Dissipation in a Pressure
Flow
Thesis directed by Assistant Professor James C. Y. Guo.
A finite element numerical model using quadratic
isoparametric elements was developed to analyze a two-
dimensional pressure flow and to predict the flow energy
losses associated with the friction and cross sectional
contraction. The unsteady flow governing equations were
converted into stream function and vorticity diffusion
equations. Initial conditions were determined by potential
flow and boundary conditions included a uniform incoming
flow and non-slip condition on the solid wall.
Flow cases studied included the opening ratios
ranging from 0.3 to 1 and Reynolds numbers varying from 50
to 200. In general, this numerical approach provides stable
computations for flow to reach its steady state when >-
Reynolds number is less than 200 and opening ratio is larger
than 0.3. Although for high Reynolds number flows the exit
velocity profile may slightly deviate from the exact
velocity profile, the predicted friction loss and orifice
loss still give good agreements with the theoretical
solutions and laboratory data.


ACKNOWLEDGMENTS
This project was supported by the Denver Center for
the Performing Arts and funded by the Helen G. Bmfils
Foundation.
I would like to thank Dr. Ronald C. Scherer, senior
scientist, Denver Center for the Performing Arts, and Dr.
James C. Y. Guo, assistant professor, University of
Colorado at Denver, for their advice and encouragement on
this project.
Also, I wish to express my thanks to thesis
committee, Dr. David W. Hubly, and Dr. John R. Mays, for
introducing me to the Fluid Mechanics and Finite Element
Method.
Finally, I would like to give my special thanks to
my parents, my wife and children, for their understanding
throughout the entire project.


CONTENTS
CHAPTER
I. INTRODUCTION ...................................... 1
Objective ....................................... 1
Scope of the study .............................. 2
II. LITERATURE REVIEW ................................. 3
Energy Losses in An incompressible Flow ...... 3
Numerical Modeling of An Incompressible Flow . 12
III. GOVERNING EQUATIONS OF FLOW MOTIONS .............. 16
Flow Motions and Descriptions .......... 16
Stream Function-Vorticity
Equations of Flow Motion .................... 19
Nondimensionalization
of the Governing Equations ................... 23
Boundary Conditions ............................ 25
IV. NUMERICAL MODELING BY FINITE ELEMENT SCHEME .... 29
Basic Concepts ................................ 29
Derivation oif Finite Element Equations
of Flow Governing Equations ......... 35
Element Shape Function ........................ 40
Numerical Initial and Boundary Conditions .... 49
Computational Flow Chart ...................... 53
V. CASE STUDIES AND RESULTS ......................... 57
Case One : Study of Flow Patterns ............ 63


vi
Case Two : Prediction of Energy Losses ....... 77
VI. CONCLUSIONS ...................................... 99
REFERENCE ................................ . ........ 102
APPENDICES
A. THE DERIVATION OF ELEMENT EQUATION BY
GAUSS-GREEN THEOREM ....................... 105
B. THE INTEGRATION OF SHAPE FUNCTIONS .............. 109
C. DEFINITION OF SYMBOLS ........................... 115
D. THE FORTRAN SOURCE CODE ......................... 117


TABLES
Table
1. Computational Flow Chart .......................... 55
2. Summary of the Flow Parameter Used in
the Case Studies ................................ 58
3. Computed Velocity Distribution for
Opening Ratio = 0.7 and Reynolds number = 50 .. 73
4. Comparison of Peak Vorticity....................... 76
5. Summary of Computational Time Intervals
and the Required Times for Flow to
Reach Steady State .............................. 78
6. Computed Velocity Distribution
in a Poiseuille Flow
for Reynolds Number = 50 ........................ 81
7. Computed Velocity Distribution
in a Poiseuille Flow
for Reynolds Number = 100 ....................... 81
8. Comparison Between the Computed Exit Velocity
Profile at X/B = 8.7 and The Analytical Profile
for Poiseuille Flow............................. 84
9. Comparison of Friction Loss Between
Computational and Analytical Predictions ........ 87
10. Comparison Between The Computed Exit Velocity
Profile and The Analytical Profile
for the cases with a contraction .............. 94
11. Distribution of Flow Pressure
with a Contraction .............................. 95
Comparison of Discharge Coefficients Between
the Computed Results and Experimental Data .... 97
12.


FIGURES
Figure
1. Illustration of Sudden Expansion ................ 7
2. Illustration of
the Flow Domain with Orifice Plate ........... 10
3. Rotational and Irrotational Flows
between Two Infinite Plates .................... 18
4. Boundary Conditions of Poiseuille Flow ........... 27
5. Illustration of One-Dimensional
Linear Finite Element .......................... 30
6. The Example of Eight-Node
Quadratic Quadrilateral Elements ............... 34
7. The Eight-Node Quadratic Quadrilateral
Isoparametric Element with
Natural Coordinates ............................ 43
8. Finite Element Net Work for
the Case with Opening Ratio of 1 ............. 59
9. Finite Element Net Work for
the Case with Opening Ratio of 0.7 ........... 60
10. Finite Element Net Work for
the Case with Opening Ratio of 0.5 ........... 61
11. Finite Element Net Work for
the Case with Opening Ratio of 0.3 ........... 62
12. The Steady State Distribution of Stream
Function and Contours of Vorticity for
Opening Ratio = 0.7 and
Reynolds Number = 50 .......................... 64
13. Illustration the Evolution of Flow Patterns
Changing From Potential Flow to Its
Steady State for Opening Ratio = 0.7
and Reynolds Number =50 ........................
66


ix
14. The Variation of Velocity Profile along the
Flow Direction for Opening Ratio = 0.7
and Reynolds number = 50 ..................... 72
15. The Distribution of the Wall Vorticity for
Opening Ratio = 0.7 and
Reynolds number = 50 ......................... 75
16. The Steady-State Velocity Profiles
In A Poiseuille Flow for
Reynolds Number = 50 ......................... 79
17. The Steady-State Velocity Profiles
In A Poiseuille Flow for
Reynolds Number = 100 ........................ 80
18. The Variation of Velocity Profile Along the
Flow Direction for Opening Ratio = l
and Reynolds number = 50 ..................... 82
19. Comparison Between The Computed Exit Velocity
Profile and The Analytical Profile
for Poiseuille Flow with X/B = 8.7 ........... 83
20. Comparison of Friction Loss Between
Computational and Analytical Predictions
for Poiseuille Flow with X/B = 8.7 ........... 88
21. The Variation of Pressure Along the Flow
Direction for Poiseuille Flow
with X/B = 8.7 ............................... 89
22. Comparison Between The Computed Exit Velocity
Profile and The Analytical Profile for
Opening Ratio = 0.7 .......................... 91
23. Comparison Between The Computed Exit Velocity
Profile and The Analytical Profile for
Opening Ratio = 0.5 .......................... 92
24. Comparison Between The Computed Exit Velocity
Profile and The Analytical Profile for
Opening Ratio = 0.3 .......................... 93
25. Computed Pressure Drops......................... 96
26. Comparison of Discharge Coefficients Between
Computed Results and Experimental Data ....... 98
27. Coordinate Transformation .........*............ Ill


CHAPTER I
INTRODUCTION
Objective
Energy dissipation is an important subject to many
professions. Examples include crude oil delivered in pipes
water supply systems, air flow through the larynx and blood
flow in arteries. Numerous studies have been conducted to
develop analytical and empirical equations to estimate
energy losses associated with the shear friction and eddy
motions. Recently, energy loss in the larynx flow has
attracted much attention because of the lack understanding
of flows near the glottis. The Denver Center for the
Performing Arts, in association with the University of Iowa
has been engaged in a laboratory investigation of this
phenomena. Since the energy dissipation in a flow is
related to both flow and fluid properties, laboratory
arrangements may not easily cover wide ranges of flow
variables such as Reynolds number. As a result, numerical
modeling is considered a possible helpful alternative.
The purpose of this study is to develop a numerical
model to analyze a two-dimensional pressure flow and to
predict the flow energy losses associated with the wall


2
friction and cross sectional contraction. In this work,
plate orifices are considered in order to establish the
finite element technique before applying the technique to
laryngeal flow problems.
Scope of the Study
After a review of numerical methods, the finite
element method was selected to convert the nonlinear
governing equations of unsteady flow into their numerical
forms and to include the complicated boundary conditions in
the numerical integrals. Solutions in terms of the stream
function, vorticity, velocity and pressure drop at each time
step were obtained from the numerical computations for
unsteady flow. The steady state solutions were achieved by
allowing the successive time-dependent solutions to
converge. It is found that for low Reynolds number flows,
comparisons between the predicted flow variables and
analytical solutions are in good agreement.
Although this developed numerical scheme is only for
two-dimensional flow, it is believed that the approach is
also suitable for a similar three-dimensional flow problem.


CHAPTER TWO
LITERATURE REVIEW
Energy Losses in an Incompressible Flow
The main objective of this study is to develop a
numerical model for the prediction of flow patterns and
energy losses in a two dimensional pressure flow. In
general, there are two kinds of energy losses involved in a
pressure flow. They are friction loss and minor loss.
Friction in a flow is caused by the non-uniform velocity
distribution and fluid viscosity. Friction results in
energy dissipation that converts mechanical energy into
thermal energy. This energy loss to the environment is
irreversible.
Friction loss is determined by the length of the
conduits, flow Reynolds number, and wall roughness. Using
dimensional analysis, the friction loss equation in a
pressure flow may be derived as follows (Daugherty,
Franzini, and Finnemore, 1985):
where, hT= the head loss due to friction, L= the length of
the conduit, R^ = hydraulic radius, V= cross sectional
average flow velocity, g = acceleration of gravity, and C^=


Darcy friction factor.
4
For a pipe flow, substituting Rh=D/4, in which D is
the pipe diameter, into Eq. 2.1, the Darcy-Weisbach equation
is obtained:
h
L
f
L
V
D 2g
(Eq. 2.2)
where D = the pipe diameter, and f = 4C^.
To determine the pipe friction factor f the pipe
flow experiments were conducted as early as the mid-1800s.
For high Reynolds number flow, flow in rough pipe resulted
in a higher energy loss than in smooth pipe. In 1933, J.
Nikuradse experimented using artificially roughened pipes
and obtained reliable information for friction factor
(Gerhart and Gross, 1985). In 1944, the Moody Chart for
Friction Factor f versus Reynolds number and relative
roughness was published and, since then, it has been widely
adopted by engineers.
To expand the friction loss, Hagen and Poiseuille
successfully derived analytical solutions in 1840 for
laminar flows in a circular pipe and between two plates.
They stated that for laminar flow in a circular pipe, the
loss of head in friction is proportional to the first power
of the cross sectional average velocity and is given by
IL L
h = 32------------= V (Eq. 2.3)
L <5 D2
where ii = absolute viscosity, and <: = the density of fluid.
By equating Eq. 2.2 and 2.3, we have


5
64 M>
f = ---------
Since CDV/n is flow Reynolds number,
expressed as
(Eq. 2.4)
Eq. 2.4 can be further
64
f = ------ (Eq. 2.5)
Re
in which Re= flow Reynolds number.
From Eq. 2.5, it is found that in laminar flow the
friction loss is only a function of Reynolds number and
independent of the roughness of the pipe wall. On the Moody
chart, the friction factor f for laminar flow is a straight
line.
This fact is not only true for a laminar pipe flow
but also for a laminar flow in a noncircular duct. Jones
(1976) expanded Eq. 2.4 to duct flow with the circular
diameter D replaced by the hydraulic radius, D^, and
conducted a series of experiments for various duct
configurations. He suggested that for the laminar flow in a
rectangular duct with width, a, and height, b, the friction
factor may be described by
k
f = ------ (Eq. 2.6)
Reh
where k =
64
2 11 b
3 24 a
(2
(Eq. 2.7)
Reh -
___!i_
and Dh = hydraulic diameter.


6
For the laminar flow between two infinite plates,
the value of variable a, in Eq. 2.7, is infinite. As a
result, k is equal to 96 when is equal to 4B (Blevins
1984, and Gerhart and Gross i985).
Losses due to the change of local pipe
configurations are termed minor losses. The flow in a pipe
system may need to pass through fittings, bends, or abrupt
changes in the pipe configuration. These kinds of changes
in local pipe geometries generate a great deal of disorderly
motion and result in dissipation of flow momentum. Whenever
the flow velocity is altered, eddy currents are induced and
a loss of energy in addition to friction loss is created.
Although minor losses are usually confined to a short length
of flow path, the effects may not disappear for a
considerable distance downstream.
The minor loss may be expressed as
v2
hT = K ----- (Eq. 2.8)
L 2g
where K = loss coefficient, which must be determined
experimentally for each situation.
The only loss coefficient which can be calculated by
simple analytical methods is that the axisymmetric sudden -
expansion as shown in Figure 1. This problem was first
treated by Borda and Carnot. Using continuity, linear
momentum, and energy equations and neglecting the shear
stress along the wall, the head loss for sudden expansion


Figure 1

Illustration of sudden expansion


8
can be expressed as
hL = ( 1 -
)2 (
ZL
2g
(Eq. 2.9)
where A1= the cross sectional area at position 1, A2= the
cross sectional area at position 2, and V1= the average
velocity at position 2.
Comparing Eq. 2.9 with Eq. 2.8 yields
A1 2
K = ( 1-----p- )z (Eq. 2.10)
A2
Although we may be uneasy about the assumptions of
the Borda-Carnot approach and feel that Eq. 2.10 is too
simple to evaluate the accurate energy loss for sudden
expansion, the laboratory experiments showed that the
difference between observed and predicted values of K in Eq.
2.10 were within a few percent (Vennard 1969). Based on
this fact, it can be assumed that minor losses may be
considered for practical purposes to occur in a much more
localized flow region.
As to the minor losses other than sudden expansion,
they all have to be determined by experiments. In 1885,
Weisbach followed the approach used by Borda-Carnot to
observe the energy losses of abrupt contraction and a fair -
estimate of loss coefficients was obtained. In 1952,
Gibson's studies on diffusers of conical form grasped the
trend of the change of loss coefficient, which are defined
as a function of cone angle and area ratio. In 1942,


9
Langhaar successfully calculated the loss of head at
entrance by using boundary-layer theory. As to the losses
in bends and elbows, there were studied by Beij and Hofmann
in 1933. Hofmann gave loss coefficients for wide ranges of
bend shape, roughness, and Reynolds number.
Nigro (1978) and Pao(1961) indicated that for a flow
with Reynolds number below 104, the discharge coefficient
becomes increasingly sensitive to Reynolds number, first
rising and then falling as Reynolds number decreases. The
same phenomenon as indicated by Negro and Pao was observed
by Scherer (1981) in his experiments on laryngeal flow for a
convergent-divergent glottis and opening ratios d/B less
than 0.05, where d = the glottal diameter, and B = the
subglottal (tracheal) lateral diameter.
The discharge coefficient for the pipe flow used by
Nigro is
Ud "
2 A p
h
[ i (-f-r ]
(Eg. 2.11)
where C = discharge coefficient, 11^= the mean flow velocity
through the orifice throat, "a = the cross sectional area of
orifice, A = the cross sectional area of the pipe, Ap = the
pressure drop between the flange taps as shown in Fig 2,
and = fluid density.
Considering a unit width of a two-dimensional flow,
Eq. 2.11 can further be simplified to evaluate the flow
velocity between two infinite plates.


Figure 2 Illustration of the flow domain
with orifice plate


11
ud =
2 A p
h
[ 1 (-f-r ]
(Eq. 2.12)
where d/B = the opening ratio between two plates, d = half
of the opening width at the orifice, B = half of the width
between the two plates.
As indicated by Nigro, for laminar flow through a
thin square-edged orifice with Reynolds number between 100
and 200, the discharge coefficient for flange taps is
approximately .70 to .80.
The values of loss coefficients obtained by the
pioneer researchers mentioned above were continually
modified by researchers. There are presently many handbooks
and technical papers published with experimental data
(Idel'chek 1960, and ASHRAE Handbook of Fundamentals).
Unfortunately, owing to different experiments conducted
under different experimental conditions, data values for
minor losses may be different for the same flow
configuration. Therefore, the most representative one is
difficult to determine. Furthermore, experiments for the
evaluation of minor losses are common for turbulent flows,
whereas the empirical data for laminar flows are less
common. These are the main reasons why this numerical study
for estimating the energy losses for a laminar flow is
needed.
In this study, the laminar flow patterns and energy
losses associated with the plate orifice as shown in Figure


12
2, are numerically calculated. Before the selection of
numerical methods, the following literature review was
performed.
Numerical Modeling of An Incompressible Flow
There are three approaches by which we can solve
fluid problems, namely
1. Experimental
2. Theoretical? and
3. Numerical
The experimental approach has the capability of
producing the most realistic answers for many flow problems.
However, its cost is the highest one among these three
methods.
In theoretical approaches, assumptions should be
made in order to make the problem tractable. Therefore, the
problems that can be solved by theoretical approaches are
restricted to simple boundary conditions and usually limited
to linear problems.
However, the shortcomings of experimental and
theoretical methods do not exist in numerical approach. In
the numerical approach, a limited number of assumptions are
made and a digital computer is used to solve the governing
eguations of a flow problem. There is no restriction to
flow linearity and even complicated geometry can be treated.
Furthermore, the computational costs are always much less


13
than that of laboratory equipment for experiments. It is
why the numerical method has become more important when the
computing facilities become more efficient and accessible.
Most people attribute the first definite work of
numerical approach to Richardson (1910) who introduced point
iterative schemes to numerically solve Laplace's equation.
His scheme used data available from the previous iteration
to update each value of the unknown and was called the
finite difference method (FDM). Richardson actually carried
out calculations for the stress distributions in a masonry
dam with great success.
Up to the late 1930's, FDM was used for solving flow
problems. Allen and Southwell (1955) applied relaxation
scheme to solve the incompressible, viscous flow around a
cylinder. Although they were not the first scientists to
solve these kinds of problems numerically, their results did
capture the flow phenomenon and had good accuracy. At the
same time, John von Neumann developed his method for
evaluating the stability of numerical methods for solving
time-dependent problems. Douglas and Rachford (1950)
developed an implicit method for solving parabolic and
elliptic partial differential equations. From then on, the
basic knowledge for solving fluid dynamics problems by FDM
was complete.
Before the finite element method (FEM) was fully
developed, FDM was used to solve flow problems. The


14
advantage of FEM is that its rate of convergence is much
faster than that of FDM.
The finite element method was originally used by
aircraft structural engineers in the 1950's to analyze a
large system of structural elements in the aircraft.
Usually, there are two ways to formulate element approximate
equations: one is the variational method, which has a close
relationship to the classical variational concepts of the
Rayleigh-Ritz method (Reyleigh 1877;Ritz 1909); the other
one is the weighted residual method, which is modeled after
the well-known method of Galerkin (1915). In the weighted
residual method, the error in the approximate satisfaction
of the governing equations is not set to zero, but instead
its integral with respect to the selected "weights" is
required to vanish.
The application of FEM to nonstructural problems
such as fluid flows was initiated by Zienkiewicz and Cheung
(1965). Pioneering research has employed weighted residual
methods to formulate approximate element equations. Oden
(1972) was among the first to derive the basic theoretical
analog for the Navier-Stokes equations. Baker (1971, 1973)
published his results for the analysis of incompressible
flow using the Galerkin method. Olson (1974) used a pseudo-
variational finite-element algorithm for solving the
biharmonic equation for stream functions. Lynn(1974)
utilized the least-squares methods to develop a finite-


15
element algorithm for laminar boundary-layer flow. In
1978, Backer successfully extended his studies from 2-
dimensional laminar flow to 3-dimensional turbulent flow.
From that time on, FEM was widely adopted for solving fluid
dynamics problems.
In this study, the numerical model for solving
laminar, incompressible flow will be derived using FEM
because of its flexibility of element shape and efficient
computation. Galerkin's method is used to weight residuals
in the approximate governing equations of flow. The
procedure of formulation will be discussed in Chapter Four.


CHAPTER THREE
GOVERNING EQUATIONS OF FLOW MOTIONS
Flow Motions and Descriptions
Using viscosity as a criterion, fluid may be
classified as real fluid and ideal fluid. For a real fluid,
flow motion results in shear forces between fluid layers.
But for an ideal flow, the shear force is assumed to be zero
everywhere. Although this is an idealized situation, it
does simplify the complexity of flow problems.
The shear stress in a real flow is termed internal
friction stress and results from the cohesion and momentum
interchange between fluid layers (Daugherty, Franzini and
Finnemore, 1985). The magnitude of shear stress is
determined by fluid viscosity and the flow velocity gradient
which is usually expressed by flow rotationality.
For a viscous flow between two infinite flat plates,
the relationship between viscosity and shear stress may be
expressed as
__ du
r = (i ---- C Eg 3.1)
dy
where T = shear stress, y. = the absolute viscosity, and
du/dy = the time rate of angular deformation or the velocity
y
gradient.
From the above equation, it is noted that whenever


17
and wherever a flow velocity gradient exits, shear stress
occurs correspondingly. Based on laboratory observations,
the fluid in direct contact with a solid boundary has the
same velocity as the boundary itself. Therefore, when the
boundary is fixed, the no-slip condition is applied, i.e.,
the velocity of the fluid on the solid surface must be zero.
As a result, the farther the fluid from the solid wall, the
faster the flow velocity can be. This fact explains the
reason why the existence of the wall results in flow
velocity gradient and shear stress. As shown in Figure
3(a), the cross sectional velocity profile between two
plates is parabolic for a laminar flow. This is a typical
motion of a low Reynolds number viscous flow. Considering
element 1 in Figure 3(a), the uneven velocity gradients in
the x- and y-directions cause a rotation. This kind of flow
motion is called rotational flow (Janna 1987). It means that
the shape and orientation of any fluid element change from
time to time and from one position to another as well.
When the flow is assumed to be ideal, as shown in
Figure 3(b), the flow velocity profile between two plates at
any cross section is a straight line. As a result, the shape
of a fluid element remains identical. It means that the
fluid element does not rotate. This is so called
irrotational flow (Janna 1987).
Applying the concept of control volume to a two-
dimensional, viscous flow moving on the xy-plane, the


I
18
'ZA'//////S////////////V//////////;//////JS///////////;//////,///,,
4- ^
7n7f/f/jj)//7Tn'nnnu/ui/ii////m/wwj/r/T777////M/wvM .<
(a) Rotational flow ^
////////// /// ////////////l/milLL
W ri
H. 4

(b) irrotational flow
Figure 3 Rotational and irrotational flows
between two infinite plates
\
V
\i


19
rotational velocity of a fluid element in the z axis may be
expressed as
1 dv ou
ez = ---- (------------) (Eq. 3.2)
2 dx 6y
in which, z= rotational velocity in the z-direction, u =
velocity component in the x-direction, and v = velocity
component in the y-direction.
The vorticity e used to express the rotational
characteristic of flow, is defined as two times the
rotational velocity, namely
6v 6u
ez = 20z -------------- (Eq- 3.3)
6x 6y
in which ez= vorticity in the z-direction.
Although Eq's 3.2 and 3.3 are derived from
mathematic concepts, they offer a feasible way to measure
flow motion. The following section will further explain how
to use vorticity to form the governing equations of the
motion.
Stream Function-Vorticity Equations
of Flow Motion
The basic governing equations for a two-dimensional
incompressible ,viscous flow can be described by the
principles of continuity and momentum.
Conservation of momentum states:
ou ou du dp d2u d2u
ot ox oy dx dx*4 dy44


20
and,
ou 6v 6v op
0 (--- + u --- + v ---) = - + p,
ot ox oy 6y
The continuity Equation states:
(Eq. 3.4)
62v 62v
ox2 6y2
(Eq. 3.5)
6u ov
---- + ----- = 0 (Eq. 3.6)
6x 6y
where There are three different formulations used to solve
Eq's 3.4 through 3.6 for flow problems. They are:
1. The stream function formulation combines all
governing equations together and expresses all flow
variables by stream function. As a result, a fourth order
partial differential equation of stream function needs to be
solved. Olson (1975) first applied this approach using
finite element method to solve the flow around a circular
cylinder. The results, including the separation phenomenon
behind the cylinder, were found to be accurate.
2. The velocity-pressure formulation directly solves
momentum and continuity equations with the pressure and
velocity as unknowns. Kawahara et. al. (1974) used this
approach to solve several steady state flow problems with
Reynolds number up to 1200.
3. The stream function-vorticity formulation uses


21
stream function, and pressure. The advantage of this approach is to reduce
unknowns from three variables, u, v, and p, to two
variables, approach to solve flows in a constricted channel with
Reynolds numbers less than 100.
Using the ip-e formulation, Eq's 3.4 and 3.5 are
reduced to a diffusion equation of vorticity and Eq. 3.6
becomes a poisson equation of the stream function. In this
study, the combination of ip and e is used to describe flow.
The following section further presents the detailed
derivations.
From the principle of conservation of mass, the
velocity components u and v are related to the stream
function by the following equations,
6 u = ---- (Eq. 3.7)
oy
v = -
dip
ox
(Eq. 3.8)
Substituting Eq's 3.7 and 3.8 into Eq. 3.3 yields
ip = e
(Eq. 3.9)
62 62
in which =
. 2 . '2
6x oy
For an irrotational flow, Eq. 3.9 is further simplified to
t tp = o
(Eq. 3.10)
Taking the derivative 6/oy of Eq. 3.1 and the


derivative -o/dx of Eq. 3.2 and adding the resultant
equations together, we obtain
22
M-
d
or
62e
dx'
d2e
oy
de
---- + u
ot
de
dx
de
+ v ----
dy
M* 2 De
--- e = ----- (Eq. 3.11)
d Dt
Eq. 3.11 is a two-dimensional diffusion equation of
vorticity. Eq's 3.9 and 3.11 represent two coupled second
order equations governing incompressible viscous flows. The
continuity equation is implicitly satisfied by the
relationship between stream function and velocity.
In order to predict energy loss, we need to
calculate the pressure distribution through out the flow
field. Taking the derivative of Eq. 3.4 with respect to x
and the derivative of Eq. 3.5 with respect to y and then
summing the resultant equations together, we have
du 2 ov 2 du ov
* P = 0 [ ( ---- ) + ( ----- ) + 2----------]
dx dy dy dx
(Eq. 3.12)
Substituting Eq. 3.6 into Eq. 3.12 yields
du dv du dv
2p = 2d [----------------------] (Eq. 3.13)
dx dy dy dx
Eq. 3.13 may further be expressed by stream function, namely


23
r2p = 2C [ -
6y dx
) ] (Eg. 3.14)
2
Eq's 3.9, 3.11, and 3.14 are governing equations
used in this study. It is noted that Eq's 3.9 and 3.14 take
the same form as the Poisson equation provided that the
terms on the right-handed side of each equation can be
treated as constants.
may become more generalized when they are expressed in
nondimensional forms. Without nondimensionalization,
solutions of a numerical model may only be suitable for the
specific conditions studied. In the process of
nondimensionalization, certain flow dominating parameters
like Reynolds number may be found.
flow quantities have to be chosen as characteristics of the
flow system (Daily and Harleman, 1966). For a Poiseuille
flow, the inflow velocity, U, may be chosen as the
characteristic velocity and a half of the distance between
two infinite plates, B, may be chosen as the characteristic
length. A set of dimensionless flow quantities may be
defined as follows:
Nondimensionalization of
the Governing Equations
In fluid mechanics, the solution of flow problems
To perform nondimensionalization, some reference
0
x
0
Y
x
y
B
B


24
u
0
t
0
0
u
U ,
tu
UB ,
0
V
0
P
V
u
p

eB
U
(Eq. 3.15)
Substitution of Eq. 3.15 into Eq. 3.9, the equation of
stream function becomes
62( (p ) 62(

6 ( x0)2 6 ( y0)2
or

2
( (Eq. 3.16)
Similarly, Eq 3.11 becomes
6(e) 0 6(e) 0 6(e)
6(t)
+ u
' / 0.
o(x )
+ V
6(y)
CUB
(
o2(e) 62(e)
'/ 0X2 A/ 0.2
o(x ) 6(y )
)
(Eq. 3.17)
The term of in the above equation is the reciprocal of
flow Reynolds number which delineates flow regimes, i.e.,
laminar or turbulent.
CUB
Re = ----- (Eq. 3.18)
Reynolds number physically represents the ratio of inertia
force of flow to viscous force. For a viscous flow between
infinite plates, gravity does not affect the flow pattern.


25
It is also obvious that capillarity and surface tension are
of no practical importance. Hence the significant forces
affecting flow motions are inertia force and viscous force.
Substituting Eg. 3.18 into Eg. 3.16 yields
1
Re
*2(e)
D(e )
D(t )
(Eg. 3.19)
Based on the similarity theorem, Eg's 3.17 and 3.19
shall depict the same flow dynamic behavior as long as the
boundary conditions and Reynolds number remain the same.
The dimensionless form of the pressure distribution
eguation is derived as follows:
. 2 . 0. i.2. 0. . 0. , . 0.
6 (p ) 6 (p ) 6(u ) 6(v )
(--------- + ---------) = -2 (----------------
a/ 0 .2 .. 0 .2 w 0. .. 0.
6(x ) 6(y ) o(y ) 6(x )
6(u) 6(v)
---------------)
6(x) o(y)
or
t2(P) = -2
. 0 N . 0.
6(u ) o(v )
. 0\ / 0 x
6(y ) 6(x )
'/ a/ o.
o(u ) 6(v )
6(x) 6(y)
(Eg. 3.20)
Comparing Eg. 3.20 with Eg.3.14, it can be seen that the
density <: of fluid is eliminated by nondimensionalization.
In the following sections, the flow problem will be
solved using the nondimensional eguations, and the
superscript 0's are neglected.
Boundary Conditions
When solving a partial differential eguation, the
boundary conditions play a determining role. In the study


26
of the flow between two infinite plates, as shown in Figure
2, flow is symmetric to its centerline, therefore, only a
half of the flow field is used to define the computational
domain. The inflow is assumed to have a uniform velocity
profile with a constant horizontal velocity U. The boundary
conditions associated with each governing equation of flow
are shown in Figure 4 and summarized below :
On the boundary AB :
Flow enters with a uniform velocity profile,
therefore, we have
v = 0
u = 1

e = 0 and
the value of pressure is set to be a selected datum
value such as 100.
On the boundary BC :
Applying the non-slip condition to this boundary, we
have

To estimate the vorticities generated from the wall,
it is necessary to apply the Taylor expansion to Eq.
3.3. The detailed derivation will be discussed
later.
On the centerline AD :
Based on the continuity and Eq. 3.15, the stream


p = Selected
datum value
Figure 4 Boundary conditions of Poiseuille flow
to


28
function at the centerline must be

Due to the occurrence of the maximum horizontal
velocity, we can write
6u/6y = 0, (Eq. 3.21)
Considering the symmetry of the flow about the
centerline, we conclude that
v = 0, and e = 0, on this boundary.
On the boundary DC :
This is the exit of the computational domain. Cheng
(1972) suggested to introduce a velocity profile
defined by a cubic equation of distance measured
from the wall to calculate stream functions aiong
the line DC, to ensure that Poiseuille flow is
generated. This suggestion forces the velocity
profile on the line DC to agree with the analytic
velocity profile, i.e., a parabolic curve. In this
study, it is found that as long as the computational
domain is long enough, the exit velocity profile
agrees with the analytic solution quite well.
Therefore, the condition that the vertical velocity
vanishes on the boundary DC is adopted and used in
the computation.


CHAPTER FOUR
NUMERICAL MODELING BY FINITE ELEMENT SCHEME
Basic Concepts
The finite element method is a numerical scheme and
computational procedure for obtaining approximate solutions
to many engineering problems (Segerlind, 1984). In general,
the application of the method takes the following five
steps:
1. Discretization and Selection of Element Configuration
This step involves dividing the entire computational
domain into a suitable number of "small" bodies called
finite elements. The shapes of the elements may be
triangular or rectangular, which is determined by the user.
2. Determination of the Approximation Eguation
For each element, it is reguired to have an
approximation equation for solutions in terms of the unknown
nodal values. For example, a one-dimensional linear element
AB, as shown in Figure 5, has an approximation equation
expressed by a linear equation.

The coefficients a^ and a2 can be determined by using the
nodal conditions, namely
6 II H 0) ft x = X1



30
Figure 5 ^Illustration of one-dimensional
linear finite element


31
Substituting Eq. 4.2 into Eq. 4.1 gives
X_-x x-X1

L 1 L Z
Eq. 4.3 is a typical form for approximation equations used
in the finite element method. The nodal values are
multiplied by linear functions of x, which are called shape
functions or interpolation functions. To be more
generalized, Eq. 4.3 may be rewritten as

where = the unknown nodal values, = shape functions.
In the above linear case, the shape functions become
N = , and N = ------- (Eq. 4.5)
1 h z h
The subscript i in Eq. 4.4 is an index that indicates how
many nodes are involved in an element. For a linear
triangular element, the value of n equals 3. As to a
quadratic quadrilateral element, the value of n equals 8,
3. Derivation of Element Equations
There are several ways to derive the element
equations in the finite element method such as the
variational method, and the weighted residual method. In
this study, Galerkin's method is used, which is a kind of
weighted residual method.
The weighted residual method is used to minimize the
residual resulting from an approximate solution that is
substituted into the differential equations governing a


32
field problem (Desaif1979). Usually, it requires that
J Wi(x)R(x) dx = 0 (Eq. 4.6)
J0
where R(x) = the residual, W^(x) the ith nodal weighting
function, and L = the length of integral.
There are a number of schemes developed for the
weighted residual method to determine the weighting
functions, among which are collocation, subdomain, least
squares and Galerkin's methods.
When Galerkin's method is used, Eq. 4.6 is evaluated
using the same functions for W^(x) as were used in the
approximate solution (Segerlind 1984). Thus, for any
element, the weightihg functions used for deriving element
equation are exactly the same as its shape functions. When
there are n nodes involved in an element, there are n shape
functions that can be used as weighting functions to form n
equations as Eq. 4.6 for solving n unknown nodal values
For a two-dimensional field problem, Eq 4.6 may be
rewritten in its matrix form
I [N ]T R(x) dA = 0 (i = 1 to n) (Eq. 4.7)
J A
T
where [N] = row vector containing the element shape
functions.
After integration, we have n simultaneous equations
with n unknowns of These simultaneous equations may be
expressed as


33
[K]{$) = {F} (Eq. 4.8)
where [K] = stiffness matrix or element property matrix, {$}
= the unknown nodal values, and {F} = force matrix or vector
of nodal forces, which may be a zero vector.
4. Assembly of Element Equations to Obtain Global Equations
and Introduction of Boundary Conditions
The final aim is to obtain the global equations that
approximately define the behavior of the entire system. As
shown in Figure 6, when the eight-node quadratic
quadrilateral elements are used ( for instance, the element
one has nodes 1, 10, 15, 16, 17, 11, 3, and 2), the element
equation can be written as
(1)
t K'*'] { = {
(1)
F ^1 ^ }
where the superscripts represent the element number, and
{ *
r
$
$
$
$
$
$
$
$
1
10
15
16?
17
11
3
2
The sequence of node numbers takes the
counterclockwise rule. Since the different elements have
different nodes, we have to, based on the node numbers,
incorporate all element equations into the global equation.'
After the boundary conditions are introduced into
the computation, the modified global equation similar to Eq.
4.8 may further be obtained.


Figure 6 The example of eight-node
quadratic quadrilateral elements


35
5. Computation and Solution Finding of the Global Equation
The global equations are a set of linear or
nonlinear simultaneous algebraic equations. They can be
solved using the well-known Gaussian elimination method to
obtain the value at each node.
In the following sections, a numerical model of two-
dimensional incompressible flow will be derived according to
the steps stated above.
Derivation of Finite Element Equations
of Flow Governing Equations
To formulate element equations, a general form of
the dimensional field equation can be used (Segerlind 1984).
o

D ----r + D ------ Gip + Q = 0 (Eq. 4.11)
X 6x^ y
Comparing Eq. 4.11 with Eq. 3.9, we conclude that for the
calculation of stream function, Dx = Dy = 1, G = 0, and Q =
e.
According to the Galerkin's method, the residual
integral equation for any element e is given by
6x 6y
(Eq. 4.12)
where [N] is the row vector containing the element shape
functions and is regarded as the weighting function in the
Galerkin's method.
To integrate the second-derivative terms in Eq.


36
4.12, the Green-Gauss theorem may be applied. After
integration, Eq. 4.12 becomes
{R(e)} = + [k(e)]{$(e)} {f(e)} = 0 (Eq. 4.13)
where
{I
tip
cos + D ------- sin )dr
Y 6y
(Eq. 4.14)
{k
D.
T
dfN]1
OX
6[N]
6x
+ D.
T
[N] A
6y
6[N]
oy
)dA
+
dA
(Eq. 4.15)
and
(Eq. 4.16)
As to two-dimensional time dependent problems, the
general form of their governing equations may be expressed
as
x 2
ox
+ a. o + q =
6t
y 2
1 oy
Eq. 4.17 can be rewritten as follows:
, , 6

where x oxz Y tyz
The residual integral is
(R(e)) = -| [N]T ( (Eq. 4.17)
(Eq. 4.18)
(Eq. 4.19)


37
Similarly, after integration
{R(e)} = {I(e)} + [k(e)
= 0
where {I^}, {k^e)},
expressed in Eg. 4.13,
we have
{$(e)j _ + [c^ ] }
(Eq. 4.20)
( g }
and {f1 } are the same as
[C(e)] = f [N]T[N] dA
JA
/ e)
and ;) = the nodal values of d

(Eq. 4.21)
{
(e) }T of {4(e)} is
6$. 6$. 6$
[ 1 2 n
dt ot dt
]
(Eq. 4.22)
The detailed derivation of Eq's 4.13 and 4.20 are presented
in Appendix A.
Next we apply the general forms to the flow problem.
The governing equations for a two-dimensional,
incompressible flow as derived in chapter two are:

£
1
Re
e
De
Dt
(Eq. 3.16)
(Eq. 3.19)
_ 6u 6v 6u dv
* p = 2 [-----------------------] (Eq. 3.20)
dx 6y dy dx
in which Eq's 3.16 and 3.19 are the two coupled governing
equations of flow motion. Eq. 3.19 may be further expressed
by stream function



38
1 2 6e b

- e = ------ +----------------------- (Eq. 4.23)
Re 6t 6y 6x 6x 6y
Because both

partial differential equation. To overcome the nonlinearity
of Eq. 4.23, some researchers (Chehg, 1972) suggested that
an unsteady flow problem can be regarded as linear in stream
function tp and vorticity e at each time step. We may
consider two solutions ( Vn £n ) at the nth time step and (
solutions are achieved by allowing the successive time-
dependent solutions to converge. Therefore, the governing
equations of flow motion become
tp 2 = e (Eq. 4.24)
Mi+1 n
and
1
Re
2e
n+1
6t
6y
6en
bx
dtp
n+1
bx
^fn
by
(Eq. 4.25)
It is noted that, at each time step, the in Eq. 4.25 are no longer unknowns. Thus Eq* 4.25 becomes
solvable numerically.
After the steady state is reached, the velocity
components u and v at every node can be obtained by taking
the first derivative of the stream function with respect to
y and x ^respectively. Then the velocity gradients in the
right-hand side of Eq. 3.20 can further be computed. As a
result, Eq, 3.20 may be regarded as a Poisson equation.


39
If the shape functions used for calculating
vorticities are the same as those used for stream functions,
.... !
applying Eq. 4.13 to Eq's 4.24 and 3.20, and Eq. 4.20 to Eq.
4.25, the element equation can be obtained as
{I(e)} + []Je)] = {f(e)}
1 ip 1 1 ip J l^n+l (Eq. 4.26)
] (cn+i> + [c(e)]{fe(n)> =
(Eq. 4.27)
+ (k^^ ] {p (e)} = {f} 1 P (Eq. 4.28)
(Eq. 4.29)
(1 (Eq. 4.30)
<^> = [ [N]t( j r op COS0 + 6p sin )dr
ox oy
(Eq. 4.31)
{k(.!}} r 6[N]T = ( 6(N] + T o[N]1 6[N] ) dA (Eq. 4.32)
A 6x ox 6y 6y
{k
1
Re
(k<^)
(Eq. 4.33)
{k

(Eq. 4.34)
I


40
N]T[N] dA (Eq. 4.35)
[N]T en) dA (Eq. 4.36)
[N]T ( -------------------------- )dA
oy 6x 6x 6y
(Eq. 4.37)
_ 6u 6v 6u ov
[N]1 (----------------------) dA
6x 6y 6y 6x
(Eq. 4.38)
Eq#s 4.26 through 4.38 are the element equations for solving
a two-dimensional, incompressible flow problem.
For an irrotational flow, the governing equation is
2

Then, the element equation is reduced to
] {tp^} = 0 (Eq. 4.26)
where {I^e^} and [k^e^] can be referred to Eq's 4.29 and
4.32 respectively.
In the following sections, the element shape
functions will further be discussed and determined.
Element Shape Function
To apply the finite element method to a flow
problem, we need to decide what kind of elements and element
equations should be used. The elements may be triangular,
rectangular, and other shapes. The element equations may be
[C(}] = [ [
JA
- JA <
(f'^) --2 J


41
linear, quadratic, or cubic.
There are two requirements that should be satisfied
in the determination of element shape function. The first
requirement is the continuity between adjacent elements. In
Eq. 4.6, R(x) represents the governing equation. Suppose
that R(x) contains derivatives up to the r-th order. To
avoid the integral becoming infinite, the compatibility
requirement (Huebner, 1975) requires that its (r-l)-th order
derivative must be continuous. This ensures that the
piecewise continuity is maintained at element interfaces.
In this study, the linear element shape function may be
chosen. Because the governing equations of flow motion
contain the second-derivative terms, a linear element with
continuity of first order can sufficiently satisfy this
requirement.
The second requirement, criterion of completeness,
results from the approximation theorem ( Zienkiewicz, 1977).
The finite element method has an implicit assumption that
convergence occurs as the size of each element decreases. To
ensure that this assumption is satisfied well, at least a
constant value can be obtained in a local element when the
size of any element tends to zero. Therefore, if a
derivative of order r exists in the governing equation, it
is necessary for the element shape function to be at least
of the order r so that, in the limit, such a constant value
can be acquired. Based on this criterion, a quadratic


42
element shape function was chosen.
It is noted that the criterion of completeness is
not as necessary condition as the requirement of
compatibility, for determining the Order of element shape
function. But, in this study, since a time-dependent
equation of vorticity (Eq. 3.19) needs to be dealt with, it
is better to use a quadratic element shape function to
ensure an efficient convergence. In addition, considering
the application of this numerical model to a flow problem
with curved boundary, the quadratic quadrilateral
isoparametric elements were selected and used in this study.
Isoparametric elements are useful in modeling a
field problem with the curved boundary and in grading a mesh
from coarse to fine ( Cook, 1981). They make it possible to
use nonrectangular quadrilateral elements for solving a
field problem with the curved boundary. The interpolation ,
equation for the eight-node quadratic quadrilateral
isoparametric element shown in Figure 7 is

1 2 J 4 D b / o
(Eq. 4.39)
with nodal conditions :
a = -1 and B = -1 at node 1/
a = 0 and B = -1 at node 2,
a = 1 and B = -1 at node 3,
a ~ 1 and B = 0 at node 4,
a = 1 and B = 1 at node 5,


43
Figure 7 The eight-node quadratic quadrilateral ,
isoparametric element with natural coordinates


44
a = 0 and B = 1 at node 6,
a = -1 and fl = 1 at node 7, and
ct - -1 and B = 0 at node 8.
in which a and B are natural coordinates.
The element shape functions can be obtained by applying
their basic characteristic to Eq. 4.39. The basic
characteristic of the shape functions is that each shape
function has a value of one at its own node and zero at the
other node. Applying this characteristic to Eq. 4.39 gives
a standard form as follows:

where = the unknown nodal values, and = shape
functions, in which
Nx = \(1-a)(1-fl)(-a-B-l)
N2 = ^(1-a2)(1-B)
N3 = *(l+a)(l-fl)(a-B-l)
N4 = \(1-B2)(1+a)
N5 = \(1+a)(1+fl)(a+fl-1)
N6 = \(1-a2)(1+B)
Ny = %(l-a)(1+B)(-a+6-1)
N8 = ^(1-B2)(1-a) (Eq. 4.41)
Applying Eq 4.41 to Eq's 4.29 through 4.38, the
integration can be performed. There are two types of
integrations:
1. The integrations occur on the boundary, such as


45
Eq's 4.29, 4.30, and 4,31. This kind of derivative boundary
condition will be discussed in section 4.4.
2. The integrations occur within the elements. They
can further be subdivided into two different types: the
integrations of the shape functions such as Eq's 4.35, and
4.36, and the integrations containing partial derivative
terms such as Eq's 4.32, 4.33, and 4k34, that involve the
first derivatives of the shape functions> and Eq's 4.37 and
4.38, that include the first derivatives of stream
functions, vorticities, and velocities.
For the purpose of integration, the change of
variables is required. For a two-dimensional integral, the
change of variables equation is
| | f(x,y)dA =
A
[ [ f[x(a,B),y(a,B)] | det[J] | dadB
J-1 J-1
(Eq. 4.42)
in which f(x,y) = the function expressed by global
coordinates, f[x(a,B),y(a,B)] the function expressed by
natural coordinates, and [J] = Jacobian matrix, which is
defined as
d(x,y)
[J] = ---------
6(a,B)
6x 6y
oa 6a
6x 6y
6B 6B
(Eq. 4.43)
and det[J] = the determinant of the Jacobian matrix.
Applying Eq. 4.43 to Eq. 4.35 yields


46
[C
[N(a,B)]T[N(a,B)]
det[J]| dadfi
(Eq. 4.44)
Similarly, Eq. 4.36 becomes
{f
a,B)]T[N(a,B)] {e })
|det[J]| dadfi
(Eq. 4.45)
Eq's 4.44 and 4.45 can be evaluated by Newton-Cotes or
Gauss-Legendre quadrature method.
For those integrals containing the partial
derivatives of the shape functions, we need to change the
local derivatives in terms of a and B to global derivatives.
The relationships between local derivatives and global
derivatives are
6N.
dni
6y
(a,B)
-1
= [J(a,B)]
where [J] = Jacobian matrix, and [J]
6Ni
6a
bN^
""6B
-1
(a,B)
(Eq. 4.46)
(a,B)
= inverse of [J].
However, in Eq. 4.46 an explicit form of [J]-
cannot be obtained (Zienkiewicz,1977). Numerical schemes may
be used to evaluate element integrals. In this study, Gauss-
Legendre quadrature is adopted. The Gauss-Legendre
quadrature selects the sampling points between the integral
intervals to achieve the best accuracy. For double
integrals, it states that


47
ll pi in n
R= f(a,B) dadB ~ S S [W.W. f(q.,B.)]
J -1J -1 i=l j=i 1 ? J 1
(Eq. 4.47)
where m is the number of the selected sampling points with
respect to B, which can be determined by equating (2m-l) to
the highest power of B existing in the function f(a,B). If
the resulting m is not an integer, the next larger integer
is used.
n is the number of the selected sampling points with
respect to a, which can be determined the same way as that
used for determining m.
is the weighting factor for the sampling points with
respect to B, and Wj is the weighting factor for the
sampling points with respect to a.
The detailed derivation of Eq's 4.42 and 4.46 and
the further explanations of Gauss-Legendre quadrature are
presented in Appendix B.
How many sampling points should be chosen when using
Gauss-Legendre quadrature to maintain the accuracy ? In this
study, referred to Eq. 4.44, the highest power of variables
. T 4 4
appearing in [N] [N] are a and B Therefore, the value of
3 is used for both m and n. It means that the order of
quadrature is 3 and there are 3x3 points used for
computation. Cook (1974) suggested for plane elements,
either linear or quadratic, the 2x2 points are usually the
best, except for a markedly nonrectangular or elongated


48
quadratic element, where 3x3 points may give better
accuracy. In this study, the numerical model developed may
further be applied to the analysis of laryngeal flow.
Considering the irregularity of larynx geometries, the
nonrectangular and elongated elements were expected.
Therefore, 3x3 points are considered the most reasonable
choice.
Applying Eq. 4.47 to Eq. 4.44, it gives
3 3
[C^e)] = 2 E {W W [N(a ,fl )]T[N(a ,B )]} |det[J]|
6 i=i j=i 1 i i J 1
(Eq. 4.48)
The local coordinates of each sampling point and
its corresponding weighting factor in the above equation are
presented in Appendix B. After the coordinates of each
sampling point were determined, det[J] at each point can be
obtained using Eq. B.3. Then Eq. 4.48 can further be
evaluated numerically. The same procedure is applied to Eq.
4.45.
As to Eq 4.32, which contains the first derivatives
of the shape functions, with the aid of Eq's 4.46 and 4.47,
it becomes
, . 3 3
{k(e)} = E E W. W
* i=l j=l 1 ^
d[N(ot.,S.)] o[N(a.,B.)]
( __________J -1- _________J_______
-1 6x 6x
6[N(a.,fl.)] d[N(a.,B.)]
+ ------Ji---- -------J----- ) | det [ J ] | .
6y oy
(Eq. 4.49)
Similar approach to the forementioned can be applied to Eq#s


49
4.33 and 4.34.
The term of [ (6 in Eq 4.37, may be evaluated by applying Eq. 4.46 to the
center-point of each element. For any element,
[ (6 By applying Gauss-Legendre quadrature, Eq. 4.37 becomes
(eK _
3 3
(f V > = s s {W.W. [N(a ,B.)] } det[J] Cc
i=l i=l 1 3 3 1
where Cc = (
6

mi+1 n
dip ,, be
^n+1 n
bx
6y
(Eq. 4.50)
) evaluated at
by bx
the center of each element.
Eq. 4.38 can be solved by the same way as used for
Eq. 4.37. The term of (?u ov
6u S- )' in Eq- 4-37
ox by by bx
are also evaluated at the center of each element.
Numerical Initial and Boundary Conditions
Since Eq. 3.19 is a time-dependent equation, the
initial conditions have to be determined. In this study, the
initial flow field was considered to be a potential flow
(ideal flow) with a uniform inflow of velocity U at the
entrance. As shown in Figure 4, the initial vorticity is
zero everywhere in the computational domain. To determine
the distribution of stream function, the governing equations
for the irrotational flow were used.
There are two kinds of boundary conditions involved
in this study. They are, Dirichlet boundary conditions,


50
that are the values of (p's specified on the boundaries, and
Neumann boundary conditions, that are the values of
derivatives 6(p/6x and 6(p/6y specified on the boundaries.
Recall that in the derivation of element equations, the
Neumann boundary conditions have already included as a part
of element equations (Eq's 4.29, 4.30, 4.31) and take the
following form
( D,
6 cose + d
6 sine ) ,
(Eq. 4.51)
" 6x * 6y
where e = the angle to the outward normal to the boundaries
of each element.
If Dx = Dy = D, then Eq. 4.51 becomes
d o
6n
where 6(p/6n is the derivative normal to the boundary.
To evaluate Eq. 4.29, in which the Neumann boundary
conditions of the stream functions were involved, the
derivative boundary conditions of stream functions in Figure
4 are listed below:
along line AB, 6

along line BC, 6(p/6n = 0,
along line AD, 6(p/6n = 0, and 9 = 90
along line DC, 6

(Eq. 4.52)
Substituting Eq. 4.52 into Eq. 4.29 yields
(Eq. 4.53)
Eq. 4.53 implies that the Neumann boundary conditions fbr
>


51
stream functions are self-consistent in the element
equations. If boundaries AB and DC in Figure 4 are vertical
straight lines, the contribution of the derivative boundary
conditions of the stream functions to the element equations
automatically vanishes. Therefore, for the stream function,
only Dirichlet boundary conditions need to be considered. As
a result, Eq. 4.29 can be removed from the element
equations.
The Dirichlet boundary conditions of the stream
functions as shown in Figure 4 are listed below:
Along line BC,

along line AD, tp = 1.
How to impose Dirichlet boundary conditions by
modifying the global equations can be found in Reference 7,
9, and 29.
To consider the numerical boundary conditions of
vorticity, a linear combination of e and its normal
derivative is required to be specified throughout the entire
computational domain (Backer, 1983). On the fixed wall, BC,
the value of vorticity is a priori unknown. To evaluate the
wall vorticity, the no-slip boundary condition alone is
insufficient. From Eq. 4.24, at a point (xQ,y0) on the wall,
the vorticity may be calculated by
e
n
d2 6B
n+-l
2
(
y0 >
(Eq. 4.54)
where B is the coordinate normal to the wall. To evaluate


52
2 2
6 (Xi,yi) which is located a small distance from the wall
along the fl direction.
n+i B2 2n+1
Vn^xl/yl^ ,pn+l^xO/yO^+fl n+11 0'* 0
6B 0, 0 2 6B
(Eq. 4.55)
Since the no-slip condition indicates that
6 oB
( x0 Y0 ) = 0
we have
2
n (Eq. 4.56)
Wall vorticity is then expressed in terms of the stream
functions on the wall and the node at a small distance from
the wall. Because & is the coordinate normal to the wall,
the construction of the grids for computing wall vorticity
do need special attention.
As to the derivative boundary condition of
vorticity, the values of 6en+1/6x along boundaries AB, BC,
CD and den+1/6y along boundary BC are undeterminable. Eq.
4.30 may be considered as a negligible term which is further
ignored from the global equations. In this study, it is
shown that neglecting the derivative boundary conditions of
vorticity can still give acceptable results which agree to
analytic solution well. The same situation may be applied


53
to the derivative boundary conditions of pressure. By
adding cos times (Eq.3.4) to sin times (Eq.3.5) and using
the appropriate boundary conditions, we have
r 6 u 6 u
[cos( j + 2 )+sin(
6x^ oy
(Eq. 4.57)
Substituting Eq. 4.57 into Eq* 4.31, it is found that the
integrations of the second derivative of velocity along the
element boundaries are difficult to evaluate. Like Eq. 4.30,
Eq. 4.31 is considered as a numerically negligible term and
ignored from the global equation.
a two-dimensional, incompressible, viscous flow between two
infinite plates, the initial conditions are determined using
the governing equations for ideal flow and then Eq's 4.26
and 4.27 are used to calculate stream functions and
vorticities for successive time steps. Neglecting the
derivative boundary conditions, Eq's 4.26 and 4.27 are
reduced to
Computational Flow Chart
To develop a digital computer program for analyzing
(Eq. 4.58)
(Eq. 4.59)
in which
(Eq. 4.60)


54
Substituting Eq. 4.60 into Eq. 4.59, the recurrence relation
is
{ itfk^mc^] H£n+1} = Cc(6)](e(n),+ At
(Eq. 4.61)
At each time step, Eq. 4.58 is first used to calculate Substituting new values of corresponding wall vorticities can be obtained. Using Eq.
4.61, the vorticity can be updated for the next
computational time step.
The criterion for determining the steady state was
defined to be when the largest value of 6e/ot in the entire
computational domain was less than 0.04 (Backer, 1983).
After the steady state is reached, using Eq. 4.28, the
corresponding pressure distribution is further computed.
Comparing pressure drops between downstream and upstream
flows, we can determine the energy dissipation.
A flow chart illustrated the computational
procedures is presented in Table 1. The source code and
FORTRAN program can be found in Appendix D.


FLOW CHART A
(Continued)
Table 1
Computational Flow Chart


Flow Chart B
End


CHAPTER FIVE
CASE STUDIES AND RESULTS
There are four different contraction geometries used
in this study. They can be distinguished by the opening
ratio, r, which is defined as
r = -f (Eq. 5.1)
where d = a half width of the opening, and B = a half width
of the distance between two plates. When r = 1, flow does
not have any contraction. Results of these cases present
predictions of friction losses only. When r < 1, flow is
subject to minor losses due to contraction. The opening
ratios used in this study were 0.3* 0.5, 0.7, and 1.0.
Table 2 summarizes the total number of nodal points,
elements, and Reynolds numbers used in the studies of these
four cases. The corresponding element network for each case
is shown in Figures 8 through 11. Since the abrupt changes
of stream function and vorticity are expected to occur near
the orifice, those cases with a smaller opening should use a
larger number of smaller elements around the orifice. Thus,
the number of nodal points and elements around the orifice
used in the case with an opening ratio of 0.3 has as nearly
1.5 times nodal points and elements as used in the case with
an opening ratio of 0.7.
The total length of the computational domain also


Table 2 Summary of the Flow Parameters used in the Case Studies.
Opening ratio Reynolds number Number of nodal points Number of element Length L/B Separation eddy
1.0 50 413 120 8.7 No
1.0 75 413 120 8.7 No
1.0 100 413 120 8.7 No
1.0 150 413 120 8.7 No
1.0 200 413 120 8.7 No
0.7 50 235 64 6.5 Yes
0.7 100 235 64 6.5 Yes
0.7 200 235 64 6.5 Yes
0.5 50 309 86 8.7 Yes
0.5 100 309 86 8.7 Yes
0.5 200 309 86 8.7 Yes
0.3 50 305 84 9.3 Yes
0.3 100 305 84 9.3 Yes
0.3 150 305 84 9.3 Yes
Ul
00


FINITE ELEMENT GRID NETWORK FDR R=1,0
(Solid Wall) <8.7, 0.0)
X/B
Figure 8 Finite element grid net work for
the case with opening ratio of 1
ui


FINITE ELEMENT GRID NETWORK FOR R=0.?'
(0.0, 1.0)
y/B
(Centerline of Flou)
(6.5, 1.0)

\ /

(0.0 ,0.0)
Orifice
(Solid Wall) (6.5, 0.0)
X/B
Figure 9 Finite element grid net work for
the case with opening ratio of 0.7
&
o


FINITE ELEMENT GRID NETWORK FOR R=0.5
(O.O, 1.0)
y/B
(CENTERLINE OF FLOW)
(8.7, 1.0)




(0.0 ,0.0)
RIFICE
(SOLID WALL)
(8.7, 0.0)
X/B
Figure 10 Finite element grid net work for *
the case with opening ratio of 0.5


FINITE ELEMENT GRID NETVDRK FDR R=0.3
>-<>' 1-0) .(CENTERLINE OF FLOW)
(9.3, 1.0)
y/B
(0.0
\ \ y 7 z

rO.O) Orifice- (SOLID WALL)
x/B
Figure 11 Finite element grid net work for
the case with opening,ratio of 0.3
7

<
62


63
needs special attention. After the flow passing the
orifice, it takes a longer distance for the flow with more
contraction to return to its parabolic velocity profile
than that with a smaller contraction. When the
computational domain downstream of the orifice is not long
enough, the numerical procedure may become unstable and
results are not reliable. In this study, it was found that
the range of computational domain, in terms of L/B varied
from 6.5 to 9.3. Predictions using these computational
domains give good agreements to the laboratory and
analytical results.
Predictions and comparisons for each case study are
presented as follow:
Case One : Study of Flow Patterns
As may be expected, eddies are generated behind the
orifice plate. Table 2 indicates that as long as the
opening ratio is less than unity, eddies always occur. For
high Reynolds number flows, eddies immediately downstream of
the orifice plate may further shed out. Figure 12 presents
steady state flow with an opening ratio of 0.7 and Reynolds
number of 50. The computational time interval, At, was
0.01. It took 640 time steps, i.e. t = 6.4, for flow to
reach the steady state. The upper half of Figure 12 shows
the distribution of stream function and the lower half
presents the contours of vorticity. It is noted that a
reverse flow occurs immediately downstream of the


Figure 12 The steady state distribution of stream function
and contours of vorticity for Opening Ratio =0,7
and Reynolds Number = 50


65
orifice plate and the highest vorticity is produced near
the orifice.
On this figure, the standing eddies are clearly
shown. The general flow patterns predicted in this study
are very similar to that reported by White (1974), Chung
(1978), and, Kwon and Pletcher (1981), although the
numerical approaches were different.
Evolutions of flow patterns with respect to time
changing from potential to steady state flow are shown in
Figure 13. It can be seen that the eddies are growing in
size as time goes on. Using the steady state criterion
mentioned in Chapter four, flow becomes steady when t =
6.4. However, it is found that when t = 3.2, the variations
of streamline patterns have already been numerically
negligible. This may imply that the criterion for steady
state suggested by Baker (1983) seems to be conservative.
Due to the inadequate accuracy of software for plotting
streamlines and the coarseness of elements used, some flow
patterns in Figure 13 are not shown to be separating from
the edge of the orifice.
Figure 14 displays the plot of velocity profiles
along the flow direction at x/B = 0., 1.1, 2.2, 2.5, 2.9,
4.7, and 6.5. Near the entrance, the velocity profile shows
that the flow has not been fully developed. Behind the
orifice, negative velocity occurs. This means that eddy
motions in the area of wake produce a reverse flow. Table 3


66
t o.o



_ ^ .


'I'll l 1 1 l 1 l 'I ' r r~ti | i i i i [ rm | i i i11 | i r
, 0 1 2 3 H 5 X/B G





1 f-1 1 | i i i i | r i | i i i i | i i i i | i i i i | i i
0 1 -2 3 H- 5 G
X/B
t 0.6
Figure 13 Illustration the evolution of flow patterns
changing from potential flow to its steady state
for Opening Ratio = 0.7 and Reynolds Number = 50


67
0.... 1 2 ' 3 5 X/B 6
v (continued)


6a
l
y/B
0
t 2.0

0
i i~r~f~[~ i i i'i | r i i ~t~| i i ~i r | i itt
1 2 3 . H
iir
S.
X/B
i"|T r
G
(continued)


69
(continued)


70
(
(contInued)


71


Figure 14 The variation of velocity profile along
the flow direction for Opening Ratio = 0.7
and Reynolds number = 50
-j
ts>


73
Table 3 Computed Velocity Distribution
For opening ratio= 0.7, and Reynolds number= 50
0.00 1.10 X/B 2.50 2.90 4.70 6.50 2.20
Distance Velocity Distribution At Orifice
y/B u/Uo u/Uo
0.00 1.00 0.00 0.00 0.00 0.00 0.00
0.05 1.00 0.18 -0.09 -0.08 0.08 0.14
0.10 1.00 0.36 -0.17 -0.13 0.15 0.27
0.20 1.00 0.65 0.03 0.05 0.36 0.54
0.30 1.00 0.92 0.33 0.28 0.58 0.77 0.00
0.45 1.00 1.11 1.03 0.92 0.96 1.06 1.27
0.60 1.00 1.27 1.67 1.61 1.37 1.29 1.50
0.80 1.00 1.30 1.74 1.80 1.58 1.45 1.57
1.00 1.00 1.33 1.80 1.96 1.77 1.57 1.64
Exit Profile Analytical Profile
Distance Y/B Velocity Profile u/Uc u/Uc
0.00 0.00 0.00
0.05 0.09 0.10
0.10 0.18 0.19
0.20 0.34 0.36
0.30 0.49 0.51
0.45 0.68 0.70
0.60 0.82 0.84
0.80 0.92 0.96
1.00 1.00 1.00
Note: X= horizontal distance. y= distance from solid wall.
B= half width between two plates.
Uo= entrance velocity. Uc= centerline velocity


74
presents velocity profiles at different locations and a
comparison between the exit velocity profile and the
analytic parabolic velocity profile. It is noticed that
although the existence of the orifice does disturb velocity
profile, flow still returns to Poiseuille profile far
downstream. In this study, the velocity profile for
Poiseuille flow can be expressed as (Gerhart and Gross,
1895)
u/Uc = 1 ( 1-y2)
where u= velocity in x-direction, Uc= u at centerline, y=
distance from the solid wall.
The wall vorticity is also worthy of studying
because the wall shear stress is directly proportional to
its value. Figure 15 depicts the distribution of the wall
vorticity for the flow with an opening ratio of 0.7 and
Reynolds number of 50. It shows that the peak wall
vorticity occurs at the contraction section and the wall
shear stress changes its direction corresponding to the flow
direction near the wall. Table 4 presents the values of
peak vorticity for various cases under different Reynolds
numbers. It is noted that for the case with an opening
ratio = 1, the peak vorticity, highest shear stress, occurs
at the entrance because the flow experiences the steepest
velocity gradient at the beginning of the entrance length
and its value is proportional to the Reynolds numbers. On
the contrary, for those cases with an opening ratio less


Woll Vorticity
22
X/B
Figure 15 The distribution of the wall vorticity for
Opening Ratio = 0.7 and Reynolds number *= 50
Ul


Table 4 Comparison of Peak Vorticity
The Values
Opening ratio Reynolds number of Peak Vorticity Occuring Position
1.0 50 16.62 Near Entrance
1.0 75 17.04 Near Entrance
1.0 100 17.36 Near Entrance
1.0 150 17.84 Near Entrance
1.0 200 18.19 Near Entrance
0.7 50 20.78 Near Orifice
0.7 100 23.49 Near Orifice
0.7 200 26.43 Near Orifice
0.5 50 34.15 Near Orifice
0.5 100 38.98 Near Orifice
0.5 200 42.40 Near Orifice
0.3 50 60.89 Near Orifice
0.3 100 66.80 Near Orifice
0.3 150 73.77 Near Orifice


77
than one, the highest shear stress occurs at the orifice
and their values are a function of both Reynolds number and
the opening ratio. With the same Reynolds number, the value
of peak vorticity with an opening ratio of 0.3 is as nearly
three times as that with an opening ratio of 0.7.
Computational time intervals and the required times
for flow to reach the steady state are presented in Table 5.
It is noted that the higher the Reynolds number is, the more
time steps are required.
Case Two : Prediction of Energy Losses
Energy dissipation in a pressure flow includes
friction loss and minor loss. In this portion, flows are
studied with and without an orifice. When the opening ratio
is unity, the energy loss is essentially a friction loss.
Figures 16 and 17 present the steady-state velocity profiles
at various locations with Reynolds number equal to 50, and
100 respectively. The corresponding data can be found in
Tables 6 and 7. Flow starts with a uniform velocity profile
at the entrance and gradually changes to a fully developed
flow through the entrance length. Figure 18 depicts the
variations of velocity profiles along the flow direction
with Reynolds number of 50. Figure 19 and Table 8 show the
comparison between the exit velocity profiles under Reynolds
numbers equal to 50, 75, 100, 150, and 200 with the analytic
velocity profile of a fully developed flow.
It is noted that the exit velocity profiles with


Table 5 Summary of Computational Time Intervals and the
Required Times for Flow to Reach Steady State.
Time Interval Time Steps
Opening Reynolds Used for for Reaching
ratio number Computation Steady State
1.0 50 0.0100 960
1.0 75 0.0100 1200
1.0 100 0.0100 1200
1.0 150 0.0100 1900
1.0 200 0.0050 2400
0.7 50 0.0200 320
0.7 100 0.0100 780
0.7 200 0.0050 2015
0.5 50 0.0100 633
0.5 100 0.0100 754
0.5 200 0.0050 1992
0.3 50 0.0100 567
0.3 100 0.0100 737
0.3 150 0.0025 3376


Velocity u/Uo


Velocity u/Uo
Figure 17 The steady-state velocity profiles in
a Poiseuille flow for Reynolds Number = 100
03
o


81
Table 6 Computed Velocity Distribution in a Poiseuille Flow
Reynolds Number= 50.
X/B Exit Analytical
0.00 2.00 3.90 6.30 8.70 Profile Profile
Distance Velocity distribution Velocity Profile
Y/B u/Uo u/Uc u/Uc
0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00
0.05 1.00 0.18 0.16 0.15 0.15 0.15 0.98
0.10 1.00 0.35 0.32 0.30 0.30 0.30 0.19
0.15 1.00 0.49 0.45 0.43 0.42 0.42 0.28
0.20 1.00 0.64 0.59 0.57 0.56 0.56 0.36
0.40 1.00 1.07 1.01 0.99 0.98 0.98 0.64
0.60 1.00 1.27 1.27 1.27 1.27 1.27 0.84
0.80 1.00 1.33 1.39 1.43 1.44 1.44 0.96
1.00 1.00 1.33 1.43 1.48 1.49 1.49 1.00
Table 7 Computed Velocity Distribution in a Poiseuille Flow
Reynolds Number= 100
X/B Exit Analytical
0.00 2.00 3.90 6.30 8.70 Profile Profile
Distance Velocity Distribution Velocity Profile
y/B U/UO u/Uc u/Uc
0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00
0.05 1.00 0.21 0.18 0.17 0.16 0.11 0.01
0.10 1.00 0.40 0.36 0.33 0.32 0.23 0.19
0.15 1.00 0.56 0.50 0.47 0.46 0.32 0.28
0.20 1.00 0.74 0.66 0.61 0.60 0.42 0.36
0.40 1.00 1,13 1.08 1.03 1.02 0.72 0.64
0.60 1.00 1.24 1.27 1.27 1.27 0.89 0.84
0.80 1.00 1.24 1.32 1.37 1.39 0.98 0.96
1.00 1.00 1.24 1.32 1.39 1.42 1.00 1.00
. __
Note: Uo = entrance veloctiy. Uc= centerline velocity.
X = horizontal distance,
y = distance from the plate.
B = half width between two plates.


Figure 18 The variation of velocity profile along
the flow direction for Opening Ratio = l
and Reynolds number = 50


Velocity u/Uc
Figure 19 Comparison between the computed exit velocity
profile and the analytical profile for
Poiseuille flow with X/B = 8,7
t '
83


Table 8 Comparison between the Computed Exit Velocity
Profile at x/B = 8.7 and the Analytical
Velocity Profile for Poiseuille Flow
Reynolds Number 50 75 100 150 200 Analytical Velocity Profile
Distance Y/B Exit Velocity u/Uc at X/B=8.7 u/Uc
0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.05 0.10 0.11 0.11 0.13 0.14 0.01
0.10 0.20 0.21 0.23 0.24 0.27 0.19
0.15 0.28 0.29 0.32 0.35 0.38 0.28
0.20 0.38 0.38 0.42 0.45 0.51 0.36
0.40 0.66 0.67 0.72 0.76 0.82 0.64
0.60 0.85 0.86 0.89 0.92 0.96 0.84
0.80 0.97 0.96 0.98 0.98 0.99 0.96
1.00 1.00 1.00 1.00 1.00 1.00 1.00
Note: Uo = entrance veloctiy. Uc= centerline velocity.
X = horizontal distance,
y = distance from the plate.
B = half width between two plates.


85
Reynolds number less than 100 give a better agreement to
the exact solution. The discrepancy between the predicted
exit velocity profiles and analytic solution for high
Reynolds number flows may be due to the length of
computational flow domain. When it is not long enough, flow
may not become
fully developed at the exit.
Theoretically, an infinite distance is required for
a flow with uniform entrance velocity to become fully
developed, but it has been proved both by theory and
observation that for a pipe flow, the maximum velocity at
the centerline of the pipe may reach 99 percent of its
ultimate value in the distance determined as follows:
L = 0.058 RePD (Eq. 5.2)
where Rep is the Reynolds number for pipe flow and defined
as
UD
ReP = -----
in which U = average flow velocity, and D = the diameter of
the pipe (Daugherty, Franzini and Finnemore).
Applying Eq 2.7 to Eq. 5.2 and using Reynolds number
determined by hydraulic radius we have
L = .348 Re B (Eq. 5.3)
For Re = 100, the required length to reach 99% of fully
developed flow is 35 B. This is impractical to adopt this
computational domain because of the elongation of


86
computational time. After many tests, the length of 8.7 B
is used in this study based on the judgement of the
agreement between the predicted exit velocity profile and
the exact velocity profile for flows with Reynolds number
less than 100.
As to the prediction of friction loss, computational
domain of L/B=8.7 can still give fair accurate predictions
even for the flows with higher Reynolds number. Numerical
predictions of friction loss expressed by flow pressure
drop are presented in Table 9 and Figure 20. The analytical
pressure drop in Table 9 was computed by Eq 5.3 that can be
derived by applying Eq. 2.7 to Eq. 2.1 and using hydraulic
radius for Reynolds number.
3 L
a p =---------- (Eq. 5.3)
Re B
where L/B = 8.7 in this study. It can be seen that the
computed friction loss is in a good agreement to the
analytical solutions.
Figure 21 presents the variations of flow pressure
along the flow direction with various Reynolds numbers. It
shows that independent of Reynolds number, kinetic effects
dominate over the first 30 percent of the entrance length.
Thereafter, viscous effects gradually increase in importance
and become the sole influence at the end of the duct, where
the pressure gradient is a constant.
To predict the energy losses due to the contraction,


Table 9 Comparison of Friction Loss between
Computational and Analytical Predictions
Reynold Number Computed Pressure Analytical Drop Pressure Drop
50 0.50 0.52
75 0.36 0.35
100 0.25 0.26
150 0.18 0.17
200 0.14 0.13
00
N]


Friction
Reynolds Number
Computed + Analytical
Figure 20 Comparison of friction loss between
computational and analytical predictions
for Poiseuille flow with X/B =8.7
oo
oo


Pressure
Distance x/B
Figure 21 The variation of pressure along the flow
direction for Poiseuille flow with X/B = 8.7


90
opening ratios of 0.3, 0.5 and 0.7 are used. The exit
velocity profiles for different cases are shown in Figures
22, 23, and 24 and Table 10. It seems that the
computational domain was not long enough for the flow with r
= 0.3. Table 11 presents the computed pressure drop. Figure
25 plots the computed pressure drops for the flows with the
opening ratios of 1, 0.7, 0.5, and 0.3.
Using Eg. 2.11, the orifice pressure drop can be
converted to discharge coefficient which can further be
compared with the experimental data observed by Nigro (1978)
for pipe flow with a circular orifice. It is noted that the
orifice pressure drop measured by Nigro is the pressure
difference between two flange taps that were located at one
inch on the both sides of circular orifice. In this study,
the nearest nodal points at one inch on both sides of the
plate orifice were chosen to calculate orifice pressure
drops.
Comparison between the predicted pressure drop and
experimental data is shown in Table 12 and Figure 26.
Although energy dissipation of a plate orifice may be
somewhat different from a circular orifice, these two sets
of pressure drops capture a similar trend and even numerical
deviations are fairly small.


Velocity u/Uo
Distance y/B
Re = 50 + Re=100 O Re=200 A Analytical
Figure 22 Comparison between the computed exit velocity vd
profile and the analytical Profile for ~
Opening Ratio =0.7


Full Text

PAGE 1

FINITE ELEMENT MODELING FOR ENERGY DISSIPATION IN A PRESSURE FLOW by Chwen-geng Guo B.S., National Cheng-Kung University, Taiwan, 1978 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of civil Engineering 1988

PAGE 2

This thesis for the Master of Science Degree by Chwen-geng Guo has been approved for the Department of civil Engineering by John R. Mays Ronald c. Scherer b. I PtPtft

PAGE 3

Guo, Chwen-geng (M.S. civil Engineering) Finite Element Modeling for Energy Dissipation in a Pressure Flow Thesis directed by Assistant Professor James c. Y. Guo. A finite element numerical model using quadratic isoparametric elements was developed to analyze a twodimensional pressure flow and to predict the flow energy losses associated with the friction and cross sectional contraction. The unsteady flow governing equations were converted into stream function and vorticity diffusion equations. Initial conditions were determined by potential flow and boundary conditions included a uniform incoming flow and non-slip condition on the solid wall. Flow cases studied included the opening ratios ranging from 0.3 to 1 and Reynolds numbers varying from 50 to 200. In general, this numerical approach provides stable computations for flow to reach its steady state when Reynolds number is less than 200 and opening ratio is larger than 0.3. Although for high Reynolds number flows the exit velocity profile may slightly deviate from the exact velocity profile, the predicted friction loss and orifice loss still give good agreements with the theoretical solutions and laboratory data.

PAGE 4

ACKNOWLEDGMENTS This project was supported by the Denver Center for the Performing Arts and funded by the Helen G. Bmfils Foundation. I would like to thank Dr. Ronald c. Scherer, senior scientist, Denver Center for the Performing Arts, and Dr. James c. Y. Guo, assistant professor, University of Colorado at Denver, for their advice and encouragement on this project. Also, I wish to express my thanks to thesis committee, Dr. David w. Hubly, and Dr. John R. Mays, for introducing me to the Fluid Mechanics and Finite Element Method. Finally, I would like to give my special thanks to my parents, my wife and children, for their understanding throughout the entire project.

PAGE 5

CONTENTS CHAPTER I. INTRODUCTION . . . . . . . . . . . . . . . . . 1 Objective . . . . . . . . . . . . . . . . . . 1 Scope of the study . 2 II. LITERATURE REVIEW 3 Energy Losses in An Incompressible Flow . . . 3 Numerical Modeling of An Incompressible Flow 12 III. GOVERNING EQUATIONS OF FLOW MOTIONS . 16 Flow Motions and Descriptions . . . . . . . . 16 Stream Function-Vorticity Equations of Flow Motion . . 19 Nondimensionalization of the Governing Equations . 23 Boundary conditions . . . . . . . . . . . . . .25 IV. NUMERICAL MODELING BY FINITE ELEMENT SCHEME .. 29 Basic Concepts ....... 29 Derivation of Finite Element Equations of Flow Governing Equations ... 35 Element Shape Function ... 40 Numerical Initial and Boundary Conditions 49 Computational Flow Chart . 53 V. CASE STUDIES AND RESULTS 57 Case One : Study of Flow Patterns 63

PAGE 6

vi Case Two: Prediction of Energy Losses 77 VI. CONCLUSIONS 99 ........ 102 APPENDICES A. THE DERIVATION OF ELEMENT EQUATION BY GAUSS-GREEN THEOREM 105 B. THE INTEGRATION OF SHAPE FUNCTIONS 109 C. DEFINITION OF SYMBOLS 115 D. THE FORTRAN SOURCE CODE 117

PAGE 7

TABLES Table 1. Computational Flow Chart . 55 2. Summary of the Flow Parameter Used in the case studies .......................... 58 3. Computed Velocity Distribution for Opening Ratio= 0.7 and Reynolds number= 50 73 4. Comparison of Peak Vorticity 76 5. Summary of Computational Time Intervals and the Required Times for Flow to Reach Steady State . 78 6. Computed Velocity Distribution in a Poiseuille Flow for Reynolds Number= 50 81 7. Computed Velocity Distribution in a Poiseuille Flow for Reynolds Number= 100 . 81 8. Comparison Between the Computed Exit Velocity Profile at X/B = 8.7 and The Analytical Profile for Poiseuille Flow . .. 84 9. Comparison of Friction Loss Between Computational and Analytical Predictions 10. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile 87 for the cases with a contraction .. 94 11. Distribution of Flow Pressure with a Contraction . 95 12. Comparison of Discharge Coefficients Between the Computed Results and Experimental Data 97

PAGE 8

FIGURES Figure 1. Illustration of Sudden Expansion . 7 2. Illustration of the Flow Domain with Orifice Plate . . . . . 10 3. Rotational and Irrotational Flows between Two Infinite Plates 18 4. Boundary Conditions of Poiseuille Flow . 27 5. Illustration of one-Dimensional Linear Finite Element ........ 30 6. The Example of Eight-Node Quadratic Quadrilateral Elements ..... 34 7. The Eight-Node Quadratic Quadrilateral Isoparametric Element with Natural Coordinates ... 43 B. Finite Element Net Work for the case with Opening Ratio of 1 . 59 9. Finite Element Net Work for the Case with Opening Ratio of 0.7 ....... 60 10. Finite Element Net Work for the Case with Opening Ratio of 0.5 ...... 61 11. Finite Element Net Work for the case with Opening Ratio of 0.3 .... 62 12. The steady State Distribution of Stream Function and Contours of Vorticity for Opening Ratio = 0.7 and Reynolds Number = 50 ....... 64 13. Illustration the Evolution of Flow Patterns Changing From Potential Flow to Its Steady State for Opening Ratio = 0.7 and Reynolds Number= 50 ....... 66

PAGE 9

14. The Variation of Velocity Profile along the Flow Direction for Opening Ratio = 0.7 ix and Reynolds number= 50 ........ 72 15. The Distribution of the Wall Vorticity for Opening Ratio= 0.7 and number= 50 . ." ... 75 16. The steady-state Velocity Profiles In A Poiseuille Flow for Reynolds Number= 50 79 17. The Steady-state Velocity Profiles In A Poiseuille Flow for Reynolds Number= 100 80 18. The Variation of Velocity Profile Along the Flow Direction for Opening Ratio = 1 and Reynolds number= 50 ... 82 19. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Poiseuille Flow with X/B = 8.7 83 20. Comparison of Friction Loss Between Computational and Analytical Predictions for Poiseuille Flow with X/B = 8.7 88 21. The Variation of Pressure Along the Flow Direction for Poiseuille Flow with X/B = 8.7 89 22. Comparison Between The Computed Exit Velocity and.The Analytical Profile for Open1ng Rat1o = 0.7 .. 91 23. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Opening Ratio = 0 5 9 2 24. Comparison Between The computed Exit Velocity Profile and The Analytical Profile for Opening Ratio= 0.3 93 '25. computed Pressure Drops .. 96 26. comparison of Discharge Coefficients Between Computed Results and Experimental Data 98 27. Coordinate Transformation 111

PAGE 10

CHAPTER I INTRODUCTION Objective Energy dissipation is an important subject to many professions. Examples include crude oil delivered in pipes, water supply systems, air flow through the larynx and blood flow in arteries. Numerous studies have been conducted to develop analytical and empirical equations to estimate energy losses associated with the shear friction and eddy motions. Recently, energy loss in the larynx flow has attracted much attention because of the lack understanding of flows near the glottis. The Denver Center for the Performing Arts, in association with the University of Iowa, has been engaged in a laboratory investigation of this phenomena. Since the energy dissipation in a flow is related to both flow and fluid properties, laboratory arrangements may not easily cover wide ranges of flow variables such as Reynolds number. As a result, numerical modeling is considered a possible helpful alternative. The purpose of this study is to develop a numerical model to analyze a two-dimensional pressure flow and to predict the flow energy losses associated with the wall

PAGE 11

friction and cross sectional-contraction. In this work, plate orifices are considered in order to establish the finite element technique before applying the technique to laryngeal flow problems. Scope of the study 2 After a review of numerical methods, the finite element method was selected to convert the nonlinear governing equations of unsteady flow into their numerical forms and to include the complicated boundary conditions in the numerical integrals. Solutions in terms of the stream function, vorticity, velocity and pressure drop at each time step were obtained from the numerical computations for unsteady flow. The steady state solutions were achieved by allowing the successive time-dependent solutions to converge. It is found that for low Reynolds number flows, comparisons between the predicted flow variables and analytical solutions are in good agreement. Although this developed numerical scheme is only for two-dimensional fl.ow, it is believed that the approach is also suitable for a similar three-dimensional flow problem.

PAGE 12

CHAPTER TWO LITERATURE REVIEW Energy Losses in an Incompressible Flow The main objective of this study is to develop a numerical model for the prediction of flow patterns and energy losses in a two dimensional pressure flow. In general, there are two kinds of energy losses involved in a pressure flow. They are friction loss and minor loss. Friction in a flow is caused by the non-uniform velocity distribution and fluid viscosity. Friction results in energy dissipation that converts mechanical energy into thermal energy. This energy loss to the environment is irreversible. Friction loss is determined by the length of the conduits, flow Reynolds number, and wall roughness. Using dimensional analysis, the friction loss equation in a pressure flow may be derived as follows (Daugherty, Franzini, and Finn.emore, 1985): ( Eq. 2.1) 2g where, hL= the head loss due to friction, L= the length of the conduit, Rh = hydraulic radius, V= cross sectional average flow velocity, g = acceleration of gravity, and Cf=

PAGE 13

4 Darcy friction factor. For a pipe flow, substituting Rh=D/4, in which D is the pipe diameter, into Eq. 2.1, the Darcy-Weisbach equation is obtained: ( Eq. 2. 2) D 2g where D = the pipe diameter, and f = 4Cf. To determine the pipe friction factor f the pipe flow experiments were conducted as early as the mid-1800s. For high Reynolds number flow, flow in rough pipe resulted in a higher energy loss than in smooth pipe. In 1933, J. Nikuradse experimented using artificially roughened pipes and obtained reliable information for friction factor (Gerhart and Gross, 1985). In 1944, the Moody Chart for Friction Factor f versus Reynolds number and relative roughness was published and, since then, it has been widely adopted by engineers. To expand the friction loss, Hagen and Poiseuille successfully derived analytical solutions in 1840 for laminar flows in a circular pipe and between two plates. They stated that for laminar flow in a circular pipe, the loss of head in friction is proportional to the first power of the cross sectional average velocity and is given by 1.1. L ( Eq. 2. 3) where 1.1. = absolute viscosity, and = the density of fluid. By equating Eq. 2.2 and 2.3, we have

PAGE 14

5 64 f = ----DV ( Eq. 2. 4) Since is flow Reynolds number, Eq. 2.4 can be further expressed as 64 f =---( Eq. 2. 5) Re .in which Re= flow Reynolds number. From Eq. 2.5, it is found that in laminar flow the friction loss is only a function of Reynolds number and independent of the roughness of the pipe wall. On the Moody chart, the friction factor f for laminar flow is a straight line. This fact is not only true for a laminar pipe flow but also for a laminar flow in a noncircular duct. Jones (1976) expanded Eq. 2.4 to duct flow with the circular diameter D replaced by the hydraulic radius, Dh' and conducted a series of experiments for various duct configurations. He suggested that for the laminar flow in a rectangular duct with width, a, and height, b, the friction factor may be described by k f =---( Eq. 2. 6) 64 where k = (Eq. 2.7) 2 + 11 _E._ (2 __ ) -3-24 a a and Dh = hydraulic diameter.

PAGE 15

For the laminar flow between two infinite plates, the value of variable a, in Eq. 2.7, is infinite. As a result, k is equal to 96 when Dh is equal to 4B (Blevins 1984, and Gerhart and Gross i985). Losses due to the change of local pipe configurations are termed minor losses. The flow in a pipe 6 system may need to pass through fittings, bends, or abrupt changes in the pipe configuration. These kinds of changes in local pipe geometries generate a great deal of disorderly motion and result in dissipation of flow momentum. Whenever the flow velocity is altered, eddy currents are induced and a loss of energy in addition to friction loss is created. Although minor losses are usually confined to a short length of flow path, the effects may not disappear for a considerable distance downstream. The minor loss may be expressed as h = K L 2g (Eq. 2.8) where K = loss coefficient, which must be determined experimentally for each situation. The only loss coefficient which can be calculated by simple analytical methods is that the axisymmetric sudden expansion as shown in Figure 1. This problem was first treated by Borda and Carnot. Using continuity, linear momentum, and energy equations and neglecting the shear stress along the wall, the head loss for sudden

PAGE 16

7 Figure 1 Illustration of sudden expansion '.

PAGE 17

can be expressed as 2g ) (Eq. 2.9) where A 1= the cross sectional area at position 1, A 2= the cross sectional area at position 2, and v1= the average velocity at position 2. Eq. 2.9 with Eq. 2.8 yields 8 K A1 2 =(1---) A2 (Eq. 2.10) Although we may be uneasy about the assumptions of the Borda-carnot approach and feel that Eq. 2.10 is too simple to evaluate the accurate energy loss for sudden expansion, the laboratory experiments showed that the difference between observed and predicted values of K in Eq. 2.10 were within a few percent (Vennard 1969). Based on this fact, it can be assumed that minor losses may be considered for practical purposes to occur in a much more localized flow region. As to the minor losses other than sudden expansion, they all have to be determined by experiments. In 1885, Weisbach followed the approach used by Borda-Carnot to observe the energy losses of abrupt contraction arid a fair estimate of loss coefficients was obtained. In 1952, Gibson's studies on diffusers of conical form grasped the trend of the change of loss coefficient, which are defined as a function of cone angle and area ratio. In

PAGE 18

Langhaar successfully calculated the loss of head at entrance by using boundary-layer theory. As to the losses in bends and elbows, there were studied by Beij and Hofmann in 1933. Hofmann gave loss coefficients for wide ranges of bend shape, roughness, and Reynolds number. 9 Nigro (1978) and Pao(1961) indicated that for a flow with Reynolds number below 104, the discharge coefficient becomes increasingly sensitive to Reynolds number, first rising and then falling as Reynolds number decreases. The same phenomenon as indicated by Negro and Pao was observed by Scherer (1981) in his experiments on laryngeal flow for a convergent-divergent glottis and opening ratios d/B less than 0.05, where d = the glottal diameter, and B = the subglottal (tracheal) lateral diameter. Nigro is The discharge coefficient for the pipe flow used by c -( [ 1 ] A 2 4 p ) (Eq. 2.11) where C = discharge coefficient, Ud= the mean flow velocity_ through the orifice throat, a = the cross sectional area of orifice, A = the cross sectional area of the pipe, p = the pressure drop between the flange taps as shown in Fig 2, and = fluid density. Considering a unit width of a two-dimensional flow, Eq. 2.11 can further be simplified to evaluate the flow velocity between two infinite plates.

PAGE 19

'lo .. ... ..... .. I I __J i'l-I -Flange tap ...... -., ... _.. ..-: .( 'a,. ...... ,..,. ...... .. ,; -, _,. r.,.r .. ,., ... ,.-_.,. ........... ., _., .... -, ,!,; plate . __ .. :----:---Centerli!le ______ ... __ 2B 2d j_ If :: :-: .. .... -:. .. """11., ... _., .. ' ,, ,. ,. ,, tC' .... ........ ,. --. ;. ;. Solid wall Figure 2 Illustration of the flow domain with orifice plate Outflow ---..:,. _\ .. 1 ) / .. / ......... o

PAGE 20

11 c ( [ 1 (_Q_)2 ] 2 .. p ) ( Eq. 2.12) B where d/B = the opening ratio between two plates, d = half of the opening width at the orifice, B = half of the width between the two plates. As indicated by Nigrq, for laminar flow through a thin square-edged orifice with Reynolds number between 100 and 200, the discharge coefficient for flange taps is approximately .70 to .ao. The values of loss coefficients obtained by the pioneer researchers mentioned above were continually modified by researchers. There are presently many handbooks and technical papers published with experimental data (Idel'chek 1960, and ASHRAE Handbook of Fundamentals). Unfortunately, owing to different experiments conducted under different experimental conditions, data values for minor losses may be different for the same flow configuration. Therefore, the most representative one is difficult to determine. Furthermore, experiments for the evaluation of minor losses are common for turbulent flows, whereas the empirical data for laminar flows are less common. These are the main reasons why this numerical study for estimating the energy losses for a laminar flow is needed. In this study, the laminar flow patterns and energy losses associated with the plate orifice as shown in Figure

PAGE 21

2, are numerically calculated. Before the selection of numerical methods, the following literature review was performed. Numerical Modeling of An Incompressible Flow There are three approaches by which we can solve fluid problems, namely 1. Experimental 2. Theoretical; and 3. Numerical 12 The experimental approach has the capability of producing the most realistic answers for many flow problems. However, its cost is the highest one among these three methods. In theoretical approaches, assumptions should be made in order to make the problem tractable. Therefore, the problems that can be solved by theoretical approaches are restricted to simple boundary conditions and usually limited to linear problems. However, the shortcomings of experimental and theoretical methods do not exist in numerical approach. In the numerical approach, a limited number of assumptions are made and a digital computer is used to solve the governing equations of a flow problem. There is no restriction to flow linearity and even-complicated geometry can be treated. Furthermore, the computational costs are always much less

PAGE 22

13 than that of laboratory equipment for experiments. It is why the numerical method has become more important when the computing facilities become more efficient and accessible. Most people attribute the first definite work of numerical approach to Richardson (1910) who introduced point iterative schemes to numerically solve Laplace's equation. His scheme used data available from the previous iteration to update each value of the unknown and was called the finite difference method (FDM). Richardson actually carried out calculations for the stress distributions in a masonry dam with great success. Up to the late 1930's, FDM was used for solving flow problems. Allen and Southwell (1955) applied relaxation scheme to solve the incompressible, viscous flow around a cylinder. Although they were not the first scientists to solve these kinds of problems numerically, their results did capture the flow phenomenon and had good accuracy. At the same time, John von Neumann developed his method for evaluating the stability of numerical methods for solving time-dependent problems. Douglas and Rachford (1950) developed an implicit method for solving parabolic and elliptic partial differential equations. From then on, the basic knowledge for solving fluid dynamics problems by FDM was complete. Before the finite element method (FEM) was fully developed, FDM was used to solve flow problems. The

PAGE 23

advantage of FEM is that its rate of convergence is much faster than that of FDM. 14 The finite element method was originally used by aircraft structural engineers in the 1950's to analyze a large system of structural elements in the aircraft. Usually, there are two ways to formulate element approximate equations: one is the variational method, which ha.s a close relationship to the classical variational concepts of the Rayleigh-Ritz method 1909); the other one is the weighted residual method, which is modeled after the well-known method of Galerkin (1915). In weighted residual method, the error in the approximate satisfaction of the governing equations is not set to zero, but instead its integral with respect to the selected 11weights11 is required to vanish. The application of FEM to nonstructural problems such as fluid flows was initiated by Zienkiewicz and Cheung (1965). Pioneering research has employed weighted residual methods to formulate approximate element equations. Oden (1972) was among the first to derive the basic theoretical analog for the Navier-Stokes equations. Baker (1971, 1973) published his results for the analysis of incompressible flow using the Galerkin method. Olson (1974) used a pseudovariational finite-element algorithm for solving the biharmonic equation for stream functions. Lynn(1974) utilized the least-squares methods to develop a finite-

PAGE 24

15 element algorithm for laminar boundary-layer flow. In 1978, Backer successfully extended his studies from 2-dimensional laminar flow to 3-dimensional turbulent flow. From that time on, FEM was widely adopted for solving fluid dynamics problems. In this study, the numerical model for solving laminar, incompressible flow will be derived using FEM because of its flexibility of element shape and efficient computation. Galerkin's method is used to weight residuals in the approximate governing equations of flow. The procedure of formulation will be discussed in Chapter Four.

PAGE 25

CHAPTER THREE GOVERNING EQUATIONS OF FLOW MOTIONS Flow Motions and Descriptions Using viscosity as a criterion, fluid may be classified as real fluid and ideal fluid. For a-real fluid, flow motion results in shear forces between fluid layers. But for an ideal flow, the shear force is assumed to be zero everywhere. Although this is an idealized situation, it does simplify the complexity of flow problems. The shear stress in a real flow is termed internal friction stress and results from the cohesion and momentum interchange between fluid layers (Daugherty, Franzini and Finnemore, 1985). The magnitude of shear stress is determined by fluid viscosity and the flow velocity gradient which is usually expressed by flow rotationality. For a viscous flow between two infinite flat plates, the between viscosity and shear stress may be expressed as du r = (Eq. 3.1) dy shear stress, = the absolute viscosity, and dujdy = the time rate of angular deformation or the velocity I gradient. From the above equation, it is noted that whenever

PAGE 26

17 and wherever a flow velocity gradient exits, shear stress occurs correspondingly. Based on laboratory observations, the fluid in direct contact with a solid boundary has the same velocity as the boundary itself. Therefore, when the boundary is fixed, the no-slip condition is applied, i.e., the velocity of the fluid on the solid surface must be zero. As a result, the farther the fluid from the solid wall, the faster the flow velocity can be. This fact explains the reason why the existence of the wall results in flow velocity gradient and shear stress. As shown in Figure 3(a), the cross sectional velocity profile between two plates is parabolic for a laminar flow. This is a typical motion of a low Reynolds number viscous flow. Considering element 1 in Figure 3(a), the uneven velocity gradients in the x-and y-directions cause a rotation. This kind of flow motion is called rotational flow (Janna 1987). It means that the shape and orientation of any fluid element change from time to time and from one position to another as well. When the flow is assumed to be ideal, as shown in Figure 3(b), the flow velocity profile between two plates at any cross section is a straight line. As a result, the shape of a fluid element remains identical. It means that the fluid element does not rotate. This is so called irrotational flow (Janna 1987). Applying the concept of control volume to a twodimensional, viscous flow moving on the xy-plane, the

PAGE 27

(a) Rotational flow I .(b) flow .-Figure 3 Rotational and irrotational flows between two infinite plates 18 0 .............. i i L. S.. :{ -..,,,'

PAGE 28

19 rotational velocity of a fluid element in the z axis may be expressed as 1 ov ou az = ( ---) (Eq. 3.2) 2 ox oy in which, e = rotational velocity in the z-direction, u z velocity component in the x-direction, and v = velocity component in the y-direction. The vorticity used to express the rotational characteristic of flow, is defined as two times the rotational velocity, namely ov ou = = 29 = z z ( Eq. 3. 3) ox oy in which z= vorticity in the z-direction. Although Eq's 3.2 and 3.3 are derived from mathematic concepts; they offer a feasible way to measure flow motion. The following section will further explain how to use vorticity to form the governing equations of the motion. Stream Function-Vorticity Eguations of Flow Motion The basic governing equations for a two-dimensional incompressible ,viscous flow can be described by the principles of continuity and momentum. Conservation of momentum states: ou ou ou op o 2 u o 2 u (-+ u + v -) = ---+ IJ. ( + ot ox oy ox ox2 oy2 )

PAGE 29

20 (Eq. 3.4) and, ou ov ov op o 2 v o 2 v (-+ u + v -) = --+ ( + ) ot ox oy oy ox2 oy2 (Eq. 3.5) The continuity Equation states: ou OV --+ = 0 ( Eq. 3. 6) ox oy where = the fluid density, and = the viscosity. There are three different formulations used to solve Eq's 3.4 through 3.6 for flow problems. They are: 1. The stream function formulation combines all governing equations together and expresses all flow variables by stream function. As a result, a fourth order partial differential equation of stream function needs to be solved. Olson (1975) first applied this approach using finite element method to solve the flow around a circular cylinder. The results, including the separation phenomenon behind the cylinder, were found to be accurate. 2. The velocity-pressure formulation directly solves momentum and continuity equations with the pressure and velocity as unknowns. Kawahara et. al. (1974) used this approach to solve several steady state flow problems with Reynolds number up to 1200. 3. The stream function-vorticity formulation uses

PAGE 30

stream function, and vorticity, e, to replace velocity and pressure. The advantage of this approach is to reduce unknowns from three variables, u, v, and p, to two 21 variables, and e. Cheng (1972) successfully applied this approach to solve flows in a constricted channel with Reynolds numbers less than 100. Using the formulation, Eq's 3.4 and 3.5 are reduced to a diffusion equation of vorticity and Eq. 3.6 becomes a poisson equation of the stream function. In this study, the combination of and e is used to describe flow. The following section further presents the detailed derivations. From the principle of conservation of mass, the velocity components u and v are related to the stream function by the following equations, u = (Eq. oy v = -----(Eq. 3.8) ox Substituting Eq's 3.7 and 3.8 into Eq. 3.3 yields 2 = e (Eq. 3.9) in which + For an irrotational flow, Eq. 3.9 is further simplified to (Eq. 3.10) Taking the derivative ojoy of Eq. 3.1 and the

PAGE 31

derivative -ojox of Eq. 3.2 and adding the resultant equations together, we obtain or < + 0 ) = + v ot ox oy 22 .,.2 = (Eq. 3.11) Dt Eq. 3.11 is a two-dimensional diffusion equation of vorticity. Eq's 3.9 and 3.11 represent two coupled second order equations governing incompressible viscous flows. The continuity equation is implicitly satisfied by the relationship between stream function and velocity. In order to predict energy loss, we need to calculate the pressure distribution through out the flow field. Taking the derivative of Eq. 3.4 with respect to x and the derivative of Eq. 3.5 with respect to y and then summing the resultant equations together, we have .,.2p = ou 2 OV 2 ou ov [ ( ) + ( ) + 2 ----] ox oy oy ox (Eq. 3.12) Substituting Eq. 3.6 into Eq. ;3.12 yields ou ov ou ov ----] (Eq. 3.13) OX oy oy ox Eq. 3.13 may further be expressed by stream function, namely

PAGE 32

23 2 + -( ) ] (Eq. 3.14) Eq's 3.9, 3.11, and 3.14 are governing equations used in this study. It is noted that Eq's 3.9 and 3.14 take the same form as the Poisson equation provided that the terms on the right-handed side of each equation can be as constants. Nondimensionalization of the Governing Equations In fluid mechanics, the solution of flow problems may become more generalized when they are expressed in nondimensional forms. Without nondimensionalization, solutions of a numerical model may only be suitable for the specific conditions studied. In the process of nondimensionalization, certain flow dominating parameters like Reynolds number may be found. To perform nondimensionalization, some reference flow,quantities have to be chosen as characteristics of the flow system (Daily and Harleman, 1966). For a flow, the inflow velocity, u, may be chosen as the characteristic velocity and a half of the distance between two infinite plates, B, may be chosen as the characteristic length. A set of dimensionless flow quantities may be defined as follows: X y B B

PAGE 33

24 0 u 0 v u = v = u u to tu 0 p = p = B u2 cp 0= cp 0 B = (Eq. 3.15) UB u Substitution of Eq. 3.15 into Eq. 3.9, the equation of stream function becomes 02( 0 ) 02( 0 ) cp cp + = 0 ( x0)2 0 ( y0)2 or 'f' 2 ( cpo ) = (Eq. 3.16) similarly, Eq 3.11 becomes + u + v = ( + ) o(yo) ua o(xo)2 o(yo)2 (Eq. 3.17) The term of in the above equation is the reciprocal of flow Reynolds number which delineates flow regimes, i.e., laminar or turbulent. UB Re = ---(Eq. 3.18) Reynolds number physically represents the ratio of inertia force of flow to viscous force. For a viscous flow between infinite plates, gravity does not affect the flow pattern.

PAGE 34

25 It is also obvious that capillarity and surface tension are of no practical importance. Hence the significant forces affecting flow motions are inertia force and viscous force. Substituting Eq. 3.18 into Eq. 3.16 yields 1 (Eq. 3.19) Re Based on the similarity theorem, Eq's 3.17 and 3.19 shall depict the same flow dynamic behavior as long as the boundary conditions and Reynolds number remain the same. The dimensionless form of the pressure distribution equation is derived as follows: = -2 (-----------------) or o(u0 ) o(v0 ) o(u0l o(v0 ) = -2 ( ) (Eq. 3.20) o(yo) o(x0 ) o(x0 ) o(yo) Comparing Eq. 3.20 with Eq.3.14, it can be seen that the density of fluid is eliminated by nondimensionalization. In the fol.lowing sections, the flow problem will be solved using the nondimensional equations, and the superscript O's are neglected. Boundary Conditions When solving a partial differential equation, the boundary conditions play a determining role. In the study

PAGE 35

26 of the flow between two infinite plates, as shown in Figure 2, flow is symmetric to its centerline, therefore, only a half of the flow field is used to define the computational domain. The inflow is assumed to have a uniform velocity profile with a constant horizontal velocity u. The boundary conditions associated with each governing equation of flow are shown in Figure 4 and summarized below : On the boundary AB : Flow enters with a uniform velocity profile, therefore, we have v = 0 u = 1 = uy = Y = 0 and the value of pressure is set to be a selected datum value such as 100. on the boundary BC Applying the non-slip condition to this boundary, we have = 0, and. v = u = 0 To estimate the vorticities generated from the wall, it is necessary to apply the Taylor expansion to Eq. 3.3. The detailed derivation will be discussed later. On the centerline AD : Based on the continuity and Eq. 3.15, the stream

PAGE 36

p = Selected datum value .-. u=l centerline D : Y/B v=O cp=l' =0 ---u=l I v=O cp=y =0 c ...... .-"-_._.; ,.::-__:-:. ::. 6> ,: :-.. "'!-',:; .. ;, ._.._ :,,.-,.'--' 'J, . :.;._..._,....,,. -.--..--;t, I I : / __ ___. . -:;:--;::..-;.. r. Solid wall u=v=O cp=O f=fw Figure 4 Boundary conditions of Poiseuille flow X/B N -.J

PAGE 37

function at the centerline must be = 1 Due to the occurrence of the maximum horizontal velocity, we can write 28 oujoy = o, (Eq. 3.21) Considering the symmetry of the flow about the centerline, we conclude that v = 0, and f = o, on this boundary. On the boundary DC This is the exit of the computational domain. Cheng (1972) suggested to introduce a velocity profile defined by a cubic equation of distance measured from the wall to calculate stream functions the line DC, to ensure that Poiseuille flow is generated. This suggestion forces the velocity profile on the line DC to agree with the analytic velocity profile, i.e., a parabolic curve. In this study, it is found that as long as the computational domain is long enough, the exit velocity profile agrees with the analytic solution quite well. Therefore, the condition that the vertical velocity vanishes on the boundary DC is adopted and used in the computation.

PAGE 38

CHAPTER FOUR NUMERICAL MODELING BY FINITE ELEMENT.SCHEME Basic Concepts The finite element method is a numerical scheme and computational procedure for obtaining approximate solutions to many engineering problems (Segerlind, 1984). In general, the application of the method takes the following five steps: 1. Discretization and Selection of Element Configuration This step involves dividing the entire computational domain into a suitable number of "small" bodies called finite elements. The shapes of the elements may be triangular or rectangular, which is determined by the user. 2. Determination of the Approximatiqn Equation For each element, it is required to have an approximation equation for solutions in terms of the unknown nodal values. For example, a one-dimensional linear element AB, as shown in Figure 5, has an approximation equation expressed by a linear equation. The coefficients a1 and a2 can be determined by using the nodal conditions, namely cp = 1 cp = 2 at at X= X 1 X= X 2 (Eq. 4. 2)

PAGE 39

cp = al. + a 2 x Tt2 tl 1 2 rxl L x2 Figure 5 Illustration of one-dimensional linear finite element 30 I X

PAGE 40

Substituting Eq. X -x lp = L 1 4.2 into Eq. x-x + L 2 4.1 gives { Eq. 4. 3) Eq. 4.3 is a typical form for approximation equations used in the finite element method. The nodal values are 31 multiplied by linear functions of x, which are called shape functions or interpolation functions. To be more Eq. 4.3 may be rewritten as cp = N. l. l. ( i = 1 to n ) {Eq. 4.4) .. where = the unknown nodal values, Ni = shape functions. In the above linear case, the shape functions become X -x 2 N =--1 L and x-x 1 N =--2 L (Eq. 4.5) The subscript i in Eq. 4.4 is an index that indicates how many nodes are involved in an element. For a linear triangular element, the value of n equals 3. As to a quadratic quadrilateral element, the value of n equals a. 3. Derivation of Element Equations There are several ways to derive the element equations in the finite element method such as the variational method, and the weighted residual method. In this study, Galerkin's method is used, which is a kind of weighted residual method. The weighted residual method is used to minimize the residual resulting from.an approximate solution that is substituted into the differential equations governing a

PAGE 41

32 field problem (Desai,l979). Usuc;llly, it requires that L I Wi(x)R(x) dx = o 0 ( Eq. 4. 6) where R(x) =the residual, Wi(x) =the ith nodal weighting function, and L = the length of integral. There are a number of schemes developed for the weighted residual method to determine the weighting functions, among which are collocation, subdomain, least squares and Galerkin's methods. When Galerkin's method is used, Eq. 4.6 is evaluated using the same functions for Wi(x) as were used in the approximate solution 1984). Thus, for any element, the weightihg functions for deriving element equation are exactly the same as its shape functions. When there are n nodes involved in an element, there are n shape fqnctions that can be used as weighting functions to form n equations as Eq. 4.6 for solving n unknown nodal values 'n For a two-dimensional field problem, Eq 4.6 may be rewritten in its matrix form I [N ]T R(x) dA = 0 A (i = 1 to n) (Eq. 4. 7) where [N]T = row vector containing the element shape functions. After integration, we have n simultaneous equations with n unknowns of ti. These simultaneous equations may be expressed as

PAGE 42

33 = {F} (Eq. 4.8) where [K] = stiffness matrix or element property matrix, = the unknown nodal values, and {F} = force matrix or vector of nodal forces! which may be a zero vector. 4. Assembly of Element Equations to Obtain Global Equations and Introduction of Boundary Conditions The final aim is to obtain the global equations that approximately define the behavior of the entire system. As shown in Figure 6, when the eight-node quadratic quadrilateral elements are used ( for instance, the element one has nodes 1, 10, 15, 16, 17, 11, 3, and 2), the element equation can be written as where the superscripts represent the element number, and The sequence of node numbers takes the counterclockwise rule. Since the different elements have different nodes, we have to, based on the node numbers, incorporate all element equations into the global After the boundary conditions are introduced into the computation, the modified global equation similar to Eq. 4.8 may further be obtained.

PAGE 43

4 3 2 1 (2) (6) 11 25 ... ( 1) (5) 10 15 24 Figure 6 The example of eight-node quadrilateral -elements l 34 ..

PAGE 44

5. Computation and Solution Finding of the Global Equation The global equations are a set of linear or nonlinear simultaneous algebraic equations. They can be solved using the well-known Gaussian elimination method to obtain the value at each node. 35 In the following sections, a numerical model of twodimensional incompressible flow will be derived according to the steps stated above. Derivation of Finite Element Equations of Flow Governing Eguations To formulate element equations, a general form of the dimensional field equation can be used (Segerlind 1984). 02'P 02'P D --+ D Gcp + Q = 0 X OX2 y oy2 (Eq. 4.11) Comparing Eq. 4.11 with Eq. 3.9, we conclude that for the calculation of stream function, Dx = Dy = 1, G = o, and Q = e According to the Galerkin's method, the residual integral equation for any element e is given by ( ) I T o2m o2m {R e } = -[N] (Dx + D Gcp + Q)dA ax oy (Eq. 4 .12) where [N] is the row vector containing the element shape functions and is regarded as the weighting function in the Galerkin's method. To integrate the second-derivative terms in Eq.

PAGE 45

36 4.12, the Green-Gauss theorem may be applied. After integration, Eq. 4.12 becomes where I = r [N]T( D ---cosa + D sine )dr x ox Y oy (Eq. 4.14) {k(e)} IA< ox o[N]T o[NJ o[N]T o[NJ = + Dy )dA OX ox oy oy + IAG[N]T[N] dA (Eq. 4.15) and (Eq. 4.16) As to two-dimensional time dependent problems, the general form of their governing equations may be expressed as o2n on Dx --=.-L + D --"' + Q = -"'-ox2 Y oy2 ot (Eq. 4.17) Eq. 4.17 can be rewritten as follows: + Q = ( Eq. 4.18) where = D --+ D --x ox2 Y oy2 The residual integral is
PAGE 46

37 similarly, after integration, we have = 0 ( Eq. 4. 20) where {I(e)}, {k(e)}, and {f(e)} are the same as expressed in Eq. 4.13, and {t(e)} =the nodal values of The transpose {i(e)}T of {t(e)} is {t(e)}T = [ at n at ( Eq. 4. 21) ] (Eq. 4.22) The detailed derivation of Eq's 4.13 and 4.20 are presented in Appendix A. Next we apply the general forms to the flow problem. The governing equations for a two-dimensional, incompressible flow as derived in chapter two are: ., 2 = e: ( Eq. 3 .16) 1 .,2 De: e: = (Eq. 3 .19) Re Dt 2 ., p = 2 [ ou OV au ov -----] (Eq. 3.20) ox oy oy ox in which Eq's 3.16 and 3.19 are the two coupled governing equations of flow motion. Eq. 3.19 may be further expressed by stream function as follows:

PAGE 47

38 1 +----(Eq. 4.23) Re ot oy ox ox oy Because both and e are unknowns, Eq. 4.23 is a nonlinear partial differential equation. To overcome the nonlinearity of Eq. 4.23, some researchers (Cheng, 1972) suggested that an unsteady flow problem can be regarded as linear in stream function and vorticity e at each time step. We may consider two solutions ( en ) at the nth time step and ( en+ 1 ) at a time increment t later. The steady-state solutions are achieved by allowing the successive timedependent solutions to converge. Therefore, the governing equations of flow motion become 2 =e n (Eq. 4. 24) and 1 ,2 oen OE'n E'n+1 = + Re ot oy ox OE'n ox oy (Eq. 4.25) It is noted that, at each time step, the and en in Eq. 4.25 are no longer unknowns. Thus Eq. 4.25 becomes solvable numerically. After the steady state is reached, the velocity components u and v at every node can be obtained by taking the first derivative of the stream function with respect to y and x .,respectiveiy. Then thevelocity g+adients in the right-hand side of Eq. can further be computed. As a result, Eq, 3.20 may be regarded as a Poisson equation.

PAGE 48

39 If the shape functions used for calculating vorticities are the same as those used for stream functions, applying Eq. 4.13 to Eq's 4.24 and 3.20, and Eq. 4.20 to Eq. 4.25, the element equation can be obtained as where {I(e)} = {k(e)} = {k(e)} p (Eq. 4.26) ( Eq. 4. 27) ( Eq. 4. 28) --case + sine ) dr ox oy ( Eq. 4. 29) oEn+l OEn+l --'------'-case + sine ) dr ox oy (Eq. 4.30) op op ---case + sine )dr ox + / oy o[NJ oy ( Eq 4. 31) )dA ( Eq. 4. 32) (Eq. 4.33) ( Eq. 4. 34) i .. ( ..

PAGE 49

40 [C(e)] = IA [N]T[N] dA (Eq. 4.35) {f(e)} = I ([N]T n) dA (Eq. 4.36) 'P .A {f(e)} I [N] T ( 6 'Pn+l On On )dA = oy ox ox oy A (Eq. 4.37) {f(e)} -2 I [N]T au ov au ov = (-------) dA p ox oy oy ox (Eq. 4.38) Eq's 4.26 through 4.38 are the element equations for solving a two-dimensional, incompressible flow problem. For an irrotational flow, the governing equation is y 2 'P = 0 (Eq. 3.10) Then, the element equation is reduced to ( Eq. 4. 26) where {I(e)} and [k(e)] can be referred to Eq1s 4.29 and 4.32 respectively. In the sections, the element shape functions will further be discussed and determined. Element Shape Function To apply the finite element method to a flow problem, we need to decide what.kind of elements and element equations should be used. The elements may be triangular, rectangular, and other shapes. The element equations may be

PAGE 50

!'i1 linear, quadratic, or cubic. There are two requirements that should be satisfied in the determination of element shape function. The first requirement is the continuity between adjacent elements. In Eq. 4.6, R(x) represents the governing equation. suppose that R(x) contains derivatives up to the r-th order. To avoid the integral becoming infinite, the compatibility requirement (Huebner, 1975) requires that its (r-1)-th order derivative must be continuous. This ensures that the piecewise continuity is maintained at element interfaces. In this study, the linear element shape function may be chosen. ;Because the governing equations of flow motion contain the second-derivative terms, a linear element with continuity of first order can sufficiently satisfy this requirement. The second requirement, criterion of completeness, results from the approximation theorem ( Zienkiewicz, 1977). finite element method has an implicit assumption that convergence occurs as the size of each element decreases. To ensure that this assumption is satisfied well, at least a constant value can be obtained in a local element when the size of any element tends to zero. Therefore, if a derivative of order r exists in the governing equation, it is necessary for the element shape function to be at least of the order r so that, in the limi.t, such a constant value can be acquired. Based on this criterion, a quadratic

PAGE 51

42 element shape function was chosen. It is noted that the criterion of completeness is not as necessary condition as the requirement of compatibility, for determining the order of element shape function. But, in this study, since a time-dependent equation of vorticity (Eq. 3.19) needs to be dealt with, it is better to use a quadratic element shape funqtion to ensure an efficient convergence. In addition, considering the application of this numerical model to a flow problem with curved boundary, the quadratic quadrilateral isoparametric elements were selected and Qsed in this study. Isoparametric elements are useful in modeling a field problem with the curved boundary and in grading a mesh from coarse to fine ( Cook, 1981). They make it possible to use nonrectangular quadrilateral elements for solving a field problem with the curved boundary. The interpolation equation for the eight-node quadratic quadrilateral isoparametric element shown in Figure 7 is = a1 + a2 a + a3 B + a4 aB + a5 a 2 + a6s2 + a 7 a 2 B + a8aA2 (Eq. 4. 39) with nodal conditions : a = -1 and B = -1 at node 1, a= 0 and B = -1 at node 2, a = 1 and B = -1 at node 3, a= 1 and B = 0 at node 4, a= 1 and B = 1 at node 5,

PAGE 52

a J Figure 7 The eight-node quadratic quadrilateral r isoparametric element with natural coordinates 43

PAGE 53

44 Q = 0 and B = 1 at node 6, Q = -1 and B = 1 at node 7, and Q = -1 and B = 0 at node a. (Eq. 4.40) in which a and B are natural coordinates. The element shape functions can be obtained by applying their basic characteristic to Eq. 4.39. The basic characteristic of the shape functions is that each shape function has a value of one at its own node and zero at the .other node. Applying this characteristic to Eq. 4.39 gives a standard form as follows: cp = N.t. 1 1 ( i = 1 to a ) where t. = the unknown nodal values, 1 functions, in which N1 = N2 = N3 = \(1+a)(1-B)(a-B-1) N4 = N5 = \{1+a)(1+B)(a+B-1) N6 = N 7 = \(1-a)(1+B)(-a+B-1) N 8 = (Eq. 4.4) and N. 1 = shape (Eq. 4.41) Applying Eq 4.41 to Eq's 4.29 through 4.38, the integration can be performed. There are two types of integrations: 1. The integrations occur on the boundary, such as

PAGE 54

45 Eq's 4.29, 4.30, and 4,31. This kind of derivative boundary condition will be discussed in section 4.4. 2. The integrations occur within the elements. They can further be subdivided into two different types: the integrations of the shape functions such as Eq's 4.35, and 4.36, and the containing partial derivative terms such as Eq's 4.32, 4.33, and 4.34, that involve the first derivatives of the shape functions, and Eq's 4.37 and 4.38, that include the first derivatives of stream functions, vorticities, and velocities. For the purpose of integration, the change of variables is required. For a two-dimensional integral, the change of variables equation is I I A 1 1 f{x,y)dA =I I f[x{a,B),y{a,B)] I det[J] I dadS -1 -1 (Eq. 4. 42 )" in which f(x,y) = the function expressed by global coordinates, f[x{a,B),y{a,B)] =the function expressed by natural coordinates, and [J] = Jacobian matrix, which is defined as [J] = o(x,y) o(a,B) = ox oa ox 68 {Eq. 4.43) and det[J] = the determinant of the Jacobian matrix. Applying Eq. 4.43 to Eq. 4.35 yields

PAGE 55

46 (Eq. 4.44) Similarly, Eq. 4.36 becomes 1 {f(:)} = JJ_1([N(a,B)]T[N(a,B)] {e }) ldet[JJI dadS ( Eq. 4. 45) Eq's 4.44 and 4.45 can be evaluated by Newton-Cotes or Gauss-Legendre quadrature method. For those integrals containing the partial derivatives of the shape functions, we need to change the local derivatives in terms of a and B to global derivatives. The relationships between local derivatives and global derivatives are oN. oNi 1 (a,B) (a,B) ox -1 oa = [J(a,B)] (Eq. 4.46) oN. oN. 1 (a,B) l (a,B) oy os where [J] =Jacobian matrix, and =inverse of [J]. However, in Eq. 4.46 an explicit form of [J]-1 cannot be obtained (Zienkiewicz,1977). Numerical schemes may be used to evaluate element integrals. In this study, GaussLegendre quadrature is adopted. The Gauss-Legendre quadrature selects the sampling points between the integral intervals to achieve the best accuracy. For double it states that

PAGE 56

47 R = J l Jl m -l 1f(a,B) dadB (Eq. 4. 47) where m is the number of the selected sampling points with respect to B, which can be determined by equating (2m-1) to the highest power of B existing in the function f(a,B). If the resulting m is not an integer, the next larger integer is used. n is the number of the selected sampling points with respect to a, which can be determined the same way as that used for determining m. wi is the weighting factor for the sampling points with respect to B, and Wj is the weighting factor for the sampling points with respect to a. The detailed derivation of Eq's 4.42 and 4.46 and the further explanations of Gauss-Legendre quadrature are presented in Appendix B. How many sampling points should be chosen when using Gauss-Legendre quadrature to maintain the ? In this study, referred to Eq. 4.44, the highest power of variables. appearing in [N]T[N] are a 4 and s4 Therefore, the value of 3 is used for both m and n. It means that the order of quadrature is 3 and there are 3x3 po1nts used for computation. Cook (1974) suggested for plane elements, either linear or quadratic, the 2x2 points are usually the best, except for a markedly nonrectangular or elongated

PAGE 57

48 quadratic element, where 3x3 points may give better accuracy. In this study, the numerical model developed may further be applied to the analysis of laryngeal flow. Considering the irregularity of larynx geometries, the nonrectangular and elongated elements were expected. Therefore, 3x3 points are considered the most reasonable choice. Applying Eq. 4.47 to Eq. 4.44, it gives [C(e)] 3 3 )]T[N(a = I: I: {W W [N(a ,B ,B )]} 1 det[JJ 1 i=l j=l i j j i j i (Eq. 4.48) The local coordinates of each sampling point and its corresponding weighting factor in the apove equation are presented in Appendix B. After the coordinates of each sampling point were determined, det[J] at each point can be obtained using Eq. B.3. Then Eq. 4.48 can further be evaluated numerically. The same procedure is applied to Eq. 4.45. As to Eq 4.32, which contains the first derivatives of the shape functions, with the aid of Eq's 4.46 and.4.47, it becomes {k (e)} 3 3 o[N(aj,Bi)J o[N(a.,B.)J = I: I: w. w. ( J ]. cp i=l j=l ]. J ox ox o[N(a.,B.)] o[N(aj,Bi)] + J ]. ) 1 det[JJ I -oy oy (Eq. 4.49) Similar approach to the forementioned can be applied to Eq's

PAGE 58

49 4.33 and 4.34. The term of in Eq 4.37, may be evaluated by applying Eq. 4.46 to the center-point of each element. For any element, equals a constant. By applying Gauss-Legendre quadrature, Eq. 4.37 becomes {f(e)} 3 3 T = :E :E (w.w. [N(aj,Bi)] } 1 det[JJ I Cc i=l j=l 1 J (Eq. 4.50) 0 0 where Cc = ( n n ) evaluated at oy OX ox oy the center of each element. Eq. 4.38 can be solved by the same way as used for ( ou ov ou ov Eq. 4.37. The term of ox oy ox ), 1n Eq. 4.37 are also evaluated at the center of each element. Numerical Initial and Boundary Conditions Since Eq. 3.19 is a time-dependent equation, the initial conditions have to be determined. In this study, the initial flow field was considered to be a potential flow (ideal flow) with a uniform inflow of velocity u at the entrance. As shown in Figure 4, the initial vorticity is zero everywhere in the computational domain. To determine the distribution of stream function, the governing equations for the irrotational flow were used. There are two kinds of boundary conditions involved in this study. They are, Dirichlet boundary conditions, ..

PAGE 59

50 that are the values of specified on the boundaries, and Neumann boundary conditions, that are the values of derivatives and specified on the boundaries. Recall that in the derivation of element equations, the boundary conditions have already included as a part of element equations (Eq's 4. 29, 4. 30, 4. 31) and take the following form ( D ---cose + D Y sine ) x ox oy (Eq. 4.51) where e = the angle to the outward normal to the boundaries of each element. If Dx = Dy = D, then Eq. 4.51 becomes D = 0 on where is the derivative normal to the boundary. To evaluate Eq. 4.29, in which theNeumann boundary conditions of the stream functions were involved, the derivative boundary conditions of stream functions in Figure 4 are listed below: along line AB, = o, and e = 1ao0 along line BC, = o, along line AD, = o, and e = 90, along line DC, = o, and e = oo, Substituting Eq. 4.52 into Eq. 4.29 yields {I (e) } = 0 (Eq. 4. 52) (Eq. 4. 53) Eq. 4.53 implies that the Neumann boundary conditions for

PAGE 60

51 stream functions are self-consistent in the element equations. If boundaries AB and DC in Figure 4 are vertical straight lines, the contribution of the derivative boundary conditions of the stream functions to the element equations automatically vanishes. Therefore, for the stream function, only Dirichlet boundary conditions need to be considered. As a result, Eq. 4.29 can be removed from the element equations. The Dirichlet boundary conditions of the stream functions as shown inFigure 4 are listed below: Along line BC, = o, and along line AD, = l How to impose Dirichlet boundary conditions by modifying the global equations can be found in Reference 7, 9, and 29. To consider the numerical boundary conditions of vorticity, a linear combination of E and its normal derivative is required to be specified throughout the entire computational domain (Backer, 1983). On the fixed wall, BC, the value of vorticity is a priori unknown. To evaluate.the wall vorticity, the no-slip boundary condition alone is insufficient. From Eq. 4.24, at a point (x0,y0 ) on the wall, the vorticity may be by (Eq. 4. 54) where B is the coordinate normal to the wall. To evaluate

PAGE 61

52 joB2 we can apply Taylor series expansion to a point (x1,y1 ) which is located a small distance from the wall along the B direction. ( Eq. 4. 55) since the no-slip condition indicates that OB ( xo Yo ) = o we have 2 en(xo,Yo) = 8 2 (Eq. 4 . 56) Wall vorticity is then expressed in terms of the stream functions on the wall and the node at a small distance from the wall. Because B is the coordinate normal to the wall, the construction of the grids for computing wall vorticity do need special attention. J As to the derivative boundary condition of vorticity, the values of oen+l/ox along boundaries AB, BC, CD and oen+l/oy along boundary BC are undeterminable. Eq. 4.30 may be considered as a negligible term which is further ignored from the global equations . In this study, it is shown that neglecting the derivative boundary conditions of vorticity can still give acceptable results which agree to analytic solution well. The same situation may be applied

PAGE 62

53 to the derivative boundary conditions of pressure. By adding case times (Eq.3.4) to sine times (Eq.3.5) and using the appropriate boundary conditions, we have ( Eq. 4. 57) Substituting Eq. 4.57 into Eq. 4.31, it is found that the integrations of the second derivative of velocity along the element boundaries are difficult to evaluate. Like Eq. 4.30, Eq. 4.31 is considered as a numerically negligible term and ignored from the global equation. Computational Flow Chart To develop a digital computer program for analyzing a two-dimensional, incompressible, viscous flow between two infinite plates, the initial conditions are determined using the governing equations for ideal flow and then Eq's 4.26 and 4.27 are used to calculate stream functions and vorticities for successive time steps. Neglecting the derivative boundary conditions, Eq's 4.2 and 4.27 are reduced to [k(e)] cp {cp(e)} = n+1 {f(e)} cp ( Eq. 4. 58) [k(e)]{f:(e)} + [c(e)]{;(e)} = {f(e)} n+1 n (Eq. 4.59). in which d fn+1 n {f:n } = { dt } = { .&t } (Eq. 4.60)

PAGE 63

54 Substituting Eq. 4.60 into Eq. 4.59, the recurrence relation is (Eq. 4.61) At each time step, Eq. 4.58 is first used to calculate Substituting new values of to Eq. 4.58, the corresponding wall vorticities can be obtained. Using Eq. 4.61, the vorticity can be updated for the next computational time step. The criterion for determining the steady state was defined to be whem the largest value of OE/Ot in the entire computational domain was less than 0.04 (Backer, 1983). After the steady state is reached, using Eq. 4.28, the corresponding pressure distribution is further computed. Comparing pressure drops between downstream and upstream flows, we can determine the energy dissipation. A flow chart illustrated the computational procedures is presented in Table 1. The source code and FORTRAN program can be found in Appendix D.

PAGE 64

FLOW CHART A I Define Flow Field 1 Discretize and Select the Element configuration Specify the Approximation Equation The Isoparametric Quadratic Quadrilateral Elements are Used Set Up Initial condition .,2 cp= 0 The Boundary Conditions for Initial State Are Shown in Figure 4. (Continued) Table 1 computational Flow Chart. 55

PAGE 65

56 Flow Chart B t = 0 .,.2 cp = 0 and = 0 l J t = t + &t J 'l l Determination of wall vorticity .. 2 ... f = -(cp cp ] , .... w 2 w w+l B l 1 .,.2 D Re = Dt l .... 2 cp = l Check steadiness no 0/0t < 0.04 I yes Calculation of Pressure "2P 2 [ ou ov ou OV ] = CSX6Y-6Y6X I End I

PAGE 66

CHAPTER FIVE CASE STUDIES AND RESULTS There are four different contraction geometries used in this study. They can be distinguished by the opening ratio, r, which is defined as r = d B (Eq. 5.1) where d = a.half width of the opening, and B =a half width of the distance between two plates. When r = 1, flow does not have any contraction. Results of these cases present predictions of friction losses only. When r < 1, flow is subject to minor losses due to contraction. The opening ratios used in this study were 0.3; 0.5, 0.7, and 1.0. Table 2 summarizes the total number of nodal elements, and Reynolds numbers used in the studies of these four cases. The corresponding element network for each case is shown in Figures 8 through 11. Since the abrupt changes of stream function and vorticity are expected to occur near the orifice, those cases with a smaller opening should use a larger number of smalier elements around the orifice. Thus, the number of nodal points and elements around the orifice used in the case with an opening ratio of 0.3 has as nearly 1.5 times nodal points and elements as used in the case with an opening ratio of 0.7. The total length of the computational domain also

PAGE 67

Table 2 Summary of the Flow Parameters used in the case Studies. Number of Number Opening Reynolds nodal of Length Separation ratio number points element L/B eddy ==--=================--=====--==================== 1.0 50 413 120 8.7 No 1.0. 75 413 120 8.7 No 1.0 100 413 120 8.7 No 1.0 150 413 120 8.7 No 1.0 200 413 120 8.7 No 0.7 50 235 64 6.5 Yes o. 7 100 235 64 6.5 Yes .7 200 235 64 6.5 Yes 0.5 50 309 86 8.7 Yes 0.5 100. 309 86 8.7 Yes 0.5 20" 0 309 86 8.7 Yes 0.3 5 0 305 84 9.3 Yes 0.3 100 305 84 9.3 Yes 0.3 150 305 84 9 3 Yes ====--== ... U1 Q)

PAGE 68

FINITE (0.0, 1.0) y/9 (0.0 ,0.0)" ELEMENT GRID FOR
PAGE 69

FINITE ELEMENT GRID NET'w'DRK FOR R=0.7 -. y/B I I II kk Orifice
PAGE 70

._ .. .l. :: .. . FINITE ELEMENT GRID NETV/DRK FOR R=0.5 (CENTERLINE OF (0.0, 1.0) y/B 1111 DNJ( (0.0 ,0.0) -ORIFICE .. (SOLID VIALL) x/B. Figure 10 Finite element grid net work for the case with opening ratio of o.s \ \ (8.7, 1.0) (8.7, 0.0)" --_j ..r.;,'-'! ) l 0'1 .... ..... \:; ,., ',. ...., ___ ,..

PAGE 71

FINITE ELEMENT GRID NET\JORK FOR 'R=0.3 co.o, 1.0) CCENTERLINE OF FLO"'w'). .y/B II II ll if (0.0 ,0.0) .Orifice (SOLID 'w'ALL) Figure 11 Finite element grid net work for the case with opening_ratio of o.J x/B. (9.3, 1.0) <9.3, o.o> 0\ N /-,-. i ...

PAGE 72

63 needs special attention. After the flow passing the orifice, it takes a longer distance for the flow with more contraction to return to its parabolic velocity profile than that with a smaller contraction. When the computational domain downstream of the orifice is not long enough, the numerical procedure may become unstable and results are not reliable. In this study, it was found that the range of computational domain, in terms of L/B varied from 6.5 to 9.3. Predictions using these computational domains give good agreements to the laboratory and analytical results. Predictions and comparisons for each case study are presented as follow: Case One : Study of Flow Patterns As may be expected, eddies are generated behind the orifice plate. Table 2 indicates that as long as the opening ratio is less than unity, eddies always occur. For high Reynolds number flows, eddies immediately downstream of the orifice plate may further shed out. Figure 12 presents steady state flow with an opening ratio of 0.7 and Reynolds number of 50. The computational time interval, At, was 0.01. It took 640 time steps, i.e. t 0 = 6.4, for flow to reach the steady state. The upper half of Figure 12 shows the distribution of stream function and the lower half presents the contours of vorticity. It is noted that a reverse flow occurs immediately downstream of the

PAGE 73

L------:----====()-------------J --cp=. -20--------------_J o.S Figure 12 'The steady state distribution of stream function and contours of vorticity for Opening Ratio = 0.7 and Reynolds Number = 50 .... .. / 0\

PAGE 74

orifice plate and the highest vorticity is produced near the orifice. On this figure, the standing eddies are clearly shown. The general flow patterns predicted in this study are very similar to that reported by White (1974), Chung (1978), and, Kwon and Pletcher (1981), although the numerical approaches were different. 65 Evolutions of flow patterns with respect to time changing from potential to steady state flow are shown in Figure 13. It can be seen that the eddies are growing in size as time goes on. Using the steady state criterion mentioned in Chapter four, flow becomes steady when t 0 = 6.4. However, it is found that when t 0 = 3.2, the variations of streamline patterns have already been numerically negligible. This may imply that the criterion for steady state suggested by Baker (1983) seems to be conservative. Due to the inadequate accuracy of software for plotting streamlines and the coarseness of elements used, some flow patterns in Figure 13 are not shown to be separating from the edge of the orifice. Figure 14 displays the plot of velocity profiles along the flow direction at X/B = o., 1.1, 2.2, 2.5, 2.9, 4.7, and 6.5. Near the entrance, the velocity profile shows that the flow has not been fully developed. Behind the. orifice, negative velocity occurs. This means that eddy motions in the area of wake produce a reverse flow. Table 3

PAGE 75

.:..: .. (. .. t o.o 1 Y/B -:::::::: 61 I I I I I I I I I I I I I 0 1 2 3 Lt. 5 X/B i 0 ... t 0.2 1 Y/B i :::::=::::: 6$. I I I I I ' i fJ 1 2 3 Lt 5 X/B to 0.4 1 } Y/B t I I I I I I I I I I 0 1 't 5 X/B to 0.6 .. 1 Y/B (continued> Figure 13 Illustration the evolution of flow patterns changing from potential flow to its steady state for Opening Ratio 0.7 and Reynolds Number= 50 ) --66 I I 6 I G ' ... I 6

PAGE 76

67 ______ _.::,_: J ,0 . . : 1. 2 . : : _. 3 .. 't . 5 X/B 6 .. .. --.--.-. -. I I : 'l-c=================================::::::' to 1 o -.. .. ,. Y/B 2 X/8 ::-.. 0 t l. 2 l 1 3 X/B .. tL-----:-___ ---::===:---------.0 I I I I I I I I 1 .. 1 II I I I 1 . 2 I . 3 't 5 X/B ... ---.

PAGE 77

I. 0 t 1.6 . . ... I( 1 2 3 5 I I I G t0 1.a l'-c==========:======================== y/B 0 t 2.0 Y/B 1 2 3 Y/B 2 3" ___ 5 ._X/B G ... r-

PAGE 78

69 0 t 2.4 ..:;;::'. y/B f:---:::pr l \-2._ -.... .. ,,, 0 1 2 3 tt 5 X/B G 0 t II 2.6 I' I I -1., I e 1 2 3 tt 5 x;a G t0 11 2.a r?... 1 . Y/B C--.--:------:---======-------. .!: .. 0 1 2 3 lt. 5 X/B I I I I : . G 0 t II 3o0 .. 1-c==================================== Y/B I 1 = I I I I I I I .. ---... -.. ---. .:. -----. -----------. 5 X/B. G
PAGE 79

70 0 t 3.2 Y/B :r::: -:-I I I I I If I I I I I 0 1 2 3 ... '+ to 3 .4 I I I I I I I 5 X;B 6 . 1--c:=============================== Y/B 6 t::----?f&:_ -, 1''' 'I 0 1 2 3 '+ 0 t 3.6 I I I I I I 5 X/B I G Y/B .-r 1 2 3 Lf 5 6 0 .t 3.8 1-c=================================== Y/B 6 I 0 a 1 2 3 ( '+ I I :: I I I I I I i I 5 X/B 6

PAGE 80

D 1 2 3 o. t 4.8 I I I .lt I ' 5 xis I I 6 71' y/B. '. 1 2 3 0 : t .. 5.2 : I j I 4 I I I I I 5 6 .X/B Y/B i---------:-: 0 I 0 1 2 3 4 0 t .4 I I I I 5 X/B I I I G Y/B I I I I I I i I I I I I I 0 1 2 3 4 5 G X/B ."' 0 ,' j I : ..

PAGE 81

Y/B _f//// 4u=l 2 Figure 14 The variation of velocity profile along the flow direction for Opening Ratio = 0.7 and Reynolds number = 50 X/B ....,J N

PAGE 82

73 Table 3 Computed Velocity Distribution For opening ratio= 0.7, and Reynolds number= 50 ============================================================ o.oo X/B 1.10 2.50 2.90 4.70 6.50 2.20 ============================================================ Distance Velocity Distribution At orifice yjB ujUo ujUo ============================================================ 0.00 1.00 0.00 0.00 0.00 0.00 o.oo 0.05 1.00 0.18 -0.09 -o.o8 0.08 0.14 0.10 1.00 0.36 -0.17 -0.13 0.15 0.27 0.20 1.00 0.65 0.03 0.05 0.36 0.54 0.30 1.00 0.92 0.33 0.28 0.58 0.77 o.oo 0.45 1.00 1.11 1.03 0.92 0.96 1.06 1.27 0.60 1.00 1.27 1.67 1.61 1.37 1.29 1.50 0.80 1.00 1.30 1.74 1.80 1.58 1.45 1.57 1.00 1.00 1.33 1.80 1.96 1.77 1.57 1.64 ============================================================ ----------------------------------------------------------------------------Exit Profile Analytical Profile ====================================== Distance Y/B Velocity Profile ujUc U/UC ====================================== o.oo 0.05 0.10 0.20 0.30 0.45 0.60 0.80 1.00 o.oo 0.09 0.18 0.34 0.49 0.68 0.82 0.92 1.00 0.00 0.10 0.19 0.36 0.51 0.70 0.84 0.96 1.00 Note: X= horizontal distance. y= distance from solid wall. B= half width between two plates. Uo= entrance velocity. Uc= centerline velocity

PAGE 83

74 presents velocity profiles at different locations and a comparison between the exit velocity profile and the analytic parabolic velocity profile. It is noticed that although the existence of the orifice does disturb velocity profile, flow still returns to Poiseuille profile far downstream. In this study, the velocity profile for Poiseuille flow can be expressed as (Gerhart and Gross, 1895) U/Uc = 1 -( 1-y2 ) where u= velocity in x-direction, Uc= u at centerline, y= distance from the solid wall. The wall vorticity is also worthy of studying because the wall shear stress is directly proportional to its value. Figure 15 depicts the distribution of the wall vorticity for the flow with an opening ratio of 0.7 and Reynolds number of 50. It shows that the peak wall vorticity occurs at the contraction section and the wall shear stress changes its direction corresponding to the flow direction near the wall. Table 4 presents the values of peak vorticity for various cases under different Reynolds numbers. It is noted that for the with an opening ratio = 1, the peak vorticity, highest shear stress, occurs at the entrance because the flow experiences the steepest velocity gradient at the beginning of the entrance length and its value is proportional to the Reynolds numbers. On the contrary, for those cases with an opening ratio less

PAGE 84

22 20 18 16 14 .a 12 ] 1:! 10 0 8 6 4 2 0 -2 0 2 4 6 X/B Figure 15. The distribution of the wall vorticity for opening Ratio= 0.7 and Reynolds number 50 Ul

PAGE 85

Table 4 Comparison of Peak Vorticity. ============================================================= The Values Opening. Reynolds of Peak occuring ratio number Vorticity Position ============================================================= 1.0 50 16.62 Near Entrance 1.0 75 17.04 Near Entrance 1.0 100 17.36 Near Entrance 1.0 150 17.84 Near Entrance 1.0 200 18.19 Near Entrance 0.7 50 20.78 Near Orifice 0.7 100 23.49 Near Orifice 0.7 200 26.43 Near Orifice 0.5 50 34.15 Near Orifice 0.5 100 .38. 98 Near Orifice 0.5 200 42.40 Near Orifice 0.3 50 60.89 Near Orifice 0.3 100 66.80 --. Near_Orifice ---o-:3-150 73.77 Near Orifice ============================================================= -.J O'o

PAGE 86

77 than one, the highest shear stress occurs at the orifice and their values are a function of both Reynolds number and the opening ratio. With the same Reynolds number, the value of peak vorticity with an opening ratio of 0.3 is as nearly three times as that with an opening ratio of 0.7. Computational time intervals and the required times for flow to reach the steady state are presented in Table 5. It is noted that the higher the Reynolds number is, the more time steps are required. Case Two : Prediction of Energy Losses Energy dissipation in a pressure flow includes friction loss and minor loss. In this portion, flows are studied with and without an orifice. When the opening ratio is unity, the energy loss is essentially a friction Figures 16 and 17 present the steady-state velocity profiles at various locations with Reynolds number equal to 50, and 100 respectively. The corresponding data can.be found in Tables 6 and 7. Flow starts with a uniform velocity profile at the entrance anq gradually changes to a fully developed flow through the entrance length. Figure 18 depicts the variations of velocity profiles along the flow direction with Reynolds number of 50. Figure 19 and Table 8 show comparison between the exit velocity profiles under Reynolds. numbers equal to 50, 75, 100, 150, and 200 with the analytic velocity profile of a fully developed flow. It is noted that the exit velocity profiles with

PAGE 87

Table 5 summary of Computational Time Intervals and the Required Times for Flow to Reach steady State. ========================================================== Time Interval Opening. Reynolds Used for ratio number Computation Time Steps for Reaching Steady State ========================================================== 1.0 50 0.0100 960 1.0 75 0.0100 1200 1.0 100 0.0100 1200 1.0 150 0.0100 1900 1.0 200 0.0050 2400 0.7 50 0.0200 320 0.7 100 0.0100 780 0.7 200 0.0050 2015 0.5 50 0.0100 633 0.5 100 0.0100 754 0.5 200 0.0050 1992 0.3 50 0.0100 567 0.3 100 0.0100 7 .37 0.3 150 0.0025 3376 '.J (%)

PAGE 88

0 :J :J ,... u 0 cu > 1.5 1.4 -, _m .:t 2.0 0 0 v ).9 ----+ 1.3 1.2 1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 A 6.3 X s. 7 0.2 0.4 Distance y /B 0.6 0.8 Figure 16 The steady-state velocity profiles ln a Poiseuilie flow for Reynolds Number 50 1.0 -.J \D

PAGE 89

2 1.9 1.8 1.7 ..6 1.5 1.4 1.3 0 1.2 :::> 1.1 :::J p 1 'ij 0.9 0 Ill 0.8 > 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.00 X/B o o + 0 3.9 .1::.. 6. 3 v a. 1 I I--0.20 0.40 Distance y/B 0.60 0.80 Figure 17 The steady-state velocity profiles in a Poiseuille flow for Reynolds Number.= 100 .. 1.00 CXI 0

PAGE 90

81 Table 6 Computed Velocity Distribution in a Poiseuille Flow Reynolds Number= so. ============================================================ X/B o.oo 2.00 3.90 6.30 8.70 Exit Analytical Profile Profile ============================================================ Distance Velocity distribution Velocity Profile yjB U/UO u;uc U/UC ==================.========================================== o.oo 1.00 o.oo o.oo o.oo o.oo 0.00 0.00 0.05 1.00 0.18 0.16 0.15 0.15 0.15 0.98 0.10 1.00 0.35 0.32 0.30 0.30 0.30 0.19 0.15 1.00 0.49 0.45 0.43 0.42 0.42 0.28 0.20 1.00 0.64 0.59 0.57 0.56 0.56 0.36 0.40 1.00 1.07 1.01 0.99 0.98 0.98 0.64 0.60-1.00 1.27 1.27 1.27 1.27 1.27 0.84 0.80 1.00 1.33 1.39 1.43 1.44 1.44"0.96 1.00 1.00 1. 33 1.43 1.48 1.49 1.49 1.00 ============================================================ Table 7 Computed Velocity Distribution in a Poiseuille Flow Reynolds Number= 100 ============================================================ X/B Exit Analytical 0.00 2.00 3.90 6.30 8.70 Profile Profile ============================================================ Distance Velocity Distribution Velocity Profile Y/B U/UO U/UC ujUc ============================================================ 0.00 1.00 0.00 o.oo 0.00 0.00 o.oo o.oo 0.05 1.00 0. 21-0.18 0.17 0.16 0.11 0.01 0.10 1.00 0.40 0.36 0.33 0.32 0.23 0.19 0.15 1.00 0.56 0.50 0.47 0.46 0.32 0.28 0.20 1.00 0.74 0.66 0.61 0.60 0.42 0.36 0.40 1.00 1.13 1.08 1.03 1.02 0.72 0.64 0.60 1.00 1.24 1.27 1.27 1.27 0.89 0.84 0.80 1.00 1. 24 1.32 1.37 1.39 0.98 0.96 1.00 1.00 1.24 1.32 1.39 1.42 1.00 1.00 ============================================================ Note: Uo = entrance veloctiy. Uc= centerline velocity. X horizontal distance. y = distance from the plate. B = half width between two piates.

PAGE 91

_y)s 1 g u=1 ... 0 .. I I I I I I _1 I w f "1'' h' U.I.IU.I.Ii1.111.1/.l.l 'I 11//// .1/r.I.IIN F /N F/1'11" f.l /.'/.II' .I /NF'f/'-I'H".IUF 2 4-6 Figure 18 The variation of velocity profile along the flow direction for Opening Ratio = 1 and Reynolds number = 50 I I J I r r .. ,, I """"" 8 .. Q) liJ

PAGE 92

2 1.9 1.8 1.7 1.6 1.5 1.4 1.3 u 1.2 :::> 1.1 ::J .0 1 u 0.9 0 Ill 0.8 > 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 .. 0.00 Re 0 50 + 75 <> 100 )(. 150 /:). 200 ."'\1 Analytical 0.20 0.40 0.60. 0.80 Distance y/B Figure 19 -Comparison between the computed exit velocity profile and the analytical profile for Poiseuille flow with X/B = 8.7 .(--_ 1.00 co .w :...

PAGE 93

Table 8 Comparison between the Computed Exit Velocity Profile at X/B = 8.7 and the Analytical Velocity Profile for Poiseuille Flow ============================================================ Reynolds Number 50 75 100 150 200 Analytical Velocity Profile ============================================================ Distance Exit Velocity at X/B=8.7 Y/B u;uc U/UC ============================================================ o.oo 0.00 o.oo 0.00 0.00 o.oo o.oo 0.05 0.10 0.11 0.11 0.13 0.14 0.01 0.10 0.20 0.21 0.23 0.24 0.27 0.19 0.15 0.28 0.29 0.32 0.35 0.38 0.28 0.20 0.38 0.38 0.42 0.45 0.51 0.36 0.40 0.66 0.67 0.72 0.76 0.82 0.64 0.60 0.85 0.86 0.89 0.92 0.96 0.84 0.80 0.97 0.96 0.98 0.98 0.99 0.96 1.00 1.00 1.00 1.00 1.00 1.00 1.00 ============================================================ Note: Uo = entrance veloctiy. Uc= centerline velocity. X = horizontal distance. y = distance from the plate. B = half width between two plates. Q)

PAGE 94

85 Reynolds number less than 100 give a better agreement to the exact solution. The discrepancy between the predicted exit velocity profiles and analytic solution for high Reynolds number flows may be due to the length of computational flow domain. When it is hot long enough, flow may not become fully developed at the exit. Theoretically, an infinite distance is required for a flow with uniform entrance velocity to become fully developed, but it has been proved both by theory and observation that for a pipe flow, the maximum velocity at the centerline of the pipe may reach 99 percent of its ultimate value in the distance determined as follows: L = 0.058 RePo (Eq. 5.2) where Rep is the Reynolds number for pipe flow and defined as UD IJ. in which U =average flow velocity, and D =the diameterof the pipe (Daugherty, Franzini and Finnemore). Applying Eq 2.7 to Eq. 5.2 and using Reynolds numper determined by hydraulic radius we have L = .348 Re B (Eq. 5.3) For Re = 100, the required length to reach 99% of fully developed flow is 35 B. This is impractical to adopt this computational domain because of the elongation of

PAGE 95

computational time. After many tests, the length of 8.7 B is used in this study based on the judgement of the agreement between the predicted exit velocity profile and the exact velocity profile for flows with Reynolds number less than 100. 86 As to the prediction of friction loss, computational domain of L/B=8.7 can still give fair accurate predictions even for the flows with higher Reynolds number. Numerical predictions of friction loss expressed by flow pressure drop are presented in Table 9 and Figure 20. The analytical pressure drop in Table 9 was computed by Eq 5.3 that can be derived by applying Eq. 2.7 to Eq. 2.1 and using hydraulic radius for Reynolds number. 3 L 4 p = ( Eq. 5. 3) Re B where L/B = 8.7 in this study. It can be seen that the computed friction loss is in a good agreement to the analytical solutions. Figure 21 presents the variations of flow pressure along the flow direction with various Reynolds numbers. It shows that independent of Reynolds number, kinetic effects dominate over the first 30 percent of the entrance length. Thereafter, viscous effects gradually increase in importance and become the sole influence at the end of the duct, where the pressure gradient is a constant. To predict the energy losses due to the contraction,

PAGE 96

Table 9 Comparison of Friction Loss between Computational and Analytical Predictions ================================================ Reynold Number Computed Analytical Pressure Drop Pressure Drop 50 75 100 150 200 0.50 0.36 0.25 0.18 0.14 0.52 0.35 0.26 0.17 0.13 ================================================ Q) '-.J

PAGE 97

0.55 0.5 0.45 0.4 Ill '1 _3 0.35 c 0 :;::; 0 0.3 .:: LL. 0.25 0.2 0.15 0.1 50 70 90 110 130 150 Reynolds Number [] Computed + Analytical .Figure 20 Comparison of friction loss between comput"ational and analytical predictions for Poiseuille flow with X/B = s.7 170 190 (X) (X)

PAGE 98

p uo """2 QJ Ill Ill QJ 0.8 ..... ----' ------' ,, ------------------., ,____ ----,_. ...... .......... ----------------------------!: .. 0. 6 0 2 '+' 6 Distance X/B Figure 21. The variation of pressure along the flow directiort for Poiseuille flow with X/B = 8 Re 200 150 100 75 50 .> (X) \0

PAGE 99

90 opening ratios of 0.3, 0.5 and 0.7 are used. The exit velocity profiles for different cases are shown in Figures 22, 23, and 24 and Table 10. It seems that the computational domain was not long enough for the flow with r = 0.3. Table 11 presents the computed pressure drop. Figure 25 plots the computed pressure drops for the flows with the opening ratios of 1, 0.7, 0.5, and 0.3. Using Eq. 2.11, the orifice pressure drop can be converted to discharge coefficient which can be compared with the experimental data observed by Nigro (1978) for pipe flow with a circular orifice. It is noted that the orifice pressure drop measured by Nigro is the pressure difference between two flange taps that were located at one inch on both sides of circular orifice. In this the nearest nodal points at one inch on both sides of the plate orifice were chosen to calculate orifice pressure drops. Comparison between the predicted pressure drop and experimental data is shown in Table 12 and Figure 26. Although energy dissipation of a plate orifice may be somewhat different from a circular orifice, these two sets of pressure drops capture a similar trend and even numerical deviations are fairly small.

PAGE 100

0.9 0.8 0.7 0 0.6 :::> :J >. 0 5 u 0 Ill 0.4 > 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 Distance y/B 0 Re=50 + Re=100 <> Re=200 1::. Analytical Figure 2"2 comparison between the computed exit. velocity profile and the analytical Profile for \D -Opening Ratio = 0.7 -. .

PAGE 101

1 0.9 0.8 0.7 0 0.6 ::::> "-.. ::J >. 0.5 u 0 (I) 0.4 > 0.3 0.2 0.1 0 0 0 Re=SO I 0.2 0.4 0.6 0.8 Distance y /8 + Re=100 <> Re=200 h. Analytical Figure 23 between the computed exit velocity profile and the analytical profile opening Ratio = 0.5 1 U) N

PAGE 102

0.9 0.8 0.7 0 0.6 :J ......... :J >. 0.5 u 0 Gl 0.4 > 0.3 0.2 0.1 0 ,. 0 0 Re=50 0.2 0.4 Distance y /B + Re=100 <> Re=150 0.6 0.8 A Analytical \D w

PAGE 103

94 Table. 10 Comparison between the Computed Velocity Profiles and the Analytical Profile for cases with Contraction. =========================================================== Exit Velocity Profile ujUc Analytical d/B = 0.7 d/B = 0.5 ujUc =========================================================== Distance Reynolds number yjB 50 100 200 50 100 200 =========================================================== o.oo o.oo 0.00 o.oo 0.00 o.oo o.oo o.oo 0.05 0.09 0.08 o.o8 0.09 0.09 0.12 0.10 0.10 0.17 0.17 0.16 0.19 0.19 0.23 0.19 0.20 0.34 0.33 0.32 0.37 0.37 0.40 0.36 0.30 0.49 0.48 0.48 0.52 0.52 0.51 0.51 0.45 0.68 0.66 0.66 0.69 0.69 0.69 0.70 0.60 0.82 0.84 0.84 0.84 0.84 0.84 0.84 0.80 0.92 0.93 0.94 0.96 0.96 0.95 0.96 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 =========================================================== ================================================== Exit Velocity Profile u;uc d/B = 0.3 Analytical Profile ujUc ================================================== Distance Reynolds Number yjB 50 100 150 ================================================== o.oo o.oo 0.00 0.00 0.00 0.05 0.12 0.11 0.08 o.1o 0.10 0.28 0.24 0.22 0.19 0.25 0.52 0.44 0.44 0.44 0.40 0.70 0.57 0.60 0.64 0.70 0.87 0.85 0.84 0.91 1.00 1.00 1.00 1.00 1.00 ================================================== Note: Uo = entrance veloctiy. Uc= centerline velocity. X = horizontal distance. y = distance from the B = half width between two plates.

PAGE 104

Table 11 Pressure Drop of the Flow with a contraction. ============================================== b/B Ratio Reynolds Number 0 0.3 0.5 0.7 ============================================== 50 75 100 150 200 0.50 0.36 0.25 0.18 0.14 10.31 9.75 10.03 --2.94 0.98 2.74 0.84 2.89 0.81 ============================================== Note: --= not computed d/B = opening ratio between two plates computed in this study \0 ur

PAGE 105

p uo a. e Q "QI Ill, Ill Gl 10 5 d/8 = O.J d/8 = o.s )( )( d/8 1:1: 0.7 1 d/8 = 1 ---Reynolds Number Figure Computed pressure drops U) 0\

PAGE 106

Table 12 Comparison of Discharge Coefficients Between the Computed Results and Experimental Data. ============================================================ d/B Ratio d/D Ratio Reynolds Number 0.3 0.5 0.7 0.3 0.5 0.7 ============================================================ 50 100 150 200 0.70 0.72 0.71 0.71 0.74 0.72 0.73 0.78 0.80 0.70 0.71 0.70 0.69 0.71 0.73 0.72 0.71 0.74 0.79 0.81 0.82 ============================================================ Note: --= not computed d/B = opening ratio between two plates computed in this study. d/D = opening ratio in a circular pipe observed in lab. \0 -.,J

PAGE 107

0.9 ..., ; 0.8 ..... u ..... Ql 0 u CD 00 as .c: o. 7 ..... Q X A + .. d/B 0.3 o.s 0.7 experimental data d/8 .. 0.7 d/8 .. 0.5 d/B 0.3 10 20 40" 60 100 200 400 600 Reynolds Number Figure 26 Comparison of.discharge coefficients between computed results and experimental \0 co

PAGE 108

CHAPTER SIX CONCLUSIONS A finite element numerical model using quadratic isoparametric elements was developed to analyze a Poiseuille flow with and without a contraction. In this study, the unsteady flow governing were converted into stream function and vorticity transport equations. A general numerical form of the Poisson equation was developed to integrate finite element equations. Initial conditions were determined by potential flow and boundary conditions included a uniform incoming flow and non-slip condition on the plate. The range of computational dimensionless time increment was from 0.0025 to 0.01. Flow cases studied include the opening ratios ranging from 0.3 to 1 and Reynolds number varying from 50 to 200. In general, this numerical approach provides stable computations for flow to reach its steady when Reynolds number is less than 200. It was found that the computational domain of L/B =. 8.7 gives, in general, reasonably good predictions compared with analytical solutions. Although for high Reynolds number flows the exit velocity profile may slightly deviate from the exact velocity profile, the predicted friction loss and

PAGE 109

orifice loss still give good agreements to theoretical solutions and laboratory data. From the comparisons and discussions above, this numerical model is capable of computing flow patterns and predicting energy losses for a two dimensional flow with Reynolds number less than 200 and orifice opening ratio 100 larger than 0.3. This range may cover most blood flows in the artery and oil movements in lubrication devices. To extend its application to high Reynolds number flows, some improvements are recommended as follows: 1. Efficiency of Computation: In this study, for the case with the opening ratio of 0.7 and Reynolds number of 50, it took 4 hours of computer time to reach the steady state. Requirement of such a longcomputation time may be due to the use of quadratic isoparametric elements. If the linear triangular elements were used, it might reduce the computer time to one half of that used in this study. In addition, the values of wall vorticities play an importance role in the determination of numerical convergence. The. value of the wall .vorticity, as shown in Eq. 4.54, is proportional to the reciprocal to B 2 At each time step, the wall vorticity has to be updated first. If B is too large, the value of the wall vorticity becomes too small. It decelerates the rate of convergence; as a result, it requires more computer time to reach the steady state. Conversely, if B is too small, the numerical procedure may

PAGE 110

101 become unstable due to the abrupt increment in the values of the wall vorticity. Therefore, how to select the size of wall elements and determine the value of 8 are subjects for the future studies. 2. Expansion of Its Usefulness: It has been found that when Reynolds number is greater than 200 or the opening ratio is less than 0.3, this numerical model becomes unstable and the steady state can never be reached. As indicated in the beginning of this study, the development of this model was aimed at the application of the finite element scheme to the numerical modeling of larynx flows. How to extend the application of this numerical model from lower Reynolds number flows to turbulent flows may present another challenge to the future study.

PAGE 111

REFERENCE 1. Anderson, D. A., Tannehill, J. c., and Pletcher, R. H., Computational Fluid Mechanics and Heat Transfer, McGrawHill, N.Y., 19a4. : Chaps 1, and 6. 2. Baker, A. J., Finite Element Computational Fluid Mechanics, Hemisphere Publishing, N.Y., 19a3. : Chaps 1, 2, 3, 4, and 5. 3 .. Baker, A. J. ''Finite Element Solution Algorithm for Incompressible Fluid Dynamics" in Finite Elements in Fluids 2. Ed. Gallagher, Oden, Taylor, and Zienkiewicz. John Wiley, London, 1975 : 67:a2. 4. Bathe, K. J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, N.J, 19a2. : Chap 5. 5. Blevins, R. D., Applied Fluid Dynamics Handbook, Nostrand Reinhold, N.Y., i9a4. : Chaps 5, 6, and 7. 6. Cheng, R. T., "Numerical Solution of Navier-stokes Equations by Finite Element Method" The Phy. of Fluids 12, Dec. 1972. : 209a:2105. 7. Chung, T. J., Finite Element Analysis in Fluid Dynamics, McGraw-Hill, 197a. : Chaps 1, 2, 3, 4, and 5. a. Connor, J. J. and Brebbia, c. A., Finite Element Techniques for Fluid Flow, Newnes-Butterworths, LOndon, 1976. : Chap 9. 9. Cook, R. D., concepts And Applications of Finite Elenient. Analysis, John Wiley and Sons, N.Y., 19a1. : Chap 5. 10. Daily, J.W. and D.R.F., Flyid Dynamics, Addison-Wesley, Reading, Mass., 1966 : 137:143. 11. Daugherty, R. L., Franzini, J. B., and Finnemore, Fluid Mechanics with Engineering Applications, McGraw-Hill, N.Y., 19a5. : Chaps 1, 7, a, and 12. Desai, c. s., Elementary Finite Element Method, Prentice-Hall, Englewood N.J., 1979. :Chap 1. 13. Fox, R. w. and McDonald, A. T., Intro. to Fluid Mechanics, John Wiley and Sons, N.Y., 19a5 :Chaps 6, and a.

PAGE 112

14. Gerald, c. F., Applied Numerical Analysis, AddisonWesley, Reading, Mass., 1970 : Chaps 7, and 9. 103 15. Gerhart, P. M. and Gross, R. J., Fundamentals of Fluid Mechanics, Addison-Wesley, Reading, Mass., 1985. : Chap 7. 16. Huebner, K. H., The Finite Element Method for Engineers, John Wiley and Sons, N.Y., 1975. : Chap 9. 17. Idel'chek, I. E., Handbook of Hydraulics ResistanceCoefficients of Local Resistance and of Friction, u.s National Technical Information Service, 1960. 18. Janna, W. s., Intra. to Fluid Mechanics, PWS Publishers, Boston, Mass., 1987 : 554:557. 19. Jones, o. c., "An Improvement in the Calculation of Turbulent Friction in Rectangular Ducts" ASME J. of Fluids Eng., June 1976. 20. Kawahara, M., Yoshimura, N., and Nakagawa, K., "Analysis of steady Incompressible Viscous Flow" in Finite Element Methods in Flow Problems. Ed. Oden, Zienkiewicz, Gallagher, and Taylor. The University of Alabama Press, Huntsville, Ala, 1974 : 107:120. 21. Kersten, R. D., Engineering Differential Systems, McGraw-Hill, N.Y., 1969 : Chaps 3, and 7. 22. Kreyszig, E., Advanced Engineering Mathematics, John Wiley and Sons, N.Y., 1983 : 412:447. 23. Kwon, o. K. and Pletchcher, R. H. "Prediction of the Incompressible Separated Boundary Layers Including Viscous-Inviscid Interaction" J. Fluids Eng. 101,1979.: 466 :472. 24. Longwell, P. A., Mechanics of Fluid Flow, McGraw-Hill, N.Y., 1980 : Chap. 25. Nigro, F. E. B., et al., 11A numerical study of the Laminar Viscous Incompressible Flow Through a Pipe Orifice" J. of Fluid Eng. 100, 1978. : 467:472. 26. Olson, M.D. "Variational Finite Element Methoc;is for TwoDimensional and Axisymmetric Navier-Stokes Equations" in Finite Elements in Fluids 1. Ed. Gallagher, Oden, Taylor, and Zienkiewicz. John Wiley, London, 1975 : 57:72.

PAGE 113

27. Pao, R. H. F. Fluid Mechanics, John-Wiley and Sons, N.Y., 1961. 28. Pin Tong, "On the Solution of the Navier-Stokes Equations in Two-Dimensional and Axial-Symmetric Problems" in Finite Element Methods in Flow Problems. Ed. Oden, Zienkiewicz, Gallagher, and Taylor. The University of Alabama Press, Huntsville, Ala, 1974 57:66. 104 29. Segerlind, L. J., Applied Finite Element Analysis, John Wiley and Sons, Canada, 1984. 30. Scherer, R., Laryngeal Fluid Mechanics:Steady Flow Considerations Using Static Models, University of Iowa, diss., 1981. 31. Streeter, v. L., Fluid Dynamics, McGraw-Hill, N.Y.,1948 : Chap 1, 2, 10, 11, and 12. 32. Vennard, J. K., "One-Dimensional Flow" in Handbook of Fluid Dynamics. Ed. Streeter. McGraw-Hill, N.Y., 1969 : Chap 3. 33. White, F. M., Viscous Fluid Flow, McGraw-Hill, N.Y., 1974 : Chap 4. 34. Zienkiewics, o. c., The Finite Element Method, McGrawHill, London, 1977 : Chaps 1,3,7, and 8. 35. Zienkiewics, o. c. and Cheung, Y. K., "Finite Element in the Solution of Field Problems", The Engineer, 1965 : 507:510.

PAGE 114

APPENDIX A THE DERIVATION OF ELEMENT EQUATION BY GAUSS-GREEN THEOREM The Gauss-Green divergence theorem in a xy-plane states (Kreyszig, 1983): Let R be a closed bounded region in the xy-plane whose boundary r consists of finitely many smooth curves. Let U(x,y) be a vector function which is continuous and has continuous first partial derivatives in some domain containing R. Then, we can write I I R di v u dxdy = I r u n dr (Eq. A.l) where div = o +-o-6X oy and n = the outer unit normal vector. If we write u and n in terms of unit components, namely ( Eq. A. 2) n = cose i + sine j (Eq. A 3) Substituting Eq's .A.2 and A.3 into Eq. A.l gives + (Eq. A. 4) The integral equation we have to deal with is ( ) I T 02'" (R e } = -[N] (Ox _,.._2 + Dy Gp + Q)dA ox oy

PAGE 115

Using the product rule for differentiation in the above equation gives 106 { R (e) } = -J [ 0 _!:____ ([ N] T ocp ). + 0 ( [ N] T ocp ) ] dA A x ox ox Y oy oy o[NJT ocp o[NJT ocp + JA( 0 -+ 0 +G[N]Tcp + Q[N]T)dA x ox ox Y oy oy ( Eq. A. 5) Applying Eq. A.4 to the first integral of the right-hand side of Eq. A.5, it becomes = ocp ocp Jr( o [N]T -cose + oy [NJT -sine )dr x ox oy (Eq. A. 6) Eq. A.6 can be transformed to its final form using the following relationship. (Eq. A. 7) Substitution and rearranging yield J T ocp ocp {R(e)} =-r.( ox [NJ -cose + oy [N]T -.-sine )dr ox oy JAc ox o[N]T o[NJ o[N]T o[NJ )dA ] {t(e)) + [ + 0 ox ox y oy oy + ( JAG[N]T[N] dA ) { t (e)} -JAQ[N]T dA (Eq. A. 8)

PAGE 116

107 As to a two-dimensional time-dependent problem, the residual integral is {R(e)} = -J [N]T ( cp(x,y) + Qgf )dA (Eq. 4.19) D + o2tn where cp(x,y) = D Gcp x ox2 Y oy2 Comparing Eq. 4.19 with Eq. 4.13, the only additional term is I ( [N]T )dA ot ( Eq. A. 9) Assume that ocpjot varies quadratically within any element. ocpjot can be expressed as = N ;. (e) ot i 'i!i (Eq. A.10). where t = nodal values of ocpjot. Substituting Eq. A.10 into Eq. A.9 gives (Eq. A.11) Adding Eq. A.ll to Eq. A.S, the final form of the residual integral of a time-dependent field problem is {R(e)} ocp ocp = -J ( D [N]T -case + Dy [N]T -sine )dr r x ox oy + r JAc o[N]T o[NJ o[N]T o[Nl {t(e)} D + D )dA ] X ox ox y oy oy

PAGE 117

+ ( JAG[N]T[N] dA ) {t(e)} JAQ[N]T dA + ( J [N]T[N] dA ) {;(e)} 108 A.12)

PAGE 118

APPENDIX B THE INTEGRATION OF SHAPE FUNCTIONS Since the shape functions are expressed in terms of natural coordinates, to perform integration, the change of variables is required. In case of a double integral J I f(x,y)dxdy R the new variables of integration a,B can be introduced by setting x = x(a,B), and y = y(a,B) where the functions x(a,B) and y(a,B) are continuous and have continuous first derivatives in some region R* in the aB-plane so that each point (a,B) in R* corresponds to a point [x(a,B), y(a,B)] in Rand conversely. Kreyszig (1983) indicated that I I f(x,y)dxdy =I I f[x(a,B),y(a,B)] R R* o(x,y) I dadB o(a,B) (Eq. B.l) where o(x,y)jo(a,B) is the Jacobian matrix, which can further be expressed as [Jl = o(x,y) o(a,B) = ox OQ ox 613 (Eq. B. 2) To prove Eq. B.l, the two coordinate systems as

PAGE 119

110 shown in Fig 27 were considered (Chung; 1978). The directions of the cartesian coordinates x-y and the arbitrary nonorthogonal (possibly curvilinear ) isoparametric coordinates a-B are expressed by qnit vectors i1 i2 and the tangent vectors t1 t2 respectively. The relationship between i and t are t = 2 ox oy --oa ox oB The differential area (shaded) is ox _l!y oa oa o dadB or dxdy i3 = I det[J] I dadB i3 Thus, Eq. B.l can be obtained. ox ..21.. oB oB 0 (Eq. B.J) For an isoparametric quadratic quadrilateral element, the transformation equations are given by y = N. (a,B) Y. l. l. where N1 are shape functions, x.,Y. are the global l. l. coordinates of the nodes. The Jacobian matrix can be written as the matrix product as

PAGE 120

i \ 1 1 1 i__ y X Figure 27 Coordinate transformation

PAGE 121

112 x1 y2 oN1 oN2 oN3 oN4 oN5 oN6 oN7 oN8 x2 y2 oa oa oa oa oa oa oa oa XJ YJ [J] = oN1 oN2 oN3 oN4 oN5 oN6 oN7 oN8 x4 y4 oB oB oB os oB as a a as xs y5 x6 y6 x7 y7 xa Ya (Eq. B.J) in which oN1 joa = \(1-B)(2a+B), oN2 ;oa = and so on. In Eq. B.J, the terms of oNifoa and oNi/68 can be evaluated. The global coordinates of nodes, i.e., X. and l. Yi, are already known. Therefore, the Jacobian matrix for a element can be obtained, which is a 2x2 matrix. Substituting the determinant of the Jacobian matrix into the integrate equations, i.e., Eq's 4.35 and 4.36, the integral without the partial derivative terms can be performed. For those integral containing the partial derivative terms, we need to change the local derivative in terms of a and B to global derivatives in terms of x and y. Applying the chain rule for differentiation toa shape function Ni(a,B)t it gives oNi(a,B) oNi(a,B) 6x oNi(a,B) oy = --+ oa ox oa oy oa oNi(a,B) oNi(a,B) ox oNi(a,B) oy oB = ox ""68 + oy 08 or in a matrix form '. ;

PAGE 122

113 i oN. ox oy oN. l. (a,B) l. (a, B) oa 6a ox = oN. ox oy oN. l. (a, B) l. (a,B) oB 68 68 oy (Eq. B.4) The coefficient matrix is the Japobian matrix [J]. To find the global derivatives we invert [J] and write oN. oN. l. (a,B) l. (a,B) 6x -1 oa = [J(a,B)] (Eq. B.S) oN. oN. l. (a,B) l. (a,B) oy oB However, an explicit form of [J]-1 pannot be obtained (Zienkiewicz 1977). The numerical methods have to be used. In this study, the Gauss method is recommended because it has been proved to be in finite element method. Gauss-Legendre Quadrature To approximate the integral 1 R = J f(x) dx -1 (Eq. B. 6) we can select f(x) at the midpoint, i.e., f(O) and multipty by the length of the interval, i.e., 2. Thus we find R 2f(O). This result is an exact solution if the function y is a straight line. Generalization of Eq. B.6 yields 1 R = J f ( x) dx :E ( w i f (xi ) ) -1 ( Eq. B. 7) where = number of points selected. Thus, to approximate

PAGE 123

114 R, we evaluate f(xi) at each location xi, multiply f(xi) by an appropriate weighting value Wi' and sum the resultants up. The Gauss-Legendre quadrature locates the sampling points to achieve the greatest accuracy. Sampling points are located symmetrically with respect to the center of the interval. Symmetrically paired points have the same weighting value Wi. In general, a polynomial of degree 2n-1 is integrated exactly by n sampling points. further R = Applied to two-dimensional be expressed as 1 1 1 J_1J-1f(a,B) dadB I_1 problems, Eq. B.7 may n [,:E(Wj f (aj,B))] dB ]=1 m n :E :E [W.W.f(a.,B.)] i=1 j=1 1 J J 1 (Eq. B.8) The locations of sampling points and weighting coefficients can be found in Reference 4, a, 9, and 29.

PAGE 124

Symbol A a a B b c g hL [J] k K L N. l. p Re APPENDIX C DEFINITION OF SYMBOLS Definition Cross-sectional area The width of a rectangular duct The cross sectional area of orifice Half of the width between two plates The height of a rectangular duct Discharge coefficient Darcy friction factor Pipe diameter Hydraulic diameter Half of the opening width at the orifice Darcy friction factor for pipe flow Acceleration of gravity Head loss Jacobian matrix Friction coefficient ( k = f Re ) Loss coefficient Length of conduit Shape function Pressure Pressure difference Flow Reynolds number

PAGE 125

r t u u v v X y w. 1 Q B r r 1 e a Reynolds number for the flow in a rectangular duct Hydraulic radius Opening ratio Time 116 Inflow velocity for a Poiseuille flow Mean flow velocity through the orifice Velocity in x-direction Cross sectional average flow velocity Velocity in y-direction Cartesian coordinate, usually parallel to main flow direction Cartesian coordinate, usually perpendicular to main flow Weighting function Natural coordinate Natural coordinate Element boundary Shear stress Absolute viscosity Density Nodal value The angle to the outward normal Rotatoional velocity Stream function Vorticity

PAGE 126

APPENDIX D THE FORTRAN SOURCE CODE {

PAGE 127

c program TDICF common/elmatx/esmC8,8) ,efCB> ,phiCB) ,nsCS) ,wsmC8,8) ,p(200) common/matl/ge,qe,ivor,vk,wvC600) common/pdxy/vnCB),pnxC8l,pnyC8>,xxC8),yy(8),xd,yd,det common/crd/xc(600),ycC600),velx(600),velyC600) common/tle/title(20) common/av/aC40000),jqf,jqsm,np,nbw dimension nel(200,8),ick(600),w(8),q(200),sv(600) dimension iwc(l50,2),1pc(20),u(8),v(8) dimension ibc(l50,2),bc(l50),cc(40000) data in/60/,io/61/ openCin,f1le='FLOW.INP',status='old',form='formatted', $ access='sequential') open( io,file= 'FLOW. OUT' ,status= 'unknown' ,form=' formatted', $ access='sequential') c INPUT OF THE TITLE CARD AND PARAMBI'ERS c c writeC",l) 1 format(/lx,'Plaese enter Reynolds number :') readC",">re vk=l./re wr1teC",2) 2 format(/lx,'T1me interval for every computational step?') readC","Jdeltat readC1n,5)t1tle 5 format(20a4) read(in,")np,ne,nbc,nwc,npc read(in,")Cxc(i),i=l,np) read(in,")(yc(i),1=l,np) writeCio,8)title,np,ne,nbc,nwc,npc,re 8 format(////10x,20a4//lOx,5hNP =,16/lOx,ShNE =,16/ $ lOx,ShNBC =,i6/10x,ShNWC =,i6 $ /lOx,ShNPC =,16//lOx,ShRE =,16) do 9 kk::l,ne read(in,")n,(ne1Cn,i),1=1,8) 9 continue read(in,") (1bc(i,l),ibcC1,2),bc(1),1=l,nbc) read(in,"> CiwcCi,l),iwcCi,2),1=l,nwc) read(in,") (ipc(i),i=l,npc) c CALCULATION OF BANDWIDTH c nbw::Q do 20 kk,;.l,ne do 22 1=1,8 22 ns(i)=nel(kk,i) lk=7 do 21 1=l,lk ij=i+l do nj=ij,8 nb=iabs(ns(i)-ns(j)) ifCnb.le.nbw) qoto 21 inbw=kk nhw=nh 21 continue 20 continue nhw=nhw+l writeCio,24) nhw,inbw 24 format(//10x,l2hBANDWIDTH IS,I4,11H IN ELEMENT,I4) 118

PAGE 128

c c INITIALIZATION OF THE COLUMN VECTOR A() c -jgf=np jgsm=jgf+np j1=jend-jgf if(jend.gt.40000)then 25 format(l0x,30hDIMENSION OF A VECTOR EXCEEDED/lOX, $ 20HEXECUTION TERMINATED) goto 100 end if c -, c GENERATION OF THE SYSTEM OF EQUATIONS c c do 26 i=l,np 26 wvCi)=O. istep=O ivor=O istop=O 30 continue do 27 i=l,jend a(i)=O. cc(i)=O. 27 continue do 32 kk=l,ne do 31 1=1,8 ns(i)=nel(kk,i) 31 j=ns(i) ge=O. qe=O. ifCivor.eq.l)then qe=qCkk) ge=l. end if iftivor.eq.2)then ge=O. qe=pCkk) end if CALL ELSTMX C KK) ----c DIRECT STIFFNESS PROCEDURE c c:: do 33 i=1,B ii=ns(i) aCjgf+ii)=aCjgf+ii)+ef(i) do 34 j=l,B jj=ns(j)+l-11 if(jj.1e.O) goto 34 cc(ji)=cc(ji)+wsm(i,j) 34 continue 33 continue 32 continue c:: MODIFICATION AND SOLUTION OF THE SYSTEM OF EQUATIONS c:: if(ivor.ne.l)then if(ivor.eq.2)goto 75 119

PAGE 129

do 40 i=l,nbc if(ibc(i,2).eq.O)then ib=ibc(i,l) bv=bc(i) end if CALL MODIFYCib,bv) 40 continue CALL DCMPBD CALL SLVBD do 56 i=l,np 56 sv(i)=a(i) 98 continue do 76 i=l,np velx(i)=O. 76 vely(i)=O. do 70 kk=l,ne do 71 1=1 ,8 ns(i)=nel(kk,i) j=ns(i) ,. 11 phi( i> =a( j) ELGRAD ( KK) 70 continue 81 do 80 kk=l,ne do 81 1=1,8 ns(i)=nel(kk,i) j=ns(i) xx(i)=xc(j) yy(i)=yc(j) u(i)=velx(j) v(i)=vely(j) continue gdxu=O. gdxv=O. gdyu=O. gdyv=O. call pderv(O.,O.) do 83 j=l,8 83 gdyv=gdyv+pny(j)*v(j) continue 80 continue ivor=2 goto 30 75 continue do 54 i=l,npc ib=ipcCi) bv=lOO. CALL MODIFY{ib,bv) 54 continue CALL DCMPBD CALL SLVBD .. '.'. 44 format(//lOx,'step =',f8.4/) write(io,45) 45 formatCllx,'NO.',lx,'Streamline',4x,'Vorticity' ,5x'Pressure', $ 9x, 'Vel-x' ,6x, 'Vel-y' I) 120

PAGE 130

write(io,43) (i,sv(i),wv(i),a(i),velx(i),velyCil,i=l,np) 43 format(l0x,13,2x,fl0.5,2x,fl0.5,2x,fl2.2,4x,f9.5,2x,f9.5) end if if(istep.qt.3000)qoto 100 if(istop.eq.l)goto 100 do 47 i=l,nwc if(iwc(i,2).ne.Olthen j=iwc(1,1) k=iwc(1,2) & end if 47 continue do 51 kk=l,ne do 52 1=1,8 ns(i)=nelCkk,i) j=ns(i) xxC1)=xc(j) yy(i)=yc:(j) w(i)=wv(j) 52 phi(1)=sv(j) qdxphi=O. qdyphi=O. qdxw=O. qdyw=O. PDERV(O.,O.) do 53 j=l,8 53 continue 51 continue ivor=l istep=istep+l goto 30 end if do 60 i=l,np sum=o. ilimi t=nbw+i-1 if(i11mit.gt.np)11imit=np do 61 j=l,111mit if(i.le.j)then k=j-1 else k=i-j+1 ij=iij+jgsm ifCiij.qt.nnc)goto 61 end if 61 continue aCjgf+i)=aCjgf+i)+sum/deltat 60 continue do 6.2 1=1 ,nne 121

PAGE 131

c j=i+jgsm a(j)=a(j)+cc(j)/de1tat 62 continue do. 63 i=1,nwc if(iwc(i,2).eq.O)then ib=iwc(i,1) hv=o. else ib=iwcC i, 1) bv=wv(ib) endif CALL MODIFY(ib,bv) 63 continue CALL DCMPBD CALL SLVBD do 90 i=1,np ick(i)=l def=aCi)-wv(i) deft=def/de1tat if(deft.lt.0.04) ick(i)=O 90 continue do 99 i=l,np if(ick(i).eq.1)then do 64 i=1,np 64 wv(i)=a(i) ivor=O goto 30 end if 99 continue istop=l goto 98 100 continue stop end SUBROUTINE ELSTMXCKK) c EVALUATES THE ELEMENT STIFFNESS MATRIX AND ELEMENT FORCE VECTOR c common/e1matx/esm(8,8),ef(8),phi(8),ns(8),wsm(8,8),p(200) common/matl/ge,qe,ivor,vk,wv(600) common/pdxy/vn(B),pnx(B),pny(B),xx(B),yyCB),xd,yd,det common/ipts/vxC9r,vyC9) ,wc(9) common/crd/xc(600),yc(600),velx(600),ve1y(600) dimension w(B),ff(8,8) do 1 1=1,8 ef(i)=O.O do 1 j=l,8 ff(i,j)=O. wsm(i,j)=O.O 1 esm(i,j)=O.O do 10 1=1,8 j=ns(i) w( i) =wv( j) xx(i)=xc(j) 10 yy(i)=yc(j) 122

PAGE 132

CALL INGPTS do 3 ii=l,9 call pderv(vx(ii),vy(ii)) do 2 i=l,B do 2 j=l,B if ( 1 vor. eq. 0) ff ( i, j) =ff ( i, j) +vn( 1) j) ( i1) 2 continue 3 continue if(ivor.eq.l)then. do 4 i=l,B do 4 j=l,B 4 continue end if if(ivor.eq.O)then do 5 i=l,B do 5 j=l,B 5 continue end if return end SUBROUTINE ELGRADCKK) comrnon/elrnatx/esrn(8,B),ef(B),phi(8),ns(B),wsm(B,8),p(200) common/pdxy/vn(8),pnx(8),pny(8),xx(B),yy(8),xd,yd,det cornmon/ipts/vx(9),vy(9),wc(9) cornmon/crd/xc(600),yc(600),velx(600),vely(600) dimension xq(9),yq(9),gdx(9),gdy(9) xq/-l.,O.,l.,l.,l.,0.,-1.,-1.,0./ data yq/-l.,-l.,-l.,O.,l.,l.,l.,O.,O./ data io/611 do 1 i=l,B j=ns(i) xx(i)=xc(j) yy(i)=yc(j) 1 continue do 10 1=1,9 gdx(i)=O. gdy(i)=O. CALL PDERV
PAGE 133

c velx(j)=Cvelx(j)+qdx(1))/2. endif ifCvely(j).eq.O.)then vely(j)=qdy(i) else vely(j)=(vely(j)+qdy(1))/2. end if 15 continue return end SUBROUTINE INGPTS common/ipts/vx(9),vy(9),wc(9) dimension a(3),b(3) data data b/5.,8.,5./ c GENERATION OF THE NINE INTERGATION POINTS c c n=O do 1 1=1,3 do 1 j=l,3 n=n+l vxCn)=aC1) vyCn)=aCj) 1 return end SUBROUTINE PDERV(Xl,X2) C EVALUATES THE JOCOBIAN TRANSFORMATION MATRIX C AND USES THE INVERSE JOCOBIAN MATRIX TO OBTAINTHE DERIVATIVES C OF THE SHAPE FUNCTIONS WITH TO X AND Y. c common/pdxy/vnn(8),pnxC8),pny(8),x(8),y(8),xd,yd,det common/derv/vn(8),pnsCB),pne(8) real jocb(2,2) do 1 1=1,2 do 1 j=l,2 1 jocb(i,j)=O.O xd=O.O "yd=O.O CALL QDSHFNCXl,X2) do 2 1=1,8 2 .. jocbC 2, 2 > =joc:bc 2, 2 > +pnec 1 > b=jocb(l ,1) jocb(l,l)=jocb(2,2)/a jocb(l,2)=-jocb(l,2)/a jocb(2,1)=-jocb(2,1)/a jocbC 2, 2) =b/a det=abs(a) 124

PAGE 134

c do 3 1=1,8 vnn(i)=vn(i) pnx(i)=jocb(l,l)*pns(i)+jocb(1,2)*pne(i) pny(i)=jocbC2,1)*pns(i)+jocbC2,2)*pne(i) return end SUBROUTINE QDSHFNCSI,ETA) 125 c CALCULATES THE VALUE OF THE SHAPE FUNCTIONS AND THEIR DERIVATIVES c common/derv/vn(B),pns(B),pne(B) dimension sq(B),eqCB) data sq/-1.,0.,1.,1.,1.,0.,-1.,-1./ data eq/-1.,-1.,-1.,0.,1.,1.,1.,0./ do 5 1=1,7,2 vn(i)=.25*C1.+si*sq(i))*C1.+eta*eq(i))*(si*sq(i)+etalleq(i)-1.) pns(i)=.25*C1.+eta,i:eqCi))ll(2.,r,:si+eta*eq(i)*sqCi)) 5 pne( i) =. 25* C l.+si*sq( i)) *C 2. *eta+si*sq( i) lleq( i) )_ do 6 1=2,6,4 vn(i)=.5*(1.-s1**2>*C1.+eta*eq(i)) pns(i)=-si*C1.+eta*eq(i)) 6 pneCi)=.5*eqCi)*C1.-si**2) do 7 1=4,8,4 vnC i) =. 5* ( 1. +si*sq( i)) *< 1. -eta**2) pns(i)=.5*sq(i)*(1.-eta**2) 7 pne(i)=-eta*C1.+si*sq(i)) return end .. SUBROUTINE MODIFY(ib,bv) common/av/a(40000),jqf,jqsm,np,nbw data inl60/,io/61/ k=ib-1 do. 211 j=2,nbw m=ib+j-1 if(m.qt.np) qoto 210 -aCjgf+m)=aCjqf+m)-a(ij)*bv aCij)=O.O 210 if(k.le.O) qoto 211 kj=jgsm+(j-l)*np+k-Cj-l)*Cj-2)/2 a(jgf+k)=aCjqf+k)-a(kj)*bv a(kj)=O.O k=k-1 211 continue a(jqf+ib)=a(jqsm+ib)*bv 221 continue return end -SUBROUTINE DCMPBD common/av/a(40000),jqf,jgsm,np,nbw io=61 npl=np-1 do 226 i=1,npl mj=i+nbw-1 if(mj.qt.np) mj=np

PAGE 135

nj=i+1 mk=nbw if((np-i+1).1t.nbw) mk=np-i+1 nd=O do 225 j=nj,mj mk=mk-1 nd=nd+1 n1=nd+1 do 225 k=1,mk nk=nd+k ii=jgsm+i 225 a(jk)=a(jk)-a(in1)"a(ink)/a(ii) 226 continue return end SUBROUTINE SLVBD common/av/a(40000),jgf,jqsm,np,nbw npl=np-1 do 250 i=1,np1 mj=i+nbw-1 if(mj.gt.np) mj=np nj=i+1 1=1 do 250 j=nj,mj 1=1+1 250 continue a(npl=a(jqf+np)/a(jgsm+np) do 252 k=1,np1 i=np-k mj=nbw if((i+nbw-1).gt.np) mj=np-1+1 sum=O.O do 251 j=2,mj n=i+j-1 251 252 a(i)=(a(jgf+i)-sum)/a(jgsm+i) return end 126