HYDRAULIC MODELING OF OPEN CHANNEL FLOWS

OVER AN ARBITRARY 3-D SURFACE

AND ITS APPLICATIONS

IN AMENITY HYDRAULIC ENGINEERING

TRAN NGOC ANH

August, 2006

Acknowledgements

The research work presented in this manuscript was conducted in River System

Engineering Laboratory, Department of Urban Management, Kyoto University, Kyoto,

Japan.

First of all, I would like to convey my deepest gratitude and sincere thanks to Professor

Dr. Takashi Hosoda who suggested me this research topic, and provided guidance,

constant and kind advices, encouragement throughout the research, and above all, giving

me a chance to study and work at a World-leading university as Kyoto University.

I also wish to thank Dr. Shinchiro Onda for his kind assistance, useful advices especially

in the first days of my research life in Kyoto. His efforts were helping me to put the first

stones to build up my background in the field of computational fluid dynamics.

My special thanks should go to Professor Toda Keiichi and Associate Professor Gotoh

Hitoshi for their valuable commences and discussions that improved much this

manuscript.

I am very very grateful to my best foreign friend, Prosper Mgaya from Tanzania, for all

of his helps, discussions and strong encouragements since October, 2003.

In addition, my heartfelt gratitude is extended to all of my Vietnamese friends in Japan,

Kansai Football Club members, who helped me forget the seduced life in Vietnam,

particular Nguyen Hoang Long, Le Huy Chuan and Le Minh Nhat.

Last but not least, the most deserving of my gratitude is to my wife, Ha Thanh An, and

my family, parents and younger brother. This work might not be completed without their

constant support and encouragement. I am feeling lucky because my wife, my parents and

my younger brother are always by my side, and this work is therefore dedicated to them.

iii

Abstract

Two-dimensional (2D) description of the flow is commonly sufficient to analyze

successfully the flows in most of open channels when the width-to-depth ratio is large

and the vertical variation of the mean-flow quantities is not significant. Based on

coordinate criteria, the depth-averaged models can be classified into two groups namely:

the depth-averaged models in Cartesian coordinate system and the depth-averaged models

in generalized curvilinear coordinate system. The basic assumption in deriving these

models is that the vertical pressure distribution is hydrostatic; consequently, they possess

the advantage of reduction in computational cost while maintaining the accuracy when

applied to flow in a channel with linear or almost linear bottom/bed. But indeed, in many

cases, water flows over very irregular bed surfaces such as flows over stepped chute,

cascade, spillway, etc and the alike. In such cases, these models can not reproduce the

effects of the bottom topography (e.g., centrifugal force due to bottom curvature).

In this study therefore, a depth-averaged model for the open channel flows over an

arbitrary 3D surface in a generalized curvilinear coordinate system was proposed. This

model is the inception for a new class of the depth-averaged models, which was classified

by the criterion of coordinate system. In conventional depth-averaged models, the

coordinate systems are set based on the horizontal plane, then the equations are obtained

by integration of the 3D flow equations over the depth from the bottom to free surface

with respect to vertical axis. In contrary the depth-averaged equations derived in this

study are derived via integration processes over the depth with respect to the axis

iv

perpendicular to the bottom. The pressure distribution along this axis is derived from one

of the momentum equations as a combination of hydrostatic pressure and the effect of

centrifugal force caused by the bottom curvature. This implies that the developed model

can therefore be applied for the flow over highly curved surface. Thereafter the model

was applied to simulate flows in several hydraulic structures this included: (i) flow into a

vertical intake with air-core vortex and (ii) flows over a circular surface.

The water surface profile of flows into vertical intake was analyzed by using 1D steady

equations system and the calculated results were compared with an existing empirical

formula. The comparison showed that the model can estimate accurately the critical

submergence of the intake without any limitation of Froude number, a problem that most

of existing models cannot escape. The 2D unsteady (equations) model was also applied to

simulate the water surface profile into vertical intake. In this regard, the model showed its

applicability in computing the flow into intake with air-entrainment.

The model was also applied to investigate the flow over bottom surface with highly

curvature (i.e., flows over circular surface). A hydraulic experiment was conducted in

laboratory to verify the calculated results. For relatively small discharge the flow

remained stable (i.e., no flow fluctuations of the water surface were observed). The model

showed good agreement with the observations for both steady and unsteady calculations.

When discharge is increased, the water surface at the circular vicinity and its downstream

becomes unstable (i.e., flow flactuations were observed). In this case, the model could

reproduce the fluctuations in term of the period of the oscillation, but some discrepancies

could be still observed in terms of the oscillation’s amplitude.

In order to increase the range of applicability of the model into a general terrain, the

model was refined by using an arbitrary axis not always perpendicular to the bottom

surface. The mathematical equation set has been derived and some simple examples of

v

dam-break flows in horizontal and slopping channels were presented to verify the model.

The model’s results showed the good agreement with the conventional model’s one.

vi

Preface

The depth-averaged model has a wide range of applicability in hydraulic engineering,

especially in flow applications having the depth much smaller compare to the flow width.

In this approach the vertical variation is negligible and the hydraulic variables are

averaged integrating from bed channel to the free surface with respect to vertical axis. In

deriving the governing equations, the merely pure hydrostatic pressure is assumed that is

not really valid in case of flows over highly curved bed and cannot describe the

consequences of bed curvature. Therefore, this work is devoted to derive a new

generation of depth-averaged equations in a body-fitted generalized curvilinear

coordinate system attached to an arbitrary 3D bottom surface which can take into account

of bottom curvature effects.

This manuscript is presented as a monograph that includes the contents of the following

published and/or accepted journal and conference papers:

1. Anh T. N. and Hosoda T.: Depth-Averaged model of open channel flows over an

arbitrary 3D surface and its applications to analysis of water surface profile. Journal

of Hydraulic Engineering, ASCE (accepted on May 12, 2006).

2. Anh T. N. and Hosoda T.: Oscillation induced by the centrifugal force in open

channel

flows over circular

surface. 7th

International Conference on

Hydroinformatics (HIC 2006), Nice, France, 4~8 September, 2006 (accepted on April

21, 2006)

3. Anh T. N. and Hosoda T.: Steady free surface profile of flows with air-core vortex at

vii

vertical intake. XXXI IAHR Congress, Seoul, Korea, pp 601-612 (paper A13-1),

11~16 September, 2005.

4. Anh T.N and Hosoda T.: Water surface profile analysis of open channel flows over a

circular surface. Journal of Applied Mechanics, JSCE, Vol. 8, pp 847-854, 2005.

5. Anh T. N. and Hosoda T.: Free surface profile analysis of flows with air-core vortex.

Journal of Applied Mechanics, JSCE, Vol. 7, pp 1061-1068, 2004.

viii

Table of contents

Acknowledgment

iii

Abstract

iv

Preface

vii

List of Figures

xi

List of Tables

xv

1

Chapter 1. INTRODUCTION

1.1 Classification of depth-averaged modeling

2

1.2 Depth-averaged model in curvilinear coordinates

3

1.3 Objectives of study

4

1.4 Scope of study

5

1.5 References

6

Chapter 2. LITERATURE REVIEW

7

2.1 Depth-average modeling

7

2.2 Depth-average model in generalized curvilinear coordinate system

10

2.3 Effect of bottom curvature

13

2.4 Motivation of study

16

2.5 References

16

Chapter 3. MATHEMATICAL MODEL

20

3.1 Coordinate setting

20

3.2 Kinetic boundary condition at the water surface

23

3.3 Depth-averaged continuity and momentum equations

24

Chapter 4. STEADY ANALYSIS OF WATER SURFACE PROFILE OF FLOWS

WITH AIR-CORE VORTEX AT VERTICAL INTAKE

30

30

4.1 Introduction

35

4.2 Governing equation

47

4.3 Results and discussions

54

4.4 Summary

ix

4.4 References

54

Chapter 5. UNSTEADY PLANE-2D ANALYSIS OF FLOWS WITH

AIR-CORE VORTEX

56

5.1 Governing equation

56

5.2 Numerical method

59

5.3 Results and discussions

62

5.4 Summary

64

5.5 References

65

Chapter 6. WATER SURFACE PROFILE ANALYSIS OF FLOWS OVER

CIRCULAR SURFACE

66

6.1 Preliminary

66

6.2. Hydraulic experiment

67

6.3 Steady analysis of water surface profile

74

6.4 Unsteady characteristics of the flows

81

6.5 2D simulation

94

6.6 Summary

94

6.7 References

99

Chapter 7. MODEL REFINEMENT

100

7.1 Preliminary

100

7.2 Non-orthogonal coordinate system

101

7.3 Application

105

Chapter 8. CONCLUSIONS

111

x

List of Figures

Chapter 2

Figure 2.1 Definition sketch for variables used in depth-averaged model…….. 8

Figure 2.2 Definition of terms in curvilinear system…………………………...11

Figure 2.3 Definition sketch by Sivakumaran et al. (1983)……………………..14

Chapter 3

Figure 3.1 Definition sketch for new generalized coordinate system…………..21

Figure 3.2 Kinetic boundary condition at water surface………………………..23

Chapter 4

Figure 4.1 An example of free surface air-vortex………………………………31

Figure 4.2 Various stages of development of air-entraining vortex:

S1>S2>S3>S4 (Jain et al, 1978)……………………………………31

Figure 4.3 The inflow to and circulation round a closed path

in a flow field (Townson 1991)……………………………………..33

Figure 4.4 The concept of simple Rankine vortex that including

two parts: free vortex in outer zone and forced vortex

in inner zone (Townson 1991)………………………………………33

Figure 4.5 Definition of coordinate components……………………………….36

Figure 4.6 An example of computed water surface profile with

quasi-normal depth line and critical depth line……………………..45

Figure 4.7

The effect of circulation on water surface profile and

discharge at the intake with same water head………………………48

Figure 4.8 Variation of intake discharge with circulation

(a=0.025m, b=10-5 m2, water head=0.5m)………………………….49

Figure 4.9 Different water surface profiles with different values

of circulation while maintaining the constant intake discharge……..49

Figure 4.10 Changing of water surface profile with different shape of the intake..51

xi

Figure 4.11 The effects of b on discharge (17a) and submergence

(17b) at an intake…………………………………………………51

Figure 4.12 Definition sketch of critical submergence………………………..52

Figure 4.13 Comparison of computed critical submergence by the

model (Eq. 47) and by Odgaard’s equation (51)………………….52

Figure 4.14 The variation of critical submergence wit different values of b …53

Chapter 5

Figure 5.1 Definition sketch of the new coordinates…………………………57

Figure 5.2

Illustration of the computational grid……………………………..59

Figure 5.3 Definition sketch of cell-centered staggered grid in

2D calculation……………………………………………………..60

Figure 5.4

Illustration for the discretization scheme in momentum

equations…………………………………………………………..61

Figure 5.5 Water surface of flow with different discharges at the intake…….63 Figure 5.6 Water surface of flow with different velocity at the outer-zone

boundary…………………………………………………………..63

Figure 5.7 Water surface of flow with different shape of the intake………….64

Chapter 6

Figure 6.1 Side view of the experimental facility ……………………………68

Figure 6.2

Experimental site………………………………………………….68

Figure 6.3

Schematic of sensor connection…………………………………..71

Figure 6.4 Sensor calibration…………………………………………………71

Figure 6.5

Time history of the free surface at four locations in different

experiments:

a) Exp-1; b) Exp-2; c) Exp-3; d) Exp-4;…………………72

Figure 6.6 The oscillation density at four locations in circular region………..73

Figure 6.7 Curvilinear coordinates attached to the bottom…………………....75

Figure 6.8 Illustration of computed water surface profile with

quasi-normal and critical depth lines………………………………78

Figure 6.9 Steady water surface profile with conditions of Exp-1…………….79

Figure 6.10 Steady water surface profile with conditions of Exp-2…………….79

Figure 6.11 Steady water surface profile with conditions of Exp-5…………….80

Figure 6.12 Steady water surface profile with conditions of Exp-6…………….80

xii

Figure 6.13 Illustration of staggered grid………………………………………..81

Figure 6.14 Computed water surface profile in Exp-1………………………….84

Figure 6.15 Computed water surface profile in Exp-2………………………….84

Figure 6.16 Computed water surface profile in Exp-5………………………….85

Figure 6.17 Computed water surface profile in Exp-6………………………….85

Figure 6.18 Computed water surface profile in Exp-3………………………….86

Figure 6.19 Computed water surface profile in Exp-4………………………….86

Figure 6.20 Computed water surface profile in Exp-7………………………….87

Figure 6.21 Computed water surface profile in Exp-8………………………….87

Figure 6.22

Power spectrum of water surface displacement at point 3 in Exp-3…88

Figure 6.23

Power spectrum of water surface displacement at point 4 in Exp-3…88

Figure 6.24 Comparison of calculated and experimental results at point 3

in Exp-3………………………………………………………………89

Figure 6.25 Comparison of calculated and experimental results at point 4

in Exp-3………………………………………………………………89

Figure 6.26

Power spectrum of water surface displacement at point 3 in Exp-4…90

Figure 6.27

Power spectrum of water surface displacement at point 4 in Exp-4…90

Figure 6.28 Comparison of calculated and experimental results at point 3

in Exp-4………………………………………………………………91

Figure 6.29 Comparison of calculated and experimental results at point 4

in Exp-4………………………………………………………………91

Figure 6.30

Power spectrum of water surface displacement at point 3 in Exp-8…92

Figure 6.31

Power spectrum of water surface displacement at point 4 in Exp-8…92

Figure 6.32 Comparison of calculated and experimental results at point 3

in Exp-8………………………………………………………………93

Figure 6.33 Comparison of calculated and experimental results at point 4

in Exp-8………………………………………………………………93

Figure 6.34

Carpet plot of water surface in 2D simulation of Exp-1 …………….95

Figure 6.35

Carpet plot of water surface in 2D simulation of Exp-2……………..95

Figure 6.36

Carpet plot of water surface in 2D simulation of Exp-3……………..96

Figure 6.37

Carpet plot of water surface in 2D simulation of Exp-4……………..96

Figure 6.38

Carpet plot of water surface in 2D simulation of Exp-5……………..97

Figure 6.39

Carpet plot of water surface in 2D simulation of Exp-6……………..97

Figure 6.40

Carpet plot of water surface in 2D simulation of Exp-7……………..98

Figure 6.41

Carpet plot of water surface in 2D simulation of Exp-8……………..98

xiii

Chapter 7

Figure 7.1

Illustration for limitation of the model in concave topography………101

Figure 7.2

Definition sketch of new generalized coordinate system…………….105

Figure 7.3 Calculated water surface profile at different time steps of

dam break flow in a dried-bed sloping channel:

0

0

;

;

………………………………108

0.1=

Hini

m

60=α

1

2

30=α Figure 7.4 Calculated water surface profile at different time steps of

dam break flow in a dried-bed horizontal channel: 0

0

;

;

………………………………..108

0.1=

Hini

m

60=α

1

2

0=α Figure 7.5 Comparison of water surface profile for dried horizontal

channel at different times: T = 0.4s; 1.0s; and 1.6s; ………………..109

Figure 7.6 Comparison of water surface profile for wetted horizontal

channel at different times: T = 0.4s; 0.6s; and 0.8s; ..............................109

Figure 7.7 Calculated water surface profile at different time steps of

dam break flow in a dried-bed sloping channel: 0

0

;

;

…………………………………110

Hini

0.1=

m

30=α

60=α

2

1

Figure 7.8 Calculated water surface profile at different time steps of

dam break flow in a dried-bed horizontal channel:

0

0

=

=

;

;

……………110

Hini

;0.1 m

Hini

5.0

m

30=α

60=α

1

2

up

down

xiv

List of Tables

Table 4.1

Parameters used in the calculations of results in Figure 4.7……..48

Table 4.2

Parameters used in the calculations of results in Figure 4.9……..49

Table 4.3

Parameters used in the calculations of results in Figure 4.10……51

Table 6.1 Experiment conditions……………………………………………73

xv

Chapter 1

INTRODUCTION

The advent of modern computers has had a profound effect in all branches of engineering

especially in hydraulics. The recent development of numerical methods and the

capabilities of modern machines has changed the situation in which many problems were,

up to recently, considered unsuited for numerical solution can now be solved without any

difficulties (Brebbia and Ferrante 1983).

The most open-channel flows of practical relevance in civil engineering are always

strictly three-dimensional (3D); however, this feature is often of secondary importance,

especially when the width-to-depth ratio is large and the vertical variation of the

mean-flow quantities is not significant due to strong vertical mixing induced by the

bottom shear stress. Based on these facts, a two-dimensional (2D) description of the flow

is sufficient to successfully analyze the flows in most of open channels using the

depth-averaged equations of motion. The depth averaging process used to derive these

equations sacrifices flow details over the vertical dimension for simplicity and

substantially reduces computational effort (Steffler and Jin 1993).

1

1.1 Classification of depth-averaged models:

In spite of the variation of numerical methods applied in solving the governing equations

in different practical problems, the depth-averaged models can be classified using several

criteria such as:

(1) Time dependence:

a. Steady,

b. Unsteady.

(2) Spatial integral or spatial dimension:

a. Integrate over a cross-section to get 1-D equations,

b. Integrate the 3D equations from bottom to water surface (i.e. depth averaged

model) to get 2D equations.

(3) Pressure distribution:

a. Hydrostatic pressure,

b. Consideration of vertical acceleration (Boussinesq eq.).

(5) Velocity distribution and evaluation of bottom shear stresses:

a. Uniform velocity distribution or self-similarity of distribution,

b. Modeling of local change of velocity distribution (secondary currents caused by

stream-line curvature, velocity distribution with irrotational condition, etc.).

(6) Turbulence model:

a. 0-equation model (eddy viscosity proportional to depth multiplied by friction

velocity),

b. Depth averaged

model.

ε−k

(7) Single layer model or two layered model:

A multi-layered model which has more than two layers is classified as 3-D model.

(8) Open channel flow or partially full pressurized flow:

2

a. Fully free water surface,

b. Co-existence of open channel flows and pressurized flows in underground

channels such as sewer networks.

(9) Coordinate system:

a. Cartesian coordinate set on a horizontal plane,

b. (moving) Generalized curvilinear coordinate on a horizontal plane,

c. Generalized curvilinear coordinate on an arbitrary 3-D surface.

Based on their characteristics, one model can be classified as a combination of the above

sub-classes as follows: Unsteady 2D model for open channel flows in a curvilinear

coordinate system using the assumption of hydrostatic pressure and 0-equation closure for

turbulent model.

In fact, to select the most appropriate model to solve a practical problem is usually the

most important step for the hydraulic engineers that can satisfy the required accuracy and

meet the reasonable computational cost.

1.2 Depth-averaged models and Coordinate system

The original depth-averaged model was derived in Cartesian coordinate system (Kuipers

and Vreugdenhill 1973) that was then applied extensively in various practical problems

mostly related with the regular boundaries.

To apply the depth-averaged equations to calculate the flows in the natural rivers or

meandering channels, the boundary-fitted curvilinear coordinates were introduced. The

use of curvilinear coordinates can overcome the problems of irregular boundaries

consequently it could decrease the computational cost effectively, spreads out widely and

increases the range of the applicability of depth-averaged equations.

Up to now, most of hydraulic depth-averaged models have been developed based on

3

either Cartesian coordinates or a generalized curvilinear coordinates that set on a

horizontal plane. The basic assumption in deriving these models is that the vertical

pressure distribution is hydrostatic; hence they possess the advantage of reduction in

computational cost while maintaining the accuracy when applied to flow in a channel

with linear or almost linear bottom/bed. But indeed, in many cases, water flows over very

irregular bed surfaces such as flows over stepped chute, cascade, spillway, etc and the

alike. In such cases, these models can not reproduce the effects of the bottom topography

(e.g., centrifugal force due to bottom curvature). For this reason, the present research

work is developed at proposing a new class of models that use the generalized curvilinear

coordinates based on an arbitrary surface which can conform and reflect better the

variation of bed surface; that is classified as class 9(c) in Section 1.1.

1.3 Objectives of Study

The objectives of this study are:

- To develop the new class of general depth-averaged equations in a generalized

curvilinear coordinate system set based on an arbitrary 3D surface that could be

applied to simulate the flows over a complex channel bed’s topography with highly

curvatures.

- To apply the equations to analysis of water surface profile of flows in Amenity

Hydraulic Structures.

A general mathematical model based on a new coordinate system attached to 3D arbitrary

