A three-dimensional fixed grid model for shallow-water flow

I

I

I

I

I

I

I

I

J

I

I

I

I

I

rÀ'h~,

TU

Delft

I

Delft Universityof Technology

Department of Civil Engineering

Hydraulicand GeotechnicalEngineering Division HydromechanicsSection

I

I

I

I

I

I

I

I

I

I

-

I

I

I

I

I

I

A three-dimensional fixed grid model for

shallow-water flow

M.D.J.P. Bijvelds report no. 6-98

1998

The work reported herein has been made possible by financial support of RIKZ, Rijkswaterstaat

~f~

'

,I

,

~

t

~.

TU Delft

Ph.D. student, Hydromechanies Section, Faculty of Civil Engineering and Geosciences, Delft University of Technology, P.O. Box 5048, 2600 GA, the Netherlands. E-mail: M.Bijvelds@ct.tudelft.nl

I

~

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

,

J

I

I

I

I

I

I

I

I

Abstract

I

I

I

,

I

In this report the implementation and testing of a numerical model that is based on a Carte-sian fixed grid in vertical direction is described. The model uses the shallow-waterequations

and accounts for effects of stratification. In stratified environments, the terrain-following0

-transformation, which is often used in shallow-water models,is known to introduce numerical

errors in computing the baroclinic pressure gradient,especially in regions with varying depth. Although quite some research was carried out in the past on the improvement of the imple -mentation of the baroclinic pressure gradient, it still is not entirely clear to what extent the

e

rrors related to this problem harm computational results.

In order to compareboth approaches, i.e. the Cartesian fixed grid and the co-ordinate

system based on the o-transformation, in a consistent way,identical discretizations of time

and space were employed in both models, This allows to draw conclusions about the errors

that are related to the co-ordinate system used,without interference of other uncertainties.

It is shown that for relatively simple flow configurations,both models yield comparable

results. In case of stratified flow over a rapidly varying bathymetry, significant differences

between the models wereobserved.

I

I

I

I

I

I

I

I

I

I

I

I

Contents

I

I

1 Introduction

2

2

Definition of the model

6

2.1

Governingequations 6

2.2

Numerical model . . . 7

2.2.1

Grid definition 7

2.2.2

Solution strategy 9

2

.

2.3

Discretization of pde

'

s

12

3

Testcases 15

3.1

Steady channel flow

15

3.2

Exchange flow. . . .

18

3.3 Standing wave in basin.

20

3.4

Flow over a trench

..

21

4

Preliminary conclusions

26

References

28

List of figures

30

List of symbols

32

A Spatial discretization of pde's

34

I

I

I

,

I

I

I

I

I

I

I

1

I

I

I

I

I

I

Chapter

1

Introduction

I

I

I

The application of large three-dimensional numerical models for free surface flows in civil

engineering is common practice nowadays. These models often cover areas such as coastal seas and estuaries, where the bathymetry may vary rapidly. Moreover, due to tidal influences

the position of the free water surface may vary significantly as a function of time resulting in

drying and flooding of tidal flats. Hence

,

the local water depth depends both on the horizontal

position and time.

In these mode Is the o-transformation, generally ascribed to Phillips (Phillips, 1957), is often employed. By introducing this dimensionless vertical co-ordinate, the number of grid layers becomes fixed and the relative layer thickness is independent of time and space. A transformation that is widely used is given by

I

I

I

I

z

