|
UFL/COEL-90/016
TWO-DIMENSIONAL FINITE-DIFFERENCE MODEL
FOR ESTUARINE AND COASTAL FRONTS
by
Yu Luo
Thesis
1990
TWO-DIMENSIONAL FINITE-DIFFERENCE MODEL
FOR ESTUARINE AND COASTAL FRONTS
By
YU LUO
A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
UNIVERSITY OF FLORIDA
1990
ACKNOWLEDGEMENT
I would like to express my gratitude to my advisor and
committee chairman Professor Y. Peter Sheng for his
initiating, motivating, directing and financially
supporting my work and for his criticism and patience
during the writing process. I would also like to extend my
thanks and appreciation to my thesis committee members,
Professor Donald M. Sheppard, Professor Wei Shyy and
Professor Robert J. Thieke, for their comments about this
work and patience in reviewing this thesis.
I must thank Mr. Hongjun Li, my husband, for his
numerous important and valuable suggestions about this work
and for his general help on computer programming and during
the writing process.
My gratitude must also be extended to Dr. James
O'Donnell for his useful suggestions on this work.
The support of all my fellow graduate assistants and
postdoctoral associate is also deeply appreciated. I would
like to thank Dr. Donald Eliason, Mr. Hyekeun Lee, Mr.
Steve Peene, Mr. Xingjian Chen and Mr. Jei Choi. My
appreciation must also be extended to Mr. Suburna Malakar
ii
for his help in dealing with the VAX computers.
Finally, I would like to thank my family in China for
their encouragement throughout this study.
iii
TABLE OF CONTENTS
page
ACKNOWLEDGMENTS .......................................... ii
LIST OF FIGURES ........................................ vi
ABSTRACT .............................................. x
CHAPTERS
1 INTRODUCTION ..................... ................. 1
2 NUMERICAL SCHEME REVIEW ........................... 13
2.1 First-Order Upwind Scheme ................... 15
2.2 Flux-Corrected Transport (FCT) Scheme ...... 15
2.2.1 The First-Order Scheme ................ 18
2.2.2 The Second-Order Scheme ............... 18
2.2.3 The Combination-Hybrid Scheme ......... 19
2.3 Central-Difference Scheme .................... 20
2.4 Total Variation Diminishing (TVD) Scheme .... 21
2.5 Quadratic Upstream Interpolation for Convective
Kinematics with Estimated Streaming Terms
(QUICKEST) Method ........................... 22
2.6 Time Marching Method ...................... 24
2.7 Summary ........................................ 26
3 NUMERICAL MODEL ..................................... 27
3.1 Governing Equation .......................... 27
3.2 Numerical Model ................................ 30
3.2.1 Computation Domain ..................... 31
3.2.2 Descretization of Equations .......... 31
3.3 Finite-Difference Equation and Procedure ...... 34
3.3.1 Find i+lby Solving Continuity Equation.. 34
3.3.2 Solve u ?1 and v+1 ...................... 35
3.3.3 Solve S + ............................... 40
3.4 Stability ...................................... 42
4 ADAPTIVE GRID TECHNIQUE .......................... 44
4.1 Grid Adaptive Algorithm ..................... 45
4.2 Weight Function ................................ 47
4.3 Steady-State vs. Transient Problems ........ 50
4.4 Temporal Coupling ........................... 52
4.5 Verification of The Adaptive Grid System .....
5 NUMERICAL EXPERIMENTS ............................
53
60
5.1 "Breaking Dam" Problem ........................ 60
5.1.1 Physical Problems ....................... 60
5.1.2 Numerical Results ...................... 62
5.2 Flow with Linear Density Distribution ....... 64
5.2.1 Physical Problems ......................... 64
5.2.2 Numerical Results ....................... 69
5.3 Sudden Release Problem .......................... 74
5.3.1 Physical Problems ......................... 74
5.3.2 Numerical Results ...................... 77
5.4 Water Discharge Problem ......................... 79
5.4.1 Physical Problems ......................... 79
5.4.2 Numerical Results ........................ 91
5.5 The Application of Adaptive Grid Technique ... 100
6 CONCLUSION .......................................
APPEND IX ...............................................
REFERENCES ............................................
BIOGRAPHICAL SKETCH .....................................
116
119
124
128
LIST OF FIGURES
Estuarine front with foam line observed just
north of the lightering area inside
Delaware Bay from Bowman and Iverson(1977) .....
Thermal imagery of Point Beach power plant
from Scarpace et al.(1975) .....................
South Pass effluent into ambient water masses
during flooding and ebbing tides
from Wright and Coleman(1974) .......
The comparison of numerical schemes
at CFL=0.8 .........................
The comparison of numerical schemes
at CFL=0.2 .........................
QUICKEST method ....................
Results of QUICKEST method
(a) Linear equation ................
(b) Nonlinear equation ...........
Computation domain and grid system .
Initial mass distribution .........
Momentum balance principle ........
Case (1) on adaptive grid points 20
Case (2) on adaptive grid points 20
Case (3) on adaptive grid points 20 ...........
The comparison of solutions on uniform grid and
adaptive grid
(a) Case (1) on uniform grid 20 and 80
(b) Case (3) on uniform grid 20 and 40 .........
Initial State of Dam ..........................
1.2
1.3
2.1
2.2
2.3
2.4
........... 12
........... 16
........... 17
........... 23
3.1
4.1
4.2
4.3
4.4
4.5
4.6
5.1
...........
............
............
...........
5.2 "Breaking dam" problem
(a) Analytic solutions at unit intervals of time
(b) Numerical solutions corresponding to (a)... 65
5.3 The numerical solutions for the layer depth with
rotation effects in (a) t=10 and (b) t=25 ..... 66
.5.4 Nonratating "breaking dam" problem
at t=l and t=6 ................................. 67
5.5 With rotating effects, the gravitive
adjustment at t=l and t=3 ........................ 68
5.6 Density Distribution ......................... 69
5.7 Numerical time-sequence showing the isopycnals r
with p= .058x10-3g/cm3. a=15xl0-cm-.
at t=8.Os, 12.3s. 16.7s. 20.9s and 24.9s ....... 71
5.8 Corresponding to Figure 5.2-1, a time-sequence
showing the motion of the isopycnals
in laboratory experiments from
Simpson and Linden(1988) ...................... 72
5.9 A log-log plot of the angle 9 of inclination
of the isopycnals to the vertical against the
dimensionless t'==(Iga) /'2t ..................... 73
5.10 Laboratory Experiment Case 1 ..................... 75
5.11 Laboratory Experiment Case 2 .... ................ 75
5.12 Large Scale Bottom Plume ........................ 76
5.13 Large Scale Surface Plume ...................... 76
5.14 Shadowgraphs of a volume of salt water
with ho=7cm, H=49cm, xo=50cm, g'=47cms-2
collapsing into fresh water at 2.2s, 3.8s
and 6.0s after release ...... .. ... .... ........ 80
5.15 Shadowgraphs of a volume of salt water with
ho/H=1, ho=7ctm, x0=50cm, g'=47cms-2
collapsing into fresh water at 4.7s,
7.7s and 10.7s after release ..... ............. 81
5.16 Numerical solutions corresponding
to Figure 5.14 ................................ 82
vii
sta' En"eee-n, .-AIr^ li
3.17 Numerical solutions corresponding
to Figure 5.15 .................. ............. 83
5.18 Numerical solutions of Case 3 at 20hrs, 40hrs,
60hrs and 80hrs ................................. 84
5.19 Numerical solutions of Case 4 at 20hrs. 40hrs.
60hrs and 0Ohrs .................................... 85
5.20 Numerical solution for surface plume at 20hrs.
4Ohrs, 60hrs and SOhrs by Wang(1984) ............ 86
5.21 Numerical solution for bottom plume at 20hrs,
40hrs. 60hrs and 80hrs by Wang(1984) ............. 87
5.22 Numerical solution for surface plume at 20hrs.
40hrs, 60hrs and SO/rs by current model ......... 88
5.23 Numerical Solution for Bottom Plume at 20hrs,
40hrs, 60hrs and 80hrs Coresponding tp 5.3-12 ... 89
5.24 Water Discharge Problem ...................... 90
5.25 Water discharge problem with Fr=0.04
for overflow case at 2hrs, 4hrs. 6hrs and 8hrs
(a) Velocity vector ............................ 92
(b) Salinity contour ........................... 93
5.26 Water discharge problem with Fr=0.22 for
overflow case at 20hrs 40ihrs. 60hrs and 80hrs
(a) Velocity vector ........................... 94
(b) Salinity contour .............................. 95
5.27 Water discharge problem with Fr=0.057 for
underflow case at 4h/rs. 8hrs, 12hrs and 16hrs
(a) Velocity vector ............................. 96
(b) Salinity contour ........................... 97
5.28 Water discharge problem with Fr=0.64 for
underflow case at 4hrs, 8hrs, 12hrs and 16hrs
(a) Velocity vector ............................. 98
(b) Salinity contour ........................... 99
5.29 Two regular grids with grid size 80x10 and
40x10 for sudden release problem Case 4 ....... 102
5.30 Adaptive grid at 20hrs. 40hrs, 60hrs and 80hrs
with grid size 40x10 ......................... 103
viii
5.31 The salinity contour obtained on uniform grid
80x10 at 20hrs. 40hrs. 60/irs and SOhrs ......... 104
5.32 The salinity contour obtained on uniform grid
40x10 at 20hrs. 40hrs. 60hrs and SOhrs ......... 105
5.33 The salinity contour obtained on adaptive grid
40x10 at 20hrs, 40/hs. 60hrs and 80hrs ......... 106
5.34 Two uniform grids with grid size 50x10 and
25x10 for overflow water discharge problem
at Fr=0.22 ...... .... ... ......................... 107
5.35 Adaptive grid at 2hrs. 4-hrs. 6hrs and Shrs
with grid size 25x10 ............................. 108
5.36 The velocity vector plots at 2hrs. 4hrs, 6hrs and
Shrs on uniform grid 25x10 ..................... 109
5.37 The velocity vector plots at 2hrs, 4hrs, 6hrs and
8hrs on adaptive grid 25x10 .................... 110
5.38 The salinity contour obtained on uniform grid
50x10 at 2hrs, 4hrs, 6hrs and Shrs ............. 111
5.39 The salinity contour obtained on uniform grid
25x10 at 2hrs, 4hrs, 6hrs and Shrs ............. 112
5.40 The salinity contour obtained on adaptive grid
25x10 at 2hrs, 4hrs, 6hrs and Shrs .............. 113
5.41 The salinity contour obtained on uniform grid
25x10 at 2hrs. 4hrs. Chrs and Shrs by Upwind scheme. 114
5.42 The velocity vector plots at 2hrs, 4hrs, 6hrs and
8hrs on adaptive grid 25x10 by Upwind scheme ... 115
A.1 Sketch of characteristic method ............... 119
Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science
TWO-DIMENSIONAL FINITE-DIFFERENCE MODEL
FOR ESTUARINE AND COASTAL FRONTS
By
Yu Luo
December 1990
Chairman: Professor Y. Peter Sheng
Major Department: Coastal and Oceanographic Engineering
A two-dimensional numerical model is developed to
predict the properties of gravity currents and resolve the
sharp density gradient at the frontal zone in estuarine and
coastal waters. The governing equations, which are
vertically integrated Navier-Stokes equations of fluid
motion are solved by using a finite-difference scheme with
a staggered grid in space and the Euler forward-time scheme
for time integration. In order to capture the sharp
variation density gradient at the front zone accurately and
efficiently, this study employs a high resolution TVD
(Total Variation Diminishing) scheme, which is simplified
by "freezing" the local constants for complicated system
equations, and an adaptive grid method with weight function
correction. A series of analytical solutions and laboratory
experiment results is used to verify the accuracy of the
model. The model is then employed to investigate the
dynamics of fresh water or cool water discharge into a
large water body.
CHAPTER 1
INTRODUCTION
The dynamics of coastal and estuarine are often
locally dominated by shallow layers of buoyant water or
plumes. These layers are often established by the discharge
of fresh water from a river or bay into estuarine or
coastal water, and as another example, the cooling water or
wastewater from a power plant into a river or estuary. As
the discharging and receiving waters with different
densities (due to different salinities or different
temperatures) meet, they are often separated by a sharp
front, i.e. a narrow zone of sharp density gradient.
Fronts occur over a wide range of spatial scales. In
an estuary, a bottom salinity front (salt wedge) separates
the incoming ocean water from the estuarine water. Near the
river mouth, a surface salinity front (river plume)
separates the spreading, brackish water from the ocean
water. On shallow continental shelves, fronts are formed at
the transition between well-mixed and stratified waters.
And along the shelf break, a semi-permanent shelf-slope
front is formed during the winter months, separating cold,
fresh shelf water from warm, saline open-ocean water. Other
fronts are found on the continental shelf, separating a
zone of coastal water from oceanic water. Fronts also occur
on a larger scale in the deep ocean, between water masses
of different properties (Bowden, 1983).
There have been many observations of plumes and
fronts. The Mississippi River's discharge at South Pass
produces a plume into the Gulf of Mexico. The behaviour of
the plume. with its load of suspended sediment, has been
the subject of a number of investigations, including that
of Wright and Coleman(1974). The spreading river water has
a considerable influence on the coastal band of water out
to the edge of the continental shelf, as marked by the 200m
contour which, in the region of the delta, is about 50 km
from the coast.
The Columbia River flows into the Pacific Ocean as the
boundary between the states of Washington and Oregon, USA.
Duxbury(1965) found that its estuary is partially mixed as
a result of a high river discharge averaging 7,300 m3s-,
combined with strong tidal currents. The outflowing upper
layer forms a plume of relatively low salinity water which
can be distinguished over a large area of the adjacent
ocean. The behaviour of the plume depends largely on the
wind and wind-driven currents and shows a seasonal
variation.
Another example of river plume is from Bowman(1978).
The Hudson River plume, which is largely determined by
wind conditions, flows into the New York Bight.
Szekielda and Mitchell (1972) found heavy metal
concentrations three orders of magnitude higher than
background values in frontal zones along Delaware Bay for
chromium, copper, lead, mercury, silver, and zinc. Klemas
and Polis(1977) graphically illustrated how convergent
tidal currents in Delaware Bay captured and advected oil
slicks over period as short as one hour.
Scarpace et al.(1975) presented a study of the
structure of the thermal plume of the Point Beach power
plant in Lake Michigan. Again a major feature of the plume
was the sharp frontal boundaries in the near field.
Secondary fronts were also observed in the far field. These
fronts appeared as boundaries between areas of slightly
different temperatures.
Garvine and Monk(1974) reported the results of
observation of the hydrographic fields and the circulation
in the neighborhood of the Connecticut River front. They
found that the front essentially separates the brackish
plume water from the ambient Long Island Sound water and
that mixing between them takes place within 50m of the
front.
Since the gravity current is an interesting and
important problem of concern to scientists,
4
environmentalists and geographiers. there have been- many
attempts to model plumes and fronts. For example, Kao et
al.(1978) and Valentine and Kao(1983) examined the detailed
structure of the flowfield produced by a buoyant inflow
into a quiescent receiving basin. They reduced the
primitive equations into stream function and vorticity
equations. Which are then solved by using a conventional
over-relaxation method with the stream function for the
point source potential flow solution as an initial guess.
In their discussion of accuracy related to the advection
term, they defined a false diffusion term which resulted
in truncation error. They analyzed the effects of the false
diffusion by comparing the advective term, the physical
diffusion term and the false diffusion term at various
areas in plume region including frontal zone. They
concluded that the advective term is dominant over the
false and physical diffusion in the headwave region.
However, in the main established region behind the
headwave, the false diffusion becomes important. This
computational diffusion is probably large enough to -mask
the physical diffusion and viscous effect. Their results
showed that the currents reached a steady speed of
propagation and the layer behind the front appeared to
reach a steady finite depth. They discussed the influence
of the densimetric Froude number for the stratified gravity
currents in near field and reported that their results are
consistent with the observation that the frontal zone forms
a strong surface convergence. Furthermore, their model
prediction agreed quantitatively with data from the
Connecticut River Plume provided by Garvine and Monk(1974).
Fischer(1976) developed a numerical model for density
currents in estuaries. He pointed out that the main
difficulty for the numerical model lies in the numerical
treatment of the convection terms. The usual methods for
treating the convection terms (e.g., Forward Time/Centered
Space, Upstream Difference, Lax-Wendroff, ADI, Leapfrog or
Galerkin Finite-Element Method) break down when the
convection of fairly sharp 'density jumps is to be
calculated, as is the case for salt wedges. He described a
new method "Hermite-interpolation method," in which the
spatial density distribution is represented by IIermite
interpolation functions and the density value p and its
ap
local gradient are defined as time dependent variables.
ax
In his model, the advection term is neglected in the
momentum equation which is solved by an explicit leap-frog
method with an implicit treatment of the viscosity term. He
concluded that Hermite-interpolation method gave
satisfactory results in every case.
O'Donnell and Garvine(1983) and O'Donnell(1986)
developed an approach based on the assumption that the
plume may be divided into two dynamically distinct domains.
the inviscid region or plume body, and the frontal zone.
The plume body is considered to be the upper layer of a two
layered flow between which turbulent exchange with the
ambient fluid is negligible. The dynamics of the flow in
this region is dominated by gravitational spreading across
which flow properties are related by algebraic jump
conditions. In their two-dimensional model, they used the
"front patch technique." which employed the
bicharacteristic relationship and the characteristic cones
of the systems in the frontal zone. In the interior of the
plume, the two-step Lax-Wendroff scheme was employed. They
discussed and calculated the 'growth of a plume from a
radial source in a steady crossflow. It was found that the
nearshore areas of the plume are most sensitive to the
effects of rotation. The effect of the oscillation tidal
crossflow on the Connecticut River plume indicated that the
plume is diluted by vertical mixing soon after a change in
the direction of the crossflow.
Wang(1984) developed a primitive equation model to
examine the two-dimensional mutual intrusion of a gravity
current in vertical plain and the formation of a density
front. He used a staggered-grid, central-difference scheme
in space and a leapfrog scheme in time, except for the
advection terms in the salinity equation which were treated
with an upwind scheme. A mode-spliting technique (Sheng et
al.. 1978) was used to achieve computational efficiency.
Later Wang(1984) changed the upwind scheme in the salinity
equation to a diffusion-corrected advection scheme
(Smolarkiewicz, 1983. 1984) for solving the advection term.
He pointed out that the upwind scheme can generate large
numerical diffusion which can be significantly reduced by
using a dffusion-corrected scheme. His solution gave a
realistic description of the frontal structure and the mean
southward current by comparing with Garvin and Monk's(1974)
observations from the New England shelf-slope front.
Wang(1984) extended his two-dimensional model to a three-
dimensional one to study gravity currents produced by
instantaneous release of a buoyant fluid in a rectangular
channel with and without the Coriolis effect. His results
were compared with laboratory experiments provided by
Rottman and Simpson(1983).
Chao and Boicourt(1986) developed a three-dimensional,
primitive-equation model to study the onset of estuarine
plumes. Their numerical solution essentially followed the
scheme of Semtner(1974). Equations and boundary conditions
are finite-differenced in space to second-order accuracies
and leap-frogged in time with lagged friction and
diffusion. In case of outflow boundary condition, upwind
finite differencing is employed for the advection term. The
8
finite-difference scheme for the advection terms conserves
mass. energy, salt and salt variance. At each time step,
the vertically-averaged current and the deviation from the
mean are found separately. The former was obtained by a
five-point, sucessive over-relaxation of the vorticity of
the depth-averaged flow predicted via a vorticity equation.
The latter was found by integrating the momentum equations
in time while removing the barotropic components. Regular
averaging of two consecutive time steps was carried out to [
avoid the splitting of the developing solutions at even and [
odd time steps. If the ocean becomes gravitationally
unstable, instantaneous vertical mixing was carried out to
achieve a stable stratification. Their results showed that
the plume expands in the direction of propagation of the
coastally trapped waves after the freshwater release. The
seaward expansion of the bulge decreases, and the
undercurrent leaks out of the bulge and propagates behind
the nose of the bore. The three-dimensional structures of
the density-driven estuarine circulation and also of the
bore intrusion along the shelf have been identified.
All these models tried to obtain good results by using
a variety of techniques aimed to improve the accuracy and
efficiency of the numerical solution in vicinity of a
front. This thesis will present a study to further improve
the accuracy and efficiency of the numerical solution of
frontal problems. A two-dimensional numerical model is
developed by using a finite-difference method in which a
high resolution TVD (Total Variation Dissipation) scheme
for simulating advection terms is employed in both the
momentum and the salinity equations and an adaptive grid
technique is applied. To verify the accuracy of this model,
numerical solutions are compared with analytical solutions
and some laboratory experiment results. The capabilities of
this model in solving the realistic problems are
demonstrated with further numerical experiments.
This thesis consists of 6 chapters. An introduction is
given in Chapter 1. In Chapter 2. a brief review of several
numerical schemes for solving advection terms are
presented. A model equation (Bergur's Equation) is used to
verify those schemes separately and the comparisons of
solutions are presented. In Chapter 3, the derivation of
the model from the shallow water equations is presented.
The computational domain and the computational procedure
are explained. Chapter 4 describes an adaptive grid
technique. The effect of weight functions is discussed and
different test cases are evaluated to assess the
effectiveness and efficiency of the numerical grid
adaption. In Chapter 5. results of a series of numerical
experiments are presented. Comparisons with analytical
solutions and laboratory experiment results are presented.
10
Some more realistic problems are also solved by the model.
Chapter 6 gives the conclusion and suggestions.
:'-A,;.
- 5.* .
I _* 0
Figure 1.1 Estuarine front with foam line observed just north
of the lightering area inside Delaware Bay from
Bowman and Iverson(1977)
Figure 1.2 Thermal imagery of Point Beach power plant
from Scarpace et al(1975)
12
Gulf water -
S 00 1000 0 00 1000
met*'$ Momrs
Figure 1.3 South Pass effluent into ambient water masses
during flooding and ebbing tides from Wright and
Coleman(1974)
ii'
oor Ae r
FJ~ure1.3Souh Pss fflen int Ben w~e mtse
duringI odn n bbn =ds rm rgh n
Co-len (174
CHAPTER 2
REVIEW OF NUMERICAL SCHEMES
The mathematical model of gravity currents is governed
by Navier-Stokes equations. However, the numerical
solutions of Navier-Stokes equations are plagued with
difficulties caused by the numerical treatment of
advection terms. An ideal numerical scheme should have the
following properties:
1. Accuracy: Second order accuracy is usually
required.
2. Capibility to capture the discontinuities.
3. Efficiency.
4. Robustness: Able to handle different problems.
5. Convenience: Easy to program.
Unfortunately, no scheme can meet all of the above
requirements presently. Many of the schemes work well only
for the special application for which they have been
designed. Wang(1984) used an upwind scheme and a diffusion-
corrected scheme for the advection terms. O'Donnell(1986)
used the Lax-Wendroff scheme for the advections terms to
solve the two layer front problem.
In this chapter, a few schemes which are usually
employed to solve hyperbolic equations are disscused based
on a model equation, i.e., a one-dimensional advection-
diffusion equation:
atO+=F x (2.1)
where
u Dependent variables
1"
F u
y Viscosity coefficient
t Time
x Spacial coordinate
For inviscid flow, p=0. The equation becomes
Sat> (2.2)
In this equation, the characteristic speed is aF=u.
au
The initial conditions are
1 x=O, t=O
u= (2.3)
0 x>O. t=O
and the boundary conditions are
1 x=O
u=f (2.4)
0 x=l
All numerical computations are completed on a grid of
20 points. The Courant number, CFL=uy, is chosen 0.8 and
0.2.
2.1 First-Order Upwind Scheme
The first-order upwind scheme is written as
u('+=u,-A(OF',++( -2o0)F'-( 1 -)F'1!) (2.5)
0 ui>O
Ax
Aj A= x (2.6)
1 ui
The results in Figure 2.1(a) and 2.2(a) with different
CFL numbers demonstrates that a non-oscillating shock is
captured, but the shock is smeared due to the high
dissipativity of this scheme.
2.2 Flux-Corrected Transport (FCT) Scheme
McDonald(1984) made an improvement to the first-order
upwind scheme by introducing a flux-correction term which
results in a second order FCT scheme. His method is
separated into two parts, the first-order scheme and the
second-order scheme. The combination of these two parts
1.5
EXACT SOL.
1.0 O O O Oe9
S0 UPHINO
0.5
0.0 --~-e-o0-9-O-e-e-
-0.50.0 0.1 o.2 0.3 0o. 0.5 0.6 0.7-O.1 0.9 1.0
x
(a.)
I.5
1.0 i ----,--,--C-,.----- -,-----------
EXACT SOL.
SFCI
0.5
0.0 C -e- --e-
.0 .-. -.2 -.3 .4 o.- 0.6 0. 0.8 0.9 1.0
x
(b)
1.5 1.5
EXACI SOL. _EXnCT SOL.
1.0 --e---e-e--e-e -e-o-- -e- -e-e-e-eO-e--- e
E) LAX-UENOROFF IVO
.S 0.5 -
0.0 -e-e -e- -e- e- o.o -e-o-e-o- oe--e-
-0.5o.0 0. o2 0.3 O. 0.5 0.16 0.7 0.8 0.9 1.0 -0.50.0 0.1 0.2 0.3 0.1 o.S 06 0.7 0.8 0.9 1.0
X X
(c) (cl)
Figiu re 2.1 IThe complariison of numerical schemes at CIF'l=0.8
-1
0L
0a
'pQ
i
a, a
n
S oP
)(
Uj 3
SI)
C )
0
1
3
0
0,
0
n,
0>x .
0
0,
'I
0
In 0 1 0 41
0, 0 0, a
a
0 i
0
411 0 In 0 In
- ~ 0
1)
-J
V
U-
x c
u, Z
lo
0
0
0
0
0
5
Q
makes a hybrid form by the method of flux-correction
(Zalesak, 1979).
2.2.1 The First-Order Scheme
The numerical difference equation is written as
(2.7)
0 is the numerical flux given by
At F
Ax -
wi+1/2>0
(2.8)
LFX i+l
where w is a variable which is used to determine the upwind
direction. It can be defined as
wi+l/2=-(Fi+l-Fi)(ui+i-ui) (2.9)
2.2.2 The Second-Order Scheme
From the relation
F ,,au
t -^t
a finite difference equation can be written as
U =i -i .+l/2+? i--1/2
(2.10)
Oi+1/2=
Fn+ l/=F"- i (2.11)
The one-sided second-order differences in flux conservative
form is defined as
31"n+l/2 1n+1/2
si "2 i-1 W+1/2 0
Ai+t/= (2.12)
3 n+1/2 1- n+1/2
Fi+l i+2 Wi+1/2<0
2.2.3 The Combination-Hybrid Scheme
u,=ut--i+1/2+0i-1/2 (2.13)
Where, uj represents a first-order update. The second-order
flux-correction is difined as
64i +l/2==
Because of non-monotonicity of this flux-correction, a
filter, i.e., a flux-limiter must be introduced (Boris and
Book, 1973).
645i+1/2=-sgn(6i+1/2)max O, min(16i+1/21 sgn(6 ,i+i/2)(ui+2-u +I),
sgn(6i+,/2)(u$--u;_I)] (2.15)
Then the hybrid scheme is completed as follows
u2+'1= ui-6 4i+12+6i-1/2' (2.16)
The numerical results are shown in Figure 2.1(b) and
Figure 2.2(b). As can be seen from the figures, the shock
is sharper than the first order results and it is also free
of oscillation. However, this scheme requires much more CPU
time in each time step. The extension of this scheme to a
system of hyperbolic equations is quite complicated because
the direction of the characteristic velocity of the system
has to be determined at each grid point.
2.3 Central-Difference Scheme
In comparison with the upwind schemes, the central-
difference schemes require less CPU time and are easier to
program. However, classical second order central difference
schemes are not able to give smooth solutions in the
neighbhood of discontinuities. Generally all central
difference schemes can be expressed as the sum of a flux
term and a high order diffusion term.
uj ,=u-1 ( F --1 (2.17)
where F,4i/3=1(F,+ +F,), F,_1/2=.(Fi,_+Fi). Fi+/2-Fi-1/2 is the flux
term. Different ways to construct the diffusion term D will
result in different schemes.
If Di is expressed as
=i=( A-)2[(Ai+/2(F -F1)-Ai-_/2(F -Fit)]
AX 1+ 9 r I-
(.2.18)
Ai+l/2=-- 1(ui+ui+1)
Ai -/2-=-(Ui+ui-t)
the well-known Lax-Wendroff scheme (Lax and Wendroff, 1960)
is obtained. It can be seen from Figure 2.1(c) and 2.2(c)
that oscillations in the solutions are presented obviously.
2.4 Total Variation Diminishing (TVD) Scheme
The high resolution oscillation-free TVD scheme (Yee
and Harten, 1985) is constructed. If diffusion term I( works
as a flux-correction term. As treated in FCT scheme, all
oscillations are "limited" by a function called "limiter".
The TVD diffusion term is given by
i -= i+1/2-- i-1/2 (2.19)
The diffusion function 0 can be written as
1i+1/2=---b(ai+l/2)[Aui+/2 -Qi-1/2] (2.20)
where 0 is a correction term to prevent zero value of a.
( a)= 6=10-4 (2.21)
a'+62 <6
26
a=a=u (2.22)
Q is the limiter function
Q+1/2=minmod(Aui_l/2 ,I i+1/2 Ui+3/2) (2.23)
ui-1/2=ui--ui-1 (2.24)
The minmod function is defined as
minmod(x, y)=sgn(x)max{ O, min[ljx, ysgn(x)]} (2.25)
As shown in Figure 2.1(d) and 2.2(d), TVD scheme gives
a high resolution discontinuity. There is no oscillation in
the solution. TVD scheme also needs more CPU time than
upwind and Lax-Wendroff schemes do, but it is easier to
program than FCT scheme.
2.5 Quadratic Upstream Interpolation for Convective Kinematics
with Estimated Streaming Terms (QUICKEST) Method
The QUICKEST method is a third order scheme based on
the quadratic upstream interpolation (Leonard, 1979). This
method involves three-point upstream-weighted quadratic
polynomial interpolation for estimating control-volume face
values and gradients.
The finite difference equation can be written as
un+ n=+C,1U-Cru
U i = U7+! ut-- rtar
(2.26)
where C=uAt is the local
Ax
and r denote i-1 and i+1
2) 2'
Courant number. The subscripts I
i-1 i i+1
Figure 2.3 QUICKEST method
u;= (uC+u)- (1-C)CURV- ~Cr1GRAD
u;2CC 1 2
(2.27)
The curvature term, .which must be upstream-weighted
for stability, is defined as
uR-2uC+uL
CURV={
u,>O
(2.28)
uRR 2uR+UC Ur
where subscripts RR, R, C, L and LL denote i+2, i+1, i, i-1 and
i-2.
The gradient term is defined as
GRAD=ui+-u, (2.29)
u" can be found with corresponding formulas (shifted
one index to the left).
Leonard(1979) gave a result for linear model equation
as shown in Figure 2.4(a). For nonlinear equation, when CFL
is larger than 0.4, 'the numerical results present serious
oscillation. Even when the CFL reduces to 0.2, the result
shown in Figure 2.4(b) is still not good as other schemes
discussed above.
2.6 Time Marching Method
An Euler forward explicit scheme is used for time
marching. Even though implicit scheme can guarantee large
time steps due to unconditional stability, it needs much
more CPU time in each time step. .On the other hand, the
explicit method is more accurate because of the smaller
truncation error O(At). Another reason to chose explicit
method is that the explicit method is much easier to be
implemented than the implicit method for any kind of
computational domain.
1.5
1.0
S0.5
0.0
'0. 0o 0 0.1 0.2 0.3 0.11 0.s 0.16 0.7 0.8 0.9 1.0
x
(a) Linear equation
0
0
0 EXACT SOL.
00 0 0QUICKESI
-
-0.So 0.1 0.2 0.3 0. 9 0 .6 o.l? 0.L--. 1.0-
(b) Non- iin;r t (eq(iiLt. i ill
Figure 2.4 Ilesul-ts of QUICKEST meitl.hod
S0.5
S EXACT SOL.
E 0 OUICKESI
Eo E 6 9 -
0
_ iC~___;i_; ~ ___
26
2.7 Summary
From the comparison of these five schemes shown in
Figure 2.1, Figure 2.2 and Figure 2.4, only the solutions
solved by FCT and TVD have the sharper gradient variations.
However, the FCT method is more complicated than TVD
method. The programing of FCT method is also more difficult
than that of TVD method. Therefore, TVD method is chosen as
the numerical scheme for advection terms in the current
model .
CHAPTER 3
NUMERICAL MODEL
3.1 Governing Equations
The two-dimensional dynamics in the vertical plane
across the density front may be mathematically described in
space and time by means of a set of five equations together
with appropriate boundary and initial conditions. Four of
the equations represent the conservation of volume,
momentum and salt. The fifth relates the density of water
to salinity.
In Cartesian coordinates,
Conservation of mass:
auaw( )
=0 (3.1.1)
dx Oz
Conservation of momentum:
au a(uu) 8(uw)
9t Ox Oz
P= P (A u,+ (A, ) (3.1.2)
Poax 57 z aZtz OX OX
Suva vw (AA- )+ (A v) (3.1.3)
at _x Oz Oz Oz x dx
(3.1.4)
Conservation of salt:
as O(us) +(ws)
-t x Oz
a as 8 O Os
Equation of state
P=po++3s
(3.1.5)
(3.1.6)
where
x,z Subscripts denoting partial derivatives in
the longitudinal and vertical directions
u, w Longitudinal and vertical velocity components
v Latitude velocity component. It exists only when
the Corioli effect is considered
p Pressure
p Density
s Salinity
f Coriolis parameter
Ah,AV Horizontal and vertical eddy viscosity
Kh ,KV Horizontal and vertical eddy diffusivity
0 An empirical constant : 7.8 x 10-4
In equations (3.1.1) (3.1.6), the pressure is
Op
7- P g
assumed to be hydrostatic and the flow is assumed to be
incompressible and isothermal. Also, the Boussinesq
approximation is used.
The boundary conditions are
(1) At the sea surface there is no surface stress.
(3.1.7)
1K'Osso
C)
(2) At the ocean bottom there is bottom friction, and no
flow coming through the bottom
(3.1.8)
w=O
where A is a linear drag coefficient.
(3) At the coast, the normal flow and flux vanishes
u=O
(2.1.9)
all av)=
all az
L9 Z
AU,_x V=O
K.J =0
i8x
(4) At the open ocean end, the ambient condition are
specified as
7r=0
A,(u, v, )=O (3.1.10)
(''x ax Dx
K(as)=O
where 7 = Water surface level
=0 mean water level
>0 above mean water level
3.2 Numerical Model
A longitudinal and vertical two-dimensional model for
nonrotating water model is developed in this section.
Equations (3.1.1) through (3.1.5) are solved with a
two-time-level and central-difference scheme in a staggered
grid. In order to reduce the numerical oscillation and to
capture the sharp density variation in front zone, the high
resolution TVD scheme and adaptive grid technique are
employed in the horizontal direction.
3.2.1 Computational Domain
Figure 3.1 shows the computational domain and the
staggered grid. Computation begins from the bottom of the
water body to the surface for the continuity equation. The
momentum equations are then computed from the surface to
the bottom. The i-index is counted from the closed
boundary(left) to the open boundary(right). The k-index is
counted from the bottom to the surface.
3.2.2 Discretization of Equations
The first step in this numerical model is to divide
the water body into horizontal layers of thickness Az.
Equation (3.1.1) (3.1.5) are then applied to each layer
and integrated vertically across the layer thickness
resulting in a set of hybrid differential-finite-difference
equations (Cerco, 1982).
After integration, the continuity equation (3.1.1) is
rewriten on surface layer as
SW 8(uh) (3.2.1)
t = x .
Below the free surface, the equation (3.1.1) becomes
= (uh) (3.2.2)
wt = W x
2 3 4 5 '6
Figure 3.1 Computation domain and grid system
where
u = Vertically averaged value of velocity in each layer
h = Layer thickness
Az+1 surface layer
h=
Az otherwise
t,b = Subscripts denoting parameter evaluated at the
top or bottom of a layer
Rewrite Eq. (3.1.2) and Eq. (3.1.3) in the conservation
form
(3.2.3)
at +Fax Oz
Q^t 9x az
where
- uuh
uvh]
(3.2.4)
(3.2.5)
vh -u,w, +ubwb -P A +( h)+rT,-r
-fuh-vw, +ubWb+ 0(Ah)+ry,-ry,
Ox Ox
v = Vertically averaged value of velocity
in each layer
V= Iuhl
vhJ
From Eq. (3.1.5)
8(sh) O(u h)
-+ -+SteW-Sbbw
Ot Ox
=-(K(h )+(K) )-(K0 )b
where
s = Vertical averaged value of salinity in each layer
P = Integrated pressure over the layer thickness
r7,ry= Shear stresses
r,,,r,-= (Avu A ),
(A"z, A', ),' otherwise
T7,b 'yb -I
(-Au, -Av)b
estuary bottom
3.3 Finite-Difference Equation and Procedure
3.3.1 Find 1?+' by Solving Continuity Equation
Computation begins from the bottom and ends at the
surface.
From equation (3.2.2) and equation (3.2.1)
z ( 3 3 1)Z( U
,j Ax Zj-l il- '.f
(3.2.6)
(3.3.1)
At the bottom. w'=0. i.e. w'1=0
(3.3.2)
Sn _(u ~2, .n u -1
i.2-- =_.(U i u -" i-1.2)
At the surface, wt= Thus,
h Tl +Atwu
(3.3.3)
3.3.2 Solve u1' and v+1
The finite-difference equation for equation (3.2.3) is
t' =?+t(FLUX(W +R) (T?) (3.3.4)
where FLUX term can be separated into two terms
FLUX()v') = FLUXC(W )+FLUXT(W )
(3.3.5)
where FLUXC(W ) are the net fluxes solved by central
difference schemes and FLUXT(W ) are the TVD type diffusion
fluxes.
QF.-- F i+i,-F i-i/2
FLUXC(Wi)=---
(3.3.6)
F.+F,,
where F i+1/2--
F i_- +F
F i-l/2
In the current development, the concept of TVD scheme
is only defined for scalar conservation laws or constant
coefficient hyperbolic systems. It is impossible to fully
extend the theory of nolinear scalar TVD scheme to the
nonlinear systems. What we can do is to extend the scalar
TVD scheme to the system so that the resulting scheme is
TVD for the "locally frozen" constant coefficient system.
(Yee, 1987; Li and Kroll, 1989) For our two-dimensional
problems, this procedure can be greatly simplified.
The TVD type diffusion term FLUXT(Wi) is expressed as
FLUXT(W,)= (Ti;1/2Ji+l/*>-Ti-,/.9-1/2 (3.3.7)
where T is the right-eigenvector matrix of the Jacobian.
From equation (3.2.4)
Suu u\. (3.3.8)
F LuvhJ =iuV14
For a the "frozen" constant coefficient system, the
Jacobian matrix can be written as
A u (3.3.9)
Obviously, the two eigenvalues of A are
At=u A,=u
and the corresponding right-eigenvector matrix and its
invers are
T = 1 T- 1 0
T 0 1 0 1
(3.3.10)
In this case, the TVD FLUXT can be reduced to
(3.3.11)
The diffusion function 0 can be expressed as
(3.3.12)
Si+1/2=- i+1/(Ali/2 i+/2-0 i+1/2)
AW j+i/2 =T/12( i+i -WW)=' i+-W i
(3.3.13)
I'i+/2l
(3.3.14)
i+1/2=
i+1/21>10-4
lI.+l/21 < 10-4
where
a2i+1/2+6
26
FLUXT(W i)=. .Zi+/2-0 '-1/2-
6 is a small positive value.
ai+1/2= 2 -" U+1/2
The limiter function Q can be written as follows
Q i+,/2=minmnod(AW l;_1/, AW i+1/2' i+3/2)
(3.3.15)
In these calculations. no matrix is involved, so the
computation of FLUXT is much faster.
R(W1) consists of many terms. Writing equation (3.2.5)
again
1 0P 8 8ub
fvh-u,w,+ubwb Po- x + (A+ -')+r, ,- rb
-fuh-vw,+uIbwb+') (Ah Vh)+ ry-Tyb
O N 1N
where
t= .(Ui,j +ui,j+1)
vt= ~(ij+i,j +l)
Wt= (W i,j+wi+l,j)
U =. (uij-. +ui, )
2+ U
Vb= 1("ij- Vij)
b = (WJi,j-l+ iW+tj-t)
constants, the stress terms can be
(3.2.5)
(3.3.16)
Assuming Ah and Av are
39
written as
r = I U ij+= Ui.j
9v A =v vi.j+ --vi j
rt --A = AZ
r2 = '^ z~' A A' Az-
(3.3.17)
rb = Av- lb Avuj zj-
baz Az
Vi -V i.j-1
T yb -A b -= Av .Z- 1
j; b \z
Hence
S_ 9uh A ui+ i+/2 -2ujhi-/2.j +ui-l,jhi-3/2,
'kh-'-A )=Ah A2
xA ax Ax2
(3.3.18)
The pressure grident -lP is calculated from the
PoOx
surface to the bottom. At surface layer, from hydrostatic
equation (3.1.4)
(3.3.19)
p(z)=pg(7- z)
ap= 2- g- z)
ax 8r9x
Ox jax JI-ax
(3.3.20)
(3.3.21)
By Liebnitz's rule,
api = ,0 20P;
=P9g(i+Az) +g(Azd + xx
(3.3.22)
=pg( + z) I+ i- g(i+ z)2Pi+ j--Pi-l.J
=pg(?+Az ) "2Ax 4 +' Ax
At the subsurfaces,
S P IP .P Aza
ij.i--- ij4+1 gAz i,+t g (Pi+,j-Pi-1,i)
j=JM-1, ..... 1
(3.3.23)
3.3.3 Solve S'+1
Once the velocity components u, v are obtained, the
salinity equations can be solved.
From equation (3.2.6). the finite-difference equation
for salinity equation with two flux terms can be written as
n+ =( hi FLUXC(sh)+FLUXT(sh)+RHS(sh)3.3.4)
t,j h,,,
where h=Az+q for surface, and h=Az for subsurface.
Again, FLUXC(sh) is the fluxe solved by central
difference schemes. FLUXT(sh) is the TVD type diffusion
flux.
FLUXC(sh)=-_ x- sWt-SBwb)
1 (sh)i+l.,+(sh)i.j (sh)i.j+(sh)i-, i ,
= -A(x -2-u -2 u2-1'
(sw,-sbwb) (3.3.25)
where
st = (s ij l+s;
wl =wi,j
s=.-(si,j +s
b i j-i j
W"b= ;i~j--1
Since the salinity equation is solved independently,
the FLUXT(sh) can be calculate in the same way as solved in
scale equation discussed in Chapter 2.
FLUXT(sh)= (.i+1/2 _i-/2)
(3.3.26)
where 0 can be solved by substituting A(sh) into Au in
equation (2.22).
(3.3.27)
and the limiter function can be given as
(3.3.28)
The diffusion term of salinity equation is solved with
the assumptions of constant K,, and K1,.
Qi+,/2=minmod(A(sh)i-1/2., A(sh)i+1/2, A(sh) i+3/2)
A(sh)i+/2=. ((sh)i+(sh)i+,)
0H 1h sh= sh
RHS(sh) =Kh -+v
ox" 0z2
S2 =K 1 ((sh)i+ -2(sh)i,+(sh)i-1j)
(h -=K ((sh)i.+-,-2(sh)i,j+(sh);ij-
3.4 Stability
Equations (3.2.3) and (3.2.4) are solved explicitly in
the longitudinal direction. The computations have to be
subjected to several stability and convergence constraints.
Cerco(1982) discussed the following inequalities.
(Ah,Kh)>
(Ah, ,)h
(A uK2a)
(AAKh) 2
(3.4.1)
(3.4.2)
(3.4.3)
Inequality (3.4.1) is restrictive to the selection of
Ah and Kh. uAx/kh is well known as "cell Reynolds number."
If the "cell Reynolds number" is equal to 2, the equations
(3.4.2) and (3.4.3) satisfy the stability criterion CFL
where
(3.3.28)
(3.3.29)
(3.3.30)
If this number is greater than 2. equation (3.4.2) will not
satisfy CFL condition. In this case, the computation
remains stable but oscillations develop in the solutions.-'
Inequality (3.4.2) arises from the explicit treatment of
the longitudinal dispersion term and is more restrictive
than inequality (3.4.3) which is due to the explicit
treatment of the advection term. If At satisfies the
inequality (3.4.2), the numerical calculations are stable
and convergent. Inequality (3.4.3) is the dead limitation.
The violation of this constraint results in ultimate
instabilities in numerical computation.
Finally, the most restrictive stability criterion
associated with the wave propagation on the free surface
must be satisfied.
At < Ax (3.4.4)
~-gh
Actually, this is the dominant stability criterion to
limit the time step At.
CHAPTER 4
ADAPTIVE GRID TECHNIQUE
Adaptive grid technique is a method by which the mesh
points can be self-adjusted as the numerical solution
evolves in time such that grid resolution increases in
regions with large gradients of solutions but decraeases in
regions with little or no gradients. This adjustment can
reduce the truncation error and improve the numerical
accuracy in regions with sharp gradients such as an
estuarine or coastal front. In addition, grid points will
be reduced in regions with little or no gradients, thus
increasing the numerical efficiency.
Due to its potential for improving the accuracy and
efficiency of the numerical solutions, adaptive grid method
has have rapidly emerged as an important tool in numerical
simulations of flow problem. With adaptive grid, each grid
point acts like an internal observer whose purpose is to
monitor the solution. The grid points distribution over the
computational domain is correspondingly adjusted
dynamically to concentrate grid points in regions of high
solution variability. The improved accuracy is achieved by
the redistribution of the grid points to increase the
45
number of grid points in the region where the solution has
large gradients. This reduces the truncation error in the
numerical solution. The ability of the adaptive grid to
redistribute the grid points into the regions requiring
finer grid resolution on an as needed basis reduces the
number of grid points required to compute the solution.
Thus requiring less storage and computational time for
performing a numerical simulation.
In the front zone, sharp density gradient often exists
over a spatial scale of tens of meters. Since the frontal
dynamics are affected by the large-scale estuarine and
coastal circulations, it is essential to model both the
large-scale circulation and the front-scale dynamics. If a
non-adaptive grid is used. very fine grid resolution in the
entire computational domain is necessary to obtain high
accuracy solution, thus requiring large computer storage
and huge CPU time. In thus cases, the adaptive grid
technique can be employed to refine the grid locally while
keeping the total grid point number unchanged, to obtain
accurate solutions at the front zone. This should lead to
savings in computer storage and CPU time.
4.1 Grid Adaptive Algorithm
The adaptive grid scheme developed here is developed
based on the idea of "minimal moment" of the moment-
balancing principle according to Connett et al.(1988). It
belongs to the category of equidistribution adaptive grid
method (Thompson et al., 1985).
Let a segment L in a one-demensional grid have a mass
V. then the density P of this segment (not the density of a
fluid) is W/L (see Figure 4.1). The movement of a grid
point, i.e. the center point in the segment L is determined
by the two adjacent grid points and the densities assigned
to the two segments (Ax, and Ax,l,). The algorithm is set up
to move a single grid point at one time. Once the local
algorithm is defined, the entire grid is traversed. The
traversals are repeated until a satisfactory adaption of
overall grid is achieved. The local algorithm is discussed
in the following.
W=PL
Pi-1/2 Pi+1/2
xi_- Ax, x, Axi+i xi+l
I L
Figure 4.1 Initial Mass Distribution
If the density P is uniform over the entire segment L,
the mass center is on the segment center. Then the moment
of the segment L is minimized.
When the density on either side of x, i.e. -x. and
-x.. is changed. the mass center should be chanced. If
P,. on -:<. is larger than P,./\. the mass center should be
moved to the left to a ne' position x:,) ( Figure 4.2).
N-
new mass center (x.) new
old mass center (x,),od
'-1/2 I Wi+l/2
I:,
Figure 4.2 Moment Balance Principle
As shown in Figure 4.2, w.i;/ (wi-1/2= Pil/ 2xi)
w i/ (wi+1/2-=Pi+/12) are weight functions on Ax, and -Ax+i .
and
Let us define
10o=11+1,= (Ax,+Axxi)
(4.1.1)
According to moment-balance principle,
wi-l/2 11= "i+ /2 12-
(4.1.2)
SWil!/2 I ',-k/'
Wi+ /2
-Wi+ (' ." .-i) (4. 1- 3)
Thus.
S- -) (4.1.4)
-1 2 w-/*+
Therefore. the mass center will be shifted from (xi)o, to
(x) new
(xi) neUw=x .i-I+1- + l
=(Xi)o.1d -A'NSi i+2 :>:>
=(x,)old xi + li (4.1.5)
The adaption will be advanced from i=l to i=Imar. The
adjustment of each point is based on the previous point
positions. Then repeating the same procedure to find 14 and.
(x,)ne. until I|(Ax,)new-(.Ax.i)o dIe (here e is a small value).
4.2 Weight Function
The main idea of the adaptive grid method is to have
the grid points move as the physical solution developed. A
sensor S has to be defined on the grid points, and the grid
is driven to respond to the magnitude of the gradients of
S, so as to obtain small cells in the presence of the
larger gradients of S and large cells in the presence of
the small gradients. Usually, there are many quantities to
be solved in a physical problem. The quantity which has
discontinuty is chosen as the sensor S(r). For the front
problem discussed in this paper. salinity is chosen as the
sensor S.
As described in section 4.1, the weight function is
the product of the density and the length of the interval
Ax,.
wi-1/2=Pi,-1/2- xi (4.2.1)
Where, the density function P is some kind of combination
of the first derivative of the sensor S(z), S, and the second
derivative of the sensor S(x). Si, so that P can present the
large variation of the salinity in solutions. The selection
of P for the adaptive grid method is somewhat arbitrary and
subjective. Two types of P are used in this study:
P=(1 +eYts, js.l) (4.2.2)
p l +aIlS I+3|SI (4.2.3)
1+|ISrI+ISI
The arbitrary constants a and 3 are the adaptation
parameters. When a increases, the sensitivity of P to Szz
will increase. When 3 increases, the sensitivity of P. to S,
will increase. There is no theoretical way to chose a and
3. The determination of their values are based on numerical
experiments.
Another influence on the sensitivity of P is the
formular of P selected. If P is too sensitive with the
variation of ISzzj and ISZi like the exponential form Eq.
(4.2.2). the adaption may not converge and may sometimes be
unstable. Usually, the formula of P is chosen as polynomial
type (Thompson et al., 1985. Shyy, 1987), Eq. (4.2.3), so
that it is not too sensitive to select the adaption
parameters a and 0.
Satisfactory adaption generally can not be achieved in
one step. During the period of repeated adaptions, new sets
of grid system continue to be built up, by using new weight
functions which are interpolated from the old values.
4.3 Steady-State vs. Transient Problems
Equidistribution adaptive method is very effective for
solving steady state problem. For unsteady case, the grid
system always follow the solution. Unfortunately, since
the adaption is performed based on knowledge of the
previous solutions, the region where large gradient of the
solution occurs moves forward with the development of time.
Thus, the refined grid region always lags behind the new
large-gradient region. Therefore the adaptive effects are
not obvious for this unsteady flow problems.
If the adaption can capture the movement of the large
gradient or the refined grid region can catch up with
large-gradient region, this grid system will be able to
simulate the large variation of the solution solved in next
time step. Thus, an effective adaption can be achieved. To
catch the movement of large gradient, the moving grid
method (Eiseman and Bockelie. 1989) can be used. However,
the moving grid method is more complicated and requires
more CPU time than the equidistribution method. In this
work, an alternative method. which is basically a modified
equidistribution method. is developed to refine the grid
points in large-gradients region.
In general, weight function includes the influence of
the first derivative and the second derivative of the
sensor S. However, at the downstream region of the front,
there is no influence from the upstream region, hence the
two derivatives are very small or even zero, i.e. the weight
functions are very small. In this case, the adaptive grid
cells are always coarse and the fine grid cells lag behind
the sharp "head" zone of the front. Therefore, there must
be a weight function correction to transfer the influence
of the two derivative in front "head" zone forward to the
downstream ambient area. The following weight function at
point i includes the influence of weight functions at two
previous points i-i and i-2.
w,= (Aw,_2+Bwi_+Cwi) (4.3.1)
where A, B and C are arbitrarily chosen constants.
4.4 Temporal Coupling
Between every adaption, interpolations have to be
made. If adaption advances through the entire computation
period, the solution should be satisfactory. However,
interpolation error will always exist whenever adaption is
employed. In order to reduce the interpolation error,
adaption will not be performed at every time step. Starting
from a particular moment, the solutions obtained on a fixed
grid are used as previous knowledge and the adaption
begins. Over a certain time interval, adaption is performed
every time step. At every step, a new grid is developed
based on known solutions. The known solutions are then
interpolated onto the new grid for marching into the next
time step. By successively repeating this procedure, the
forward matching front is captured on the adaptive grid
system until the desired final time instant.
For front problems, the layer in which the salinity
has a typical horizontal gradient distribution representing
33
a front, is chosen as the reference layer. The sensitive
variable S is then taken from this layer. After adaption.
all variables are interpolated onto the newly established
grid system.
The grid cells will become more and more clustered as
the adaption period increases in time. In order to avoid
the development of infinitesimally small grid cells, a
minimum grid size is chosen to be approximately 10% of the
initial fixed srid -size. and 3 can also affect the
smoothness of the grid cells.
4.5 Verification of The Adaptive Grid System
The one-dimensional Burgers equation is used to verify
the adaptive grid method in capturing the discontinuity.
Three cases are provided here to present the steady state
and transient solutions. The transient case has the same
properties as the front problem.
Consider the Burgers equation.
au 4.ILu (4.5 .1)
The initial conditions are
Case (1)
1 x=0
u = linear O
-1 x=1
Case (2) u = cos[(2k+1),rxl O
Case (3) 1 x=O
u = ^(4.5.4)
0 x>O
In case (1) and case (2). a steady-state solution can
be obtained, whereas case (3) is essentially a transient
case.
Figure 4.3 through 4.6 present the comparison between
the solutions on adaptive grid and fixed grid. For cases
(1) and (2), the weight functions are exponential types,
and solutions are given in Figures 4.3 and 4.4.
respectively. The discontinuity captured in the adaptive
grid with 20 points is as sharp as that in the fixed grid
with 80 points (See Figure 4.6(a)). Figure 4.3 through 4.6
also show the evolution of the solution from the transient
to a steady state. The grid points move with the variation
of the solution until the final steady state is arrived.
Figure 4.5 presents the marching discontinuity. A
polynomial type weight function is employed in this
adaptive grid. The discontinuity captured in adaptive grid
with 20 points is as sharper as that in fixed grid with 40
points shown in Figure 4.6(b).
2.00
1.50
1.00
0.50
S0.00
-0.50
-1.00
-1.50
-2.00
2.00
1.50
1.00
0.50
0.00
-0.50
-1.00
-1.50
-2.00
2.00
1.50
1.00
0.50
0.00
-0.50
-1.00
-1.50
-2.00
-1.00
-1.50
-2.00
II ME 0.1193
-I__fl.-LLULLLUU LLJ Jlll__-
x
I I L I). IJ iUJl
LLLJ A .AI-LIII LL-Lt L- L
l'i gi;re -1.3 (Cas (1) ,in iLda.pt i ve yi icl I itinl.s 20()
TIE- 0.0000
I I 11111I1I1 II 11
2.00
1.50
1.00
0.50
, 0.00
-0.50
TIHE- 0.25011
I I I I I I IIIII I I I I I
I IHE- 0. 19'16
1.001 --
0.50
0.00
-0.50
-1.00
-1.50
-2.00
2.00
1.50
1.00
0.50
30.00
-0.50
-1.00
-1.50
-2.00
LILIIJL__LL__ LLLLI__.LIJL-_
x
I*'igis' p -i.1 Ca~s. c (2) oni adiaplt-. iv<: grl Id"| iints.. "-0
2.00
1.50
1.00
0.50
S0.00
-0.50
-1.00
-1.50
-2.00
T IHE-0. 11167
TIME-0. 7631
I I I 1 I I I
-
- IJ_ _J _J _L J .. _
2.00
1.50
1.00
0.50
: 0.00
-0.50
-I .n
-1.50
-z.oo LlIIIIII I I
2.00,
1.50
1.00
0.50
: 0.00
-0.50
-1.00
-1.50
-2.00
I IHE- 0.372'1
--LILLL11 IIIILLII I
I'igi're '1. 5 Case .(3) oil adaptLive' g i(I point.H 20
T IME-0.0000
I I I I I I I I I I I I I I I I I
2.00
1.50
1.00
0.50
S0.00
-0.50
-1.00
-1.50
-2.00
2.00
1.50
1.00
0.50
0.00
-0.50
-1.00
-1.50
-2.00
IIHE-0. 1101
F-
I I I I
iTi E- 0.2216
1 I li l i t IIllIll I I I I
,
-------
I I I
tx 21.
-I I 1 I LLL___LIL__L_L LLI
IX 'iI .
(aL) Case (1) on uniform grid 20 land 80
(b) Ca u (3:) on uni orm grid 20 "ni( -10
'a'lT u coijv e riid l ion ( 1 ii i in i ,gri
adaptive grid
2.00
1.50
1.00
0.50
o 0.00
-0.50
-1.00
-1.50
-2.00
2.00
1.50
1.00
0.50
0.00
-0.50
-1.00
-1.50
-2.00
2.00
1.50
1.00
0.50
S0.00
-0.50
-1.00
-1.50
-2.00
2.00
1.50
1.00
0.50
0.00
-0.50
-1.00
-1.50
NX- 1.
1I1111111111 II II 11 111111 II ill lll III IIIIIIlII1111 1 II Ill1
Figure I4.0(i
CHAPTER 5
NUMERICAL EXPERIMENTS
Front problems are very complicated. In general, there
are no theoretical solutions for the real physical
problems. Some theoretical solutions and laboratory
experiment results have been obtained by simplifying the
real problems. In this chapter, four kinds of physical
problems will be solved by the numerical model described
earlier in this thesis. Three of the problems contain
simplifying assumptions while the fouth one is more
realistic. Whenever possible, the solutions are compared
with theoretical solutions, laboratory experiment results
and field observations.
5.1 "Breaking Dam" Problems
5.1.1 Physical Problems
Two fluids with different densities can be adjusted
under gravity from initial unbalanced state to final
equilibrium. This has been an old but a very typical topic
discussed by ocean and atmosphere scientists since the
beginning of 1900's.
Water spreading out from a suddenly broken dam becomes a
plume with the water front at its boundary. Due to its
similar physical feature to the suddenly released fluid
flowing into the surrounding ambient fluid with a
vertically varying density structure, the "breaking dam"
problem can be regarded as a simplified front problem.
x=O
Figure 5.1 Initial State of Dam
Assumming the "'breaking dam" flow is one-dimensional
i.e. along the horizontal direction, and neglecting viscosity
and diffusion terms. equation (3.2.1) and (3.2.2) become
8uD+D =0
x Oax
a D uuD -fD= 0D
at +x =O
+vD1uvD+fuD=
at dx
(5.1.1)
(5.1.2)
(5.1.3)
where
62
x
x>O. t=O
If Coriolis effect
equations are reduced to
is neglected,
the governing
OD, u=uD_
Ot ax
uD+9 uuDl 8D2
Ot Ox 2 .O
For these simplified equations,
solutions are given in the Appendix.
the theoretical
5.1.2 Numerical Results
In the "Breaking Dam" problem. the finite difference
equations are solved with spatial increament of Ax=0.025.
The time steps are limited by the following condition:
(5.1.7)
At < CFL( Ax )
nazxu+cl
where c=IgD and the Courant number CFL=0.6.
Figure 5.2(a) presents the analytic solutions at 10
time instants. The numerical solutions. of equation (5.1.5)
(5.1.4)
(5.1.5)
(5.1.6)
and (5.1.6) are shown in Figure 5.2(b). As can be seen from
the figures, the numerical solutions are in excellent
agreement with the theoretical solutions.
Figure 5.2 shows that immediately after the release of
the initial discontinuity, a hydraulic jump forms. The wave
fronts move away from the initial discontinuity at the
constant speed. Due to the hydraulic balance and nonlinear
effects, the discontinuity propagates like a compression
wave towards lower pressure region. and a rarefaction wave
travels in the opposite direction into the high pressure
region.
Figures 5.3(a) and 5.3(b) demonstrate the effect of
Coriolis terms on the numerical solutions at 10 time
instants and 25 time instants respectively. Figures 5.4 and
5.5 corresponding to the Figures 5.2 and 5.3 present the
numerical solutions at 6 time instants and 3 time instants
respectively. As pointed out by Rossby(1938), the Coriolis
force simply balances the pressure gradient. Based on this
geostrophic balance, the gravitational adjustment process
leads ultimately to a steady state.
Rotation effects contribute an additional dispersive
term to the original one. i.e.
W2= J +k'c (5.1.8)
The wave with this dispersion relation is reffered to as
the "Poincare Wave" (Stoker. 1957). The wave front moving
away from the initial step is followed by a "wake" of waves
that trail behind because of the dispersion. The decay of
the rightward propagating jump is caused by a rarefaction
of Poincare waves behind the compression. The leftward
propagating jump forms also by the trailing wake of
Poincare waves which increase the pressure after the
passage uf the front. At t=6. the second jump appears
because of the nonlinear steepening of the Poincare wave.
From the solutions, it can be seen that there is
energy loss between the initial state and the final
geostrophic equilibrium. The wave fronts carry energy, so
for any finite region the energy is lost through the sides
by "radiation" of Poincare waves until the only energy left
is that to hold the final balance. Cahn(1945), Gill(1982)
and O'Donnell(1988) discussed the linear and nonlinear
rotation adjustment problems by theoretical and numerical
calculations.
5.2 Flow With Linear Density Distribution
5.2.1 Physical Problem
When a horizontal density gradient exists in a fluid,
vertical gravitational field vorticity is generated, and
109 ----- -__T9 .-.r-- ...
0 --- U ....... .....
S------ ------ -
3 U ------------ -
3 -- -- 3 -
0 | .---- _______ -- .--- ---- -- -- -L--- --
-12 -10 -8 -6 -4 -2 0 2 i 6 8 10 12 -12 -t0 8 -6 -4 -2 0 2 1
x x
(u) (l,)
Figure 5.2 "nreakinig dla" problems
(a) Ana.lytlic -olutionii at unit. int.urval o( I' t.iIiie
(1b) N ntiiurical rol utiioni c:o -r,<; pi)",, i ngu t,< (i.
' I I I I
S .
_-.6 I 10 12_ I_
6 0] 10 12
0-
a)
'I I- I S *I 5 V11(~ 1rgyr
-JN
-1 I -10 -8 -L -
-12 -10 -8 -6 -4 -2
~L_1_
0
2 1 6 10. I
x
(a)
20 ------
10 --- -- -
.I.
6 __--------'--
I -------- --- --
S _--__...._.. ... . .
0- -- --- -~___ -- -
o
CTJ
03
- 1 -- 1 --I --1- 1 J-I 1.._ __.l .... -. I._ l._.l 1
12 -12 -10 -8 -6 -4 -2 0 2 I 6 0 10 12
x
( ')
Figure 5.3 'The 1numerical solutions F'or 1.1t,- layt.: (lt:lt).i itl.l
rIo,.i.ta ion clffct;s in (a) 1.:z 10 i(I (I>) t.=2-
Figure 5.4
Nonratating "breaking dam" problem at t=l and t=6
Figure 5.5
With rotating effects, the gravitive adjustment
at t=1 and t=3
flow is induced no matter how small the density gradients.
As an example, an uniform horizontal density gradient
distribution producing fluid motion in a rectangular tank
is discussed.
X
Figure 5.6 Density Distribution
As shown in Figure 5.2-1, the density distribution
follows p=p(1--ax), where a=!p is a measure of the
horizontal density gradient and ~ is the value of the
largest density in the tank. The length and depth of the
tank are 360cm and 15cm, respectively. Since the flow is
two-dimensional, the width of the tank is eliminated from
the equation. Simpson and Linden(1988) presented their
theoretical and laboratory results under the same
condition.
5.2.2 Numerical Results
The numerical experiment is solved with spatial
increment Ax=10cm and Az=1.5cm. The time step is
At=0.05s.
Figure 5.7 is the numerical contour plots of salinity
obtained at 8.0s, 12.3s. 16.Ts, 20.9s and 24.9s after
initial state with =1.058 x 10-3g/cm-3. The vorticity
generated by the horizontal density produces a flow towards
the right in the lower half of the tank and towards the
left in the upper half. The vertical stratification
increases with the development of time. The isopycnals
remain straight throughout the motion. Corresponding to the
numerical experiment. Figure 5.3 shows a time sequence of
photographs of laboratory experiment provided by Simpson
and Linden(1988). As can be seen, the numerical results and
laboratory experiment have good agreement. In figure 5.3,
the angle of inclination 9 'of the isopycnals to the
vertical for theoretical, numerical and laboratory results
is plotted as a function of time t. In this log-log plot,
the time has been non-dimensionalized by (.ga)2 According
to Simpson and Linden's(1988) theoretical computation, the
initial angle of inclination is given by tanO=,. The
laboratory experiment is realized by suddenly rotating the
tank from vertical position to horizontal position so that
the motion begins with an initial 0. In the numerical
experiments, the numerical results are recorded from a
initial 0 which is the same value as that in laboratory
experiment. Two sets of numerical results are obtained with
a=6.5x10-scm-l, p=I .058x10-g/cm3 and a=15x10-5cm-, p=1 .0
r~ ----=--=---~==-L~=_~-===---,-_-
IFigiire 5.7 'Nimtecicxl t;imI-sec(tionce showing t.Iicisociias
Wit.ii f= 1.058x 103gl/cn 3, a= I5X 1 0-S CI -
atl. t.=8.0s" 12.3s, 16.7s, 20.9s and 2,1.9s.!)s
~1
r (
1;~~~--~----------------
~"=-~==~"=~-~===----------'
II. U I1Y*.)S ** &
Figure 5.8
Corresponding to Figure 5.7 a t.ime-scquc:
showing the motion of the isopyca: ls in
laboratory experiments from Simpson iULIl
Lini(Ien(1988). Only the right-hIn tl I( l" f1' tIl O i
tank is visible and the vertical I i nes are iLt.
10cm intervals
'Theory
=-6.5xl0-scrl-l
.- ---. r=15x10-cm-I
o o 0 I.ab 1
X x x Lab 2
-. 0 0
O
-
0
- -- 0
Sx
xx
-0.60 -0.20 0.20 0.60
log(' )
F igture 5 .
A log-log plot of the angle 0 of iicn:l ii t.iion (,1
the isopycnals to lhlie ver.tl i:ic l againsitl t.ill:
cl imensionless t =(ly/2) ti
2.0
1.6
1.2
0.8 -
0.4 -
0.0
1-.00
1.0(
g/cm3 are shown in the Figure 5.2-4. rhe numerical results
are identical to the laboratory results and theoretical
computation. In the numerical experiments, the top layer of
the tank is treated as a free surface whereas it is rigid
boundary in actual laboratory experiment. Therefore, some
error is brought into the numerical solutions. For rigid
boundary, the forms of governing equations will be changed.
Thus, there are must be the changes in the numerical model.
5.3 Sudden Release Problem
5.3.1 Physical Problem
When a fluid is suddenly released into the ambient
fluid which has a vertically varying density structure, a
gravity front is created. In the laboratory experiment, the
gravity front is formed by releasing fluid of a different
density from behind lock-gates. Since the volume of the
fluid behind the gate is usually much smaller than that of
the other fluid, the experiment approximates the release of
a finite volume of fluid into another fluid of infinite
volume. Rottman and Simpson(1983) presented their
laboratory results. In order to compare with the laboratory
results, four test cases are simulated. The first two cases
are in the same conditions as in the laboratory
experiments, which represent bottom plumes.
Figure 5.10 Laborat
Where xo=50 cm
ho=7 cm
Po=1 .048
ory Experiment Case 1
L=150 cm
H=49 cm
g/cm p1,= g/cm3
Figure 5.11 Laboratory Experiment Case 2
Where xo=50 cm L=150 cm
h0=7 cm H=7 cm
po=1.048 g/cm pl=l g/cm3
For realistic large-scale river and ocean flows, there
Case 1:
Case 2:
are not only bottom plumes, but also surface bouyant
plumes. Therefore. the next t-'o cases describe the Large
scale bottom and surface gravity currents.
Case 3:.
L
H !
lI
Figure 5.12
Large Scale Bottom Plume
Where x0=3 km L=40 km
ho=7 m H=10 m
Po=1.000546 g/crn
p,=l g/cm3
Case 4:
L
Pi H
XFi
Figure 5.13 Large Scale Surface Plume
where xo=3 kmn L=40 km
ho=T m H=10 m
Po=1 g/cm
p =1.000546 g/cmn3
5.3.2 Numerical Results
For the sudden release problem, the numerical
solutions are given in Figure 5.16 to Figure 5.19. From the
figures, it appears that the two-dimensional non-rotating
gravity current produced by an instantaneous release passes
through two distinct phases. There is an initial adjustment
phase during which the initial conditions are important,
and an eventual self-similar phase. In the initial phase
the "nose" velocity, i.e. the velocity of gravity-current
front, is about constant. In the self-similar phase, the
nose velocity decreases. The transition from the first to
the second phase occurs when a disturbance generated at the
end wall overtakes the front.
Figure 5.14 and Figure 5.15 present Rottman and
Simpson's(1983) laboratory experimental results, which
correspond to Figure 5.10 and 5.11. The numerical solutions
in Figures 5.16 and 5.17 show the same density front
formation as obtained in the laboratory experiments. The
bores, surges and hydraulic jumps exist during the front
7$
propagation. The comparison of computed and measured front
velocities are shown in Table 1.
For realistic large-scale coastal or river density
plume. the numerical solutions shown in Figure 5.18 and
5.19 also can give the same results about the front
formation.
The solution shown in Figures 5.18 and 5.19 demonstrate
that the hydraulic jumps occur in both bottom and surface
currents and the two phases of fronts are obvious. The
front velocities are given in Table 2 in comparison with
the laboratory experiment and theoretical solutions. The
theoretical solutions which are solved by characteristic
method are from Rottman and' Simpson(1983). The front
velocities follow the relation
u =F(I) (5.3.1)
Theoretically, the bottom front has larger velocity
than surface front due to the reduced acceleration of
gravity g' g. This can be found from table 2 which gives
both bottom and surface front velocities at the same
initial conditions.
Front Velocities
Table 2. Bottom
and
Surface Front velocity y
Wang(1984) investigated both surface and bottom front
problems by his numerical model. His results are shown in
Figure 5.20 and 5.21. In comparison with his results,
solutions obtained by using current model are shown in
Figure 5.22 and 5.23.
5.4 Water Discharge Problem
5.4.1 Physical Problem
The discharge water of a different density into
Case Theory Lab Numerical
1 11.79 12.24 12.5
2 8.16 5.S 6.9
Case Theory Lab Numerical
3 8.709 9.68 10.68
4 8.707 9.67 8.47
Table 1.
Bottom
Sf1r
I
figure 5. 14' Shado(wgralphs of a. vol ei of Isal 1t. wa ;.t:r w i li
Iho=7cm, II=49cm,, xo=50cm, g' =17cms- col li apsi i ntl
freshly wat.t'r at 2.2s, 3.8s anld 6.() a; t.(:r I.':l:;at.-i
The dlLtel v'ert.ical line indicat. v th.lie ipsition
oC t*lhe I ock galte, tlie uIe wal I his al. f;ir left., i;nld
-th th in v;:rt. ical linlies .ar at. I Ol ().n ii n4ti e val'i
(llotlrnani and Simpson (1983))
Figre 5. 15 Shadowgraphs of a voliae of sait. wtt r Wi Li*
. /11 x* .0 17c 'n .
'. .' *
fl I
Figure 5.15 Shadowgraphs of a volume of salt waltr wiitlh
ho/|i/=l ho==7cm, xo=50cm, g'=417cus-" col lnplsing
into fresh water at 4.7s, 7.7s anmd 10.7s a Tt cr
release (Ilottmlan arid Simpson (1983))
1 _/ ii ___i iii~iii~ __ ~_~__C_____~__I_____I _~ ____i
39.00
32.50
26.00
19.50
13.00
6.500
i'igure 5.16 NPumerical solutions correspol tlnding to lF'i gllic 5 .1 1
v I I Ir E 0
.*,,-- ------
L. .. ... .].. . . .
. .
2 .................
.I... .... .
....... .......
.. .. .
r-~i~i: ..P~r........ .... .............. ..
...... .... ..
;ii F'!I"."M--~-~1'1M02-121% M,'----
50.83
32.50
23.33
14.17
Figure 5.17 Numerical solutions corresponding to Figure 5.15
i 7-7 1 -T i E
---- --, ------- -
0 40 km
___________________ ________- ______ -- !t I t sii :
-10.m ... .
.. "2 ,-: -------,-,--
0.5540 0.5080 0. i620 0.4160
Figure 5.13 Numerical solutions of Case 3
at 20hrs. 40hrs, 60hrs and 80hrs
- - - - -
b7- .. .... .......
SI I t J~ I I 1 I I t i ii I t L I [ I i -
I- -
tI.
::-: 4
]1
I C [ I. It i~I i I. i I- I~ .i.. i L .i i i t. i. i. '-..
l_ i t Fl~~r LrE I F t I I I It I i
...... ............. ....... ......
.- : . : : : : . . .. . . . . . .... .. . . . . ....~-~~~:~. ~~
.......... ...... ...~:~- ~:
. . .
. . . .
0.1860
0.1520
0.1180
0.8400E-01
Figure 5.19 Numerical solutions of Case 4 at 20hrs,40ihvs, 60hrs
and 80hrs
8*
o 10 20 30
0o o 20 so
aI
to 20 10
Figure 5.20 Numerical Solution for Surface Plume
at 20hrs. 40hrs, 60hrs and 80hrs by Wang(1984)
0 t0 20 i0
.4
E1
.r..
0.
0
2
0.
E
.u~'
0.
0
2
a
0.
S
0
a.
a.
CL
2 0 30
a 10 20 30
Figure 5.21 Numerical Solution for Bottom Plume
at 20hrs. 40hrs, 60hrs and 80hrs by Wang(1984)
Ii
Figure 5.22 Numerical .Solution for Surface Plume at 20hrs,
40hrs, 60hrs and SOhrs Coresponding to 5.3-11
|
RSS
HIDE