bottom surface with an axis perpendicular to it, is developed to solve the 2D

depth-averaged equations. In this approach, the assumption of shallow water is utilized

and the internal turbulent stresses are neglected. Firstly, the initial 3D equations are

transformed into generalized curvilinear coordinate system, then, the continuity and

4

momentum equations are integrated over the depth from bottom to free water surface with

respect to the axis perpendicular to the bottom.

The model is then applied to several flows to analyze the free water surface profiles and

the results are compared with the experimental data or/and with an existing empirical

formula to show the applicability and effectiveness of the model.

1.4 Scope of Study

The manuscript consists of eight chapters including the Introduction and the Conclusions.

Chapter 1 describes the classification of depth-averaged models, difference of Coordinate

systems and clarifies the objectives of the study.

Chapter 2 reviews the related studies in the past, including studies on depth-averaged

models in Cartesian coordinates, generalized curvilinear coordinates and studies on the

effect of bottom curvature. It analyzes these methods and their limitation in applying in

flow over complex bed topography consequently it emphasizes the motivation of the

study. In Chapter 3, the derivation of new generalized depth-averaged equations from the

original Reynold Averaged Navie-Stokes (RAN) equations is described. The flow over a

vertical intake with the air-core vortex is investigated in Chapter 4 and Chapter 5, in

which the steady water surface profiles are examined firstly in Chapter 4, and then

Chapter 5 is devoted for the 2D characteristics using the unsteady equations. The same

model is applied to simulate the steady and unsteady characteristics of flows over circular

surface in Chapter 6.

Although the new model has shown its applicability in several problems, it has some

limitations as well in some convex bottoms; therefore Chapter 7 is dedicated to model

refinements. This improvement is illustrated by a simple application of dam-break flows

in horizontal and sloping channels. Finally, the Conclusions - Chapter 8 summarizes all

5

the contributions of this study.

1.5 References

1. Brebia C. A. and Ferrante A.: Computational Hydraulics. Butterworths Press, 1983.

2. Kuipers, J. and Vreugdenhill, C. B. 1973. Calculations of two-dimensional horizontal

flow. Rep. S163, Part 1, Delft Hydraulics Laboratory, Delft, The Netherlands.

3. Steffler M. P. and Jin Y. C. 1993 Depth averaged and moment equations for

moderately shallow free surface flow. J. Hydr. Res. 31 (1), 5-17.

6

Chapter 2

LITERATURE REVIEW

2.1 Depth-averaged models

Kuipers and Vreugdenhill (1973) developed one of the first mathematical models capable

of solving the two-dimensional (2D) depth averaged equations using the finite-difference

scheme developed by Leendertse (1967) (Molls and Chaudhry, 1995). Four assumptions

were the basis of the equations set used in this study: (i) water is incompressible; (ii)

vertical velocities and accelerations are negligible; (iii) wind stresses and geostrophic

effects are negligible; and (iv) average values are sufficient to describe flow properties

which vary over the flow depth.

The equations were derived by integrating the Reynold-Averaged Navier Stokes

equations from the bed to free surface with respect to vertical direction in plugging with

the kinematic boundary condition at the free surface. Therefore, this equations system is

often called as the depth-averaged or shallow water equations.

7

z

h w u

zb

x

Continuity equation:

)

)

+

+

0=

(2.1)

∂ h ∂ t

( hu ∂ x

( hv ∂ y

Momentum equation:

(

)

(

)

xy

+

+

−=

+

+

u

v

g

(2.2)

τ bx

∂ u ∂ t

∂ u ∂ x

∂ u ∂ y

∂ z s ∂ x

1 ρ h

1 ρ h

hT xx ∂ x

1 ρ h

hT ∂ y

(

)

(

)

yx

yy

+

+

+

−=

+

u

v

g

(2.3)

τ by

∂ v ∂ t

∂ v ∂ x

∂ v ∂ y

∂ z s ∂ y

1 ρ h

1 ρ h

hT ∂ x

1 ρ h

hT ∂ y

in which

vu,

: depth-averaged velocities

Figure 2.1 Definition sketch for variables used in depth-averaged model

t : time

: Cartesian coordinates

yx,

g : gravitational acceleration

z

+= h

z

h : flow depth

sz : water surface elevation (

s

b

bz : bed elevation

8

)

: bottom shear stresses ττ , by bx

TTT ,

xy

xx

yy

z

s

2

=

T

2

ρ u

2'

( ρ u

) dz

u

(2.4)

xx

1 h

∂ u ρυ ∂ x

⎡ ⎢ ⎣

⎤ ⎥ ⎦

z

b

z

s

=

=

+

ρ

T

T

' vu

'

( ρ u

)( vu

v

) dz

(2.5)

xy

yx

1 h

∂ u ∂ y

∂ v ∂ x

⎛ ρυ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

z

⎤ ⎥ ⎦

⎡ ⎢ ⎣

b

z

s

2

=

T

2

ρ v

2'

( ρ v

) dz

v

(2.6)

yy

1 h

∂ v ρυ ∂ y

⎤ ⎥ ⎦

⎡ ⎢ ⎣

z

b

,' vu

'

: turbulent velocity fluctuations (Ponce and

where υ: kinematic viscosity; and

Yabusaki 1981).

Since then, many researchers have developed and applied the 2D depth-averaged models

in variety of hydraulic engineering problems as follows:

McGuirk and Rodi (1978) studied the side discharge into an open channel and its

re-circulating flow field by using a depth-averaged model. Kalkwijk and De Vriend

(1980) calculated flow pattern in a shallow bend, in which the depth gradually approaches

to zero toward the banks by mean of a depth-averaged model. Ponce and Yabusaki (1981)

applied the shallow water model to simulate the horizontal circulation of the flows in

channel-pool configuration. Vreugdenhil and Wijbenga (1982) studied the flow pattern in

rivers with a depth-averaged model in Cartesian coordinates and introduced the curved

river boundary to the model through steps. Tingsanchali and Mahesawaran (1990)

presented a mathematical depth-averaged model for prediction of flow pattern around

groynes.

Jin and Steffler (1993) developed a depth-averaged model and included the effects of

velocity distribution in depth by using empirical velocity distribution equations. They

9

, : effective shear stresses defined as follows:

then applied their model to a 270o bend and compared the results with experimental data

for water surface profile and velocity distribution with satisfactory results. Steffler and Jin

(1993) also used the classical depth-averaged equations for shallow free surface flow

extended to treat the problem with non-hydrostatic pressure and non-uniform velocity

distributions. Khan and Steffler (1996) analyzed the momentum conservation within a

hydraulic jump utilized the depth-averaged equations; with the velocity distribution was

evaluated using a moment of longitudinal momentum equation, coupled with a simple

linear velocity distribution (Steffler and Jin 1993).

Kimura and Hosoda (1997) investigated the properties of flows in open channels with

dead zone by mean of depth averaged model in a variable grid system. The fairly good

agreement was found between the computed results and experimental data. More recently,

Molls et al. (1998) applied the 2D shallow water equations to examine the effect of

sidewall friction in analyzing a backwater profile in a straight rectangular channel.

2.2 Depth-averaged Model in Generalized Curvilinear Coordinate System

In solving these equations for practical problems, the natural rivers are almost

meandering and have the irregular boundaries, a “stair stepped” approximation is

commonly used by the conventional finite-difference method and finite volume method

(Vreugdenhil and Wijbenga 1982). This can result in either poor resolution due to coarse

grid distribution near the boundary or increase computational cost by refining the grid

system to achieve a better representation of the physical boundary. The boundary-fitted

system of generalized curvilinear coordinate technique has been developed to overcome

this deficiency that Thompson (1980) was one of the pioneers.

)ηξ,

)yx,ξξ= (

Firstly, the boundary-fitted curvilinear coordinates (

are introduced as

,

(Figure 2.2), and then the governing equations (2.1-2.3) are transformed

)yx,ηη= (

10

η

η

y

ξ

ξ

x

a)

b)

Figure 2.2 Definition of terms in curvilinear system

a) Physical mesh b) Computational mesh

) , →yx

)ηξ, (

from Cartesian to non-orthogonal curvilinear coordinates, (

, as follows:

Continuity equation:

∂ ∂ t

h J

Uh J

Vh J

⎛ ⎜ ⎝

⎞ +⎟ ⎠

∂ ⎛ ⎜ ξ ∂ ⎝

⎞ +⎟ ⎠

∂ ⎛ ⎜ η ∂ ⎝

⎞ 0=⎟ ⎠

Momentum equation:

+

gh

∂ z s ξ ∂

∂ ∂ t

M J

UM J

VM J

ξ x J

η x J

τ bx ρ J

⎛ ⎜ ⎝

⎞ +⎟ ⎠

∂ ⎛ ⎜ ξ ∂ ⎝

⎞ +⎟ ⎠

∂ ⎛ ⎜ η ∂ ⎝

⎞ −=⎟ ⎠

⎛ ⎜⎜ ⎝

∂ z ⎞ s ⎟⎟ η ∂ ⎠

(

)

)

(

)

+

+

+

+

(2.7)

2 ' hu

' hvu '

2 ' hu

)hvu ' '

( ∂ − η ∂

∂ η ∂

( ∂ − ξ ∂

∂ ξ ∂

ξ x J

ξ y J

η x J

η y J

+

gh

∂ z s ∂ ξ

τ by ρ J

∂ ∂ t

N J

UN J

VN J

ξ y J

η y J

⎛ ⎜ ⎝

⎞ +⎟ ⎠

∂ ⎛ ⎜ ∂ ξ ⎝

⎞ +⎟ ⎠

∂ ⎛ ⎜ ∂ η ⎝

⎞ −=⎟ ⎠

⎛ ⎜⎜ ⎝

∂ ⎞ z s ⎟⎟ ∂ η ⎠

2

)

(

)

)

(

+

+

+

+

(2.8)

'

hvu ' '

2 hv '

hvu ' '

)hv

( ∂ − ∂ ξ

∂ ∂ ξ

( ∂ − ∂ η

∂ ∂ η

ξ x J

ξ y J

η x J

η y J

and

where

M =

Uh

N =

Vh

: contravariant velocities that are perpendicular to ξη,

curvilinear

VU ,

coordinates, respectively

11

=

=

,

(2.9)

U

v

V

v

ξξ + u y

x

ηη + u y

x

J : Jacobian matrix defined by

J

= 1

( yx ηξ

)ξη yx

Using

the same

idea, Younus and Chaudhry (1994)

transformed

the original

depth-averaged equations from the physical coordinates (

)tyx , ,

to the computational

,ηξ

)t,

then analyzed (i) the developed uniform flow in a straight

coordinates (

rectangular channel, (ii) hydraulic jump in a diverging channel, (iii) supercritical flow in a

diverging channel as well as (iv) circular hydraulic jump. Molls and Chaudhry (1995)

presented a depth-averaged model in a general coordinate system with an explicit

algorithm and constant eddy viscosity. This model was used for both subcritical and

supercritical flows in a contraction, hydraulic jump, spur dike, and a channel with an 180o

bend. A two-dimensional, vertically averaged, unsteady circulation model using a

non-orthogonal boundary-fitted technique was employed in spherical coordinates for

predicting sea level, currents in estuarine and shelf water (Muin and Spaulding 1996). Ye

and McCorquodale (1997) developed a depth-averaged model in a non-orthogonal

curvilinear system with collocated grid arrangement. Comparison of their results with

available experimental data for side discharge into an open channel, flow in meandering

channel, and flow in Parshall flume with supercritical outflow was satisfactory. Zhou and

Goodwill (1997) presented a depth-averaged model and compared its predictions for a

water surface profile along the 180o bend channel. The bend-flow was also considered

using depth-averaged model in a generalized curvilinear coordinate by Lien et al. (1999),

Hsieh and Yang (2003), Lackey and Sotiropoulos (2005), etc. The similar equations were

also employed to investigate the fundamental characteristics of high velocity flows in a

sinuous open channel by Hosoda and Nishihama (2002). The calculated results of water

surface profile were visualized and compared with the experimental data to verify the

12

numerical model.

In a more recent paper, Zarrati et al. (2005) developed a two-dimensional depth-averaged

model in a nonorthogonal curvilinear coordinates for prediction of flow pattern in

meandering channel with 60o and 90o bend, and also with compound meandering channel.

In this study, the Cartesian velocity components were selected as the dependent variables

to avoid curvature sensitive terms (Christoffel symbols). The comparison showed that the

model predicted the water surface profile and velocity distribution well in simple

channels, and the predictions of the model in main channel of compound meandering

channel were also in general agreement with the experimental results.

2.3 Effect of Bottom Curvature

Despite of the efficiency and reasonable accuracy of above conventional depth-averaged

models, they are limited in use due to the assumption of moderate slope of the channel

bed. In all above models, the governing equations in Cartesian coordinate are firstly

integrated over the depth from the bottom to the water surface in vertical direction, (i.e.

bz to

sz which respect to a horizontal datum plane), and then, transformed into a 2D

generalized curvilinear coordinate. The basic assumption in deriving equations mentioned

above is that the vertical accelerations of a particle are negligible, (i.e. assuming the

vertical pressure distribution to be merely hydrostatic). When the channel bottom is

almost linear either horizontal or with only a small inclination, these models yield a good

solutions for shallow flows, but in complex terrains they can not describe the significant

effect of the bottom curvature.

Dressler (1978) derived a set of one-dimensional shallow water equations that included

the effect of bottom curvature. New equations were derived which is the generalization of

the nonlinear unsteady shallow-flow partial differential equations, usually called the

13

Saint-Vernant equations. The new equations possess terms containing explicitly the

curvature and the derivative of the curvature. The resulting streamlines were therefore

curved, and the pressure expression contained terms, in addition to the hydrostatic term,

which describe the effect of streamline curvatures (Equation 2.15).

)

=

+

( 1

0

(2.10)

∂ h ∂ t

∂ q ∂ s

=

+

0

(2.11)

1 g

∂ u 0 ∂ t

∂ E ∂ s

where the flow per unit width q and energy head E are

)

(2.12)

( , tsq

( 1ln

)h κ

u 0 κ

)

) 2

−≡ −

(2.13)

( tsE ,

h

cos

κ h

p h ρ g

u 2

2 ( 0 1 g

n -axis (upward normal to the bed)

),( hs

θ

tsu ),(0

h

θcos

q

s -axis

z

Radius of curvature at

κ1=P

x

Centre of curvature of the bed profile at P

Figure 2.3 Definition sketch by Sivakumaran et al (1983)

14

+≡ ζ + θ + −

and the velocity components and pressure are given by

)

( , tnsu ,

(2.14)

u 0 κ− n

1

2

)

)

)

]2

)

=

ρ