-(

a=d+( (1.1)

I

I

where z is the Cartesian vertical co-ordinate, ((x, y,

t)

denotes the position of the free surface

as a function of horizontal co-ordinates x and y, and time

t,

and d(x, y) is the depth with respect to a reference level, see Figure 2.1. Hence, the total water depth is given by H

=

d+(.

This terrain-following vertical co-ordinate system has several advantages compared to the Cartesian co-ordinate system 1 (for the lat ter see e.g. (Leendertse & Liu, 1975), (Casulli &

Cheng, 1992)). Since the computational domain is fixed when the SCM is applied, treat-ment of kinematic boundary conditions is very simple in this case. In the FGM the free surface boundary is not fixed which makes it elaborate to treat. Furthermore, the stepwise representation of the bed topography can distort the bed shear stress and thereby effect the flow (Sheng, 1983). The staircase representation of the bed and the free surface, also leads to a 10cal degradation of the truncation error of the spatial discretization, since only first-order

I

I

I

1In the remainder of this report the zr-co-ordinate model and the Cartesian fixed grid model will be referred to as"SCM" and "FGM", respectively.

I

I

I

I

I

approximations can be used at these artificial boundaries. Besides,it makes implementation of the discretized equations substantially more cumbersome. Finally, in the FGM the shallow

areas may be under-resolved unIess a very fine discretization is used.

The o-co-ordinate transformation has some inherent disadvantages as weIl. Due to this transformation, the number of nodes in the vertical direction must be constant over the entire computational domain. This lack of flexibility leads to unnecessarily high grid resolution in shallow areas and unsufficient grid resolution in deep parts of the computational domain. In very shallow drying areas, the mapping may even become singular. Further numerical problems may arise near steep bed slopes combined with vertical stratification of the density. For stratified flows over continental slopes, the currents induced by discretization errors may be as large as the typical physical veloeities (Walters & Foreman, 1992).

The horizont al derivatives in o-co-ordinates can be expressed as

I

I

I

I

I

(1.2)

I

I

where x" is the horizontal co-ordinate along a plane where O"=constant. As can be seen in eq. (1.2) in terrain-following co-ordinates the pressure gradient is the sum of two terms, which

may both be relatively large near steep topography and of opposite signs. Small truncation

errors in the approximation of both terms may result in a large error in the total pressure gradient, causing artificial flow and hence unrealistic mixing in case of stratification. If we consider the discretization of eq. (1.2), employing central differences and assume a rigid lid approximation (i.e.

H ~ d)

8if> ,....,if>i+l,k - if>i-l,k 0" 8d if>i,k+l - if>i,k-l

8x ,...., 2.6.x -

d

8x* 2.6.0" (1.3)

I

I

where if>is concentration or velocity,for example,the truncation error becomes

(1.4)

I

I

This error can be reduced by subtracting a reference state density field before computing the pressure gradient terms (Gary, 1973). However, this may be oflittle help in large areas where density can deviate significantly from the reference state. Related to this problem, several authors (see e.g. (Janjié, 1977), (Mesinger, 1982)) have identified a problem of "hydrostatic consistency". In (Haney, 1991) this consistency relation is given by

I

I

I

1 0"

8HI

- - .6.x

<

.6.0" H 80" (1.5) 3

I

I

I

I

I

I

This condition guarantees that the a-surface immediately below (above) a given o-surface remains below (above) the given a-surface within a horizontal distance of one grid interval. If this condition is not satisfied, the numerical scheme may be non-convergent. The first term on the right hand side of eq. (1.4) shows that unless the hydrostatic consistency is fulfilled, the error caused by the b..x term on the right-hand-side is dominant and increasing the vertical resolution will have little effect on the error. In drying regions it is hardly possible to fulfill condition (1.5), which means that convergence might become impossible.

Another problem related to the a-transformation is due to the non-conformality of the mapping. The transformation is not angle preserving which leads to additional terms that are proportional to

8H/8x

and

8H/8y.

For the diffusion terms this leads to cross-derivatives terms which are often neglected in order to simplify implementation and to anticipate pro-blems concerning accuracy, stability and monotonicity. In the presence of steep beds this may lead to significant errors as shown in (Paul, 1994) and (podber & Bedford, 1994).

Stelling and Van Kester (Stelling & Van Kester, 1994) proposed a method for compu-ting horizont al diffusive fluxes. They approximate the horizont al diffusive fluxes based upon a finite volume approach. This allows treatment of these fluxes in Cartesian co-ordinates although a a-grid is used. Fluxes through cell edges are determined by using a non-linear filter. H, instead of a filter some interpolation formula is applied wiggles may occur and a small persistent artificial vertical diffusion will be present (except for linear vertical density distribution). The non-linear filter guarantees at least a first-order approximation of the ho-rizontal diffusive fluxes, depending on total water depth variation and density variation. It is also shown in this artiele that the min-max condition (no spurious oscillations in the solution) is fulfilled which is important to get physically realist ie solutions. It turns out that this con-dit ion is also sufficient for stability in case of explicit time integration. The same procedure is also used for calculating the horizontal pressure gradients. Three testcases show that arti-ficial flow due to a-transformation is minimized in case of stratified situations in combination with steep bed slopes. This finite volume approach does not have to satisfy the hydrostatic consistency condition, which is a positive property of this method, especially when estuaries are considered.

In this report, we focus on the development of a numerical model that uses Cartesian co-ordinates in the vertical direction. Although quite some attention has been paid to the problems related to the a-transformation, very little consistent comparison of both approaches is available from the literature. Quite often, it is unclear whether differences between results from numerical simulations and observations can be imputed to the schematization of the geometry, the turbulence model that is used or the co-ordinate system, for example. The

FGM is developed to compare results of numeri cal experiments to those obtained with the SCM. Similar discretizations for time and space are used in both models, in order to prevent

differences stemming from the method of integration of the partial differential equations (pde's). Such a comparison allows us to draw more definite conclusions about the performance of the two co-ordinate systems, without interference of other uncertainties.

I

I

I

I

I

I

I

I

I

I

I

I

4

I

I

I

I

To test the implementation of the FGM first flow configurations in which the SCM reduces

to the FGM are considered, i.e. in which gradients of the water depth are small. Both unstratified (wide channel flow) and stratified (exchange flow) are considered. The treatment of the moving free surface in the FGM is tested by simulation of a standing wave.

After testing the implementation of the FGM we focus on situations in which deviations between the two approaches are expected to occur. Flow in a wide channel with a sloping bed gives rise to a stair-case boundary at the bed in the FGM, and is known to introduce errors. Furthermore, we focus on a schematic representation of stratified flow in estuaries. As a first step towards practical computations, a stratified tidal flow over a rapidly changing bathymetry is considered, a situation in which errors due to the a transformation are most pronounced.

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

5

I

I

I

I

I

I

Chapter

2

,

I

_{Definition of the model}

I

I

I

I

2

.

1

Governing equations

I

For the

applications that we consider

,

i.e.

flowsin coastal areas,

a characteristic length

scale

of the flow in the vertieal direetion is mueh smaller than a eharaeteristie length seale in the

horizontal directions. In addition, the characteristic turbulent stresses of the fluid are small

compared to the characteristic hydraulic pressure gradients (Vreugdenhil, 1994). Under these circumstances, the full three dimensional Navier-Stokes equations may be reduced to the

boundary layer form and the equation for the vertieal veloeity w reduees to the hydrostatic

pressure assumption:

I

8

p

-

=-pg

8

z

(2.1)

I

I

where p is the pressure and 9 is the gravitational acceleration and p is the density of the

Huid. Furthermore weuse the Boussinesq approximation which assumes that the influence of density variations is only envisaged in the buoyancy term. The momenturn equations then read

(p

o

is a referenee density):

Du

=

_

g

8( _!!_

{ç

8p d

z

'

+ _I_

(8T

XX

+

8T

xy

+

8T

XZ) _

J

v

(2.2a)

Dt

8

x

po}

z

8

x

Po

8x

8y

àz

I

I

I

I

D

v

=

_g

8(

.:s.

{ç

8p d

z

'

+

_I_

(8T

yx

+

8

T

yy

+

8T

yZ)

+

J

u

Dt

8y

p

o

}

z

8y

Po

8

x

8y

8

z

(2.2b) with D. 8. 8. 8. 8.

-=

-+u-+

v

-

+w-Dt

8t

8

x

8

y

8

z

(2.3) 6

I

I

I

I

I

Here, u, v and ware veloeities in x- y- and z-directions, T is the Reynolds stress tensor and

f

is the Coriolis parameter. The non-conservative formulation of eqs. (2.2a-2.2b) results from the assumption that the Huid is incompressible. The continuity equation may then be written as

I

ou

+

o

v

+

ow

= 0

ox

ay

a

z

The above equations (2.1-2.17) are generally referred to as the shallow-water equations. The Reynolds stress tensor T in the momentum equations is modelled by introducing the

Boussinesq hypothesis: (2.4)

I

I

I

I

I

Ti

j

( OUi OUj)

- ==V«

-+-p

aXj

aXi

Here tensor notation is used for convenience sake. A large amount of models are available for computing the eddy viscosity

Vt

.

Among the most popular turbulence models in civil engineering are the simple algebraic turbulence model and the more advanced k-é turbulence model. For more detailed information on turbulence modelling the reader is referred to the

literature (e.g. (Rodi, 1980), (Wilcox,

1993))

(2.5)

Transport of scalars is represented by

De _ ~

(Vt

oe)

+ ~

(Vt

oe)

+

i_

(Vt

oe)

+

S

Dt -

OX

at

OX

OX

at

ay

oz

at

az

(2.6)

I

where e is a scalar quantity,

at

is the PrandtljSchmidt number and

S

is a souree term.

I

2.2

2.2.1

Numerical model

Grid definition

I

In contrast with the a-co-ordinate model, the number of active grid layers in the FGMmay be

a function of time and space, see Figure 2.1. Given a number of horizontallevels that do not change with time and horizont al co-ordinates, with vertical co-ordinates between Zkmin and Zkmax, the active vertical grid is computed every half time step. The vertical grid spacing in a (-point is computed by:

I

I

(2.7)

.

1

.

I

I

In this report

~Xi

(~z in eq. (2.7)) represents a grid spacing in the xi-direction, the indicesi,j, k are discrete representations of the X-, y- and z-co-ordinates respectively and the superscript n is a counter for time administration. A computational cell is set active whenever ~Z~j,k

>

O.

In contrast with models that are based on fixed grids in the vertical direction reported in the 7

I

I

I

I

Zkmax ,.---,---:----,

I

_{z,a, k}

I

I

1

k+1

ç

... . ... ...; ~ : . •••••• ,•• ••••••• ••1 • - + .; k ···;···1···:··· k-1 d k=kmax Zkmin+l .

I

I

I

I

a=-ll.---Zkmin '--- __ ~ __ ~ ~_~_~_~---'

I

I

I

I

I

I

I

I

---X,Y,t,}

Figure 2.1: Grid definition in case of fgm model and scm model, respectively.

literature, the free water surface is aUowed to cross ceU faces. No restrietion is imposed for

the minimal vertical grid size. For small grid spacings, the coefficientmatrix of the system of equations that needs to be solved may become badly conditioned. These matrices are solved with the double-sweep algorithm. To prevent the double sweep algorithm from producing too large rounding errors, it might be necessary to introduce a threshold (an equivalent of the threshold for drying and flooding in the SCM) under which a computational cell is set dry. However, in the present study, numerical experiments did not suffer from this type of behaviour. To control the drying and flooding process properly, the introduetion of three -dimensional mask arrays proved to be imperative.

To ensure positive water levels, the total water depth in a velocity point (u or v) is

determined by an upwind procedure (Stelling & van Kester, 1996):

.6.Z~j,k =min

((

+

) -

max (-d+,Zk-l) (2.8) where ifUi+l/2,j

2:

0 ifUi+l/2,j

<

O. (2.9) and (2.10) 8

I

I

I

I

i,j+1/2,k

I

I

i-1/2j,k ij,k

+

- i+1/2j,k

I

:v - :u I

I

I ij-l/2,k

I

I

I

I

Figure 2.2: Definition of the staggered grid in the horizontal plane

The depth in a v-point is computed by replacing the subscript ibyj and u byv in the above

formula. A staggered grid (C-grid of Arakawa) for the position of the variables is used, see

Figure 2.2.

Sincethe grid spacing near the bed and free surface may vary as a function of iand j,

velocity points on the staggered grid of two adjacent grid cells may be situated at different

vertical

l

evels. However

,

t

h

e free s

ur

face leve

l

genera

ll

y va

r

ies slowly in time and space and

bed topography is smooth in most areas. Therefore, in the present approach, no correction for this effect is taken into account.

I

2.2.2

Solu

t

ion strateg

y

A coupled set of equations for the veloeities (u and v) and the water surface elevation ( is obtained by manipulation of the continuity equation. We start with the kinematic boundary condition at the free surface, which requires that a partiele must follow the motion of the free

surface,

I

I

I

I

o

;

a(

o;

w

=

-

+

u-

+

v- at z

= (

at

ax

ay

(2.11)

Weintegrate eq. (2.4) with respect to the water depth,

j

(

(

au

_-+-+-

av

aw

)

d

_z-

- 0

-d

a

x

ay

a

z

(2.12)

and apply Leibnitz' rule to obtain

I

I

I

a

(j

(

)

a

(

a

(j

(

)

o

;

-

_a

_x

ud

z

-uk-+-

vd

z

-vk-+wk-wl

-

d=O

-d

ax

ay

-d

ay

(2.13) 9

I

I

I

I

Inserting the free surface kinematic boundary condition into eq. (2.13) gives

I

~~+

:x (/_~ u

dZ)

+

:y (/_~

v

dZ)

=0 (2.14)

I

Eq. (2.14) and (2.2a-2.2b) form the set of equations that are solved numerically and yield theunknowns (, u and v. Before discussing the numerical details of the model, the schematic procedure of the solution algorithm is presented in table 1.

The shallow-water equations are solved by a two-step procedure in time (ADI) using finite differences. In the first half time step n

+

1/2, the water surface elevation ( and the velocity in the x-direction, u, are computed implicitly. The velocity in the y-direction, v, is taken from the previous time step n and thus integrated explicitly. In the vertical direction, both adveetion and diffusion are computed implicitly in order to avoid excessive small time steps imposed by the relatively small grid spacing in this direction. All other terms are computed explicitly. In the second half time step n

+

1 where v and ( are coupled, the same procedure as in the first half time step is followed.

Transport of matter is computed using a finite volume approach in order to ensure conser-vation of mass. Horizontal derivatives are treated explicitly and vertical transport of matter

is treated implicitly, again to avoid small time steps.

I

I

I

I

I

I

I

I

I

I

I

10

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

initialize variables

for all time steps (n=l, nmax) do

stage 1

- compute water depth and in u-points using upwind approach, compute index upper grid layer in u-points and (-points _ compute (n+l/2, un+1

- compute new water depth and in u-points using upwind approach, compute new index upper grid layer in u-points and (-points - check for dryingjflooding in u-points

if necessary (dryingj flooding is true) repeat stage 1

k n+l/2 /

c

n

- adjust un+1 to fulfill ""_L,.,kmmmax (un+1ó.zn+1 2) =

f

_-d un+1dz

_ compute cj;n+l/2

- compute bed shear stress in v-points

stage 2

- compute water depth and in v-points using upwind approach, compute index upper grid layer in v-points and (-points _compute (n+l, vn+1

-compute new water depth and inv-points using upwind approach, compute new index upper grid layer in v-points and (-points if necessary (dryingj flooding is true) repeat stage 2

_adjust vn+1 to fulfill""kmaxn+1 (vn+1ó.zn+1) = fCn+1/2 vn+1dz

L,.,kmm -d

_ compute cj;n+l

- compute bed shear stress in u-points _ compute pn+1 (cj;)

_compute kn+l

- compute [n+l

od

Table 1: schematic outline of the numerical procedure 11

I

I

I

I

2.2.3

Discretizati

o

n

o

f pde's

Shallow-water equations

Asmentioned in the previous section,in the first half timestep the surface elevation( and the

velocity u are computed implicitly. In order to do so,we start from the momentum equation in the x-direction per grid layer k,

I

I

I

I

I

I

I

I

un+_i+l1_/2,j,k - un_Hl/2,j,k 8rn+1_..,

₊

_n 8un+_H11/2_,J._,k

~

t

+

g~ wi

+

1/2,j,k 8z _ ~ ( n 8U~:11/2,j,k) _ RHSn 8z l/i+l/2,j,k

B

z

-

Hl/2,j,k (2.15)

The right-hand-side (RHS) ofeq. (2.15) contains horizontal advectiveterms, the baroclinic pressure gradient, horizontal viscosity and the Coriolis force. The horizontal adveetion of

momentum is approximated by a first-order streamwiseupwind method. Verticaladveetion and diffusion isdiscretized by meansof first-order upwindand second-ordercentral differences,

respectively. Detailson spatial derivativescan befound in appendix A.

B

y

int

e

grating eq

.

(2.15) ov

e

r the water depth

H

~-H

/

2

,j'

a linearized expres

s

ion for th

e

depth averaged momentum flux through the cell faceat (i

+

1/2, j) is obtained:

1

(i+

1/2,j

- Twind

+

Tbed

+

_o

RHS

~I'2

_{I ,}J k_,dz

(2.16)

I

Thematrix-coefficientsfor the system ofequations thusobtained arecoupled inthe horizontal direction via the barotropic pressure gradient, and in the vertical direct ion via the vertical

adveetion and diffusion terms. In order to solvethis five-points stencilefficiently,the system isdecoupled in the vertical direction by applying a double-sweep procedure.

Using central differences for the horizontal derivatives, the discrete continuity equation,

is given by

I

I

I

C+_i,j1/2 _- (n_i,j

₊

(:::-:-11_U +1Hn) _{Hl/2,j -} (:::-:-11+1_U Hn) _i-l/2,j

~

~t

~x

+ (

vnHnkj+l/2 - (vnHnkj_l/2 =0

b.

y

(2.17)

I

I

I

By substituting eq. (2.16) into eq. (2.17) weget atridiagonal system in(n+l/2, that is solved by applying double-sweep recursion. Back-substitution of the new water level(n+l/2 into

I

I

I

I

the momentum equation per layer yields the velocity u~+l, which in two-dimensional cases is

equal to the depth averaged velocity

un

+

1.

Accurate integration of the explicit part of the time derivative in eq. (2.16), isofcourse

important in order to be consistent with the numerical method proposed and hencewith the SCM. A natural candidate to use for this explicit part of the time derivative would be the

dischargeper unit width of a cell from the previous time step, QH-l/2,j,k' that directly follows from substitution of the new water levels in eq (2.15). Division by the vertical grid spacing

and the time step is necessary to obtain an acceleration term. Sincethe vertical grid spacing

may approach zero,this gives riseto instabilities in the numericalexperiments, and therefore the velocity itself divided by the time step is used. However,the maximum number ofgrid layers mayvary in time at a certain position

(

i,

j)

and hence correct integration of

J

~~

u

n+1

d

z

isnot ensured. To fulfill this constraint, it isnecessary to adjust the velocity profilefor the

horizontal velocities.

In the present approach, fluxes at cell faces of computational cells that do not have an

activeneighbouring cell,maybeunequalto zero. Forcingthe veloeitiesat the cell facestobe equal to zero,would lead to additional damping of the system, not related to the physicsbut

to the treatment of the free surfaceboundaries. In order to fulfill the continuity equation, the

m

ass

fiux

es

through the cell face

s

need to be added to the neighbouring column of

c

ells. This

means that in the determination of the coefficientsrelated to the explicit part ofeq. (2.4), a

larger number ofgrid layersis needed than for the implicit part of the equation.

Dueto the complicated tracking of the maximum number of active cells at a given position

(i,

j)

and the problems related to this tracking, it isconvenient to introduce arrays in which information of the activegrid of the previous time step isstored.

In the second half time step, the same procedure isrepeated in the y-direction in order to compute en+! and v~+l.

It is noted that in contrast with the numerical procedure followed by Delft3D-FLOW (Delft Hydraulics, 1998) and TRIWAQ(Zijlema,M., 1997),noequation forvn+1/2 and un+1/2 is solved in the first and secondstages, respectively. In combination with the explicit tre at-ment of the advective terms in the momentum equations, this yields a considerablespeed-up

of the computations. From experiencewith this numerical procedure it was found that this approach may reduce the amount of work by a factor 4 compared to the "standard" pro -cedure. This semi-implicit treatment allows for an practically acceptable magnitude of the

time step. The restrietion of this time step is imposed by the CFL stability criterion for the advective terms in horizontal directions only (if weneglect horizontal viscosity). In civil

engineering applications, horizontal grid spacings are usually large compared to the water depth, leading to practically acceptable maximum time steps, the moresoif time dependent

processes in irregular geometries are considered. In this casetime steps must be sufficiently small to model physical processes accurately.

I

I

I

I

I

I

I

I

I

I

.

1

-I

I

13

I

I

I

I

Transport

equation

The transport equation, eq. (2.6), is solved by a finite volume method on a structured grid.

This means that the value ofq;n+l is computed by considering the fluxes across the six ce ll-faces:

I

I

q;n+l q;n

i,j,k _f::l-_t i,j,k =F,_t._-l /2'_,),k

+

p._1,)-,1/2k

+

P._1,.₎k_{, -}1/2 -

P

_t+1/2_,),.k - p._t,)+,1/2 k -

P

_1,

.

₎k_,+1/2

(2.18)

I

Here, sourees or sinks are left out of consideration. The fluxes P, consist of transport due to adveetion and diffusion. Analogous to the procedure with the momentum equation, tran s-port in horizontal directions (x- or y-direction) is treated explicitly and transport in vertical direction (z-direction) is integrated implicitly. As a consequence

I

I

I

{

p

n+1

P=

k pn k if k=K+l/2, K=O..kmax-l if k=K+l, K=O..kmax-l (2.19)

I

Implementation ofa mass conserving solver for the sealar transport equation is more eumber -some than in the SCM. Due to the variation of the index of the top eell,kmax, and since cells

ofneighbouring columns (i+l ori-I in the x-direction forexample) may be situated above or under the eell under consideration, it is most eonvenient to start integration at the bed. At this closed boundary the flux of matter normal to the boundary is set equal to zero. If the

index of the upper grid layer of a position (i,j) has not changed during a half time step, the

new concentration in this top eell is eomputed by summation of the fluxes through the ce ll-faees, including the flux through the water surface. The latter, Fi,j,kmax+1/2 sterns from the summation ofthe mass fluxes of neighbouring cells with a higher top grid layer index kmax,

whieh are treated explicitly. In case the index of the top grid layer has increased compared

to the previous half time step, the fluxes of neighbouring cells are redistributed equally over

the eell ranging from kmax" to kmax'"!". If the free surfaee of neighbouring eells is below the free surface of the cell under consideration, this flux is set equal to zero. This way a mass eonservative implementation of the transport equation is obtained.

I

I

I

I

I

I

I

I

14

I

I

I

I

I

I

Chapter

3

I

_Testcases

I

I

I

I

I

In this chapter the performance of the FGM is tested against the SCM. The flows that we consider will generally be controlled by bed friction, barotropic pressure gradients and baro-clinic pressure gradients. In the test cases we focus on these physical phenomena. We start with two situations in which water depth variations are much smaller than the local water

depth. In these situations both models should yield (almost) identical results. After that, the

treatment of the free surface boundary and the bed friction in case of varying water depth is studied. In order to show some differences between both models, we conclude this chapter with computations that are representative of the stratified flow in areas with rapidly varying bathymetry, like estuaries.

3.

1

S

tea

dy channel flow

I

I

I

In case of steady and frictionless two-dimensional channel flow with horizontal bed, the SCM

reduces to the FGM. As bed friction is absent, there will be no gradient in surface elevation and hence the grid cells in the SCM coincide with a Cartesian grid. In this situation results of both models should converge since the time and spatial discretization are dealt with in the same manner in both models. However, this trivial test case is not adequate for testing the model since important features of open channel flow, such as bed friction and diffusion of momentum, are not included. Therefore,a channel with a Chezy coefficientequal to 55

m

1/2

js

was used for testing. The theoretical slope of the free surface can be obtained by reduction of the shallow-water equations. Assuming a steady uniform flow, which is a good approximation for small water surface gradients, the vertical viscosity term balances the barotropic pressure gradient:

I

I

I

I

!_

(v

t

ou)

=

gHo(

OZ OZ

ox

(3.1)

15

I

I

I

I

((x) u(z) 0.04 0 0.03 -0.2

.§.

_0.02 ...!... -0.4 ~ _-0.6 "Jo

--0.01 N _-0.8 0 -1 2500 5000 0 5 10 15 20 25 x

[mJ

ulu.

[

-

J

k(z) & ê(Z) l/t(z) 0 0 ~ -0.2 -0.2 ...!... -0.4 k ...!... -0.4 ~ _-0.6 ~ _-0.6

--

N _{-0.8 E}

--

N _-0.8 -1 -1 0 1 2 3 0 0.025 0.05 0.075 0.1

klu;, êHlu~

[-J

1/t!(Hu.)

[-J

I

I

I

I

I

I

I

Fi

g

ur

e

3.1

:

Water elevation as a function of space and non

-

dimensional velocity

,

turbulent

kinetic energy (TKE), dissipation of TKE and eddy viscosity as a function of depth for the

FGM (solid lines) and SCM (dashed lines).

After integration over the water depth, and using u.lu =.;gIG and 'Tbed =pu;, we get:

8(

u

i ui

=

8

x

-

HG2 (3.2)

I

This expression is known as Chezy's formula. For a 10 m deep channel with a mean flow velocity,

u,

of 0.5

mis

and a Chezy coefficient of 55 m1/2

I

s,

this yields a slope of about

8.3 x

10

-

6. For a 5000 m. long channel this leads to a water level drop of .04 m over the

channel. Since the slope of the water surface is small, differencesbetween the SCM model and

the FGM ought to be negligible, seeeq. (1.4).

The standard k-e turbulence model was used for the computations of the channel flow. A non-equidistant vertical grid, with decreasing mesh size towards the bed, was used for the

numerical experiment. The number of grid layers in the model was equal to 25. In Figure

3.1, the free surface elevation is shown as a function of the streamwise co-ordinate x. The

velocity profile, the profiles of k and

e

and eddy viscosity are plotted in the same figure as

function of depth z, No difference between the models can be observed, and the computed free surface slope agrees well with the theoretical one, given by eq. (3.2).

In case of a sloping channel bed, it is expected that differences between the SCM and

FGM do result. For this type of problems, the SCM is more suitable than the FGM. Due to

I

I

I

I

I

I

16

the stair-case representation of the bed in the FGM the flow locally behaves like flow over a backward facing step since in the numerical model the advective terms

u8ul8x

are switched off near a closed boundary. Furthermore, the grid spacings in the vertical direction of the two models do not coincide anymore as in the case with a horizontal bed.

In our experiment, we use a channel with a sloping bed of 1x10-3. The depth ranges

linearly from 5 m near the entrance and 10 m at the outflow boundary, the length L of the channel is 5000 m. The standard k-e turbulence model is used for the computations. Since two transport equations are solved for the determination of the eddy viscosity,these solution of these equations also suffer from the stair-case boundary at the bed. Deviations between the two models are likely to be most pronounced in this case. We use 10 grid layers of equal (relative) thickness. This means that a step is present every 10 grid cells in the FGM. The mean flow velocity is equal to 1

mis

at the inflowboundary. In Figure 3.2, veloeities and bed shear stress near x] L=0.8, just behind a step, are shown as a function of depth. It is observed that veloeities of the FGM agree weUwith those of the SCM model. Near the bed, differences

can be explained from the exclusion of the advective term

u8ul8x

near the step. Itis noted that veloeities are plotted at the (-points which means that veloeities are interpolated. The

I

I

I

I

1

I

I

I

I

I

I

I

u(z)

0

...---.

--...---.,

~i _)

-10 L.__--'-- __ __L___ __j

o

I

0.25 0.5

u [mis]

u(z)

0.75

I

_~i

o .----..,----...,..---,

_'

_{_)}

_/

-10 '-_--' __ ____L __ __.J

o

I

I

I

I

I

I

0.25 0.5

u [mis]

0.75

0.01

,....-

--..-

-.---,

0.0075 0.005 0.0025

o

'---~---~

2500 x

[m]

Tbed(x) 5000 0.01.---r---, 0.0075 0.005 0.0025

O

'--

-~

---__j

2500 x

[m]

5000

Figure 3.2: Velocity as a function of depth z at x=0.8L and bed shear stress as a function of

x. Solid lines represent the FGM and dashed lines represent the SCM. In the lower two figures the bed shear stress was determined in the second grid ceUabove the bed in case ofFGM.

bed shear stress fluctuates highly due to the strong variation of the vertical grid spacing and

I

I

I

I

velocity in this area. Although the velocity profile seems not to be affected too severely, it is an undesirable side-effect when morphological computations are executed. In these models

the sediment transport is strongly related to the bed shear stress. As a first improvement of this deficiency, the peaks in the bed shear stress profile can be suppressed relatively easy

by computing the bed shear stress in the second grid cell above the bed. This increases the magnitude of velocity somewhat near the bed since the bed shear stress has decreased,

and consequently decreases the velocity near the free surface. More sophisticated manners of determining the bed shear should be applied to improve this result.

I

I

I

I

I

I

I

3.2

Exchange flow

We consider a closed basin with length L=1000 mand depth H =10 m which is divided into

two areas with different concentrations of salt:

c =

{~

5

kg/rn:' if0::;x

< ~

if ~ ::; x ::;L (3.3)

I

I

I

Due to the baroclinic pressure gradient, gravity currents develop at the bed and the free surface; the lighter water tends to intrude under the more dense fluid resulting in two fronts

moving in opposite directions. Besides, an external wave,with very small amplitude compared to the water depth, is generated. The speed of the front of the gravity current can be deduced

fromconsidering the energy budget of the system. At t=O the total amount of potential energy is equal to ~gH(PI

+

P2). Assuming a frictionless flow, with no mixing between the layers,

the system must consist of two layers with different density and of equal thickness (the layer

of lighter water (pI), ofcourse, overlays the layer of heavier water (P2)). The potential energy of the transient state equals ~gH PI

+

tgH P2. Because of conservation of energy, the net change in potential energy, hH!::lp must be equal to the total kinetic energy, PQu;ront' hence

I

I

I

I

I

%0"=

~t~p

It is shown in (Kranenburg, 1993) that for frictionless exchange flow the celerity of the front issomewhat higher than given by eq. (3.4). The concentration distribution eq. (3.7), yields a

density difference of 19.4 kg/rn", leading to a theoretical celerity of the front of 0.69

mis.

The external wave generated by the intruding current has an amplitude which is much smaller, a

factor b..p

I

p, than the amplitude of the "internal wave", and therefore both models should lead to very similar results again. In Figure 3.3,the velocity profiles of u(z) and the vertical density distribution p(z) are plotted at a monitoring station at

t

L

for four different times. The front reaches this position at the same time for both models and the "final" distribution of velocity and density are almost identical at t=540 s. However differences can be observed

(3.4)

I

I

I

I

12.5 25 - _~ - ....!.... - ~

---- N 1 0 0 -0.2 ....!.... -0.4

::t1

_-0.6 c -<; N _-0.8 -1 -1 0 0 -0.2 ....!.... -0.4 >-~ ~ _-0.6

---

N _-0.8 -1 -1 -0.5 0 0.5 1 u [mfs];c [kgfm3] I u

I

I

-0.5 0 0.5 u [mfs];c [kgfm3] 12.5 25

I

I

I

I

o

o

-0.2 -0.4 -0.6 -0.8 -1 -1 12.5 25

,)

-,

_, 1\ " -0.5 0 0.5 1 u [mfs];c [kgfm3] 12.5 25

o

O

r--

-

-

-

~

-

~~~

-0.2 -0.4 -0.6 -0.8 -1 ~~L- __ ~ __ ~ __~ -1 -0.5 0 0.5 1 u [mfs];c

[kg/rn"]

Figure 3

.

3

:

Velocity

(u)

and concentration

as

a function of relative depth at four different times

at x=~. The Figures correspond with t=450 s, t=480 s,t=510 s and t=540 s,respectively.

Adveetion of momentum is discretized in non-conservative form. The solid lines represent the

FGM and the dashed lines represent the SCM.

I

I

I

at t=480 s and t=510 s which may be explained by considering the implementation of the

barotropic pressure gradient in the SCM. Using eq.

(1.2),

we may write the pressure gradient as

1 ap 9

en

1

0 I 9 a

1

0 I pg (a( aH)

--=-- pda+-H- pda+-

-+a-Po ax Po ax a Po ax a Po ax ax

Leaving out the water surface gradient, which is integrated implicitly,and rearranging terms weget

1apRHS 9

en

(1

0dl) gH a

1

0 d I

---=-- pa +pa

+--

pa

Po ax Po ax o Po ax a

This expression is used in the SCM code. In case of zero depth gradient, this expression is

equal to the baroclinic pressure gradient in the FGM. However,near the front the largest water surface gradients occur,especially in case a hydrostatic code is used. This may explain the differences observed in Figure 3.3.

The computed speed of the front is approximately 0.52 mfs, significantly slower than the

theoretical front speed. A significant part of this discrepancy can be attributed to the

disere-I

I

I

I

I

I

(3.5) (3.6)

19

I

I

I

I

I

I

tization of the advective terms in the moment urn equation. Using a conservative formulation for the advection, i.e.

~ou2lox

instead of

uoulox,

the front speed increases to about 0.62 mis, see Figure 3A. It is noted that ot her approaches for treatment of adveetion are also possible, allleading to different adveetion approximations in the vicinity of discontinuities.

12.5

o

o

.-

---~----

--

~

-0.2 f--OA I-c -0.6 I--0.8 I--1 '---'---'-'---'---__j -1 u

I

I

-0.5 0 0.5 u

[m/s]

;«

[kg/m3] 12.5 25

I

I

I

o

o

.-

---

~

--~--~

-0.2 -0.4 -0.6

-0

.

8

-1to..____ -'-- __ -'- __ ---'- __ ~ -1 -0.5 0 0.5 1 u [m/s];c [kg/m3] 25

o

0

..--

-

-

.,...,---,

-0.2 -OA -0.6 -0.8 -1 '---="'--'---'---__j -1 12.5 25 -1 -0.5 0 0.5 u

[m/s]

;«

[kg/m3]

o

12.5 25

o

-0.2 -OA -0.6

-0.8

1 -1~ ___l_____ __L___ ___L___ _J -1 -0.5 0 0.5 1

u

[m/s]:

«

[kg/m"]

Figure 3A: Velocity (u) and concentration c as a function of relative depth at four different times at x=t. The Figures correspond with t=360 s,t=390 s, t=420 s and t=450 s, resp ec-tively. Adveetion of momentum is discretized in conservative form. The solid lines represent

the FGM and the dashed lines represent the SCM.

I

I

I

3

.3

S

ta

nding wave

in

ba

s

i

n

An important difference between the SCMand the FGM is the treatment of the boundary con

-ditions at the free surface and bed. The staircase representation of the free surface makes this boundary elaborate to treat in the FGM. As mentioned in section 2.2.3, inaccurate treatment

of the free surface boundary leads to excessive damping, which was clearly demonstrated du-ring the development of the FGM. In order to test the free surface administration properly, all terms except the barotropic pressure gradient and horizont al advective terms were excluded during testing. In Figure 3.5, the water surface elevation as a function of time and space

are plotted. No significant differences between the models can be observed, indicating that

in the FGM the free surface is dealt with in a correct manner. Moreover, it indicates that

I

I

I

I

I

20

I

I

I

I

I

errors introduced by the misalignment of two adjacent velocitypoints are onlyminor. The computed period of the oscillation agrees well with the thereotical period, which is computed byT=2Ljv9B ~ 100 s.

I

0 0 fgm-- ·.t=Os fgm --scm _._.- scm _._. --0.2 -0.2 -0.4 -0.4

..ê..

-0.6

8

-0.6 _{t=27 s} <:» <» -0.8 -0.8 -1 -1 -1.2 -1.2

0

100 200 300 400 500

0

100 200 300 400 500

t [sj x [mJ

I

I

I

I

I

.

1

Figure 3.5: Water elevation ( for the FGM and SCM model as a function of time (near z=Il)

and space, respectively. Horizontal lines represent vertical celledges of the computational

grid incaseof the FGM.

I

I

3

.4

Flow over a

trench

In this section, wefocus on the performance of the SCM and the FGM in case of a flow over

a trench. Two different situations are considered that are both schematizations of processes

that take placein estuaries with rapid variation of batymetry. The first one concernsa fresh water flow over a trench which is partially filled with moredense water. The concentration inthe trench is given by:

c= {25 kgjm3 if 0:::;z

< ~

o

if

!f

:::;

x:::;

H

(3.7)

I

eDiffurrorssion of mclearly.oThmeentum and mattbaroclinic pressure term in ther areneglected in the SCMemodmodeel in ordlisapproximerto brinated in two diffg out numericael -rent waysduring the computations. Besides eq. (3.6), the anti-creep method of (Stelling &

I

I

I

I

I

I

I

Van Kester, 1994) is used. The numerical model consists of a 1 km long channel. In vertical direct ion the domain is divided into 20 equidistant grid layers. In horizontal direction a grid spacing of 10 m is used. Since diffusion is neglected, transport of salt out of the trench is entirely due to the start up of the model (which generates free surface waves) and numerical errors. Initially, when the flow in the main stream has not developed yet, the anti-creep method is superior to the standard method, eq. (3.6). This can be deduced from the water surface elevation at t=O+s plotted in Figure 3.6. The elevation near the open boundary is

I

I

0.005,---,---,,---,---,

-0.0075

0

I

I

I

I

El V> -0.0025 ~\ ,~ I, ti I, " 0.0025~ : "f \ I" :I I' / I I -,...._... I I I

o

~----~

,

--~

,

---~

,

--~---

~

, I I ,'" \ I I I \ I I I \ I I I \ I J I \I II \ I I/ 11 11 11 11 \: ~' fgm , scm eq. (3.6)

==

scm anti-creep ---- --0.005 250 500 750 1000 x

[m]

Figure 3.6: Water surface elevation at t=O+ s.

I

I

I

I

caused by the boundary condition whereas the elevation in the inner area is solely due to numerical errors in the baroclinic pressure term eq. (3.6). In Figure 3.8, concentration pro-files at several positions are plotted. The salt water is washed out much faster compared to the FGM in Figure 3.7. Application of the anti-creep method for the baroclinic pressure does not improve this behaviour. Apparently, large errors are not due to the barotropic pressure discretization but are rat her introduced by the treatment of advection. This error stems from the fact that the grid is not aligned with the streamlines of the flow. Although the FGMsuffers from numerical diffusion,it is clear that results are much more satisfactory than the results obtained with the SCM.

In the next situation, the same geometry is used to simulate the propagation of a density front over a trench filled with fresh water. The error related to the treatment of the advective terms is even more clear in this case. In Figure 3.9, again velocity profiles are plotted at several positions in space. At t=300 s,the front has reached the trench. In contrast with the FGM the first grid cell in the trench is filled with salt water by adveetion from its neighbouring

cell situated just outside the trench. This type of direct transfer of matter from one grid cell to a neighbouring cell that is situated relatively far above or below this cell, is not possible

I

I

I

22

I

I

I

I

,

I

in region with rapid variations of bathymetry in the FGM. A physically more realistic flow

pattern is generated by the FGM in such cases.

0

~

h h h h

~

S

_-5 N -10

I

400 500 600 700 x [m] 0

N

l

L

l

l

~

L L

S

_-5 N -10 400 500 600 700 x [m] 0

~

l

l l l

~[

l

L

S

_-5

N -10

I

400 500 600 700 x

[m]

I

I

I

I

I

I

I

Figure 3.7: Vertical profiles of concentration obtained with the SCM, at different positions

in space at three different times; t=O+ s, t=900 s and t=1800 s, respectively. The solid line represents the anti-creep approach, the dashed line represents eq. (3.6) for the baroclinic pressure term. The maximum value at t=O+ s corresponds with a concentration of 25 kgJm3

I

I

,

I

I

I

23

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

.E..

-:5

L-l L-l

hl

L-lH2

~ -10 I 400 500 600 700 x rml

s

-:

5

h

h

_l

lld

~ -10 I 400 500 600 700 x

[m]

s

-:

5

L

•

l

lN!

~ -10

400

500

600

700

x [m]

Figure 3.8: Vertical profiles of concentration obtained with the FGM at different positions in

space at three different times; t=O+ s, t=900 s and t=1800 s, respectively. The maximum

value at t=O+ s corresponds with a concentration of 25 kg/m3

I

I

I

I

I

I

I

I

0

NI

~ ~

S

_-5 N -10 I 400 500 600 700 x [m] 0

51

~h

~ ~ S _-5 N

j"

-10 400 500 600 700 x

[m]

_~:

~'t\

s

~,

k k

~~

h h

N 400 500 600 700 x

[m]

I

I

I

I

I

I

Figure 3.9: Vertical profiles of concentration obtained with the FGMmodel (solid lines) and

SCM with anti-creep (dashed lines), at different positions in space at three different times;

t=150 s, t=300 s and t=450 s, respectively. The maximum value at t=150 s corresponds with

a concentration of 25 kg/m3

I

I

I

I

25

I

I

I

I

I

I

I

Chapter

4

Preliminary

conclusions

I

I

In this report a numerical model for shallow-water flow that is based on a fixed grid in the vertical (z-) direct ion is described (FGM). In contrast with numerical models that are based on the well-known a-transformation (SCM), the number of active layers is a function of the total water depth and thereby a function of time in this model. This means that the free

water surface is allowed to cross cell faces which does not impose arestriction

for a minimum

grid size near the surface. Depending on the local water depth, the water column is covered with a certain number of active grid cells. This flexibility may decrease computational effort significantly compared to SCMmodels in regions with large variations in the bathymetry due to the presence of tidal flats and deep channels. On the other hand, this flexibility increases the complexity of the model. Variation of the number of active grid layers makes implementation of the model more cumbersome than the SCM.

Several flow configurations were used as test cases to substantiate the proper implemen-tation of the model. Situations in which the SCM reduces to the FGM, i.e. in situations where the water depth hardly varies, both models yielded the same results. Moreover, the FGM

proved to produce results, that were very similar to the SCM results, in case of a frictionless free surface wave.

The representation of a sloping bed by a stair-case in the FGM model is known to introduce errors. Although this statement was supported in a numerical experiment, the differences in velocity compared to results from the SCM seemed to be acceptable. Significant errors were observed in the computed bed stress in the FGM. A vast reduction of the errors in the bed stress was obtained by computing the bed stress in the second grid cell above the bed instead of the first one. However,it is noted that other flowconfigurations may weaken this conclusion and therefore, additional simulations should show the errors related to this problem in more complex flows.

For simulation of exchange flowof a stratified fluid, in which gradients in the water depth are negligible, some discrepancies between the two models were observed near the head of

I

I

I

I

I

I

I

I

I

I

26

I

I

I

I

the density current. This discrepancy is ascribed to the numerical treatment of the baroclinic pressuregradient in the SCM.Since these differencesare restricted to a relatively small region only,they are of minor importance. It was found that the treatment of the advectiveterms in themomentum equation highly influence the position of the front in time. Using a conservative formulation of the advective terms in the momentum equations leads to an increasing front

celerity compared to a non-conservativeformulation. It is therefore recommended that other formulations of the advectiveterms are tested in situations werediscontinuities may playan important role.

Stratified flow in regionswith rapid variation of the bathymetry revealed significant dif -ferences between the SCMand FGM.The major contribution to these differences is not due

to the treatment of the baroclinic pressure gradient but is ascribed to errors in the adve c-tion term. This enhances the transport of salt water into or out of a trench. Therefore, the

FGMmay be a promising tool in computing stratified flowsin regions with rapidly varying

bathymetry.

As a general conclusion it can be stated that both approaches havetheir inherent dra w-backs; the SCMgivesrise to errors in the pressure term whereas in the FGMerrors are in

-troduced by the boundary treatment near the bed and free surface. Accurate simulations

require detailed and elaborate numerical methods in both approaches, just to minimize the

erros. Nevertheless, still significant differences between the two models remain, in case of

a stratified flow overa trench, for example. Further testing of the FGMrelated to practical problems, such as the stratified flow in the Haringvliet in the Netherlands, is imperative to

substantiate above contentions.

I

I

I

I

I

I

I

I

I

I

I

I

I

27

I

I

I

I

I

I

·

1

References

I

I

I

CASULLI,

v.

,

& CHENG, R.T. 1992. A semi-implicit finite difference model for three -dimensional tidal circulation. Pages 620-631 of: SPAULDING,M.L., BEDFORD,K.,

BLUMBERG,A., CHENG,R., & SWANSON,C. (eds), Estuarine and coastal modeling,

proc. 2"d conf. New York: American Society of Civil Engineers.

DELFT HYDRAULICS.1998. Delft3D-flow manual. Delft Hydraulics.

GARY,J.M. 1973. Estimate of truncation error in transformed coordinate, primitive equation

atmospheric mod

e

ls. J.

Atmo

s

. Sc

i

.

,

30

,

223

-

233.

I

HANEY,R.L. 1991. On the pressure gradient forceover steep topography insigma coordinate

ocean models. J. Phys. Oceanography,21, 610-619.

I

I

JANJlé, I.Z. 1977. Pressure gradient force and adveetion scheme used for forecasting with

steep and small scale topography. Contrib. Atmos. Phys., 50,186-199.

KRANENBURG,C. 1993. Unsteady gravity currents advancing along a horizontal surface. J.

of Hydr. Res., 31(1), 49-59.

I

I

LEENDERTSE,J.J.,

&

Lnr

,

S.K. 1975. A three-dimensional model for estuaries and coastal seas. 11 Aspects of computation. Tech. rept. R-1764-0WRT. Rand Corporation, Santa

Monica.

MESINGER,F. 1982. On the convergence and error problems of the calculation ofthe pressure

gradient force insigma coordinate models. Geophys. Astrophys. Fluid Dyn., 19, 105-117.

I

I

PAUL,J.F. 1994. Observations related to the use of the sigma coordinate transformation for

estuarine and coastal modeling studies. Pages 336-350 of: SPAULDINGM.L, ., BEDFORD,

K., BLUMBERG,A., CHENG,R., & SWANSON,C. (eds), Estuarine and coastal modeling

lIl, proc. conf. New York: American Society of Civil Engineers.

PHILLIPS,N.A. 1957. A coordinate system having some special advantages for numerical forecasting. J. Meteorology, 14, 184-185.

I

I

I

I

I

I

PODBER,D.P., & BEDFORD, K.W. 1994. Tributary loading with a terrain following coor-dinate system. Pages 475-488 of: SPAULDING,M.L., BEDFORD,K., BLUMBERG,A., CHENG,R., & SWANSON,C. (eds), Estuarine and coastal modeling 111,proc. conf. New York: American Society of Civil Engineers.

I

I

I

RODI, W. 1980. Turbuletice models and their applications in hydraulics - A state of the art review. Delft, The Netherlands: International Association for Hydraulic Research.

I

SHENG,Y.P. 1983. Mathematical modeling of three-dimensional coastal curretiis and sediment dispersion. Tech. rept. CERC-33-2. U.S. Army Corps of Engineers.

STELLING,G.S.,

&

VAN KESTER, J.A.TH.M. 1994. On the approximation of horizontal gradients in sigma co-ordinates for bathymetry with steep slopes. J. lor Num. Meth. in Fluids, 18,915-935.

STELLING,G.S., & VANKESTER,J.A.TH.M. 1996. A non-hydrostatic flow model in Car-tesian coordinates. Tech. rept. Z901-10. Delft Hydraulics Laboratory.

I

I

VAN LEER

and eonservation eombined in a seeond order seheme

, B. 1974. Towards the ultimate conservative difference scheme Il. Monotonicity

.

J.

Comp. Phys

,

32

,

101

-

136

.

I

VREUGDENHILtherlands: Kluwer Academie Publishers.,C.B. 1994. Numerical methods lor shallow-water flow. Dordrecht, The

Ne-I

WALTERS,R.A., & FOREMAN,M.G.G. 1992. A 3D, finite element model for baroclinic circulation on the Vancouver Island continental shelf. J. Marine Sys., 3.

I

WILCOX,D.C. 1993. Turbulence modeling for CFD. California: DCW Industries Inc. ZIJLEMA,M. 1997. TRIWAQ-three-dimensional incrompessible shallow flow model. National

Institute for Coastal and Marine ManagementjRIKZ. version l.O.

I

I

I

I

I

I

29

I

I

I

I

I

I

List of Figures

I

I

2.1 Grid definition in case of fgm model and scm model, respectively. 2.2 Definition of the staggered grid in the horizontal plane . . . .

8

9

I

I

3.1 Water elevation as a function of space and non-dimensional velocity, turbulent kinetic energy (TKE), dissipation of TKE and eddy viscosity as a function of depth for the FGM (solid lines) and SCM (dashed lines). . . .. 16 3.2 Velocity as a function of depth z at x=0.8L and bed shear stress as a function

of x. Solid lines represent the FGM and dashed lines represent the SCM. In the

lower two figures the bed shear stress was deterrnined in the seeond grid cell

above the bed in case of FGM. 17

3.3 Velocity (u) and concentration as a function of relative depth at four diffe-rent times at x=t. The Figures correspond with t=450 s, t=480 s,t=510

s and t=540 s,respectively. Adveetion of moment urn is discretized in non-conservative form. The solid lines represent the FGM and the dashed lines represent the SCM. . . .. 19 3.4 Velocity (u) and concentration cas a function of relative depth at four different

times at x=t. The Figures correspond with t=360 s,t=390 s, t=420 s and t=450 s, respectively. Adveetion of moment urn is discretized in conservative

form. The solid lines represent the FGM and the dashed lines represent the SCM. 20 3.5 Water elevation ( for the FGM and SCM model as a function of time (near z =Ü]

and space, respectively. Horizontal lines represent vertical cell edges of the computational grid in case of the FGM. . . .. 21 3.6 Water surface elevation at t=O+ s. . . .. 22 3.7 Vertical profiles of concentration obtained with the SCM, at different positions

in space at three different times; t=O+s, t=900 s and t=1800 s, respectively. The solid line represents the anti-creep approach, the dashed line represents eq. (3.6) for the baroclinic pressure term. The maximum value at t=O+ s corresponds with a concentration of 25 kg/rn" . . . .. 23

I

I

I

I

I

I

I

I

I

30