UFL/COEL2000/014
A THREEDIMENSIONAL CONSERVATIVE EULERIANLAGRANGIAN MODEL FOR COASTAL AND ESTUARINE CIRCULATION
by
Jun Lee
Thesis
2000
Copyright 2000
by
Jun Lee
ACKNOWLEDGMENTS
First of all, I would like to thank my advisor, Dr. Sheng, for his guidance and financial support throughout my graduate research. He encouraged me when I was frustrated doing my research. In addition, I would like to thank my committee members, Dr. Dean and Dr. Thieke, for their careful review of this thesis.
I would like to thank all of colleagues in my research group and in particular :Justin and Sun, for their help with
discussing, developing, debugging and verifying the model; Dave and Du, for doing field work together; Chenxia, Vadim and Jack, for studying and working together; and my Korean friends Kij in and Yeonshik, for studying and relaxing together.
I would like to thank everyone in the administrative offices, Becky, Lucy, and Sandra. Special thanks are sent to
Helen who helped me to f ind many references. Grateful thanks go to everyone in the Coastal Engineering Laboratory and particularly Sidney and Vic who are my diving buddies.
Finally, I would like to thank my parents, my wife and my son, and God, for without their love and endless support, I would never be where I am.
TABLE OF CONTENTS
. . . . . . . ii
ACKNOWLEDGMENTS LIST OF TABLES LIST OF FIGURES ABSTRACT ....
1 INTRODUCTION
1.1 Objectives ........ ..................
1.2 Organization of This Study .... ..........
2 GOVERNING EQUATIONS AND BOUNDARY CONDITIONS IN THREEDIMENSIONAL CARTESIAN COORDINATES .... ..........
2.1 Introduction ....... .................
2.2 Governing Equations in Cartesian Coordinates 2.3 Boundary Conditions in Cartesian Coordinates
3 A THREEDIMENSIONAL NUMERICAL MODEL ... ..........
3.1 Introduction ....... .................
3.2 Finite Difference Form of Momentum Equations
3.3 Solution Algorithm ..... ..............
3.4 Conjugate Gradient Method ..........
3.5 Salinity Scheme . . . . . . .
* V
* vi
* xvi
3
5
4 EULERIANLAGRANGIAN (SEMILAGRANGIAN) METHOD ...... 34
4.1 Introduction ...... ................. ..34
4.2 Review of EulerianLagrangian Method (ELM) .... 34
4.3 Derivation of EulerianLagrangian Method (ELM) 37
4.4 Conservative EulerianLagrangian Method (ELM) 43
4.5 Numerical Experiments ..... .............. 51
4.5.1 Linear Advection Case ... ........... 51
4.5.2 Conservation Test .... ............. 59
4.6 Momentum Advection Equation ... ........... 61
4.6.1 NonLinear Advection Case .......... 61
4.6.2 Momentum Conservation Check ......... 70
4.7 TwoDimensional Conservation Test..........71
4.8 ThreeDimensional Conservation Test .........76
4.9 Summary.......................81
5 MODEL VERIFICATION WITH ANALYTICAL SOLUTION..........82
5.1 Introduction......................82
5.2 Wind Setup......................83
5.3 Tidal Propagation..................86
5.4 Tidal Propagation with NonLinear Advection ... 98 5.5 Tidal propagation with Coriolis Effect ........104
5.5.1 Kelvin Wave Propagation............104
5.5.2 The Governing Equation for Kelvin Waves ..105
5.5.3 Comparison between Numerical and Theoretical
Solutions...................109
5.6 Tidal Propagation with Bottom Friction Effect 117
5.7 DensityDriven Circulation.............125
5.7.1 Horizontal Diffusion Test............125
5.7.2 Density and Tidally Driven Circulation 129
5.8 Wetting and Drying Test over Tidal Flats .......132
5.8.1 Determining Flooding and Drying Cells . 133 5.8.2 Comparison with Analytical Solution . . 135
6 COMPARISON WITH CH3D....................139
6.1 Introduction...................................139
6.2 Comparison with CH3D in a Rectangular Basin . 140
7 APPLICATION OF MODEL TO LAKE OKEECHOBEE...........150
8 CONCLUSIONS AND FUTURE STUDY................161
8.1 Conclusions....................161
8.2 Future Study.....................162
APPENDIX A NUMERICAL SCHEMES FOR LINEAR ADVECTION EQUATION 165
APPENDIX B DESCRIPTION OF PROCEDURES OF THE PRESENT MODEL 171
LIST OF REFERENCES.....................175
BIOGRAPHICAL SKETCH.....................182
LIST OF TABLES
Table oacre
Table 5.1 Summary of test for analytical solutions and model simulations........................83
Table 5.2 Comparison between analytical and simulated wind setup...........................84
LIST OF FIGURES
Figure page
Figure 3.1 Schematic diagram of stencil .. ......... 15
Figure 4.1 Characteristic lines and their feet ..... 35
Figure 4.2 An Eulerian Lagrangian Mesh .. ......... 38
Figure 4.3 Analytical solution of linear advection equation
at T = 50 Sec ........ .................... 52
Figure 4.4 Analytical solution of linear advection equation
at T = 250 Sec ....... ................... 52
Figure 4.5 Analytical solution of linear advection equation
at T = 500 Sec ....... ................... 52
Figure 4.6 Results of the 1D linear advection equation
obtained with upwind method at T = 50 Sec ...... 53
Figure 4.7 Results of the 1D linear advection equation
obtained with upwind method at T = 250 Sec ..... 53
Figure 4.8 Results of the ID linear advection equation
obtained with upwind method at T = 500 Sec ..... 53
Figure 4.9 Results of the 1D linear advection equation
obtained with the Leith's method at T = 50 Sec . 54
Figure 4.10 Results of the 1D linear advection equation
obtained with the Leith's method at T = 250 Sec . 54
Figure 4.11 Results of the 1D linear advection equation
obtained with the Leith's method at T = 500 Sec . 54
Figure 4.12 Results of the 1D linear advection equation
obtained with the TVD method at T = 50 Sec ...... 55
Figure 4.13 Results of the 1D linear advection equation
obtained with the TVD method at T = 250 Sec. . . 55
Figure 4.14 Results of the 1D linear advection equation obtained with the TVD method at T = 500 Sec. . . 55
Figure 4.15 Results of the 1D linear advection equation obtained with the ELM with linear interpolation at
T = 55 sec ........ ..................... 56
Figure 4.16 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at
T = 253 sec ......... ..................... 56
Figure 4.17 Results of the ID linear advection equation
obtained with the ELM with linear interpolation at
T = 506 Sec ......... ..................... 56
Figure 4.18 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at
T = 55 Sec ......... ..................... 57
Figure 4.19 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at
T = 253 Sec ......... ..................... 57
Figure 4.20 Results of the ID linear advection equation
obtained with the ELM with linear interpolation at
T = 506 Sec ......... ..................... 57
Figure 4.21 Results of the 1D linear advection equation
obtained with the Conservative QMSL at T = 55 Sec . 58
Figure 4.22 Results of the 1D linear advection equation
obtained with the Conservative QMSL at T = 253 Sec 58
Figure 4.23 Results of the ID linear advection equation
obtained with the Conservative QMSL at T = 506 Sec 58
Figure 4.24 Computational cell used in conservation test 59
Figure 4.25 Conservation error of the conventional ELM 61
Figure 4.26 Conservation error of the QMSL ....... 61
Figure 4.27 Results of the nonconservative upwind method
at T = 0.5 sec ....... ................... 63
Figure 4.28 Results of the nonconservative upwind method
at T = 1.0 sec ....... ................... 63
vii
Figure 4.29 Results of the nonconservative upwind method
at T = 1.5 sec ....... ................... 63
Figure 4.30 Results of the conservative upwind method
at T = 0.5 sec ....... ................... 64
Figure 4.31 Results of the conservative upwind method
at T = 1.0 sec ....... ................... 64
Figure 4.32 Results of the conservative upwind method
at T = 1.5 sec ....... ................... 64
Figure 4.33 Results of the Leith's method for nonlinear
advection equation at T = 0.5 sec .. .......... 65
Figure 4.34 Results of the Leith's method for nonlinear
advection equation at T = 1.0 sec .. .......... 65
Figure 4.35 Results of the Leith's method for nonlinear
advection equation at T = 1.5 sec .. .......... 65
Figure 4.36 Results of the Godunov method for nonlinear
advection equation at T = 0.5 sec .. .......... 66
Figure 4.37 Results of the Godunov method for nonlinear
advection equation at T = 1.0 sec .. .......... 66
Figure 4.38 Results of the Godunov method for nonlinear
advection equation at T = 1.5 sec .. .......... 66
Figure 4.39 Results of the nonlinear advection equation of
ELM with linear interpolation at T = 0.5 sec . . 67
Figure 4.40 Results of the nonlinear advection equation of
ELM with linear interpolation at T = 1.0 sec . . 67
Figure 4.41 Results of the nonlinear advection equation of
ELM with linear interpolation at T = 3.25 sec . 67
Figure 4.42 Results of the nonlinear advection equation of
ELM with quadratic interpolation at T = 0.5 sec . 68
Figure 4.43 Results of the nonlinear advection equation of
ELM with quadratic interpolation at T = 1.0 sec . 68
Figure 4.44 Results of the nonlinear advection equation of
ELM with quadratic interpolation at T = 3.25 sec . 68
viii
Figure 4.45 Results of the nonlinear advection equation of conservative QMSL method at T = 0.5 sec.........69 Figure 4.46 Results of the nonlinear advection equation of conservative QMSL method at T = 1.0 sec..........69
Figure 4.47 Results of the nonlinear advection equation of conservative QMSL method at T = 3.25 sec.........69 Figure 4.48 Conservation error of the conservative QMSL 70
Figure 4.49 Numerical solution of 2D nonconservative ELM 72 Figure 4.50 Conservation error of 2D nonconservative ELM 72 Figure 4.51 Contour map of numerical solution of nonconservative ELM......................73
Figure 4.52 Numerical solution of 2D conservative ELM 74 Figure 4.53 Conservation error of 2D conservative ELM 74 Figure 4.54 Contour map of numerical solution of nonconservative ELM.....................75
Figure 4.55 Numerical solution of 3D nonconservative ELM 77 Figure 4.56 Conservation error of 3D nonconservative ELM 77 Figure 4.57 Numerical solution of 3D nonconservative ELM in
XZ plane.........................78
Figure 4.58 Numerical solution of 3D nonconservative ELM in
YZ plane..........................78
Figure 4.59 Numerical solution of 3D conservative ELM .. 79 Figure 4.60 Conservation error of 3D conservative ELM .. 79 Figure 4.61 Numerical solution of 3D conservative ELM in XZ
plane..........................80
Figure 4.62 Numerical solution of 3D conservative ELM in YZ
plane..........................80
Figure 5.1 Comparison between analytical and numerical
solutions of water level for wind setup at three locations
within the domain(2D case)................85
Figure 5.2 Comparison between analytical and numerical solutions of water level for wind setup at three locations
within the domain(3D case) 85
Figure 5.3 Schematic diagram of a rectangular basin 87
Figure 5.4 Computational grid system for tidal propagation
test 88
Figure 5.5 Twodimensional simulation of water surface
elevation of tidal propagation at three locations with
O= 1.
Line : Analytical Solution
Circle : Numerical Solution 90
Figure 5.6 Twodimensional simulation of water velocity
of tidal propagation at three locations with O= 1.
Line : Analytical Solution
Circle : Numerical Solution 91
Figure 5.7 Threedimensional simulation of water surface
elevation of tidal propagation at three locations with
O= 1.
Line : Analytical Solution
Circle : Numerical Solution 92
Figure 5.8 Threedimensional simulation of water velocity
of tidal propagation at three locations with O= 1.
Line : Analytical Solution
Circle : Numerical Solution 93
Figure 5.9 Twodimensional simulation of water surface
elevation of tidal propagation at three locations with
O= 0.501.
Line : Analytical Solution
Circle : Numerical Solution 94
Figure 5.10 Twodimensional simulation of water velocity
of tidal propagation at three locations with 0 = 0.501.
Line : Analytical Solution
Circle : Numerical Solution 95
Figure 5.11 Threedimensional simulation of water surface
elevation of tidal propagation three locations with
O= 0.501.
Line : Analytical Solution
Circle :Numerical Solution
Figure 5.12 Threedimensional simulation of water velocity
of tidal propagation at three locations with 0 =0.501.
Line Analytical Solution
Circle Numerical Solution...................97
Figure 5.13 Twodimensional simulation of water surface
elevation of nonlinear advection at three locations with
0 =1 .
Line Analytical Solution
Circle Numerical Solution...................100
Figure 5.14 Twodimensional simulation of water velocity
of nonlinear advection at three locations with 0 = 1
Line :Analytical Solution
Circle :Numerical Solution...................101
Figure 5.15 Threedimensional simulation of water surface
elevation of nonlinear advection at three locations with
0 =1 .
Line Analytical Solution
Circle Numerical Solution...................102
Figure 5.16 Twodimensional simulation of water velocity
of nonlinear advection at three locations with 0 = 1
Line :Analytical Solution
Circle :Numerical Solution...................103
Figure 5.17 Twodimensional model results of water surface
elevation at three locations with Coriolis effect when
0 =1 .
Line Analytical Solution
Circle Numerical Solution...................112
Figure 5.18 Twodimensional model results of water velocity in
X direction at three locations with Coriolis effect when
0 =1 .
Line :Analytical Solution
Circle :Numerical Solution...................113
Figure 5.19 Threedimensional model results of water surface
elevation at three locations with Coriolis effect when
0 =1I.
Line Analytical Solution
Circle :Numerical Solution14
Figure 5.20 Threedimensional model results of water velocity
in X direction at three locations with Coriolis effect when
06=1 .
Line :Analytical Solution
Cicl Nue.a Solution .. .. 115
Figure 5.21(a) water level(contours) and velocity field vector
after 48 hours of simulation with Coriolis effect when
0 =1 .
Contour lines represent water surface elevation . 116
Figure 5.21(b) water level(contours) and velocity field vector
after 60 hours of simulation with Coriolis effect when
0 =1.Contour lines represent water surface elevation . 116 Figure 5.22 Twodimensional model results of water surface
elevation at three locations with bottom friction when
06=1I.
Line :Analytical Solution Circle :Numerical Solution12
Figure 5.23 Twodimensional model results of water velocity
in X direction at three locations with bottom friction when
0 =1I.
Line :Analytical Solution Circle :Numerical Solution
Figure 5.24 Threedimensional model results of water surface
elevation at three locations with bottom friction when
06=1I.
Line Analytical Solution
Circle Numerical Solution
. . 122
. . . 123
Figure 5.25 Threedimensional model results of water velocity
in X direction at three locations with bottom friction when
06=1I.
Line :Analytical Solution Circle :Numerical Solution
Figure 5.26 Comparison between analytical and numerical
solutions of diffusion equation in X direction.
Figure 5.27 Comparison between analytical and numerical
. . . 124
127
xii
. 114
. 121
solutions of diffusion equation in Y direction . 127
Figure 5.28 Numerical result of temperature distribution in
rectangular plate ....... .................. 128
Figure 5.29 Comparison between the numerical and analytical
solution of the density driven currents.
Line Analytical solution
Circle Numerical solution .... ............. 131
Figure 5.30 Definition of total computational depth . 134
Figure 5.31 Wave propagation on a linearly sloping beach 135
Figure 5.32 Nondimensional comparison between wave profiles
as predicted by theory and numerical model of Wetting and
drying(time=0 )T/2) ...... ................ 137
Figure 5.33 Nondimensional comparison between wave profiles
as predicted by theory and numerical model of Wetting and
drying(time = 2z/3 g ) ..... .............. ..138
Figure 6.1 Comparison to CH3D with wind forcing ..... 140
Figure 6.2 Water surface elevation of propagation test with
CH3D when 0 =1.0
Line CH3D
Circle Present model ...... ................ 141
Figure 6.3 Water surface elevation of propagation test with
CH3D when 0=0.501
Line CH3D
Circle Present model ...... ................ 142
Figure 6.4 Depthaveraged water velocity of propagation test
with CH3D when 0 =0.501
Line CH3D
Circle Present model ...... ................ 143
Figure 6.5 Water surface elevation of nonlinear test
with CH3D when 6 =1.0
Line : CH3D
Circle : Present model ...... ................ 144
Figure 6.6 Water surface elevation of Nonlinear test
with CH3D when 6 =0.501
Line : CH3D
xiii
Circle :Present model14
Figure 6.7 Depthaveraged water velocity of nonlinear test
with CH3D when 06=0. 501
Line :CH3D
Circle :Present model.................146
Figure 6.8 The test of Coriolis effect with CH3D when 06=1.0
Line :CH3D
Circle :Present model.................147
Figure 6.9 Comparison of water surface elevation with
Coriolis ef fect with CH3D when 6 =0. 501
Line CH3D
Circle Present model .............
Figure 6.10 comparison of depthaveraged velocity with
Coriolis ef fect with CH3D when 0 =0. 501
Line CH3D
Circle Present model.................149
Figure 7.1 The Bathymetry of Lake Okeechobee
(From the 1987 contour map of Lake Okeechobee) . . 151
.148
Figure 7.2 Computational grid
. . . 151
Figure 7.3 Numerical solution of total water depth as a
function of time at location A, B, C, D, E, and F in Lake
Okeechobee subject to an impulsive wind .........153
Figure 7.4 Numerical solution of water surface elevation
with time at location A, B, C, D, E, and F. .....154
Figure 7.5 Velocity fields and dry regions at time 1 and 3
Hours. Blank cells indicate the dry cells........155
Figure 7.6 Velocity fields and dry regions at time 6 and 9
Hours. Blank cells indicate the dry cells ........156
Figure 7.7 Velocity fields and dry regions at time 12 and 15
Hours. Blank cells indicate the dry cells ........157
Figure 7.8 Velocity fields and dry regions at time 18 and 24
Hours. Blank cells indicate the dry cells.........158
Figure 7.9 Velocity fields and dry regions at time 30 and 36
Hours. Blank cells indicate the dry cells........159
xiv
. 145
Figure 7.10 Mass conservation error . ..
. .160
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
A THREEDIMENSIONAL CONSERVATIVE EULERIANLAGRANGIAN MODEL
FOR COASTAL AND ESTUARINE CIRCULATION
By
Jun Lee
August 2000
Chairperson: Y. Peter Sheng
Major Department: Civil and Coastal Engineering
This study presents a threedimensional conservative Eulerian Lagrangian hydrodynamic model with wetting and drying which is significantly enhanced over the nonconservative EulerianLagrangian model developed by Casulli and Cheng. Semiimplicit finite difference method has been used for the numerical solution of threedimensional shallow water equations. The advection terms and horizontal diffusion terms in the momentum equations are discretized with the explicit conservative EulerianLagrangian Method (ELM). The barotropic pressure gradient terms in the momentum equations and the velocities in the vertically integrated continuity equation are discretized with the escheme, while the baroclinic
xvi
pressure gradient terms in the momentum equations are discretized with the explicit scheme. The vertical diffusion terms are discretized with the implicit numerical scheme.
Conservation is an important property when studying the transport of salt, heat, and other conservative species in coastal and estuarine numerical modeling. If a numerical model does not have conservative property, the model is unable to capture a front correctly and the total amount of conservative substances will not remain unchanged. The conventional
EulerianLagrangian numerical schemes do not have the conservative property and hence cannot be used for practical longterm simulations. The present model possesses the global conservative property because of the use of the quasiconservative EulerianLagrangian method.
The Coriolis terms in the momentum equations are discretized explicitly in the X direction and implicitly in the Y direction, as opposed to the fully explicit treatment in Casulli and Cheng. This method is conditionally stable but the time increment is much bigger (about 3 hours) than the other time step limits. If the Coriolis terms are not considered, the present model is composed of a large fivediagonal system of linear equations. The matrix is symmetric and positive definite. This fivediagonal system of equations is solved efficiently by the preconditioned conjugate gradient method.
xvii
If the Coriolis terms are considered, the ninediagonal
matrix is positive definite but not symmetric, and is solved by the successive overrelaxation (SOR) method.
In addition this study compares several advective schemes versus analytical solutions. "Wetting and drying" features, which are important for coastal modeling, have been
implemented in the model and compared with analytical solution as well. Several additional analytical tests were conducted to demonstrate the accuracy of the model.
After successful comparisons with analytical solutions, the model is applied to the winddriven circulation of Lake Okeechobee, Florida. The results demonstrate that the present
model represents the wetting and drying dynamics well and possesses the conservative property for water mass in Lake Okeechobee.
xviii
CHAPTER 1
INTRODUCTION
A large number of numerical models of circulation in coastal, ocean, estuaries, and lakes have been developed over the last few decades. The earlier models are generally onedimensional or twodimensional. With the advancement in computer technology, more recent models are threedimensional and timedependent.
One pioneering threedimensional coastal and estuarine model was developed by Leendertse (1967, 1975). The model uses a timeexplicit method to solve the threedimensional hydrostatic equation of motion in Cartesian coordinates. Sheng (1983) developed a threedimensional model of coastal and estuarine circulation in horizontally Cartesian and vertically stretched(a grid) coordinates. They solved the external mode(water level and vertically integrated velocities which govern the surface gravity waves) using a factorized implicit scheme and the internal mode(vertical structure of threedimensional velocities and temperature and salinity) using a verticallyimplicit scheme.
Blumberg and Mellor(1983, 1987) developed a three
2
dimensional estuarine and coastal Ocean circulation model using both Cartesian and orthogonal curvilinear coordinates in the horizontal plane and U grid in the vertical direction. They solved the external mode using an explicit method and the internal mode with a verticallyimplicit method. More recently they used an implicit scheme developed by Casulli(1992) to solve the external mode.
CH3D (Curvilinear Hydrodynamic Three Dimensional) model developed by Sheng(1986) is a stateofart finite difference numerical model which uses boundaryfitted curvilinear coordinates in the horizontal directions and U coordinate system in the vertical direction. The CH3D model also includes a sediment transport model and water quality model.
TRIM (Tidal Residual Intertidal Mudflat Model) was developed by Casulli and Cheng(1992). A Cartesian coordinate system is used in both horizontal and vertical directions. The model does not separate the external mode from the internal mode but solves them together with a timeimplicit method. However it uses a nonconservative EulerianLagrangian method(ELM) for the advection terms in the momentum equations and the linear advectiondiffusion equations. Hence for long term simulation, there will not be conservation of mass, momentum, heat, and salinity. In addition their model uses an explicit scheme for the Coriolis terms which could render the
3
solution unstable. The model allows wetting and drying in shallow waters.
1.1 Obiectives
The objectives of this study include the following
* Produce a robust threedimensional model by improving the
Casulli and Cheng(1992) model by implementing a conservative EulerianLagrangian Method(ELM) and adding
baroclinic terms and salinity equation,
9 Improve the EulerianLagrangian method by implementing an
explicit scheme in Xmomentum equation and an implicit
scheme in Ymomentum equation for the Coriolis terms,
0 Test the model with analytical solutions, and
* Compare the model results with those of CH3D model.
1.2 Organization of This Study
In Chapter 2, the governing threedimensional shallow equations which are momentum, continuity and salinity transport equations are written for a threedimensional Cartesian coordinates with boundary and initial conditions.
In Chapter 3, the governing shallow equations and the transport equation will be written by finite difference equation with boundary conditions. And a preconditioned conjugate gradient method to solve system of equations
In Chapter 4, the EulerianLagrangian method (ELM), or
4
semiLagrangian numerical method to solve linear and nonlinear advection equations will be discussed in detail. Model results are compared to analytical solution.
In Chapter 5, several analytical problems are developed and compared to model results. Wind test, tidal propagation, nonlinear term, Coriolis effect, bottom friction effect, diffusion test, baroclinic effect and wetting and drying test over tidal flats for both two and three dimensional cases.
In Chapter 6, the present model results will be compared to those of CH3D model results.
In Chapter 7, the present model will be applied to simulate winddriven circulation in Lake Okeechobee, Florida.
Chapter 8 will present a summary of the model and conclusions.
CHAPTER 2
GOVERNING EQUATIONS AND BOUNDARY CONDITIONS IN THREEDIMENSIONAL CARTESIAN COORDINATES
2.1 Introduction
In this chapter, the basic threedimensional equations
and boundary conditions for estuarine and coastal circulation are described with respect to Cartesian (x, y, z) coordinate system.
2.2 Governing Equations in Cartesian Coordinates.
The four basic assumptions used for deriving the
governing equations are : (1) the flow is incompressible so that the flow is nondivergent, (2) the vertical accelerations are negligible compared tc gravity, thus, the vertical pressure distributions satisfies hydrostatic approximation,
(3) density is approximated by its mean value except when multiplied by gravity, i.e., the Boussinesq approximation is valid, and (4) the eddy viscosity concept is used with constant turbulent eddy viscosity and diffusivity.
Using a righthanded Cartesian coordinate system (x, y, z) with the origin at the water surface and z measured upward, the equations are written as following nonconservative form (Lamb,1945) :
Continuity equation, du dv dw
+ + 0 (2.1)
dx dy dz
Xmomentum equation,
du d u du du 1 dp d 2u d 2u
X A(T 2)
dt dx dy dz po dx dx + y
d du
+ (AV )+ (2.2)
Ymomentum equation,
dv v + v dv 1 dp + 2v 2v
y+ u + v + wz + AH( + )
dt dx dy dz po dy dx2 2
d dv
+ z(Av)fu (2.3)
Zmomentum equation,
p= pg (2.4)
dz
Additionally, by using the Leibnitz integration rule, the hydrostatic pressure gradient can be split into the barotropic and the baroclinic components as follows
Vph = gp(xy,,t)V (xyt)+ gV p(xy,,t)d7
VPh = gp(X,Yliqt)Vq1(Xlylt)+ gJVp(x,,,t)1
(2.5)
Finally the momentum equations in x and y directions are
du du v+ wdy dz
g d +
fo d + po dx
d 2 2 2H 2T
A ( + 2 )
dv
dt
dv Ux+
dx
dv vd+ dy
g pd 2v d 2V Pd +AH(x2 +2 po dy x y
Salinity equation,
dS aS +u +
V S +wS = 2S d 2S v +wz =1(T + 2) dy dz dx2 y
Equation of state,
Where u(x,y,z,t), v(x,y,z,t), and w(x,y,z,t) are
the velocity
components in the horizontal x, y and vertical z direction, t is the time, q(x,y,t) is the free surface elevation, g is the
du
dt St
du U+
dx
8 du + (A, )+ fv
dz vdz
dv 89 w g y
(2.6)
d dv
+ (AV ) fu dz dz
(2.7)
d dS
+(D ) dz v dz
(2.8)
p = p(T,S)
(2.9)
dq
gx
8
gravitational acceleration and A,, DH and A., D. are the horizontal and vertical turbulent eddy coefficients, f is the Coriolis parameter defined as 2Q sino where Q is the rotational speed of the earth and 0 is the latitude, T is temperature, S is salinity, p0 is the density of freshwater, p is the density, p is pressure, respectively. Various forms of the equation of state can be used. Eckart (1958) used
p = P / (a + 0.698P)
P= 5890+ 38T 0.375T2 +3S (2.10)
a = 1779.5 + 11.25T 0.0745T2 (3.8 + 0.1OT)S
where T is in 'C, S is in psu and p is in g/cm3. CH3D uses equation (2.10), while Blumberg(1978) used p = po(1+ aS) (2.11)
where a is a function of temperature, and for a temperature of 15C, a is 7.57xi04. The value of a can be derived from U, tables (Bowden, 1967). In the present model, Blumberg's equation is used because the present model doesn't have a temperature equation.
2.3 Boundary Conditions in Cartesian Coordinates
The boundary conditions at the free surface are specified
w w
by the prescribed wind stresses Tx and Tw
2 2
pATz T x pCda Uw +Vw (2.12)
dv w 122
pAvl PCdaVw w +WVw (2.13)
Wz w
where Txw and Ty are the wind stresses in x and y directions
y
at the free surface, uw and vw are the components of wind speed measured at some distance above the free surface and C, is the drag coefficient. The drag coefficient is normally a function of the roughness of the sea surface and the wind speed at some height above the water surface. For this study, the empirical relationship developed by Garratt (1977) is used. Garratt defined the drag coefficient as a linear function of wind speed measured at 10 meters above the water surface, FT :
Cda = 0.001(0.75 + 0.067W) (2.14)
where W is the wind speed in meters per second.
At the free surface, the kinematic condition states that
Dij d) diq d) l)(215 w= +u +v (2.15)
Dt dt dx dy dt
Integrating the continuity equation over the depth and using the kinematic condition at the free surface leads to the following free surface equation :
+ vdz 0 (2.16)
dt dx dy h
where h(x,y) is the water depth measured from the undisturbed
water surface and H(x,y,t) is the total water depth, given by H(x,y,t) = h(x,y)+ ql(x,y,t).
The boundary conditions at the sedimentwater interface are given by specifying the bottom stress in a form of the ChezyManning formula
du b
pAVz = x yu (2.17)
Iv
d pAz Y (2.18)
ii
g U2 + V2
where u and v are velocities at the bottom
layer, and C2 is the Chezy friction coefficient which can be formulated as
C, = 4.64 (2.19)
n
where R is the hydraulic radius given in centimeters and n is Manning's n In shallow estuaries, the hydraulic radius can be approximated by the total depth. In CH3D model, the bottom stresses are calculated differently (rx, ,)= PCdb(112 + V2 (1,,V) (2.20)
where Cdb is the bottom friction coefficient, and tj and V1 represent the velocity components at the first grid point above the bottom. Taking z, as half of the bottom layer thickness (which starts at the bottom roughness height, z0), Cdb, for a hydraulic rough flow, is given by (Sheng, 1983)
Cdb = k In (2.21)
L
where k is the von Karman constant.
12
The open boundary condition for the surface elevation is prescribed by
q(x,y,t)= An cosr + 0 (2.22)
where An, Tn and On are the amplitude, period and phase angle of the tidal constituents. When open boundary conditions are given interms of the surface elevation, the normal velocity component is assumed to be of zero slope while the tangential velocity component may be either (1) zero, or (2) zero slope or (3) computed from the momentum equations. In the present study, it is assume that the velocity gradients are zero at the open boundary. At the fixed boundary, no flux through the boundary is allowed.
Salinity is prescribed as a function of time and depth at the open boundary or salinity is calculated from the one dimensional advection equation. In this study, the salinity at the open boundary is prescribed. At the water surface and the bottom the normal salinity flux vanishes, i.e.
as
D, = 0 (2.23)
If the estuary is connected to freshwater river flow then the boundary conditions at the interface are given as
S(x, y,z,t) = 0
(2.24)
13
u(x yz, t) = (2.25)
(h + q)b
where q is the river flow, usually a constant, (m3/s) and b is width of the corresponding river cells. At the open boundary the salinity is prescribed as S(x,y,z,t)= SO (2.26)
where So is the prescribed salinity value in psu.
CHAPTER 3
A THREEDIMENSIONAL NUMERICAL MODEL
3.1 Introduction
The governing equations and boundary conditions shown in Chapter 2 are discretized using finite difference method similar to that described in Casulli and Cheng (1992) A simple stability analysis of the twodimensional, vertically averaged shallow water equations showed that the celerity term lH in the equation affects the stability of the finite difference equations. Results of this analysis led to the development of an implicit method (Casulli and Cheng, 1992) and a semiimplicit method (Casulli and Cattani, 1994) which have been proven to be unconditionally stable.
When the implicit method is used, there is undesirable numerical damping in the numerical solution. To improve the accuracy of the implicit method, an implicitness parameter G is introduced such that for 0 = 1, this method is fully implicit and for 0 = 1/2, this method is semiimplicit and remains unconditionally stable and doesn't have undesirable numerical wave damping.
15
In this chapter, the threedimensional semiimplicit circulation model is developed by using the so called 8scheme.
3.2 Finite difference form of momentum ecquations
The horizontal gradients of surface elevation in the momentum equations (2.6), and (2.7) and the velocity gradients in the free surface equation (2.16) will be discretized by the escheme, and the vertical mixing terms will be discretized implicitly because the vertical length scale is much smaller than the horizontal one.
Figure 3.1 Schematic diagram of stencil
As shown in Figure 3.1, the spatial mesh consists of rectangular boxes of length Ax, width Ay and height Azk,
16
i.e., a cellcentered grid system is used. Each box is numbered at its center with indices i, j, and k. The discrete u velocity is defined at half integer i and integers j and k; velocity v is defined at integers i, k and half integer j; vertical velocity w is defined at integers i, j and half integer k. Water surface elevation q is defined at the center of each column i and j. The density of water is defined at the center of each box i, j and k. The water depth h(x,y) is specified at the u and v horizontal points.
Then, parametrized semiimplicit discretization for the momentum equations (2.6) (2.7) and for the free surface equation (2.15) takes the following form.
n+ j A t n+1 n+l in ,j in
g At L Az (p",j,, pi, ) Az Pi, ,m Pi ,m
poAX m= k iy2, j i t+ ,j 2
n+1 n+1 Un+1 n+1
A i+ '2,j,k+1 i+yjk i+ j,k i+yj,k
Vk+J4 A k)E AZ+k2
A V k + 2 i + j ,k + Y 2A 2 A i + j ,k Y
+At
(+ ,j,k
vn
+ f,+ j,k
(3.1)
n+1 = nAt [ n ++l )+qn ]
i,+Y,k i,= Fvj+Y2,k I ,+ 0 )X ,";+ _,")
A t n pn.
M '
 g t Az"2 (p"2i p" ,,) Az