( tnsp , ,

p

( − nhg

cos

[ ( 1

κ h

( −− 1

κ n

(2.15)

h

2 0

1 ρθ u 2

where κ is the curvature of the bed, s and n are the respective coordinates along and

normal to the bed, t is the time; the height ζ above some datum defines the bed

profile, and the angle θ measures bed slope,

hp is the constant atmospheric pressure

(often set equal to zero) and ρ is the fluid density as shown in Figure 2.3 (Simakuvaran

1983).

Sivakumaran et al. (1981, 1983) applied these equations to analyze the steady state

shallow flows over curved beds such as spillway. These studies were based on the

procedure of asymptotic approximation proposed by Friedrichs (1948) and extended by

Keller (1948). However, this approach is efficient only when the flow is sufficiently

shallow and gradually varied. For this reason, these models appear to be inapplicable for

the case of highly curved surface.

Berger and Carey (1998a, b) generalized the derivation to a two-dimensional surface

based on an orthogonal curvilinear co-ordinate system, which improved the Dressler’s

equations by including the vorticity features. In fact, it is difficult to set an orthogonal

curvilinear coordinate on an arbitrary surface.

On the other hand, to apply the 2D shallow water equations for prediction of dam-break

flows over complex bed topography, Zhou et al. (2004) developed a surface gradient

method for the accurate treatment of the sources term. The ideas were originally proposed

by Zhou et al. (2001) and extended for treating bed topography involving a vertical step

(Zhou et al. 2002). However, the technique used in these studies that involved simply

15

= + + −

adding the head loss in the source term of the equations is similar to that of suddenly

contracted channel. This is not applicable in general geometry and did not reflect the

effect of the bottom curvature.

2.4 Motivation of Study

A huge number of studies have been performed in the past regarding to the conventional

depth-averaged model, however some difficulties as mentioned earlier still remain

unsolved entirely. In order to reach to a compromise that could take the advantages of

shallow water models while could represent the effect of bottom curvature; this study is

therefore devoted to develop a new general form of depth-averaged equation that can be

applied to simulate the flow over complex topography.

2.5 References

1. Berger, R. C. and Carey, G. F. 1998. Free-surface flow over curved surfaces, Part I:

Perturbation analysis. Int. J. Numer. Meth. Fluids. Vol. 28, pp. 191-200.

2. Berger, R.C. and Carey, G. F 1998. Free-surface flow over curved surfaces, Part II:

Computational model. Int. J. Numer. Meth. Fluids. Vol. 28, pp. 201-213.

3. Brebia C. A. and Ferrante A. 1983. Computational Hydraulics. Butterworths Press.

150pp.

4. Dressler, R. F. 1978. New nonlinear shallow-flow equations with curvature. J.

Hydraul. Res., Vol. 16, pp. 205-222.

5. Friedrichs, K.O. 1948. Appendix: On the derivation of the shallow water theory.

Comm. Appl. Math., Vol. 1, pp. 81-87.

6. Hosoda T. and Nishihama M. 2002. Water surface profile of high velocity flows in a

sinuous open channel. The 10th International Symposium on Flow Visualization,

August 26-29, Kyoto, Japan.

16

7. Hsieh T. and Yang J. C. 2003. Investigation on the Suitability of Two-Dimensional

Depth-Averaged Models for Bend-Flow Simulation. J. Hydr. Engrg., ASCE, Vol.

129(8), pp. 597-612.

8. Jin Y. C. and Steffler P. M. 1993. Predicting flow in curved open channels by

depth-averaged method. J. Hydr. Engrg., ASCE, Vol. 119(1), pp. 109-124.

9. Kalkwijk J. P. T. and De Vriend H. J. 1980. Computational of the flow in shallow

river bends. J. Hydraul. Res., Vol 18(4), pp. 327-342.

10. Keller, J. B. 1948. The solitary wave and periodic waves in shallow water. Comm.

Appl. Math., Vol. 1, pp. 323-339.

11. Khan A. A. and Steffler P. M. 1996. Physically based hydraulic jump model for

depth-averaged computations. J. Hydr. Engrg., ASCE, Vol. 122(10), pp. 540-548.

12. Kimura I. And Hosoda T. 1997. Fundamental properties of flows in open channel

with dead zone. J. Hydr. Engrg., ASCE, Vol. 123(2), pp. 98-107.

13. Kuipers, J. and Vreugdenhill, C. B. 1973. Calculations of two-dimensional horizontal

flow. Rep. S163, Part 1, Delft Hydraulics Laboratory, Delft, The Netherlands.

14. Lackey T. C. and Sotiropoulos F. 2005. Role of artificial dissipation scaling and

multigrid acceleration in numerical solutions of the depth-averaged free-surface flow

equations. J. Hydr. Engrg., ASCE, Vol. 131(6), pp. 476-487.

15. Leendertse, J. J. 1967. Aspects of a computational model for long period water-wave

propagation. Memo RM-5294-PR, Rand Corporation.

16. Lien H. C., Hsieh T. Y., Yang J. C. and Yeh K. C. 1999. Bend-flow simulation using

2-D depth-averaged model. J. Hydr. Engrg., ASCE, Vol. 125(10), pp. 1097-1108.

17. McGuirk, J. J. and Rodi, W. 1978. A depth-averaged mathematical model for the near

field of side discharge into open-channel flow. J. Fluid Mech. Vol. 86(4), pp.

761-781.

17

18. Molls T. and Chaudhry M. H. 1995. Depth-averaged open-channel flow model. J.

Hydr. Engrg., ASCE, Vol. 121(6), pp. 453-465.

19. Molls T., Zhao G. and Molls F. 1998. Friction slope in depth-averaged flow. J. Hydr.

Engrg., ASCE, Vol. 124 (1), pp. 81-85.

20. Muin M. and Spaulding M. 1996. Two dimensional boundary-fitted circulation model

in spherical coordinates. J. Hydr. Engrg., ASCE, Vol. 122(9), pp. 512-521.

21. Ponce V. M. and Yabusaki S. B. 1981. Modeling circulation in depth-averaged flow. J.

Hydr. Div., ASCE, Vol. 107(11), pp. 1501-1518.

22. Sivakumaran N. S, Hosking, R. J. and Tingsanchali, T. 1981. Steady shallow flow

over a spillway. J. Fluid Mech., Vol. 111, pp.411-420.

23. Sivakumaran N. S, Tingsanchali, T. and Hosking, R. J. 1983. Steady shallow flow

over curved bed. J. Fluid Mech., Vol. 128, pp.469-487.

24. Steffler M. P. and Jin Y. C. 1993. Depth averaged and moment equations for

moderately shallow free surface flow. J. Hydr. Res., Vol. 31(1), pp. 5-17.

25. Thompson J. F. 1980. Numerical solution of flow problems using body-fitted

coordinate systems. Computational fluid dynamics. W. Kollmann, ed., Hemisphere

Publishing Corp., Bristol, PA.

26. Tingsanchali T. And Maheswaran S. 1990. 2-D depth-averaged flow computation

near groyne. J. Hydr. Engrg., ASCE, Vol. 116(1), pp. 71-86.

27. Vreugdenhill C. B. and Wijbenga J. 1982. Computation of flow pattern in rivers. J.

Hydr. Div., ASCE, Vol. 108(11), pp. 1296-1310.

28. Ye J. and McCorquodale J. A. 1997. Depth-averaged hydrodynamic model in

curvilinear collocated grid. J. Hydr. Engrg., ASCE, 123(5), pp. 380-388.

29. Younus M. And Chaudhry M. H. 1994. A depth-averaged

turbulence model

ε−k

for the computation of free-surface flow. J. Hydr. Res., Vol. 32(3), pp. 415-444.

18

30. Zarrati A. R., Tamai N. and Jin Y. C. 2005. Mathematical modeling of meandering

channels with generalized depth averaged model. J. Hydr. Engrg., ACSE, Vol. 131(6),

pp. 467-475.

31. Zhou J. G., Causon D. M., Mingham C. G. and Ingram D. M. 2004. Numerical

prediction of dam-break flow in general geometries with complex bed topography. J.

Hydr. Engrg., ACSE, Vol. 130(4), pp. 332-340.

32. Zhou J. G., Causon D. M., Ingram D. M. and Mingham C. G. 2002. Numerical

solutions of the shallow water equations with discontinuous bed topography. Int. J.

Numer. Methods Fluids, Vol. 38, pp. 769-788.

33. Zhou J. G., Causon D. M., Mingham C. G. and Ingram D. M. 2001. The surface

gradient method for the treatment of source terms in the shallow water equations. J.

Comput. Phys., Vol. 168, pp. 1-25.

34. Zhou J. G. and Goodwill I. M. 1997. A finite volume method for state 2D shallow

water flows. Int. J. Numer. Methods Heat Fluid Flow, Vol. 7(1), pp. 4-23.

19

Chapter 3

MATHEMATICAL MODEL

3.1 Coordinate setting

A co-ordinate system composed of three axes (

)ζηξ , ,

is used to define the

)ηξ,

are set on a 3D arbitrary curved bottom

three-dimensional flow domain, where (

surface as a body fitted coordinates and ζ is a straight line perpendicular to the surface

as shown in Figure 3.1. The surface is expressed mathematically by equation

(

) 0 =zyxζ , ,

, which implies that

( ,ζ

) const

=,

zyx

defines a co-ordinate surface above

the bottom.

In order to introduce the notation to be used, we first consider a transformation from the

i

(

(

to the curvilinear coordinates

)zyxx i ,

,

)ζηξξ ,

,

. We define

Cartesian coordinates

the transformation matrix:

=

(3.1)

c

i j

j

∂ ξ i ∂ x

20

ζ

Arbitrary surface

η

ξ

Figure 3.1 Definition sketch for new generalized coordinate system

and its inverse

i

=

i j

c

(3.2)

∂ x ξ∂ j

x

y

x

x

y

,

,

,

,

) are denoted

The derivatives of spatial independent variables (

y , ζηξζηξ

by the following equations:

− ζηζη y

y

z

z

=

(3.3a)

x

ξ

J

− ξζξζ y

y

z

z

=

(3.3b)

x η

J

− ηξηξ y

y

z

z

=

(3.3c)

x ζ

J

− ζηζη z

x

x

z

=

y

(3.4a)

ξ

J

− ξζξζ z

x

z

z

=

(3.4b)

y η

J

21

K,

− ηξηξ z

x

x

z

(3.4c)

=y ζ

J

− ζηζη x

y

x

y

=

z

(3.5a)

ξ

J

− ξζξζ x

y

y

x

=

(3.5b)

z η

J

− ηξηξ x

y

y

x

(3.5c)

=z ζ

J

, zyx ,

are the Cartesian co-ordinates and the Jacobian of transformation ( J ) is

where:

defined as:

x

=

(3.6)

J

x

ξξξ z y ηηη y z ζζζ y z

x

The transformation relationship between the contravariant velocity vector

WVUV ,

(

,

)

in curvilinear coordinate and its counterpart

in Cartesian coordinate system is

( )wvu , ,v

j

j

V = i

v = i

i c V j

(3.7)

, or

i c v j

The conditions for the axis ζ to be perpendicular to the surface (i.e. the axis ζ

perpendicular to both axes ξ and η) are expressed as follows:

)

( ζ

( ξ

grad

grad

x,y,z

(3.8)

, , zyx

) 0 =

)

( ζ

( η

(3.9)

) 0 =

grad

zyx , ,

grad

x,y,z

or:

+

+

0=

(3.10)

ξζξζξζ y z

x

y

x

z

+

+

0=

(3.11)

ηζηζηζ y z

x

y

x

z

And condition for ζ axis to remain a straight line with the scale factor equal 1 is

described as

22

2

2

2

+

+

=

y

z

1

(3.12)

x ζ

ζ

ζ

Plugging (3.10, 3.11) with (3.1, 3.2), it yields:

+

+

0=

(3.13)

xx ζη

yy ζη

zz ζη

+

+

0=

(3.14)

xx ζξ

yy ζξ

zz ζξ

=

=

=

From the definition in Equations (3.3-3.5) we can derive ( zyJ

( zyJ

( zyJ

(3.15)

zy ηζ

zy ζξ

ξζ

ζη

ηξ

)ξη zy

) ζ , x

) η , x

ξ x

=

=

=

( xzJ

( xzJ

( xzJ

(3.16)

ζη

xz ηζ

ξζ

xz ζξ

ηξ

)ξη xz

ξ y

) η , y

) ζ , x

=

=

=

( yxJ

( yxJ

( yxJ

(3.17)

ζη

yx ηζ

ξζ

yx ζξ

ηξ

)ξη yx

ξ z

) η , z

) ζ , z

Using Equations (3.15-3.17), it easily gives

( i

)zyx ,

,=

(3.18)

∂ ζ ∂

ξ i J

η i J

ζ i J

∂ ⎛ ⎜ ξ ∂ ⎝

∂ ⎛ ⎜ η ∂ ⎝

⎞ +⎟ ⎠ ⎞ +⎟ ⎠ ⎛ ⎜ ⎝ ⎞ 0=⎟ ⎠

3.2 Kinematic boundary condition at the water surface

ζ

h

ξ

Figure 3.2 Kinetic boundary condition at water surface

are

The free water surface (Figure 3.2) at a time t and

t ∆+ t

=

h ηξζ

,(

), t

0

−∆+

∆+

∆+

h

t

ηηξξζζ ,

(

,

=∆+ ) t

0

23

thus

=

+

+

(3.19)

∂ h ξ ∂

∂ h η ∂

ζ ∆ ∆ t

∂ h ∂ t

ξ ∆ ∆ t

η ∆ ∆ t

where

=

+

+

=

+

+

=

u

v

w

U

(3.20)

ξ x

ξ y

ξ z

ξ x

s

ξ y

s

ξ z

s

s

ξ ∆ ∆ t

∆ x ∆ t

∆ y ∆ t

∆ z ∆ t

and similarly,

=

(3.21)

sV

∆η ∆ t

=

(3.22)

sW

∆ζ ∆ t

are contravariant components of velocity vector at free water

in which

,

)

, ( WVU s

s

s

surface.

Applying Equations (3.20-22) into (3.19), the kinetic boundary condition at the free

surface is obtained

=

+

+

W

U

V

(3.23)

s

s

s

∂ h ξ ∂

∂ h η ∂

∂ h ∂ t

3.3 Depth-averaged continuity and momentum equations

The fundamental equations in Cartesian co-ordinates are firstly transformed into an

arbitrary generalized curvilinear coordinate system as follows:

Continuity equation:

(3.24)

∂ ξ ∂

∂ ζ ∂

U J

V J

W J

⎛ ⎜ ⎝

⎞ +⎟ ⎠

∂ ⎛ ⎜ η ∂ ⎝

⎞ +⎟ ⎠

⎛ ⎜ ⎝

⎞ 0=⎟ ⎠

Momentum equation:

24

+

+

+

J

J

J

x

x

x

+

+

+

+

+

+

Jv

J

J

J

y

y

∂ ζ ∂

∂ ξ ∂

∂ ∂ t

+

+

+

Jw

( Uu ( Uv ( Uw

) ρξ p ) ρξ p y ) ρξ p

J

( Vu ( Vv ( Vw

) ρη p ) ρη p ) ρη p

J

( Wu ( Wv ( Ww

) ρζ p ) ρζ p ) ρζ p

J

Ju ⎡ ⎢ ⎢ ⎢ ⎣

z

z

z

+

+

+

+

g

J

z

x

y

xz

xx

xx

xz

x

y

z

x

+

+

+

+

+

+

=

g

J

y

z

x

y

yx

yz

yx

yz

y

x

z

∂ ξ ∂

) ) )

) ) )

+

+

+

+

g

z

zz

x

y

z

⎤ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ⎢ ⎢ η ∂ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎦

zx

zz

y

x

z

+

xx

xz

y

x

z

+

+

+

(3.25)

yx

yz

x

y

z

∂ ζ ∂

( ρτξτξτξ J xy ( ρτξτξτξ J yy ( ρτξτξτξ J zy ) ) )

+

+

