- Email: [email protected]
A comparison of hydrostatic and nonhydrostatic pressure components in seiche oscillations
MATHEMATICAL AND COMPUTER MODELLING
Available online at www.sciencedirect.com
8CIENCl~DIRECTe
ELSEVIER
Mathematical and Computer Modelling 41 (2005) 887-902 www.elsevier.com/locate/mcm
A Comparison of Hydrostatic and Nonhydrostatic Pressure Components in Seiche Oscillations XINJIAN CHEN Resource Conservation and Development Department Southwest Florida Water Management District 7601 Highway 301 N o r t h , T a m p a , F L 33637, U.S.A. xinj ian. chen©swfwmd, state, fl. us
(Received April 2004; revised and and accepted August 2004)
A b s t r a c t - - T h e hydrostatic pressure assumption has been widely used in studying water movements in shallow waters where the water d e p t h (D) is much smaller t h a n the horizontal length scale (L). It is assumed t h a t the pressure field in a shallow water flow is hydrostatic because the nonhydrostatic pressure component (Pn) is much smaller t h a n the hydrostatic pressure component (Ph) and thus be safely neglected in model simulations. W i t h this assumption, a three-dimensional, hydrodynamic model can be significantly simplified because the m o m e n t u m equation in the vertical direction is reduced to an expression of hydrostatic pressure. While such a simplification is valid in many shallow water applications, there exist no rigorous guidelines as to what extent can the hydrostatic assumption be safely used without any problems. In an effort to get a better understanding of the effect of nonhydrostatic pressure on free-surface flows, this paper presents a comparison of hydrostatic and nonhydrostatic pressure components and their horizontal gradients in two seiches. T h e first seiche represents a deepwater wave oscillation, while the second one is a shallow water wave oscillation. A nonhydrostatic model was used to simulate the two oscillations, and simulated hydrostatic and nonhydrostatic pressure components were compared. A scale analysis was also carried out for the two pressure components in seiching. It is found t h a t although the nonhydrostatic component was about 2 to 3 orders of magnitude smaller t h a n the hydrostatic component for the two oscillations, the horizontal gradient of the nonhydrostatic component is not necessarily smaller t h a n t h a t of hydrostatic pressure by 2 to 3 orders of magnitude. On the other hand, while the horizontal gradient of nonhydrostatic pressure may b e in the same order of magnitude as t h a t of hydrostatic pressure, it does not necessarily require nonhydrostatic pressure itself to be in the same order of magnitude as hydrostatic pressure. Model simulations also show t h a t the nonhydrostatic portion of the horizontal pressure gradient affects the velocity field during the course of the oscillation and its effects are not negligible even when t h e aspect ratio is less t h a n 0.05. ~) 2005 Elsevier Ltd. All rights reserved.
Keywords--Hydrostatic pressure, Nonhydrostatic pressure, Hydrostatic pressure assumption, Seiche oscillation, Nonhydrostatic model.
1. I N T R O D U C T I O N The hydrostatic pressure assumption has been widely used in studying shallow rivers, lakes, estuaries, and continental shelves. This assumption is valid in most cases and has been successfully 0895-7177/05/$ - see front m a t t e r (~) 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2004.08.005
Typeset by ~4A/tS-~I~_tX
888
X. CHEN
used in many three-dimensional, hydrodynamic models (e.g., [1-5]). However, there are also many cases where this assumption may be questionable [6-8]. One example is the flow induced by a very large horizontal density gradient. Other examples include short waves, flows through structures, the near field of brine disposal from a desalination plant, and flows near an uptake point of a water withdrawal system. In these cases, nonhydrostatic pressure effects may be comparable to those of hydrostatic pressure, and thus cannot be neglected in model simulations. The basic idea of the hydrostatic pressure assumption comes from the fact that the acceleration and eddy viscosity terms in the momentum equation for the vertical velocity component are much smaller than the gravitational acceleration term and can thus be neglected [9,10], because vertical length scales in shallow waters are general much smaller than horizontal length scales. The momentum equations for water movement in three directions can be expressed as follows
p-~x + -~x Ah
d--[ =
dv dt -
+ ~Y k, h ~y J -Jc--~Z \
fu_lOp 0 (AhOV h 0 AhOy ~z\ p -Oyy+ -Oxx\ OxJ + -~y ÷
dw
lop
dt-g-
+Tx
0 fA Ow~
hTx)+
0 fAhOW~
(1)
OZ / '
Oz]
,
(2)
0 fAvOW~
\
(3)
az]'
where x, y, and z are Cartesian coordinates (x is from west to east, y is from south to north, and z is vertical pointing upward); u, v, and w are velocities in the x-, y-, and, z-directions, respectively; t, f, g, and p denote time, the Coriolis parameter, the gravitational acceleration, and pressure, respectively; Ah and Av represent horizontal and vertical eddy viscosities, respectively; and p is density which is a function of temperature and salinity. Equation (3) can be rewritten in the following form
lop
dw
O fAhOW~
p Oz -- g - --~ + -~x \
O fAhOW~
Ox ] + ~y \
0 fAvOW~
Oy ] + -~z \
Oz ] = - g + a3,
(4)
where a3 includes the vertical acceleration of a water particle and the vertical viscosity term. Integrating the above equation from an elevation z to the free surface (~7), we have P - Pa = g
//
p d~-
S/
pa3 de,
(5)
where Pa represents pressure at the water surface (~) and is assumed to be zero throughout this paper. The first integration on the right-hand side of equation (5) is hydrostatic pressure (Ph) and the second integration is nonhydrostatie pressure (P,0. In the hydrostatic pressure assumption, the magnitude of the second integral on the righthand side of the above equation is assumed to be much smaller than the first one because the aspect ratio is very small for shallow waters. Therefore, equation (5) is reduced to p= g
//
p de.
(6)
Obviously, if O(a3) << g, equation (6) is always valid. By replacing total pressure in the horizontal momentum equations with hydrostatic pressure, the number of governing equations for three-dimensional shallow water flows is reduced from four to three: two horizontal momentum equations and a continuity equation
d-~ -= dt --
- pox P-~Y
pd¢ + ~x \ pd¢ + ~ \
Ox) +-~y \ Ox/#+-~y \
ay) + Oz \
Oz]'
Oy] +-~z A . ~
(7)
,
(8)
A Comparison
Ou
Ov
889
Ow
0---~-b ~yy -F ~ z -- 0.
(9)
Many three-dimensional hydrodynamic models (e.g., [1-5]) solve the above three equations. However, this kind of treatment of the horizontal pressure gradient term is mathematical not rigorous because it implies that if the nonhydrostatic portion of pressure is much smaller than the hydrostatic portion, then the horizontal gradients of nonhydrostatic pressure would also be much smaller than those of hydrostatic pressure. In other words, if IP~ [ << Ph, then the following expression would be automatically true:
Opn 0x
OOph
Opn
Oy
<<
Oph
(10)
<<
The above logic is only true when both p~ and Ph have the same length scale in the horizontal direction. However, because hydrostatic pressure may have multiple length scales for shallow water flows, the above logic is not necessarily always correct. There could be cases where ]p~ [<
H
T
cos
(~rx)
(11)
'
where H is the wave height of the oscillation and l~ is the length of the rectangular basin. At time > 0, water in the basin undergoes oscillation under the barotropic pressure gradient caused by the initial free surface setup. In the following sections, a three-dimensional, nonhydrostatic model is briefly described, before model applications to the two oscillations, along with a detailed comparison of hydrostatic and nonhydrostatic pressure components and their horizontal gradients, are presented. Model results with and without the hydrostatic pressure assumption during the oscillations are also discussed. A scale analysis using analytical solutions for small amplitude waves is then described to show Z
D
W
X F i g u r e 1. A seiche oscillation in a r e c t a n g u l a r b a s i n is c a u s e d by a n initial free-surface setup.
890
X. CHEN
the relative scales of nonhydrostatic pressure in comparison with those of hydrostatic pressure. Conclusions drawn from the study are included in the last section of the paper.
OF
A
2. A B R I E F DESCRIPTION NONHYDRODYNAMIC MODEL
A three-dimensional nonhydrodynamic model named LESS3D (Lake and Estuary Simulation System in Three Dimensions) [7,8] was used in this study. The model solves the governing equations for free-surface flows, including equations (1)-(3) and equation (9), by performing two predictor-corrector steps. In the first predictor-corrector step, the model uses hydrostatic pressure at the previous time step as an initial estimate of the total pressure field at the new time step. Based on the estimated pressure field, an intermediate velocity field is calculated, which is then corrected by adding the nonhydrostatic component of the pressure to the estimated pressure field. By forcing the velocity field divergence-free, a Poisson equation for nonhydrostatic pressure is obtained. In a structured grid system, the Poisson equation forms a seven-diagonal matrix system that is positive definite if a uniformly distributed grid system is used. In this case, the conjugate gradient method with incomplete Cholesky preconditioning [11] is used to find the inverse of the matrix. If the grid sizes are not uniformly distributed in the computation domain, the sevendiagonal matrix is generally asymmetric. In this case, a Bi-CGSTAB method [12] is used to solve the matrix. After the nonhydrostatic pressure component is solved, the intermediate velocity field is corrected to obtain the second intermediate velocity field that is the second prediction of the final velocity field. In the second predictor-corrector step, an intermediate free-surface location at the new time step is calculated based on the second intermediate velocity field. The intermediate free surface is then adjusted to find the final free-surface location based on a solution of a free-surface correction equation. The final velocity field is finally solved after the final freesurface location is found. A detailed description of the numerical procedure for implementing the two predictor-corrector steps in the nonhydrostatic model LESS3D is presented in reference [8]. After finding the free-surface location and the velocity field at the new time step, the model solves the following transport equation of concentration oc
ouc
owc
o
(B
O---~+ --.~x + --~y + O---~.-= O-.-.~\
o
o
(Boac
Ox j + -~y ~, h ~y / + -~z \
Oz ) + S + R,
(12)
where c is concentration (can be salinity, temperature, suspended sediment concentration, nutrient species, etc.), S denotes source/sink terms, R represents reaction terms, and Bh and B, are horizontal and vertical eddy d/flus/v/ties, respectively. If the simulated concentration involves settling, w in the above equation includes the settling velocity. In this study, however, because the seiche oscillations do not involve concentration, equation (12) is not solved. The model calculates eddy viscosity and diffusivity in the horizontal directions (Ah and Bh) using Smagorinsky's approach [13]. For the vertical eddy viscosity and diffusivity (Av and B~), a turbulent kinetic energy (TKE) model [14,15] is used. In this study, however, a small constant eddy viscosity of 0.1 cm 2 sec-1 was used in the simulation. Boundary conditions in the horizontal directions are specified with either free surface elevations or velocities for open boundaries. At solid boundary, both the velocity component and the pressure gradient in the normal direction are set to zero. Boundary conditions specified in the vertical direction are shear stresses. At the free surface, wind shear stresses are used. At the bottom, the log-layer distribution of velocity is used to calculate the bottom shear stresses, or =
+
(13)
where a is the vonKarman constant (0.41), Ub and Vb are horizontal velocity calculated at a level zb near the bottom, zo --- ks~30 and k~ is the bottom roughness.
A Comparison
891
The LESS3D model can be run either with or without the nonhydrostatic effect. When the model is run without the nonhydrostatic effect, the model skips the substep of finding the inverse of the seven-diagonal matrix for nonhydrostatic pressure and the substep of correcting the first intermediate velocity field. The model uses the so-cMled z-level in the vertical direction and no transformation is made to the vertical coordinate. The model allows the free surface to travel between layers, so that there is no need to make the top layer thick enough to cover the whole range of free surface variation. To fit the bottom topography, hybrid grid cells are used in the model [16,17]. 3. N U M E R I C A L TWO SEICHE
SIMULATION OF OSCILLATIONS
This section presents applications of the LESS3D model to the two seiche oscillations described in the Introduction section. As mentioned before, the basin length and the water depth are the same in the first seiche oscillation and both equal to 10 m. In the second oscillation, however, the basin length is 44m, which is 11 times of the water depth of 4m. Both oscillations have Hydrodynamic
Pressure
20
A
E 15 o~ t- 10 o
5
U,I
-5
0 o -10 M 1= -15 i
=
-20 0 . . . . . .
......
......
a' . . . . . . . . . . . . .4. . . . . . .
s
.
.
.
.
.
.
.
.
.
.
.
.
6
; ......
; ......
10
Time (second) 20 Ix = O.lm, z f S.gml
15
A
o
5
0 0
0
-5
~-10 -lS .200
A
1
2
3
4
~ ......
5 Time (second)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7, . . . . . . . . . . 8. . . . . . . . . . 9
10
5O 4(1
10 0
-10 -20 -30 -4O -50
. . . . . .
0
I
1
. . . . . .
I
2
,
,
r
,
,
,
I
3
,
,
,
,
,
,
I
4
,
,
,
,
,
,
I
,
,
,
,
i
5 Tlme (second)
,
I
6
,
,
,
i
,
,
I
7
,
,
r
i
,
,
I
8
. . . . . .
I
. . . . . .
9
F i g u r e 2. C o m p a r i s o n of s i m u l a t e d ( s o l i d l i n e s ) s u r f a c e e l e v a t i o n s a n d v e l o c i t i e s w i t h a n a l y t i c a l s o l u t i o n s ( d a s h e d l i n e s ) for t h e first s e i c h e o s c i l l a t i o n w i t h a : 0.5. S u r f a c e e l e v a t i o n s a r e c o m p a r e d a t x = 0.1 m a n d 5.5 m , w h i l e v e l o c i t i e s a r e c o m p a r e d a t x = 0 . 1 m , z = 8 . 9 m a n d x = 5 . 5 m , z = 4 . 9 m . T h e n o n h y d r o s t a t i c effect w a s included in the model simulation.
10
892
X. CHEN Hydrostatic
20 A
Pressure
I/I/Ill'"
15 5 0
-5 -10 -15 -20
1
2
3
4
5 Tlme (second)
1
2
3
4
5 Time (second)
6
7
8
9
10
20 15 10 5 0 0
>:
0
-5 -10 -15 -20 0
10
.--
20
10
," "-'.
0
~-I0 -20
.20
.................. 0
I
. . . .
2
3
4
Ir
5 Tlme (second)
. . . . . . . . . . . . . . . . . . . . . .
6
7
8
9
10
Figure 3. Comparison of simulated (solid lines) surface elevations and velocities with analytical solutions (dashed lines) for the first seiche oscillation with a ----0.5. Surface elevations are compared at x = 0.1 m and 5.5 m, while velocities are compared at x ----0.1m, z ----8.9m and x = 5.5m, z ----4.9m. The nonhydrostatic effect was not included in the model simulation. the same wave height of 30 cm. To c o m p a r e the m a g n i t u d e of n o n h y d r o s t a t i c pressure and its horizontal gradient with those of hydrostatic pressure, b o t h pressure c o m p o n e n t s were saved as model o u t p u t s during the model run. Horizontal gradients of h y d r o s t a t i c and n o n h y d r o s t a t i c pressure c o m p o n e n t s were calculated using a post-process p r o g r a m after t h e simulation was done. W i t h the initial free-surface setup expressed in equation (11), a s t a n d i n g wave is formed in the rectangular basin with a wavelength t h a t is twice the length of t h e basin (L = 2/~). Obviously, the horizontal length scale of the oscillation is the wavelength (L) and the vertical length scale is t h e water d e p t h D. Thus, the aspect ratio (~) equals D / L . I n the first oscillation, because the water d e p t h is the same as the basin length, the D / L ratio is 0.5. It is expected t h a t the n o n h y d r o s t a t i c effect is significant for seiching in this kind of basin. In the second oscillation, the aspect ratio ~ is relatively small and has a value of 0.0455. In b o t h seiche oscillations, the wave slopes are m u c h smaller t h a n 1. As such, analytical solutions according to the linear wave t h e o r y can be o b t a i n e d as follows H ~- ~- cos kx c o s a t ,
(14)
(u,w) - H ~ s i n ~ t ( s i n k x c o s h k z , - coskxsinhkz) (15) 2 sinh kD where k --- 2~r/L -~ ~/l~ and a is related to k by the dispersion relationship cr2 --- g k t a n h kD.
A Comparison
893
Model Simulation of the First Seiche Oscillation
When the first oscillation was simulated using the LESS3D model, a uniform grid system was used with a grid size of 0.2 m in both the horizontal and vertical directions. The model was run for 10 s with a time step of 0.05 s. From the dispersion relationship, the wave period can be calculated and equals to 3.59 s. Figure 2 show the model results of surface elevations and velocities at two locations in the basin, together with the analytical solutions. One location is near the left end of the basin with the coordinates of x = 0.1 m and z = 8.9 m, while the other location is near the middle of the basin at x = 5.5m and z = 4.9m. The model was run with the consideration of nonhydrostatic pressure. The model results are plotted with solid lines and the analytical solutions axe in dashed lines. As can be seen from the figure, the model was able to reproduce the analytical solution fairly well. If the nonhydrostatic pressure effects were not included in the model simulation, the simulated surface elevation and velocity field would be totally different from analytical solutions. Figure 3 shows model results (solid lines) without the consideration of nonhydrostatic effect. Analytical solutions are also shown in the figure (dashed lines). Obviously, because the nonhydrostatic effect is significant in this oscillation case, the hydrostatic pressure assumption is not valid and the use of this assumption produces incorrect model results. Two things can be noticed from Figure 3. One is that the simulated oscillation has a much shorter wave period than the analytical solution does. The other thing is that higher mode oscillations are introduced to the model results shortly after the initial free surface setup. It is clear that forcing the horizontal pressure gradient to
100000
-9000 -7000 -5000
-3000 - 1 ~
111110 3000
5000
7000
9000
(b)
(~) .~.~:%~~.~.~i:~~::~!,~' ~
Ip,~pd
~!:il:~;I!~I~::,
dp./dx o.ooo o.oos O.OLO o.o~5 o.o-2o 0.025
II
I:::~:~1:=:::!~lii~:'~=]i~i~~ m~ ~ -~ -~ o
-ro-~e.~e-~4-~2-1o-~
(d)
(c)
ldpJdxl/IcWJdxl
dp~/dx 0.0
4.0
8.0
(e)
12.0
16.0
20.0
I I
t::
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5,0
(f)
F i g u r e 4. M o d e l r e s u l t s at t i m e = 6 s for t h e first seiche oscillation w i t h (a), (b), a n d (c) a r e s i m u l a t e d n o n h y d r o s t a t i c p r e s s u r e (p~), h y d r o s t a t i c (Ph), a n d t h e a b s o l u t e value of t h e ratio o f p n over Ph, respectively. (d) a n d horizontal g r a d i e n t s of p n a n d Ph, respectively. (f) is t h e a b s o l u t e v a l u e of of t h e two g r a d i e n t s .
(~ ~ 0.5. pressure (e) s h o w t h e ratio
894
X. CHEN 100000
100000--
~00c : :::::~::i~::i:i: ~:iti::: ::i::: ::: :i::i:i:::::j:::!i ~ :
p~ I]
I I
-9000 -7000 -5000 -3000 -1000
1000
3000
5000
7000
9000
1.0E+05
3.0E+OS
(a)
IpJpd
ill:::; ::
0.000
0.010
0.020
5.0E+05
7,0E+05
g.OE+05
(b)
0.030
dpr/dX
0.040
-17 -13
-9
(c)
-5
-1
3
7
11
1.5
2,0
15
19
23
(d)
I~ I ~
llit t
i~!~:il H~l:~i[
I
itll ~i:l!!l~i[ I
dp~dx
i i
I -27.0
Idp./dxl / Idpjdxl I I -21.0
-15.0
-9.0
-3.0
3.0
0.0
0.5
1,0
2.5
3.0
3.5
(e) (f) Figure 5. Model results at time = 7.7s for the first seiche oscillation with a -- 0.5. (a), (b), and (c) are simulated nonhydrostatic pressure (pn), hydrostatic pressure (Ph), and the absolute value of the ratio ofp~ over Ph, respectively. (d) and (e) show horizontal gradients of pn and Ph, respectively. (f) is the absolute value of the ratio of the two gradients. be the horizontal gradient of hydrostatic pressure is equivalent to forcing the phase speed of the wave to be V ~ , which corresponds to a wave period of 2.02 s for this seiching case. This is the same wave period shown in the simulated oscillation in Figure 3, which is much shorter t h a n the analytical solution of 3.59 s. The model simulation with the consideration of nonhydrostatic pressure allows both hydrostatic pressure and nonhydrostatic pressure as well as their horizontal gradients to be output and compared. Figures 4 and 5 show contour plots of simulated hydrostatic pressure and nonhydrostatic pressure at time = 6s and 7.7s, respectively, in the first oscillation. In Figures 4 and 5, (a) and (b) show simulated nonhydrostatic pressure (Pn) and hydrostatic pressure (Ph), respectively. Their corresponding horizontal gradients are shown in (d) and (e), respectively. Figures 4c and 5c are plots of the absolute value of the ratio of p~ over Ph, while Figures 4f and 5f are plots showing the absolute value of the ratio of the horizontal gradients of the two pressure components. As expected, contours of hydrostatic pressure are parallel to the free surface because there is no baroclinic effect in the oscillation. The magnitude of hydrostatic pressure varies linearly from 0dyn-cm -2 at the free surface to about 1 × 106 dyn-cm -2 near the b o t t o m . On the other hand, the nonhydrostatic pressure distribution in the basin is more complex. It can be either negative or positive, with the magnitude being as large as 9,000dyn-cm -2 during the course of the oscillation. Because the horizontal gradient of hydrostatic pressure is not a function of z, its contours are vertical straight lines. Because the horizontal gradient of the free surface approaches zero near the two ends of the basin, the horizontal gradient of hydrostatic pressure also approaches zero near
A Comparison
895
the two ends. The horizontal gradient of nonhydrostatic pressure varies not only with the time and the x-coordinate, but also with the z-coordinate. The magnitude of the horizontal gradient of nonhydrostatic pressure is almost the same as that of the horizontal gradient of hydrostatic pressure. From Figures 4 and 5, it can be seen that although nonhydrostatic pressure is about two orders of magnitude smaller than hydrostatic pressure, their horizontal gradient are comparable in magnitude for this seiching where the aspect ratio is 0.5. M o d e l S i m u l a t i o n of t h e Second Seiche Oscillation For the second seiche oscillation, a uniform grid spacing of 0.5 m was used in the horizontal directions, while the vertical grid spacing varies between 10 cm and 30 cm. The time step used was 0.1 s. From the dispersion relationship, the phase speed of the wave is 6.18m/s and the period of the first mode oscillation is about 14.24s. The phase speed calculated from the dispersion relationship is about 1.3 percent smaller than that calculated from ~ / ~ , which equals 6.26 m/s and would be the phase speed under the hydrostatic assumption.
Hydrodynamic Pressure
25 ,..2O
E 15 o
g o ,-n -5
~ -10 ~ -15
m 2O 20
411
[x = 32.25m,
30
A
30
z = 2.85m
40
70 80 Time (second)
90
100
11 0
120
130
140
I
g2o lO
_~
o
¢)
2 -10 Q
>,.20
I~-- o ys~, ~ = o . .
-30
-40
. . . . . .
I
0
. . . . . .
I
10
. . . . . .
I
20
i
. . . . . .
30
I
. . . . . .
[
40
. . . . . .
i
50
. . . . . .
i
,
,
,
,
,
,
i
. . . . . .
60 70 80 Time (second)
i
. . . . . .
90
i
. . . . . .
100
i
. . . . . .
11 0
i
. . . . . .
120
i
. . . . . .
130
140
10
~5
Ix = 0.75m, z = 0 . 9 5 m j
o
~o 4}
-5 = 32.z5,,,, ~ = 2.85,,, I mlO
. . . . . .
0
I
. . . . . .
10
I
. . . . . .
20
I
. . . . .
30
,
~ I
40
. . . . .
I
. . . . . .
50
i
. . . . . .
I
. . . . . .
i
,
60 70 80 Time (second)
[
,
,
i
, 1
.....
90
r I
. . . . . .
100
I
. . . . . .
110
I
. . . . . .
120
I
. . . . .
130
F i g u r e 6. C o m p a r i s o n of s i m u l a t e d (solid lines) s urfa c e e l e v a t i o n s a n d v e l o c i t i e s w i t h a n a l y t i c a l s o l u t i o n s ( d a s h e d lines) for t h e second seiche o s c i l l a t i o n w i t h a :- 0.0455. Surface e l e v a t i o n s are c o m p a r e d at x = 0.75 m a n d 32.25 m, w h i l e v e l o c i t i e s axe compared at x : 0.75m, z : 0.95m and x = 32.25m, z = 2.85m. The nonhydrostatic effect w a s i n c l u d e d in t h e m o d e l s i m u l a t i o n .
140
X. CHEN
896
Hydrostatic Pressure 25 ..20 E 15 10 = 5 _~ 0
i -5
.20 10
20
30
60 70 80 Time (second)
90
100
110
40
40
30
[~=32.3s~z--z.ssm[
lO
.~,
o
_~-10 >-20 -3o
I ~ = o.Tsm, ~ = o . . m
-40 0 ...... 10'...... 20' ...... 30' .....
~
I
..... 50'"
....60' ...... 70I...... 80'...... 90'...... 100l...... 11' 0...... 120'...... 130'...... 140 Time (second)
10 x = 0.75m, z -- 0 . 9 5 m 1
/h 0
A
: -10
"
-
"
'
'
M " &' '
-
Ix = 32.2s.,
0 . . . . 1'0. . . . 2'0. . . . ~
....
40 . . . .
= 2.SSm I ~
.... 6
.... iO .... ~
.... ~
.... i~
.... 110 .... li0 .... i30 ....
MO
Time (second) F i g u r e 7. C o m p a r i s o n of s i m u l a t e d (solid lines) s urfa c e e l e v a t i o n s a n d v e l o c i t i e s w i t h a n a l y t i c a l s o l u t i o n s ( d a s h e d lines) for t h e s e c ond seiche o s c i l l a t i o n w i t h ~ = 0.0455. Surface e l e v a t i o n s are c o m p a r e d at x = 0.75 m a n d 32.25 m, w h i l e v e l o c i t i e s a re comp a r e d a t x = 0 . 7 5 m , z = 0 . 9 5 m a n d x = 3 2 . 2 5 m , z ---- 2 . 8 5 m . T h e n o n h y d r o s t a t i c effect was not i n c l u d e d in t h e m o d e l s i m u l a t i o n .
The model was run for 142 s (about ten cycles of the first mode oscillation) with and without the hydrostatic pressure assumption. For the nonhydrostatic simulation, both hydrostatic and nonhydrostatic pressure fields were output. Similar to Figures 2 and 3, simulated time series of surface elevation and velocities at two locations are compared with analytical solutions in Figures 6 and 7 for simulations with and without the consideration of the nonhydrostatic effects, respectively. In the figures, solid lines are model results and dashed lines are analytical solutions. The two locations are at (x, z) -- (0.75 m, 0.95 m) and (x, z) = (32.25 m, 2.85 m). It can be seen from Figures 6 and 7 that both nonhydrostatic and hydrostatic simulations yield similar model results during the first couple of cycles because the nonhydrostatic effect is not significant at the beginning of the seiche oscillation case. However, because the error caused by the exclusion of the nonhydrostatic effect is cumulative for this initial value problem, model results of the hydrostatic simulation start to stray away from the first mode analytical solution as time goes on, while model results of the nonhydrostatic simulation still have a reasonable match with the analytical solution. Near the end of the simulation period, it is evident that model results without the nonhydrostatic effect possess a phase lead with respect to the analytical solution as the cumulative effect of the slightly larger phase speed under the hydrostatic assumption becomes significant. Clearly, model
A Comparison
~,~,~, ,~
,,~
897
' ~ff~ ~ ~7~7~#,7~i~7'>".~i~:~i~ '~e
:<~
i7~i77i:~:ii!i
p. I I -400
I -300
I......... ~ . ~ , ~ / ~ ~ -200
-100
0
100
200
P~
300
II
!:::: : : f i~t7q 7 ! : : ' i % 1 7 Z ~ /
O.OE+O0
1.0E+05
(a)
I[
0
I: I
0,001
:I,:~. l 0.002
3.0E+05
4.0E+05
(b/
~
IpdpJ
2.0E+05
h ~
0.003
~
¢
!
dp,,/dx
0,004
~
[
-0.4
i
:
-0,3
~
-0.2
~"
-0.1
0.0
~ ::~i~:~ ]
0.1
0.2
(d)
(c)
.......:::;, ::: :i.:: !!4~ : ' ~!~i!'~
dp~/dx
0.0
1.0
i
2.0
3.0
4,0
5.0
O.Q5rL 7! "!
O.(X) 0.05 0.10 0,15 0.20 0.25 0.30 0.35 0.40
6.0
(e)
(0
F i g u r e 8. M o d e l r e s u l t s at t i m e = 130.4s for t h e second seiche oscillation w i t h = 0.0455. (a), (b), a n d (c) are s i m u l a t e d n o n h y d r o s t a t i c p r e s s u r e (p~), h y d r o s t a t i c p r e s s u r e (Ph), a n d t h e a b s o l u t e value of t h e ratio ofpn over Ph, respectively. (d) a n d (e) s h o w horizontal g r a d i e n t s of pn a n d Pu, respectively. (f) is t h e a b s o l u t e value of t h e ratio of t h e two gradients.
results with the nonhydrostatic effect have a much better agreement with analytical solutions than those without the nonhydrostatic pressure effect, especially after the first couple of wave cycles. Figures 8 and 9 are simulated hydrostatic and nonhydrostatic distributions at time = 130.4s and time -- 137.6 s, respectively. Similar to Figures 4 and 5, (a) and (b) in Figures 8 and 9 show nonhydrostatic pressure and hydrostatic pressure, respectively, while (d) and (e) are plots of their horizontal gradients. The Pn over Ph ratio and the °o~ over ~ ratio are plotted in (c) and (f) of the two figures, respectively. From Figures 8 and 9, it can be seen that nonhydrostatic pressure is at least two orders of magnitude smaller than hydrostatic pressure for the second seiching case that has an aspect ratio of a - 0.0455. Again, the horizontal gradient of hydrostatic pressure is not a function of z, while the horizontal gradient of nonhydrostatic pressure varies with both x- and z-coordinates. Obviously, the horizontal gradient of nonhydrostatic pressure is no longer two orders of magnitude smaller than that of hydrostatic pressure. It is only about one order of magnitude smaller than the horizontal gradient of hydrostatic pressure.
4. A
SCALE
STUDY
To get a better understanding of the model results presented in the last section, a scale analysm for hydrostatic and nonhydrostatic pressure components and their horizontal gradients was conducted using analytical solutions expressed in equations (14) and (15). From the unsteady
898
X. CHEN
: ::: ~:: :: : ~ h ~H~
~ "
-400
-
-
~
-300
¥
:
~
-200
.~
-100
.
0
.
.
100
.
200
~ . ~
Ph
300
O.OS+O0
1.0E+05
(a)
IP~/Phl
II0
0.001
0.002
2.0E+05
3,0E+05
4.0E,P05
(b)
0.003
-0.1
0.004
-0,0
(c)
0.1
0.2
0.3
0.4
0.5
(d)
f: :i::i !:!:i.li:~'~' Idp~/dxl t IdpJdxl -5.0
-4.0
-3.0
-2.0
-1.0
0.00 0.05 0,10 0.15 0.20 0.25 0.30 0.35 0.40
0.0
(f)
(e)
Figure 9. Model results at time = 137.6s for the second seiche oscillation with cz = 0.0455. (a), (b), and (c) are simulated nonhydrostatic pressure (pn), hydrostatic pressure (Ph), and the absolute value of the ratio of p~ over Ph, respectively. (d) and (e) show horizontal gradients of p~ and Ph, respectively. (f) is the absolute value of the ratio of the two gradients. Bernoulli equation, the total pressure can be o b t a i n e d and takes the following form [19]:
p = - p g z + pg~
cosh k(D + z) cosh k D
(16)
Splitting the pressure expressed in the above equation into h y d r o s t a t i c and n o n h y d r o s t a t i c parts, we have: ph = pg(v - z)
[cosh k(D + z) ] P'~ = Pg71[ e--~shk D -1
(17)
T h e scales of Ph a n d p,~ are then: [ o (p~), o (p~)] = p9
+ D
,7
1
co2- )I,
and the ratio of the two scales is:
O(pn) H (1 0 (Ph) = 2D + H
1 ) coshkD
(19) '
As can be seen from the above equation, the ratio of t h e two pressure scales is a function of k D and H / D . Because k D = 2~rD/L, the scale of n o n h y d r o s t a t i c pressure is related to t h a t of
A Comparison
899
hydrostatic pressure by the aspect ratio and the H/D ratio. From equation (17), it can be seen that Ph has two parts. The first part is associated with ~, while the second part is proportional to the z-coordinate. Because the horizontal length scale for 7/is L, the first part of hydrostatic pressure has a horizontal length scale of L. The second part of hydrostatic pressure has an infinite horizontal length scale because it is not a function of x. Correspondingly, the hydrostatic pressure scale shown in equation (18) also has two parts. The first part is associated with the wave height, while the second part is related to the water depth. For a small amplitude wave with a wave height that is much smaller than the water depth (H << D), the second part is much larger than the second part. Therefore, we have
O(p~,)lH( O(ph)-2D
1
1 ) cosh27ra
(20) '
where c~(= D/L) is the aspect ratio. It is obvious that regardless the value of the aspect ratio, O(pn) << O(ph) is always true as long as the wave height (H) is much smaller than the water depth (D). For example, for a H/D ratio of 0.04, equation (18) suggests that O(pn) is always smaller than 2% of O(ph). For a small aspect ratio of c~ = 0.0455, O(pn) ~- O.O00790(ph) when H/D = 0.04, or the scale of nonhydrostatic pressure is more than three orders of magnitude smaller than that of hydrostatic pressure. Even for a relative large aspect ration of a = 0.5, we have O(pn) ~- O.O180(ph) when H/D = 0.04, or almost 1.74 orders of magnitude smaller than the scale of hydrostatic pressure. Now, let us look at the horizontal gradients of the two pressure components. From equation (17), we have:
(Oph Opt,) 0~7 [ coshk(D+z) Ox ' Ox = Pg"~x 1, cosh kD
] - 1
(21)
and
o\--~-z],O\ox/j=pg-~ or
0 (Op,~/Ox)0 (Oph/OX)
cosh kD cosh kD
1-
1,
1
c-~hshkZ~
1 cosh 2~ra"
'
(22)
(23)
The above equation indicates that the ratio of the two horizontal pressure gradients is a function of the aspect ratio only. Unlike equation (19), the right-hand side of equation (23) is not proportional to the H/D ratio, or the H/D ratio does not have any effect on the relative scale between the two horizontal gradients. Form equation (22), it can be seen that the horizontal gradient scale for hydrostatic pressure is only related to the wave height. The reason for the disappearance of the water depth from the expression of O ( - ~h ) is that the second part of hydrostatic pressure has a horizontal length scale that is infinite. Comparing equations (23) with equation (20), it can be seen that the right-hand side of equation (23) is about 2D/H times of that of equation (20). The ratio of the horizontal gradient of nonhydrostatic pressure over that of hydrostatic pressure approaches 1 as the water depth becomes larger (deep water). As the aspect ratio approaches zero, it approaches 0. For an aspect ratio of o = 0.0455, the right-hand side of equation (23) is about 0.04. Or, the scale of the horizontal gradient of nonhydrostatic pressure is about 1.40 orders of magnitude smaller than that of hydrostatic pressure. In this case, although nonhydrostatic pressure is more than three orders of magnitude smaller than hydrostatic pressure, their horizontal gradients are not differed by three orders of magnitude. Instead, their horizontal gradient scales are differed only by 1.40 orders of magnitude. On the other hand, for an aspect ratio of c~ -- 0.5, the right-hand side of equation (22) is about 0.91. In this case, the horizontal gradients for nonhydrostatic and hydrostatic pressure components are basically in the same order of magnitude, even when the scale of nonhydrostatic pressure is about 1.74 orders of magnitude smaller than that of hydrostatic pressure.
900
X. CHEN T a b l e 1. A c o m p a r i s o n of scales of h y d r o s t a t i c a n d n o n h y d r o s t a t i c p r e s s u r e c o m p o n e n t s as well as t h e i r h o r i z o n t a l g r a d i e n t s for H I D ---- 0.04.
(?(p,o o(pu)
o(op,~/o~) o(oph/ox)
0.0455
0.000079 < 10 - 4
0.04 ~ 10 -1"4°
0.5
0.018 = 10 -1' 74
0.91 = 10 -0.04
D L
Table 1 is a recap of the above discussion. A much smaller nonhydrostatic pressure component in comparison with the hydrostatic pressure component does not necessary guarantee its horizontal gradient to be much smaller than that of hydrostatic pressure. The scale difference is much smaller for the horizontal gradients than the scale difference between the two pressure components themselves. For a large aspect ratio of a = 0.5, the horizontal gradient of nonhydrostatic pressure is almost in the same order of magnitude as that of hydrostatic pressure without the need t h a t the scale of nonhydrostatic pressure to be in the same order of magnitude as that of hydrostatic pressure. This means that the condition that the nonhydrostatic pressure is negligible in comparison with hydrostatic pressure is not sufficient for the horizontal gradient of nonhydrostatic pressure in the horizontal momentum equation to be small enough (in comparison with the horizontal gradient of hydrostatic pressure) so that it can be taken away from the equation. For shallow water waves where a approaches zero, the right-hand side of equation (23) does approach zero, indicting that the nonhydrostatic effect in the horizontal momentum equation is indeed insignificant. In practice, if the D / L ratio is less than 0.05, the wave is treated as a shallow water wave [19] and the hydrostatic assumption is normally assumed to be valid. Nevertheless, as discussed above, for a = 0.0455, the magnitude of the horizontal gradient of nonhydrostatic pressure is only about 4% of that of hydrostatic pressure. Although the effect of nonhydrostatic pressure on the flow is indeed much smaller than that of hydrostatic pressure, such a nonhydrostatic effect cannot be neglected for an initial-value problem such as seiching in a closed basin. Model results with and without the consideration of the nonhydrostatic effect is quite different, even when the horizontal gradient of nonhydrostatic pressure is more than one order of magnitude smaller than that of hydrostatic pressure for seiching. The reason for this is that the nonhydrostatic effect is cumulative and cannot be cancelled out within the wave cycle. This can be seen from equation (21) and be explained here. Equation (21) shows t h a t both hydrostatic and nonhydrostatic pressure gradients in the horizontal direction are proportional to the gradient of surface elevation. However, their effects on the flow are opposite because they have opposite signs. The sum of the two horizontal gradients is:
Op _ Oph Op,~ cosh k(D + z) O~ Ox O~ + ~ = pg c o s h k D Ox"
(24)
If the nonhydrostatic effect on the flow is neglected, then
Op 077 cosh k(D + z) O~ O---x= Pg--~x -~ pgl cosh kD Ox' where g,
cosh kD -- gcosh k(D + z)"
(25) (26)
In other words, neglecting nonhydrostatic terms in the horizontal momentum equations is equivalent to an increase of the gravitational acceleration, which, for a = 0.0455, can be as much as 4%. From the dispersion relationship, it can be seen that an increase of the gravitational acceleration means an increase in wave frequency because the wave number k is fixed for the first mode standing wave in a closed basin. It can be proved that for a --- 0.0455, the increase
A Comparison
901
in the wave frequency using the hydrostatic assumption is about 1.5%. Although this is just a small increase, the solutions with and without nonhydrostatic effect will be quite different after a few cycles, as can be seen in Figures 6 and 7. For c~ = 0.5, the increase in wave frequency is about 77.6%. As discussed in the last section and shown in Figure 4, the wave period is reduced from 3.59 to 2.02. The solutions with and without the nonhydrostatic effect will be significantly different even in the first wave cycle. Errors induced due to the exclusion of nonhydrostatic pressure in one oscillation cycle will affect model results in the subsequent cycles.
5. C O N C L U S I O N S The hydrostatic pressure assumption has been widely used in many numerical models for free surface flows in shallow waters, because the assumption allows governing equations for the free-surface flows to be simplified. However, there are many cases where nonhydrostatic effects should be taken into account. A general rule of thumb for a safe use of the hydrostatic pressure assumption in a free-surface flow is that the length scale in the vertical direction should be much smaller that that in the horizontal direction, so that the nonhydrostatic pressure component is negligible in comparison with the hydrostatic pressure component. Nevertheless, there exist no quantitative guideline as to how small the aspect ratio should be for the nonhydrostatic effect to be negligible. In literature, it is generally believe that when the aspect ratio is less than 0.05, it is a shallow water flow and the nonhydrostatic effect can be neglected. For a better understanding of the relative magnitude of nonhydrostatic pressure in comparison with hydrostatic pressure, this study compares magnitudes of hydrostatic and nonhydrostatic pressure components and their horizontal gradients in two seiching cases. The first one is a deepwater wave oscillation with an aspect ratio of 0.5, while the second is a shallow water wave oscillation with an aspect ratio of 0.0455. Magnitudes and scales of hydrostatic and nonhydrostatic pressure components were compared and analyzed using model results and analytical solutions. Both model results and analytical solutions have shown that the nonhydrostatic component is much smaller than the hydrostatic component for small amplitude waves, even when the vertical length scale is comparable to the horizontal length scale. On the other hand, the horizontal gradient of the nonhydrostatic component can be in the same order of magnitude as that of hydrostatic pressure, even when nonhydrostatic pressure is much smaller than hydrostatic pressure. This kind of an "amplification" of the nonhydrostatic effect is due to the fact that hydrostatic pressure has a horizontal length scale that is infinite. The scale analysis shows that nonhydrostatic pressure is much smaller than hydrodynamic pressure because the ratio of the two scales is not only a function of the aspect ratio, but also proportional to the ratio of wave height to the water depth (H/D), which is normally small for small amplitude waves. However, once the horizontal derivative is taken, the ratio of the two horizontal gradient scales is no longer proportional to the D/H ratio. Comparison of model results of hydrostatic simulation with those of nonhydrostatic simulation reveals the differences between the two. Figures 2 and 3 and Figures 6 and 7 have shown that simulated velocities and surface elevations with and without the hydrostatic pressure assumption are totally different. Model results indicate that the hydrostatic assumption tends to produce higher mode oscillations and forces the first mode oscillation to have a phase speed of x/~, even for a deep-water wave.
In summary, this study has shown that the condition [Pnl << Ph does not guarantee that [°0-~I is smaller than [-~z~ ] by the same order of magnitude. Even when ~ is in the same order of magnitude as °o-~zh, the nonhydrostatic pressure component is still about two orders of magnitude smaller than the hydrostatic pressure component. As such, the condition that °0-~is comparable to ~ does not necessarily require ]Pnl ~ Ph. In other words, ]p~] ~ Ph is not the necessary condition for ~ and Oph to be comparable.
902
X. CHEN
REFERENCES 1. Y.P. Sheng, A three-dimensional mathematical model of coastal, estuarine and lake currents using boundary fitted grid, Report No. 585, A.R.A.P. Group of Titan Systems, Princeton, NJ, (1986). 2. A.F. Blumberg and C.L. Mellor A description of a three-dimensional coastal ocean circulation model, In Three-Dimensional Coastal Ocean Models, (Edited by N. S Heaps), pp. 1-16, American Geophysical Union, Washington, DC, (1987). 3. V. Casulli and R.T. Cheng, Semi-implicit finite difference methods for three-dimensional shallow water flow, Int. J. for Numer. Methods in Fluids 15, 629-648, (1992). 4. J.M. Hamrick, A three-dimensional environmental fluid dynamic computer code: Theoretical and computational aspects, Special Report in Applied Marine Science and Ocean Engineering No. 317, Virginia Institute of Marine Science, College of Willian 8z Mary, Gloucester Point, VA~ (1992). 5. X. Chen, A flee-surface correction method for simulating shallow water flows, Journal of Computational Physics 189, 557-578, (2003). 6. V. Casulli and G.S. Stelling, Numerical simulation of 3D quasi-hydrostatic free-surface flows, Journal of Hydraulic Engineering 124 (7), 678-686, (1998). 7. X. Chen, Three-dimensional, hydrostatic and non-hydrostatic modeling of seiching in a rectangular basin, In Estuarine and Coastal Modeling, Proe. of the 6th Intern. Conf., (Edited by M.L. Spaulding and H.L Butler), pp. 148-161, ASCE, New Orleans, LA, (1999). 8. X. Chen, A fully hydrodynamic model for three-dimensional, free-surface flows, Int. J. for Numer. Methods in Fluids 42, 929-952, (2003). 9. H. Lamb, Hydrodynamics, Sixth Edition, Dover Publications, New York, (1932). 10. J. Pedlosky, Geophysical Fluid Dynamics, Springer-Verlag, New York, (1979). 11. G.H. Golub and C. vanLoan, Matrix Computations, John Hopkins University Press, Baltimore, MD, (1990). 12. H.A. vanderVorst, Bi-CGSTAB: A fast and smoothly converging variant of bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13 (2), 631-644, (1992). 13. J. Smagorinsky, General circulation experiments with primitive equations, I. The basic experiment, Monthly Weather Review 91 (3), 99-164, (1963). 14. Y.P. Sheng and C. Villaret, Modeling the effect of suspended sediment stratification on bottom exchange process, J. Geophys. Res. 94, 14429-14444, (1989). 15. X. Chen, Effects of Hydrodynamics and Sediment Transport Processes on Nutrient Dynamics in Shallow Lakes and Estuaries, Ph.D. Dissertation, University of Florida, Gainesville, FL, (1994). 16. X. Chen, Responses of a hybrid z-level model to various topography treatments for a boundary-value problem and an initial-value problem, In Estuarine and Coastal Modeling, Proc. of the 7eh International Conference, (Edited by M.L. Spaulding), pp. 614-627, ASCE, St. Pete Beach, FL, (2001). 17. X. Chen, Fitting topography and shorelines in a 3-D, Cartesian-grid model for free-surface flows, In Computational Fluid and Solid Mechanics 2003, Volume 2, (Edited by K.J. Bathe), pp. 1892-1895, Elsevier Science Ltd., (2003). 18. X. Chen, Using a piecewise linear bottom to fit the bed variation in a laterally averaged, Z-coordinate hydrodynamic model, Int. J. for Numer. Methods in Fluids 44, 1185-1205, (2004). 19. R.G. Dean and R.A. Dalrymple, Water Wave Mechanics for Engineers and Scientists, Princeton-Hall, Inc., Englewood Cliffs, NJ, (1984).
Available online at www.sciencedirect.com
8CIENCl~DIRECTe
ELSEVIER
Mathematical and Computer Modelling 41 (2005) 887-902 www.elsevier.com/locate/mcm
A Comparison of Hydrostatic and Nonhydrostatic Pressure Components in Seiche Oscillations XINJIAN CHEN Resource Conservation and Development Department Southwest Florida Water Management District 7601 Highway 301 N o r t h , T a m p a , F L 33637, U.S.A. xinj ian. chen©swfwmd, state, fl. us
(Received April 2004; revised and and accepted August 2004)
A b s t r a c t - - T h e hydrostatic pressure assumption has been widely used in studying water movements in shallow waters where the water d e p t h (D) is much smaller t h a n the horizontal length scale (L). It is assumed t h a t the pressure field in a shallow water flow is hydrostatic because the nonhydrostatic pressure component (Pn) is much smaller t h a n the hydrostatic pressure component (Ph) and thus be safely neglected in model simulations. W i t h this assumption, a three-dimensional, hydrodynamic model can be significantly simplified because the m o m e n t u m equation in the vertical direction is reduced to an expression of hydrostatic pressure. While such a simplification is valid in many shallow water applications, there exist no rigorous guidelines as to what extent can the hydrostatic assumption be safely used without any problems. In an effort to get a better understanding of the effect of nonhydrostatic pressure on free-surface flows, this paper presents a comparison of hydrostatic and nonhydrostatic pressure components and their horizontal gradients in two seiches. T h e first seiche represents a deepwater wave oscillation, while the second one is a shallow water wave oscillation. A nonhydrostatic model was used to simulate the two oscillations, and simulated hydrostatic and nonhydrostatic pressure components were compared. A scale analysis was also carried out for the two pressure components in seiching. It is found t h a t although the nonhydrostatic component was about 2 to 3 orders of magnitude smaller t h a n the hydrostatic component for the two oscillations, the horizontal gradient of the nonhydrostatic component is not necessarily smaller t h a n t h a t of hydrostatic pressure by 2 to 3 orders of magnitude. On the other hand, while the horizontal gradient of nonhydrostatic pressure may b e in the same order of magnitude as t h a t of hydrostatic pressure, it does not necessarily require nonhydrostatic pressure itself to be in the same order of magnitude as hydrostatic pressure. Model simulations also show t h a t the nonhydrostatic portion of the horizontal pressure gradient affects the velocity field during the course of the oscillation and its effects are not negligible even when t h e aspect ratio is less t h a n 0.05. ~) 2005 Elsevier Ltd. All rights reserved.
Keywords--Hydrostatic pressure, Nonhydrostatic pressure, Hydrostatic pressure assumption, Seiche oscillation, Nonhydrostatic model.
1. I N T R O D U C T I O N The hydrostatic pressure assumption has been widely used in studying shallow rivers, lakes, estuaries, and continental shelves. This assumption is valid in most cases and has been successfully 0895-7177/05/$ - see front m a t t e r (~) 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2004.08.005
Typeset by ~4A/tS-~I~_tX
888
X. CHEN
used in many three-dimensional, hydrodynamic models (e.g., [1-5]). However, there are also many cases where this assumption may be questionable [6-8]. One example is the flow induced by a very large horizontal density gradient. Other examples include short waves, flows through structures, the near field of brine disposal from a desalination plant, and flows near an uptake point of a water withdrawal system. In these cases, nonhydrostatic pressure effects may be comparable to those of hydrostatic pressure, and thus cannot be neglected in model simulations. The basic idea of the hydrostatic pressure assumption comes from the fact that the acceleration and eddy viscosity terms in the momentum equation for the vertical velocity component are much smaller than the gravitational acceleration term and can thus be neglected [9,10], because vertical length scales in shallow waters are general much smaller than horizontal length scales. The momentum equations for water movement in three directions can be expressed as follows
p-~x + -~x Ah
d--[ =
dv dt -
+ ~Y k, h ~y J -Jc--~Z \
fu_lOp 0 (AhOV h 0 AhOy ~z\ p -Oyy+ -Oxx\ OxJ + -~y ÷
dw
lop
dt-g-
+Tx
0 fA Ow~
hTx)+
0 fAhOW~
(1)
OZ / '
Oz]
,
(2)
0 fAvOW~
\
(3)
az]'
where x, y, and z are Cartesian coordinates (x is from west to east, y is from south to north, and z is vertical pointing upward); u, v, and w are velocities in the x-, y-, and, z-directions, respectively; t, f, g, and p denote time, the Coriolis parameter, the gravitational acceleration, and pressure, respectively; Ah and Av represent horizontal and vertical eddy viscosities, respectively; and p is density which is a function of temperature and salinity. Equation (3) can be rewritten in the following form
lop
dw
O fAhOW~
p Oz -- g - --~ + -~x \
O fAhOW~
Ox ] + ~y \
0 fAvOW~
Oy ] + -~z \
Oz ] = - g + a3,
(4)
where a3 includes the vertical acceleration of a water particle and the vertical viscosity term. Integrating the above equation from an elevation z to the free surface (~7), we have P - Pa = g
//
p d~-
S/
pa3 de,
(5)
where Pa represents pressure at the water surface (~) and is assumed to be zero throughout this paper. The first integration on the right-hand side of equation (5) is hydrostatic pressure (Ph) and the second integration is nonhydrostatie pressure (P,0. In the hydrostatic pressure assumption, the magnitude of the second integral on the righthand side of the above equation is assumed to be much smaller than the first one because the aspect ratio is very small for shallow waters. Therefore, equation (5) is reduced to p= g
//
p de.
(6)
Obviously, if O(a3) << g, equation (6) is always valid. By replacing total pressure in the horizontal momentum equations with hydrostatic pressure, the number of governing equations for three-dimensional shallow water flows is reduced from four to three: two horizontal momentum equations and a continuity equation
d-~ -= dt --
- pox P-~Y
pd¢ + ~x \ pd¢ + ~ \
Ox) +-~y \ Ox/#+-~y \
ay) + Oz \
Oz]'
Oy] +-~z A . ~
(7)
,
(8)
A Comparison
Ou
Ov
889
Ow
0---~-b ~yy -F ~ z -- 0.
(9)
Many three-dimensional hydrodynamic models (e.g., [1-5]) solve the above three equations. However, this kind of treatment of the horizontal pressure gradient term is mathematical not rigorous because it implies that if the nonhydrostatic portion of pressure is much smaller than the hydrostatic portion, then the horizontal gradients of nonhydrostatic pressure would also be much smaller than those of hydrostatic pressure. In other words, if IP~ [ << Ph, then the following expression would be automatically true:
Opn 0x
OOph
Opn
Oy
<<
Oph
(10)
<<
The above logic is only true when both p~ and Ph have the same length scale in the horizontal direction. However, because hydrostatic pressure may have multiple length scales for shallow water flows, the above logic is not necessarily always correct. There could be cases where ]p~ [<
H
T
cos
(~rx)
(11)
'
where H is the wave height of the oscillation and l~ is the length of the rectangular basin. At time > 0, water in the basin undergoes oscillation under the barotropic pressure gradient caused by the initial free surface setup. In the following sections, a three-dimensional, nonhydrostatic model is briefly described, before model applications to the two oscillations, along with a detailed comparison of hydrostatic and nonhydrostatic pressure components and their horizontal gradients, are presented. Model results with and without the hydrostatic pressure assumption during the oscillations are also discussed. A scale analysis using analytical solutions for small amplitude waves is then described to show Z
D
W
X F i g u r e 1. A seiche oscillation in a r e c t a n g u l a r b a s i n is c a u s e d by a n initial free-surface setup.
890
X. CHEN
the relative scales of nonhydrostatic pressure in comparison with those of hydrostatic pressure. Conclusions drawn from the study are included in the last section of the paper.
OF
A
2. A B R I E F DESCRIPTION NONHYDRODYNAMIC MODEL
A three-dimensional nonhydrodynamic model named LESS3D (Lake and Estuary Simulation System in Three Dimensions) [7,8] was used in this study. The model solves the governing equations for free-surface flows, including equations (1)-(3) and equation (9), by performing two predictor-corrector steps. In the first predictor-corrector step, the model uses hydrostatic pressure at the previous time step as an initial estimate of the total pressure field at the new time step. Based on the estimated pressure field, an intermediate velocity field is calculated, which is then corrected by adding the nonhydrostatic component of the pressure to the estimated pressure field. By forcing the velocity field divergence-free, a Poisson equation for nonhydrostatic pressure is obtained. In a structured grid system, the Poisson equation forms a seven-diagonal matrix system that is positive definite if a uniformly distributed grid system is used. In this case, the conjugate gradient method with incomplete Cholesky preconditioning [11] is used to find the inverse of the matrix. If the grid sizes are not uniformly distributed in the computation domain, the sevendiagonal matrix is generally asymmetric. In this case, a Bi-CGSTAB method [12] is used to solve the matrix. After the nonhydrostatic pressure component is solved, the intermediate velocity field is corrected to obtain the second intermediate velocity field that is the second prediction of the final velocity field. In the second predictor-corrector step, an intermediate free-surface location at the new time step is calculated based on the second intermediate velocity field. The intermediate free surface is then adjusted to find the final free-surface location based on a solution of a free-surface correction equation. The final velocity field is finally solved after the final freesurface location is found. A detailed description of the numerical procedure for implementing the two predictor-corrector steps in the nonhydrostatic model LESS3D is presented in reference [8]. After finding the free-surface location and the velocity field at the new time step, the model solves the following transport equation of concentration oc
ouc
owc
o
(B
O---~+ --.~x + --~y + O---~.-= O-.-.~\
o
o
(Boac
Ox j + -~y ~, h ~y / + -~z \
Oz ) + S + R,
(12)
where c is concentration (can be salinity, temperature, suspended sediment concentration, nutrient species, etc.), S denotes source/sink terms, R represents reaction terms, and Bh and B, are horizontal and vertical eddy d/flus/v/ties, respectively. If the simulated concentration involves settling, w in the above equation includes the settling velocity. In this study, however, because the seiche oscillations do not involve concentration, equation (12) is not solved. The model calculates eddy viscosity and diffusivity in the horizontal directions (Ah and Bh) using Smagorinsky's approach [13]. For the vertical eddy viscosity and diffusivity (Av and B~), a turbulent kinetic energy (TKE) model [14,15] is used. In this study, however, a small constant eddy viscosity of 0.1 cm 2 sec-1 was used in the simulation. Boundary conditions in the horizontal directions are specified with either free surface elevations or velocities for open boundaries. At solid boundary, both the velocity component and the pressure gradient in the normal direction are set to zero. Boundary conditions specified in the vertical direction are shear stresses. At the free surface, wind shear stresses are used. At the bottom, the log-layer distribution of velocity is used to calculate the bottom shear stresses, or =
+
(13)
where a is the vonKarman constant (0.41), Ub and Vb are horizontal velocity calculated at a level zb near the bottom, zo --- ks~30 and k~ is the bottom roughness.
A Comparison
891
The LESS3D model can be run either with or without the nonhydrostatic effect. When the model is run without the nonhydrostatic effect, the model skips the substep of finding the inverse of the seven-diagonal matrix for nonhydrostatic pressure and the substep of correcting the first intermediate velocity field. The model uses the so-cMled z-level in the vertical direction and no transformation is made to the vertical coordinate. The model allows the free surface to travel between layers, so that there is no need to make the top layer thick enough to cover the whole range of free surface variation. To fit the bottom topography, hybrid grid cells are used in the model [16,17]. 3. N U M E R I C A L TWO SEICHE
SIMULATION OF OSCILLATIONS
This section presents applications of the LESS3D model to the two seiche oscillations described in the Introduction section. As mentioned before, the basin length and the water depth are the same in the first seiche oscillation and both equal to 10 m. In the second oscillation, however, the basin length is 44m, which is 11 times of the water depth of 4m. Both oscillations have Hydrodynamic
Pressure
20
A
E 15 o~ t- 10 o
5
U,I
-5
0 o -10 M 1= -15 i
=
-20 0 . . . . . .
......
......
a' . . . . . . . . . . . . .4. . . . . . .
s
.
.
.
.
.
.
.
.
.
.
.
.
6
; ......
; ......
10
Time (second) 20 Ix = O.lm, z f S.gml
15
A
o
5
0 0
0
-5
~-10 -lS .200
A
1
2
3
4
~ ......
5 Time (second)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7, . . . . . . . . . . 8. . . . . . . . . . 9
10
5O 4(1
10 0
-10 -20 -30 -4O -50
. . . . . .
0
I
1
. . . . . .
I
2
,
,
r
,
,
,
I
3
,
,
,
,
,
,
I
4
,
,
,
,
,
,
I
,
,
,
,
i
5 Tlme (second)
,
I
6
,
,
,
i
,
,
I
7
,
,
r
i
,
,
I
8
. . . . . .
I
. . . . . .
9
F i g u r e 2. C o m p a r i s o n of s i m u l a t e d ( s o l i d l i n e s ) s u r f a c e e l e v a t i o n s a n d v e l o c i t i e s w i t h a n a l y t i c a l s o l u t i o n s ( d a s h e d l i n e s ) for t h e first s e i c h e o s c i l l a t i o n w i t h a : 0.5. S u r f a c e e l e v a t i o n s a r e c o m p a r e d a t x = 0.1 m a n d 5.5 m , w h i l e v e l o c i t i e s a r e c o m p a r e d a t x = 0 . 1 m , z = 8 . 9 m a n d x = 5 . 5 m , z = 4 . 9 m . T h e n o n h y d r o s t a t i c effect w a s included in the model simulation.
10
892
X. CHEN Hydrostatic
20 A
Pressure
I/I/Ill'"
15 5 0
-5 -10 -15 -20
1
2
3
4
5 Tlme (second)
1
2
3
4
5 Time (second)
6
7
8
9
10
20 15 10 5 0 0
>:
0
-5 -10 -15 -20 0
10
.--
20
10
," "-'.
0
~-I0 -20
.20
.................. 0
I
. . . .
2
3
4
Ir
5 Tlme (second)
. . . . . . . . . . . . . . . . . . . . . .
6
7
8
9
10
Figure 3. Comparison of simulated (solid lines) surface elevations and velocities with analytical solutions (dashed lines) for the first seiche oscillation with a ----0.5. Surface elevations are compared at x = 0.1 m and 5.5 m, while velocities are compared at x ----0.1m, z ----8.9m and x = 5.5m, z ----4.9m. The nonhydrostatic effect was not included in the model simulation. the same wave height of 30 cm. To c o m p a r e the m a g n i t u d e of n o n h y d r o s t a t i c pressure and its horizontal gradient with those of hydrostatic pressure, b o t h pressure c o m p o n e n t s were saved as model o u t p u t s during the model run. Horizontal gradients of h y d r o s t a t i c and n o n h y d r o s t a t i c pressure c o m p o n e n t s were calculated using a post-process p r o g r a m after t h e simulation was done. W i t h the initial free-surface setup expressed in equation (11), a s t a n d i n g wave is formed in the rectangular basin with a wavelength t h a t is twice the length of t h e basin (L = 2/~). Obviously, the horizontal length scale of the oscillation is the wavelength (L) and the vertical length scale is t h e water d e p t h D. Thus, the aspect ratio (~) equals D / L . I n the first oscillation, because the water d e p t h is the same as the basin length, the D / L ratio is 0.5. It is expected t h a t the n o n h y d r o s t a t i c effect is significant for seiching in this kind of basin. In the second oscillation, the aspect ratio ~ is relatively small and has a value of 0.0455. In b o t h seiche oscillations, the wave slopes are m u c h smaller t h a n 1. As such, analytical solutions according to the linear wave t h e o r y can be o b t a i n e d as follows H ~- ~- cos kx c o s a t ,
(14)
(u,w) - H ~ s i n ~ t ( s i n k x c o s h k z , - coskxsinhkz) (15) 2 sinh kD where k --- 2~r/L -~ ~/l~ and a is related to k by the dispersion relationship cr2 --- g k t a n h kD.
A Comparison
893
Model Simulation of the First Seiche Oscillation
When the first oscillation was simulated using the LESS3D model, a uniform grid system was used with a grid size of 0.2 m in both the horizontal and vertical directions. The model was run for 10 s with a time step of 0.05 s. From the dispersion relationship, the wave period can be calculated and equals to 3.59 s. Figure 2 show the model results of surface elevations and velocities at two locations in the basin, together with the analytical solutions. One location is near the left end of the basin with the coordinates of x = 0.1 m and z = 8.9 m, while the other location is near the middle of the basin at x = 5.5m and z = 4.9m. The model was run with the consideration of nonhydrostatic pressure. The model results are plotted with solid lines and the analytical solutions axe in dashed lines. As can be seen from the figure, the model was able to reproduce the analytical solution fairly well. If the nonhydrostatic pressure effects were not included in the model simulation, the simulated surface elevation and velocity field would be totally different from analytical solutions. Figure 3 shows model results (solid lines) without the consideration of nonhydrostatic effect. Analytical solutions are also shown in the figure (dashed lines). Obviously, because the nonhydrostatic effect is significant in this oscillation case, the hydrostatic pressure assumption is not valid and the use of this assumption produces incorrect model results. Two things can be noticed from Figure 3. One is that the simulated oscillation has a much shorter wave period than the analytical solution does. The other thing is that higher mode oscillations are introduced to the model results shortly after the initial free surface setup. It is clear that forcing the horizontal pressure gradient to
100000
-9000 -7000 -5000
-3000 - 1 ~
111110 3000
5000
7000
9000
(b)
(~) .~.~:%~~.~.~i:~~::~!,~' ~
Ip,~pd
~!:il:~;I!~I~::,
dp./dx o.ooo o.oos O.OLO o.o~5 o.o-2o 0.025
II
I:::~:~1:=:::!~lii~:'~=]i~i~~ m~ ~ -~ -~ o
-ro-~e.~e-~4-~2-1o-~
(d)
(c)
ldpJdxl/IcWJdxl
dp~/dx 0.0
4.0
8.0
(e)
12.0
16.0
20.0
I I
t::
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5,0
(f)
F i g u r e 4. M o d e l r e s u l t s at t i m e = 6 s for t h e first seiche oscillation w i t h (a), (b), a n d (c) a r e s i m u l a t e d n o n h y d r o s t a t i c p r e s s u r e (p~), h y d r o s t a t i c (Ph), a n d t h e a b s o l u t e value of t h e ratio o f p n over Ph, respectively. (d) a n d horizontal g r a d i e n t s of p n a n d Ph, respectively. (f) is t h e a b s o l u t e v a l u e of of t h e two g r a d i e n t s .
(~ ~ 0.5. pressure (e) s h o w t h e ratio
894
X. CHEN 100000
100000--
~00c : :::::~::i~::i:i: ~:iti::: ::i::: ::: :i::i:i:::::j:::!i ~ :
p~ I]
I I
-9000 -7000 -5000 -3000 -1000
1000
3000
5000
7000
9000
1.0E+05
3.0E+OS
(a)
IpJpd
ill:::; ::
0.000
0.010
0.020
5.0E+05
7,0E+05
g.OE+05
(b)
0.030
dpr/dX
0.040
-17 -13
-9
(c)
-5
-1
3
7
11
1.5
2,0
15
19
23
(d)
I~ I ~
llit t
i~!~:il H~l:~i[
I
itll ~i:l!!l~i[ I
dp~dx
i i
I -27.0
Idp./dxl / Idpjdxl I I -21.0
-15.0
-9.0
-3.0
3.0
0.0
0.5
1,0
2.5
3.0
3.5
(e) (f) Figure 5. Model results at time = 7.7s for the first seiche oscillation with a -- 0.5. (a), (b), and (c) are simulated nonhydrostatic pressure (pn), hydrostatic pressure (Ph), and the absolute value of the ratio ofp~ over Ph, respectively. (d) and (e) show horizontal gradients of pn and Ph, respectively. (f) is the absolute value of the ratio of the two gradients. be the horizontal gradient of hydrostatic pressure is equivalent to forcing the phase speed of the wave to be V ~ , which corresponds to a wave period of 2.02 s for this seiching case. This is the same wave period shown in the simulated oscillation in Figure 3, which is much shorter t h a n the analytical solution of 3.59 s. The model simulation with the consideration of nonhydrostatic pressure allows both hydrostatic pressure and nonhydrostatic pressure as well as their horizontal gradients to be output and compared. Figures 4 and 5 show contour plots of simulated hydrostatic pressure and nonhydrostatic pressure at time = 6s and 7.7s, respectively, in the first oscillation. In Figures 4 and 5, (a) and (b) show simulated nonhydrostatic pressure (Pn) and hydrostatic pressure (Ph), respectively. Their corresponding horizontal gradients are shown in (d) and (e), respectively. Figures 4c and 5c are plots of the absolute value of the ratio of p~ over Ph, while Figures 4f and 5f are plots showing the absolute value of the ratio of the horizontal gradients of the two pressure components. As expected, contours of hydrostatic pressure are parallel to the free surface because there is no baroclinic effect in the oscillation. The magnitude of hydrostatic pressure varies linearly from 0dyn-cm -2 at the free surface to about 1 × 106 dyn-cm -2 near the b o t t o m . On the other hand, the nonhydrostatic pressure distribution in the basin is more complex. It can be either negative or positive, with the magnitude being as large as 9,000dyn-cm -2 during the course of the oscillation. Because the horizontal gradient of hydrostatic pressure is not a function of z, its contours are vertical straight lines. Because the horizontal gradient of the free surface approaches zero near the two ends of the basin, the horizontal gradient of hydrostatic pressure also approaches zero near
A Comparison
895
the two ends. The horizontal gradient of nonhydrostatic pressure varies not only with the time and the x-coordinate, but also with the z-coordinate. The magnitude of the horizontal gradient of nonhydrostatic pressure is almost the same as that of the horizontal gradient of hydrostatic pressure. From Figures 4 and 5, it can be seen that although nonhydrostatic pressure is about two orders of magnitude smaller than hydrostatic pressure, their horizontal gradient are comparable in magnitude for this seiching where the aspect ratio is 0.5. M o d e l S i m u l a t i o n of t h e Second Seiche Oscillation For the second seiche oscillation, a uniform grid spacing of 0.5 m was used in the horizontal directions, while the vertical grid spacing varies between 10 cm and 30 cm. The time step used was 0.1 s. From the dispersion relationship, the phase speed of the wave is 6.18m/s and the period of the first mode oscillation is about 14.24s. The phase speed calculated from the dispersion relationship is about 1.3 percent smaller than that calculated from ~ / ~ , which equals 6.26 m/s and would be the phase speed under the hydrostatic assumption.
Hydrodynamic Pressure
25 ,..2O
E 15 o
g o ,-n -5
~ -10 ~ -15
m 2O 20
411
[x = 32.25m,
30
A
30
z = 2.85m
40
70 80 Time (second)
90
100
11 0
120
130
140
I
g2o lO
_~
o
¢)
2 -10 Q
>,.20
I~-- o ys~, ~ = o . .
-30
-40
. . . . . .
I
0
. . . . . .
I
10
. . . . . .
I
20
i
. . . . . .
30
I
. . . . . .
[
40
. . . . . .
i
50
. . . . . .
i
,
,
,
,
,
,
i
. . . . . .
60 70 80 Time (second)
i
. . . . . .
90
i
. . . . . .
100
i
. . . . . .
11 0
i
. . . . . .
120
i
. . . . . .
130
140
10
~5
Ix = 0.75m, z = 0 . 9 5 m j
o
~o 4}
-5 = 32.z5,,,, ~ = 2.85,,, I mlO
. . . . . .
0
I
. . . . . .
10
I
. . . . . .
20
I
. . . . .
30
,
~ I
40
. . . . .
I
. . . . . .
50
i
. . . . . .
I
. . . . . .
i
,
60 70 80 Time (second)
[
,
,
i
, 1
.....
90
r I
. . . . . .
100
I
. . . . . .
110
I
. . . . . .
120
I
. . . . .
130
F i g u r e 6. C o m p a r i s o n of s i m u l a t e d (solid lines) s urfa c e e l e v a t i o n s a n d v e l o c i t i e s w i t h a n a l y t i c a l s o l u t i o n s ( d a s h e d lines) for t h e second seiche o s c i l l a t i o n w i t h a :- 0.0455. Surface e l e v a t i o n s are c o m p a r e d at x = 0.75 m a n d 32.25 m, w h i l e v e l o c i t i e s axe compared at x : 0.75m, z : 0.95m and x = 32.25m, z = 2.85m. The nonhydrostatic effect w a s i n c l u d e d in t h e m o d e l s i m u l a t i o n .
140
X. CHEN
896
Hydrostatic Pressure 25 ..20 E 15 10 = 5 _~ 0
i -5
.20 10
20
30
60 70 80 Time (second)
90
100
110
40
40
30
[~=32.3s~z--z.ssm[
lO
.~,
o
_~-10 >-20 -3o
I ~ = o.Tsm, ~ = o . . m
-40 0 ...... 10'...... 20' ...... 30' .....
~
I
..... 50'"
....60' ...... 70I...... 80'...... 90'...... 100l...... 11' 0...... 120'...... 130'...... 140 Time (second)
10 x = 0.75m, z -- 0 . 9 5 m 1
/h 0
A
: -10
"
-
"
'
'
M " &' '
-
Ix = 32.2s.,
0 . . . . 1'0. . . . 2'0. . . . ~
....
40 . . . .
= 2.SSm I ~
.... 6
.... iO .... ~
.... ~
.... i~
.... 110 .... li0 .... i30 ....
MO
Time (second) F i g u r e 7. C o m p a r i s o n of s i m u l a t e d (solid lines) s urfa c e e l e v a t i o n s a n d v e l o c i t i e s w i t h a n a l y t i c a l s o l u t i o n s ( d a s h e d lines) for t h e s e c ond seiche o s c i l l a t i o n w i t h ~ = 0.0455. Surface e l e v a t i o n s are c o m p a r e d at x = 0.75 m a n d 32.25 m, w h i l e v e l o c i t i e s a re comp a r e d a t x = 0 . 7 5 m , z = 0 . 9 5 m a n d x = 3 2 . 2 5 m , z ---- 2 . 8 5 m . T h e n o n h y d r o s t a t i c effect was not i n c l u d e d in t h e m o d e l s i m u l a t i o n .
The model was run for 142 s (about ten cycles of the first mode oscillation) with and without the hydrostatic pressure assumption. For the nonhydrostatic simulation, both hydrostatic and nonhydrostatic pressure fields were output. Similar to Figures 2 and 3, simulated time series of surface elevation and velocities at two locations are compared with analytical solutions in Figures 6 and 7 for simulations with and without the consideration of the nonhydrostatic effects, respectively. In the figures, solid lines are model results and dashed lines are analytical solutions. The two locations are at (x, z) -- (0.75 m, 0.95 m) and (x, z) = (32.25 m, 2.85 m). It can be seen from Figures 6 and 7 that both nonhydrostatic and hydrostatic simulations yield similar model results during the first couple of cycles because the nonhydrostatic effect is not significant at the beginning of the seiche oscillation case. However, because the error caused by the exclusion of the nonhydrostatic effect is cumulative for this initial value problem, model results of the hydrostatic simulation start to stray away from the first mode analytical solution as time goes on, while model results of the nonhydrostatic simulation still have a reasonable match with the analytical solution. Near the end of the simulation period, it is evident that model results without the nonhydrostatic effect possess a phase lead with respect to the analytical solution as the cumulative effect of the slightly larger phase speed under the hydrostatic assumption becomes significant. Clearly, model
A Comparison
~,~,~, ,~
,,~
897
' ~ff~ ~ ~7~7~#,7~i~7'>".~i~:~i~ '~e
:<~
i7~i77i:~:ii!i
p. I I -400
I -300
I......... ~ . ~ , ~ / ~ ~ -200
-100
0
100
200
P~
300
II
!:::: : : f i~t7q 7 ! : : ' i % 1 7 Z ~ /
O.OE+O0
1.0E+05
(a)
I[
0
I: I
0,001
:I,:~. l 0.002
3.0E+05
4.0E+05
(b/
~
IpdpJ
2.0E+05
h ~
0.003
~
¢
!
dp,,/dx
0,004
~
[
-0.4
i
:
-0,3
~
-0.2
~"
-0.1
0.0
~ ::~i~:~ ]
0.1
0.2
(d)
(c)
.......:::;, ::: :i.:: !!4~ : ' ~!~i!'~
dp~/dx
0.0
1.0
i
2.0
3.0
4,0
5.0
O.Q5rL 7! "!
O.(X) 0.05 0.10 0,15 0.20 0.25 0.30 0.35 0.40
6.0
(e)
(0
F i g u r e 8. M o d e l r e s u l t s at t i m e = 130.4s for t h e second seiche oscillation w i t h = 0.0455. (a), (b), a n d (c) are s i m u l a t e d n o n h y d r o s t a t i c p r e s s u r e (p~), h y d r o s t a t i c p r e s s u r e (Ph), a n d t h e a b s o l u t e value of t h e ratio ofpn over Ph, respectively. (d) a n d (e) s h o w horizontal g r a d i e n t s of pn a n d Pu, respectively. (f) is t h e a b s o l u t e value of t h e ratio of t h e two gradients.
results with the nonhydrostatic effect have a much better agreement with analytical solutions than those without the nonhydrostatic pressure effect, especially after the first couple of wave cycles. Figures 8 and 9 are simulated hydrostatic and nonhydrostatic distributions at time = 130.4s and time -- 137.6 s, respectively. Similar to Figures 4 and 5, (a) and (b) in Figures 8 and 9 show nonhydrostatic pressure and hydrostatic pressure, respectively, while (d) and (e) are plots of their horizontal gradients. The Pn over Ph ratio and the °o~ over ~ ratio are plotted in (c) and (f) of the two figures, respectively. From Figures 8 and 9, it can be seen that nonhydrostatic pressure is at least two orders of magnitude smaller than hydrostatic pressure for the second seiching case that has an aspect ratio of a - 0.0455. Again, the horizontal gradient of hydrostatic pressure is not a function of z, while the horizontal gradient of nonhydrostatic pressure varies with both x- and z-coordinates. Obviously, the horizontal gradient of nonhydrostatic pressure is no longer two orders of magnitude smaller than that of hydrostatic pressure. It is only about one order of magnitude smaller than the horizontal gradient of hydrostatic pressure.
4. A
SCALE
STUDY
To get a better understanding of the model results presented in the last section, a scale analysm for hydrostatic and nonhydrostatic pressure components and their horizontal gradients was conducted using analytical solutions expressed in equations (14) and (15). From the unsteady
898
X. CHEN
: ::: ~:: :: : ~ h ~H~
~ "
-400
-
-
~
-300
¥
:
~
-200
.~
-100
.
0
.
.
100
.
200
~ . ~
Ph
300
O.OS+O0
1.0E+05
(a)
IP~/Phl
II0
0.001
0.002
2.0E+05
3,0E+05
4.0E,P05
(b)
0.003
-0.1
0.004
-0,0
(c)
0.1
0.2
0.3
0.4
0.5
(d)
f: :i::i !:!:i.li:~'~' Idp~/dxl t IdpJdxl -5.0
-4.0
-3.0
-2.0
-1.0
0.00 0.05 0,10 0.15 0.20 0.25 0.30 0.35 0.40
0.0
(f)
(e)
Figure 9. Model results at time = 137.6s for the second seiche oscillation with cz = 0.0455. (a), (b), and (c) are simulated nonhydrostatic pressure (pn), hydrostatic pressure (Ph), and the absolute value of the ratio of p~ over Ph, respectively. (d) and (e) show horizontal gradients of p~ and Ph, respectively. (f) is the absolute value of the ratio of the two gradients. Bernoulli equation, the total pressure can be o b t a i n e d and takes the following form [19]:
p = - p g z + pg~
cosh k(D + z) cosh k D
(16)
Splitting the pressure expressed in the above equation into h y d r o s t a t i c and n o n h y d r o s t a t i c parts, we have: ph = pg(v - z)
[cosh k(D + z) ] P'~ = Pg71[ e--~shk D -1
(17)
T h e scales of Ph a n d p,~ are then: [ o (p~), o (p~)] = p9
+ D
,7
1
co2- )I,
and the ratio of the two scales is:
O(pn) H (1 0 (Ph) = 2D + H
1 ) coshkD
(19) '
As can be seen from the above equation, the ratio of t h e two pressure scales is a function of k D and H / D . Because k D = 2~rD/L, the scale of n o n h y d r o s t a t i c pressure is related to t h a t of
A Comparison
899
hydrostatic pressure by the aspect ratio and the H/D ratio. From equation (17), it can be seen that Ph has two parts. The first part is associated with ~, while the second part is proportional to the z-coordinate. Because the horizontal length scale for 7/is L, the first part of hydrostatic pressure has a horizontal length scale of L. The second part of hydrostatic pressure has an infinite horizontal length scale because it is not a function of x. Correspondingly, the hydrostatic pressure scale shown in equation (18) also has two parts. The first part is associated with the wave height, while the second part is related to the water depth. For a small amplitude wave with a wave height that is much smaller than the water depth (H << D), the second part is much larger than the second part. Therefore, we have
O(p~,)lH( O(ph)-2D
1
1 ) cosh27ra
(20) '
where c~(= D/L) is the aspect ratio. It is obvious that regardless the value of the aspect ratio, O(pn) << O(ph) is always true as long as the wave height (H) is much smaller than the water depth (D). For example, for a H/D ratio of 0.04, equation (18) suggests that O(pn) is always smaller than 2% of O(ph). For a small aspect ratio of c~ = 0.0455, O(pn) ~- O.O00790(ph) when H/D = 0.04, or the scale of nonhydrostatic pressure is more than three orders of magnitude smaller than that of hydrostatic pressure. Even for a relative large aspect ration of a = 0.5, we have O(pn) ~- O.O180(ph) when H/D = 0.04, or almost 1.74 orders of magnitude smaller than the scale of hydrostatic pressure. Now, let us look at the horizontal gradients of the two pressure components. From equation (17), we have:
(Oph Opt,) 0~7 [ coshk(D+z) Ox ' Ox = Pg"~x 1, cosh kD
] - 1
(21)
and
o\--~-z],O\ox/j=pg-~ or
0 (Op,~/Ox)0 (Oph/OX)
cosh kD cosh kD
1-
1,
1
c-~hshkZ~
1 cosh 2~ra"
'
(22)
(23)
The above equation indicates that the ratio of the two horizontal pressure gradients is a function of the aspect ratio only. Unlike equation (19), the right-hand side of equation (23) is not proportional to the H/D ratio, or the H/D ratio does not have any effect on the relative scale between the two horizontal gradients. Form equation (22), it can be seen that the horizontal gradient scale for hydrostatic pressure is only related to the wave height. The reason for the disappearance of the water depth from the expression of O ( - ~h ) is that the second part of hydrostatic pressure has a horizontal length scale that is infinite. Comparing equations (23) with equation (20), it can be seen that the right-hand side of equation (23) is about 2D/H times of that of equation (20). The ratio of the horizontal gradient of nonhydrostatic pressure over that of hydrostatic pressure approaches 1 as the water depth becomes larger (deep water). As the aspect ratio approaches zero, it approaches 0. For an aspect ratio of o = 0.0455, the right-hand side of equation (23) is about 0.04. Or, the scale of the horizontal gradient of nonhydrostatic pressure is about 1.40 orders of magnitude smaller than that of hydrostatic pressure. In this case, although nonhydrostatic pressure is more than three orders of magnitude smaller than hydrostatic pressure, their horizontal gradients are not differed by three orders of magnitude. Instead, their horizontal gradient scales are differed only by 1.40 orders of magnitude. On the other hand, for an aspect ratio of c~ -- 0.5, the right-hand side of equation (22) is about 0.91. In this case, the horizontal gradients for nonhydrostatic and hydrostatic pressure components are basically in the same order of magnitude, even when the scale of nonhydrostatic pressure is about 1.74 orders of magnitude smaller than that of hydrostatic pressure.
900
X. CHEN T a b l e 1. A c o m p a r i s o n of scales of h y d r o s t a t i c a n d n o n h y d r o s t a t i c p r e s s u r e c o m p o n e n t s as well as t h e i r h o r i z o n t a l g r a d i e n t s for H I D ---- 0.04.
(?(p,o o(pu)
o(op,~/o~) o(oph/ox)
0.0455
0.000079 < 10 - 4
0.04 ~ 10 -1"4°
0.5
0.018 = 10 -1' 74
0.91 = 10 -0.04
D L
Table 1 is a recap of the above discussion. A much smaller nonhydrostatic pressure component in comparison with the hydrostatic pressure component does not necessary guarantee its horizontal gradient to be much smaller than that of hydrostatic pressure. The scale difference is much smaller for the horizontal gradients than the scale difference between the two pressure components themselves. For a large aspect ratio of a = 0.5, the horizontal gradient of nonhydrostatic pressure is almost in the same order of magnitude as that of hydrostatic pressure without the need t h a t the scale of nonhydrostatic pressure to be in the same order of magnitude as that of hydrostatic pressure. This means that the condition that the nonhydrostatic pressure is negligible in comparison with hydrostatic pressure is not sufficient for the horizontal gradient of nonhydrostatic pressure in the horizontal momentum equation to be small enough (in comparison with the horizontal gradient of hydrostatic pressure) so that it can be taken away from the equation. For shallow water waves where a approaches zero, the right-hand side of equation (23) does approach zero, indicting that the nonhydrostatic effect in the horizontal momentum equation is indeed insignificant. In practice, if the D / L ratio is less than 0.05, the wave is treated as a shallow water wave [19] and the hydrostatic assumption is normally assumed to be valid. Nevertheless, as discussed above, for a = 0.0455, the magnitude of the horizontal gradient of nonhydrostatic pressure is only about 4% of that of hydrostatic pressure. Although the effect of nonhydrostatic pressure on the flow is indeed much smaller than that of hydrostatic pressure, such a nonhydrostatic effect cannot be neglected for an initial-value problem such as seiching in a closed basin. Model results with and without the consideration of the nonhydrostatic effect is quite different, even when the horizontal gradient of nonhydrostatic pressure is more than one order of magnitude smaller than that of hydrostatic pressure for seiching. The reason for this is that the nonhydrostatic effect is cumulative and cannot be cancelled out within the wave cycle. This can be seen from equation (21) and be explained here. Equation (21) shows t h a t both hydrostatic and nonhydrostatic pressure gradients in the horizontal direction are proportional to the gradient of surface elevation. However, their effects on the flow are opposite because they have opposite signs. The sum of the two horizontal gradients is:
Op _ Oph Op,~ cosh k(D + z) O~ Ox O~ + ~ = pg c o s h k D Ox"
(24)
If the nonhydrostatic effect on the flow is neglected, then
Op 077 cosh k(D + z) O~ O---x= Pg--~x -~ pgl cosh kD Ox' where g,
cosh kD -- gcosh k(D + z)"
(25) (26)
In other words, neglecting nonhydrostatic terms in the horizontal momentum equations is equivalent to an increase of the gravitational acceleration, which, for a = 0.0455, can be as much as 4%. From the dispersion relationship, it can be seen that an increase of the gravitational acceleration means an increase in wave frequency because the wave number k is fixed for the first mode standing wave in a closed basin. It can be proved that for a --- 0.0455, the increase
A Comparison
901
in the wave frequency using the hydrostatic assumption is about 1.5%. Although this is just a small increase, the solutions with and without nonhydrostatic effect will be quite different after a few cycles, as can be seen in Figures 6 and 7. For c~ = 0.5, the increase in wave frequency is about 77.6%. As discussed in the last section and shown in Figure 4, the wave period is reduced from 3.59 to 2.02. The solutions with and without the nonhydrostatic effect will be significantly different even in the first wave cycle. Errors induced due to the exclusion of nonhydrostatic pressure in one oscillation cycle will affect model results in the subsequent cycles.
5. C O N C L U S I O N S The hydrostatic pressure assumption has been widely used in many numerical models for free surface flows in shallow waters, because the assumption allows governing equations for the free-surface flows to be simplified. However, there are many cases where nonhydrostatic effects should be taken into account. A general rule of thumb for a safe use of the hydrostatic pressure assumption in a free-surface flow is that the length scale in the vertical direction should be much smaller that that in the horizontal direction, so that the nonhydrostatic pressure component is negligible in comparison with the hydrostatic pressure component. Nevertheless, there exist no quantitative guideline as to how small the aspect ratio should be for the nonhydrostatic effect to be negligible. In literature, it is generally believe that when the aspect ratio is less than 0.05, it is a shallow water flow and the nonhydrostatic effect can be neglected. For a better understanding of the relative magnitude of nonhydrostatic pressure in comparison with hydrostatic pressure, this study compares magnitudes of hydrostatic and nonhydrostatic pressure components and their horizontal gradients in two seiching cases. The first one is a deepwater wave oscillation with an aspect ratio of 0.5, while the second is a shallow water wave oscillation with an aspect ratio of 0.0455. Magnitudes and scales of hydrostatic and nonhydrostatic pressure components were compared and analyzed using model results and analytical solutions. Both model results and analytical solutions have shown that the nonhydrostatic component is much smaller than the hydrostatic component for small amplitude waves, even when the vertical length scale is comparable to the horizontal length scale. On the other hand, the horizontal gradient of the nonhydrostatic component can be in the same order of magnitude as that of hydrostatic pressure, even when nonhydrostatic pressure is much smaller than hydrostatic pressure. This kind of an "amplification" of the nonhydrostatic effect is due to the fact that hydrostatic pressure has a horizontal length scale that is infinite. The scale analysis shows that nonhydrostatic pressure is much smaller than hydrodynamic pressure because the ratio of the two scales is not only a function of the aspect ratio, but also proportional to the ratio of wave height to the water depth (H/D), which is normally small for small amplitude waves. However, once the horizontal derivative is taken, the ratio of the two horizontal gradient scales is no longer proportional to the D/H ratio. Comparison of model results of hydrostatic simulation with those of nonhydrostatic simulation reveals the differences between the two. Figures 2 and 3 and Figures 6 and 7 have shown that simulated velocities and surface elevations with and without the hydrostatic pressure assumption are totally different. Model results indicate that the hydrostatic assumption tends to produce higher mode oscillations and forces the first mode oscillation to have a phase speed of x/~, even for a deep-water wave.
In summary, this study has shown that the condition [Pnl << Ph does not guarantee that [°0-~I is smaller than [-~z~ ] by the same order of magnitude. Even when ~ is in the same order of magnitude as °o-~zh, the nonhydrostatic pressure component is still about two orders of magnitude smaller than the hydrostatic pressure component. As such, the condition that °0-~is comparable to ~ does not necessarily require ]Pnl ~ Ph. In other words, ]p~] ~ Ph is not the necessary condition for ~ and Oph to be comparable.
902
X. CHEN
REFERENCES 1. Y.P. Sheng, A three-dimensional mathematical model of coastal, estuarine and lake currents using boundary fitted grid, Report No. 585, A.R.A.P. Group of Titan Systems, Princeton, NJ, (1986). 2. A.F. Blumberg and C.L. Mellor A description of a three-dimensional coastal ocean circulation model, In Three-Dimensional Coastal Ocean Models, (Edited by N. S Heaps), pp. 1-16, American Geophysical Union, Washington, DC, (1987). 3. V. Casulli and R.T. Cheng, Semi-implicit finite difference methods for three-dimensional shallow water flow, Int. J. for Numer. Methods in Fluids 15, 629-648, (1992). 4. J.M. Hamrick, A three-dimensional environmental fluid dynamic computer code: Theoretical and computational aspects, Special Report in Applied Marine Science and Ocean Engineering No. 317, Virginia Institute of Marine Science, College of Willian 8z Mary, Gloucester Point, VA~ (1992). 5. X. Chen, A flee-surface correction method for simulating shallow water flows, Journal of Computational Physics 189, 557-578, (2003). 6. V. Casulli and G.S. Stelling, Numerical simulation of 3D quasi-hydrostatic free-surface flows, Journal of Hydraulic Engineering 124 (7), 678-686, (1998). 7. X. Chen, Three-dimensional, hydrostatic and non-hydrostatic modeling of seiching in a rectangular basin, In Estuarine and Coastal Modeling, Proe. of the 6th Intern. Conf., (Edited by M.L. Spaulding and H.L Butler), pp. 148-161, ASCE, New Orleans, LA, (1999). 8. X. Chen, A fully hydrodynamic model for three-dimensional, free-surface flows, Int. J. for Numer. Methods in Fluids 42, 929-952, (2003). 9. H. Lamb, Hydrodynamics, Sixth Edition, Dover Publications, New York, (1932). 10. J. Pedlosky, Geophysical Fluid Dynamics, Springer-Verlag, New York, (1979). 11. G.H. Golub and C. vanLoan, Matrix Computations, John Hopkins University Press, Baltimore, MD, (1990). 12. H.A. vanderVorst, Bi-CGSTAB: A fast and smoothly converging variant of bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13 (2), 631-644, (1992). 13. J. Smagorinsky, General circulation experiments with primitive equations, I. The basic experiment, Monthly Weather Review 91 (3), 99-164, (1963). 14. Y.P. Sheng and C. Villaret, Modeling the effect of suspended sediment stratification on bottom exchange process, J. Geophys. Res. 94, 14429-14444, (1989). 15. X. Chen, Effects of Hydrodynamics and Sediment Transport Processes on Nutrient Dynamics in Shallow Lakes and Estuaries, Ph.D. Dissertation, University of Florida, Gainesville, FL, (1994). 16. X. Chen, Responses of a hybrid z-level model to various topography treatments for a boundary-value problem and an initial-value problem, In Estuarine and Coastal Modeling, Proc. of the 7eh International Conference, (Edited by M.L. Spaulding), pp. 614-627, ASCE, St. Pete Beach, FL, (2001). 17. X. Chen, Fitting topography and shorelines in a 3-D, Cartesian-grid model for free-surface flows, In Computational Fluid and Solid Mechanics 2003, Volume 2, (Edited by K.J. Bathe), pp. 1892-1895, Elsevier Science Ltd., (2003). 18. X. Chen, Using a piecewise linear bottom to fit the bed variation in a laterally averaged, Z-coordinate hydrodynamic model, Int. J. for Numer. Methods in Fluids 44, 1185-1205, (2004). 19. R.G. Dean and R.A. Dalrymple, Water Wave Mechanics for Engineers and Scientists, Princeton-Hall, Inc., Englewood Cliffs, NJ, (1984).