|
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
ACIKNOW LEDGEMENT
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 NI. 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 Mir. 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.
i ii
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.T 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 i7+lby Solving Continuity Equation 34 3.3.2 Solve u and v .. ..................... 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 .......... :** * * I 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
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 Colernan(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
3.5 W'ith rotating effects, the gravitive
adjustment at t=l and t=3 ..................... 68
I
5.6 Density Distribution ......................... 69
5.7 Numerical time-sequence showing the isopycnals "
with p=1.058x10-3g/cm3. a=15x10-5cm-t
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'=(Ige) '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, x0=50cm, g=47cms-2
collapsing into fresh water at 2.2s, 3.8s
and 6.Os after release ....................... 80
5.15 Shadowgraphs of a volume of salt water with
h0/H=1, ho=7cmn, x0=50cm, g'=47cms "
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
53.17 Numerical solutions corresponding
to Figure 5.15 .................. ............ $3
5.18 Numerical solutions of Case 3 at 20hrs, 40hrs,
60hrs and SOhrs ............................... $4
5.19 Numerical solutions of Case 4 at 20hrs. 40hrs.
60hrs and SOhrs ............................... 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 8Ohrs by Wang(1984) ............ $7
5.22 Numerical solution for surface plume at 20hrs.
40hrs. 60hrs and SO/wrs by current model .......... $8
5.23 Numerical Solution for Bottom Plume at 20hrs,
40hrs, 60hrs and 0Ohrs 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, 40hrs. 60hrs and 80hrs
(a) Velocity vector .......................... 94
(b) Salinity contour .......................... 95
5.27 Water discharge problem with Fr=0.057 for
underflow case at 4hrs. 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 80xlO 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
80xO10 at 20hrs. 40hrs. GOhis and SOhrs ......... 104
5.32 The salinity contour obtained on uniform grid
40x10 at 20hrs. 40hrs. 60hrs and SO/irs ......... 105
5.33 The salinity contour obtained on adaptive grid
40x10 at 20hrs, 40/hrs. 60/rs and SOhrs ......... 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
Shrs 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, 6ihrs and Shrs ............. 112
5.40 The salinity contour obtained on adaptive grid
25x10 at 2hrs, 4hrs, 6/hrs 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
Requ irements 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 seperated 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 M3s51, 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 Nonk(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, [ao et al.(1978) and Valentine and IKao(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 91 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 Ilermite interpolation functons and the density value p and its local gradient L 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 efficency. Later Wang(19S4) 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 threedimensional one to study gravity currents produced by instantanous 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
S
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 (Bergurs 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.
I.
.- I
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
S.
EBBING TIDE Gulf watr
Tac 0 Soo00c
Figure 1.3 South Pass effluent into ambient water masses
during flooding and ebbing tides from Wright and
Coleman(1974)
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 difficulities caused by the numerical treatment of advection terms. An ideal numerical scheme should have the following properties:
1. Accuracy: Second order accuracy is usually requ i red.
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 diffusioncorrected 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
14
employed to solve hyperbolic equations are disscused based on a model equation, i.e., a one-dimensional advectiondiffusion equation:
9OU=+F 02U (2.1)
where
u Dependent variables
F -1 u
2
S- Viscosity coefficient
t Time
x Spacial coordinate
For inviscid flow, p=0. The equation becomes
- =- -(2.2)
In this equation, the characteristic speed is -F=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=1
All numerical computations are completed on a grid of 20 points. The Courant number, CFL=u is chosen 0.8 and
0.2.
2.1 First-Order Upwind Scheme
The first-order upwind scheme is writen as
u '+=u(-A( OF',++(1 -20)F'-( 1 -0)F'-) (2.5)
0 u;>O
O A= At (2.6)
Ax
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
seperated into two parts. the first-order scheme and the second-order scheme. The combination of these two parts
1.5
EXACT SOL.
1.0 9
.~0 CC C CC O UPHINDO
0.5
0.0 -- ---o-9- - -a
-0.50.0 0. 1 0.2 0 .3 q 0. 5 0!6.O 0.18 0!9 10
x
(a)
1.5
EXACT SOL.
0 FCI
0.5
0.0 -C C 9 -- -I 0.S 0.6 0. 0.8 0.9 1.0
x
(t))
1.5 1.5
EXACT SOL. EXACT SOL.
1.0 --------- .0 -e-e---e- -e--e-0 e
0E LAX-NENDROFF 0 TVU
O.S 0.5
o.o -e--e-E-e- 0o.o -e---- e--e- 0.0 0.9 0.2 0.3 o 0.5 0.6 0.7 0.8 0.9 -0.50.0 0. 0.2 0.3 0.4 s 0.6 0.7 0.6 0.9 1.0
x x
(c) (ci)
Figi re 2.1 The c |)comparIson of numerical schemes at CoF'I,=0.8
3.5 1.5 -_ _ _ _ _ _ _ _ _ _
ExACflt SOL. ExACT 501.
1.0 a a 0 0 0 e a 1.0 -0-- -O--------E UPIND10 E FCT
0.5 0.5
0.0 -------a-e- 0.0 -0-- -0-0-*
0.0 0.1 02 .3 0.1 0. 6 01 0.T 0.0 1.0 0. 0.1 0.2 0.3 0.50 0.5 0.6 0.7 0.8 0.9 1.0
X
x x
1.5 3 E.
O
0 0 XENSOL c lyn
0.5 0.5
0.0 a0.0
. o -. a : : :: : -O.S-a a-
I I -0..5O. _JnL ,. J .,- J L0.5
-O.5. 0.1 0.2 0.3 0.13 0.S 0.6 0.7 0.8 0. 1.0 -0. 0.1 0.2 0.30. 06 0.8 .9 3.0
x x
(c) (
'Figur e 2.2 The comparis' on of litnl ica iil schcnieS ai, Cl(,'l.=0.2
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
AtF
Axt;p
wi+1/2>O0
(2.8)
AFi~
/\t i+1
where w is a variable which is used to determine the upwind direction. It can be defined as
wi+1/2--=(Fi+l1-Fi)(ui+ i-ui) (2.9)
2.2.2 The Second-Order Scheme
From the relation
!2F au Lgt Bt
a finite difference equation can be written as
un+ =upn-#.+4 .
(2.10)
F / =F"-I u(+/-- /) (2.11)
The one-sided second-order differences in flux conservative form is defined as
3 n+1/2 n+1/2
2 Fi -1F-1 wi+12i+/2= t (2.12)
3 n+1/2 n+1/2
i+1 i+2 wi+1/
2.2.3 The Combination-Hybrid Scheme
u7u'-u tbi+1/>+0i-1/2 (2.13)
Where, u? represents a first-order update. The second-order flux-correction is difined as
64+1/2= <+1l2- +1/2" (2.14)
Because of non-monotonicity of this flux-correction, a filter, i.e., a flux-limiter must be introduced (Boris and Book, 1973) .
64i+,12=sgn(64i+1/2)max{O, min(164i+1/21,sgn(6 i+12)(u+2-u +I),
sgn(6 i+/2)(u u i_)]} (2.15)
Then the hybrid scheme is completed as follows
u+1='u-6 46i+ .,+60i-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 centraldifference 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.
u7+ =- u (Fi+1/2-F -1/2- (2.17)
--ui- g 2- i-l/ D
where Fji=1/(Fj + F,), F,_q,2= (F _1+Fi). Fi+1/2-Fi-1/i 2 is the flux term. Different ways to construct the diffusion term D will result in different schemes.
If Di is expressed as
= ( ~ )2 A F (
(Dj= "t.-iA+j/'2(- Il-F!1)-Aj-j/.(F- -F7- )]
(.2.18)
A
Ai+1/2-- 1(u+ui+1)
21 +1
Ai-1/2- (ugi+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 (D 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--di+1/2-i-i /2 (2.19)
The diffusion function 0 can be written as
0i+1/2= -b(ai+1/2)[Aui+i/2 -Qi-1/2] (2.20)
where 0 is a correction term to prevent zero value of a.
jai ja>6
O(a)= 6=10-4 (2.21)
a'+62 |al<6
26
a- =aFu (2.22)
Bu
Q is the limiter function
Q+1/2=minmod(Aut/2, Auui/2,Au+3/2) (2.23)
Aui-1/2=ui-ui-1 (2.24)
The minmod function is defined as
minmod(x, y)=sgn(x)max{ 0, mnin[(jxj 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
ui =u7+Cu u-Cru;
(2.26)
where C=uAt is the local
d r e Ax
and r denote i- and i+1.
Courant number. The subscripts I
i-1 i i+1
Figure 2.3 QUICKEST method
u= (uC+u)- G(1-C)CURV- CGRAD
(2.27)
The curvature term, which must be upstream-weighted
for stability, is defined as
UR-2uC+uL
CURV={
ur,>O
(2.28)
URR-2UR+UC ur<0
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+1-ui (2.29)
ul 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.
.5
i.0 0.5
0.0
" 0. 0 0.1 0.2 0.3 0.11 0. 0.6 0.7 0.8 0.9 1.0
x
(a) Linear equation
0
EO
E- EXACT SOL.
0 )0 OUICKESI
0 -- ~~0 eooo
-0. 0o 0.1 0.2 0.3 0.4 0. 0.6 U.1 0 .a 0. 9 1.0
(b) Notn-I ii( ," eq'( tlli. nill
Figure 2.4 Ilesiults of QUICI(KEST method
0.5
S EXACT SOL.
o aucs
0 QUlIKESI
E)Eut 9--
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:
au wo
. i=0(3.1 .1)
Conservation of momentum:
au a(uu) + 8(uw)
3t Ox Oz
-lap 0 + a (A, u) (3.1.2)
Vox 57X a-" t)J Z OXO
Ovt+_Ouv~vw~ a (AA '-)+ a (A. v) (3.1.3)
a-t-x 9z Oz-Oz O-hxa
(3.1.4)
Conservation of salt:
s 8(us)+ O(ws)
-t 8x Oz
a 8 Os O, Os) = :: (a h X) + ( )7
Equation of state
P=po+ s
(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,K 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
Op7z 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)
8Kz
(2) At the ocean bottom there is bottom friction, and no flow coming through the bottom
(3. ..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)
A,(9 '- )= O
az az
az a v-A(u, v
K,-s= 0 L9z
A Ox )=o
(4) At the open ocean end, the ambient condition are specified as
A,(u, v, Ow_ (3.1.10)
(K a Ox
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
- Wb 8(uh) (3. 2 .1)
Ox
Below the free surface, the equation (3.1.1) becomes
S= w a(uh) (3.2.2)
Ox
Figure 3.1 Computation domain and grid system
where
u = Vertically averaged value of velocity in each layer
h = Layer thickness
Az+ 17 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)
Ot +x 0+az
(" x a-45-
where
- -uuh F = uvh
(3.2.4)
(3.2.5)
Jvh-u,w,+ubwb P (A uh+r -r
-fuh-vtw,+uw+ (Ah)+r,-r,
v = Vertically averaged value of velocity
in each layer
- Iuh S= vh
From Eq. (3.1.5)
8(sh) O(u2h)
t+ O+SeW-Sbwb Ot dx
= (K h + (K ,-(K )b
where
s = Vertical averaged value of salinity in each layer P = Integrated pressure over the layer thickness re,ry= Shear stresses
r ,, ,r= (Av'-u, ),
(A z A,-),' otherwise
T ,b 'yb I
(-Au, -Av)b
estuary bottom
3.3 Finite-Difference Equation and Procedure
3.3.1 Find 77+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)
Wn( 3.W3,.(u) U Wi/ ,j '-'Ax i,j- l-l,j)
(3.2.6)
(3.3.1)
At the bottom. w'=0. i.e. w'1==0
(3.3.2)
w = _(u "2 Ui-1.2)
At the surface, w,= Thus,
t+ += +Atwu
lh =Nt
(3.3.3)
3.3.2 Solve u" and v+1
The finite-difference equation for equation (3.2.3) is
?=+%At(FLUX ( ?) + ( ) (3.3.4)
where FLUX term can be separated into two terms
FLUX(l ) = FLUXC(W)+FLUXT(W)
(3.3.5)
where FLUXC(W ) are the net fluxes solved by central difference schemes and FLUXT( ) are the TVD type diffusion fluxes.
) F= F ,7,2-F /2j FLUXC(W )=--
(3.3.6)
where F i+1/2 2
- Fi,+Fi
F i-1/2 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/24 i+1/*>-Ti-1/29 -1/., (3.3.7)
where T is the right-eigenvector matrix of the Jacobian.
From equation (3.2.4)
S uu= uW, (3.3.8)
F uvh = [uW14j,
For a the "frozen" constant coefficient system, the Jacobian matrix can be written as
A-u 0 (3.3.9)
Obviously, the two eigenvalues of A are
At=u Az=u
and the corresponding right-eigenvector matrix and its
invers are
S O T-=[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/2(AW i+1/2-Q0 i+1/2)
AW ',+1/2,=T 1/2(W i+I-Wi)=W i+I-W i
(3.3.13)
Ii+/.2l
(3.3.14)
i+1/2"=
i+ 1/2 >10-4 P'i+1/ 21 < 10-4
where
-X 2 +62
a +1/2+
26
FLUXT lig)= 4 2-ofi-1/2)
6 is a small positive value.
i i+1/2 4in2 c w e+1/2
The l imiter function Q can be written as follows
Qi+1/2,=minmod(AW i_,-1/2, Ai+1/2' AW 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
1O8P O u
fyh-u+w,+uw PoOx +(A x )+rct-fre
-fuh-vw, +ubwb+ ) (Ah Vh)+TY-ryb
O N 1
where
1
U=(Ui,j"+ui,j+1)
vt=!(vij+ ij+ )
t (w + +i,) 14t=-" 2 i,jm i+l,j)
=1 (ui -+uij) Vb -i,j-1+ ig)
wb= (wij-1 +wi+.j-t)
constants, the stress terms can be
(3.2.5)
(3.3.16)
Assuming Ah and A, are
39
written as
9V = "\"Vi.j+I -Vij
= Aq-- = z
(3.3.17)
reb = Aulb = Avu zi-l
Ujz Az
Tyb bAv( = Av Az' j1
j ; b \z
Hence
a hub1 hi+t/2j-2ui hi-/2.j +ui-,jhi-3/2,j
--(As 'h)=AhAx
dx Ax2
The pressure grident -1P is calculated from the PoOx
surface to the bottom. At surface layer, from hydrostatic equation (3.1.4)
(3.3.19)
p(z)= pg(q-z)
p= pg(q-z) ax 8r9
-- 0dz= pg(q-z)dz
Ox j f8x x
(3.3.20)
(3.3.21)
By Liebnitz's rule,
(3.3.18)
api d 1 2P;
=pg0,+Az) +.-g(Az +) Ox
(3.3.22)
)IN+1 i-I + I ,P lj--pi-IJ = pg(7,i+Az) +-) g(rqi+ z ")2pi+ '_P
=p(?~+z)2AX 4 'Ax
At the subsurfaces,
aP P SP BP Az
Si,j --- x ij+1 gAz = i,+t g-4--(Pi+,ji-Pi-1j)
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+l=h h FLUXC(sh)+FLUXT(sh)+RHS(sh)
s,=j s+t h 3.3.24)
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)- swt-swb)
1 ((sh)i+l.+(sh)ij (sh)i.j+(sh)i-j uii) = --x 2 ui'i 2
- (st'W'-sbwb) (3.3.25)
where
st (sij+Ij+sij)
wt =wi, j
s=.(sij +sj-1) "b i.,j-I
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)= (4i+1/2 -i-1/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 K;1.
Qi+1/,=minmod( (sh)i-/2, '(sh)j+1/2, A(sh)i+a/2 )
A(sh)i+1/2. ((sh)i+(sh)i+j)
RHS(sh) =Kh +K V
0x- 8z
.2O 1 (sh),,j-2(sh)i,j+(sh)i-1 ,j)
( ~= vK, ((sh)i+-2(sh)ij+(sh)j-1
[ (9=z,, 2 Az 2~ ~ hij .-
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.
u uAx
(AhKIh)>2
(Ah Ax2
u2At
(A,Kh) 2 (As,K~a)> u'.
(3.4.1)
(3.4.2)
(3.4.3)
Inequality (3.4.1) is restrictive to the selection of Ah and K. 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)
43
If this number is greater than 2. equation (3.4.2) will not satisfy CF'L condition. IIn 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) wh ic h is due to the explicit treatment of the advection term. If \t 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 < Axs. (3.4.4)
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 w i th 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
number of grid points in the region ,;here 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 (Axi 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 P+1/2 I
xi- Ax xi Ax+1 xi+l
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 r. the density on either side o- x, i.e. _x. an i
-x._. ,is changed. the mass center should boe changed. If P,_ on .x. is larger than P./. the mass center should be moved to the left to a ne p-osition :,),w (Figure 4.2).
Ni-
new mass center (x.) e,
old mass center (xi),od "i-1/2 w+1/2
I Wi41 :
X....
Figure 4.2 Moment Balance Principle
As shown in Figure 4.2, wi_/,(wi-1/2=Pi-t/2xi) wi-1/ (wi+1/2=--Pj+12) are weight functions on Axi and Axi+,.
and
Let us define
10=11+1 ,= (Ax +Axyt)
(4.1.1)
Acoording to moment-balance principle,
w i-1/211 = "i+ /2 12 2
(4.1.2)
Wi+l 1iL~ /2
-Wi+_- .- (4. 1 3)
Thus.
1 Wi+1/a(AXi +Axi,1)
I =- I ( + ( ) (4.1.4)
Therefore. the mass center will be shifted from (xi),Id to
(x) ,new.
(xi) new=xi-I+1 + I
= (x),ta dxgIAx,+ L1*
=(x)od- !Ax+ I 1 (4.1.5)
The adaption will be advanced from i=1 to i=Imax. The adjustment of each point is based on the previous point positions. Then repeating the same procedure to find 1 and
(x,)n..e until Ei(Axi)new-(AXi)o20dIse (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 developes. 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
discontinuity is chosen as the sensor S(z). 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 Axi.
W i-1/2.= P,_1/2..Axi (4.2.1)
Where. the density function P is some kind of combination
of the first derivative of tihe sensor S(z), Sz, and the second derivative of the sensor S(x). S,, 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=Y( + e l (4.2.2)
p=1 +aIS.I+3!SI (4.2.3)
l+ISxrI+IS:I
The arbitrary constants a and 3 are the adaptation parameters. When a increases, the sensitivity of P to Sz,
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 ISz-I and ISZj 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 f.
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. Unfourtunately, since the adaption is performed based on knowleadge 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 Dockelie. 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 1-2.
wi=!(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 maching 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 ,,rid .ize. and 3 can also affect the
smoothness of the grid cells.
4.5 Verification of The Adapti,.e 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.
The inta o(4.5 1)
The initial conditions are
Case (1)
1 x-U linear O
-1 x=1
Case (2) u = cos[(2k+I),rx] O
Case (3) 1 x=O
(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
- 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 0.00
-0.50
-1.00
-1.50
-2.00
-1.00
-1.50
-2.00
IIME- 0.1193
--- --
--_-LfLLLLLLUAALLLLL1Li
X
I It-L ). 1 JIUt)
LLLJ .iWLIIIII LL~L L.L -
"igur ia e -1.3 (ais (1 ) (In aLdapt.,ive g' il I i t. 20
- I hE- 0.0000
I I 11111I1I1 II 11
2.00 1.50 1.00 0.50
- 0.00
-0.50
T IHE- 0.25011 I I I I I I I lIIII I I I I I
I IHE- 0. 1916
1.001 --
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
-2.00
.LLOIJ_LL I L LL..LJLI
X
l*'igns-c' -1.4I u c (2) on ada.pliv: g rid lain.. t 20
2.00 1.50 .00 0.50 S0.00
-0.50
-I.00
-1.50
-2. 00
T I hE 0. 11167
I I I I I ll ia I i i 1 1 I i
TIME- 0. 7631
I 1 I lLL lJ_ .. .. L
2.00 1.50 1.00 0.50
- 0.00
-0.50
- I Ao
-1.50
-z.oo L llIJIlII I I
2.00,
1.50
1.00 0.50
0.00
-0.50
-1.00
-1.50
-2. 00
STIHE- 0.3724
.-.LLLALL ILLI
,. L LI 1 1 1I 1 I lilll I I I
Iigire '1. -5 Ca.LSe .(3) oil adap-Live' grid point.H 20()
- T IME- 0.0000
- ~ fJ~fLJ LJ Lfi
2.00 1.50 1.00
0.50 S0.00
-0.50
- I.00
-1.50
-2.00
2.00 1.50 1.00 0.50 0.00
-0.50
-1 .00
-1.50
-2.00
TIMlE- 0. 1101
F-
I I I I L
TI HE- 0.2216
Il li ll i IIll Il l I I I I
,
------
I I I I
j)IP II' I~ 1.t ) ti i0t I 0z:)It( I to!I.I Ion j! p* j, n.ggI. g u- pun o, pt.-W inan.y~ itut (i w ) OR 1 W) O, Iu' OE It.. uslof.ut uo (1) a- o (W)
I i1 T ff IT- T1TIT fTT T FII I I-TTT1
\/
\-
--FI T TT-T-T-1---1 I I T T-
I~, ~-x1J
00."
OS" 1O0"I(S 0
00*I O0 "
000
00*1 OS1I
an -z
o0 "e00,1OS'I
00 0 C: Os-I
00"1 0o'? 00 "0
TIIII lT IITIlIIITIIIITTI I IIII T n IIII II II llii l 11111111lii111111ll
"lg-XN
OS"too" l00* I 00 1OsVO
0 0 c:
00 1 OS I Do*?
0S" 100 *00* I
OS't
000 C 000
0011
Os., 00*?
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 discribed earlier in this thesis. Three of the problems contain simplifying assumptions wh i le 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 begining 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
BuD+OD =0 ax Ox
auDuuD-fD= IgOD
Ot x ON
+vD8 uvD+fuD=0
Bt x
(5.1.1) (5.1.2) (5.1.3)
where
62
xO. t=O
If Coriolis effect equations are reduced to
is neglected.
the governing
D.=OuD=0 8t x
8uD+ OuuD 12 Ot Ox 2O
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=O.025. The time steps are limited by the following condition:
(5.1.7)
At < CFL( Ax )
nmazjlu+cl
where c=.g 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 (..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 discontinuit" 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 coresponding 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 origional one. i.e.
W"= "+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
10 '-vI L-- ..2 I' -' -- -,-1-r--I -1--i- -l -1'-0 -I-s- ( -, -: :v I I' I o I" I s- I"
9. 9 -......
10 U---T
o-u
7 76
6 ..-.J LLL J-- _L J ..aL L --.-t ----..i -.-_I --- -- -- ... al_-.-J a |3 --- 3
2 2-0 --0 --
-12 -10 -8 -6 -i -2 0 2 i 6 6 10 12 -12 -tO 8 -6 -' -2 0 2 4
X x
( ) (b,)
Figiire 5.2 "[rcakiing dam" problemiii
(a) Aialyt.ic solutitiOiI at i it i ntii.urval: (1 i. i ile:
(b) Niimn erical solkiutiolio i cori p-ipd-li ig t (1.)
' I "l I I I "
I
. .. . i . . .
i-. I 8- 1 - 1-6 t0 10 12
0" OI
I I S 5 V111T1fFr 1rg1yrN
-I I 0 -.. .
-12 -10 -8 -6 -4 -2
~1
0
2 1 6 8 10.
x
(G)
22
20 18
!0 -------------.
8 --. .--
I I-4 - -10 -------~-
a-J
03
2 -2 I-L- I-- I .. I.. I.... I._ I. I 1 1
2 -12 -ID -8 -6 -4 -2 0 2 I 6 0 t10 12
x
'igtre 5.3 'T'he 11111nerical solutions For 1.1he laytI (Ilet.bI *ithl
ritionl e cffeoct;s ii (a) 1.z 1 it (I (I) t.=2
Figure 5.4
Nonratating "breaking dam" problem at t=1 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.
p
X
Figure 5.6 Density Distribution
As shown in Figure 5.2-1, the density distribution follows p=(1--ex), where a=!- is a measure of the Pax
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 p=1.058x10-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 (Jga)2 According to Simpson and Linden's(1988) theoretical computation, the initial angle of inclination is given by tanO=,r. 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-5cm-, p=I .058xO-3g/cM3 and a=15x10-cm-1, P=I .0
-I
a-.
Iigiurc 5.7 Nimine ical t inc-se(jetnce showing t.1 isopyciial-s
with f =1.058x 10()-3g/cm a=1I5xlO-'(:cmat. t;=8.0s, 12.3s, 16.7s, 20.9s and 2-1.9s.
----------
~~~4i ~ ~ ~ ~ ~ .... ...~rJwrrrwr rrY r ~ ~ r.
Figure 5.8
Corresponding to Figure 5.7 a t.iume-scqucm,,t:t showing the motion of the isoy:u Is in laboratory experiments from Simnpson uud Lin(Ien(1988). Only the right-an1ld I i I. 'f tOlF tank is visible and the vertical l ines arc aLt. 10cm intervals
Theory
a=6.5x10-5crm-1
----- -a=15x10-6cm-'
o o I.ab 1 XX I Lab 2
- -~-0
o00
- --% -- 0
X X
xx
-0.60 -0.20 0.20 0.60()
log(tL" )
Figure 5.9
A log-log plot of the angle 0 oF incl itnti.iont (,o the isopyc nals to Ih e vertical against. Ith. climensionless tI'=(ya)/2 t
2.0
1.6
1.2
0.8
0.4
0.0
1.00
1.00
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 P-hysical 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 experimen t, the gravity front is formed by releasing fluid of a different density from behind lock-gates. Since the volum e 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 pg=1 g/cm3
H= h
Figure 5.11 Laboratory Experiment Case 2
Where xo=50 cm L=15O cm
ho=7 cm H=7 cm
Po=1.048 g/cm Pt=l g/cm3
For realistic large-scale river and ocean flows, there
Case 1:
Case 2:
L
Ixo
are nor only bottom plumes, but also surf-ace bouyant plumes. Therefore. the next t-wo cases describe zhe large scale bottom and sur-ace gravity currents.
Case 3:.
8 L
Figure 5.12 Large Scale Bottom Plume
Where .xo=3 km L=40 km
h0=7 m H=10 m
po=1.000546 g/cm
p,=1 g/cm3
Case 4:
L
Pl H
.Figure 5.13 Large Scale Surface Plume
Figure 5.13 Large Scale Surface Plume
where xo=3 kin L=40 kn
h0=7 n H=10 n
Po=1 g/cm
pl .000546 g/cn3
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.16 and .5.19 also can give the same results about the front
formation.
The solution shown in Figures 5.1S and 5.19 demostrate 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
-f (5.3.1)
Theoretically, the bottom front has larger velocity than surface front due to the reduced acceleration of gravity gl=Pg. 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 \'elocity
Wang(1984) investegated 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
I
Case Theory Lab Numerical
1 11.79 12.24 12.5
2 8.16 5.S 6.9
Case Theory I Lab Numerical
3 8.709 9.68 10.68
4 8.707 9.67 8.47
Table 1.
Bottom
(IsI
Ii
Vgi fe 5. 14 Si(tiWgraphs of a vo I n of sal t. wate:r wri l
lio=7cm, II=49cm, xo=50cu, g) =17cus co I laps i n to
fresh water at 2.2s, 3.8s and 6.Os afte relvase 'IThe dLted vertical line inldiclate.s the position
of tf 1e I ock gate, time cimldwal I is at. fr left., iil
the thiin v rIical ines are t. IO n ii e a l
(Illtiian and Simpson(1983))
Figre 5.15 Shadowgraphs of a vole ait t i I
. 1
ho/11=1 1n-7cm, x0=50cm, g' -17cnSs col ig into fresh water at 1.7s, 7.7s and 10.7s a.tcr release (Ilottman and Simnpson(1983))
I I I I I. I- I I I I -I I - I I I I I
. .... . . . .. . .. .. .
", .: 1 ." .. . .2 ... .. . . . .
... . . "a ] ~ . .. . ..-... .. . ... ..... . ..=== = == = == = .. . . .. . .
. .. . . . . . . . . .
. .I . . . .. .
... .. ........,. ..,. ., .| a S | . .. ,, ,. :... .. ..... . . .. . . . . ."
.'" . -. . .. : :- . . ..... ... .q . . ...... . .: : ::. : : : . ... . . .
.. .. . ,. .. . . . . . . . . . . ..,,,, J . . . '.' ". . . . . . .
: : ,,: ." ." . ., ....... . . .: ... ... ... :...... . .
I I I I
I I I I
I I I I I IWWJ
LLLJ LL LL~
.. . . .. .. . . ... . ... .. .. . .. . ....:. .. . .. ..
. ... .. .. .. .. . . ,. . .. + , .. .. . . .. . .. . .. . .
...: ... .. .. ................... . .! ., . . . . ........ ...,.- ........... .. -. ..... ... ... . . . . . . . .
: "= t ,:' : ....... ...-. ........ : :.'::' :: ..... ..:: :: ......." :: : :' : . .
Ntime r i c al so I ii t i o1S co 1 II cslPO)I(I g L)
39.00
32.500 26.0
13.00
6.500
ii r'e 5. 16
- . .. ".. . . .. i 1 a I I I I
Filir,,'. 'r. l'l
. . . . . . . . . . . . . . . . . . . . . . . .
. -g : -: '! i : : : : : : . .: : : : .: : : : : : .: : : .. . . . ..w
. .. . . . .-.. . . . . .. .... . . . : . . . . . .
:: :.: k :'F :::> k :::::::: :::: :: :: :. : .: : : :. : ........... .. .. ... ..... .
- -- -"4 : 'I' |" L :!i;.ii.; ;:..!i i................ . ... .- .. ... .:: .... .. I : : : : i
"M HMI I
F....................
.:.....................
':1.. .. . . .. .
. . . . . .
L ... ........ .
. .__ _ _ ._ _ ... ...
. ~ ~ ~ ~ ~ ~ ~ ........
50.83
32.50
233
14. 17
Figure 5.17 Numerical1 solutions corresponding to Figure 3.15
0 40km
I I I
_ _.. __.......
am ...........O.5540 0.5080 0. 20 0. 4160
Figure 5.1S Numerical solutions of Case 3
at 20hrs. 40hrs, 60hrs and 80hrs
I I t ~ I I I I t i I t L I I
I-
~I.
-4
]
I I r I L H
E 1. 1 i I. i t 11 1 E I L i 1 1 i I. 1 L I, f
I. ~ U i~ !1
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . I . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . .. . . . .. .. . . I . . . . . . . .
. . . . . . . . . I . . . . . . . . I . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. ... ......... ...................
. . . . . . . . . . . . .
...... ................. ..
............................
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . I . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . I . . . . . . . I .
. I I .
. : .. .
. I I I .
. .
. .
. .
. ..... ... .
. .
. .......................... .
. .
. .
. .
. . . . .
. ....... .. .. .
. ................... ............
.. ........
.........
.............. ...
. . . . . .......... ...........
0 40 km
I E i. I L I- L 1. -1 1: .4:.. 1- 1: t I 1.: .1:.: V.
. . . . . . . . .
. .
. .
. ..
. .
. .
. .
. .
. .
. .
. .
. I I .
. I *
. .
. .
. .
. M .
. .
. .
* .
. I .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. .... .. .
. ......... .
. .
. .
. . . . . . . . . . . . . . . . . . .
. ........
. ............. ..........
. ...........................
........................... .
. . . . . .
.... ...... .. ........................................
. ........................ ...............
..................
............ 0 .......
-10M - - - -- --- -- ---- ----
0. 1850
0. 1520
0.1180
0. 840E-0 I
Figure 5.19 Numerical solutions of Case 4 at 2Ohrs,40hivs, 6Olirs
and 80hrs
. ..... .........
............. .. ... .
... ........... ... ...... .
. . . . . . . . . . .
8*'
II
0 2030
0 o 20o
aI
Ia 20 10
Figure 5.20 Numerical Solution for Surface Plume
at 20hrs, 40hrs, 60hrs and 80hOrs by Wang(1984)
0 to 20 30
.4
~1
.~ ..0.
0
2
0.
E
.~ ~'.
0.
0
2
E.
SI
0
eL
I.
20 30
S10 20 O30
Figure 5.21 Numerical Solution for Bottom Plume
at 20hrs. 4Ohrs, 60hrs and SOOhrs by Wang(1984)
i
Figure 5.22 Numerical .Solution -for Surface Plume at 20hrs,
40hrs, 60hrs and 80hrs Coresponding to 5.3-11
|