⎤ ⎥ ⎥ ⎥ ⎦ ( ρτητητη J xy ( ρτητητη J yy ( ρτητητη J zy ⎡ ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ⎢ ⎢ η ∂ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎦

zx

zz

x

y

z

Based on the shallow water condition, the transforming Jacobian is approximately

homogenous with respect to ζ-axis and is set equal to the value of Jacobian at the

bottom

0J . Other flow quantities such as U and V are uniform in ζ -direction.

Integrating Equation (3.24) over the depth from the curved bottom to the water surface

with respect to ζ-axis and applying the kinetic boundary condition at the water surface

(Equation 3.23), the depth-averaged continuity equation is derived:

+

+

=

0

(3.26)

∂ ξ ∂

∂ η ∂

1 J

∂ h ∂ t

M J

N J

0

0

0

M =

Uh

N =

Vh

and

. From now onwards, the subscript zero “0” denotes the

where

values right on the bottom surface.

Substituting the relations described in Equation (3.7) into Equation (3.25), the ξ

component of momentum equation in 3D space is recast

2

2

Γ

+

( U

UV

UW

VU

ξ +Γ ξξ

ξ +Γ ξη

ξ +Γ ξζ

ξ ηξ

∂ ζ ∂

∂ ∂ t

U J

U J

VU J

WU J

1 J

⎤ ⎥ ⎥ ⎥ J ⎦ zx ( ρτζτζτζ + J xy ( ρτζτζτζ J yy ( ρτζτζτζ J zy ⎡ ⎢ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎥ ⎦

∂ ⎛ ⎜ η ∂ ⎝

∂ ⎛ ⎜ ⎜ ξ ∂ ⎝

ξ

ξξ

2

2

+

Γ

Γ

V

VW

WU

WV

W

-

(pressure

term)

) =

+

ξ +Γ ηη

ξ +Γ ηζ

ξ +Γ ζξ

ξ + ζη

ξ ζζ

τ ρ J

G J

∂ ⎛ ⎜ ⎜ ξ ∂ ⎝

⎛ ⎜ ⎝ ⎞ +⎟ ⎠ ⎞ +⎟ ⎠ ⎛ ⎜ ⎝ ⎞ +⎟ ⎠ ⎞ ⎟ ⎟ ⎠

25

⎞ ⎟ ⎟ ⎠

ξη

ξζ

ξξ

ξη

ξζ

ηξ

ηη

ηζ

+

( τ

τ

τ

τ

τ

τ

+

+

ξ +Γ ηη

ξ +Γ ξξ

ξ +Γ ξη

ξ +Γ ξζ

ξ +Γ ηξ

ξ Γ ηζ

∂ ζ ∂

τ ρ J

τ ρ J

1 ρ J

∂ ⎛ ⎜ ⎜ η ∂ ⎝

ζξ

ζη

ζζ

+

τ

τ

τ

Γ

(3.27)

ξ +Γ ζξ

ξ +Γ ζη

ζζ

where

)

)

+

+

+

+

+

(pressure

= term)

( ξξξξξξ y z

x

x

y

z

( ηξηξηξ y z

x

y

x

z

1 J

1 J

∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρη ∂ ⎠ ⎝

∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρξ ∂ ⎠ ⎝

)

+

+

+

+

y

y

x

x

z

z

ξ x

∂ ζ ∂

( ζξζξζξ1 J

ξ x J

η x J

ζ x J

⎛ ⎜ ⎝

⎞ ⎟ ⎠

∂ ⎛ ⎜ η ∂ ⎝

⎞ +⎟ ⎠

⎞ +⎟ ⎠

∂ ⎛ ⎞ p ⎟⎟ ⎜⎜ ρζ ∂ ⎝ ⎠

⎫ ⎬ ⎭

∂ ⎧ p ⎛ ⎜ ⎨ ξρ ∂ ⎝ ⎩

+

+

+

+

ξ y

ξ z

∂ ∂ ζ

∂ ∂ ζ

ξ y J

η y J

ζ y J

ξ z J

η z J

ζ z J

⎛ ⎜ ⎝

⎞ ⎟ ⎠

∂ ⎛ ⎜ ∂ η ⎝

⎞ +⎟ ⎠

⎞ +⎟ ⎠

⎫ ⎬ ⎭

∂ ⎧ p ⎛ ⎜ ⎨ ∂ ξρ ⎝ ⎩

⎛ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ η ⎝

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

⎫ ⎪ ⎬ ⎪⎭

∂ ⎧ ⎛ p ⎪ ⎜⎜ ⎨ ∂ ξρ ⎪⎩ ⎝

(3.28)

Applying the perpendicular condition (Equations 3.10, 3.11) and relation in Equation

(3.18) to Equation (3.28) gives:

2

2

2

)

)

+

+

+

+

+

(pressure

= term)

( ξξξ y z

x

( ηξηξηξ y z

x

y

x

z

1 J

1 J

∂ ⎛ ⎞ p ⎜⎜ ⎟⎟ ρη ∂ ⎝ ⎠

∂ ⎛ ⎞ p ⎜⎜ ⎟⎟ ρξ ∂ ⎝ ⎠

(3.29)

Substituting pressure term in Equation (3.29) into (3.27) and taking integration over the

depth from the curved bottom to the water surface with respect to ζ -axis, the

depth-averaged momentum equation in ξ direction is obtained:

2

2

+

+

+

+

+

+

( M

MN

NM

N

ξ ξξ

ξ ξη

ξ ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ∂ t

M J

UM J

VM J

0

0

0

1 hJ 0

⎛ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ ξ ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ η ⎝

⎞ ⎟⎟ ⎠

h

h

+

+

+

+

ξ

0

0

0

0

0

x

x

y

z

2 x 0

=

ζ

G

ξ τ b ρ J

ηξηξηξ 0 y z J

ξξξ 2 2 y z 0 0 J

h J

0

0

0

0

∂ p ∫ d ∂ ρη 0

∂ p ζ ∫ d ∂ ρξ 0

h

h

τ

τ

+

+

ζ

+

+

(3.30)

d

ζ d

y

x

x

z

( ζξζξζξ y z 0

0

0

0

0

0

ξξ ρ

ξη ρ

∂ ξ ∂

∂ η ∂

1 J

1 J

1 J

0

) ⎛ ⎜⎜ ⎝

⎞ p ⎟⎟ ρ ⎠

00

00

b

26

⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠

Using the same procedures the η-component of depth-averaged momentum equation is

derived:

2

2

+

+

+

+

+

+

( M

MN

NM

N

η ξξ

η ξη

η ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ∂ t

N J

UN J

VN J

0

0

0

1 hJ 0

⎛ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ξ ∂ ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ η ∂ ⎝

⎞ ⎟⎟ ⎠

h

h

+

+

+

+

η

0

0

0

0

0

x

x

y

z

2 x 0

ζ

=

G

η τ b ρ J

ηξηξηξ 0 y z J

h J

ηηη 2 2 y z 0 0 J

0

0

0

0

∂ p ζ ∫ d ∂ ρξ 0

∂ p ∫ d ρη ∂ 0

h

h

τ

τ

+

+

+

+

ζ

d

ζ d

0

0

0

0

0

( ζηζηζη y z 0

y

x

x

z

∂ ∂ ξ

ηξ ρ

∂ ∂ η

ηη ρ

1 J

1 J

1 J

0

) ⎛ ⎜⎜ ⎝

⎞ p ⎟⎟ ρ ⎠

00

00

b

(3.31)

where

ξG , ηG , ξτ

b and ητ

b are the contravariant components of gravitational vector

and shear stress acting on the bottom surface;

Γ is the Riemann-Christoffel symbols:

i jk

m

−=Γ

,

(m =

1,2,3)

(3.32)

i jk

m

∂ ∂ x ∂ ∂ ξξ k j

∂ ξ i ∂ x

⎞ ⎟ ⎟ ⎠

⎛ ⎜ ⎜ ⎝

Substituting the conditions of orthogonality (Equations 3.10 and 3.11) into Equations

(3.30, 3.31), and neglecting the effect of Reynold stress gives:

2

2

+

+

+

+

+

+

( M

MN

NM

N

ξ ξξ

ξ ξη

ξ ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ∂ ξ

∂ ∂ t

M J

UM J

VM J

0

0

0

1 hJ 0

⎛ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

⎛ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ η ⎝

⎞ ⎟⎟ ⎠

h

h

+

+

+

+

ξ

x

x

y

z

0

0

0

0

0

2 0 x

=

ζ

G

h J

ξξξ 2 2 0 0 y z J

ηξηξηξ 0 y z J

ξ τ b ρ J

0

0

0

0

∂ p ∫ d ∂ ρη 0

∂ p ζ ∫ d ∂ ρξ 0

(3.33)

2

2

+

+

+

+

+

+

( M

MN

NM

N

η ξξ

η ξη

η ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ξ ∂

∂ ∂ t

N J

UN J

VN J

0

0

0

1 hJ 0

⎛ ⎜ ⎜ ⎝

⎛ ⎜ ⎜ ⎝

∂ ⎛ ⎜ ⎜ η ∂ ⎝

⎞ ⎟ ⎟ ⎠

⎞ ⎟ ⎟ ⎠

⎞ ⎟ ⎟ ⎠

h

h

+

+

+

+

η

x

x

y

z

2 0 x

0

0

0

0

0

ζ

=

ζ d

d

G

η τ b ρ J

ηξηξηξ 0 y z J

h J

ηηη 2 2 0 0 y z J

0

0

0

0

∂ p ∫ ρη ∂ 0

∂ p ∫ ρξ ∂ 0

(3.34)

27

Pressure term in the momentum equation in ζ-direction is

)

)

+

+

+

+

+

(pressure

term)

=

( ζξζξζξ y z

x

x

y

z

( ζηζηζη y z

x

x

y

z

1 J

1 J

∂ ⎛ ⎞ p ⎟⎟ ⎜⎜ ρξ ∂ ⎝ ⎠

∂ ⎛ ⎞ p ⎟⎟ ⎜⎜ ρη ∂ ⎝ ⎠

)

+

+

+

( ζζζ 2 2 y z

2 x

1 J

∂ ⎛ ⎞ p ⎟⎟ ⎜⎜ ∂ ρη ⎝ ⎠

Applying the orthogonal conditions, it is reduced as:

)

+

+

(pressure

= term)

( ζζζ 2 2 y z

2 x

1 J

∂ ⎛ ⎞ p ⎜⎜ ⎟⎟ ∂ ρη ⎝ ⎠

To estimate the pressure terms in Equations (3.33), (3.34), the momentum equation in

ζ -direction of Equation (3.25) is considered. If the acceleration and shear stress

components are negligible, the equation can be reduced as:

+

ζζζζ + 2 z

2 x

2 z

)

Γ

=

UV

VU

VV

( UU

(3.35)

ζ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

∂ p ρζ ∂

G J

J

1 J

0

0

0

therefore,

p

)

=

[ −

Γ

+

( UU

UV

VU

VV

G

ζ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

+

+

∂ ∂ ζζζρζ 2 2 x z

1 2 z

+

+

=

ζ is a straight axis consequently,

1

,

ζζζ 2 2 z z

2 x

ζ

)

−=

Γ

+

(3.36)

( UU

UV

VU

VV

G

ζ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

∂ p ∂ ρζ

The expression for pressure distribution is obtained by integrating Equation (3.36) with

respect to ζ. It follows that:

2

ζ

}

=

+

+

+

{ UU

VU

UV

VV

G

(3.37)

ζ ξη

ζ ηη

Γ 0

Γ 0

ζ Γ ξξ 0

ζ Γ ξη 0

ph ζ ∫ d ρ

h 2

0

Equations (3.26, 3.33, 3.34) and (3.37) are the depth-averaged equations in a general form

used in calculating water surface profile of flows in the following applications.

28

Notations

t : time

, zyx ,

: Cartesian coordinates

ζηξ ,

,

: generalized curvilinear coordinates

,

,

,

,

: components of transformation matrix

ηηηξξξ , z z

y

x

y

x

J : transformation Jacobian

: velocity components in Cartesian coordinate

)wvu , ( ,

: contravariant components of velocity vector

,

)WVU , (

)

,

: contravariant components of velocity vector at water surface

( , WVU s

s

s

h : water depth in ζ-direction

p : pressure

ρ: water density

g : gravitational acceleration

: contravariant components of gravitational vector

ζηξ ,

GGG

,

τττττ

,

,

,

,

: effective shear stress

xx

xy

yx

xz

K,

: contravariant of shear stress acting on the bottom

ηξ ττ b , b

Γ : Christoffel symbol

i jk

29

Kyy

Chapter 4

STEADY ANALYSIS OF WATER SURFACE PROFILE OF FLOWS WITH AIR-CORE VORTEX AT VERTICAL INTAKE

4.1 Introduction

Air-core vortex formation at intakes is a significant hydraulic engineering problem in

many situations such as intakes for irrigation, drainage system, hydropower generation in

which the water is normally drawn from rivers, channels or reservoir,… as well as in

water art-works or sculptures. It occurs typically whenever the submergence is less than a

critical value and causes some detrimental effects as reduction in intake discharge,

resulting vibrations and noises as well as operational difficulties. Figure 4.1 shows an

example of free surface air-core vortex occurred in laboratory.

When the flow in a large body of slowly moving water is diverted and locally accelerated

or drawn off, any associated vortex tube is extended and its rotation is thereby increased.

Higher velocities incur lower pressures and, if a free surface exists, it becomes locally

30

Figure 4.1 An example of free surface air-vortex

S1

S2

Q

Q

b. Formation of a dimple

a. No depression on the surface at large submergence

S3

S4

Q

c. Air core extends deeper as the vortex becomes

d. Air-entraining vortex

stronger

Figure 4.2 Various stages of development of air-entraining vortex: S1>S2>S3>S4

(Jain et al., 1978)

31

depressed. Thus in hydraulic structures, where flow negotiates a vertical shaft, submerged

conduit or gated outlet, the vortex may be sufficiently intense to form a hollow core,

which could transmit air. This may lead to the undesirable effects already mentioned, and

considerable efforts have been made to predict and control such phenomena. Indeed, the

deliberated introduction of a vortex at the vertical intake has been described in Jain et al.

(1978).

Figure 4.2 shows the various stages in the development of an air-entraining vortex when

the water depth is decreased gradually. At the first stage, when the submergence is large,

there is no depression on the surface (Figure 4.2a), but together with decreasing water

depth, a dimple forms as in Figure 4.2b. If the submergence is further decreased the

air-core vortex occurs (Figure 4.2c) and Figure 4.2d shows critical condition under which

a vortex just tends to entrain air (Jain et al. 1978).

Suppose that steady flow occurs in the horizontal

y

x − plane and that closed path, S ,

is traced out in the fluid (Figure 4.3). The path encloses an area A by linking adjacent

flow particles and subsequently moves with the flow. At a point on the path, one particle

may have velocity components

nV and

tV normal and tangential to the path. For the

is

path to remain unbroken and assuming the flow is incompressible, the sum ∑ ∆SVn

the net inflow into the closed area A and must always be zero. The corresponding sum

is not necessarily zero and is known as the circulation Γ

along the path ∑ ∆SVt

(being defined as positive anticlockwise). The Kenvin’s theorem states that the circulation

remains constant with time unless an external shear stress exists along S , as described in

Townson (1991).

32

S

A

nV

S∆

tV

Figure 4.3 The inflow to and circulation round a closed path in a flow field (Townson 1991)

r

∞→r 0→ tV

H∆

H

h

free vortex

forced vortex

ideal fluid

real fluid

Figure 4.4 The concept of simple Rankine vortex including two parts: free vortex in

outer zone and forced vortex in inner zone (Townson 1991)

33

If we consider an infinitely wide reservoir of constant depth is drawn off at a central point

with rate q per unit depth. At some distance from the center, a circular vortex tube may

be defined with circulation Γ about a vertical axis through the drawn-off point. As the

tube contracts to radius r under the influence of drawn-off, the tangential velocity

. The radial velocity is

π2/= r

q

, and both clearly accelerate, producing a

π2/Γ= r

Vt

Vr

spiral flow towards the center. Applying the Kelvin’s theorem, the circulation is

conserved therefore the velocity at the center is infinite. If the pressures remain

hydrostatic, the depression of the free surface is also infinite. But the viscosity present in

real fluids prevents this condition arising, and a zone in center of the flow rotates as a

solid mass (Townson 1991). This central part is known as a forced vortex and the

composite system as Rankine’s vortex (Figure 4.4).

Several approaches have been presented in the literature to deal with the problem of

determination and prediction of critical submergence serving in design works. These

approaches basically can be labeled as analytical models (Jain 1984; Odgaard 1986; Hite

and Mih 1994) and physical models (Anwar et al. 1978; Jain et al. 1978; Yildirim and

Kocabas 1995; Yildirim and Kocabas 1997; Yildirim and Kocabas 1998).

Many analytical attempts have been presented in the literature in order to attain a

theoretical view of the far-field velocity; in fact the flow representation has not been

defined so far by any comprehensive analytical analysis. The concept of simple Rankine

vortex normally used in the basic equations (Odgaard 1986; Hite and Mih 1994),

therefore this approach obviously could not be applied for the case of air-entraining

vortex and moreover, the Kelvin’s theorem is invalid for the central region.

Trivellato et al. (1999) set the water surface equal to the stationary headwater while other

experimental works only focused on the critical submergence. Consequently, these

approaches could not be used to predict the water surface profile of flow with air-core

34

vortex.

In this study, the water surface profile of a steady air core vortex flow into a vertical

intake was derived through out a depth-averaged model of open channel flows over the

3-D curvilinear bottom using a generalized and body fitted coordinate system. The

assumption of fully free air-core vortex in the new coordinate allowed us to use the

Kelvin’s theorem of the conservation of circulation for the whole flow field. The vortex

was assumed axisymmetric and steady. The assumption of shallow water and kinetic

boundary condition at water surface were also used. The equation describing water

surface profile was derived and calculated results were compared to the formula

introduced by Orgaard (1986).

The application’s results showed us the ability of the model in analyzing the water surface

profile and can be improved to simulate the flow structure of an air-core vortex.

4.2 Governing equations

Coordinate setting

To consider the vortex occurring at a cavity on the bottom surface (Figure 4.5), the

where (

position of any point say P , is defined by three coordinates (

)ζηξ , ,

)ηξ,

define the position of

0P (projection of P ) on the bottom, and ζ is the distance from

point P to that bottom surface.

Assuming that the shape of bottom surface (i.e at

0=ζ ) has the form of

=

z

,

(4.1)

0

− b −

a

)

r ( 0

where a and b are the coefficients which define the shape of the bottom;

0r is the

distance from any point on the bottom to the z -axis

2

2

=

+

x

y

(4.2)

r 0

0

0

35

r

y

z • l

ξ

a

x

η

a

a−

x

O

• ζ P

P0

r

d

a

2=

Figure 4.5 Definition of coordinate components

36

−=

y

Based on Figure 4.5, we can derive the relations

ηcos

and

.

x = 0

r 0

0

ηsin0 r

Then, the bottom surface is expressed by the following equation:

Φ

=

(

,

,

)

(

0

x

y

z

za )

=+ b

(4.3)

0

0

0

r 0

0

)(n normal to the bottom surface is derived as:

The unit vector

a

1

1

=

+

+

=

r n

r i

r j

r k

2

2

2

Φ Φ

2

2

2

grad grad

+

+

r 0 +

zx 00 r 0

zy 00 r 0

z

a

z

a

z

a

(

)

(

)

(

)

0

r 0

0

r 0

0

r 0

At ζ the relations between

), ,( zyx

and

ζηξ ) ,( ,

can be expressed as:

η

b

1

cos

ζ

η −

=

+

(4.4a)

x

x

ζ =

cos

r 0

0

4

2

2

2

+

+

zx 00 r 0

a

b

(

)

z

a

(

)

r 0

0

r 0

η

b

1

sin

ζ

ζ

=

+

−=

η +

(4.4b)

y

y

sin

r 0

0

4

2

2

2

+

+

zy 00 r 0

a

b

(

)

z

a

(

)

r 0

0

r 0

2

a

a

(

)

=

+

−=

+

ζ

ζ

z

z

(4.4c)

0

4

2

2

b −

2

a

(

)

r 0 −

+

r 0 +

r 0

a

b

(

)

z

a

(

)

r 0

0

r 0

Then the covariant base vector components on the bottom surface are denoted as follows:

η

=

+

+

−=

+

j

e

y

j

z

k

η r i

k

(4.5a)

ξ

ξ

v ix ξ 0

0

ξ 0

2

b −

a

)

cos p 1

sin p 1

rp ( 1 0

=

+

+

−=

y

j

z

k

r η i

j

sin

cos

(4.5b)

e η

η η

r ix η 0

η 0

η 0

r 0

r 0

η

η

=

+

+

−=

+

x

r i

y

j

z

k

r i

r j

e ζ

ζ 0

ζ 0

ζ 0

4

2

4

2

cos b −

+

sin b −

+

a

b

a

b

)

(

)

(

r 0

r 0

2

a

(

)

+

r k

(4.5c)

4

2

r 0 −

+

a

b

(

)

r 0

Let’s denote that

2

+

+

)

=

(4.6)

p 1

2

( rl 0 2

) a 2

2

+

a +

br 0 −

+

a

l

a

a

b

(

)

(

( rb 0 − 2 )

(2 lb

)

r 0

r 0

r 0

r 0

37

in which l is defined in Figure 4.5.

From Figure 4.5, it could be easily seen that

+

l

b −

l

z

a

0

=

=

tan ξ =

(4.7)

) a −

+ b )a

r 0 r 0

r 0

( rl 0 ( rr 0 0

then

2

+

+

a

a

)

)

( rl 0

br 0

−=

=

(4.8)

2

+ b )

) a −

a

( rb 0 − a

)

(

∂ ξ tan ∂ r 0

∂ ∂ r 0

( rl 0 ( rr 0 0

2 r 0

r 0

⎞ ⎟⎟ ⎠

⎛ ⎜⎜ ⎝

From (4.6) we can obtain:

2

2

2

2

)

)

+

+

ξ

a

l

a

b

2 lb

2

( r 0

2 r 0

( r 0

=

=

ξ

+= 1

tan

tan ξ ∂

ξ

+ )2

) a −

1 2 cos

a

( r 0 2 r 0

( r 0

(4.9)

but

ξ

ξ

ξ

tan ξ ∂

=

=

ξ

tan ξ ∂

∂ r 0 ξ ∂

∂ r 0 ξ ∂

tan ∂ r 0

tan ∂ r 0

(4.10)

Thus, substituting (4.8) and (4.9) into (4.10) gives

2

2

2

2

2

+

+

+

(

a

)

l

a

)

a

)

b

r 0

r 0

−=

−=

∂ r 0 ∂ ξ

r ( 0 + 2

(2 lb +

)

a

a

)

1 p 1

rb ( 0

r 0 br 0

rl ( 0

(4.11)

From (4.4a):

3

cos

=

+

η

cos

(4.12)

4

) a ] ζ 232

( η r 0 ) +

∂ x ∂ r 0

b

a

b 2 [ ( r 0

Using the relations in (4.11), (4.12), it can be obtained:

3

3

η

cos

cos

−=

=

+

η

x

cos

4

4

∂ x = ξξ ∂

( η r 0 )

) a ] ζ 232

) a ] 232

( η r 0 ) +

+

cos p 1

1 p 1

⎛ −⋅ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

a

b

b

a

b 2 [ ( r 0

2 b [ ( rp 0 1

⎡ ⎢ ⎢ ⎣

⎤ ζ ⎥ ⎥ ⎦

(4.13)

38

Similarly,

3

η

sin

=

+

y

(4.14)

ξ

4

) a ] ζ 232

( η r 0 )

+

sin p 1

a

b

b 2 [ ( rp 1 0

a

−=

z

(4.15)

ξ

2

4

) 232

)

b −

a

] ζ

( 2 rb 2 0 )

+

( rp 0 1

a

b

[ ( rp 0 1

and

η

b

ζ

−=

η

+

(4.16)

sin

x η

r 0

4

2

sin )

+

a

b

( r 0

η

b

ζ

−=

η

+

(4.17)

cos

y η

r 0

4

2

cos )

+

a

b

( r 0

(4.18)

0=ηz

η

b

−=

x

(4.19)

ζ

4

2

cos )

+

a

b

( r 0

η

b

=

y

(4.20)

ζ

4

2

sin )

+

a

b

( r 0

2

)

z

(4.21)

4

2

a )

( r 0 −

+

a

b

( r 0

and the Jacobian right on the bottom surface is:

ζηξ 0 0

0

4

2

ζηξ 0 0

0

2

0

= ζ

0

ζηξ 0 0

0

x x x + − ) b ( 1 r 0 = = = (4.22) y y y a − 1 J J ( a ) 1 p 1 r 0 r 0 z z z

The gravitational vector can be expressed in the new coordinate system as follows:

ξ eGG

η eG

ζ eG

ξ

η

ζ

= + + (4.23)

39

Substituting the base vector components from (4.5) into (4.23)

(

)

(

)

(

)

ξ ixGG ξ

η ixG η

ζ ixG ζ

)

= + + + + + + + + jy ξ kz ξ jy η kz η jy ζ kz ζ

( ξ xG

η xG

) ζ ixG

( ξ yG

η yG

ζ yG

η zG

ζ zG

( ξ zG

) k

η

ξ

ζ

η

ζ

ξ

η

ζ

ξ

= + + + + + + + + j

G −=

kg

but

ξ

consequently,

η

ξ

ζ

ξ

G x x x ζηξ 0 ⋅ = y y G (4.24) y η 0 − g z z G z ζηξ

2

ξ

)

hence, the contravariant components of gravitational vector are derived

( yxgJ

ζη

2

( r 0 −

) +

( r 0

η

−=

−= − = (4.25) G yx ηζ bgp 1 a 4 − ) a b

G

( yxgJ

) 0=

ζξ

yx ξζ

2

)

(4.26)

ξ

)

( yxgJ

ηξ

4

2

( r 0 −

( r 0

− −= − −= (4.27) G g yx ξη a ) + a b

4

Using Equations (4.13-4.22) in plugging with (3.15-3.17), it gives:

ξ x

2

( r 0 −

) +

η ( r 0

4

−= (4.28) p 1 cos a 4 − ) a b

ξ y

2

( r 0 −

) +

η ( r 0

2

= (4.29) p 1 sin a 4 − ) a b

z

4

2

) a +

( rb 0 − a

( r 0

η

−=ξ (4.30) − ) b

η −= x

sin r 0

η

(4.32)

η −= y

cos r 0

40

(4.31)

η 0=

(4.33)

z

η

b

−=

(4.34)

ζ x

4

2

cos )

+

a

b

( r 0

η

b

=

(4.35)

ζ x

4

2

sin )

+

a

b

( r 0

2

)

(4.36)

=ζ z

4

2

a )

( r 0 −

+

a

b

( r 0

Equation of water surface profile

Dealing with the flows at a vertical intake as in Figure 4.5, the depth-averaged equations

(3.26, 3.33 and 3.34) are applied. With the assumptions of steady vortex (

0=

) and

∂ ∂ t

axisymmetric flow, i.e. the flow is homogenous in η-direction (

), those

0=

∂ η ∂

equations become:

Continuity equation:

=

=

=

0

or

Q

const

.

(4.37)

0

d ξ d

M J

M J

0

0

Momentum equations:

ξ-direction:

ξ

2

2

)

+

+

=

( M

MN

NM

N

G

ξ +Γ ξξ 0

ξ Γ ξη 0

ξ +Γ ηξ 0

ξ Γ ηη 0

UM J

h J

ξ τ b ρ J

0

0

0

1 hJ 0

⎛ d ⎜ ⎜ ξ d ⎝

⎞ ⎟ ⎟ ⎠

2

+

+

ζ

2 x 0

({

)

+

+

+

(4.38)

UU

VU

UV

VV

G

ζ ξη

ζ ηη

Γ 0

Γ 0

ζ Γ ξξ 0

ζ Γ ξη 0

h 2

ξξξ 2 2 z y 0 0 J

d ξ d

0

⎫ ⎬ ⎭

η-direction:

η

2

2

)

+

+

=

G

( M

MN

NM

N

η Γ ηη 0

η +Γ ξξ 0

η Γ ξη 0

η +Γ ηξ 0

η τ b ρ J

h J

UN J

0

0

0

1 hJ 0

⎞ ⎟ ⎟ ⎠

⎛ d ⎜ ⎜ ξ d ⎝

41

2

+

+

ζ

0

0

0

0

0

x

x

y

z

)

+

+

+

({ UU

VU

UV

VV

G

ζ ξη

ζ ηη

Γ 0

Γ 0

ζ Γ ξη 0

ζ Γ ξξ 0

h 2

ηξηξηξ y z 0 J

d ξ d

0

⎫ ⎬ ⎭

(4.39)

Using the definition in Equation (3.32) in plugging with Equations (4.28 – 4.36), the

Christoffel symbols are hereafter calculated:

2

=

+

−=

x

y

z

(4.40)

ξ

ξ

ξ

ξ ξξ

Γ 0

4

ξ ∂ z ξ ∂

ξ ∂ x ξ ∂

ξ ∂ y ξ ∂

)

]2

+

dp 1 dr 0

1 2 p 1

a

b

b 2 [ )( ra 0

( rp 1 0

=

0

(4.41)

ξ ξη

ξ ηξ

Γ 0

Γ= 0

4

=

(4.42)

Γ ξ ηη 0

p 1

4

2

− )

) a +

( rr 0 0 − a

b

( r 0

=

0

(4.43)

η ξξ

η ηη

Γ 0

Γ= 0

−=

(4.44)

η ξη

η ηξ

Γ 0

Γ= 0

1 rp 01

−=

(4.45)

Γ ζ ξξ 0

2

4

2

)

+

a

a

b

p 1

( r 0

b 2 ) ( r 0

=

0

(4.46)

ζ ξη

ζ ηξ

Γ 0

Γ= 0

−=

(4.47)

Γ ζ ηη 0

4

2

br 0 ) −

+

a

b

( r 0

4

)

2

2

2

+

+

=

(4.48)

0

0

0

ξ x

ξ y

ξ z

4

2

− )

a +

( 2 rp 0 1 − a

b

( r 0

2

2

+

+

=

(4.49)

0

0

0

2 ηη y

x

η z

1 2

r 0

ξ and η introduced in Figure (4.5) are obviously perpendicular with each other, thus it

is given that

+

+

=

0

0

0

0

0

0

ξηξηξη y z 0

x

y

x

z

42

Therefore, Equation (4.39) becomes:

η

2

2

)

+

+

=

G

( M

MN

NM

N

(4.50)

η +Γ ξξ 0

η Γ ξη 0

η +Γ ηξ 0

η Γ ηη 0

h J

η τ b ρ J

UN J

0

1 hJ 0

0

0

⎛ d ⎜ ⎜ ξ d ⎝

⎞ ⎟ ⎟ ⎠

By substituting Equations (4.26, 4.43-4.44) and neglecting the effect of friction in

η-direction, Equation (4.50) is reduced as:

=

+

0

(4.51)

MV J

0

MN hJ 0

2 rp 01

⎞ ⎟ ⎟ ⎠

⎛ − ⎜ ⎜ ⎝

⎞ ⎟ ⎟ ⎠

⎛ d ⎜ ⎜ ξ d ⎝

Applying Equation (4.37) into (4.51) gives

=

Q

( ) + QV

0

0

0

d ξ d

N h

2 rp 01

⎞ ⎟⎟ ⎠

⎛ − ⎜⎜ ⎝

=

or

−= 2

(4.52)

dV ξ d

dV V

V 2 rp 01

dr 0 r 0

Integrating Equation (4.52) between two arbitrary points A and B with distance

Ar0 and

respectively from z -axis, the following relation is obtained:

0Br

B

B

(4.53)

−= 2

dV V

dr 0 r 0

V =⇒ B V A

A

A

2 r A 0 2 r B 0

are the contravariant components of velocity vector at A and B

where

AV and

BV

respectively.

Based on equation (4.53), the depth averaged tangential velocity component at

0r is

defined as:

B

B

A

B

A

2 r 0 A 2 r 0 B

= = = (4.54) ⇒ ⇒ r = eVv η Vr 0 rv 0 A rv 0 B

A

v r 0 v r 0

Eq. (4.54) is identical to the free vortex condition which allows us to assume the

vr02π=Γ

02π=Γ

43

for the whole flow regime, thus Vr 2 and conservation of circulation

from which:

= V (4.55) Γ π 2 02 r

The evaluation of the contravariant component of shear stress acting on the bottom is

2

approximated using the following relation:

= = = UefUf f

r V

r ξ

2

ξ τ b ρ

ξ τ b ρ

0

2 Q 0 hr 0

(4.56) ⇒ J

Using the relations described in Equations (4.40-4.49) and (4.45-4.46), Equation (4.38)

2

2

ξ

becomes

ξ ξξ

ρ

ξ τ b− J

0

0

0

2

2

2

2

2

2

)

+ = G Γ 0 M hJ h J M hJ 0 ⎞ ⎟ ⎟ ⎠ ∂ ⎛ ⎜ ⎜ ξ ∂ ⎝

ζ ξξ

( ξ x

ξ y

ξ z

ζ Γ ηη 0

0

0

0

ζ hG 2

0

− + + + − (4.57) Γ 0 ∂ ξ ∂ M 2 N 2 1 J ⎡ ⎢ ⎣ ⎤ ⎥ ⎦

2

ξ

2

2

2

2

0

)

Substituting (4.37) into (4.57) gives:

+

+

ξ ξξ

0

( ξ x

0

ξ y

0

ξ z

0

ρ

]ζ [ ∂ Γ ξξξ 0 ∂

Q 0 2 J

2 JQ 0 h

ξ τ b− J

0

0

0

0

2

2

2

2

2

2

2

)

)

]2

+ = Q G Γ 0 1 hJ h J ∂ ⎛ ⎜⎜ ξ ∂ ⎝ ⎞ ⎟⎟ ⎠

ζ ηη

( ξ x

0

ξ y

0

ξ z

0

( ξ x

ξ y

ξ z

0

0

0

[ ∂ ζ hG ξ ∂

0

4 r 0

− + + + + + Γ 0 ∂ ξ ∂ Γ π 2 8 J 1 J 02 ⎡ ⎢ ⎢ ⎣ ⎤ ⎥ ⎥ ⎦

Then, using the relations (4.56) and after some manipulation, we obtain the equation for

water surface profile in the form of

dh = ξ d

)

rhf ) ,( 0 1 ,( rhfp 0 21

(4.58)

44

in which

0.2

)

m

0.0

( z

Water surface Critical depth line

Quasi-normal depth line Bottom

-0.2

-0.1

0.0

0.1

0.2

x (m)

45

Figure 4.6 An example of computed water surface profile with quasi-normal depth line and critical depth line

4

2

3

2

]

[ (3

2

4

4

4

π 2 8

]

2 r 0 −

[ (

2 r 0

2

2

4

3

(2

a

)

(3

a

)

[ (

b

2

r 0

r 0

r 0

+

+

h

2 bQ 0

2 bQ 0

4

a ) 252

4

232

]

]

+

+

a

)

b

[ (

[ (

a

)

b

r 0

r 0

r 0

⎧ ⎪ ⎨ ⎪⎩

] ⎫ ⎪ ⎬ ⎪⎭

212

4

2

]

− + + − − Γ ( a b a ) r 0 r 0 = ) gb − b h rhf ,( 1 0 ) 232 r r (2 0 0 ] 232 − + + a ) b a ) [ ( a ) b r 0 r 0 ⎫ ⎪ ⎬ ⎪⎭ ⎧ ⎪ ⎨ ⎪⎩

[ (

3

2 0

2 0

2

2 r 0 −

π 2 4

− + Γ a r 0 + − − h - h Qf . (4.59) b 2 ) − ( a ) a ) ( 1 r 0 Q r 0 r 0 r 0 ⎧ ⎪ gb ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

2

2

and

3

3

2 0

4

212

4

212

π 2 4

]

]

2 r 0 −

[ (

[ (

− Γ ( a ) 1 r 0 = + − (4.60) ) gh bh Q rhf ,( 2 0 + − + a ) b a ) b r 0 r 0 r 0

Method of Calculation

For convenience of simulating the water surface profile, the transformation of special

0r is needed. Substitute Equation (4.11) into (4.58) it follows

coordinates from ξ to

that

= −= = . (4.61) dh ξ d dr 0 ξ d ) ) ) dh dr 0 1 p 1 dh dr 0 rhf ,( ) 0 1 ,( rhfp 0 21 dh −=⇒ dr 0 rhf ,( 0 1 ,( rhf 0 2

A common method of analysis including singular point often used in hydraulic

)

rhf ,( 0 1

engineering is applied. The singular point is the point at which both functions

in Equation (4.58) are equal to zero. The equation ,( ) 0 and ) rhf ,( 2 0 =rhf 0 1

,(

0)

=rhf 2 0

expresses critical depth line expresses the quasi-normal depth line whereas

as shown in Figure 4.6. The water surface profile is derived from Equation (4.58) by

using the fourth-order Runge-Kutta scheme with the initial slope near the singular point

46

defined by the following equation:

2

S

S

S

S

S

S

⎛ ⎜ ⎜ ⎝

⎞ ⎟ ⎟ ⎠

⎞ ⎟ ⎟ ⎠

⎛ ⎜ ⎜ ⎝

S

S

− + + − ± 4 ∂ f 1 ∂ h ∂ f 1 ∂ h ∂ f 2 ∂ h ∂ f 1 ∂ r 0 ∂ f 2 ∂ r 0 = (4.62) dh dr 0 2 ∂ f 2 ∂ r 0 ∂ f 2 ∂ h

(The subscript “s” denotes the derivatives at singular point)

4.3 Results and discussions

Posey and Hsu (1950) and Jain et al. (1978) have reported that a large reduction in intake

discharge was due to formation of vortices, especially in case of flow entering the intake

with entrapped air. The model herein derived has been applied with different imposed

circulations (i.e. different strengths of vortex) and the computed results show reduction of

discharge when the circulation increases while keeping the same submergence (water

head) (Figure 4.7). It is observed that the larger the circulation the larger the air core of

the vortex extends into the intake, and in case of air coming into intake the intake

− 310

× discharge decrease significantly (from 31.0 sm 3 , in case TC01 to

− 310

× 06.0 sm 3 , in case TC04 as in Table 4.1). This phenomenon is in agreement with

the previous experiments (Jain et al. 1978; Posey and Hsu 1950).

Intake discharge is inversely proportional to circulation as shown in Figure 4.8 where the

circle marks are the calculated results and solid line is the trend line calculated by the

least squares method.

The results also show that for small decrease in submergence, the air core becomes larger

(Figure 4.9). Similar results were reported by Anwar et al. (1978). In Figure 4.9, it is

noted that when the circulation increases, the submergence also would increase to

maintain the same intake discharge (see Table 4.2 for reference). Based on these results, it

is evident that the model is capable of representing the effects of circulation on the flow

47

through vertical intake.

0.15

0.10

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

0.05

+

+

)

+

+

+

+

m

+

+

( z

+ +

0.00

+

+

+

-0.05

+ TC 01 TC 02 TC 03 TC 04 Bottom

-0.05

0.05

0.10

0.00 x (m)

Figure 4.7 The effect of circulation on water surface profile and discharge at the

intake with same water head

Table 4.1 Parameters used in the calculations of results in figure 4.7

Case f a (m)

Γ (m2/s) 0.182 0.200 0.225 0.250

48

TC 01 TC 02 TC 03 TC 04 0.025 0.025 0.025 0.025 b (m2) 10-4 10-4 10-4 10-4 Q0 (10-3m3/s) 0.310 0.230 0.135 0.060 0.015 0.015 0.015 0.015

0.004

0.003

) s / 3 m

(

0.002

Q

0.001

0 0.3

0.4

0.5

0.6

Γ (m2/s)

0.2

0.1

)

m

(

z

0.0

TC 12 TC 13 TC 14 Bottom

-0.1

-0.1

0.1

0.2

0.0 x (m)

Figure 4.8 Variation of intake discharge with circulation (a=0.025m, b=10-5 m2, water head=0.5m)

Figure 4.9 Different water surface profiles with different values of circulation while maintaining the constant intake discharge

Table 4.2 Parameters used in the calculations of results in figure 4.9

f Case

Γ (m2/s) 0.30 0.25 0.20

49

TC 12 TC 13 TC 14 a (m) 0.02 0.02 0.02 b (m2) 10-4 10-4 10-4 Q0 (10-3m3/s) 0.15 0.15 0.15 0.015 0.015 0.015

' b in '

The effects of the shape of intake have been examined by changing the value of

− 4

− 5

2

equation (4.1). In this case, three simulations were considered with values of

510 m−

'a '

10

× 10 5 ;

=b

×

and while maintaining the radius of intake (i.e value of

− m210

2 . It can be seen from Figure 4.10 that the water depth in Equation 4.1) at

increases when the intake entrance becomes sharper (decreasing value of b as shown in

Table 4.3).

For more clarity, the computation has been done with the same intake and circulation to

test the effects of value of b on intake discharge without changing the submergence

'

(Figure 4.11a) and the effects on submergence with same discharge (Figure 4.11b). This

' b and it also

phenomenon is consistent with the physical meaning of the coefficient

improves the performance of bell-mouth intake.

So far we have been able to verify the model qualitatively based on the fact that no

previous studies are available which describe the variation in water surface profile at the

intake. However, a larger quantity of experimental data has been used by Odgaard (1986)

to verify the equation representing the critical submergence, which is defined as the

submergence when the tip of air-core vortex just reaches the intake (Figure 4.12), in the

5.0

=

[ N 6.5

] 5.0)( Γ

absence of surface tension (Eq. 18 in Odgaard, 1986) presented as follows:

(

Fr (4.63) S d R e 2.23

) 21gd

N

QS

Γ=Γ )(

0

, Froude number Fr = V and Reynold where circulation function

dQ υ=Re . number

This equation can be applied for the case of vortex with turbulent core, in which υ is

Γ+ kυ

×≈k 6

− 510

replaced by . From calibration the value of k was then estimated to be

50

.

0.20

0.10

)

m

(

z

0.00

TC 15 TC 16 TC 17 Bottom

-0.10

-0.1

0.0

0.1

0.2

Figure 4.10 Changing of water surface profile with different shape of the intake

Table 4.3 Parameters used in the calculations of results in figure 4.10

Case f

Γ (m2/s) 0.225 0.225 0.225

25

0.0020

20

0.0015

15

/

0.0010

) s / 3 m

d H

(

Q

10

0.0005

5

0.0000

0

0.00000

0.00012

0.00000

0.00004

0.00008

0.00012

0.00006 b (m2)

b (m2)

TC 15 TC 16 TC 17 a (m) 0.025 0.025 0.025 b (m2) 1.x10-4 5.x10-5 1.x10-5 Q0 (10-3m3/s) 0.3 0.3 0.3 0.015 0.015 0.015

Figure 4.11 The effects of b on discharge (17a) and submergence (17b) at an intake

51

S

z

d

x

102

/

) d S

( e c n e g r e m b u S

101

l

a c i t i r

'

C s d r a a g d O

100

100

102

101 Computed Critical Submergence (S/d)

Figure 4.12 Definition sketch of critical submergence

52

Figure 4.13 Comparison of computed critical submergence by the model (Equation 4.58) and by Odgaard’s equation (4.63)

For qualitative verification, the model is applied to predict the submergence at the intake,

2

=

and then comparisons are made with the Odgaard’s equation. A total of 12 simulations

b

− 510 m

01.0

m2.0

were made with and the range of intake radius between and .

Comparison of computed critical submergences and those computed using Odgaard’s

equation is shown in Figure 4.13. The horizontal coordinate is the computed critical

045 line. It can be seen that

submergences by the present model and vertical coordinate is by Odgaard’s. The perfect

fit between the model and Odgaard’s equation is on the

there exists a fairly good agreement. These comparisons show that the model can be used

12

to predict critical submergence of a vertical intake.

x +

10

b=0.000005

+

) d / S

8

b=0.00001

x

6

( e c n e g r e m b u S

x +

l

4

x +

x +

x +

a c i t i r

x +

+ x

C

+ x

+ x

+ x

2

0

0

0.04

0.08

0.12

0.16

0.2

0.24

d (m)

Figure 4.14 The variation of critical submergence wit different values of b

To verify whether the comparisons above are independent of the effect of intake shape, in

' b , the variation of critical submergence with different value of

'

this case the parameter

b while both the intake’s diameter and discharge are kept constant are shown in Figure

53

4.14. It is observed that the critical submergence is not significantly affected by a small

2

510 m−

'

change of b with order of . Based on this observation we strongly believe that

' b on the critical submergence is negligible when

the effect of ' b is small enough. '

4.4 Summary

In this application, a theoretical model based on the depth-averaged and momentum

equations has been developed to calculate the variation of water depth of flow with

air-core vortex at an intake. Hence, the water surface profile can be directly derived from

the model. The use of curvilinear coordinate avoids invalidation of the Kelvin’s theorem

in the core of vortex, therefore the conservation of circulation can be applied for whole

flow.

The comparisons of calculated data using the model and an empirical equation show that

the proposed model yields reliable results in predicting the critical submergence of the

intake without any limitation of Froude number, a problem that most of existing models

cannot escape.

4.5 References

1. Anwar O. H., Weller A. J. and Amphlett B. M. 1978. Similarity of free-vortex at

horizontal intake. J. Hydr. Res., Vol. 16 (2). pp. 95-105.

2. Hite J. E. and Mih C. W. 1994. Velocity of air core vortices at hydraulic intakes. J.

Hydr. Eng. ASCE, Vol. 120(3), pp.284-297.

3. Hosoda T. 2003. Depth averaged model of open channel flows over an arbitrary

surface. Proceeding of the International symposium on Shallow flows, Part II, June

16-18, pp. 267-271.

4. Jain A. K., Raju K. G. R. and Garde R. J. 1978. Vortex formation at vertical pipe

intake. J. Hydr. Div. ASCE, Vol. 104 (10), pp.1429-1445.

54

5. Jain S. C. 1984. Tangential vortex-inlet. J. Hydr. Eng. ASCE, Vol. 110(12),

pp.1693-1699.

6. Odgaard A. J. 1986. Free-surface air core vortex. J. Hydr. Eng. ASCE, Vol. 112(7),

7. Posey C. J. and Hsu H. 1950. How the vortex affects orifice discharge. Engineering

pp.610-620.

News Record, Vol. 144, p. 30.

8. Steffler M. P. and Jin Y. C. 1993. Depth averaged and moment equations for

moderately shallow free surface flow. J. Hydr. Res., Vol. 31(1), pp. 5-17.

9. Townson J. M. 1991. Free-surface hydraulics. London Unwin Hyman.

10. Trivellato F., Bertolazzi abd Firmani B. 1999. Finite volume modeling of free

surface draining vortices. J. Comp. App. Mech., Vol. 103, pp. 175-185.

11. Yıldırım N. and Kocabas F. 1995. Critical submergence for intakes in open channel

flow. J. Hydr. Eng., ASCE, Vol. 121(12), pp 900–905.

12. Yıldırım N. and Kocabas F. 1997. Closure to ‘Critical submergence for intakes in

open channel flow. J. Hydr. Eng., ASCE, Vol. 123(6), pp 589–590.

13. Yıldırım N. and Kocabas F. 1998. Critical submergence for intakes in still-water

55

reservoir. J. Hydr. Eng., ASCE, Vol. 124(1), pp. 103–104.

Chapter 5

UNSTEADY PLANE-2D ANALYSIS OF FLOWS WITH AIR-CORE VORTEX

In most of the studies in the past, because of the complexity of 2D equation the variation

of the flow quantities in the tangential direction at the same radius from the center of the

vortex was neglected with the assumption of homogeneity of flow in circular direction. In

this study, with the advantages of the 2D depth-averaged model, it is possible to estimate

the variation in 2D flow.

5.1 Governing equations

Coordinate setting:

In order to apply the depth-averaged equations (3.26, 3.33 & 3.34) to simulate the 2D

)ζ perpendicular with that

flows with air-core vortex, a new coordinate system is introduced as in Figure 5.1. Two

)

coordinates ,( ηξ are set on the bottom and the other ( )

0

1

) 1R is

56

surface. Parameters ( a b R R define the geometry of the bottom in which ( , , ,

)0R refers to the curvature of the edge of the intake entrance,

radius of the intake, (

while a and b are radii of the flow approach area and the length of intake,

respectively.

ζ

η

z

0R

a

ξ

b y

1R

x

Figure 5.1 Definition sketch of the new coordinates

Governing equations:

The original depth-averaged equations are used:

+

+

=

0

Continuity equation:

∂ ∂ ξ

∂ ∂ η

1 J

∂ h ∂ t

M J

N J

0

0

0

(5.1)

ξ-direction

2

2

+

+

+

+

+

+

( M

Momentum equations:

ξ ξξ

ξ ξη

ξ ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ξ ∂

∂ ∂ t

0

0

0

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

⎛ ⎜⎜ ⎝

⎛ ⎜⎜ ⎝

∂ ⎛ ⎜⎜ η ∂ ⎝

h

h

+

+

+

+

ξ

0

0

0

0

0

2 0 x

x

x

y

z

=

ζ

MN NM N M J UM J VM J 1 hJ 0

ξ τ b ρ J

ηξηξηξ 0 y z J

ξξξ 2 2 0 0 y z J

0

0

0

0

∂ p ∫ d ρξ ∂ 0

∂ p ζ ∫ d ρη ∂ 0

57

(5.2) G h J

η-direction

2

2

+

+

+

+

+

+

( M

η ξξ

η ξη

η ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ξ ∂

∂ ∂ t

0

0

0

⎞ ⎟ ⎟ ⎠

⎞ ⎟ ⎟ ⎠

⎞ ⎟ ⎟ ⎠

⎛ ⎜ ⎜ ⎝

⎛ ⎜ ⎜ ⎝

∂ ⎛ ⎜ ⎜ η ∂ ⎝

h

h

+

+

+

+

η

0

0

0

0

0

2 0 x

x

x

y

z

=

ζ

MN NM N N J UN J VN J 1 hJ 0

ζ d

ηηη 2 2 0 0 y z J

ηξηξηξ 0 y z J

η τ b ρ J

0

0

0

0

∂ p ∫ ρη ∂ 0

∂ p ∫ ρξ ∂ 0

G d h J

(5.3)

In this case, it can easily be seen that the coordinates setting on the bottom surface with

ξ 0=Γ=Γ ξη

ξ ηξ

η ηη

η 0=Γ=Γ ξξ

ζ 0=Γ=Γ ξη

ζ ηξ

+

+

=

0

ηξηξηξ 0 y z

0

0

0

0

y

x

x

z

0

and

0=ηG

As a result, Equations (5.2-5.3) are reduced to

2

2

+

+

+

+

( M

two axis ) ,( ηξ are orthogonal system, therefore it was proved that:

ξ ξξ

ηη

Γ 0

Γ 0

∂ ∂ t

0

0

0

1 hJ 0

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

⎛ ⎜⎜ ⎝

∂ ⎛ ⎜⎜ ξ ∂ ⎝

∂ ⎛ ⎜⎜ η ∂ ⎝

2

2

2

+

+

ξ

2 x 0

=

(5.4)

N M J UM J VM J

ζ −Γ µη

ζ +Γ ξξ

∂ ξ ∂

ζ hG 2

ξξξ 2 2 y z 0 0 J

ξ τ b ρ J

0

0

0

⎡ ⎢ ⎣

⎤ −⎥ ⎦

+

+

+

+

( MN

G M 2 N 2 h J

η ξη

ηξ

Γ 0

Γ 0

∂ ∂ t

0

0

0

1 hJ 0

⎛ ⎜⎜ ⎝

∂ ⎛ ⎜⎜ ξ ∂ ⎝

∂ ⎛ ⎜⎜ η ∂ ⎝

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

2

2

2

+

+

2 x 0

−=

(5.5)

ζ −Γ µη

ζ +Γ ξξ

NM N J UN J VN J

ζ hG 2

ηηη 2 2 z y 0 0 J

η τ b ρ J

0

0

∂ ⎡ ⎢ η ∂ ⎣

⎤ −⎥ ⎦

These equations together with continuity Equation (5.1) are then applied to calculate 2D

58

M 2 N 2

water surface profile of flow with air-core vortex at a vertical intake in the following

sections.

j=3

j=2

j=1 i=1

i=2

Figure 5.2 Illustration of the computational grid

5.2 Numerical method

Based on the generalized curvilinear coordinate system set on the bottom surface (Figure

5.1), the computation mesh is generated as shown in Figure 5.2. Then the set of equations

(5.1, 5.4-5.5) is discretized by using the cell-centered staggered grid (Prayet and Taylor,

1983), in which the hydraulic variables are defined at different positions as shown in

Figure 5.3 (Fletcher, 1991).

The Finite Difference Method and the Finite Volume Method are used in dicretization

process, in which the First-Upwind scheme is employed for convective terms while the

Two-point Forward scheme is utilized for time derivatives. It follows that:

Continuity equation:

M

M

+

+

+

+

+

+ n 1 h + i ,21

j

21

n h + i

,21

j

21

,1

i

j

21

ji ,

21

1

+

1 ∆ ξ

J

∆ t

J

J

+

+

+

+ ,210

i

j

21

+ ,10 i

j

21

,0 ji

21

⎡ ⎢ ⎢ ⎣

⎤ ⎥ ⎥ ⎦

59

N, V

j+1

M, U

M, U

h

j

i

N, V

i+1

Figure 5.3 Definition sketch of cell-centered staggered grid in the 2D calculation

N

N

+

+

,21

+ 1

,21

i

j

i

j

+

=

0

(5.6)

J

J

+ ,210

+ 1

i

j

+ ,210

i

j

⎤ ⎥ ⎥ ⎦

⎡ 1 ⎢ η ∆ ⎢ ⎣

Momentum equations:

ξ-direction

M

M

U

M

M

U

+

+

+

+

+

+

21

21

,

21

21

+− 1

21

+ 1 n + , 21 ji

n , ji

i

j

+ jai

i

, jb

i

1

+

1 ξ ∆

J

∆ t

,21 J

− ,21 j J

+

+

+

21

+ ,210

21

,0 ji

i

j

− ,210

i

j

21

⎤ ⎥ ⎥ ⎦

⎡ ⎢ ⎢ ⎣

+

21

− ;21

, ji

c

MV , ji

, ji

d

1

]

+

+

+

=

[ 2 hU

2 hV

ξ ξξ

ξ ηη

Γ 0

Γ 0

+

ji .

21

MV + 1 , ji J

J

J

+

+ 1

21

,0 ji

,0 ji

.0 ji

⎡ 1 ⎢ ∆ η ⎢ ⎣

⎤ ⎥ ⎥ ⎦

+

+

+

h ji .

21

2 x 0

[

=

Γ

G

M

+

+

2 + i

.21

j

21

+

+

ξ . ji

21

ζ ξξ i

.21

j

21

ξ τ b ρ

ξ

1 ∆

J

ξξξ 2 2 z 0 y 0 J

2

+

+

0

.0 ji

21

ji .

21

+

ji .

21

Γ

+

Γ

M

N

N

+

+

+

21

.21

21

21

2 − i .21

j

2 + i

j

2 − i .21

j

+

+

+

+

.21

21

− .21

21

ζ ηη i

j

ζ Γ ηη i

j

ζ ξξ i

j

− .21

21

]

+

G

G

(5.7)

+

+

+

+

.21

21

21

j

2 h + i

2 h − i .21

j

.21

21

21

ζ + i

j

ζ − i .21

j

60

N, V

N, V

j+1

M, U

h

h

N, V

N, V j

i-1/2

i

i+1/

M, U

M, U

h

j+1/

N, V

j

M, U

M, U

j-1/2

h

i

i+1

Figure 5.4 Illustration of the discretization scheme for the momentum equations

61

η-direction

N

N

U

N

+

+

+

+

,21

21

− 21

+ n 1 + i ,21

j

n + i

j

i

i

ja ,

NU ji ,

i

jb ,

1

+

1 ∆ ξ

J

∆ t

j ,1 J

J

+ ,210

i

j

+ ,10 i

j

,0 ji

⎡ ⎢ ⎢ ⎣

⎤ ⎥ ⎥ ⎦

V

N

N

V

+

+

+

+

+

21

,21

− 21

,21

+− 1

i

j

i

+ cj

i

j

d

i

]

+

+

[ Γ+Γ

=

η ξη

η ηξ

+

,21

i

j

,21 J

j ,21 J

UVh J

+

0

+ ,210

21

+ ,210

− 21

i

j

i

j

⎤ ⎥ ⎥ ⎦

⎡ 1 ⎢ η ∆ ⎢ ⎣

+

+

0

2 x 0

[

−=

Γ

Γ

M

M

+

,21

21

,21

− 21

2 + i

j

2 + i

j

+

+

+

,21

21

,21

− 21

ζ ξξ i

j

ζ ξξ i

j

τ b η ρ

1 ∆ η

ηηη 2 2 y z 0 J

2

+

+

0

,21

i

j

,21

i

j

+

Γ

Γ

N

N

G

+

+

,21

21

,21

− 21

2 + i

j

2 + i

j

+

+

+

+

+

,21

21

ζ + i

j

,21

21

,21

− 21

,21

21

ζ ηη i

j

ζ ηη i

j

2 h ηη i

j

]

+

G

(5.8)

+

,21

− 21

ζ + i

j

,21

− 21

2 h ηη i

j

To close the above equation system, the boundary conditions were given with the

discharge at the upstream end and Neumann condition at the downstream boundary. The

= JY

j

computational domain is chosen from

1=j

to

1+

, therefore the hydraulic

= JY

j

variables at

1+

is used for boundary condition when calculating variables at

1=j

and vice versa.

5.3 Results and discussions

)

)

With different discharges (

0V given at the outer-zone,

0M and tangential velocities (

some results of the surface profile were obtained as shown in Figures (5.5-5.7).

Similarly as the phenomenon monitored in Chapter 4, the water surface and submergence

of the flow with air-core vortex strongly depend on the flow condition such as discharge

at the intake, the circulation (which in this simulation is defined from the velocity at

boundary of outer zone) and also the geometry of the intake.

When the discharge at the intake is increased, the submergence is also increased as shown

in Figure 5.5 with the same circulation. If the discharge at the intake is maintained, when

62

a)

b)

0.5

0.5

0.4

0.4

0.3

0.3

)

)

m

m

( z

( z

0.2

0.2

0.1

0.1

30

30

0.15

0.15

0.00

0.00

0.0

0.0

-0.15

-0.15

y(m)

y(m)

0.3

0.3

0.2

0.2

-0.30

-0.30

-0.1

-0.1

-0.3

-0.3

0.0 x (m)

0.0 x (m)

Figure 5.5

Water surface of flow with different discharges at the intake

=

=

=

=

QM

0.35;

10.;

0.03 ;

0.05 ; m

a)

0

V 0

R 0

m R 1

=

=

=

=

QM

0.5;

10.;

0.03 ;

0.05 ; m

b)

0

V 0

R 0

m R 1

a)

b)

0.5

0.5

0.4

0.4

0.3

0.3

)

)

m

m

( z

( z

0.2

0.2

0.1

0.1

30

30

0.15

0.15

0.00

0.00

0.0

0.0

-0.15

-0.15

y(m)

y(m)

0.3

0.3

0.2

0.2

-0.30

-0.30

-0.1

-0.1

-0.3

-0.3

0.0 x (m)

0.0 x (m)

Figure 5.6

Water surface of flow with different velocity at the outer-zone boundary

=

=

=

=

QM

10.;

0.5;

0.05 ;

0.05 ; m

a)

V 0

R 0

m R 1

0

=

=

=

=

QM

15.;

0.5;

0.05 ;

0.05 ; m

b)

0

V 0

R 0

m R 1

63

velocity at boundary varies the water surface also varies. The higher circular velocity

)

( 0V , the higher submergence at the intake, or in other words, the circular velocity

obstruct or impede the water drainage at the intake. Consequently to maintain the same

discharge, water head would be increased with increased circular velocity (Figure 5.6).

In order to investigate the effect of the intake geometry to water surface, the parameter

0R is used to control the shape of the join between intake and bottom. Figure (5.7) shows

that the submergence increases when

0R increases while keeping other parameters

constant.

a)

b)

0.5

0.5

0.4

0.4

0.3

0.3

)

)

m

m

( z

( z

0.2

0.2

0.1

0.1

30

30

0.15

0.15

0.00

0.00

0.0

0.0

-0.15

-0.15

y(m)

y(m)

0.3

0.3

0.2

0.2

-0.30

-0.30

-0.1

-0.1

-0.3

-0.3

0.0 x (m)

0.0 x (m)

Figure 5.7

Water surface of flow with different shape of the intake

=

=

=

=

a)

QM

10.;

0.5;

0.05 ;

0.05 ; m

V 0

R 0

m R 1

0

=

=

=

=

QM

10.;

0.5;

0.03 ;

0.05 ; m

b)

0

V 0

R 0

m R 1

5.4 Summary

Based on the simulation results, it can be concluded that, not only 1D water profile as

described in Chapter 4 can be estimated by using analytical method, but also the 2D

unsteady plane water surface can be reasonably simulated using the depth-averaged

equations without any difficulty. Consequently, it also validates the broad perspective of

the model’s applicability.

64

5.5 References

2 Peyret R. and Taylor T. D. 1983. Computational methods for fluid flow, Springer Ser.

Comput. Phys. (Springer, Berlin, Heidelberg).

3 Fletcher C. A. J. 1991. Computational techniques for fluid dynamics. Springer Ser.

Comput. Phys. (Springer Verlag).

65

Chapter 6

WATER SURFACE PROFILE ANALYSIS OF FLOWS OVER CIRCULAR SURFACE

6.1 Preliminary

In a chute spillway, the stability of the free surface is strongly dependent on the geometry

of the approach channel as well as that of spillway. If there is an abrupt drop of the

channel bed, a free overfall or a vertical fully aerated drop may occur. This phenomenon

has attracted considerable attentions of investigators (e.g., Davis et al. 1998; Dey 1998,

2002, 2003, 2005; Ahmad 2005; Guo 2005; Pal and Goel 2006) since the pioneering work

of Rouse (1936) due to its practical engineering use in simple flow-measuring devices

(e.g., Dey 1998; Guo 2005). But their approaches could not be applied in case of an

abrupt change of channel’s bed with circular shape, accompanied with small discharge

that cannot result in free overfall.

In this chapter a mathematical model based on depth-averaged equations with a

66

generalized curvilinear coordinate system attached to the bottom surface is developed

from the 2D depth-averaged model described in Chapter 3. The developed model was

applied to calculate the water surface profile in open channel with circular surface at the

end. To verify the model, an experiment was conducted to measure the water surface

profile, where the model results of both steady and unsteady analysis were compared with

the observations.

6.2 Hydraulic Experiment

6.2.1 Experimental Setup

The experiments were conducted in an open-channel main flume consisted of two straight

channels joined perpendicular to each other by a step ‘A’ at the end as described in Figure

6.1. The flume was 135cm long, 10cm wide and 20cm deep (Figure 6.1a) , made of a

steel frame with glass in all side walls and bed which allowed the convenience for visual

observations.

In this study two detachable circular parts with radius of 15cm and 5cm respectively were

used (Figure 6.1b). These parts were made such that, it is easily to be attached and

detached to and from the main flume at step ‘A’, without altering the other parts of

experiment facility to form a circular channel bottom at the joint (Figure 6.1c). These

parts (B and C) were purposely constructed to investigate the effect of different

curvatures to flows. For simplicity of discussion, from now onward, the curvature with

15cm and 5cm radius will be referred to “large case” and “small case” respectively.

The water surface elevations were measured by the ultrasonic sensor Keyence UD-500

accompanied with the receiver and converter Keyence NR-2000 (Figure 6.2). For each

case the measurements were focused in the region with the curvature of the bottom at four

points 0o (point 1), 30o (point 2), 60o (point 3) and 90o (point 4) to the vertical direction

(Figure 6.1d).

67

135

20

A

20

20

75

100

c) Connection of main flume and parts B, C

20

Inflow

a) Main flume

20

Point 1 Point 2

Outflow

Point 3

R = 15

R = 5

Point 4

d) Sensors’ arrangement

b) Part B and part C

Figure 6.1 Side view of the experimental facility

(All units are in centimeter)

Sensor UD-500

Circular surface

Digitizer NR-2000

Figure 6.2 Experimental site

68

The ultrasonic sensor Keyence UD-500 measures the distance from the sensor to the

surface by calculating the reflected ultrasonic signal. Actually, the sensor is combined

both the transmitter and receiver of the ultrasonic signals, and the output is voltage of the

direct electric current, therefore it needs a receiver and digitizer to convert the analog

signal into digital signal (Keyence NR-2000) for out-putting to the PC. Thus, the sensor is

connected with the receiver and converter as in Figure 6.3.

6.2.2 Experimental procedures and results:

Before using the sensor and converter, it is required to calibrate the sensor, and establish a

correlation line between the voltage and the distance from the sensor to the surface. The

calibration process was done by using the sensor to measure the known distances (Figure

6.4) and the correlation equation was found:

=

(6.1)

35.52

*77.40

D t

V t

where:

tD is the distance (in mm ) and

tV is the measured voltage (in Voltage V ) at

time t .

)

At the first stage of experiment, the distance from the fixed sensor to the bed was measured in advance (

bD and then the time series of distance from the sensor to the

)

water surface were saved (

sD , thereafter water depth was easily obtained by

subtracting

sD from

bD .

Water was drawn from a constant head tank being filled by a turbine pump to ensure the

steady in-flow rate. The flow rate was controlled by a control valve and discharge

measurements were carried out at the end of the flume for different flow conditions.

At beginning, the discharge was set at very small value, then increased with small

intervals. After each increase, the flows were kept for some time to get the equilibrium

condition before taking measurements. Upon attainment of steady conditions at the

69

upstream end of the main flume, the water surface level in the region within circular

bottom was measured by the sensor at four locations as in Figure 6.1d. The samples of

time series results of water depth are shown in Figure 6.5. And the time-averaged values

with the flow conditions are listed in Table 1.

At the first stage, with relatively small discharge, the flows over the circular part were

observed with very slightly fluctuation and considered to be stable (Figure 6.5a &6.5b).

In the second stage, when the discharge was continuously increased, the oscillation

occurs at the circular region and become significant with increasing flow rate while

remaining other parameters as the same. With each subsequent increase of the discharge,

the amplitudes of oscillation become larger (Figure 6.5c & 6.5d). The intensities of

oscillation can be estimated by Equation (6.2) and the results are shown in Figure 6.6.

2

n

h

h i

1σ =

(6.2)

∑ n 1

h

⎛ ⎜ ⎜ ⎝

⎞ ⎟ ⎟ ⎠

From this figure, it is observed that, in both small- and large-case, the bigger intensities

were observed in downstream end of the flow field, or in other words, the oscillation

intensity increases from point 1 to point 4.

70

Transmitter &

Sensor UD-500

Receiver

Digitizer NR - 2000

PC

Figure 6.3 Schematic of sensor connection

Figure 6.4 Sensor calibration

71

a)

15

)