o 7 m= t~gm +1, ij~r A ;gm
n+1 n+1
vv
i,j+ Y2,k+1 i,j+ 2,k
Ak+Y Az ij+ ,k+Y _
+At
 Un+k
i, j,+ ,k
n n
ij+,m Pi,jm
2
n+1 n+1
v v
i,j+ 2,k i,j+ 2,k1
Si,j,+ 2k12
Azl,+>k
(3.2)
n+1 n t M n+1 M n+ 1
 kn ++ ,n k +1 ,k k= z Un j'i
7, = rl Ax 0.L m A ~jk i+Ak , i), ,Ii2,Ikj
AX [k=n i kiM 2
 t(6 Az. u"kV / Az u" ki~/,
At L M n+1
A I Z Az+)~12,k tij+ 12k Zijy, A11 Y2I
Ay t[k=,,I ,I/2 2, k=m
A t (1 Az "
A = Az ij+ y2,k i+y,jk k=m ij 2,jk i,jk
(10 k=in
At n
(1 0)[ A z. ,ki+ Ui,k I Az, , it" ]
Ay /___ki 2'k i' 2
(3.3)
where m and M denote the limit of kindex representing the bottom and the top finite difference stencil, respectively. In equation (3.1) and (3.2), Fu and Fv are explicit, nonlinear
18
finite difference operators which include the contribution arising from the momentum advection terms which are discretized using EulerianLagrangian method and the horizontal diffusion terms.
Here, it is important to discuss the treatment of Coriolis terms in the equations. The Coriolis term in the Xmomentum equation is discretized explicitly, while the term in the Ymomentum equation is discretized implicitly (Sheng, 1983). Although the Coriolis term in the Ymomentum equation is discretized implicitly, the actual computation of the term is explicit because the velocity Un+1 is already obtained from the Xmomentum equation. Stability analysis of this method shows that the time step of the scheme should be less than 1 which is about 3 hours and much bigger than the other
time step limits, e.g. the horizontal diffusion limit. Since the numerical grid is staggered, v is not defined at the u grid point and u is not defined at the v grid point. A simple way to define the undefined velocity at a grid point is to take the average of the velocities at four neighboring
points.
The fictitious values of u and v above the free surface and below the bottom in (3.1) and (3.2) can be eliminated by using the surface and bottom boundary conditions as shown in
(2.17) and (2.18). The difference forms of (2.17) and (2.18) are :
A M+ A i.%.+Y2 X
A V .'j ,M
n+l un+l
ui+Y2,j,M+1 u+ 2,j,M
Az.
i+ 2,j,M+12
n+l n+l
V ,J +12M Y2ij= y.
AzijM+ 1
n+1 n+1
SJm u ym1 n+ / n+1l
Az i+ 2,j+,,m 2 i+2,,m i+ 2,,m
(3.4) (3.5)
(3.6)
(3.7)
n+1 n+1
iY +2,m ij+ = n+2, n+1
Avidj+gm ij+ 2,m i,j+ 2,m
i,1+ 2,m 2
n+2 = (n+ 12 n+ ( 12
n+2 n At n
i+ 2,j,m i+ 2,j,n 2Ax qi+l,i
n+g2 Vn At ( .
i,j+ /2,m i,j+ 2,m 2Ay "ij
and
where
and
(3.8)
(3.9.)
(3.10)
qij
By using the boundary conditions (3.4) through (3.7), the fictitious velocities in equation (3.1) through (3.3) will be replaced by values of u and v defined within the computational domain.
3.3 Solution algorithm
To solve equations (3.1), (3.2) and (3.3) computationally, the equations are first written in matrix form as
At r ,+
A" U"+' = G" gO(LJ + ")AZ.+ + AtfV (3
i+ ,j i+ ,j ,i+ ,j AXj i ,j i+/,j f +2/,j
A "+ = G"n A ( n+I )AZ fUn+ (3.12)
.".= G+ [ (Z U AZ" ) U+
.1j J A i+/, i+Y2,j i ,Y
At
 AZ" V. 1 A (jy V" (3.13)
we U Ay id i aj jwhere U, V, AZ, G, S AndA are defined as follows
Un+I ,n+I Un+I ui+,/nl
Sn+1 i,j+Y2, M
Vn+1 iQ+Y2,M1
SAZ=
n+1
, n+1 i,+Yr ,m
AzM AzM1
.Azm )
A ,t 0 ) ,, ]
AzM Fut" +M A) + C" + BC" + At
F At 1
Az Fu," g x(1 )( i I' ,)+ C" +BC".M
[ I At ./ MI]
AZM_1 [g.,+ ...,~~ g~ o IV27+., ,J)+c ,M..,+ l + .1.+Y g/
Az..i[ FU" Az,,, Fu"'
I Az 2,jm
Az,, [Fit",4
At
At 1 ") + C" + BC, 1 Ax +Y
At
Ax g )(," ,) + C+B
At
AzMI FVI, g (1 O)(q,"j+, qi,j) + BC M+ AIt t
A, 1,j+Y2 m
At
Az,,f1 Fvn g(1 ) i) + BC"
[ At
Az,+, Fvn + g 0 ( 0)(qfi+1 G i)+ BC m+
[~ At
Azm Fv",; g0(10)(q";+1 i) +BC', + m
where C's are Coriolis terms in the Xmomentum equations
(3.14)
(3.15)
(3.16)
vn+
i,j+ Y2=
G"Y2 =
Gn
I, =
Cn1, =
AtfV+Y,j,M AtfV ,M
A v ,", +,.,, l Ay te+ jm+
and BC's are the baroclinic terms
Az (p' p,, ) Az
i+A,j,k i+1,j,k ,k i+ ,jM
At
 Z Y2 (PinIjk P1,jk) j
At z (pil P,,J) Az" 2
 o k= M1 i2 ,jk i+I,j,k j,k) i+)j, M1
PA k= +1 + ,2~ j, k n~ +2, jlm+1
AtFm
gT
PoXk=.
nAz+ ki k Pnjk AZn
Az+y2,j,k (Pi,+l,,,k ij,k) AZ+y2,j,,.
n1
Pi+I,j,M Pi,j,M
2
n n
Plj, Pj,MI
2
n n 1
Pi+l,j,m+l Pi,j,m+l
2 n
,l n.
Pi+l,j,m Pi,],n
2
(3.18)
(3.17)
AtF
g[ k=MA
BC"
i+Y2,jlk
At [
 o k
p0AY kA
AZij+ ,k Pi,j+l,k i,j,k) ij+_,nM
g P'At [k=Mm+ A i,j+2,k(Pi,j+l,kAt m
g 1m Az.J"k (Pj+lk
poAy L k=M i j+ ,k inj+1,k
At[ M
oA L k i+1 j+ ,k j+1,k gmA j+ k( i j+1,k
ij,k) A+2.P,,j,k) Azj+,+ i j,k j+,m
n np
Pij+1,M Pi,,M
2
n pm
PiJ+IM1 Pi,j,M2
,n ii
9ij+1,m+1  Pi.j,n+1
2
n n
Piylm PijmI
2
(3.19)
)[ Z T+ l2,j TU ,
itj i+ 2,j
At F
(1 ) Az )
Ayij+
 (Az U"j T ]
a ,j,j
T
V" Az,
i,j+ 1,1 2
SAzAI + a,_ a Y
 a f AzA, + amY + aMY
 a My
 am+y Az, + am+y + rAt)
(3.21)
k At
with ak = Azk
BC =
i, j+Y2,k
At
3". = (11,j by A
and
(3.20)
V"
1.1Y
Formal substitution of Equations (3.1) and (3.2)into Equation (3.3) yields :
qn+1 (n
,fj ,j
At [AAZTG g A 8(i,7 I 71+ )]A'Az} ]
S AzTY AG} g (1 ,n+l A Zn
A ~ ~ g__ yn A [ O(77i,.+,_ _,,,I A 1
 Iz A g A (n J ]) A 2
 At [[A {A 'G}t g 0( ,+, ,,l) {A Az}i+ ,j
+ L j+{A2 AG} g [ ( +",1 4,j+ )]A A'z 2
+ AI { G} g [(I +lj q+l) A AAz
+ i G g ( "j )]A i,j, { AzY2
At rO ,.+1 ++I n+1Ii
 jzy 1 AlG g ( ) Az
 [A i {A G}n g ( i +1 )]A A i2,j+
+ AAt [AG g (n 1+ 1 A'Az1
I+ A7 {A G}i2 g)[e(l, l, A jJ)] A i+'Az ,
+ A A'G g [( n+.1 +1 1 'AAz }
(3.22)
25
Neglecting the Coriolis terms and rearranging the Equation (3.22) by moving "n+l" terms to the lefthandside, Equation (3.22) becomes
, n+1 g 2 A[ [(AZ)A Z]i (iq, n+l [(A ) ),_,j~
 (Az A'Az]i (n,"i
 y At 2 2TA IAz (,j ', ',j )
Ayj + 1(3.23) [(Az) A Az.. Qnl, + ;,)]
 Ii Y2 i ( lf qli,j .i,.=" At 0 [(A Z) AG] [(Az)A'GA ]
A x Li+ ,j i ,j
At 0~ [[(Az)T A'G.y (Az)TA1G] ]
Equation (3.22) can be rewritten into a compact form :
d ,"f+' S n+1 S n+1 S .n+1 S qn+l n
i j i Y2,j i+,j iY2,jilj iJ+ 2i,J+ iJ 2iJ1 = qi,j
(3.24)
where
S" At2 n T 1 n
",= g Ax [(Az)TA 'AzYj
12n
Sn = g 9 2[(Az)T AzA+n
i, 2
d, = I+ S"n +S +S"n +n S(3.25)
j +Y ,j Y ,j i+ +Y (3 .25
q". = 8.n. At [[(Az)T AG] (Az) A'G ]
Ax ,i+ ,J i2,J
At 0[(Az) AG y [(Az)TA'G],
Ay Q+ i j 2
26
If Coriolis terms are considered, Equation (3.24) becomes more complicated :
C +n 1 1j1n l j n l in1n+1 n n l
G in1,I q q,Ij q i
7nl +l,~ j(3.18)
_, ~ C ",7 ~ _ n q n + l  ( n n + 1 l? ~ S 1 ni+l,j1 8 i+l,j +,j+1 = RHS"
where
C" = FIPn P"n
Q i,2 iY2,j1
C" = CM" + FIP" P" FIP" P1
2 iYj i,j+ i Y2,j ij2 i2,j
C=FIpn +pn"
3 i, j+ 2 i 2,j+)
C= HMn + FIP" P" FIPn" 2Pn
4 i, jY i, jY i+2, j1 ijY i2, j1/
C = 1+ CMn + CM" +HM" +HM"
C5 L i+Y2,j iY2,j+H idj+Y2+ i,jY
 FIP Pn FIPn pn FIP Pn FIP P1
i,j+Y2 i+2,j i,j+Y2 i2,j ij2 i+ Yj ij2 i2,j]
6 1,j+ i,j+Y2 i+y2, j+1 ij+ y i2,j+1
C7 = FIP" P"
7 i,j 2 i+ y2,J 1
C"(= CM" + FIP P FIP" P"
C" = FIP +P" .
\i, j+)/ it y, j+1
(3.27)
with
At2 At2 2 gA2
C gx2 H= g 2 2
SAy AxAy (3.28)
P= (AzA'AZ), M= (AzA'AZ), F= At
4
and
RHS" = 6," 6I 0(Az) A'G]". [(AZ) AG]i
lAX k 2J "2
 At [(Az) A'G (Az)T A_G]n]
Ay ij+ i,
+ {A A{:zTA} j+j 'G}>;+, + +{AzA} j+ { AG}1,+ ]
+ Atf At O[{AzTA}1 y {A'G} + {AzTA}' 2{A_'G},J
4 Ay
+{AzTAI'A A G,_ + {Az A} { A'G} I
 ~Q 8AzA)_A 2G AzA 1, AY2(
+A zAZ _ A 'Gn +AzT A'G1
(3.29)
Since the coefficient matrix A is symmetric and positive definite, A' is also positive definite, and therefore
(Az)TAIAz is a nonnegative number. Moreover the coefficient matrix A is a tridiagonal matrix, and the inversion of the matrix A is straightforward and efficient by using Thomas algorithm. Finally Equation (3.23) constitutes a linear five
28
diagonal system of equations for new water surface elevation
' and the associated matrix is symmetric and strictly diagonally dominant with positive elements on the main diagonal and negative ones elsewhere. Thus the system is positive definite and has a unique solution. If the Coriolis terms are considered, Equation (3.26) consists of a ninediagonal system of equation, i.e., a block tridiagonal matrix. The ninediagonal matrix is positive definite but not symmetric. The fivediagonal system of equations can be solved very efficiently by a preconditioned conjugate gradient method. For the ninediagonal matrix case, the block tridiagonal algorithm or successive overrelaxation(SOR) method can be applied. Present model uses SOR method. After obtaining the new water surface elevation, the new horizontal velocities U and v can be calculated by using Equation (3.11) and (3.12).
Finally, from the continuity equation, the new vertical velocity w can be obtained by setting wn1 =0 and
= Ax
,n+1 n+1 1+ ,j,k i+ xjk i j,k i ,j,k
i,j,k+y ij,kY2A
Azn v AZn Vn
i,j+ ,k i, j+ y,k i,j(,k i,j3,k
Ay (3.27)
k= m,m+ 1,...,M
3.4 Conlugate gradient method
Equation (3.24) constitutes a very large system of equations, which regime a large amount of computing time for its solution. Since the system of equations (3.24) has
symmetric and positive definite sparse matrix, the equation can be solved efficiently by the i.e., conjugate gradient method.
Conjugate gradient method ( Beckman, 1959. Casulli, 1992) is the most papular iterative method for solving large system of equation Ax= b where the matrix A is square, symmetric, and positive definite (or positive indefinite).
Consider a quadratic function of a vector with the form
I
f(x)= xT AxbTx+ c (3.28)
2
where A is a matrix, x and b are vectors, and c is a scalar constant. If the matrix A is symmetric and positive definite, f(x) is minimized by the solution to Ax= b The definition
of positive definite matrix is that if xTAx>O, for every nonzero vector x, then the matrix A is said to be positive definite matrix.
The gradient of a quadratic form (J. S. Ahewchuk, 1994) is defined to be
a
d f (x)
d
f'(x) Ax2f (x) (3.29)
d f (x)
With equation (3.29), the gradient of equation(3.28) is
1 1
f'(x) = ATX+ I Ax b (3.30)
2 2
If A is symmetric, this equation reduces to f'(x) = Axb (3.31)
Setting the gradient to zero, we obtain Ax =b which we want to solve. Therefore, the solution to Ax=b is a critical point of the function f(x) If A is positive definite as well as symmetric, this solution is a minimum of f(x), so
Ax = b can be solved by finding an x that minimizes f(x). It should be noted that if A is not symmetric, the conjugate gradient method finds the solution to the system (J. S. Ahewchuk, 1994)
](AT+ A)x= b (3.32)
2
31
To solve Equation (3.24) using the conjugate gradient method, the equation can be written in the normalized form (Casulli and Cheng, 1992) :
sn Sn
n+I i+12,] an n+1 S 2,J n n+l
" d d+i,j
S"1 Sn a
Sd d
Ili 1+111 "'J 2, I n n+ _' .+l S 2 d n n+1 qi
Id:d:+ aw'J' + d ,ad,,:_, ,j~
,Jn nliJI n j n+ 1
which, by changing the variable 1 to e,j is equivalent to
ij a+ e+l, a e,_,j a ye,,+l a,,j e,,j_ = b, (3.34) where
Sn Sn
'2,,J ai,j+2 qi,,
a i+ij /iji+l ia,j+ n d~d~+ bi'j n,
ityd d" d Vd d"jd
Now the equation (3.34) is applied to the conjugate gradient method to obtain the unknown value of ei.j The
conjugate gradient algorithm to solve the system of equations (3.26) takes the following steps.
(1) Guess (0)
i,J
(2) Set
(0) (0) e(0) (0) (0)
ij ij ij i+ 2,j i+ ,j 2,j 1,j
(0) (0) b
 a i,+Y ,j+1 a1,j_ i,j ij
(3.35)
(3) Then for k=0,1,2, .. and until (r(k),r(k)) < calculate
()(k) (k)
where a (k) = ((),Mp(k))
(kt), Mp(k))
(k+1) (k) (k) Mp(k)ij
r.. ri. (k (li ,()..
p(k+.) (k+1) + p(k)() where (k) ij =j + f ij h
((k+1) (k+1))
(k),r(k)
In equation (3.37) each element of the vector Mp which is preconditioner is simply given by
Mp k.) = p(k) ak)a (k)
1,j = "i, ai+,l2,]l"Pi+l,] i 'l2,j pil, ..
3.3 9)
_ (k) ( a () i,j+ 21i,j+1 j 2 i,j1
3.5 Salinity Scheme
The finite difference equation (2.8) is similar to the
e(k+1) (k) (k) (k)
i,yj = i,j a( Pi,]
(3.36)
(3.37)
(3.38)
33
hydrodynamic finite difference equation. The vertical diffusion term is treated implicitly, while the horizontal diffusion and advection terms are treated explicitly. The advection term will be solved by the conservative ELM method.
The finite difference equation for the salinity equation (2.8) has the form as
.n+ At DSijk+1 i,j,k ijk i,,k1
i,j,k AZ vij,k2+ AZ i,j,k2 Azn k
Ij,k i,j,k+ ijk=FS,"jk + DHAzik S1,,k 2S,jk + Si+ljk (3.40)
D+ t Azi k ( . )
DHt \i,I1,kn 2S~ik + S Ilk
+2 i jI,k 2 i,j,k i, j+I,k
where FSi, k is an explicit discretization of the substantial
derivatives of convective term by ELM. The left hand side of Equation (3.40) constitutes a tridiagonal matrix and the solution procedure for the unknown new value of S, is
exactly the same as that for the momentum equation (2.3).
CHAPTER 4
EULERIANLAGRANGIAN (SEMILAGRANGIAN) METHOD
4.1 Introduction
One major challenge in numerical modeling of free surface flows is the accurate and efficient treatment of the advection /convection terms in the governing equations. There are numerous numerical schemes for solving the advection/ convection terms. In this chapter, numerical treatment of the advection/convection equations by the EulerianLagrangian Method (ELM) will be discussed. Other methods such as upwind, Leith's method and total variation diminishing (TVD) method are briefly presented in Appendix A.
4.2 Review of EulerianLagrancian Method (ELM)
The EulerianLagrangian Method uses the Lagrangian form of governing equations in an Eulerian computational grid system. The method (Starniforth and Cote, 1991) is generally referred to as semiLagrangian method in the community of numerical weather prediction, but is called EulerianLagrangian Method (ELM) by estuarine and coastal modelers,
e.g., Baptista (1987), Casulli and Cheng (1992).
The ELM solves the equation for a variable e.g., a velocity component or concentration by following a fluid particle along the characteristic backward in time and then
interpolating at its foot (departure point of the particle) to obtain the values of the variables.
Figure 4.1 Characteristic lines and their feet
The ELM was also used by Tang and Adams (1995) with finite difference method and by Baptista (1984) with finite element method.
The major advantage of the ELM is that the method is unconditionally stable although it is explicit. In addition,
it is transportive and accurate (Casulli, 1994) Despite these advantages there are several disadvantages of this method. Although the method is accurate, its accuracy is dependent on
the interpolation function used. For instance, if a linear interpolation function (first order accuracy) is used, there is excessive numerical diffusion in the numerical solution,
with the maximum diffusion on the order of Ax2/(8At) in the ID case, although the diffusion can be controlled by reducing
the spatial increment or by increasing the time step size. The second problem is that if higher order (second order of higher) interpolation is used the numerical solution contains unphysical numerical oscillations. Godunov (1959) was the first to show that any approximation methods of order higher
than one is unable to deliver a solution free of oscillations. This result, which is known as Godunov's theorem, is a major
stumbling block for the development of integration schemes for hyperbolic problems.
The most significant problem of ELM is that the scheme
doesn't have the conservative property which is that momentum, mass, and any species should be conserved locally and globally throughout the numerical simulation period. Roache (1972) indicated that all methods of characteristics (including ELM) which solve the nonconservative form of momentum or transport equations are nonconservative. If the governing equation is written in conservative form then many numerical schemes can have conservative property.
There are several ELM methods which yield solutions free of oscillations and satisfying conservative property
(Priestley 1993 and 1994, Gravel and Staniforth 1994, Rasch 1994, Tang and Adams 1995, Leslie and Purser 1995, Lin and
37
Rood 1996 and 1997). Although the methods have different names the basic idea is the same. To control numerical oscillations these methods enforce the solution to be "monotone" i.e. free of oscillations. To achieve the conservative property, the conservative ELM uses FCT (Flux Corrected Transport)like method and checks for global conservation at every time step.
4.3 Derivation of EulerianLa rangian Method (ELM)
Consider the diffusionfree nonconservative advection equation
d+ uVc= 0 (4.1)
dt
where the c(x,t) represents any scalar quantities, u is
velocity field, and V is a gradient operator. Most of the fundamental equations in fluid dynamics can be derived from first principles in either Lagrangian form or Eulerian form. Lagrangian equations describe the evolution of the flow that would be observed following the motion of an individual parcel of fluid. While Eulerian equations describe the evolution that would be observed at a fixed point in space. The onedimensional Eulerian advection equation (4.1) can be expressed in Lagrangian form as
Dc
 =0 (4.2)
Dt
The mathematical equivalence of Equations (4.1) and (4.2)
follows from the definition of total derivative,
D 8 dx8 Dt dt dt dx
(4.3)
and the definition of the velocity,
dx
 U
(4.4)
When a Lagrangian numerical treatment is applied to Equation (4.2) the computational grid will be continuously deforming in the general case when u is nonconstant. For operational advantage, Equation (4.2) will be discretized on a fixed Eulerian grid system. Figure 4.2 shows the Eulerian Lagrangian mesh system.
#
Figure 4.2 An Eulerian Lagrangian Mesh
A finite difference scheme for Equation (4.2) is simply
Ia in
uAt
39
n+1= cfl (4.5)
where C~a = c[(ia)Ax,nAt] and n is previous time step, i is the index at a grid point and a is the CFL number. In general the CFL number is not an integer, therefore (i a ) is not the index of a grid point and a proper interpolation formula must be used to define Cina" The stability, numerical diffusion and unphysical oscillations of (4.5) depend on the interpolation formula chosen. If a linear interpolation between (ii) and
(i) is used to estimate Cia' we obtain the first order upwind scheme. If a quadratic polynomial fit is used to interpolate between (i1),(i), and (i+1), one obtains the Leith's method.
The ELM uses a generalization of the interpolation concept of cina between two or more grid points which do not necessarily include the point (i) Consider first the case that C7a is taken to be a linear interpolation using one node upstream and one downstream. For a given a 0 let n be the integer part of a and p the decimal part, then a = n+p,
0 p<1. In this case Equation (4.5) becomes
40
k+1 k k k
C. C p( c__) (4.6)
or, equivalently
Cik+1 I p)in + PCik 1 (4. 7)
Note that if a<1, then n= 0, p=a and finite difference equation of Equation (4.6) and Equation(4.7) reduce to the first order upwind method, since 0p
Since the velocity u is generally nonuniform, the correct value of a can be found from the solution of the ordinary differential equation (4.4) using any backward trajectory computation. The velocity u is known only at time level tn, hence it will be assumed to be constant over one time step. Then, at each grid point (i), Equation (4.4) will be integrated numerically from t" to tn+. In this study the backward Euler method is used for computation of the trajectory. To compute the velocity, the time step At is divided into N equal increments, T At/N, and Equation (4.4) is discretized backward as
xS1 = xs_ ruk (xs), XN = Xi
s NN, (4.8)
s= N,N 1,N 2,...,2,1
where uk(xs) is interpolated with any interpolation formula.
x. x
Then, at x, a can defined by aAx
The finite difference equation (4.6) possesses the minmax property (Casulli, 1987), thus this method is entirely free of oscillations. Moreover the minmax principle also implies stability, hence the ELM shown in (4.20), though explicit, is unconditionally stable. Rigorous stability analysis can be performed using Von Neumann method. Assuming a solution to equation (4.1) as k A e (4.9)
where A f is the amplification factor, I is the grid point index and i=". Substituting Equation (4.9) into (4.1)
Ak = [1 a(I ei(k'))]ei(knAx) (4.10)
The magnitude of Equation (4.10) is
IAk12 = 1 2a(1 a)(1 coskAx) (4.11)
If the grid CFL number satisfies 0 a 1, the scheme is stable. Because the departure point lies within the
interpolation interval (in, ini), the magnitude of (4.11) is always less than one, consequently the scheme is
unconditionally stable (Starniforth and Cote, 1991).
To examine the artificial diffusion introduced by Equation(4.6), each term in (4.6) is expanded into Taylor series to produce the following :
Dc AX2 82C
S p(p 1)T+ HOT (4.12)
Dt 2Az'
where HOT stands for higher order term. Since p<1, the least upper bound for the artificial
Ax2
diffusion coefficient is which is the same as the upwind
8At
method. However, since the ELM is unconditionally stable and the value of At can be arbitrarily large, the maximum numerical diffusion can be controlled either by reducing the spatial increment or by increasing the time step.
The ELM with linear interpolation function has no spurious numerical oscillation and little numerical diffusion, but it is only first order accurate. If three point
interpolation function is used, higher order ELM can be obtained.
Let a=n+p, where n is integer, IP
this case Equation (4.5) becomes (Casulli, 1987)
k+1 i(p + p)c,+1 + (1 p2) 1 2
Ci 2 ( +A ' f2
Note that, if jai: 1 then (4.13) reduces to the Leith's
method and it is of second order accuracy. However, when p # 0 one gets a negative coefficient on the right hand side of Equation (4.13), consequently this method does not posses the minmax property and hence it will introduce spurious oscillation.
4.4 Conservative EulerianLarangian Method (ELM)
EulerianLagrangian Methods described above do not possess the conservation property because the finite difference equations are obtained from the nonconservative form of the differential equation. To achieve conservation, some special treatments are needed. Several conservative ELMs were introduced in the past decade, inducing Priestley (1993,1994), Gravel and Starniforth (1992,1994), Rasch (1994), Tang and Adams (1995), Leslie and Purser (1995), and Lin and Rood (1996,1997). The conservative ELM presented here is based on the QuasiMonotone SemiLagrangian (QMSL) scheme (Bermejo and Staniforth, 1992) and QuasiConservative SemiLagrangian (QCSL) scheme (Priestley 1992, Garcia and Priestley 1994).
The conservative ELM consists of two steps. The first step obtains a low order solution and a monotone higher order
44
solution, and the second step combines the two solutions. As seen in the previous section, the higher order solution does not have numerical diffusion but has numerical oscillation. The low order solution is a monotone scheme.
Consider the scalar, linear advection equation for u(x,t) u+ aVu= 0 (4.14)
To obtain a higher order monotone numerical scheme, the FluxCorrectedTransport (FCT) method (Boris and Book 1973, Zalesak 1979) is used. Following the FCT method, the solution is u n+I Tn+I pn+1 4 5
L,k + F k(4.15)
where U+I is the value of the approximate solution at the
grid point xk at tn+l i.e., foot or departure point of the
particle which is obtained by Euler backward method and L k
is the low order approximation at the departure point at t, and Fn+1 is a higher order correction term such that it approximates the exact solution with an order higher than that of the low order solution. The higher order approximation is any solution which is higher than order of one. The higher order solution used here is a quadratic second order solution.
The overall procedure is :
1) Obtain a monotone low order approximate solution U',,
l'n+l
2) Compute a higher order approximate solution UH,k by a standard ELM of order higher than 1, 3) At any instance t"+l set
Fk nl n+1 (U U )n+ ' k (t"H,k "L,k
(4.16)
where (k is a grid point coefficient that controls the amount of correction to be added to the low order approximation in order for the numerical solution to be monotone. The procedure of determining the coefficient which is called a limiting procedure is as follows :
4.1) Define U+ as the maximum U value of the
neighboring points i.e.
U = max(U ',, Uk2, U3, Uk4)
(4.17)
and U as the minimum U value of the neighboring points
i.e.
U = min(Ul, Uk2, U;3, U;4 )
(4.18)
4.2) Set
+ = U+ Un+
Q L,k
Q= UUn+ (4.19)
P = U n+ U n+
"H,k "L,k
4.3) If P>O,
C = min(1, )
If P<0,
k = min(1, (4.20)
'P
If P= 0,
U n+1 T]n+l qn+
U = Ui ,i. e, = 0
Hk UL,k 'k 0
Once the low order solution and the monotone higher order solution are obtained, the next step is to recover a
conservative solution. The final quasiconservative solution can be written as
U = akUHk + (1 ak)ULk (4.21)
where ak (0 ak < 1) is a constant to be determined by the recovering process.
The procedure to recover conservation is to compare local flux with the overall flux of the quantity which is related to ak. The procedure to define ak is as follows (Priestley, 1993
and 1994) 1) Define
f U(x)dx = Uodx = C
, = (UH,k UL,k)Sk (4.22)
where Sk is the area associated with node k 2) Compute
( u s U 1usc (4.23)
Zak kHUk)Sk UL,kSk=C (4.23)
Assume that
kmax (4.24)
If this is not the case, then the definitions are changed to Ak A (4.25)
C* C*
3) determine 'k
Step 1
If <0, then a nax. Otherwise k=
If A:5 0 then a,= a i .Otherwise a k= 0
Step 2 : Define surplus
SURPLUS = C* Zakk
Step 3 : Define the average value of a as :
SURPL US
aa vg
Step 4
max max
If aavg< akm then ak :avg Otherwise ak =ak
The lower order interpolation functions include linear(1D), bilinear(2D), and trilinear(3D) interpolations. The higher order interpolation functions include quadratic interpolation function(lD), biquadratic(2D), and triquadratic(3D) interpolation functions. The lower and higher order interpolation function for ID, 2D, and 3D cases are : 1) One dimensional case Linear interpolation : 2 points used
U'= (1 p)U_ + pU1_,,_ (4.26)
Quadratic Interpolation : 3 points used
U'+ = 0.5(p2 + p)Ui_,_z + (1 p2)U1_, + 0.5(p2 p)U_, (4.27) 2) Two dimensional Case Bilinear Interpolation : 4 points used
U!.' = (1 p)[(1 q)Utn,jm + qUtn,jm_1]
+ p[(1q)Ulnl,jm + qUn1,_m] (4 28)
Biquadratic Interpolation : 9 points used
U t+U
I'
(1+ p)[( (1+ q) n,jm + (1 q )Ui,l (1 q)Ui n ,j"z+l]
+(1 p2)[ (1+ q)Un, + (1 q)U,' ,jm q) in,jm+ ]
2_injm 2
_ (1 p)[(2(1 + q)U_ in+,,nI + (1q 2) ,,1+,,,, ( q) U,+,, ,,, +l]
(4.29)
3) Three dimensional case TriLinear Interpolation : 8 points used
U" =
i,j,k
(1 r){(1 p)[(l q)Utl,jm,kn + qUt1,jm1,kn
+ p[(1 q) Ull1,jm,kn + q U11,jm1,kn }
+r{(1 p)[(1 q)Ul,jm,kn1 + qUl1,jm1,kn1
+ p[(1 q)U,1,jm,kn1+ qUlljm1l,kn11l]} (4.30)
TriQuadratic Interpolation : 27 points used
UI = (p2 + p)(q2 + q)@
11
,j,k 4
2 2U~~~~~~~_ 1 2I
2 r 2 r U,'11, m1,n1" Ui11, jm1,kn 2 r 2F r Ui'11,jm1,kn+1}
1
+ 2 +(p +p)(1q
1 1
(r2 + + U(11jm,kn1 r2)U l 1,jm,kn + 2 r) U'I,jm,kn+I
2 2
1
+( + p)(q q).
4
1 1
2) _I,jm+I,,_n + 2( U ,jj ,n+,k+I}
{2(2 r)Ui',jm+,kn + (1 r2)Uil1m+kn r)U
1
+ (1 p2)(q + q)
2
1 1
2 (r2 + r)U11,k + ( 2)U, 1,k,,_ + (r2 r)Ul,j,,,,k,+
2 2
+ (1 p2)(1 q2).
1 1
2 (r2 + r)Uil,jm,kn1 2)Uil,jm,kn + I(2 r)Uitl,jm,kn+1I
+ r)Ugl,,,,kfl + (1 2U
2 in,kn 2~ )~~lklI
1
+ (1 p2)(q2 q).
2
1 1
{(r2 + r)Uijm+knl_ + (1 r2 U_ + 2(r + r)Ui,jm+,kn1 2 U,'_1,jm+1,kn 2 2 r)Uil,jm+I,kn+1
1
+ (p2 p)(q + q).
4
lr 1
2(r2 +rUl+1,jm1,kn1+(I r 2)Ul+1,jm1,kn + 2 r)U l+1,jm1,kn+1
1
+ (p2
+ (p2 p)(1 q2)
2
1 1
(2 2+ r)U+,j,k 2)Uil+1,jmI,kn + I(2 r)U l+,jm,kn+1
2+ 2(p2(2
2
2 rUl+ +,k1 + 2)Ul+, jm+1,kn 2 )Ul+,jm+,kn+
(4.31)
where U can be any quantity such as velocity, salinity or
concentration.
The ELM can be applied to the nonlinear momentum equation
51
with nonlinear advection. The procedure is exactly the same as that for the linear advection equation.
4.5 Numerical Experiments
4.5.1 Linear Advection Case
To test the various algorithms which are described above and Appendix A, various numerical experiments were performed. The computational domain is 1000 units in length divided into 1000 computational cells. The initial and boundary conditions are :
Initial condition
c( 0 < x < 50, 0) = 0.0
c(50 < x < 150, 0) = 1.0 (4.32)
c(x < 150, 0) = 0.0 that is a square wave. Boundary condition :
c(0,t) = 0.0
c(l,t) = 0.0 (4.33)
The advection velocity u = 1.0. The analytical solution to Equation (4.1) is given by c(x,t + At)= c(x,t) (4.34)
All the numerical experiments were performed by setting the CFL number to 0.5 because all the schemes discussed in
52
Appendix A except the ELM are unstable when the CFL number is greater than 1.0. In the ELM test, the CFL number was set to 5.5. The numerical solutions are compared to the analytical solution at time T = 50, 250, and 500 seconds.
Figures 4.3, 4.4, and 4.5 show the analytical solution to Equation (4.1) when time T = 50, 250 and 500 sec.
Figure 4.3 Analytical solution of linear advection
Figure 4.4 Analytical solution of linear advection
eqato atT 25 e
2o7
Figure 4.5 Analytical solution of linear advection
equation at T = 500 Sec
53
Figures 4.6, 4.7, and 4.8 show the analytical and upwind
solutions of the 1D linear advection equation (4.1) at time
t=50, 250, and 500 sec. It is evident that there is numerical
diffusion in the numerical solution.
0 .0 .0 .6
0.2 :
0 7
0 003
x
Figure 4.6 Results of the 1D linear advection equation
obtained with upwind method at T = 50 Sec
.6 0 i i a p i, 7I.
1:4
0 .0 0 .2 0 I
00 000 0O
Figure 4.7 Results of the 1D linear advection equation
obtained with upwind method at T 250 Sec
..... .m ...... .a u t 0 .... ....
0 9
0 a 0,7 0 .0
Figure 4.8 Results of the 1D linear advection equation
obtained with upwind method at T = 500 Sec
Figures 4.9, 4.10, and 4.11 show the results of the Leith's method. Although there is no numerical diffusion in the solution,. excessive oscillation is evident.
00
0.2 0 7 0
Figure 4.9 Results of the 1D linear advection equation
obtained with the Leith's method at T = 50 Sec
0 0
0 o0
00 0 .00
02 0
00 10 (0 00
Figure 4.10 Results of the ID linear advection equation
obtained with the Leith's method at T = 250 Sec
0 .7 0 .6 220 5 .0
0 .5
.4
0
C2 250 500 760 1000
Figure 4.11 Results of the 1D linear advection equation
obtained with the Leith's method at T = 500 Sec
Figures 4.12, 4.13, and 4.14 show the results of the TVD method. The numerical and analytical solutions are in very good agreement.
Figure 4.12 Results of the 1D linear advection equation
obtained with the TVD method at T = 50 Sec
A " 4.... I 14 Io Io
o 5
Figure 4.13 Results of the 1D linear advection equation
obtained with the TVD method at T =250 Sec
A.
0 .s 00
.. . '77', B . .
Figure 4.14 Results of the lD linear advection equation
obtained with the TVD method at T = 250 Sec
A
7j
o .
o .5
. 3 0 2
Figure 4.14 Results of: the 1D linear advection equation
obtained with the TVD method at T = 500 Sec
56
Figures 4.15, 4.16, and 4.17 show the numerical results
obtained with the ELM with linear interpolation and a CFL
number of 5.5. The numerical and analytical solutions are
clearly in good agreement.
* LS .II II.
t,a
I:
10. .
0.s
2. 0 7 01 0
.tx
Figure 4.15 Results of the 1D linear
obtained with the ELM with linear at T = 55 sec
advection equation interpolation
.. ......... .... ..
o .
3
Figure 4.16 Results of the ID linear advection equation
obtained with the ELM with linear interpolation at T = 253 sec
E .. P d:9
0 0.0
0 .4
S070.00
Figure 4.17 Results of the 1D linear advection equation
obtained with the ELM with linear interpolation at T = 506 Sec
57
Figure 4.18, 4.19, and 4.20 show the numerical solution of the Equation (4.5) obtained with ELM with quadratic interpolation and a CFL number of 5.5. The numerical solutions show wiggles in the vicinity of sharp gradients.
.....o, ...
~o .s o
.0.1
250 10 750 000
Figure 4.18 Results of the ID
obtained with the ELM with at T = 55
linear advection equation linear interpolation Sec
1 2LU I I I I d. 0
v o< .8% ,. : .
Io
0 4
o .5
o. 7 A8 1
Figure 4.19 Results of the ID linear advection equation
obtained with the ELM with linear interpolation at T 253 Sec
0 .4
o I
42
o ,
Figure 4.20 Results of the ID linear advection equation
obtained with the ELM with linear interpolation at T = 506 Sec
58
Figures 4.21, 4.22, and 4.23 show the numerical solution of the conservative QMSL method. The numerical solution has no spurious numerical oscillation and agrees well with the analytical solution.
o. o ..
o .7
1 0 .0
Figure 4.21 Results of the ID linear advection equation
obtained with the Conservative QMSL at T = 55 Sec
0 ...... .n ll.I ...0
U0
0
0 0
0 4
o20 000 000 100
Figure 4.22 Results of the ID linear advection equation
obtained with the Conservative QMSL at T = 253 Sec
U0o
0
0) 0.0
~~ .5 o'i
0 I 0 000I i 00 ? 0 1000
Figure 4.23 Results of the ID linear advection equation
obtained with the Conservative QMSL at T = 506 Sec
4.5.2 Conservation Test
If one wants to simulate some conservative quantities numerically, it is necessary to use numerical schemes that possess the conservative property. Conservative property of a numerical scheme is that conservative species satisfy the integrated form of the conservative equation i.e.
o,+A x=L +1 dxdt = o
(4.35)
Thus :
J x=L(+ c j) +][(uc)Rc(c .d
fx=O _C c)dx+ (uc)Lldt = 0
(4.36)
where subscript R, and L are right and left face of computational cells as shown in Figure 4.24.
Figure 4.24 Computational cell used in conservation test
S0Lcq
60
If the convective velocity is constant, the convective equation is the same in either conservative or nonconservative form. Hence, to check the conservation of
numerical scheme, the convective velocity should be nonconstant.
The test problem for the lD linear advection equation is the same as that discussed in section 4.4.1. However, the convective velocity is nonconstant
U= 1+ (4.37)
L
where L is the length of the computational domain.
Figure 4.25 shows the conservation error of the
conventional ELM with linear interpolation, while Figure 4.26 shows the result of the conservative QMSL method. It is seen
that the conservation error of QMSL fluctuates around zero with the amplitude on the order of 1i13 while that for the conventional ELM does not converge to zero. In both figures, conservation error is defined as:
Conservation Error = Total Initial Mass Computed Mass
Figure 4.25 Conservation error of the conventional ELM
Figure 4.26 Conservation error of the QMSL
4.6 Momentum advection equation 4.6.1 NonLinear Advection Case
The EulerianLagrangian Method described in previous sections have also been applied to strictly nonlinear problems. The invicid Burgers equation is used for model tests. The nonconservative form of the governing equation is
du du
+ uk 0
whereas the conservative form of governing equation is
1
S u2
+ 0
(4.38)
(4.39)
o
5E5512
1E11. T
Elo 2 d0 300 400
E psd Tim.
The initial and boundary conditions are assume to be :
u(x,0) = (1 x) 0 x< 1
u(x,O) = 0 1< x 3 (4.40)
u(0,t) = 1 t> 0
The analytical solution of the problem, for t> 0, is given by :
u(x,t) = 1 for x t
u(x,t)= [(1x)/(1 t)] for t< x< 1 (4.41)
u(x,t) = 0 for 15 x
At time t=l a shock develops in x=l and propagates to the right at speed 0.5. Thus the analytical solution for t_ 1 is given by :
u(x, t) = 1 for x < 1+0.5(t 1)
u(x,t)= 0 for x> 1+ 0.5(t1) (4.42)
The spatial domain 0:5x 3 has been divided into 150 intervals of Ax = 0.02. Nonconservative upwind, conservative upwind, Leith's method, Godunov method, nonconservative and conservative ELM have been implemented and compared to the analytical solution.
The first scheme is the nonconservative upwind scheme. The stability of the scheme is restricted by the CFL number. The time step is set to 0.01 sec, with the max CFL number of
0.5.
63
Figures 4.27, 4.28 and 4.29 show comparison of the numerical and analytical solutions when elapsed time T=0.5, 1, and 1.5 second. The results show that the scheme is accurate but there is some numerical diffusion in the numerical
solution. This scheme does not capture the shock, i.e., the shock does not propagate after the shock develops.
Figure 4.27 Results of the nonconservative
upwind method at T 0.5 sec
0 3
0 2
Figure 4.28 Results of the nonconservative
upwind method at T 1.0 sec
'0A
o 7
0 6
0
a 3
V
Figure 4.29 Results of the nonconservative
upwind method at T = 1.5 sec
64
Figures 4.30, 4.31 and 4.32 show numerical solution of
the conservative upwind numerical scheme. The solution is accurate and capture the shock, but the propagation speed of
the shock is incorrect. The numerical solution is correct before the shock is developed.
0 3
0 2
or
Figure 4.30 Results of the conservative
upwind method at T 0.5 sec
0 .6 S 0.;
0 4
0 .2
Figure 4.31 Results of the conservative
upwind method at T 1.0 sec
0 3
0 .2 0 .1
0
Figure 4.32 Results of the conservative
upwind method at T = 1.5 sec
Figures 4.33, 4.34 and 4.35 show the results of the Leith's method for the nonlinear advection equation. This method is accurate until the shock develops but does not capture the shock. Moreover the solution has unphysical oscillation.
0. 6 0.8
.6
W 0
Figure 4.33 Results of the Leith's method for nonlinear advection equation at T = 0.5 sec
Figure 4.34 Results of the Leith's method for nonlinear advection equation at T = 1.0 sec
o .2
i o ..
Figure 4.35 Results of the Leith's method for nonlinear advection equation at T = 1.5 sec
Figures 4.36, 4.37 and 4.38 show the results of the
Godunov method, which is a kind of TVD. It captures the shock
correctly, and the solution is accurate.
0 ....< . .. .. .
D.3
0.2
.1
0o
x
Figure 4.36 Results of the Godunov method for nonlinear advection equation at T = 0.5 sec
0 .4 04 3
.3
0.2 o I
Figure 4.37 Results of the Godunov method for nonlinear advection equation at T = 1.0 sec
0 .4 0.3
0 .2 0 .1 =On 2
oX
Figure 4.38 Results of the Godunov method for nonlinear advection equation at T = 1.5 sec
67
Figures 4.39, 4.40 and 4.41 show the results of conventional nonconservative ELM with linear interpolation. The maximum CFL number is 2.5. It captures the shock correctly, and solution is accurate.
0 0 5.. ,. .
UO 4
Figure 4.39 Results of the nonlinear advection equation
of ELM with linear interpolation at T =0.5 sec
.0 E
0 0
Figure 4.41 Results of the nonlinear advection equation
of ELM with linear interpolation at T = 3.25 sec
68
Figures 4.42, 4.43 and 4.44 show the results of the conventional ELM with quadratic interpolation. It captures the shock correctly, but there are some numerical oscillation near the shock.
o .7
Figure 4.43 Results of the nonlinear advection equation
of ELM with quadratic interpolation at T 0 .0 sec
0 .2
0 .0
Figure 4.44 Results of the nonlinear advection equation
of ELM with quadratic interpolation at T = 3.5 sec
o .7 v T ,c 1L s.
a 5 0.4 o .o
Figure 4.44 Results of the nonlinear advection equation
of ELM with quadratic interpolation at T = 3.25 sec
69
Figures 4.45, 4.46 and 4.47 show the results of the conservative QMSL. It captures the shock correctly, and its results are accurate.
07 0.0
o
o 4 1 I Ia
Figure 4.45 Results of the nonlinear advection equation of conservative QMSL method at T : 0.5 sec
I 4 00
.o
.sx
Figure 4.46 Results of
of conservative
the nonlinear advection equation QMSL method at T = 1.0 sec
Figure 4.47 Results of the nonlinear advection equation
of conservative QMSL method at T = 3.25 sec
.4
.o
o .2 o I
a
4.6.2 Momentum Conservation Check
To test the momentum flux conservation, the following integral form of equation and conservation form of finite difference equation are used.
I+AI 01l 1 X =L 011t
adt + l= dx=0 (4.41)
Flux difference in time should be Zero.
( Icells+ I n+ ( Icells+1 )n+1 Icells+1 n Icells+1 n i=0 U i+1) i=ii= i+1 i=0 i = 0
(4.42)
The conservation error is defined as :
Conservation Error = Total Flux Computed Flux = 0
Figure 4.48 shows the conservation error of the conservative QMSL algorithm. The error converges to zero with time.
Conservation ofCons.rvaIv, ELM
1.0005
1 .0 0 3 i .00 02 1 .0 0 0 1
Figure 4.48 Conservation error of the conservative QMSL
4.7 TwoDimensional Conservation Test
The conservative and nonconservative EulerianLagrangian Method (ELM) were tested in twodimensional space. The physical domain is 200 km in X direction and 200 km in Y direction. The computational domain is divided into 1 Km in X and Y direction, i.e., 200x200 cells. The initial and boundary conditions are :
Initial conditions
c (5 km < x < 25 km, 0) = 1.0
c (5 km < y < 25 km, 0) = 1.0 (4.43)
c = 0, elsewhere
Boundary conditions
c(OOt) = 0.0
c(200 km, 200 km, t) = 0.0 (4.44)
The advection velocity is uniform, i.e., u = v = 30.5 cm/s. Figure 4.49 is numerical solution of the nonconservative ELM with time increment At = 1 hours, i.e, CFL number is 1.098 for 5 days simulation and Figure 4.50 shows the global conservation error of the nonconservative ELM. The global conservation error is defined by
Computed concentration Total concentration
Error(%) = x 10
Totol concentration
(4.45)
z
cF~nnevtv101 ELM Y
[ X
0.0
0.6 0.4
X(Krn)
Figure 4.49 Numerical solution of 2D nonconservative ELM
2 E 1 2 1 E 1 2
0
1 E 1 2
2 E 1 2
Figure 4.50 Conservation error of 2D nonconservative ELM
2(0 40 60 80 100 12
Tim e(H ours)
73
Figure 4.51 shows the contour map of numerical solution of nonconservative ELM. There are some numerical damping in the numerical solution.
200 175 150
125
x100
75 50 25
0o
2D NonConservative
iF!
2 days
1 lday ]Inital Condition 1:11 I.II
5 days after
I I I I I
50 100 150 200
X(km)
Figure 4.51 Contour map of numerical
nonconservative ELM
solution of
Figure 4.52 shows the numerical solution of the conservative ELM and Figure 4.53 shows the conservation error of the conservative ELM. Figure 4.52 indicates the numerical solution of conservative ELM has less damping and more accurate than nonconservative scheme.
SELM
I
0.8 o 0.6 0.4 0.2
0
z
200
X(Krn)
Figure 4.52 Numerical solution of 2D conservative ELM
2 E 12 S1E12
o
0 0
0
1 E 1 2
2 E 12 ......
Figure 4.53 Conservation error of 2D conservative ELM
75
Figure 4.54 shows the contour map of numerical solution of nonconservative ELM. The numerical damping of the
conservative ELM is much less than that of nonconservative ELM.
)fiHlI
175 150
125
EIO
.Sd 100
75 50 25
U 50 100 150 200
X(km)
Figure 4.54 Contour map of numerical solution of conservative ELM
2D Conservative ELM
5 days after
4 days
S3 days
:0 2 days
0 1 day
Initial Condition I I 1 I I ,
4.8 ThreeDimensional Conservation Test The conservative and nonconservative EulerianLagrangian Method (ELM) were tested in threedimensional rectangular basin. The physical domain is 200 km in X direction,200 km in Y direction and depth is 150 meters. The computational domain is divided into 1 Km in X and Y direction and 15 meters in Z direction, i.e., 100xlOOxlO cells. The initial and boundary conditions are :
Initial conditions
c (10 km < x < 50 km, 0) = 1.0
c (10 km < y < 50 km, 0) = 1.0 (4.46)
c(20 m < z < 70 m, 0) = 1.0
C = 0, elsewhere
Boundary conditions
c(0,0,0,t) = 0.0
C (200km,200km,150m, t) = 0.0 (4.44)
The advection velocity is uniform, i.e., u = v = 30.5 cm/s and w = 0.012 cm/s. Figure 4.55 shows the numerical solution of the nonconservative ELM with time increment At = 2 hours, i.e, CFL numbers in X and Y direction are 1.098 and CFL number in Z direction is 0.0576 for 5 days simulation and Figure 4.56 shows the global conservation error of the nonconservative ELM. The numerical solutions of 3D conservative and nonconservative ELM is very close, but the conservation error of
77
the two different schemes is significant. Figure 4.56 shows
the mass of concentration in nonconservative numerical scheme lose its original concentration with time.
3D NonConservative ELM
z
6yx
Figure 4.55 Numerical solution of 3D nonconservative ELM
0
*.0 .5
1i 5
.2
2.5 7
3 ,. . . . . .
Figure 4.56 Conservation error of 3D nonconservative ELM
Figure 4.57 and 4.58 show the numerical solution of the nonconservative ELM in XZ and YZ plane, respectively.
3D NonConservative ELM
1
4 days a fte r
2 day after
In a Itia. C.ndltipn 7
IIi ii. ___l_____.....'JHH '....."....i~/ "... , ,
0 S0 100 150 200
X (km)
X
50 25
00
N
Figure 4.57 Numerical
solution of 3D nonconservative ELM in XZ plane
3D NonConservative ELM
z
Y
Y(km)
Figure 4.58 Numerical solution of 3D nonconservative ELM in YZ plane
79
Figure 4.59 shows the numerical solution of the conservative ELM. Figure 4.60 Shows the conservation error of the conservative ELM. Conservation error of the conservative ELM is almost zero.
3D Conservative ELM
z
6yx
Figure 4.59 Numerical solution of 3D conservative ELM
2 E .1 2 1 .8 E 1 2 1 .6 E 1 2 1 .4 E .1 20
1 2 3 4
T im e (D a ys )
Figure 4.60 Conservation error of 3D conservative ELM
              
Figure 4.61 and 4.62 show the numerical solution of the nonconservative ELM in XZ and YZ plane, respectively. Figures show conservative ELM has less damping than nonconservative ELM.
3D Conservative ELM
150 125
100
50 25
X (km)
Figure 4.61 Numerical
solution of 3D nonconservative ELM in XZ plane
3D Conservative ELM
Y 50
25 D0 N f
Y (kin)
Figure 4.62 Numerical solution of 3D nonconservative ELM in YZ plane
4.9 Summary
All the above presented numerical schemes for the 1, 2 and 3dimensional linear and nonlinear advection equations are relatively accurate. The conventional ELM is accurate but it does not possess the conservative property. If a higherorder interpolation function is used, there is spurious unphysical oscillation in the numerical solutions. The total variation diminishing (TVD) method produced the best results in terms of both accuracy and conservation. However, the TVD method is restricted by the CFL condition. The conservative quasimonotone semiLagrangian (QMSL) algorithm is relatively accurate, has the conservative property, and has no numerical oscillations with the help of FCT algorithm. In addition, the conservative QMSL is unconditionally stable.
Because of its numerous advantages, the conservative QMSL method will be used for the present threedimensional model discussed in chapter 3.