m m

Point 1 Point 2 Point 3 Point 4

10

( h t p e D

5

0

0

1

2

3

4

5

Time (s)

b)

15

)

m m

10

( h t p e D

5

0

0

1

2

3

4

5

Time (s)

c)

15

)

m m

10

( h t p e D

5

0

0

1

2

3

4

5

Time (s)

d)

15

)

m m

10

( h

t p e D

5

0

0

1

2

3

4

5

Time (s)

Figure 6.5 Time history of the free surface at four locations in different experiments:

a) Exp1; b) Exp2; c) Exp3; d) Exp4;

72

Table 6.1 Experiment conditions

Time averaged water depth (mm)

Discharge

Experiment

Radius of circular part

(l/s)

Point1

Point2

Point3

Point4

Small-case

Exp 1

5 cm

0.1704

6.2

3.8

2.2

1.4

Small-case

Exp 2

5 cm

0.2441

8.6

4.8

3.1

1.7

Small-case

Exp 3

5 cm

0.6350

11.1

7.8

5.0

3.5

Small-case

Exp 4

5 cm

0.7754

14.4

10.5

7.8

5.6

Large-case

Exp 5

15 cm

0.4330

10.8

5.1

2.2

1.8

Large-case

Exp 6

15 cm

0.5421

12.7

5.7

3.9

3.1

Large-case

Exp 7

15 cm

0.8295

15.8

8.6

5.1

3.8

Large-case

Exp 8

15 cm

0.9936

18.6

10.8

7.3

5.1

0.30

Exp-1

Exp-2

0.25

Exp-3

Exp-4

0.20

Exp-5

Exp-6

σ

0.15

Exp-7

Exp-8

0.10

0.05

0.00

Point1

Point2

Point3

Point4

Figure 6.6 The oscillation density at four locations in circular region

73

6.3 Steady analysis of water surface profile

Firstly, the new generalized curvilinear coordinate is introduced. The coordinate system is

)ηξ−

generated based on the bottom surface with two axes (

attached to the surface and

the other axis (

)ζ is a straight line perpendicular to it (Figure 6.7). In Figure 6.7, η is

)

xOz or is identical with y -axis. The surface is

the straight axis normal to plane (

=zyxζ

expressed mathematically by the equation

0),

,(

.

The depth-averaged equations in this new coordinate system are used (Equations 3.26,

3.33 and 3.34):

Continuity equation:

+

+

=

0

(6.3)

∂ ξ ∂

∂ η ∂

1 J

∂ h ∂ t

M J

N J

0

0

0

Momentum equation in ξ-direction:

2

2

+

+

+

+

+

+

( M

MN

NM

N

ξ ξξ

ξ ξη

ξ ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ∂ t

M J

UM J

VM J

0

0

0

1 hJ 0

⎛ ⎜⎜ ⎝

∂ ⎛ ⎜⎜ ξ ∂ ⎝

∂ ⎛ ⎜⎜ η ∂ ⎝

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

⎞ ⎟⎟ ⎠

h

h

+

+

+

+

ξ

0

0

0

0

0

x

x

y

z

2 x 0

=

ζ

G

h J

ξξξ 2 2 z y 0 0 J

ξ τ b ρ J

ηξηξηξ 0 y z J

0

0

0

0

∂ p ζ ∫ d ∂ ρξ 0

∂ p ∫ d ∂ ρη 0

(6.4)

Momentum equation in η-direction:

2

2

+

+

+

+

+

+

( M

MN

NM

N

η ξξ

η ξη

η ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

VN J

∂ ∂ t

N J

UN J

0

0

0

1 hJ 0

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ η ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ ξ ⎝

⎞ ⎟⎟ ⎠

⎛ ⎜⎜ ⎝

h

h

+

+

+

+

η

x

x

y

z

0

0

0

0

0

2 x 0

ζ

ζ

=

G

h J

ηξηξηξ y z 0 J

η τ b ρ J

ηηη 2 2 z y 0 0 J

0

0

0

0

∂ p ∫ d ρη ∂ 0

∂ p ∫ d ρξ ∂ 0

(6.5)

74

ζ

ζer

ηer

ξer

W

COVARIANT BASE VECTORS

U

z

x

O

ξ

Figure 6.7 Curvilinear coordinates attached to the bottom

75

Neglecting transversal transport, then it follows that:

+

=

0

(6.6)

∂ ξ ∂

1 J

∂ h ∂ t

M J

0

0

2

2

ξ

2

ξξ + 2 2 x z 0 0

+

+

=

M

G

ξ ξξ

ζ ξξ

Γ 0

Γ 0

∂ ξ ∂

∂ ∂ t

M J

UM J

h J

J

ξ τ b ρ J

M 2

ζ hG 2

0

0

0

0

1 hJ 0

0

⎛ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ξ ∂ ⎝

⎞ ⎟⎟ ⎠

⎡ ⎢ ⎣

⎤ −⎥ ⎦

(6.7)

Using the system of equations (6.6-6.7), we can analyze the water surface profile of the

flow over a circular surface as depicted in Figure 6.7. Consider the steady state, equations

(6.6-6.7) can be recast as follows:

Continuity equation:

=⇔=

JQM

Q

=⇔=

0

const.

(6.8)

0

0

0

d ξ d

M J

M J

0

0

0

=

Uh

=⇔ U

or

(6.9)

JQ 0

0

JQ 0 h

Momentum equation:

2

2

ξ

2

+ ξξ 2 2 z x 0 0

+

=

(6.10)

G

M

ξ ξξ

Γ 0

ζ Γ ξξ 0

h J

d ξ d

UM J

J

ξ τ b ρ J

M 2

ζ hG 2

0

1 hJ 0

0

0

0

⎛ ⎜⎜ ⎝

⎞ ⎟⎟ ⎠

⎡ d ⎢ ξ d ⎣

⎤ −⎥ ⎦

Substitute Equations (6.8), (6.9) into (6.10), after some manipulations, the equation for

water surface profile can be obtained in the form of Equation (6.11):

−=

(6.11)

dh ξ d

ξ ) ξ )

,( hf 1 ,( hf 2

with

ζ

d

2

4

3

+ ξξ 2 2 x z 0 0

+ ξξ 2 2 x z 0 0

=

+

h

h

ξ hG

J

ξ )

ζ ξξ

2 JQ 0

0

0

Γ+ 2 0

,( hf 1

0 ξ

ξ

dJ d

ζ Γ ξξ 0 ξ d

dG d

2

2

⎡ ⎢ ⎢ ⎣

⎤ ⎥ ⎥ ⎦

+

J

− fQh

J

(6.12)

ξ ξξ

2 JQ 0

0

Γ 0

0

2 0

2 0

0 ξ

dJ d

⎡ ⎢ ⎣

⎤ ⎥ ⎦

76

and

(

=

) + 3 JQhG

)

(6.13)

,( hf 2

2 0

2 0

ζξξξ + 2 0 x

2 0 z

The common method of analysis including singular point analysis similar as in Chapter 4

is applied to calculate the water surface profile in both the upstream and downstream

direction from this point. The singular point is defined as the point at which both

functions

and

in Equation (6.11) are equal zero. The definition of

ξhf ) ,(1

ξhf ,(2 )

singular and an example of critical and quasi-normal depth lines can be found in Figure

6.8.

Figures 6.9 to 6.12 show the comparison of water surface profile between the 1D model

results under steady state condition and the experiment for Exp-01, Exp-02, Exp-05 and

Exp-06. It is observed from the figures that the calculated results are consistent with the

experimental data in both large and small cases. Note that during the experiment, because

of limitation of sensors the measurements were only focused in the vicinity of the circular

surface.

77

0.55

)

m

(

Water surface

z

Quasi-normal depth line

Critical line

Bottom

0.50

0.95

1.00

1.05

1.10

x (m)

Figure 6.8 Illustration of computed water surface profile with quasi-normal and critical depth lines

78

0.60

0.55

)

m

(

Exp.

z

0.50

Cal.

Bottom

0.45

0.90

0.95

1.05

1.10

1.00 x (m)

Figure 6.9 Steady water surface profile with conditions of Exp1

0.60

0.55

)

m

(

Exp.

z

0.50

Cal.

Bottom

0.45

0.90

0.95

1.05

1.10

1.00 x (m)

Figure 6.10 Steady water surface profile with conditions of Exp2

79

0.70

0.60

)

m

(

Exp.

z

Cal.

Bottom

0.50

0.90

1.00

1.10

1.20

x (m)

Figure 6.11 Steady water surface profile with conditions of Exp5

0.70

0.60

)

m

(

Exp.

z

Cal.

Bottom

0.50

0.90

1.00

1.10

1.20

x (m)

Figure 6.12 Steady water surface profile with conditions of Exp6

80

6.4 Unsteady characteristics of the flows

In order to study the unsteady characteristics of flows, the numerical method for unsteady

analysis is developed by applying the common method used in the field of the numerical

hydraulics.

:

UM ,

: h

i-1 i-1/2 i i+1/2 i+1

Figure 6.13 Illustration of staggered grid

Equations (6.6-6.7) and the same curvilinear coordinate with previous steady analysis is

also used, and for solving those equations numerically, the two-point forward scheme is

applied for time integration, on the other hand the First Upwind scheme is employed for

the convective term based on cell-centered staggered grid (Figure 6.13) for spatial

discretization as in the following Equations (6.14-6.15).

Continuity equation:

h

21

+ n 1 h + i 21

n + i

1

n i

n + i 1

+

=

(6.14)

0

1 ∆ ξ

J

∆ t

M J

M J

+ 21 0

i 0

+ i 1 0

i

⎡ ⎢ ⎢ ⎣

⎤ ⎥ ⎥ ⎦

Momentum equation:

+ 1

M

M

+− 1

i

+ ai

b

i

i

n i

n i

+

Γ 0

+ 2 hU i i

= Gh i

ξ i

ξ ξξ i

J ∆

i 0 ξ

− ∆ t

MU + 21 J

MU − 21 J

i 0

1-i 0

⎡ ⎢ ⎣

⎤ ⎥ ⎦

i

i

) [

]

) [

]

+

M

M

21

21

Γ 0

2 i

− Γ 2 i 0 1

ζ 2 hG + i i

ζ 2 hG − − i i 1

− 1

ζ ξξ i

ζ ξξ i

ξ τ b ρ

( ξξ + 2 2 x z 0 0 ξ 2

( + ξξ 2 2 x z 0 0 ξ 2

⎛ ⎜ ⎜ ⎝

⎞ ⎟ ⎟ ⎠

i

(6.15)

81

where:

+

i

21

− 21

i

+

i

21

i

21

≥ ≥ U U 0 = = a 0 b ; < < U 0 U 0 ⎧ if 0 ⎪ ⎨ if 1 ⎪⎩ ⎧ if 0 ⎪ ⎨ if 1 ⎪⎩

n + 1 i

n i

n + i

21

+ U U = and U 2

r

The bottom shear stress was evaluated using the following equation:

r ξ

ξ τ b ρ

= = 2UefUUf (6.16)

where f is the resistant coefficient. In this study, since all walls were made of smooth

6.7.

The model has been applied in estimating the water surface profile of flows over a

circular surface (Figure 6.7), based on the flow conditions identical with that of the

experiments (Table 6.1). The inflow discharge is given at the upstream and the Neumann

derivatives are provided at the downstream end to complete the boundary conditions.

In order to determine the unsteady characteristics of the flow over circular surface, the

same procedures as in experiment were applied for simulation. Firstly the small constant

discharge was used as upstream boundary (Exp -1, Exp -2, Exp – 5 and Exp - 6). With

constant initial depth, the equilibrium state of flow was reached after around 20 seconds.

glass, the value of f = 0.01 was used. ξer is covariant base vector as shown in Figure

Then, comparisons were made between the numerical and experimental water surface

profile as shown in Figure 6.14 – 6.15 for “small case” and Figure 6.16 – 6.17 for “large

case”. Note that, with these relative small discharges the water surface profiles become

stable as can be seen from Figures 6.14 – 6.17, and the good agreements between

the model in simulating the flows over a circular structure. These agreements are similar

82

calculated and observed results can be observed from the figures implying the ability of

with the steady analysis in section 6.3.

In the second stage, the slightly increased upstream discharge were applied. Figure

6.18-6.21 shows water surface profile with different discharge in which the bold line is

the bottom plane and the faint lines are the water surface. Note that the water surface

profiles shown in the same figure are based on different time. It is observed that with

increase of discharge the water surface oscillation started to appear (Figure 6.18, 6.20).

Further increase of discharge resulted in increase of amplitude of variation (Fig. 6.19,

6.21). Moreover, in Exp - 3 and Exp - 4, the amplitude of fluctuations also increases

The time history of water level was investigated to determine the frequency of fluctuation.

The power spectrum analysis was done for the water depth data. Figures 6.22 – 6.23 show

the measured oscillation spectra at point 3 and point 4 in Exp - 3. It can be easily seen

from these experimental results that there are some significant contributions of the

fluctuation with the periods approximately of 0.065~0.066s. These periods are

consistent with the calculated results shown in Figure 6.24 and 6.25. In these graphs, the

period of fluctuation in the calculated results is approximately 0.06s.

The same picture is also monitored from the calculation for Exp – 4 in Figures 6.26-6.29.

In observation, the significant periods are 0.065-0.069(s), while the computed one is

0.06(s). The response pattern in the calculated results is similar to the experimental data

toward downstream.

in term of frequency. It is furthermore confirmed by the comparison shown in Figures

6.30-6.33, in which the calculated results are for Exp-8. The observation periods are

0.075-0.087 (s) while the simulation one is 0.08 (s).

But, it is clearly seen from these figures that the amplitude of oscillation for the

83

calculated results is still larger than that of experiments.

0.60

0.55

)

m

(

Time averaged in Exp.

z

0.50

Water surface

Bottom

0.45

0.90

0.95

1.05

1.10

1.00 x (m)

Figure 6.14 Computed water surface profile in Exp1

0.60

0.55

)

m

(

Time averaged in Exp.

z

0.50

Water surface

Bottom

0.45

0.90

0.95

1.05

1.10

1.00 x (m)

84

Figure 6.15 Computed water surface profile in Exp2

0.70

0.65

0.60

)

m

(

Time averaged in Exp.

z

0.55

Water surface

0.50

Bottom

0.45

0.90

0.95

1.00

1.10

1.15

1.20

1.05 x (m)

Figure 6.16 Computed water surface profile in Exp5

0.70

0.65

0.60

)

m

(

Time averaged in Exp.

z

0.55

Water surface

0.50

Bottom

0.45

0.90

0.95

1.00

1.10

1.15

1.20

1.05 x (m)

85

Figure 6.17 Computed water surface profile in Exp6

0.60

0.55

)

m

Time averaged in Exp.

( z

Water surface at T=10.00s

0.50

Water surface at T=10.15s

Bottom

0.45

0.90

0.95

1.05

1.10

1.00 x (m)

Figure 6.18 Computed water surface profile in Exp3

0.60

0.55

)

m

Time averaged in Exp.

( z

Water surface at T=10.00s

0.50

Water surface at T=10.15s

Bottom

0.45

0.90

0.95

1.05

1.10

1.00 x (m)

86

Figure 6.19 Computed water surface profile in Exp4

0.70

0.65

0.60

)

m

Time averaged in Exp.

( z

0.55

Water surface

0.50

Bottom

0.45

0.90

0.95

1.00

1.10

1.15

1.20

1.05 x (m)

Figure 6.20 Computed water surface profile in Exp7

0.70

0.65

0.60

)

m

Time averaged in Exp.

( z

0.55

Water surface at T=10.00s

Water surface at T=10.05s

0.50

Bottom

0.45

0.90

0.95

1.00

1.10

1.15

1.20

1.05 x (m)

87

Figure 6.21 Computed water surface profile in Exp8

0.020

0.018

0.015

) s

0.013

2 m m

(

0.010

m u r t c e p s

0.008

s 6 6 0 . 0 = T

r e w o P

0.005

0.003

0.000

10

15

25

30

20 Frequency (hz)

Figure 6.22 Power spectrum of water surface displacement at point 3 in Exp – 3

0.020

0.018

0.015

) s

s 5 6 0 . 0 = T

0.013

2 m m

(

0.010

m u r t c e p s

0.008

r e w o P

0.005

0.003

0.000

10

15

25

30

20 Frequency (hz)

88

Figure 6.23 Power spectrum of water surface displacement at point 4 in Exp – 3

Cal.

Exp.

0.003

)

m

(

t

0.001

n e m e c a

l

0.000

p s

i

d

-0.002

e c a

f r u S

-0.003

0

0.2

0.4

0.6

0.8

1

Time (s)

Figure 6.24 Comparison of calculated and experimental results at point 3 in Exp - 3

Cal.

Exp.

0.003

)

m

(

0.001

t n e m e c a

l

0.000

p s

i

d

-0.002

e c a f r u S

-0.003

0

0.2

0.4

0.6

0.8

1

Time (s)

Figure 6.25 Comparison of calculated and experimental results at point 4 in Exp - 3

89

0.020

0.018

0.015

) s

0.013

2 m m

(

0.010

m u r t c e p s

s 9 6 0 . 0 = T

s 5 6 0 . 0 = T

0.008

r e w o P

0.005

0.003

0.000

10

15

25

30

20 Frequency (hz)

Figure 6.26 Power spectrum of water surface displacement at point 3 in Exp – 4

0.060

0.050

) s

0.040

2 m m

(

s 7 6 0 . 0 = T

s 5 6 0 . 0 = T

0.030

m u r t c e p s

0.020

r e w o P

0.010

0.000

10

15

25

30

20 Frequency (hz)

90

Figure 6.27 Power spectrum of water surface displacement at point 4 in Exp – 4

Cal.

Exp.

0.003

)

m

(

0.001

l

0.000

i

-0.002

t n e m e c a p s d e c a f r u S

-0.003

0

0.2

0.4

0.6

1

1.2

1.4

0.8 Time (s)

Figure 6.28 Comparison of calculated and experimental results at point 3 in Exp - 4

Cal.

Exp.

0.003

)

m

(

0.001

l

0.000

i

-0.002

t n e m e c a p s d e c a f r u S

-0.003

0

0.2

0.4

0.6

1

1.2

1.4

0.8 Time (s)

Figure 6.29 Comparison of calculated and experimental results at point 4 in Exp – 4

91

0.060

0.050

) s

0.040

2 m m

(

0.030

s 8 7 0 . 0 = T

m u r t c e p s

0.020

r e w o P

0.010

0.000

10

15

25

30

20 Frequency (hz)

Power spectrum of water surface displacement at point 3 in Exp – 8

Figure 6.30

0.060

0.050

) s

0.040

s 7 8 0 . 0 = T

s 1 8 0 . 0 = T

s 5 7 0 . 0 = T

2 m m

(

0.030

m u r t c e p s

0.020

r e w o P

0.010

0.000

10

15

25

30

20 Frequency (hz)

Figure 6.31

Power spectrum of water surface displacement at point 4 in Exp – 8

92

Cal.

Exp.

0.003

)

m

(

0.001

l

0.000

i

-0.002

t n e m e c a p s d e c a f r u S

-0.003

0

0.2

0.4

0.8

1

1.2

0.6 Time (s)

Figure 6.32 Comparison of calculated and experimental results at point 3 in Exp – 8

Cal.

Exp.

0.003

)

m

(

0.001

l

0.000

i

-0.002

t n e m e c a p s d e c a f r u S

-0.003

0

0.2

0.4

0.8

1

1.2

0.6 Time (s)

Figure 6.33 Comparison of calculated and experimental results at point 4 in Exp – 8

93

6.5 2D simulation

In this manuscript, in order to examine the 2D characteristics of the flow over circular

surface, the 2D equations were applied with the same discretization scheme as described

in Chapter 5.

The similar boundary conditions to that used in the previous 1D simulation were applied,

and the flows conditions were exactly as same with experiment as listed in table 6.1. After

getting equilibrium state, the results of 2D simulation are shown in Figure 6.34 – 6.41.

Generally, the 2D simulation is consistent with the 1D simulation. In Exp-1 to Exp-3 and

And in Exp-4 and Exp-8, when the discharge is larger, some oscillation at surface can be

observed.

Exp-5 to Exp7, the water surface is stable since upstream discharge is relatively small.

6.6 Summary

To investigate the oscillation of the flows over circular surface, the depth-averaged

equations based on a generalized curvilinear coordinate attached to bottom were applied.

A hydraulic experiment was conducted in laboratory. When the upstream discharge was

relatively small, the water surface was stable, and with increased discharge, the

oscillation of water surface occurred in the circular bottom region. The intensity of the

oscillation increased toward downstream. Both steady and unsteady analysis with small

discharges showed the good agreement with experimental results. Both 1D and 2D

unsteady simulations showed that, the model has been able to reproduce the phenomena

as well as the frequency of the oscillation. However, there are still some discrepancies in

94

term of oscillation amplitude.

0.25

0.20

0.0

0.15

z

0.1

(

0.10

m

)

0.2

0.05

( m )

0.3 x

0.4

0.00

0.05

0.5

0.10

0.15

y(m)

Carpet plot of water surface in 2D simulation of Exp-1

Figure 6.34

0.25

0.20

0.0

0.15

z

0.1

(

0.10

m

)

0.2

0.05

( m )

0.3 x

0.4

0.00

0.05

0.5

0.10

0.15

y(m)

Carpet plot of water surface in 2D simulation of Exp-2

95

Figure 6.35

0.25

0.20

0.0

0.15

z

0.1

(

0.10

m

)

0.2

0.05

( m )

0.3 x

0.4

0.00

0.05

0.5

0.10

0.15

y(m)

Carpet plot of water surface in 2D simulation of Exp-3

Figure 6.36

0.25

0.20

0.0

0.15

z

0.1

(

0.10

m

)

0.2

0.05

( m )

0.3 x

0.4

0.00

0.05

0.5

0.10

0.15

y(m)

Carpet plot of water surface in 2D simulation of Exp-4

96

Figure 6.37

0.30

0.20

0.0

z

0.1

(

m

0.2

)

0.10

0.3

0.4

( m )

0.00

x

0.5

0.6

0.7

0 . 1 5

0 . 0 0 0 . 0 5 0 . 1 0 y(m)

Carpet plot of water surface in 2D simulation of Exp-5

Figure 6.38

0.30

0.20

0.0

z

0.1

(

m

0.2

)

0.10

0.3

0.4

( m )

0.00

x

0.5

0.6

0.7

0 . 1 5

0 . 0 0 0 . 0 5 0 . 1 0 y(m)

97

Carpet plot of water surface in 2D simulation of Exp-6 Figure 6.39

0.30

0.20

0.0

z

0.1

(

m

0.2

)

0.10

0.3

0.4

( m )

0.00

x

0.5

0.6

0.7

0 . 1 5

0 . 0 0 0 . 0 5 0 . 1 0 y(m)

Carpet plot of water surface in 2D simulation of Exp-7

Figure 6.40

0.30

0.20

0.0

z

0.1

(

m

0.2

)

0.10

0.3

0.4

( m )

0.00

x

0.5

0.6

0.7

0 . 1 5

0 . 0 0 0 . 0 5 0 . 1 0 y(m)

98

Figure 6.41 Carpet plot of water surface in 2D simulation of Exp-8

6.7 References

1. Azmad Z. 2005. Flow measument using free overfall in inverted semi-circular

channel. Flow. Meas. Instrum., Vol 16, pp 21-26.

2. Davis, A. C., Ellett, B. G. S. and Jacob, R. P. 1998. Flow measurement in sloping

channels with rectangular free overfall. J. Hydr. Engrg., ASCE, Vol. 124 (7), pp.

760-763.

3. Dey, S. 1998. End depth in circular channels. J. Hydr. Engrg., ASCE, Vol. 124 (8), pp.

856-862.

channel flow measurement. Flow. Meas. Instrum., Vol 13, pp 209-221.

5. Dey, S. 2003. Free overfall in inverted semicircular channels. J. Hydr. Engrg., ASCE,

Vol. 129(6), pp. 438-447.

6. Dey S. 2005. End depth in U-shaped channels: a simplified approach. J. Hydr.

Engrg., ASCE, Vol. 131(6), pp. 513-516.

7. Guo, Y. 2005. Numerical modeling of free overfall. J. Hydr. Engrg., ASCE, Vol.

4. Dey, S. 2002. Free overfall in circular channels with flat base: a method of open

131(2), pp. 134-138.

8. Pal M. and Goel A. 2006. Prediction of the end-depth ratio and discharge in

semi-circular and circular shaped channels using support vector machines. Flow.

Meas. Instrum., Vol 17, pp 49-57.

9. Rouse, H. 1936. Discharge characteristic of the free overfall. Civ. Engrg. ASCE, Vol.

6(4), pp. 257-260.

99

Chapter 7

MODEL REFINEMENT

7.1 Preliminary

In Chapter 3, a 2D depth averaged model was introduced based on a generalized

curvilinear coordinate system, in which the coordinate axis ζ was always perpendicular

to the bottom surface. Therefore, the orthogonal conditions in Equation 3.10 and 3.11

were applied to obtain the simplified momentum equations as described in (3.33) and

(3.34). Nevertheless, because axis ζ is always perpendicular to the bottom surface thus,

for a concave bed there is a possibility that ζ will intersect some where above the

surface. To avoid this difficulty, the above model should be applied in a convex

topography as in Chapter 4-6 or we limit the applicability of the model in concave

surfaces with a shallow depth therefore with these limitations the condition for ζ not to

intersect within the computational domain is satisfied as shown in following figure

100

(Figure 7.1).

ζ

ξ

In order to eliminate these limitations from the previuos model, we are now at the stage to

derive a more general equation set that will not use the condition for ζ to be

perpendicular to the bottom surface, (i.e. ζ could be an arbitrary axis).

Figure 7.1 Illustration for limitation of the model in concave topography

7.2 Non-orthogonal coordinate system

Coordinate system

The new curvilinear coordinate system is set on the bottom surface similar with the

previous application but the ζ-axis is arbitrary (as shown in Figure 7.2). The kinematic

boundary condition is also derived at the free surface as in (3.23):

(7.1)

U

V

W

s

s

s

+ + = ∂ h η ∂ ∂ h ξ ∂ ∂ h ∂ t

Depth-average equations:

The transformed RAN equations in a new curvilinear coordinates are integrated through

the depth from bottom to water surface with respect to ζ-axis, without applying the

Continuity equation

101

perpendicular relations, we obtain:

+

+

=

0

∂ ∂ ξ

∂ ∂ η

1 J

∂ h ∂ t

M J

N J

0

0

0

(7.2)

2

2

+

+

+

+

+

+

( M

MN

NM

N

ξ ξξ

ξ ξη

ξ ηξ

ηη

Γ 0

Γ 0

Γ 0

Γ 0

∂ ∂ t

M J

UM J

VM J

0

0

0

1 hJ 0

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ η ⎝

⎞ ⎟⎟ ⎠

∂ ⎛ ⎜⎜ ∂ ξ ⎝

⎞ ⎟⎟ ⎠

⎛ ⎜⎜ ⎝

h

h

Momentum equations:

ξ

0

0

0

0

0

2 0 x

x

x

y

z

ζ

ξξξ 2 2 0 0 y z J

ηξηξηξ 0 y z J

ξ τ b ρ J

0

0

0

0

h

h

τ

τ

ζ

+ + + + − − − = G h J ∂ p ∫ d ρη ∂ 0 ∂ p ζ ∫ d ρξ ∂ 0

ζ d

0

0

0

0

0

( ζξζξζξ 0 y z

x

y

x

z

ξξ ρ

ξη ρ

0

) ⎛ ⎜⎜ ⎝

00

00

b

+ + + + − d (7.3) ∂ ξ ∂ ∂ η ∂ 1 J 1 J 1 J ⎞ p ⎟⎟ ρ ⎠

2

2

and

( M

MN

NM

N

η ξξ

η ξη

η ηξ

ηη

N J

UN J

VN J

0

0

0

1 hJ 0

h

h

+ + + + + + Γ 0 Γ 0 Γ 0 Γ 0 ∂ ∂ t ⎛ ⎜⎜ ⎝ ⎞ ⎟⎟ ⎠ ∂ ⎛ ⎜⎜ ξ ∂ ⎝ ⎞ ⎟⎟ ⎠ ∂ ⎛ ⎜⎜ η ∂ ⎝ ⎞ ⎟⎟ ⎠

η

0

0

0

0

0

2 0 x

x

x

y

z

ζ

G

ηηη 2 2 0 0 y z J

ηξηξηξ 0 y z J

h J

η τ b ρ J

0

0

0

0

h

h

τ

τ

+

+

+

+

ζ

(7.4)

d

ζ d

0

0

0

0

0

( ζηζηζη 0 y z

x

y

x

z

ηξ ρ

ηη ρ

∂ ∂ η

∂ ∂ ξ

1 J

1 J

1 J

0

) ⎛ ⎜⎜ ⎝

⎞ p ⎟⎟ ρ ⎠

00

00

b

In order to estimate the pressure distribution along the ζ-axis, the momentum equation

in ζ-direction is recast with the assumptions of homogeneity in ζ and shear stress are

negligible:

ζ

2

2

2

2

2

)

)

+ + + + − − − = ∂ p ∫ d ρη ∂ 0 ∂ p ζ ∫ d ρξ ∂ 0

UV

VU

V

( U

ζ Γ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

( ζ x

ζ y

ζ z

G J

1 J

1 J

)

)

= − + + ∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρζ ∂ ⎠ ⎝

( ηζηζηζ y z

x

x

y

z

( ξζξζξζ y z

y

x

x

z

− + + + + − (7.5) 1 J 1 J ∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρξ ∂ ⎠ ⎝ ∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρη ∂ ⎠ ⎝

2

2

2

)

)

)

or

ζ y

ζ z

( ηζηζηζ y z

x

x

y

z

( ξζξζξζ y z

x

x

y

z

( ζ x

102

+ + + + + + + + ∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρξ ∂ ⎠ ⎝ ∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρη ∂ ⎠ ⎝ ∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρζ ∂ ⎠ ⎝

ζ

2

2

=

Γ

G

UV

VU

V

( U

ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

(7.6)

To solve equation (7.6) by iteration method, the first estimation of pressure term is given

( ) 0

ζ

2

2

2

2

2

)

as:

( U

ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

0

0

0

( ζ x

ζ y

ζ z

= − Γ + + G UV VU V (7.7) ∂ ⎞ ⎛ p ⎟⎟ ⎜⎜ ρζ ∂ ⎠ ⎝

( ) 0

ζ

2

2

hence, taking integration from ζ to h gives

( U

[ G

]( ) h

ζ Γ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

2

2

2

)

0

0

ζ z

0

0

− − UV VU V + + p = ( ζρ x 1 ζ y

)ζ−

( )( f p h

Consequently,

( ) 0

(

)

) ζ

=

ζ ⋅−

( h

f

h

f

(7.8)

( ) 0 p

( ) 0 p

( ) ∂ 0 f p ξ ∂

( ) ∂ 0 f p ξ ∂

and similarly,

( ) 0

)

) ζ

ζ ⋅−

+ = − ∂ ∂ h = ξξ ∂ ∂ ∂ ⎛ p ⎜ ⎜ ρξ ∂ ⎝ ⎞ ⎟ ⎟ ⎠

( h

f

h

f

(7.9)

( ) 0 p

( ) 0 p

( ) ∂ 0 f p η ∂

( ∂ ∂ h = ηη ∂ ∂

( ) ∂ 0 f p η ∂

Substituting Equations (7.8) and (7.9) into (7.7) and taking integration we attain the

equation for the new estimation of pressure distribution as follows

2

ζ 2

(

)

+ = − ∂ ⎛ p ⎜ ⎜ ρη ∂ ⎝ ⎞ ⎟ ⎟ ⎠

) ζ

)( hh

f

( ξζξζξζ y z

x

x

y

z

( ) 0 p

0

( ) ∂ 0 f p ξ ∂

h 2

2

( ) 1

2

ζ 2

2

2

2

)

)

) ζ

+ − − − + + ∂ ξ ∂ ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ ⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

)( hh

( ηζηζηζ y z

y

x

x

z

( ) 0 p

( ζ x

ζ y

ζ z

0

0

( ∂ η ∂

( ) ∂ 0 f p η ∂

ζ

2

2

Γ

=

UV

VU

V

]( ) h

[ G

( U

ζ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

− + + + + − − − f p ρ h 2 2 ⎛ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠ ⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

103

as a result,

( ) 1

ζ

2

2

2

2

)

[ G

( U

]( ) h

ζ Γ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

( ζζ 2 + y

x

ζ z

0

2

ζ 2

(

)

+ −= − − UV VU V p ρ

) ζ

)( hh

( ξζξζξζ y z

y

x

x

z

( ) 0 p

0

( ) ∂ 0 f p ξ ∂

2

ζ 2

)

) ζ

− − − + + + f ∂ ξ ∂ h 2 2 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ ⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

)( hh

( ηζηζηζ y z

x

y

x

z

( ) 0 p

0

( ∂ η ∂

( ) ∂ 0 f p η ∂

+ + + − − − f (7.10) h 2 2 ⎛ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎠ ⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

( ) 1

2

ζ

2

2

2

2

2

)

)

]

Taking integration Equation (7.10) gives

ζ

[ ( U

ζ Γ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

( ζ x

ζ y

ζ z

0

ρ

0

2

3

(

)

)

+ + = − d UV VU V G ph ∫ h 2

h

f

( ξζξζξζ y z

x

y

x

z

( ) 0 p

0

( ) ∂ 0 f p ξ ∂

h 2

h 3

2

3

)

)

− + + + ∂ ξ ∂ ⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

h

f

(7.11)

( ηζηζηζ y z

y

x

x

z

( ) 0 p

0

( ) ∂ 0 f p η ∂

( ∂ η ∂

h 2

h 3

+ + + −

Equation (7.11) describe the pressure distribution along the ζ-axis. It is easily seen that,

two last terms are added in the pressure distribution and if ζ normal to the bottom

surface, the distribution reduced as the first term - a term representative for the

hydrostatic pressure and effects of centrifugal force due to bottom curvature shown in

Chapter 3. The pressure at the bottom is then obtained:

2

2

)

⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

[ ( U

UV

VU

V

] + ζ hG

ζ ηη

ζ +Γ ξξ

ζ +Γ ξη

ζ +Γ ηξ

2

2

2

)

b

( ζ x

1 ζ y

ζ z

2

(

)

)

= Γ − + + ⎞ p ⎟⎟ ρ ⎠ ⎛ ⎜⎜ ⎝

( ξζξζξζ y z

y

x

x

z

( ) 0 p

0

( ) ∂ 0 f p ξ ∂

2

)

)

+ + + − h h f ∂ ξ ∂ h 2 ⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

( ηζηζηζ y z

x

x

y

z

( ) 0 p

0

( ∂ η ∂

( ) ∂ 0 f p η ∂

104

+ + + − h h f (7.12) h 2 ⎧ ⎪ ⎨ ⎪⎩ ⎫ ⎪ ⎬ ⎪⎭

ξ

The contravariant components of gravitational vector are derived similar as in Chapter 4:

( yxgJ

ζη

)ηζ yx

η

−= − G (7.13a)

( yxgJ

ζξ

)ξζ yx

ξ

−=

G

( yxgJ

−= − G (7.13b)

ηξ

)ξη yx

(7.13c)

7.3 Application

In order to verify the validity and the applicability of the new general equations, in this

first version of manuscript, the dam-break flows in a horizontal and sloping channel bed

In this application, the 1-D flow is considered using the new coordinate system as

α and

α define the slope of the channel and

described in Figure 7.2. The parameters

1

2

the angle of the ζ-axis with the channel bottom (also means the angle between two new

α , it is easily to set the slope of channel and

axes). By changing the value

1

0

are investigated.

1

z

ζ

α 2

α 1

x

ξ

=α defines for the case of horizontal bed.

105

Figure 7.2 Definition sketch of new generalized coordinate system

ζ

2

2

2

2

2

In this research, the second and third terms in equation (7.6) are neglected, it gives

(

)

)

( U

ζ Γ + ξξ

ζ Γ + ξη

ζ Γ + ηξ

ζ ηη

ζ x

ζ y

ζ z

0

0

0

+ + = − Γ UV VU V G (7.14) ∂ ⎞ ⎛ p ⎟ ⎜ ζ ρ ∂ ⎝ ⎠

0V = ) and the calculation gives

0

ζ ξξΓ = then pressure

Consider 1D flow (i.e.

ζ

ζ

distribution in (7.14) becomes hydrostatic:

(

) ζ

2

2

2

2

2

2

(

)

(

)

0

0

0

0

0

0

ζ x

ζ z

ζ z

ζ x

= = − G h ⇒ + + + + ∂ ⎞ ⎛ p ⎟ ⎜ ζ ρ ∂ ⎠ ⎝ ⎛ ⎜ ⎝ ⎞ p ⎟ ρ ⎠ G ζ y 1 ζ y

hence, the pressure at the bottom is obtained:

2

2

2

+

+

(

)

b

0

0

0

ζ x

ζ z

and

h

2

hζ (7.15) ⎛ ⎜ ⎝ ⎞ = p ⎟ ρ ⎠ G ζ y

ζ

G

(7.16)

d

2

2

2

(

)

2

0

0

0

0

ζ x

ζ z

1 ζ y

Substituting (7.15-7.16) into (7.3) and neglecting the effect of internal stresses gives

2

h

+

+

Γ

M

ξ

2 0 x

ζ

+

+

=

G

∂ ∂ t

M J

UM J

h J

ξ ξ ξ 2 2 0 0 y z J

0

0

ξ ξξ 0 J h 0

0

0

⎛ ⎜ ⎝

⎞ ⎟ ⎠

∂ ⎛ ⎜ ∂ ξ ⎝

⎞ ⎟ ⎠

∂ p ∫ d ∂ ξ ρ 0

= + + ⎛ ⎜ ⎝ ⎞ p ⎟ ρ ⎠

)

(

0

0

0

0

0

ξ ζ ξ ζ ξ ζ 0 y z

x

y

x

z

ξτ b ρ J

1 J

0

0

b

or

+ + − − ⎛ ⎜ ⎝ ⎞ p ⎟ ρ ⎠

2 U h

ξ

ζ

2 x 0

ξ ξξ 0

2

2

2

2

ξ ξ ξ 2 2 y z 0 0 J

)

0

0

0

0

0

ζ z

0

0

0

+ + Γ + + = − G h G + + ∂ ∂ t M J UM J J h J ⎛ ⎜ ⎝ ⎞ ⎟ ⎠ ∂ ⎛ ⎜ ξ ∂ ⎝ ⎞ ⎟ ⎠ 1 ζ y ⎡ ∂ ⎢ ( ξ ζ ∂ ⎢ 2 x ⎣ ⎤ ⎥ ⎥ ⎦

(

)

x

y

z

0

0

0

0

ζ

2

2

2

0 +

ξ τ b ρ J

(

0

0

ξ ζ ξ ζ ξ ζ 0 y z ) ζ y

x ζ x

ζ z

0

0

0

+ + − − (7.17) G h + 1 J

106

Equation (7.17) is solved in plugging with continuity equation (7.2) and the geometric

variables are calculated from the mesh depicted in Figure 7.2, to investigate the

dam-break flow in horizontal and slopping channels.

Figure 7.3 and 7.4 show the dam-break flow water surface at different times in a dried

and wetted bed for horizontal channel. The bores are visually observed, and they were

compared with the conventional depth-averaged model (i.e., Kuipers and Vreugdenhill

equations). The good agreement was achieved in the bore of both cases: dried and wetted

bed , but there are still some bias at the upstream of the dam, that might be caused by the

difference in boundary conditions at upstream end due to the difference of the grid. It

always perpendicular to the bottom surface. It also was emphasized in Figure 7.5-7.6,

which show the dam-break flows in channels with slop

030 .

107

infers that, the new model can be applied in such cases that the third axis is not necessary

2.0

T = 0.00s

T = 0.20s

T = 0.40s

T = 0.60s

)

Channel bed

m

( z

0.0

2.0

4.0

6.0

8.0

x (m)

Figure 7.3

Calculated water surface profile at different time steps of dam break

0

0

;

;

0.1=

Hini

m

flow in a dried-bed horizontal channel: 30=α

60=α

1

2

2.0

T = 0.00s

T = 0.20s

T = 0.40s

T = 0.60s

)

Channel bed

m

( z

0.0

2.0

4.0

6.0

8.0

x (m)

0

0

=

=

Figure 7.4 Calculated water surface profile at different time steps of dam break

60=α

5.0

m ;0.1

Hini

Hini

m

1

2

down

up

108

; ; flow in a wetted-bed horizontal channel: 0=α

1.5

Initial depth

Conventional depth-averaged model

New model

1.0

)

m

( z

0.5

0.0

2

4

6

8

12

14

16

18

10

x (m)

Figure 7.5

Comparison of water surface profile for dried horizontal channel at

different times: T = 0.4s; 1.0s; and 1.6s;

1.5

Initial depth

Conventional depth-averaged model

New model

1.0

)

m

( z

0.5

0.0

2

4

6

8

12

14

16

18

10

x (m)

Figure 7.6 Comparison of water surface profile for wetted horizontal channel at

109

different times: T = 0.4s; 0.6s; and 0.8s;

6.0

T = 0.0s

T = 0.2s

T = 0.4s

4.0

T = 0.6s

)

m

Channel bed

( z

2.0

0.0

4.0

8.0

10.0

2.0

6.0 x (m)

Figure 7.7

Calculated water surface profile at different time steps of dam break

0

;

;

0.1=

Hini

m

flow in a dried-bed sloping channel: 30=α 0

60=α

1

2

6.0

T = 0.0s

T = 0.2s

T = 0.4s

4.0

T = 0.6s

)

Channel bed

m

( z

2.0

0.0

2.0

4.0

6.0

8.0

x (m)

Figure 7.8

Calculated water surface profile at different time steps of dam break

0

0

=

=

flow in a dried-bed horizontal channel:

5.0

30=α

60=α

Hini

m ;0.1

Hini

m

up

down

1

2

110

; ;

Chapter 8

CONCLUSIONS

In this research, a depth-averaged model of open channel flows in a generalized

curvilinear coordinate system over an arbitrary 3D surface was developed. This model is

the inception for a new class of the depth-averaged models classified by the criteria of

coordinate system. This work focused on the derivation of governing equations with

examples of application in hydraulic amenity. The major results and contributions of this

research are summarized as follows.

1. Mathematical model

- The original RAN equations were firstly transformed into the new generalized

curvilinear coordinate system, in which two axes lay on an arbitrary surface and the

other perpendicular to that surface. The depth-averaged continuity and momentum

equation were derived by integrating the transformed RAN equation and substituting

the kinematic boundary condition at the free surface.

the perpendicular axis as a combination of hydrostatic pressure and the effect of

111

- The pressure distribution was obtained by integrating the momentum equation along

centrifugal force caused by bottom curvature.

2. Flow with air-core vortex at a vertical intake

- Firstly, 1D analytical analysis was applied to calculate the water surface profile of the

flows with air-core vortex at a vertical intake. The use of curvilinear coordinate

avoids invalidation of the Kelvin’s theorem in the core of vortex, therefore the

conservation of circulation can be applied for whole flow. The comparisons of

calculated data using the model and an empirical fomula show that the proposed

model yields reliable results in predicting the critical submergence of the intake

cannot escape.

- Secondly, 2D unsteady analysis was applied to simulate the water surface profile and

also showed the applicability of the model in computing the flow into intake with

air-entrainment.

without any limitation of Froude number, a problem that most of existing models

3. Flows over circular surface

- A hydraulic experiment was conducted in the laboratory to measure the surface

profile of the flows over circular surface. It was pointed out that: with small discharge,

the free surface remained stable but some oscillations occurred when the discharge

was increased. The intensity of the oscillation was observed to increase toward the

downstream direction of the flow.

- Both 1D steady and unsteady analysis were applied to estimate the water surface

profile of the flows with small discharge. Good agreements in comparison with

measured data were achieved for both approaches. Results from 1D unsteady model

demonstrate the fluctuation at the region with circular bottom with large discharge.

comparison with observed data showed the fairly good agreement in term of

112

Spectrum analysis was applied to investigate the characteristics of the oscillation. The

frequency or period but some discrepancies in magnitude of fluctuation were still

observed. 2D simulation was conducted, and also could reproduce the phenomena as

well.

4. Model refinement

- In order to apply the model for the case with highly concave topography, the new

depth-averaged equations were derived without using the perpendicular conditions.

The improved model therefore has the ability in more general geometry of the bottom

surface

horizontal and sloping channels.

113

- A simple example was presented to verify the model with dam-break flows in