FLOWS AND OPEN-CHANNEL FLOWS

Srinivas Chippada

Rice University
Department of Mechanical Engineering and Materials Science

Submitted to the Texas Water Development Board in partial fulfillment of

the requirements for Contract No. 94-483-010.

February 1996
 Acknowledgments

I sincerely thank my thesis advisor, Dr. Bala Ramaswamy, for his guidance, advice, and

encouragement during this research. I deeply appreciate the help and advise given to me by Dr.

Mary F. Wheeler, throughout my stay at Rice. I would also like to thank Dr. John E. Akin,. Dr.

Yves Angel and Dr. Yildiz Bayazitoglu for serving on the dissertation committee.

The original rationale for the work on open-channel flows was given by Dr. Dan Tetzlaff

at Texaco, Inc. Dr. Sang Woo Joo at Wayne State University provided the inspiration for thin-

film flow research, and I thank him sincerely for the many things I learned from him.

The financial support from Rice University, Texaco Inc., and the National Science

Foundation is greatly appreciated. An important part of the work was also funded by the Texas

Water Development Board, and I gratefully acknowledge the Board's support and interest in

hydrodynamic modeling. This research wouldn't have been possible without the excellent

computational resources offered by the Department of Mechanical Engineering and Material

Science at Rice University, I learned a lot from my colleagues S. Krishnamoorthy, V. Ravi Rao,

F.J. Eaton, Uday B. Sathuvalli, and R. Moreno, with whom I had many fruitful discussions on

the subject of this research.

My heartfelt regards and deepest love to my parents and my wife for their help,

understanding, and motivation to continue through even the most difficult times.
 NUMERICAL STUDY OF THIN-FILM
FLOWS AND OPEN-CHANNEL FLOWS

Srinivas Chippada

Abstract

A numerical procedure is developed capable of simulating gravity-driven film

flows in two-dimensions. The moving boundary problem is handled through the

ALE formulation. In the case of turbulent fluid flows, the two-equations k -_ c

closure model is used to model the turbulence. A Chorin-type projection scheme is

utilized to decouple the velocity and pressure fields, and the spatial discretization

is done using the Finite Element Method.

Thin liquid films draining down vertical or inclined planes are susceptible to

long wavelength disturbances. An extensive numerical study of the surface wave

instability in isothermal thin film flows is done by solving the full-scale nonlinear

system. Temporal stability analysis of a spatially periodic disturbance reveals in-

teresting wave dynamics. The transition from nearly sinusoidal supercritical waves

to broad-banded solitary waves is found to go through a quasi-periodic regime. In

this quasi-periodic state, the fundamental mode and several of its harmonics are in

an oscillatory state, with continuous exchange of energy. An extensive parametric

search has been done to obtain the phase boundary delineating the quasi-periodic

regime. Complex wave interactions such as wave-splitting and wave-merging are

discussed. Spatial stability analysis akin to the usual experimental studies is done

and comparisons are made with the experiments.

 Ill

For the development of successful theories capable of predicting the formation

of bedforms, it is essential to understand the turbulent fluid flow on top of the bed.

To this end, mean flow and turbulence characteristics for flow over artificial stream-

wise periodic bedforms are obtained. Due to the local accelerations associated with

stream-line bending, very large velocities and stresses are found to exist at the tip

of the dune. The separation wake turbulence is found to completely dominate the

wall-generated turbulence, and the maximum turbulence intensity levels occur at a

distance, approximately equal to the dune height, away from the bed. The accuracy

of the rigid-lid approximation is determined by computing the flow field with and

without the rigid-lid approximation. The rigid-lid approximation is found to over-

predict the shear at the dune crest.

Lastly, the mean flow and turbulence characteristics in hydraulic jump are ob-

tained. In the case of flow with inlet supercritical Froude number 2.0, a small

recirculation zone is found to exist at the foot of the jump. The mixing layer tur-

bulence associated with the surface roller and the recirculation zone are found to

dominate the wall-generated turbulence.

 VI

R radius of curvature
.'j specific density of sediment

total bed load transport across a cross-section

Ss total suspended load transport across a cross-section

time

T non-dimensional transport parameter

0.t time increment

1l velocity component in the x-direction

shear velocity at the bed

U*,CT critical shear velocity for threshold motion

u mean velocity in the x-direction

v velocity component in the y-direction

mesh particle velocity in the x direction

mesh particle velocity in the y direction

Ws suspended sediment particle fall velocity

(x,y,z) Cartesian coordinates

X Material coordinates

X spatial coordinates

:X referential coordinates

z non-dimensional suspension parameter

X; Cartesian coordinate in i-th direction

~x, ~y spatial increment in x and y directions

Greek Letters

turbulent dissipation rate

 VII

von Karman const.

coefficient of viscosity

(:x, :J, gradient operator

v laminar kinematic viscosity

Vy turbulent kinematic viscosity

shape function

p fluid density

surface tension
(]' ..
I) stress tensor

T tangential direction

critical shear stress for threshold motion given by Shields

shear stress at the wall
() angle of inclination of the bed with the horizontal; also,

referred to as energy slope

Subscripts

i,j, k directions of Cartesian coordinate

t derivative with time

X derivative w. r. t. x

o:,/3,/ nodal number

partial derivative w. r. t. x;

Superscripts

n n-th time steps

refers to dimensional value
*
fluctuating quantity
 Vlll

Nondimensional Numbers

Fr Froude number. Fr = ·JgH

,u
c: Thickness parameter. C: = 2!$
v·

Re Revnolds
'
number. Re = uvH
s Surface tension parameter. S = 7
h"
3 pv 2

T J / (:3pv"f3gt!3)

We \Neber number. IV e = u
J!i

)/ote: Only the most important symbols are listed above. The symbols are defined

the first time they are used. Symbols are also subject to alteration on occasion.
 Table of Contents

Abstract ll

Acknowledgments IV

Nomenclature v
List of Figures Xlll

List of Tables XVll

1 Introduction -1

1.1 Motivation 0 1
1.2 Course of the study 2

2 Flow Over An Inclined Plane 5

201 Introduction 0 0 0 0 0 5
202 Governing Equations 5
203 Boundary Conditions 7
2.4 Nondimensionalization 9

3 Numerical Modeling of the Free Boundary Problem 11

301 Introduction o o 11

302 ALE formulation 13

303 The concept of vertical spines 20

4 Numerical Modeling of Turbulence 22

401 Introduction 0 0 0 0 0 0 0 0 22
402 Reynolds Averaged Equations 23
 X

-b.:3 Eddy viscosity hypothesis 24

-L-± Eddy viscosity models 26
4.5 k - c turbulence model 28
4.6 Boundary Conditions for k- <: Model 30
4.1 Surface Roughness ..... . 33
4.8 .\Tote on the RANS equations 35

.'5 """america! Procedure 36

5.1 Introduction . . 36
- ')
.J._ Time stepping procedure 37

.5.3 Spatial discretization 38

5.4 Mesh generation . . . 41

5.5 Miscellaneous details 43

6 Interfacial Instabilities in Thin Film Flows 46

6.1 Introduction . . . . . . . . . . . . . 46

6.2 Problem Definition and Non-dimensional Parameters 50

6.3 Literature Review . . . . . 52

6.3.1 Nusselt Film Flow 52
6.3.2 Linear stability .. 54

6.3.3 Nonlinear Stability Analysis Based On Lubrication

Approximation . . . . . . . . . . . . . . . . . . . . . . . . 58

6.3.4 Nonlinear Stability Analysis based on Boundary Layer

equations ........ . 62
6.3.5 Direct Numerical Studies . 65

6.3.6 Experimental Studies 66

6.4 Course of the Study . . . . . 67

 XI

6.5 Comparison with experiments and previous numerical simulations 69

6.6 Temporal Stability Analysis . . . . . . . . . . . . . . 81

6. 7 Spatia-temporal evolution of the thin film instability 94

6.7.1 Long domain with periodic boundary conditions 94

6. 7.2 Wave breaking in thin film flows . . . . . . . . . 103

6.7.3 Long domain with non-periodic boundary conditions 104

6.8 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . 118

7 Turbulent Fluid Flow Over Streamwise Periodic Artificial Bedforms 122

7.1 Introduction . . . . . . . . . . . . . . . . . . . 122

7.2 Problem Definition and Fluid Flow Modeling . 126

7.3 Sediment Transport Modeling 129

7.4 Results and Discussion . . . . 136

7.4.1 Mean flow and turbulence characteristics with the free

surface assumed to be rigid frictionless lid . . . . . . . . . . 138

7.4.2 Mean flow and turbulence characteristics without the

rigid-lid assumption. . 152

7.5 Concluding remarks . 161

8 Hydraulic Jump 164

8.1 Introduction 164

8.2 Literature Review . 167

8.3 Problem Definition . 170

8.4 Results and Discussion 173

8.5 Concluding remarks . 179

9 General Conclusions 195

9.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195

 XII

9.2 Future Work . 197

Bibliography 199
 List of Figures

1.1 Gravity-driven flow over an inclined plane . . . . . . . . . . . . . . . 3

3.1 Referential kinematic description . . 15

3.2 Free surface resting on vertical spines 21

4.1 The near-wall node placement.. . . . . . . . . . . . . . . . . . . . . . 31

6.1 Surface wave regimes in thin film flows. . . . . . . . . . . . . . . . . . 49

6.2 The physical configuration of a thin layer flowing down an inclined

plane. . . . . . . . . . . . . . . . . . . . . . . . . . 51

6.3 Stability diagram based on linear stability analysis . 59

6.4 Stability diagram based on weakly nonlinear stability analysis 63

6.5 Nonlinear wave evolution for the first experimental condition of

Kapitza & Kapitza (1949) . . . . . . . . . . . . . . . . . . . . . . 73

6.6 Comparison of the wave profile for the first experimental condition

of Kapitza & Kapitza (1949) . . . . . . . . . . . . . . . . . . . . . . 74

6. 7 Mesh and time step independent study for the first experimental
condition of Kapitza & Kapitza (1949) . . . . . . . . . . . . . . . . . 75

6.8 Nonlinear wave evolution for the second experimental condition of

Kapitza & Kapitza (1949) . . . . . . . . . . . . . . . . . . . . . . . . 78
 XIV

6. 9 Comparison of the wave shapes for the second experimental

condition of I\.apitza & Kapitza ( 1949) . . . . . . . . . . . . . 79

6.10 Influence of the initial condition on the nonlinear wave evolution for

the second experimental condition of Kapitza & Kapitza (1949) SO

6.11 ~onlinear wave evolution for G = .5.0 . 84

6.12 Nonlinear wave evolution for G = 25.0 86

6.13 ~onlinear wave evolution for G = 100 . 88

6.14 Free surface profiles at different instants of time for G = 100 and

k = 0.20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

6.15 Nonlinear wave evolution for G = 25 and various wavenumbers 91

6.16 Phase diagram for thin film flow obtained through full-scale

computations. . . . . . . . . . . . . . . . . . . . . . . . . . . 92

6.17 Spatia-temporal evolution in a long periodic domain for G =5 95

6.18 Spatia-temporal evolution in a long periodic domain for G = 5 with

an initial exponential pulse . . . . . . . . . . . . . . . 101

6.19 Rosenau & Oron (1989)'s wave breaking simulation . 105

6.20 Spatia-temporal evolution with f = l.5H z periodic forcing at the inlet107

6.21 Comparison with Liu & Gollub (1994)'s experiments . . . . . . . . . . 111

6.22 Spatia-temporal evolution with f = 4.5H z periodic forcing at the inlet 113

6.23 Comparison with Liu & Gollub (1994)'s experiments. 116

6.24 Comparison with Liu & Gollub (1994)'s experiments 117

6.25 Simulation of Kapitza & Kapitza ( 1949) 's experimental conditions

in a long non-periodic domain . . . . . . . . . . . . . . . . . . . . . . 119

7.1 Schematic description of fluid flow over periodic bedforms . . . . . . . 124

 XV

7.2 Typical Finite Element mesh . . . . . . . . . . . . . . . . . . 130

7.:3 Bedform shapes: (a) Type I bedform with 45° upstream and

downstream faces; (b) Type II bedform with 45° downstream slope

and a gentle upstream slope . . . . . . . . . . . . . . . . . . . . . . . . 137

7.4 Mean flow characteristics for flow over Type I bedform with rigid-lid
approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

7.5 Mean flow characteristics for flow over Type I bedform with rigid-lid

approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

7.6 Shear stress and surface pressure at the air-liquid interface for flow

over Type I bedform with rigid-lid approximation . . . . . . 142

7. 7 Turbulence characteristics for flow over Type I bedform with

rigid-lid approximation . . . . . . . . . . . . . . . . . . . . . 145

7.8 Turbulence characteristics for flow over Type I bedform with

rigid-lid approximation . . . . . . . . . . . . . . . . . . 147

7.9 Mean flow characteristics for flow over Type II bedform 148

7.10 Turbulence characteristics for flow over Type II bedform 150

7.11 Shear stress and sediment load distribution for flow over Type II

bedform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151

7.12 Mean flow characteristics for flow over Type I bedform without the

rigid-lid approximation . . . . . . . . . . . . . . . . . . . . . . . . . . 153

7.13 Turbulence characteristics for flow over Type I bedform without the

rigid-lid approximation . . . . . . . . . . . . . . . . . . . . . . . . . . 155

7.14 Effect of the rigid-lid approximation on the shear stress distribution

for flow over Type I bedform . . . . . . . . . . . . . . . . . . . . . . . 156
 XVI

7.1.5 Time history of free smface height for the case of f1m\· over Type II
bedform .. . L.58
7.16 Propagation of the surface wave in the case of flow owr Type II

bedform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
7.17 Time history of the shear stress at the bed in the case of flow over

Type II bedform without making the rigid-lid approximation. . . 160

7.18 Propagation of the surface wave in the case of flow over Type II

bedform with energy slope 0.00-!. . . . . . . . . . . . . . . . . L62

8.1 Schematic sketch of a hydraulic jump 166

8.2 Various types of hydraulic jump (adapted from Chow (1959)) . HiS

8.3 Typical Geometry for a hydraulic jump . 181

8.4 Finite Element Mesh for part of the domain 182

8.5 Velocity vector plots. . . . . 185

8.6 Mean flow field for F1 = 2(UD). 186

8.7 Mean flow field for F1 = 2(UD). 187

8.8 Turbulence characteristics for F 1 = 2(UD). 188

8.9 Turbulence kinetic energy budget for F 1 = 2(UD) 189

8.10 Mean flow for F 1 = 2(FD) ..... 190

8.11 Turbulence characteristics for F1 = 2(FD). 191

8.12 Turbulence kinetic energy budget for F 1 = 2(FD). . 192

8.13 Mean flow and turbulence characteristics for F 1 = 4(UD) 193

8.14 Variation of skin friction coefficient C1 along the jump 194

 List Of Tables

6.1 Experimental conditions of Kapitza & Kapitza(1949) and the

relevant non-dimensional parameters. . . . . . . . . . . . . . . . . . . 69

8.1 List of simulation parameters . . . . . . . . . . . . . . . . . . . . . . 183

 1

Chapter 1

Introduction

1.1 Motivation

In recent years, due to the tremendous progress achieved in computer architecture

and computational algorithms, numerical modeling is being increasingly used in the

study of many scientific and engineering problems. Computational Fluid Dynamics

(CFD) has evolved into a powerful tool and is used not only in simulating industrial

and naturally occurring fluid flows, but also in the deeper understanding of the

subject of Fluid Mechanics itself. For example, CFD models based on the Direct

Numerical Simulation (DNS) are being used in understanding the very difficult

problem of turbulence. In this dissertation, a numerical model is developed capable

of simulating a variety of gravity-driven free surface flows.

The fluid flow problems numerically simulated are the thin-film flows and open-

channel flows. Both, thin-film flows and open-channel flows are geometrically similar

and belong to the same class of gravity-driven flows, shown in Fig.l.l. Gravity is

the driving force and drag at the solid wall due to fluid viscosity is the opposing

force. The difference between the two flows, is mainly due to the vastly different

length scales. In thin-film flows, the film thickness is very small and is in the

order of millimeters. Consequently, surface tension is very important and generates

surface waves of the gravity-capillary type on the gas-liquid interface. The fluid

flow, however, is laminar in nature. In open-channel flows, the film is very thick,
usually in the order of meters, and the fluid flow is almost always turbulent in
 2

nature. Surface tension is not significant in these scales, and the interfacial waves
are primarily of the gravity-type.

The motivation for studying thin-film flows arise from their film cooling and film

coating applications. The interfacial waves on the gas-liquid interface are known to

increase heat and mass transfer in both the liquid and gaseous phases. Open-channel

flows have applications in several areas of hydraulics, geology and sedimentology. In

the vast field of open-channel flow hydraulics itself, we study in detail two fluid flow

phenomena. First is the study of turbulent fluid flow over artificial bedforms. A

detailed mean flow and turbulence field is essential for understanding the formation

of bedform features such as ripples, dunes and anti-dunes. The second problem

studied is the hydraulic jump, which is known to occur in supercritical open-channel

flows. The free surface rises abruptly and intense turbulent mixing takes place in

the jump region. Due to their energy dissipating capability, the hydraulic jumps are
widely used in dams, spillways and other outlet works.

1.2 Course of the study

The fluid flow is governed by the continuum conservation of mass and conservation

of momentum equations, and these are described in chapter 2. Also stated are the

important non-dimensional numbers arising in the gravity-driven fluid flows. The

various flow regimes pertaining to the simple flow over an inclined plane are also

discussed.

The biggest challenge in simulating these gravity-driven flows comes from the

presence of a free boundary, which in our case is the air-liquid interface. Its location

is usually not known a priori, and needs to be determined as part of the solution.

A mixed Eulerian- Lagrangian procedure called the Arbitrary Lagrangian- Eulerian

 3

Air

Inlet

Fluid /

~X

Figure 1.1: Gravity-driven flow over an inclined plane

(ALE) formulation is employed to handle the moving boundary, and the essential

aspects of this procedure are described in chapter 3.

In open-channel flows, the fluid flow is almost always turbulent in nature. Even

though it is believed by many, that an accurate and direct solution of the Navier-

Stokes equations can simulate turbulence, the present computational resources pro-

hibit such an effort for realistic fluid-flow problems. At the moment, there is no

recourse except to solve the Reynolds Averaged Navier-Stokes (RANS) equations,

and this is what is done in this work. In addition, a two-equation k- E model along

with Boussinesq Eddy Viscosity hypothesis is used, and the essential features of this

turbulence modeling are described briefly in chapter 4.

The spatial and temporal discretization aspects of the Navier-Stokes equations

are briefly outlined in chapter 5. The presence of the boundary conditions in terms of

pressure, makes the primitive variable ( u -v- p) formulation, the popular method of

choice in the numerical simulation of these moving boundary problems. A Chorin-

 4

type splitting scheme is used to decouple the pressure from the velocities, thus

permitting us to solve for u, v and p in a sequential manner. The Galerkin Finite

Element method with three-node triangular elements are employed in the spatial

discretization.

Chapters 2-5 explain the mathematical and numerical aspects of our model. The

next three chapters deal with specific applications.

Thin-film flows are of crucial importance in many heat transfer equipment and

coating applications. A variety of fascinating interfacial dynamics occur in these

flows and these are studied in detail in chapter 6. The approach used is similar

to classical hydrodynamic stability analysis. Disturbances of different wavelengths

and frequencies are imposed on the basic solution, and the interfacial stability _is

characterized based on the response of the system to these disturbances.

In open-channel flows themselves, we simulate two very important and difficult

problems. First, is the turbulent fluid flow over artificial streamwise periodic bed-

forms (chapter 7). The formation of ripples, dunes and anti-dunes on the river bed is

due to a complex interplay between the turbulent fluid flow, sediment transport and

bed dynamics. Thus, an understanding of the turbulent fluid flow over the artificial

bedforms is an important step in the understanding of the formation of bedforms.

The second problem simulated is the celebrated hydraulic jump. The formation of

hydraulic jump in open-channel flows, is marked with the presence of high gradient

regions which are very hard to simulate. The hydraulic jump is analogous to the

shock formation in compressible flows and is studied in chapter 8.

Finally, some general conclusions and future research directions are discussed in

chapter 9.
 5

Chapter 2

Flow Over An Inclined Plane

2.1 Introduction

Though, the fluid flow problems considered in this work have application in widely

different areas of engineering, the underlying physics is the same and can be modeled

as viscous fluid flowing down an inclined plane under the action of gravity (Fig.l.l).

The fluid domain is bounded below by an impermeable and rigid wall and above

by a gas-liquid interface, which is free to deform. The planar wall is inclined at an

angle B with the horizontal. x-axis is aligned in the streamwise direction and the

and y-axis is perpendicular to it. The gas-liquid interface is represented as height

h(x, t) from the x-axis. Gravity is the driving force and the drag at the solid-liquid
interface is the opposing force. However, these two primary forces are not always in

balance, thus, giving rise to additional forces such as inertial forces, pressure forces

and capillary forces. All of these forces play an important role in determining the

internal and external dynamics of the fluid flow.

2.2 Governing Equations

In nature, the fluid flow problems are always three-dimensional. However, we re-

strict ourselves to only two-dimensional flows. Thus, the span wise ( z) direction

length scale is assumed to be infinitely large and the fluid ·.·elocities and gradients

in the z direction are neglected. \Ve further, confine ourselves to incompressible,

 6

Newtonian and isothermal flows only. The governing equations are the conservation

of momentum as given by theN a vier-Stokes equations and the conservation of mass.

au au au 1 op . 2
- + U- + V- = --- + g Sill(} + V \l U (2.1)
at ox oy pox

ov ov ov 1 op
- + u - + v - = --- - g cos() + v v2v (2.2)
at ox oy P oy

(2.3)

u = u( x, y, t) and v = v( x, y, t) are the instantaneous fluid velocities in x and

y directions respectively, and p = p(x, y, t) is the instantaneous pressure. V 2 =

(;;2 + ::2) is the Laplacian operator and V = ( :x, :J is the gradient operator._ p

and v are the liquid density and kinematic viscosity, respectively. The Eqs.2.1-2.3

are for the liquid phase. A similar set of equations are also required for the gaseous

phase. However, in the applications considered in this work, the density of gaseous

phase (usually air) is in general much smaller than the density of the liquid phase

(usually water). The gaseous phase can then be assumed to be passive and the fluid

motions in that phase can be neglected, and no equations need to be solved in the

gaseous phase. The position of the gas-liquid interface is completely determined by

the fluid dynamics in the liquid phase alone, and in such situations, the gas-liquid

interface is also referred to as a free surface. Eqs.2.1-2.3 represent a set of three

coupled partial differential equations to be solved for the three unknowns u, v and

p. However, the equations do not form a closed set since the fluid domain or the

free surface height h(x, t) is not known a priori and needs to be determined also.
2.3 Boundary Conditions

To solve the system of Eqs.:2.1-:2.:l. boundary conditions are required. On the wall

boundary. the familiar no-slip and impermeability conditions are imposed as follows:

u = 0 }
at !J = 0. (2.4)
v = 0

The boundary conditions at the free surface !J = h(x, t), are more complicated and

involve the fluid stresses. The stresses in the liquid and gaseous phase are related

through the interfacial tension as follows:

iJlln)
(-p + :2f.L- ( au) cr
on I on +-R
-p+2f.l-
9
at y = h(x. t). (2.5)
Olln iJuT)
f.ll (-
OT
+an
- I =

The subscripts l and g refer to the liquid and gaseous phases. cr is the surface

tension, and n and T are the unit normal and tangential directions on the free

surface y = h(x, t). f.l is the dynamic viscosity. R is the radius of curvature at the
free surface and is calculated as follows:

1
(2.6)
R

Since, the gaseous phase is assumed to be passive, the dynamic normal and tangen-

tial stresses in the gaseous phase are neglected and the stress boundary conditions

at the interface simplify to:

} ot y d.(x, t). (2.7)

0
 8

p 9 is the pressure in the gaseous phase and it is assumed to be spatially nonvarying

and set to a datum pressure of zero. Additional simplifications can be made depend-

ing on the practical application being modeled. In the case of open channel flows

the film depth is quite large and the gravitational and inertial forces dominate the

capillary forces and the surface tension term in the free surface boundary condition

can be neglected. However, in thin film flows this is not the case and surface tension

cannot be neglected. In fact it is the interfacial tension which is responsible for the
onset of wave motion in the case of thin film flows.

The above described boundary conditions at the free surface are only the dynamic

boundary conditions. In addition to those dynamic stress boundary conditions we

also need to satisfy the kinematic condition as follows:

(2.8)

The subscript i refers to the interface. The above condition can be reformulated to

represent the gas-liquid as a material interface resulting in the following boundary

condition at the interface:

ah ah
- + u- = v at y = h(x, t). (2.9)
at ax
The above kinematic condition is the usually imposed boundary condition at the

interface along with the stress conditions given by Eq.2. 7.

A variety of boundary conditions at the the inlet and exit boundaries have been

used, and they are stated in the later chapters dealing with applications.
 9

2.4 N ondimensionalization

Considerable insight can be gained into the physics of the fluid flow through nondi-

rnensionalization. Let H be the characteristic length and U the characteristic ve-

locity. Typically, H can be taken to be the mean film thickness and U the average

velocity. The nondimensional parameters are denoted with a superscript *. The

nondimensional independent variables are defined as:

x"=r/H: .r/=y/H: t*=WjH. ( 2.10)

The nondimensional dependent variables are defined as:

u· = u/U : u" = u/U ; p* = pj(pU 2 ) ; h* = hjH. (2.U)

Based on the above defined scaling, the nondimensional governing equations are

given by:

au• + U .au·
--
at• ax•
.au•
-- + V -
ay•
8p*
= --
ax•
1 .
+ --2
Fr
Slll () + -
Re
1 (aax• + --,
2
--2
u* a u*)
ay•
2
(2.12)

av• +u"-
-
at•
av• +v*-
ax•
av• = -8p* 1 1
- - --cosB+-
ay• ay• Fr 2 Re
(aax• + aay•v*)
2
- -2
v*
2
--2 (2.13)

au• av•
-+-=0 (2.14)
ax• ay•
The boundary conditions also need to be nondimensionalized, and they are written

as:

u* = v* = 0 at y = 0; (2.15)
 10

_1
Re
(au;+ au~)
an aT 0

') "' *
-p· + -.... -uun
- at y* = h*(x*, t*). (2.16)
Re an
v·

The nondimensional numbers that fall out from the above nondimensionalization
procedure are:

• Re = (pUH)/p., the Reynolds number 1 .

• Fr = U/ -/iH, the Froude number 2 .

• We= Uf)a"f(pH), the Weber number 3 .

The nature of the film flow is largely determined by the Reynolds number, the

Weber number, and the Froude number. Similar to the flow through pipes, the

Reynolds number determines the transition form laminar to turbulent flow. The
critical Reynolds. for transition to turbulence in the case of gravity-driven flows

over inclined planes has been found to Recrit = 250 - 500 (Fulford 1964). The
Froude number is the ratio of fluid velocity to the celerity of the gravity wave,

and determines the onset of gravity waves in water films. For Froude numbers in

the range Fr = 1 - 2, gravity waves appear on the free surface (Fulford 1964).

Weber number is the ratio of fluid velocity to the celerity of the capillary wave, and
determines the onset of capillary waves in water films. The capillary waves appear

on the water surface in the range of We= 1 (Fulford 1964).

1 Alternate definition of Reynolds number based on the hydraulic radius have been used by some
researchers.
2
Some researchers define the Froude number as Fr = U 2 /(gH).
3 Many 2
alternate definitions of the Weber number exist, such as, We = t:T /(pU H), and We =
(j /(pgHz).
 11

Chapter 3

Numerical Modeling of the Free Boundary

Problem

3.1 Introduction

Boundary /Initial value problems in which the fluid domain has an unknown physical

boundary and which needs to be determined as part of the solution procedure are

called moving boundary problems or free boundary problems ( Floryan & Rasmussen

1989). Problems of this nature are of interest in open channel flows, phase change

problems, jets, flame propagation, groundwater seepage, metal forming processes,

draining film flows, and various other fields of engineering and science. The unknown

boundary is usually an interface between liquid and gas or two immiscible liquids

or fluid and solid. Moving boundary problems are highly nonlinear and analytical

solutions are extremely difficult to obtain. Hence, numerical methods are widely

used in the modeling of free boundary problems.

All numerical methods developed to handle moving boundaries can be divided

into the following two categories: interface capturing and interface tracking (Floryan

& Rasmussen 1989). In the interface capturing methods, the physical details of the
interface are resolved. No special attention need be paid to the interface with the

possible exception of use of fine mesh in the vicinity of the interface. The nonlinear-

ity in this case is material in nature with sharply changing fluid properties across

the interface. In the interface tracking procedures, the physical details of the inter-

face are not resolved and the interface is modeled as a discontinuity and explicitly
 12

tracked. The nonlinearity in this case is a geometrical one, with an unknown fluid

domain. The fluid domains of interest to us are bounded below by a solid wall

and above a moving gas-liquid interface (Fig.l.l). The interface between the liquid

and gas has a finite thickness spanning a few molecular diameters. It is extremely·

difficult to resolve this thin layer, and hence, the interface is often modeled as a sur-

face with zero thickness (discontinuity). This gives rise to an excess surface energy,
which is nothing but the interfacial tension or the surface tension. The numerical

procedure developed and used in this work is an interface tracking algorithm.

The interface tracking algorithms can further be classified as: (i) Eulerian; (ii)

Lagrangian; and (iii) mixed Eulerian-Lagrangian methods. In the Eulerian descrip-

tion of flow, the grid points are stationary in the laboratory frame of referen~e

and the fluid moves in and out of the computational cells. The Marker and Cell

(MAC) method (Harlow & Welch 1966) and the volume of Fluid (VOF) method

(Hirt & Nichols 1981) are some of the popular and successful Eulerian methods

developed to handle moving boundary problems. The Eulerian methods can handle

arbitrarily large free surface deformations without significant loss of accuracy. The

main disadvantage of these methods is that the interface is not sharply defined and

it is difficult to impose boundary conditions at the interface. In the Lagrangian

methods, the grid points move with the local fluid velocity, and the fluid within a
computational cell remains within that cell at all times. The main advantages are

that the interface is sharply defined and the boundary conditions at the interface

are easier to impose. The main disadvantage is mesh tangling and the loss of accu-

racy due to severe mesh distortion. Hirt, Cook & Butler (1970) developed a pure
Lagrangian procedure to handle incompressible fluid flows with free boundaries.

Various mixed Eulerian-Lagrangian methods have been developed to combine the

 13

advantages of pure Lagrangian and pure Eulerian methods without their disadvan-
tages. The method used in this work to handle moving boundaries belongs to the

mixed Eulerian-Lagrangian family and is called the Arbitrary Lagrangian-Eulerian

(ALE) method.

An universal method capable of handling all types of free boundary problems is

not yet available. The numerical method that is to be chosen is strongly dependent

on the physical problem that is to be modeled. Our primary interest is the study of

liquid flow down inclined planes with a physical geometry that is usually of the form

shown in Fig.l.l. The important features of this geometry are: (i) logically rectan-

gular geometry; (ii) existence of a strong mean flow; (iii) free surface deformation

normal to the mean flow; and (iv) the representation of the air-liquid interface as_a

single-valued function of the stream-wise co-ordinate x. For this type of problem it

is advantageous to let the grid points move only in the y-direction without moving

them in the x-direction. This is what is done in our ALE formulation.

3.2 ALE formulation

The ALE description of the fluid flow is also called the referential kinematic de-

scription of flow, since the governing equations are written in a frame of reference

that is moving independently of the fluid motion. The ALE method is usually at-

tributed to Amsden & Hirt (1973). A number of extensions of this work have been

proposed since then, and has been widely used in the modeling of several fluid flow

problems ( e,g. Hirt, Amsden & Cook 1974; Chan 1975; Pracht 1975; Hughes, Liu &

Zimmerman 1981; Ramaswamy & Kawahara 1987a; Huerta & Liu 1988; Soulaimani,
Fortin, Dhatt & Ouellet 1991) and fluid-structure interaction problems (e.g. Donea,

Giuliani & Halleaux 1982; Donea 1983; Belytschko & Flanagan 1982; Liu, Chang,
 14

Chen & Belytschko 1988). The ALE formulation has been derived by Donea eta!

(1982), Ramaswamy & Kawahara (1987a), Soulaimani eta! (1991), Lacroix & Garon

(1992), among many others, and is briefly described next.

Let Bo be the open region occupied by the fluid at timet= 0 (Fig.3.1). The

position vector of point P in the material domain is denoted by X = (XI, Xz);

and X are called the material coordinates. Note that, a two-dimensional space is

being assumed, and the extension to three-dimensions is straight forward. B 1 is the

open region occupied by the fluid body Bo at some later time t > 0. The position
vector of the point p in the spatial domain B 1 is denoted by x = (xi, xz); and x
are called the spatial coordinates. The mapping from the material domain to the

spatial domain is written as:

x = ,P(X, t) (3.1)

It is assumed that the above mapping is continuous, single-valued and possesses a

unique inverse Eulerian description of the form:

X= ,p-I (x, t) (3.2)

In the Lagrangian formulation, the fluid motion is described in terms of the

material coordinates X, and the frame of reference translates along with the material

point X. The material or Lagrangian velocity is defined as:

. a
,P(X,t) = ot,P(X,t), (3.3)

and the material acceleration is defined as:

. az
¢;(X, t) = ot 2 ,P(X, t) (3.4)
 15

material domain spatial domain

<1> (X i ,t)

1\
A(x,t)

referential domain

Figure 3.1: Referential kinematic description

 16

The Lagrangian description of motion is widely used in the Solid Mechanics stud-
ies. Note that, no nonlinear convective terms arise in this formulation. A pure

Lagrangian formulation in Fluid Mechanics is seldom used, since the trajectories of

the fluid particles are very complex and the numerical meshes can get very distorted
and entangled.

In the spatial or Eulerian formulation, the motion is followed in terms of the

spatial coordinates x. This amounts to a frame of reference that is fixed or moving

with a predefined velocity. The spatial or Eulerian velocity is defined as:

(3.5)

The Eulerian acceleration can be shown to be of the following form:

a(x, t)
a
=at u(,P(X, t), t)ix=q,-l(x,t) =
au
at (x, t) + u(x, t) · Vu(x, t) (3.6)

Note that in the Eulerian description of flow, nonlinear advection terms appear.

The Eulerian description of flow is widely used in the Fluid Mechanics studies.

The nonlinear advective terms u(x, t) · Vu(x, t) are the main cause for most of the

difficulties in the Fluid Mechanics studies.

In the ALE formulation, a third domain called the referential domain R is defined

as shown in Fig.3.1. At timet= 0, it occupies an open region Ro. The referential

body has its own independent motion and coincides with the material body Bt at

some later time t > 0, i.e. R 1 =B 1• The referential coordinates are written as

:X= (:X~, x2). The referential point q arrives at the spatial point pat timet through

the referential motion ,\ defined as:

x=>-(x,t). (3.7)
 17

The velocity oi the referential motion w(x. t) is given by:

8
~,\(x.t)j-_
w(x.t) = ut X-. 1
-1 1 X. tl· (3.8)

All the quantities defined on the referential domain are denoted with a superscript
hat. For example.

w ( x. t) =w ( ,\ ( x. t) , t) = w (x. t) . (3.9)

The relative motion oi the Ruid body with respect to the referential motion can be

expressed as:

x = A- 1 ( x, t) = ,\- 1 ( ¢> (X. t), t) = w (X, t) . (3.10)

The relative velocity of the fluid motion with respect to the referential motion can

be written as:

(3.11)

The material point P and referential point q arrive at the spatial point p through

independent motions, which are related as follows:

x= A(x,t) = A(l/>(X,t),t) = 4>(X,t). (3.12)

From the above relation, we can derive the following:

a
u(x, t) = u(x, t) = ot A (1/> (X, t), t) lx=>i'-1(x,t)
~~ (x, t) + F (x, t) y>i' (x, t) (3.13)

w (x, t) + F (x, t) y>i' ( x, t)

The deformation gradient F( x, t) is defined as:

 18

' x, t) = ax·
F;; ( ;::, ' I (3.14)
ux· J

The spatial acceleration in the referential coordinates can be shown to be:

a(x, t) = a(x, t) = :t u(x, t)lx=,;,-l(x.t)

(3.1.5)
:t u(x., t) + (v.P(x, t). v) u(x, t)

From Eq.3.13, V.P is replaced in the above equation resulting in:

a(x, t) = ~~ (x, t) + [(:F- 1 (u- wJ). v] u(x, t) (3.16)

The above expression, represents the fluid acceleration in the referential kinematic
-
description. Comparing the spatial acceleration (Eq.3.6) with the referential accel-

eration (Eq.3.16), the difference is that gradients are with respect to the referential
coordinates and the spatial velocity u(x, t) is replaced with the relative velocity fr.P.

The conservation of momentum equation in the spatial description is given by:

p(x, t)(
au
at (x, t) + u(x, t) · Vu(x, t)) = p(x, t)f(x, t) + \7 · o-(x, t), (3.17)

where p(x, t) is the local fluid density, f(x, t) is the body force, o-(x, t) is the Cauchy

stress tensor, and V(x, t) is the divergence operator. In the referential kinematic

description, the above momentum conservation equation becomes:

p(x, t)Jf(x, t)
(3.18)

where J detF, and P Ja-:F- 1 is the Piola-Kirchoff tensor of the first kind.

Similarly, the conservation of mass equation in referential description is given by:

 19

;)p
()t(x.tl+ [('F- 1 (u-w)
. '1
) ·V'j/J(x.tl+fJ(x.t) l
[("F- 1 v") ·u(x.t) =0.(3.19)

In the case of incompressible this reduces to:

(.F- 9). u(x, tJ =

1
o (3.20)

In the case of free boundary problems. the referential motion is related to the

fluid motion. For example. at the free boundary, the mesh points are moved normal

to the free boundary with the fluid velocity. to prevent loss or gain of fluid material.

If a discrete time-stepping procedure is used to solve the momentum and mass

conservation equations, the previous time step fluid particle velocities can be used

to approximate the referential motion. In that case, the time-stepping procedu!:'e

is only of first order accuracy in time and the governing equations in the ALE

formulation can be written as:

ou . ( x Ou ( ) Ou . l Op
-+ u-w )-+ v-wY -=gsznfJ---+v'V2 u (3.21)
at ox oy p ox

ov x Ov Y Ov 1 Op 2
-, +(u-w )-+(v-w )-=-gcosfJ---+V'v (3.22)
at ox oy p oy

ou av - 0 (3.23)
ox+ ay -

Note that if there is no referential motion, the mesh velocities wx = 0 and wY = 0,

and the ALE formulation reduces to the Eulerian or spatial formulation. If the
referential motion coincides with the fluid motion, then wx = u, wY = v, and the

ALE formulation reduces to the Lagrangian or material formulation. The relation

between the above described three formulations can be expressed as:

D 0 0

- = ot + u . v = 8t + (u - w) . v . (3.24)
Dt .___.,
Lagrangian Eulerian ALE
 20

3.3 The concept of vertical spines

The most important thing in the ALE approach 1s the choice of the referential

motion. The aim is to chose a referential motion, so that the free boundary is

followed accurately without severe mesh distortion. The type of problems we are

interested in, provide us with a simple referential motion, that is very easy to im-

plement. Physical domain generally consists of an inlet or upstream end where the

flow enters, an exit or downstream end through which the fluid leaves, a solid-liquid

interface at the bottom and a liquid-gas interface at the top (Fig.l.l). For these type

of problems it is advantageous to let the grid points move in the vertical direction

only, but restrict their motion in the horizontal direction. For this purpose, the free

surface is parameterized using vertical spines as shown in Fig.3.2. The free surface,

and consequently the numerical mesh points, slide up and down along the vertical

spines. The numerical grid points are however, fixed in the horizontal direction

(wx = 0). This type of a referential motion, lets us follow the free surface movement

accurately, without significant mesh distortion.

 21

free surface

vertical spine

Figure 3.2: Free surface resting on vertical spines

 22

Chapter 4

Numerical Modeling of Turbulence

4.1 Introduction

When the nondimensional Reynolds number Re, which represents the ratio of in-
ertial forces to the viscous forces, is over a certain cut-off value Rec, the nature of

the fluid flow undergoes a drastic change, going from a well defined simple laminar

flow to a highly complex irregular and random turbulent flow. Hinze(1975) defin~s
turbulence as

Turbulent fluid motion is an irregular condition of flow in which the vari-

ous quantities show a random variation with time and space coordinates,
so that statistically distinct average values can be discerned.

The most important thing to note is that even though the fluid flow may appear to

be irregular and random, the statistical quantities such as mean and variance retain

physical relevance. Among the applications considered in this dissertation, the open

channel flows, which refers to gravity-driven flows occurring in nature such as rivers,

streams, channels etc., are usually turbulent in nature. On the other hand, thin-film

flows, which refer to very thin liquid films draining down vertical or inclined walls,
are mostly laminar in nature.

It is believed that the 3-D Navier-Stokes equations can be used to simulate the
turbulent fluid flows. The biggest challenge to numerical modeling of turbulent fluid

flows comes from the existence of widely varying length scales requiring a very fine

mesh to resolve all the scales of motion. A turbulent fluid flow requires Re 2 mesh
 23

points in 2-D and Re 9 14 mesh points in 3-D, where Re is the Revnolds number. Even

for a fluid flow with Reynolds number 10 4 , we will need 108 mesh points in 2-D itself

and this is far beyond the capacity of the present day computers. This leaves us with

no other recourse but to model the turbulence through some approximate models.

The most common approach is the Reynolds Averaging, which is described briefly

in the next section.

4.2 Reynolds Averaged Equations

The Reynolds averaging procedure filters out the turbulent fluctuations from the

mean flow and the resulting equations are called the Reynolds Averaged N a vier-

Stokes equations (RANS). The fluid velocities and pressure are expressed as tl:te

sum of the mean (in an appropriate sense) and a turbulent fluctuation as follows:

u·
'
(4.1)
p p+p'

The mean quantities are denoted by an over bar and the fluctuating components

are denoted by the superscript "'. The averaging employed could be either time

averaging, or space averaging or ensemble averaging. All the averaging are identical

if the ergodicity hypothesis is assumed. The Reynolds averaged continuity and

momentum equations are:

(4.2)

(4.3)

a;j is the stress tensor given by:

 24

( 4.4)

11 is the dynamic viscosity of the fluid, u;,j is the velocity gradient of the fluid flow,

and D;j is the Kronecker delta. Due to the nonlinear advection terms the fluctuating

components do not disappear with Reynolds averaging but reappear as uiuj which

are called the apparent stresses or Reynolds stresses. These new quantities (six

of them) are to be determined from the mean flow field and presents us with the

classic closure problem. These terms are modeled using physical and dimensional

arguments and the unknown coefficients are tuned using experimental data.

4.3 Eddy viscosity hypothesis

As long as the continuum assumption is not validated, the Navier-Stokes equations

can be assumed to model the physics of the fluid flow accurately. Due to the

Reynolds averaging, however, information is lost irretrievably. Through closure

models we try to put at least part of the information back into the mathematical

model. No amount of complex and accurate modeling, however, can make the RANS

equations as accurate as the original Navier-Stokes equations. This i.s the penalty

we pay for resorting to averaging.

Turbulence models can be broadly classified into Boussinesq models, and non-

Boussinesq models. In the Boussinesq models, the Boussinesq Eddy Viscosity hy-

pothesis is invoked and the Reynolds stresses are written as:

(4.5)

where,

1-,-, (4.6)
k = 2U;U;,
 25

is the turbulence kinetic energy. Note that. the expression for the Reynolds stresses

is analogous to the original stress tensor proposed by Stokes. The term 2k/3 is

added to the Eq.4.5 to make it valid if a contraction is applied. The term 2k/3 can

be considered as turbulent pressure. As opposed to the laminar viscosity (f.l ), the

turbulent eddy viscosity (f.lr) is a property of the flow, and not a property of the
fluid substance. Originally, Boussinesq assumed the eddy viscosity to be a constant

scalar (Hinze 1975, pp.23). However, a constant scalar turbulent viscosity is not true
for most realistic fluid flows. The turbulent eddy viscosity varies with space and

time and is a function of the fluid flow field. It is known that, turbulence leads to

intense mixing with the attendant increase in momentum diffusion. The Boussinesq

Eddy Viscosity hypothesis, thus tries to model the turbulence as an increase in the

effective fluid viscosity or mixing coefficient. With the Boussinesq eddy viscosity

hypothesis, the stress tensor <Yij is of the form:

(4.7)

Hinze (1975) has a discussion about the various admissible forms for the Reynolds

stress -pu;uj. In the above equation a scalar eddy viscosity is assumed. However,
we could assume second-order eddy viscosity of the form:

(4.8)

Or else, we could also assume a fourth-order eddy viscosity of the form:

(4.9)

However, the above given second-order and fourth order eddy viscosity forms do not

give us a stress tensor that is invariant under all rigid body motions.
 26

An eddy viscosity hypothesis does not give us a turbulence model complete in all

respects. The other alternative is to solve for the Reynolds stress -pu'u' directlv.
' J "

This usually means, we need to solve six additional partial differential equations

for the six independent Reynolds stress components. However, even these models

are not completely rigorous, as an effort to derive equations for the second order

Reynolds stress terms results in third-order correlations of the u~uju~, which need

to be modeled. This is the classical closure problem in turbulence, with the net

result that we have more equations than unknowns. In spite of their drawbacks,

the eddy viscosity models are the most commonly used turbulence models. The

turbulence model used by us also makes the eddy viscosity hypothesis. Thus, we

model turbulence by solving the RANS along with the Boussinesq eddy viscosity

hypothesis.

4.4 Eddy viscosity models

Once we invoke the Boussinesq eddy viscosity hypothesis, the whole turbulence

modeling reduces to finding the eddy viscosity (J-tT ). The turbulence closure models

which solve the RANS along with the Boussinesq eddy hypothesis can further be

classified into:

1. Zero-equation models

2. One-equation models

3. Two-equation models

The above classification is based on the number of additional transport equations

that are to be solved in addition to the RANS to compute the eddy viscosity J-!T·
 27

The most popular of the zero-equation models is the Prandtl's mixing length
hypothesis and relates the eddy viscosity to the mean velocity gradient as follows:

(4.10)

The above relation involves a single unknown parameter, namely, the mixing length

lm, which is specified empirically. For example, in the case of boundary layer flows,
the mixing length is assumed to vary linearly with the distance from the wall.

In the above equation, vy = flT / p is the kinematic eddy viscosity. The mixing
length models are very easy to incorporate and do not increase the computational

time excessively, as no additional equations need to be solved. However, they are

employed in very simple shear and boundary layer flows, where it is easy to specify

the mixing length. For complex flows such as separating boundary layer flows, the

specification of the length scale is not straightforward and the mixing length models

are rarely used in such complex flows. The turbulent fluid flows to be simulated

in this dissertation are complex, involving boundary layer separations and both the

wall generated turbulence mechanisms and the mixing layer turbulence mechanisms

are present simultaneously. Thus, the mixing length models will be expected to

perform poorly in this instances and are not considered.

A simple dimensional arguments, will show us that the eddy viscosity is influ-

enced by a turbulent length scale and a turbulent time scale. Thus a means of

evaluating these two turbulence scales are required. The one equation models at-

tempt to compute the velocity scale through a advection-diffusion type transport

equation and specify the length scale empirically. However, since the length scale

is specified empirically, these models in general do not perform much better than

the mixing length models. The most general turbulence model requires the solution
 28

of two partial differential equations, to determine both the time and length scales
accurately, without having to specify any of them empirically.

There are several two-equation models that have been proposed to model tur-

bulence. Some of them are k- kl, k- w and k- E models. They differ from each

other in the parameters that are solved for directly. The k - l model solves for the
turbulence kinetic energy k and the turbulence length scale l. In k - w model, the

turbulent vorticity scale is obtained in addition to the k. Whereas, in the k - .:

model, the dissipation rate of the turbulence kinetic energy rate is computed in ad-

dition to k. No matter what we solve for, the idea is to put them together through

dimensional arguments to compute the eddy viscosity. In the context of hydrody-

namic flows, the two-equation k - E model has been shown to give reasonably good

results. According to the ASCE Task Committee (1988) appointed to look into the

use of turbulence models for simulation of hydraulic flows, k - E model is the one

being most widely used. In this dissertation, the two-equation k - E model is used

which is described briefly in the next section.

4.5 k- E turbulence model

k denotes the turbulence kinetic energy k and.: is rate of dissipation of the turbulence
kinetic energy. These two variable are defined in terms of the fluctuating components

as follows:

k= ~u 2 (4.11)
2 "

(4.12)
 :29

The k- E used by us is the one originally proposed by Launder & Spalding (1974).

The turbulence quantities k and E are determined through the solution of advection-
diffusion type transport equations given by:

(4.13)

(4.14)

After finding k and e, the eddy viscosity vr is computed from the following relation:

vr = cJ1.e;e (4.15)

The constants in these equations are (Launder & Spalding 1974):

cJ1. = o.o9, c1 = L44, c2 = 1.92, (Jk = 1.0, (J< = 1.3

Note that, G = vr (ii;,k + iik,i) ii;,k in the above equations represent the generation
of turbulence kinetic energy. The turbulence field derives its energy from the mean

flow which is then cascaded into the smallest eddies where it is dissipated into heat

through viscous dissipation. TheE terms represents this dissipation of the turbulence

kinetic energy. It should be noted that, the k- E model presented in this section has
been originally derived for pipe flows. The open channel flows, are in some respects

similar to the pipe flows. They are bounded below by a rigid wall and above by a

free surface, which is in many respects similar to the centerline axis of symmetry

in pipe flows. However, it has been speculated for some years that the presence of

the free surface dampens the turbulence since the free surface also acts like a wall,

albeit friction-less. However, no satisfactory and well-accepted turbulence models

 30

exist for open-channel flows even though Celik (1980), among many others tried to

modify the k - E model to be suitable for open-channel flows.

4.6 Boundary Conditions for k - E Model

The turbulence model as presented above is a high Reynolds number model and is

valid only in the fully turbulent regime. However, close to the wall there exists a

small layer where the viscous effects dominate over the turbulence effects (viscous

sublayer). The k - E turbulence model presented previously is not valid in this

region. Modifications have been proposed to the original k - E turbulence, so that

they remain valid for low Reynolds number regions also (Launder & Spalding 197 4).

However, the more popular approach has been to use the original k- E model itself,

and specify the boundary conditions at the wall using the wall function approach.

In this method, the first computational point is not placed right on the wall, but is

kept at a small distance away from the wall so that it lies in the fully turbulent log

layer as shown in Fig.4.1.

In the small layer close to the wall, the turbulence mechanism is assumed to be

in equilibrium. That is, the generation of turbulence kinetic energy is equal to the

dissipation of the turbulence kinetic energy, G = E, and there is no advection or

diffusion of the turbulence kinetic energy. In this limit, the k - E model reduces to

the mixing length hypothesis, and the universal log-law is assumed to be valid close

to the wall. The boundary conditions for k and E are imposed as follows:

(4.16)

(4.17)
 :n

( (I) ( (

( ( ( ( ( (

( Outer Layer
( ( (~ (
(Wake Layer)

000000 0 0 0 0 ooooooo 0 0 0 0 ooooooo 0( 0 Ooooooo 0 0 0 ~ooooo 0 0 0 0 _! 00 0 0 0 0 00000.1 0 0 0000000 0 0 0 0 .( ~- 0 0 0 0 oooo- ..( ooooooo 0 0 0 0 ooooooo 0 0 0 0 0000•

) ( 0 ()
Log Layer
) < N > )

) ) 00-r~ p ) ) ()
0000000 00 ooooooooo 00 oooooooo 00 0 oooooooo 0 oooooooYJ 0 0 ooooo 00 000000 00 00 00000000 00 0 ,.000000 00 000000000 0 oooooooo; 0 .. ooooooo 00 0 00000

P Viscous
t Sublayer

Figure 4ol: The near-wall node placement.

 :32

where, Yf = ypu.jv, and yp is the normal distance from the wall. The boundary
conditions for u and v at the wall are imposed as:

u.yp+ if Yf < 11.6

UF{ -u. 1n (Eyf) if Yf;:::: 11.6
(4.18)
K,

Vp = 0 ( 4.19)

Specifying the normal velocity to be zero violates the continuity requirement, but

since the wall layer is very thin, the error should be negligible. E is a function

of wall roughness and is equal to 9.0 for smooth walls. Apart from convenience,

the wall function approach is very advantageous computationally, since it bypo.sses

the high gradient region close to the wall. Moreover, information pertaining to

factors such as wall roughness, pressure gradient, curvature etc. can be introduced

through these wall functions. To determine the friction velocity u. = ;;:JP, the
approach proposed by Benim & Zinser (1985). It is assumed that UN lies on the

same logarithmic curve as up, and the friction velocity is computed as (see Fig.4.1):

if yj;, < 11.6

(4.20)
if yj;,;:::: 11.6

At the free surface the normal gradients of k and E are taken to be zero (plane

of symmetry).

ok
on ( 4.21)
OE = :} aty=h(x,t)
on
Boundary conditions at the inlet and exit are problem specific. In this disser-

tation we use periodic conditions in our modeling of the turbulent fluid flow over
 33

artificial bedforms (chapter 7). In the hydraulic jump simulation (chapter 8), the
following boundary conditions k and t: are imposed at the inlet:

10(1- y) ~
if undeveloped inlet flow

k(O,y,t) = (4.22)
Y1 ;c: if fully developed inlet flow

t:(O,y,t) = 10( )[
3
] if undeveloped inlet flow ( . )
4 23
l - .!!_ u. if fully developed inlet flow
Yl KY ( l - e-y+ /26)

The boundary conditions at the inlet imposed in the case of fully developed flow

are similar to that proposed by Alfrink & van Rijn ( 1983). In the above equations,
-
Y1 stands for the height of the free surface at the inlet. In the hydraulic jump

simulation, fully developed conditions are imposed at the exit as follows:

ak
ox (x = L, y, t) = 0 (4.24)

ot:
ox(x=L,y,t)=O (4.25)

4. 7 Surface Roughness

The surface of the bottom wall is usually very irregular and rough. In turbulent
flows it is common to represent the roughness of the bed in terms of the equivalent

sand roughness height k •. The influence of the wall roughness on laminar flows is

minimal. However, in turbulent flows, they can influence the fluid flow considerably,

by breaking up the laminar sublayer (White 1991). The relative roughness of the
wall can be expressed in terms of the parameter:

(4.26)
 31

where, vis the laminar viscosity and u. is the friction velocity. White ( 1991) defines
the following three roughness regimes:

1. hydraulically smooth wall: k; <4

2. transitional-roughness regime: 4 < k; < 60

3. fully rough flow: k; > 60.

If the roughness elements are very small and are submerged within the laminar

sublayer, then the flow is considered to be hydraulically smooth and roughness has

no effect on the velocity profile. The nature of the flow is mostly influenced by

the fluid viscosity in the hydraulically smooth wall regime. In the fully rough flo_w

regime, the laminar sublayer in broken by the roughness elements, and the fluid

flow at the wall is determined by the roughness elements and the viscous effects are

negligible in comparison. For the intermediate transitional regime, both viscosity

and roughness effects influence the fluid flow.

In the numerical model we incorporate the influence of the wall roughness using

the expressions provided by Rosten & Worrell (1988). In this model the functional

dependence of E on the wall roughness is given by:

if k: < 3.7
3. 7 < k: < 100
1
E= 112
if ( 4.27)
[a (k; /b) + (1- a)/ E!]
2

b/k: if k; > 100

In the above expression, Em 9, b 29.7, a = (1 + 2X3 - 3X 2 ), and X
 35

4.8 Note on the RANS equations

For the sake of completeness, we list below the RANS equations used by us:

au
-.-
au au 8p 1 . + (v +liT) "V u 2

at + uax- + u8y- =
pax - - - - gsmO (4.28)

av au + au = - -1 -ap. + g
-at + U- + (!I + liT) \J 2
( 4.29)
ax ay U-
pay COS () U

au au
-+-=0 (4.30)
ax ay
In the above equations, u and u are the mean fluid flow velocities. p = p + 2k./3
is the pressure, and the turbulent pressure is absorbed into the original mean pres-

sure. Note that, the RANS as stated above simplify the original stress tensor given

in Eq.4. 7. The eddy viscosity is a spatially varying quantity, and the divergence

of the stress tensor is strictly not equal to Laplacian of the velocity. However, to

retain the full stress tensor would necessitate the solution of the x-momentum and

y-momentum equations in a coupled fashion. Moreover, Lemos (1992b) found nu-

merical oscillations when he tried to retain the full stress tensor form in free surface

problems.
 36

Chapter 5

Numerical Procedure

5.1 Introduction

In the previous two chapters, the two most importa.nt elements of the numerical

procedure, namely, the ALE formulation and the k- E turbulence model, have been

described. The rest of the numerical details are briefly illustrated in this chapter.

When trying to solve the Navier~Stokes equations, one is faced with the qu~s­

tion of which primary variables to chose. The mathematical system (Eqs.2.1-2.3) as

presented in chapter 2 has u - v - p as the primary variables. There exist alternate

formulations, which employ stream function - vorticity ('ljJ- w ), or velocity- vortic-

ity as the primary variables. The biggest problem with solving the incompressible

Navier~Stokes equations comes from the pressure. Since, no time derivative for pres-

sure appears in the formulation, time-stepping for pressure is not straightforward,

as it is for velocities. The vorticity formulation, eliminates pressure by taking curl

of the momentum equations. However, in free surface problems it is advantageous

to have u - v - p (also called primitive variables) as the primary variables. This

is because, at the free surface we have a boundary condition in terms of pressure,

namely, the normal stress condition (Eq.2.7). In recent years the u- v- p method is

also gaining popularity, due to decoupling introduced by the projection schemes. In

these projection schemes, the velocity components and pressure are decoupled from

each other and can be solved in a sequential manner, thus providing significant sav-

ings in computational time and memory. In this dissertation, we use the primitive
 :n
\·ariables I, t t - 1; - p) as the primary variables and solve them using the Chorin-type

projection scheme which is de:,;cribed brietiy in the next section.

5.2 Time stepping procedure

Knowing the previous time step velocities u". v", the next time-level velocities un+I.
,,n+J are obtained as follows:

Stepl The intermediate w·locities ii.n+l. £,n+l are obtained by solving the momentum

equations dropping the pre,sure gradient terms.

l·ln+l - tl" 0 n ;:) n

- - - - - 1/ o2 ii."'"l = Jn- ( un- w") _-_u_- (vn- wY) _u_u_ (5.1)
::,.t V X ' iJx oy

vn+! - v" j) n j) n
- - . , . . . - - - - 1/ \72 z:,n+l = Jn- (un- wx) _;!-- (vn- wY) ~ (5.2)
D.t y ox fJy

J; and f; are the body forces in the x and y directions respectively. wx, wY

are the mesh velocities or the referential velocities which come from the ALE

formulation described in chapter 3. The appropriate Dirichlet and Neumann

boundary conditions are imposed in this step. Note that the diffusions terms

are treated implicitly, whereas the nonlinear advection terms are computed

explicitly, making our numerical procedure semi-implicit in time.

Step 2 The pressure is determined from the Poisson equation which is obtained

by projecting the intermediate velocities onto the divergence free spaces. The

boundary conditions for pressure are the homogeneous Neumann condition

P,n = 0 on all boundaries except the free surface where a Dirichlet bound-

ary condition determined from the normal stress continuity is imposed. If

 38

the domain is periodic in the streamwise direction, then periodic boundary

conditions are imposed for pressure.

1 2 1 (aun+I 8vn+I)
-v
P
p=- - - + - -
6t 8x 8y
(5.3)

Step 3 The final velocities un+I, vn+I are calculated from:

(5.4)
6t p8x

vn+l - vn+l 1 8p
--- (5.5)
6t p8y

The philosophy of the Fractional Step Method is to first determine the approx-

imate velocities un+l' vn+l satisfying the boundary conditions but which are not

necessarily divergence-free. Theses approximate velocities are then projected onto

divergence free spaces resulting in final divergence free velocities un+l, vn+l. This

method of discretizing Navier-Stokes equations in time is widely used (e.g., Chorin

1967; Temam 1971; Patankar 1980; Pironneau 1982; Donea, Giuliani & Laval 1982;
Mizukami & Tsuchiya 1984; Kim & Moin 1985; Ramaswamy & Kawahara 1987b;

Van Kan 1986; Glowinski 1986; Bell, Colella & Glaz 1989; Gresho 1990; Gresho &
Chan 1990; Le & Moin 1991; Oden 1992; Finlayson 1992; Ramaswamy, Jue & Akin

1992).

5.3 Spatial discretization

The spatial discretization of the fluid domain is done using three-node triangular

elements. In two-dimensions, the simplest element we could have are the triangular
 39

elements, hence thev are also called as simplex elements. The interpolation functions

are linear in these <"lernents, and the \·ariables are taken to be residing at the vertices

of the element. E4ual-order approximation, that is the same level of approximation,

which in our case is linear approximation. is used for both velocities and pressure.

The equal order approximation (also called non-staggered grids) for the velocity and

pressure has been used by several researchers in the numerical solution of the Navier-

Stokes equations (e.g. Schneider, Raithby & Yovanovich 1978; Rice & Schnipke 1986;

Ramaswamy & Kawahara 1987c: Shaw 1991: Zienkiewicz & Wu 1991: Zienkiewicz

& 'vVu 1992: Behr. Franca ,\c Tezcluyar 1992: Gresho. Chan. Christon & Hindmarsh

1994).

The velocities and pressure are approximated by the linear interpolation function

<I>" in each local triangular element as:

(5.6)

(5.7)

In the conventional Galerkin-Bubnov finite element method, the weighting functions

are same as the interpolating functions. In linear triangular elements, the spatial

gradients are constant, and the three-node triangular elements are also called as

Constant Strain Triangles in the Solid Mechanics literature.

Let us define the following matrices:

(5.8)

(5.9)

(5.10)
 40

(5.11)

Ca{Ji = inn <I>a<I>{J,i df/. (5.12)

'

(.5.13)

Mcx{J is the mass matrix, Da{J is the diffusion matrix, Aa{J,j is the advection matrix,

La{J is the Laplacian matrix, Ca{Ji is the gradient matrix, and Fai is the body force

vector. Note that the transpose of the gradient matrix C'j;{Ji is the divergence matrix.

The subscripts a and f3 refer to the local node numbers, and take the values 1, 2, 3,
since we are using three-node triangular elements. D.~ is the domain of the local

triangular element, with the subscript e referring to the fact that it is an element

level quantity, and the superscript n refers to the time level of the domain. In

addition to the above matrices, we define a lumped mass matrix M~{J which is

obtained by lumping the off-diagonal elements of the consistent mass matrix Mcx{J

on to the diagonal.

The Fractional Step Method in the Finite Element formulation become:

Step 1

Step 2

(.5.1.5)

Step 3
 41

(5.16)

:\ ote that in this step. the lumped mass matrix is used.

The above element level equations are assembled onto a global matrix, and solved
for the global velocity and pressure field.

5.4 Mesh generation

The first step in the mesh generation process is the placement of the vertical spines

along the streamwise direction. The placement of the vertical spines is dictated

by the physics of the fluid flow being investigated. Usually we require fine spacin_g

for the vertical spines in high gradient regions. In the turbulent flow over artificial

bedforms simulation, the vertical spines are close-packed near the peak of the dune,

since very high gradients exist in this region and the boundary layer separates just

downstream of the dune. In the simulation of hydraulic jump, we find that the free

surface rises abruptly across the jump region. Hence large spatial gradients occur in

the jump region, and fine spacing is required in the jump region. Unfortunately, we

do not a priori know the location of the jump, and hence place the vertical spines

in an equi-distant manner and try to capture the high gradient regions accurately

by taking sufficiently large number of spines. In the case of thin film flows, a

uniform distribution of vertical spines is sufficient and in a sense appropriate since

the traveling wave phenomena is being investigated.

A general mapping strategy described by Anderson, Tannehill & Pletcher (1984)

is used to achieve arbitrary spine placement.

x = B+ ~sinh- 1 [(:c -1) sinh(rB)], 0< r < oo (5.17)

 42

where

B=-lln [l+(eT-l)(xc/L)]
2T 1 + (e-T -l)(xc/L) .

Uniform spacing in the transformed space i produces unequal spacing in the x plane,

depending on the values of Xc and T. Xc corresponds to the point about which we

desire fine mesh spacing. In the case of turbulent flow over dunes, Xc is chosen to be

the dune crest. The fine spacing of the spines at the point Xc is controlled through

the T parameter. Large values of T would give finer mesh spacing near the point

Xc. In the turbulent flow simulations to be presented in chapter 7 T = 5.0 has been

used.

After the placement of the vertical spines, it remains to assign nodes to each

of the spines. Each of the spines are assigned a fixed number of nodes, and their

placement along the spines is specified a priori. In the case of thin-film flows, the

fluid flow is laminar, it suffices to place the nodes in an equi-distant manner along

the vertical spines. However, in turbulent flow simulations, fine mesh spacing is

desired near the wall to resolve the large gradients known to exist close to the

wall. A logarithmic mapping is obtained using the mapping strategy proposed by

Anderson, Tannehill & Pletcher (1984).

- - 1-
[n {[/3 + 1 - y-b(x)
h(x,t)-b(x)
l / [/3 -
1+ y-b(x)
h(x,t)-b(x)
l} ' 1< f3 < (X)
(5.18)
y- ln [((3 + 1)/(/3 -1)]

An equal spacmg m the f) space results in fine mesh spacing near the bed b( x)

depending on the value of (3. For large values of f3 we obtain uniform spacing,
whereas for values of f3 --.. 1, fine mesh spacing near the wall is obtained. In the
turbulent flow simulation reported in this dissertation, the values of f3 are chosen in
the range 1.001 - 1.01.
 To summanze, a uniform spacmg m the transformed plane ( i, fj) results in

nonuniform mesh spacing in the physical plane ( .r, y) depending on the values of

5.5 Miscellaneous details

In the turbulent open-channel flow simulation, in addition to the RANS equations,

advection-diffusion transport equations are to be solved for k and c. These are very

similar to the Eqs.5.1&.5.2 which solve for the intermediate velocities un+l, vn+l. The

diffusion term is treated implicitly using the Euler Backward scheme and advection

terms are treated explicitly using the Euler Forward scheme. The only difference is

the presence of highly nonlinear forcing functions in the case of k and E equations.

A simple explicit treatment of these highly nonlinear terms was found to result in
numerical oscillations, Therefore the following time discretization for the k and E

equations is used:

(5.19)

(5.20)

The explicit first-order treatment of advection in Eqs.(5.1,5.2,5.19,5.20) results

in negative diffusion being added to the system (Dukowicz & Ramshaw 1979). This

is removed using the Balancing Tensor Diffusivity procedure. In this procedure the

viscosity tensor is modified to be:

(5.21)
 44

The variable ¢ represents u, v, k or E, and V¢, the corresponding physical viscosity.

According to Rannacher (1992), the Chrorin's projection method introduces

pressure stabilization in the numerical solution of the Navier-Stokes equations and

is similar to the Petrov-Galerkin procedure of Hughes, Franca & Balestra (1986).

Unlike the solution of the Navier-Stokes equations, the k and E equations have

no implicit stabilization. A procedure called Selective Mass Lumping is used to

introduce just enough numerical diffusion to stabilize the scheme. The idea is to do

mass lumping to different extents on the left- and right- hand sides of the equations.

The stabilization introduced by selective mass lumping is in some respect similar

to the Lax-Friedrichs scheme in Finite Difference Method (Kawahara, Hirano &

Tsubota 1982; Chippada, Ramaswamy & Wheeler 1994).

Sherr (1992) has been able to show that un+I and un+l are weakly first order ap-

proximations to u(tn+I) and that pn+I is weakly order 1/2 approximation to p(tn+d-

The free surface is updated through the solution of the kinematic free surface

equation. An alternative formulation of the kinematic free surface equation can be

derived:

8h 8q- 0 (5.22)
8t +ox-

where, q(x, t) = Jb~S)' ) is the flow rate. The above equation can be considered to
1

be the conservative form of the original kinematic free surface equation. In the

problems with periodic boundary conditions, the above equation is solved using the

Fourier Spectral Method along with leap-frog time-stepping procedure. However,

in the simulation of turbulent flow over artificial bedforms, even though we have

periodic boundary conditions, we cannot use the Fourier Spectral Method since the

mesh spacing is not uniform. Second-order Adams-Bashforth scheme is used in such

 45

instances. Second-order Adams-Bashforth scheme is also used in the non-periodic

domain simulations to be reported in the chapter on interfacial instabilities in thin-

film flows (chapter 6). In the simulation of hydraulic jump, very large gradients

exist in the jump region. since the free surface rises abruptly. Certain amount of

numerical diffusion is needed to control the numerical oscillations, and this is done

using the Selective :Vlass Lumping procedure (Chippada eta! 1994).

Since, advection terms are treated implicitly, the numerical procedure is only

conditionally stable and we need to satisfy the Courant-Friedrichs-Levy (CFL) con-

dition. The time step is chosen automatically based on the stability condition in

the following manner:

(5.23)

Thus, the time step is chosen such a way that the Courant number in each of

the elements is less than one. Even though theoretically, we could take a Courant

number of one and still have stability, this result is strictly valid only for linear

problem only. The admissible Courant number is usually less than one for nonlinear

problems, and the parameter C F Lin the above equation lets us specify the Courant

number. The C F L number has been specified in the range 0.05 - 0.40, depending

on the problem. Note that, in the above expression used to find the maximum

admissible time step, ..jgh cos(} is the celerity of a gravity wave and J /ph is the
G"

celerity of a capillary wave.

The system of linear equations are solved using the Conjugate Gradient Method

using the tridiagonal preconditioner.

 46

Chapter 6

Interfacial Instabilities in Thin Film Flows

6.1 Introduction

Thin liquid films flowing down inclined walls, pipes and tubes occur widely in many

engineering applications, e.g. film cooling, gas adsorption, condensers, evaporators,

slot coating and dip painting. Thin film flows are unstable to long wavelength

disturbances, and the gas-liquid interface tends to become corrugated, resulting in

various surface wave motions. Traveling waves on the liquid-gas interface lead to

increased heat and mass transfer both in the liquid and gaseous phases (Dukler

1972). In film coating applications, however, a constant thickness liquid film is

desired. Thus, due to its obvious engineering importance, thin liquid film flows have

been extensively studied since the seminal work of Kapitza & Kapitza (1949).

The thin film flow study can be classified into the following subgroups (Dukler

1972):

1. Study of the conditions under which waves exist. These studies try to explain

why we observe waves on the interface, and if there is more than one type of

wave, the ranges of parameters in which these different types of waves exist.

2. Study of the quantitative aspects of the surface waves, such as, the wave shape,

the wave speed and other statistical features.

3. Study of the influence of these surface waves on the heat, mass and momentum

transfer, both in the liquid and the gaseous phases.

 47

The focus of our work is on the first two aspects, namely, the qualitative and quan-
titative study of the formation of surface waves on the gas-liquid interface.

The state of the art related to the thin film flow modeling has been reviewed

recently by Chang ( 1994). He identifies four distinct wave regimes as shown in

Fig.6.1:

Region I This is the wave inception region. Disturbances at the inlet are amplified

downstream and form a monochromatic wave at the end of the region. If the

flow entering is not artificially excited (natural noise) or if it is excited with

a broad band noise, the disturbance is selectively amplified and at the end of

inception region, monochromatic waves with frequency same as the maximum

growth rate frequency predicted by the linear stability analysis are observed.

However, if a single frequency forcing with small amplitude is applied to the

incoming flow, then waves with the inlet frequency are amplified and at the

end of region I monochromatic waves with frequency same as the inlet fre-

quency are observed. The amplitude of the disturbance grows exponentially

in this region. The waves are monochromatic (sinusoidal) and infinitesimal in

amplitude and linear theories can be used to study this wave regime.

Region II In this region, however, the exponential growth rate is arrested due to
nonlinear interactions and the amplitude of the wave saturates. Part of the

energy of the fundamental mode gets channeled into its first harmonic due to

the nonlinear interaction, and the wave shape is no longer sinusoidal.

Region III The primary surface wave instability exiting Region II travels for some
distance downstream without visible change in wave shape, before succumbing

to secondary side-band and subharmonic instabilities. Due to the subharmonic

 48

instability, neighboring peaks coalesce, and due to the sideband instability

longwave modulation of the wave amplitude takes place. These instabilities

lead to very localized waveforms, namely, the tear drop humps. These tear

drop humps are also referred to as solitary waves or solitary humps and have a

very long gently trailing tail and a steep front preceded downstream by small
capillary ripples. The Fourier content of this solitary wave is broad banded,

with the fundamental and several of its harmonics possessing significant en-
ergy. The spacing between the humps is nonuniform and time varying even

though the shape of the solitary wave itself may be nonvarying in time.

Region IV The solitary waves evolving from the region III finally succumb to
-
transverse disturbances, and the subsequent waves are three-dimensional, ir-
regular and random.

It should be pointed out, that the above described wave regions are only qual-

itative in nature, and based on observations. Not all of the aspects have been

rigorously proved theoretically. Especially, the nonlinear wave regimes II, III &

IV have so far been studied using simplified approximate nonlinear theories. The

present study limits itself to nonlinear 2-D waves found to occur in regions II & III,

and a full-scale nonlinear study based on the direct solution of the Navier-Stokes
equations without any a priori approximations or simplifications is done. To date,

only one other detailed direct numerical study of the thin film instability has been
reported (Salamon, Armstrong & Brown 1994). Salamon eta! (1994) assume a priori

the existence of stationary waveforms and compute their shapes by solving the 2-D

Navier-Stokes equations using FEM. The waveforms obtained in this manner are

not always stable. Their relative stability can only be ascertained through the time
integration of the initial value problem, and this is what is done in this study.
 49

air-liquid interface

Regime- I
wave inception, infmitesimal waves

Regime - II (\j\
finite amplitude permanent waves

Regime- III
solitary waves

Regime- IV
irregular, random, three-dimensional solitary waves

Figure 6.1: Surface wave regimes in thin liquid film draining down vertical
walL The figure is sketched based on the description given by Chang
(1994).
 50

6.2 Problem Definition and Non-dimensional Parameters

The model problem being considered is a thin liquid film draining down a straight

wall inclined at an angle (3 to the horizontal as shown in Fig.6.2. x-axis is along

the wall and is also referred to as the stream-wise direction. y-axis is perpendicular

to the wall, positive upwards. h(x, t) is the film thickness measured as the height
of the gas-liquid interface from the wall. The fluid flow is assumed to be laminar,

Newtonian, isothermal and incompressible. Further, it is assumed that the fluid

flow is two-dimensional, i.e., no transverse or span-wise (z-direction) variations.

The dependent and independent variables can be non-dimensionalized in several

ways and we chose to use the unperturbed film thickness h 0 as the characteristic

length and v / h0 as the characteristic velocity, which is also called the viscous scaling.
The non-dimensional governing equations are as follows:

au au au ap
-+u-+v- = --+Gsmf3+-2 + -2
. au au
2 2
(6.1)
at ax ay ax ax ay
(6.2)

(6.3)

The boundary conditions at the wall are the no-slip and impermeability conditions:

u = 0} at y = 0. (6.4)
v = 0

Neglecting the motions in the gas phase, the boundary conditions at the free surface

are as follows:
 51

free surface

Figure 6.2: The physical configuration of a thin layer flowing down an

inclined plane.

Oun ouT)
( or +on 0

3S at y == h(x, t). (6.5)

Pa-R
ah
at
+ u ax
a;, v

The first two conditions impose the continuity of the stress across the gas-liquid
interface and the last condition represents the fact that the gas-liquid interface is

a material surface and gets advected with the fluid particle velocity. As already
explained in Chapters 2 & 5, only the two stress continuity conditions are imposed

in the solution of the fluid flow given by Eq.(6.1-6.3) and the kinematic condition is
used to update the free surface location every time step.

The non-dimensional numbers that appear due to our non-dimensionalization

are:

• G == gh6j /1 2 , measure of film thickness 1 .

1 Non-dimensional number G is referred to as Galileo number by some authors

 .52

• S = ah 0 /(3pv 2 ), measure of surface tension.

The non-dimensional parameter S is not completely independent of G, since it is

dependent on unperturbed film thickness h0 . We rewrite S as follows: S = TG 113 ,

where T = a /(3pv 4 13 g 113 ). The non-dimensional parameter T 2 is a function of the

fluid properties alone and fixing T in effect fixes the fluid and the film thickness is

controlled by varying G. T for pure water at 15°C is approximately 950 (Chang

1994 ).

Some authors use the horizontal velocity at the interface U0 as the reference

velocity resulting in the following non-dimensional parameters:

• Re = U0 h0 /v, Reynolds number

• We= a/(pgh~sin/3), Weber number

The Reynolds number and the Weber number are related to the non-dimensional

parameters G and S as follows:

Re = Gsin/3 3S
We= G . a
2 Sill f./

6.3 Literature Review

6.3.1 Nusselt Film Flow

If the flow is assumed to be laminar, smooth, two-dimensional and fully developed,

the mathematical system given by Eqs.(6.1-6.3), reduce to the following simplified

system (in dimensional co-ordinates):

d?u gsinf3
-+
dy2 1/
=0 (6.6)

2
1 = 3T = 0' f(pv4/3 9 tf3) is called the Kapitza number by some researchers.
 53

dp
- = -pg cos 3 (607)
dy

The boundary conditions given by Eqso(6..±o6o5) also simplify to:

u = 0 at y = 0 (608)

-
du
dy
p
-
-
=
0

0
l at !J = ho (609)

Eqso(606.607) can be solved subject to the boundary conditions given by Eqso(6o8,6o9),

resulting in:
2
U(y) gsinJ(h oY--
y )
lJ 2
v 0 (6010)

P(y) pgcosj3(h 0 -y)

The above solution is usually attributed to Nusselt (1916) even though it is a special

case of the series solution obtained by Hopf (1910) for flow over open, inclined, two-

dimensional rectangular channels (Churchill 1988; Fulford 1964)0 The velocity at

the interface Uo is given by:

gh6 0 j3 (6011)
Uo = z; sm
2

and is 1.5 times the mean velocity given by:

gh6 j3
Uav =-Sill
0

(6012)
3v

The simple Nusselt flow given by Eqo(6o10) is a solution to the 2-D Navier-

Stokes equations (Eqs0601-6o3) and satisfies the boundary conditions (Eqso6.4,6o5)o

However, this basic flat-film solution is not observed in nature, since it is an unstable
 54

solution for most ranges of parameters of interest. Hence, much of the theoretical
interest has been to study the stability of the N usselt flat film flow and to find the

range of parameters for which this basic solution is unstable.

6.3.2 Linear stability

The classical hydrodynamic stability analysis involves imposing an infinitesimally

small perturbation on the basic flow which in this case is the Nusselt film flow given

by Eq.(6.10), and study its growth in either time or space, as the case may be. An

infinitesimally small perturbation on the flat film is imposed as follows:

h ~1
1
h(x,t)=1+h 1 (x,t), (6.13)

This perturbation of the interface results in small deviations of the velocity and

pressure field from the base flow:

u(x,y,t) U(y)+u 1(x,y,t)

v(x,y,t) v1(x, y, t) (6.14)

p(x, y, t) P(y) + p (x,y,t)

1

The deviations u 1 ,v 1 and p1 are infinitesimally small, 1.e., u 1 ~ 1, v 1 ~ 1 and

p 1 ~ 1. The above perturbed flow field is substituted into the two-dimensional3

Navier-stokes equations and linearized in u 1 ,v 1 and p1 resulting in:

1 1 1

-OU + U( y ) -OU + vdU

- = -ap
I
- + "172
v u
I
(6.15)
ot ox dy ax

av~ U( ) av~ ap~ ..,2

-+
at
y -=--+ v v
ax ay
I
(6.16)

3 As per the Squire's theorem (White 1991) it is enough to study the stability of parallel flow U(y)
for a two-dimensional disturbance.
 55

( 6.17)

Pressure is eliminated from Eqs.( 6.15,6.16) by cross-differentiating and subtracting

one from the other. The continuity equation is automatically satisfied by introducing
a stream function u•' of the form:

, ow'
u =-
I ()1/
=--
()y
0" (6.18)
rJ.r

The above analysis results in a fourth-order partial differential equation for 1/J' of
the form:

(6.19)

The initial disturbance is usually assumed to be exponential as follows:

h'( x, t) 8exp[i (kx- wt)]

' 8~ 1 (6.20)
8exp[ik(x- ct)]

One of the physical motivations behind using the exponential form for the distur-

bance is that, any function can expressed as a Fourier sum. Thus,. knowing the
response of the system to each of the harmonic disturbances, the response of the

system for any general disturbance can be constructed through linear superposition.

If the film flow is stable for all spatial wave numbers k or all frequencies w, then, the
base flow given by U(y) can be said to be a stable solution. Even if there exists a

single unstable k or w, then the base solution is said to be unstable. In the harmonic

disturbance given by Eq.(6.20), k = 2;rh 0 / .\ is the non-dimensional wave number, w

is the non-dimensional angular frequency and c = w/ k is the non-dimensional phase

speed. ,\ is the wave length of the disturbance and 8 is the amplitude of the distur-
bance. If k is fixed to be strictly real, and w is allowed to be complex, the analysis
 56

is termed temporal stability analysis. The disturbance in this case is periodic in the

stream-wise direction and either decaying or growing in time depending on the sign

of imaginary part of w. The analysis is called spatial stability analysis, if w is fixed

to be strictly real and k is allowed to take complex values. The disturbance in this

case is periodic in time and either growing or decaying in the stream-wise direction
depending on the sign of the imaginary part of k. Both the analyses give the same

critical conditions. Spatial stability and temporal stability results can be converted

into one another through the Gaster transformation.

The stream function can be written as:

1/J' = 5q6(y)exp[ik(x- ct)], (6.21)

and is substituted into Eq.(6.19), resulting in the well known Orr-Sommerfeld equa-

tion.

(6.22)

The superscript '" denotes differentiation with respect to y. The four boundary

conditions needed to solve the above fourth order linear differential equation are
provided by the no-slip conditions at the wall and stress continuity equations at the

interface.

0
- } at y = 0 (6.23)
<P' = 0

<P"' + i3k 3 s
(6.24)
0
 .57

At the interface, the boundary conditions are imposed at y =l instead of y = 1 + h',

since 6 is infinitesimally small and the error introduced is 0(8 2 ). The kinematic
condition at the interface provides another relation:

c - U0 = 6 at y = l (6.25)

The temporal stability is investigated by solving the Orr-Sommerfeld equation (Eq.

6.22) subject to the four boundary conditions given by Eqs.(6.23,6.24). Knowing

o(y), the complex wave speed is calculated from the kinematic condition (Eq.6.25).

Let c = Cr + ic;. Cr is the wave speed and kc; determines the growth rate of the

disturbance. A positive value of c; implies there is exponential growth of the initial

disturbance amplitude and the system is unstable, and vice versa if e; is negative.

c; = 0 gives us the critical conditions, for which the disturbance neither grows nor
decays.

Yih(1955) and Benjamin(1957) first formulated the interfacial instability in thin

film flows in terms of the Orr-Sommerfeld equation. The Orr-Sommerfeld equation

cannot be solved in a closed form manner. However, for very long wavelength dis-

turbances ( k « 1) and small Reynolds numbers ( Re "' 1), approximate analytical

solutions for the Orr-Sommerfeld equation have been obtained by Benjamin(1957)

and Yih( 1963) using perturbation analysis. With certain simplifying assumptions,

Anshus & Goren (1964) and Krantz & Goren (1971) obtained linear stability re-

sults not restricted to small Reynolds numbers. A complete solution of the Orr-

Sommerfeld equation valid for all Reynolds numbers and all wave numbers requires

numerical methods and this was done by Whitaker (1964), Pierson & Whitaker

(1977) and Chin, Abernathy & Bertschy (1986) among many others.
 58

The important conclusions that fall out of the linear stability analysis can be

shown in the form of a phase diagram (Fig.6.3). The linear stability analysis predicts

a cut-off Reynolds number for the onset of instability as follows:

(6.26)

The experimental confirmation of the above linear stability result for Gc was pro-

vided by Liu, Paul & Gollub (1993). In the case of vertically draining thin film,

Gc = 0, implying that a vertically draining thin film is unstable for all Reynolds
numbers. Linear stability analysis provides two neutral curves, kc = 0 and kc =
kc(G,S,/3) and for G > Gc the region of instability is 0 < k < kc(G,S,/3). Short
wavelength disturbances are damped by the capillarity action. For infinitesimal

wavenumbers (k -+ 0) linear stability predicts the phase speed to be c = G. If

a non-dimensionalization based on the velocity at the interface (Uo) is used, then

wave speed for k -+ 0 is c = 2U0 . Another quantity of interest which comes out of
the linear stability analysis is the wavenumber for which the disturbance grows the

fastest, denoted by kM. If no artificial excitations are imposed, then the waves that

are most likely to be detected would be of the maximum growth rate wavenumber

kM. This was experimentally confirmed by Alekseenko, Nakoryakov & Pokusaev

(1985).

6.3.3 Nonlinear Stability Analysis Based On Lubrication Approximation

The linear stability analysis described in the previous subsection is expected to

be valid only in the wave inception region. Beyond wave inception, the wave am-

plitude grows exponentially and the nonlinearities can no longer be neglected. In

the case of very thin layers (G ~ 0(1) ), the cut-off wavenumber kc is small, so
 59

linearly stable

___... ··· ....-···

--······
.......lioe~rtlv
· "'
unstable
_ ....--······

Figure 6.3: Stability diagram based on linear stability analysis

that the nonlinear extension of the stability analysis can be accommodated by the

small-wavenumber approximation and the lubrication theory. The aim has been to
replace the two-dimensional Navier-Stokes equations with a single nonlinear partial

differential equation for the film thickness h(x, t), and examine the stability of the
thin film in terms of this nonlinear evolution equation. Benney (1966) first derived

the nonlinear evolution equation for the two-dimensional flows, which was extended
to three dimensions by Roskes (1970). There have been a number of extensions

of this as discussed by Atherton & Homsy (1976) and Lin & Wang (1985). The
methodology used for deriving the nonlinear evolution equation is briefly described

next.
Let lc be the characteristic length in the x-direction and is typically proportional
to the wavelength of the disturbance. For very thin layers •he disturbances are of
 60

very long wavelength, i.e. lc :::> h 0 . Define a perturbation expansiOn parameter

E = ho/lc, E « 1. The original system is rescaled as follows:

~ = EX ; TJ = y ; T = Ef (6. 27)

The following scaling for the non-dimensional parameters is assumed:

G ~ 0(1) ; S = St- 2 ; S ~ 0(1) (6.28)

Surface tension is a stabilizing term in thin film instability as it weakens the initial

exponential growth resulting in finite amplitude permanent waves. S ~ 0( c 2)

is easily satisfied by most liquids. Only highly viscous oils may have a relatively

smaller 5. The dependent variables u, v and pare expressed as a power series in t

as follows:

(6.29)

(6.30)

P = Po + tPt + t 2 P2 + · · · ( 6.31)

The above asymptotic expansions are substituted into the governing equations (Eqs.
6.1 - 6.3) and boundary conditions (Eqs.6.4,6.5) and solved by grouping together

equal order terms in <.. Going through the algebra one can derive the following

expressions:

(6.32)

(6.33)
 61

po = G cos ;3 ( h - ry) - :{5 he( (6.34)

( Gh~ cos J- :fShw.) (~ 17 2 - hry)

(6.35)
+~G sin Jh, Gr7 3
- h
2
ry) + ~G 2 sin 2 ;3hh( GT/ 4
- h 3
ry)

(6.36)

The leading order terms uo, Vo and Po are the lubrication theory results. The

evaluated velocity is substituted into the kinematic equation, and all terms of higher

order than 0( t:) are dropped, resulting in the following nonlinear evolution equation

for the film thickness h(x, t):

In obtaining the above equation the system has been rescaled back to the original

(x, y, t). The second term in the above equation describes the wave propagation and
steepening. The third term is a destabilizing force and is due to the mean shear flow.

The fourth and the fifth terms are stabilizing forces and are due to hydrostatic pres-
sure and surface tension force, respectively. For two-dimensional films, Siva.shinsky

& Michelson ( 1980) showed that in the large surface tension limit ( S -+ oo), the

evolution equation can be further simplified to the Kuramoto-Sivashinsky equation:

(6.38)
 62

where ¢ = h - 1 is the wave amplitude and a, b, and c are appropriate constants

and angle {3 = ~.

Much of the nonlinear stability analysis of thin film flows has been done using the

Kuramoto-Sivashinsky equation (Eq.6.38) and the longwave evolution equation of

the Benney type (Eq.6.37) due to their relative simplicity. The longwave evolution
equation is a fully nonlinear equation and cannot be solved analytically. However, a

weakly nonlinear analysis can be performed by retaining only the fundamental mode

and its first harmonic. The weakly nonlinear analysis of Lin (1969), Gjevik (1970)

and Nakaya (1975), all predict supercritical finite amplitude waves for wavenumbt;'rs

just below the linear cut-off wavenumber kc. They also predict a new phase bound-

ary ks( G, S, {3), below which no permanent waves can be observed (Fig.6.4). Pumir,

Manneville & Pomeau (1983) and Joo, Davis & Bankoff (1991b ), solved the long-

wave evolution equation numerically. They confirm the existence of supercritical

saturated permanent waves for wavenumbers slightly lower than kc. For wavenum-

bers much below the linear cut-off wavenumber kc, they observe very long waves of

the 'solitary' type.

6.3.4 Nonlinear Stability Analysis based on Boundary Layer equations

The nonlinear stability analysis described in the previous subsection is valid only for

very thin layers, i.e. G'"" 0(1). For relatively thick layers, the inertial terms can no

longer be neglected and the lubrication theory breaks down. Another approximate

nonlinear stability theory valid for relatively thicker layers, i.e. G '"" 0(100), has
been developed based on the boundary layer equations. With only the assumption

that the disturbance wavelength is much longer than the film thickness (,\. ::;};> ho),
 63

linearly stable

subcritical

Figure 6.4: Stability diagram based on weakly nonlinear stability analysis

the original Navier-Stokes equations given by Eqs.(6.1-6.3) can be simplified to the

following boundary layer equations:

au au au . a2u
at + uax + v ay = G Sill f3- G cos f3h, + 3Shxxx + ay2 (6.39)

au av
-+-=0 (6.40)
ax ay
The differences from the original Navier-Stokes equations is that pressure is taken

to be the sum of hydrostatic and capillary pressures only, and the diffusion terms in

the stream-wise direction are dropped. The boundary conditions also simplify to:

u
(6.41)
v
 64

ou
oy (6.42)
oh oh
-+u-
ot ax
The above listed Boundary Layer system is much easier to solve than the original
N avier-Stokes equations.

The boundary layer equations can further be simplified by assuming a przorz

the velocity profile and this approach is called the Integral Boundary Layer theory

(Chang 1994) and is in the same spirit as the original Von Karman - Pohlhausen

momentum integral approach. A number of different evolution equations have been

derived, depending on the velocity profile imposed. A uniform velocity profile results

in the shallow water theory of Dressler(1949). Alekseenko eta! (1985) imposed-a

parabolic profile to derive a more appropriate evolution equation, which is further

extended recently by Prokopiou, Cheng & Chang (1991) to higher-order accuracy.

In the integral boundary layer theory, two evolution equations are obtained, one for

the flow rate q =.J0hudy, and the other for the film thickness h(x,t). A parabolic

velocity profile results in the following evolution equations:

oq a h)
ot + 1.2 ox (q 2
· ( •
= h G sm j3- G cos f3hx + 3S ( R1 ) X
) q
-3 h2 (6.43)

oh oq _ (6.44)
at+ ax- 0

With some additional approximations, the above two evolution equations can be
combined into one second-order wave equation for the film thickness h( x, t) ( Alekseenko

etal1985). Lee & Mei (1995) derived an evolution equation, valid for large Reynolds
numbers and moderate Weber numbers, based on the boundary layer approxima-

tion and a priori specification of a parabolic velocity profile. Chang, Demekhin

 65

& Kopelevich (1')93) solved the longwave boundary layer equations numerically,
without a priorz specification of the velocity profile.

In majority of the nonlinear studies based on the longwave boundary layer equa-

tions, finite amplitude permanent waves are assumed a priori, and the stationary

equations are solved in a frame of reference translating with the wave speed c. Chang
eta] ( 1993) based on this type of analysis. predict slow-moving short nearly sinu-

soidal waves referred to as /I family, and fast-moving long solitary waves with one
or more primary humps referred to as the 1 2 family. Based on a detailed bifurca-

tion analysis of the third-order dynamical system resulting from the assumption of

stationary waves, Lee eta! ( 1995) find a variety of bifurcation phenomena, such as

limit cycles, heteroclininc orbits, chaotic attractors and homoclinic orbits.

6.3.5 Direct Numerical Studies

To study the complicated nonlinear flow developments without any a priori as-

sumptions such as those made in the long-wave theory and the boundary layer

theory, the complete Na vier-Stokes equations need to be solved. Due to the irreg-
ular and time-varying flow domains involved, Finite Element Method (FEM) has

been the popular method of choice (Bach & Villadsen 1984; Kheshgi & Scriven

1987; Malamataris & Papanastasiou 1991; Salamon, Armstrong & Brown 1994).

Ho & Patera (1990) used Legendre Spectral Element Method, which is a higher

order FEM. Bach & Villadsen(1984), Kheshgi & Scriven (1987) and Malamataris

& Papanastasiou (1991) use a Lagrangian Finite Element Method to handle the
moving boundary and control the mesh distortion through rezoning. Ho & Patera

(1990) use a mixed Lagrangian- Eulerian approach and Salamon etal (1994) use the

concept of vertical spines.

 66

Kheshgi & Scriven(1987) and Ho & Patera(1990) make comparisons with the
Orr-Sommerfeld linear stability predictions for the neutral wavenumbers. Ho &

Patera (1990) compared the wave profile and wave speeds with the experimental

measurements of Kapitza & Kapitza (1949). In finite length domains without pe-

riodic boundary conditions, Kheshgi & Scriven (1987) could not observe steady

traveling waves. However, Malamataris & Papanastasiou (1991) using a modified

free boundary condition show the existence of traveling waves in truncated domains.

Salamon etal (1994) did a comprehensive direct numerical study of traveling

waves in vertical thin films and compared their results against the approximate

long-wave and boundary layer theories. They assume a priori the existence of steady

traveling waves and rewrite the governing equations in a frame of reference translat-
ing with the wave speed. They thus solve the steady state N a vier-Stokes equations

and compute the flow field, free surface profile and wave speed simultaneously for

a given wavelength and Reynolds number. They found good agreement with the

long-wave theory for small amplitude waves, but found their results to qualitatively

diverge from the longwave results for large amplitude waves. The also studied the

nonlinear interactions between the waves and the secondary subharmonic bifurca-

tions to longer waves.

6.3.6 Experimental Studies

Experimental studies of falling films have been done by, among others, Kapitza

& Kapitza (1949), Krantz & Goren (1971), Portalski & Clegg (1972), Alekseenko
etal (1985), and Lacy, Sheintuck & Dukler (1991). As summarized by Alekseenko

etal(1985), in most of the experiments performed, two-dimensional regular waves

are observed only near the wave-inception line. The waves soon become three-
 67

dimensional and irregular. In order to obtain two-dimensional wavetrains. the flow

was disturbed at a fixed frequency by. for example. wire vibrations (Krantz & Goren

1971) and pulsations of flow rate il\:apitza & I\:apitza 1949; Alekseenko eta! 1985).

Lacy eta!. ( 1991) performed the experiments without artificial perturbations, and

noticed that draining films exhibit deterministic chaos.

More recently. Liu, Paul and Gollub ( 1993) using sophisticated experimental

techniques. measured the cut-off Reynolds number as a function of the angle of

inclination of the plate. They found good agreement with the linear stability pre-

dictions for the cut-off Reynolds number. growth rates and wave velocities. They

also demonstrated the primary instability to be convective in character and hence

extremely sensitive to external noise at the source. Liu & Gollub ( 1993) found the

primary surface wave instability to be susceptible to both side-band and subhar-

monic secondary instabilities but in different ranges of frequency. These secondary

instabilities are also convective in nature and are one of the routes to spatio-temporal

chaos. Liu & Gollub ( 1994) performed a systematic analysis of solitary wave dynam-

ics in two dimensions and found the velocity of the solitary wave to be proportional

to the wave height. Due to this amplitude dependence of the wave speed, the bigger

amplitude waves travel faster and coalesce with smaller waves.

6.4 Course of the Study

From the review of the literature related to the thin film flow as outlined briefly in

the previous section, we can infer the following:

l. The solution of the linear stability problem based on the Orr-Sommerfeld

equation can be considered to be complete. The predictions of the linear

stability theory have also been validated experimentally.

 68

2. The thrust of the theoretical research in the past few years has been to derive

model equations capable of studying the nonlinear evolution of the thin film

flow. The longwave theories based on the lubrication approximation and the

boundary layer approximation have been partially successful in explaining the

nonlinear flow regimes. However, these nonlinear theories have only limited
ranges of validity.

3. A complete nonlinear study of the thin film instability requires the solution of

the Navier-Stokes equations. However, due to the presence of the free bound-

ary, numerical solution of the Navier-Stokes equations is extremely difficult

and computationally intensive. With the exception of Salamon eta] (1994), all

the other numerical studies are limited in scope and primarily aim to establish

that their numerical procedure is capable of simulating the thin film flows.

The aim of this research is:

1. To develop an efficient and accurate numerical procedure capable of solving

the free boundary Navier-Stokes equations arising in the thin film flows. In

particular, the numerical method should be based on the solution of the initial

value problem, so that the most stable waveform for the simulation parame-

ters can be directly computed. The essentials of the numerical procedure are

described in Chapters 3 & 5.

2. Utilize the numerical procedure to perform a detailed temporal stability anal-

ysis. In particular, we wish to study in detail the nonlinear flow regimes,

namely the finite amplitude permanent wave regimes (Region II in Fig.6.1)

and the solitary wave regimes (Region III in Fig.6.1). This is discussed in

sections 6.5 & 6.6.

 69

:l. Study the wave interaction processes. such as. wave splitting and wave mergers,

known to occur in thin film flow. This is discussed in section 6.7.1.

l. Lastly, we wish to establish our numerical procedure as a complement to !abo-

ratory experiments. What is studied in the laboratory, is the spatial stability

analysis. A periodic disturbance is imposed at the inlet and its evolution in

the streamwise direction is measured. To this end. we consider nonperiodic

domains with absorbing boundary conditions at the exit and make compar-

isons with the experiments of Liu & Gollub ( 1994 ). This is discussed in section
6.7.3.

6.5 Comparison with experiments and previous numerical

simulations

Kapitza & Kapitza (1949) obtained two-dimensional permanent waves on thin films

draining down a vertical cylinder by artificially perturbing the flow rate at a fixed

frequency and amplitude. Ho & Patera (1990) and Salamon eta! (1994) compared

their numerical results with the experimental findings of Kapitza & Kapitza (1949)

for the two conditions listed in Table 6.1. These two experimental conditions are

numerically simulated through our full-scale numerical procedure and compared with

the experimental results of Kapitza & Kapitza ( 1949) and the numerical results of

Ho & Patera (1990) and Salamon eta! (1994).

fluid !L
o
(em")
s 2 .\(em) Q("";") G s k
alcohol at 16.8°C 29 1.77 0.123 18.2 463.7 0.07
water at 19.6°C 74 0.80 0.201 60.0 4410 0.14

Table 6.1: Experimental conditions of Kapitza & Kapitza(1949) and the

relevant non-dimensional parameters.
 70

What is investigated in experiments is the spatial stability, where a periodic

disturbance is imposed at the inlet and its evolution in streamwise direction is de-

termined. For periodic forcing at the inlet, sufficiently far downstream from the
inlet, Kapitza & Kapitza (1949) observed saturation of the disturbance and the

waves travel downstream with fixed wave speed and wavelength. Numerically, how-
ever, the equilibrium wave profile and wave speed are obtained through temporal

stability analysis. A streamwise periodic disturbance with wavelength A observed

in the experiments is imposed as follows:

h(O,t) = 1 +bcos(kx), ( 6.45)

and its evolution in time is computed through the direct numerical solution of the full

nonlinear system given by Eqs.6.1-6.3. k = 27rh 0 / A is the nondimensional wavenum-

ber, and the boundary conditions in the streamwise direction are:

<f>(x = 0) = <f>(x =A), (6.46)

where ¢> stands for h, u, v and p. Saturation of the disturbance in time implies

steady traveling wave with fixed wave speed and wave profile. These numerically

computed wave profile and wave speed are then compared with those obtained by

Kapitza & Kapitza (1949).

The mean film thickness h0 is needed for calculating the nondimensional param-

eters G, S and k. However, Kapitza & Kapitza (1949), do not provide the mean

film thickness, and instead give the flow rate. Using the experimental flow rate q,

the mean film thickness h0 is computed as:

3gvq) 1/3
ho = ( (6.4 7)
 71

This mean film thickness is then used to compute the nondimensional parameters

G. Sand k. In the numerical formulation, the mean film thickness remains fixed in

time. but the flow rate could vary with time. In fact, with the onset of wave motion,

the flow rate increases. and the final steadv flow rate in our numerical simulation

would be greater than the experimental flow rate. The results of the nonlinear flow
evolution are presented in terms of the film thickness h(x, t) at various instants of

time and the spatial spectral coefficients en ( t), defined as

n=N
h(x,t) = L Cn(t)emkx. (6.48)
n==-N

where 2N is the number of mesh divisions in the streamwise direction. The initial
amplitude of the disturbance is set to 8 = 0.05.

The nonlinear flow evolution for the experimental condition G = 18.2, S = 463. 7,

f3 = 1r /2 and k = 0.07 is shown in Fig.6.5. The experimental wavenumber k = 0.07

is smaller than the linear cut-off wavenumber kc = 0.29 (Whitaker 1964). Thus,

in agreement with the linear theory, the amplitude of the disturbance grows expo-

nentially initially (Fig.6.5a). Eventually, however, due to the stabilizing capillary

force, the growth is arrested, resulting in steady traveling waves (Fig.6.5b ). The

dynamics of the the nonlinear processes can be better quantified through the spatial

spectral coefficients en(t) (Fig.6.5c). The fundamental mode and its harmonics grow

exponentially initially. Eventually, due to the nonlinear stabilizing, they saturate in

time, resulting in finite amplitude permanent waves traveling downstream at fixed

wave speed c. The equilibrium wave profile obtained using the present method is
compared with that reported by Kapitza & Kapitza (1949), Ho & Patera (1990)

and Salamon eta! (1994) in Fig.6.6. It is in excellent agreement with the experimen-
tal profile of Kapitza & Kapitza (1949) and the numerical profiles of Ho & Patera
 72

(1990) and Salamon eta! (1994). The dominant crest has a tear drop profile with

a sharp downstream slope and a long gently sloping tail. The sharp hump is pre-

ceded immediately downstream by small capillary waves. The free surface profile

is 'solitary wave' type and corresponds to the 'single wave' observed by Kapitza

& Kapitza (1949) for wave numbers much smaller than the critical wave number
kc. The computed wave speed is 23.1cm/s which compares very well with the wave

speed of 23.5cm/s reported by Salamon eta! (1994) and 24.7cm/s reported by Ho &
Patera (1990). The experimental wave speed reported by Kapitza & Kapitza (1949)

is 19.5cm/ s which is appreciably lesser than the numerically computed wave speeds.

As pointed out by Ho & Patera (1990) and Salamon eta! (1994), the experimental

flow rate is imposed only as an initial condition in the numerical simulation, but

due to the onset of wave motion, there is an increase in the flow rate. Salamon eta!

(1994) obtained better agreement with the experimental results when they adjusted

their mean film thickness to take into account this aspect. Also shown in Fig.6.6

is the permanent wave profile obtained through the numerical solution of the inte-
gral boundary layer system based on the assumption of a parabolic velocity profile

(Eqs.6.43,6.44). We find the wave profile obtained through the Integral Boundary

Layer system to be in good agreement with that obtained experimentally and numer-

ically through full-scale computations. The results shown in Figs.6.5 are obtained

using a 65 x 11 mesh with time step control parameter C F L = 0.1. Doubling the

number of nodes in x and y directions and reducing the time step size by half re-

sulted in identical results as shown in Fig.6.7. The evolution of the harmonic modes

is exactly the same for various mesh and time step parameters. Thus the solution

computed by us is independent of the mesh size and the time step.

 73

(a) (b)
1.1 .----~-~----~

1.5

0.9L--~-~-~-~- 0.5'-----~-....__-~_....__j

o 20 40 60 80 0 20 40 60 80
X X

(c)
0.25.---,---,..-----,----,-----,-----,------,

n=+/-1

n=+l- 2

n=+l- 3

10 20 30 40 50 60
t

Figure 6.5: First experimental condition of Kapitza & Kapitza (1949),

namely, G=18.2, 5=463.7, k=0.07, j3 = ~: (a) h(x, t) fort= 0.0 tot= 2.0 in
steps of 0.1; (b) h(x,t) fort= 64.6 tot= 64.8 in steps of 0.1; (c) Evolution
of spectral coefficients (lcn(t)l versus time).
 74

(a) (b) (c) (d) (e)

Figure 6.6: Wave shapes for G = 18.2, S = 463.7, k = 0.07, f3 = f, obtained

by: (a) Kapitza & Kapitza (1949), (b) Salamon etal (1994), (c) Ho &
Patera (1990), (d) Present numerical method, and (e) integral boundary
layer theory.
 75

(a)
0.25 , - - - - , - - - - r - - - , . - - - - - - - , - - -

0.2 n =+1- 1

- 65x11, CFL=0.1
0.15
n =+1- 2 - - 65x11, CFL=O.OS
!cn(t)!
0.1 ·- ·- 65x21, CFL=0.1
· · 129x11 , CFL=O .05

10 20 30 40 50

(b)
2,------~----~----~----~~

1.5 - 65x11, CFL=0.1

..c: - - 65x11, CFL=0.05
· - - 65x21, CFL=0.1
· · · · · 129x11, CFL=O.OS

0.5 l - -_ _ _ _...J..--_ _ _ __ . __ _ _ __ . __ _ _ _
~---'

0 20 40 60 80
X

Figure 6.7: Mesh and time step independent study for G = 18.2, S = 463.7,
k = 0.07 and !3 = ~· (a) Evolution of spatial spectral coefficients (icn(t)l
versus time); (b) Free surface profile h(x,t) at t =50.
 76

The evolution of the initially imposed sinusoidal perturbation with amplitude

0.05 for the case G = 60.0, S = 4410.0, (3 = rr/2 and k = 0.14 is shown in Fig.6.8.

The wavenumber k = 0.14 is smaller than the linear cut-off wavenumber kc ~ 0.30.
From the evolution of the spatial spectral coefficients shown in Fig.6.8( a), it is

apparent that the amplitude of the wave grows exponentially initially. Unlike in the
previous case however, the harmonic modes do not saturate in time. The spectral

coefficients continuously oscillate in time, generating a quasi-periodic waveform. The

fundamental mode and its harmonics are continuously exchanging energy , without

settling to a stationary value. In fact, the total energy of the system defined as:
n=N
E(t) = "L icn(tW, (6.49)
n=-N

was found to be continuously oscillating m time. The wave profiles at the two

extremes, namely, when the fundamental mode has the minimum energy (t = 77.75)
and wheh the fundamental mode has the maximum energy (t = 78.90) are shown

in Fig.6.8(b ). These two wave profiles along with the steady wave profiles obtained

by other researchers are shown in Fig.6.9. The experimental wave profiles and the

previously reported numerical wave profiles agree better with the wave profile at t =
78.90. A numerical solution of the Integral Boundary Layer system (Eqs.6.43,6.44)
also revealed a similar quasi-periodic wave behaviour. The wave profiles obtained

using the Integral Boundary Layer theory are also shown in Fig.6.9, and are found to

be in close agreement with that obtained by us through the direct numerical solution

of the Navier-Stokes equations. The results shown in Fig.6.8 has been obtained

using a 65 x 11 mesh and time step control parameter C F L = 0.1. Numerical

experimentation with various mesh sizes and time steps revealed identical behaviour.

Starting with a different initial condition, namely an initial disturbance amplitude

 77

of b = 0.50 instead of b = 0.0.5 also resulted in a similar behaviour (Fig.6.10).

Thus. the behaviour exhibited by the system can be considered to be independent

of numerical mesh and initial condition.

To summarize. for the first experimental condition (Fig.6.5). we obtain a steady

travelling waveform with solitary hump shape, and both the wave profile and wave

speed are in good agreement with that reported previously. For the second exper-

imental condition (Fig.6.8). however. our full-scale computations predict a quasi-

periodic waveform in contrast to the stationary waveforms reported previously.

Recent experiments of Liu & Gollub ( 1994) indicates that this type of behaviour

is also exhibited. In certain ranges of frequencies they could not observe saturated

waves. This quasi-periodic behaviour is investigated in detail in the next section.

Assuming for the moment, that our numerical results are correct, it remains to be

explained why Kapitza & Kapitza (1949), Ho & Patera (1990) and Salamon eta!

(1994) did not observe this type of a behaviour. Since, the wave profiles are not

necessarily sinusoidal, Kapitza & Kapitza (1949), defined the wave amplitude a as:

(6.50)

where hmax and hmin, are respectively, the maximum film thickness and the minimum

film thickness. The wave amplitudes for different flow rates are given in Fig.13

of Kapitza & Kapitza (1949), The second experimental condition simulated by

us corresponds to Q = 0.201cm 2 / s, and for this flow rate, they give two wave

amplitudes, possibly because they did not observe well defined stationary travelling

waves. The other possibility is that their test section is not long enough for them

to be able to observe well defined quasi-periodic behaviour. From Fig.6.8(a), we

find that the system needs about 40 nondimensional time units for the onset of
 78

(a)

10 20 30 40 50 60 70 80
t
(b)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
xj>.

Figure 6.8: Nonlinear wave evolution for the second experimental condi-
tion of Kapitza & Kapitza (1949), namely, G = 60.0, S = 4410.0, k = 0.14
and (3 = 1rj2. (a) Evolution of the spatial spectral coefficients icn(t)l; and
(b) film thickness at t = 77.75 (broken line) and t = 78.90 (continuous line).
 79

~~

(a) (b) (c) (d) (e) (f) (g)

Figure 6.9: Wave shapes for G = 60.0, S = 4410, k = 0.14 and j3 = ~'
obtained by: (a) Kapitza & Kapitza (1949), (b) Salamon etal {1994), (c)
Ho & Patera (1990), (d) Present numerical simulation at t = 77.75, (e)
Present numerical simulation at t = 78.90, (f,g) integral boundary layer
theory.
 80

(a)

90

(b)

5 10 15 20 25 30 35 40 45 50
t

Figure 6.10: Nonlinear evolution of the spectral coefficients lcn(t)l with

different initial conditions for G = 60.0, S = 4410, k = 0.14 and /3 = i· (a)
Initial disturbance is a sinusoidal perturbation with amplitude b = 0.05;
(b) initial disturbance is a sinusoidal perturbation with amplitude b =
0.50.
 81

quasi-periodic behaviour. Based on our nondimensionalization this works out to

be t = l.34s, and at a wave speed of l9.7cm/ s, it would be approximately 26cm,

before well defined quasi-periodic behaviour can be observed. The length of the test

section used by Kapitza & Kapitza ( 1949) is only 17cm. Due to a priori assumption

of steady traveling waveforms. Salamon eta! (1994) could have missed this quasi-
periodic behaviour.

6.6 Temporal Stability Analysis

From the previous section, it appears that for wavenumbers smaller than the linear

cutoff wavenumber kc, we do not always obtain steady travelling permanent wave-

forms. A similar result was shown experimentally by Liu & Gollub (1994). For

excitation frequencies w closer to linear cut-off frequency We, they observed nearly

sinusoidal saturated waves. As the frequency is reduced, however, in certain ranges

of frequencies, they observed quasi-periodic evolution in the downstream direction.

The film thicknes!l is however, still periodic in time, but does not saturate in the

streamwise direction. For much smaller wavenumbers, they observed steady trav-

elling solitary waves. Thus, experimental evidence exists for quasi-periodic trav-

elling waveforms. The experimentally investigated thin-film stability is a spatial

one, where a periodic disturbances is imposed at the inlet and its evolution in the

downstream direction is measured. Our full-scale computations, from the previous

section, tell us that this type of quasi-periodic travelling waves are also observed in

the temporal stability analysis. The quasi-periodic behaviour was also observed by

Hooper & Grimshaw (1985) in their solutions of the Kuramoto-Sivashinsky (KS)

equation and by Joo & Davis (1992) in their numerical solutions of the longwave

evolution equation. Hooper & Grimshaw (1985) observed that for some wavenum-
 82

hers the fundamental mode and its first harmonic are in a 'bouncy state', with
continuous exchange of energy between them. Our own numerical solution of the

integral boundary layer equations with the assumption of parabolic velocity profile

also indicate that this type of quasi-periodic behaviour is exhibited. A detailed

temporal stability analysis of the thin film flow is thus undertaken in this section,
to obtain the phase boundaries for this quasi-periodic waveforms.

The long time behaviour of the system is uniquely determined by G, S, (3 and

k. The nondimensional number S can be written as S = TG 113 , where the new

parameter T = a/(3pv 4 13 g 113 ) depends only on the fluid properties. We set (3 = 1rj2

and T = 100, and focus on vertical fluids with moderate surface tension. Temporal

stability of the vertical thin film is investigated in the range 5 :::; G :::; 100. For

each G, initial disturbances of various wavenumbers k and amplitude 8 are imposed

and their nonlinear evolution is computed through the direct numerical solution

of the full-scale system given by Eqs.6.1-6.3, subject to the boundary conditions

given by Eqs.6.4,6.5. In addition, periodic boundary conditions are imposed in the

streamwise direction and the initial flow field is imposed based on the lubrication

approximation as follows:
2
u(x,y,O) = Gsin(J [h(x,O)y- Y ] (6.51)
2

. (8hBx lt=o) 2·
v(x,y,O) = -Gsm(3 y
2
(6.52)

The initial amplitude of the disturbance 8 is set to 0.05 for all the cases simulated
in this section. First, results for G = 5, 25&100 are presented. For each of these

Reynolds number, three wavenumbers, namely, k = kM, kM/2&kM/4 are simulated,

and the results are shown next. Where required, the linear cutoff wavenumber kc and
 83

the maximum growth rate wavenumber k~t are taken to be that given by 'vVhitaker

( 1964), who solved the linear stability Orr-Sommerfeld equation numerically.

Fig.6.ll shows the evolution of the spatial spectral coefficients icn(t)i and the

equilibrium wave form for G = .S.O, T = 100.0 and three different wave numbers,

namely, k = 0.10(~ kM ), k = 0.05(~ kAJ/2) and k = 0.025(~ kM/4). The harmonic

modes increase exponentially, in all the three cases. Due to the nonlinear interac-

tions, there is transfer of energy from the fundamental mode to its harmonics. After

some very complex interactions. the harmonics modes eventually saturate in all the

three cases, implying steady travelling waveform. For k = kM, the wave profile is

not completely sinusoidal since the lower harmonics also possess non-zero energy

(Fig.6.11 b). As the wavenumber is further decreased, i.e. for k = kM /2, several-of

the harmonics are excited, and the wave profile is 'solitary wave' like, with a tear

drop hump preceded downstream by small capillary ripples (Fig.6.11d). Further

reduction in the wavenumber results in a solitary wave type profile with two pri-

mary humps (Fig.6.11f). Also shown in Fig.6.11 are the final permanent waveforms

predicted by the longwave evolution equation of the Benney type (Eq.6.37) 4 • We

find excellent agreement between the full-scale computations and the longwave pre-

dictions. Thus, for small Reynolds number or alternatively for small thickness film

flows, the longwave evolution equation gives us reasonably good predictions.

Many of the qualitative features observed in the simulations just discussed, sup-

port what is already known about the thin film dynamics. For wavenumbers lower

than kc, we find the flat Nusselt film flow to be unstable, and a small perturbation

imposed on this solution grows exponentially as predicted by the linear stability

4 The longwave predictions have been computed using the numerical code provided to us by Dr.

S.W. Joo from Wayne State University. The details of this numerical scheme are discussed in Joo,
Davis and Bankoff (1991b).
 84

(a) (b)
0.06

lcn(t)l
0.04

0.02
v
0.00~
0 100 200 300 400
0.8
0.0 0.5 1.0
t xj>.
(c) (d)
0.04
1.2

0.8
500 1000 0.0 0.5 1.0
t xj>.
(e) m
0.04
1.2 I'

..c 1.0

0.8
500 1000 1500 2000 0.0 0.5 1.0
t xj>.

Figure 6.11: Evolution of the spectral coefficients lcnl(t) and the final
waveform for G = 5, T = 100, j3 = 1rj2 and (a,b) k = 0.10(:::::: kM)i (c,d)
k = 0.05(:::::: kM/2); (e,f) k = 0.025(:::::: kM/4). The waveform shown in broken
line is obtained using the longwave evolution equation of the Benney
type.
 85

theory. However. as the wave amplitude grows. nonlinear interactions takeover, and

the energy is channeled from the fundamental mode to its harmonics. In some way,

this provided a stabilization and arrests the exponential growth. The nonlinearity

by itselL however, is not capable of sustaining well defined permanent waveforms.

The formation of well defined steady traveling waveforms on the gas-liquid inter-

face is due to the presence of surface tension, which acts as a restoring force and

cuts off the lower harmonics and limits the energy to only the first few modes. For

wavenumbers close to kc, these permanent waves are nearly sinusoidal in shape. But,

for wavenumbers much smaller than kc, these permanent waves are of the 'solitary
wave' type.

The nonlinear evolution of the spatial harmonic coefficients en( t) and the final

permanent wave form for G = 25 and k = 0.30(~ kM), k = 0.15(~ kM/2) and k =
0.075( ~ kM /4) are shown in Fig.6.12. The initial exponential growth is followed by

equilibration of the modes for all the three wavenumbers resulting in finite amplitude

permanent waveforms travelling downstream with fixed wave speed on the gas-liquid

interface. The notable difference from the G =5 case, is the substantially larger

wave amplitude with the maximum wave height more than double the mean film

thickness for the wave number k = 0.075 (Fig.6.12f). As the wavenumber is reduced

the waveform changes from nearly sinusoidal waves to broad-banded solitary waves.

Attempts to solve the longwave evolution equation of the Benney type (Eq.6.37)

for G = 25 have resulted in numerical breakdown for all the three wavenumbers.

This is to be expected since it is valid only for G"' 0(1). The Integral Boundary

Layer system given by Eqs.6.43,6.44 are numerically solved, and the resulting finite

amplitude waveforms are also shown in Fig.6.12. The agreement between the full-

scale computations and the boundary layer predictions is seen to be excellent.

 86

(a) (b)

icn(t)i
0201
0.10
.l:
1.5

1.0

0.00 ~ 0.5
0 5 10 15 20 0.0 0.5 1.0
t xj>.
(c) (d)
2.0
1.5

0.5
10 20 30 0.0 0.5 1.0
t xj>.
(e) (D
2.0
0.20
icn(t)i .l: 1.5

0.10 1.0
0.5
20 40 0.0 0.5 1.0
xj>.

Figure 6.12: Evolution of the spectral coefficients icni(t) and the final
waveform for G = 25, T = 100, (3 = 1rj2 and (a,b) k = 0.30(~ kM)i (c,d)
k = 0.15(~ kM/2); (e,f) k = 0.075(~ kM/4). The waveform shown in broken
line has been obtained using the integral boundary layer theory.
 87

The nonlinear 1vave e1·ol ution for G = 100 and k = 0.-!0( ~ kM ), k = 0.20( ~

kv/2) and k = 0.10(~ 1..·\f/-l). is shown in Fig.6.1:3. The harmonic modes evolution

shown in Fig.6.13(a) is quite different from that seen previously for the low Reynolds

numbers G = 5 and G = 25. The initial exponential growth is not followed by sat-

uration of the modes but rather results in a fluctuating behavior. This corresponds

to quasi-steady wave profiles and we do not observe steady permanent wave form

as observed in the previous case G = 5 and G = 25. The wave form shown in

Fig.6.1:3(b) is consequently not the permanent wave profile but is a wave profile as

some given instant of time. A similar but more complex time behavior is observed

in the k = kM/2 case. This probably corresponds to a quasi-periodic wave profile

with even more number of independent frequencies than the previous case k = kM.
Surprisingly for k = kM /4 we obtain saturated solitary wave form with several small
capillary waves in the front. A rigorous mesh and time step independent studies

have been done to confirm this behavior.

To understand what is happening physically, during these quasi-periodic evolu-

tions, free surface profiles at various instants of time are shown for G = 100 and

k = 0.20 in Fig.6.14. For clarity the free surface profiles are shown over a three

wavelength domain. To particularly note is the growth of the subsidiary peak just

upstream of the primary wave and its coalescence with the primary wave. At the end

of this coalescence there are only three subsidiary peaks between any two primary

humps. Very soon a new subsidiary peak develops and starts growing and merges

with the primary maxima and the process is repeated. It is this phenomenon which

results in the fluctuations of the harmonic modes and a quasi-steady wave behavior.

Recently, Liu and Gollub (1994) demonstrated a similar wave behavior experimen-

tally. Thus, in the transition from nearly sinusoidal permanent waves to solitary
 88

(a) (b)
0.4.----------------, 1.5.-----~------,

0.5'-----~----
5 10 15 0.0 0.5 1.0
t xj>.
(c) (d)
0.4.-------------.-,

1.5
.c
1.0
0.5
'----~--~------~
0.0 0.5 1.0
t xj>.
(e) (f)
0.4

v
1/p
0.0 rA
0 5 10 15 20 0.0 0.5 1.0
t xj>.

Figure 6.13: Evolution of the spectral coefficients !en(t)! and the waveform
for G = 100, T = 100.0, (3 = 1rj2 and (a,b) k = 0.40(~ kM)i (c,d) k = 0.20(~
kM/2); (e,f) k = 0.10(~ kM/4).
 89

wave train there exists a band of wave numbers (frequencies) where the solitary

humps are too closely packed resulting in very strong interaction between the pri-

mary hump and the subsidiary peaks, and no permanent wave form is observed in
this case, but rather a quasi-steady behavior is observed.

To address the question of whether the quasi-periodic wave forms are likely even

for smaller Reynolds numbers, we simulate the case G = :25, T = 100 and {3 = ~

for wave numbers k = 0.:25, 0.23, 0.22 and 0.20. The evolution of the harmonic

modes for these cases is shown in Fig.6.15. Saturation of the harmonic modes is

observed for all wave numbers except k = 0.23 in which case we observe fluctuations

in the harmonic modes, most likely implying quasi-periodic wave forms. There is

a remarkable change in the behavior of the modes even for a small change in the

wave number as is evident from the cases k = 0.23 and k = 0.22. Thus, this

apparent quasi-periodic wave evolution is not completely a high Reynolds number

phenomena. Except, the band of wavenumbers for which this type of behavior is

exhibited increases with increasing Reynolds numbers.

An extensive parametric search in the range 5 5 G :S: 100, has been performed

to obtain the phase diagram shown in Fig.6.16. Starting from the linear cutoff

wavenumber kc, the nonlinear wave evolution has been obtained in steps of 0.01 up

to k = 0.10, for the thickness parameters G = 5,25,40,60,80&100. The promi-

nent findings of this section including the phase diagram shown in Fig.6.16 can be

summarized as follows:

l. The finite amplitude permanent waves are nearly sinusoidal for wavenumbers

slightly smaller than kc, and are broad-banded solitary humps for very small

wavenumbers. This is in agreement with the predictions of the longwave theory

based on the lubrication approximation, which is only valid for small Reynolds
 90

!= 36.76 != 37.56

!= 36.68 != 37.48

!= 36.6 != 37.4

!= 36.52 != 37.32

!= 36.44 != 37.24

!= 36.36 t= 37.16

t= 36.28 t= 37.08

t= 36.2 t= 37

t= 36.12 t= 36.92

t= 36.04 t= 36.84

t= 35.96 t= 36.76
Figure 6.14: Free surface profiles at different instants of time shown in a
three wavelength domain for G = 100, T = 100, f3 = 1r /2 and k = 0.20.
 91

(a) (b)
0.2 0.25

0.2
0.15
-0.15
~I
c
0.1 f I
()
I
- 0.1 t
-
o~u 0.05

0 0
0 10 20 30 0 50 100
t
(c) (d)
0.25 0.25

0.2

-0.15
c
I
()
- 0.1

40 60 20 30

Figure 6.15: The evolution of the spectral coefficients icn(t)l for G = 25,
T = 100, /3 = "Tr/2 and (a) k = 0.25; (b) k = 0.23; (c) k = 022; and (d) k = 0.20
 92

0.5

~ 0.4
I kM
_1-------i------·------
/ )I(

/
/
1 )I( :
)I( )I(

* I
0.1

10 20 30 40 50 60 70 80 90 100
G

Figure 6.16: Phase diagram for thin film flow over a vertical plane ob-
tained through full-scale computations. kc and kM are the cutoff and
maximum growth rate wavenumbers predicted by the linear theory. The
wavenumbers for which quasi-periodic evolution is observed are denoted
by '*'.
 93

numbers. The saturated nearly sinusoidal waves found close to the neutral

curve correspond to the slow moving short wave /I family of Chang eta[ ( 1993 ).

The long solitary waves for very small wavenumbers correspond to the fast

moving long wave 12 family of Chang eta! ( 1993). Chang eta! ( 1993) obtained

these stationary wave families based on boundary layer approximations.

2. For intermediate range of numbers, no permanent waveforms were observed.

Rather, a quasi-periodic wave evolution was observed. This type of wave

evolution was observed experimentally by Liu & Gollub (1994), in their spatial
evolution studies.

3. In the context of temporal stability, that is time evolution in spatial periodic

domains, no previous systematic study of this quasi-periodic evolutions have

been done, even though, they have been observed in the numerical simulations

of the Kuramoto-Sivashinsky equation by Hooper & Grimshaw (1985), and in

the numerical solution of the longwave evolution equation by Joo & Davis

( 1992). Our own numerical solutions indicate that this type of quasi-periodic

behaviour is also exhibited by the integral boundary layer system.

4. In the context of time-periodic evolution is space (spatial stability), Cheng

& Chang (1995) discuss the secondary instabilities in thin film flows. The
formation of the permanent waveform is considered as the primary instabil-

ity in thin film flow. This permanent wave travels downstream at a constant
speed, but succumbs to the secondary instabilities, namely, the side-band and

subharmonic instabilities. Due to these secondary instabilities, complicated

coalescence occur between neighboring peaks and the evolution towards soli-

tary waveform occurs. Ours is a space-periodic system evolving in time. Thus,

 94

the quasi-periodic behavior is a manifestation of the side-band instability, al-

beit, a temporal one.

5. The quasi-periodic region is not contiguous, and for some intermediate wavenum-

bers permanent waveforms reappear. This is due to the secondary subhar-

monic instability.

6. 7 Spatio-temporal evolution of the thin film instability

In the previous section, we studied the temporal evolution of the surface wave in-

stability in a periodic domain. In reality, the fluid domains are not restricted by

periodicity, and the surface wave evolves both spatially as well as temporally, al-

lowing complex wave interactions such as wave-mergers and wave-splitting. Spatia-

temporal evolution of the instability and the wave interaction processes are thus

studied in this section.

6.7.1 Long domain with periodic boundary conditions

The parameters are chosen as G = 5, T = 100, (3 = 1r /2, and the length of the

stream wise periodic domain is set to 20,\M, where ).M is the maximum growth

rate wavelength predicted by the linear theory. In the first 1/20th of the domain,

i.e., 0 :=:; x :=:; ).M, the equilibrium wave profile (Fig.6.11b) computed previously is

imposed as the initial condition and the rest of the domain is undisturbed. In the

absence of external forcing, the natural waves that evolve downstream of the wave

inception are of wavelength >w (Chang 1994). This type of an initial condition

helps us understand how these natural waves evolve further downstream.

The dispersion of the initially imposed wave is shown in Fig.6.17( a). The primary

surface wave instability being convective in nature, the initial wave is transported
 95

(a)
-I I I I I
I I

~
250
~

~
.-::.

200 "'::~

;:..
.:::.

150

-
100 :§-
A
~~
;;!:>
~
~..___
so ~~
:?'v
v
~
~--
0 R: I I I I I I I

0 2 4 6 8 10 12 14 16 18 20

Figure 6.17: Spatio-temporal evolution in a long periodic domain for

G = 5.0, T = 100.0 and (3 = 1rj2. The length of the domain 20,\M. The =
free surface profiles shown are in intervals of D.t = 5 from : (a) t = 0 to
t = 250; (b) t = 250 to t = 500; (c) t = 750 to t = 1000; and (d) t = 5240 to
t = 5490.
 96

downstream. In the process, however, the initial wave disperses into several capillary

waves. The front running wave quickly attains the shape of a solitary hump and

travels downstream at a constant velocity and without further generation of waves.

Each of the trailing waves as they travel downstream evolve into a solitary waves

with sharp downstream slope and a gently sloping tail (Fig.6.17b ). Thus, very far

away from the source, the stable waveform is a solitary waveform. The domain being

periodic, the wave leaving the domain at the right enters at the left. The speed of the

wave is seen to be proportional to the wave amplitude. This fact was also observed

experimentally by Liu & Gollub (1994). The larger amplitude solitary waves travel

faster and run into small capillary waves (Fig.6.17c). There doesn't appear to be

any repulsion between the two waves as they come close to each other. At the

end of the wave coalescence, a still larger amplitude wave is formed, which travels

downstream leaving behind a quiescent interface devoid of any small scale ripples.

The coalescence between the large amplitude solitary wave and the relatively smaller

amplitude capillary ripples is thus an inelastic one. This type of wave interactions

go on for a long time, until the only type of waves seen on the interface are of the

solitary-wave type (Fig.6.17d). The solitary pulses try to keep a certain distance

from each other. If two solitary pulses get too close, then they start repelling each

other, and this is the cause of the quasi-periodic behaviour observed previously in the

temporal stability analysis (§6.6). At the same time, they cannot be too far apart,

since small capillary ripples are bound to appear on the interface. These capillary

ripples travel slower than the larger amplitude solitary pulse and eventually merge

with the solitary pulse. Two solitary pulses of comparable amplitude, however, do

not merge with each other. About t = 2000, a solitary wave train with 11 solitary

pulses is formed in the domain. Even after integrating in time up to t = 5490, the
 97

(b)
I I T I I I I

500 k\.-\..._-_-_-_---,-.:-7.-<\ ~

~"~=~~\: \-( \ ~

--?-- /'"'\
1-----,-.:7 \ '::: '---
/"\ '---

450
~ ~

-
350

I I I
' '
0 2 4 6 8 10 12 14 16 18 20
x/).M
 98

(c)

0 2 4 6 8 10 12 14 16 18 20
x/ >w
 99

(d)

0 2 4 6 8 10 12 14 16 18
x/>w
 100

number of pulses is found to remain the same, namely, 11 in a domain of length

20,\M (Fig.6.17d). The only changes in the domain are the relative spacing between

the different solitary pulses. The solitary pulses are continuously interacting with

each other to arrive at an appropriate relative distance. However, the temporal

saturation of the solitary wave train is not complete in the numerical integration

time 0 ::; t ::; 5490. A much longer time of integration may have to be done, before
complete temporal saturation can be obtained.

For the same simulation parameters as above, namely, G = 5, T = 100, (3 = 7r /2

and a streamwise periodic domain with length L = 20,\M, an exponential pulse of
the form h(O, t) = 1 + 0.05e(-O.Ol(x-tfl is initially imposed, and its spatio-temporal

evolution is studied. The initial dispersion of this disturbance and the evolution into

solitary waveform is shown in Fig.6.18( a). The spatio-temporal evolution is qual-

itatively same as that seen previously, with complicated wave-splitting and wave-

mergers taking place continuously, until a solitary wave-train is formed (Fig.6.18b).

The notable difference from the previous simulation, is the presence of 13 solitary

pulses in the wave-train in a domain of size 20,\M· Just as in the previous case, after
the formation of a solitary wave train with 13 solitary pulses (around t = 2000),

further integration in time (up to t = 4355) is only found to change the relative
spacing between the pulses, without changing the number of solitary pulses in the

domain. Thus the relative spacing between the solitary pulses, or the natural non-

linear wavelength, is found to be weakly dependent on the initial conditions.

The findings of this subsection can be briefly summarized as follows:

1. Far downstream from the source, the solitary waveform is the stable waveform.

2. The wave speed is proportional to the amplitude. This has also been observed

experimentally by Liu & Gollub (1994).

 101

(a)
I -I -, I I I I
' '
250 '-::::"'
/(
::..----- ./
/
:/'..
::::.. ;r.,.
:::::..
;:....
7
?'
200

'-
~
---~
-?"
*
150
~
- ?
h
-;:;;.

~
-:.;.....~
100
:;:...,
~
~

50

0
I I I I _l_ I
' '
0 2 4 6 8 10 12 14 16 18 20

Figure 6.18: Spatio-temporal evolution in a long periodic domain for

G = 5, T = 100 and j3 = f· The length of the domain L = 20,\M· The
initial disturbance on the free surface is an exponential disturbance of
2
the form h(x,O = 1+0.05*e(-O.OI(x-tl l. The free surface profiles are shown
in intervals of b.t = 5 from: (a) t = 0 tot= 250; and (b) t = 4105 tot= 4355.
 102

(b)

0 2 4 6 8 10 12 14 16 18 20
 103

:3. When the large amplitude wave approaches the slow moving small amplitude
wave, there doesn't appear to be any repulsion between the two. At the end
of the merger, the large amplitude wave travels downstream leaving behind an
interface devoid of waves. This has also been observed experimentally by Liu
& Gollub (1994).

4. There appears to be a natural nonlinear wavelength the system chooses to

have at the end of the nonlinear wave evolution. However, the relative spacing
between the solitary pulses is also dependent on the initial conditions.

6.7.2 Wave breaking in thin film flows

The question of wave breaking in laminar thin film flows is still an open one, and

no conclusive theoretical, experimental or numerical evidence to this effect exists

to date. The tendency for wave breaking arises from the thickness dependence of
the local phase speed, with the crests travelling faster than the troughs resulting
in steepening of the wave. The wave interaction processes described in this section
could, however, inhibit the wave breaking tendency.

Theoretical interest has focused on determining if the approximate evolution

equations such as the Kuramoto-Sivashinsky equation and the longwave evolu-
tion equation of the Benney type admit wave breaking solutions. The Kuramoto-
Sivashinsky (KS) equation can accurately model evolution of the surface wave insta-

bility in thin film flows in the limit of small G and very large surface tension. The
solutions of the KS equation are, however, smooth at all times and do not exhibit
wave breaking tendency. Rosenau & Oron (1989) hypothesized that this is due to
the over prediction of the effect of surface tension in the KS equation and proposed
some a modification to the curvature term in the KS equation. They retained the
 104

higher order terms in the denominator of the curvature and called their modified

KS equation as the Regularized Kuramoto-Sivashinsky equation (RKS). Through

numerical experiments, they show that in certain range of parameters the RKS equa-

tion exhibits wave breaking tendency. Joo, Davis & Bankoff (199la) demonstrate

that the longwave evolution equation of the Benney type (6.37) also exhibits the

numerical blowup behaviour. However, both KS and longwave evolution equations

are based on the assumption of finite wave amplitude and spatial gradients and lose

their validity much before the onset of wave breaking. Thus, it would be interest-

ing to determine if the numerical blowup of these evolution equations really implies
wave breaking. The numerical parameters chosen by Rosenau & Oron (1989) are

simulated using our full scale model and the results are shown in Fig.6.19. We find

that the harmonic modes eventually saturate implying that finite amplitude per-

manent waveforms are the stable solution for these simulation parameters. This is

contrary to what has been observed by Rosenau & Oron (1989) in their solution of

the RKS equation. The RKS thus appears to have no physical relevance.

6.7.3 Long domain with non-periodic boundary conditions

Up to now, all the domains considered are periodic in streamwise direction, and what

is studied is strictly the temporal stability of the thin film flows. In the experimental

studies of thin film instability, however, a periodic disturbance either in the form of

pressure fluctuations or film thickness perturbations are imposed at the inlet and

the evolution of this disturbance in the stream-wise direction is measured (spatial

stability analysis). For infinitesimal disturbances (linear stability), temporal and

spatial stability give equivalent critical conditions. When the disturbances are of
finite amplitude (nonlinear stability) this is not always the case, and additional
 105

(a)

0.10

0.05

0.00
0 5 10 15 20 25 30 35 40 45
t
(b)
1.6

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
xj>.

Figure 6.19: Rosenau & Oron (1989)'s wave breaking simulation with
simulation parameters G = 16.0, S = 51.63, (3 = 1r /2 and k = 0.2783. The
initial condition is h(x,O) = 1 +0.5cos(4kx). (a) Evolution of the spectral
coefficients /en( t ); and (b) Final wave shape.
 106

information can be obtained by considering both approaches. In this section, the

spatial stability of the thin film flows, in a manner akin to the physical experiments,

is done numerically. For this purpose, the experimental conditions of Liu & Gollub

(1994) are simulated. Liu & Gollub (1994) studied the solitary wave dynamics of

thin film flows. They use a 54% by weight aqueous solution of glycerin and impose

periodic pressure fluctuations at the inlet. The angle of inclination f3 = 6.4°; and the

mean film thickness is ho = 0.12789cm. These experimental conditions correspond

to G = 520.32, S = 676.65, and T = 84.10. The thin film is initially fiat and

periodic disturbances are imposed at the inlet in the following manner:

h(O,t) = 1 +bsin(wt), (6.53)

where w and b are respectively, the non-dimensional angular frequency and ampli-
tude of the external periodic forcing. The velocity boundary conditions at the inlet

are of the Dirichlet type and are imposed based on the lubrication approximation

and in such a way. that the continuity equation is satisfied.

2

u(O,y,t) Gsinf3 (h(O,t)y- ~)

(6.54)
1 8hy 2
v(O,y,t) = 2 at2·
At the exit x = L, Sommerfeld Radiation Boundary Conditions in a manner sim-

ilar to that proposed by Orlanski (1976) are imposed to let the waves leave the

computational domain with minimum reflection.

First we present our results corresponding to a 1.5 Hz forcing frequency (Fig. 3

in Liu & Gollub(1994)). The amplitude of the disturbance is set to b = 0.05. The

film thickness h(x, t) at various instants of time starting from t = O.Os tot= 2.6s in
steps of 0.104 is shown in Fig.6.20(a). The primary instability being convective, the
 107

(a)

~-
I I (l
_,(
I~
T
' I r

2.5 ~~· ~r..l' A~

~I k --_; ',.'.'
1
A
1--- ' A ~·· -::/\li,,
~ ~ /~ :~~:'"·
~ :;
~,
I' /}f..~'''
,'II'
2.0
~~· ~:''
~' )ln:'
L-------.. ' _/In ,Y r
f.-..._
....... 1.5 ~-
- -;:)" -
u ~
Q)
(/)
._..
~=-v-
Q) ~v -
E ~·
:;::; 1.0 ;:__
~---
----
~

~-----

0.5 ~-
f-

0.0
I I I I I
' '
0 20 40 ~ 00 100 1~ 140 1~
Downstream distance (em)

Figure 6.20: Spatio-temporal evolution with f = 1.5H z periodic forcing

at the inlet. G = 520.32, S = 676.65, (3 = 6.4deg and corresponds to the
experimental conditions of Liu & Gollub(1994) (Fig.3). A sinusoidal per-
turbation of the form h(O, t) = 1 + 0.05 * sin(2.4546t) is imposed at the inlet.
The free surface profiles shown are shown in intervals of 0.104.s from: (a)
t = O.O.s to t = 2.6.s; (b) t = 2.6 to t = 5.2.s; and (c) t = 5.2.s to t = 7.8.s.
 108

disturbance is transported downstream by the mean flow. According to the linear

stability theory the cut-off frequency for the onset of instability is we = 37.05. Since
the imposed frequency is much smaller than the cut-off frequency, the most likely
waveform is a solitary waveform, and that is what is observed. The amplitude of the

disturbance quickly grows downstream and the wave becomes a symmetrical, with
a sharp upstream slope and a gently trailing tail. Several small amplitude capillary

ripples are formed downstream of the primary hump. In Figs.6.20(b) & 6.20( c) the
subsequent free surface profiles in steps of 0.1.04s are shown. The wave traveling in

the front grows in amplitude, travels faster and ultimately leaves the computational
domain. After the initial growth, the wave profile doesn't appear to change much in

its journey downstream. Thus the gas-liquid interface is filled with solitary humps

having tear drop profile and there appears to be spatial saturation in the wave form.

Further downstream however, spatio-temporal chaos and transverse 3-D instabilities

will set iri. The free surface profile given by Liu & Gollub(1994) is compared against

numerically computed free surface profile in Fig.6.21. The agreement between the
experimental and numerical profiles is not very good in the first half of the domain.

This is most probably due to a difference in the amplitude of the inlet excitation. Liu

& Gollub (1994) do not mention the amplitude of the disturbance imposed at the
inlet and we set it arbitrarily to be 8 = 0.05. The saturated wave profiles as given

in the later part of the domain (Fig.6.2lc), however, show very good agreement.

The significant difference is that the amplitude of the wave as predicted by the
numerical calculations is slightly more than the experimental value. The capillary

ripples downstream of the primary hump are captured quite well by the numerical

solution.
 109

(b)
I I I I I I I

20 40 60 80 100 120 140 160

Downstream distance (em)
 110

(c)
I I I I I

5.0 L___ 1 _ __...JIL___ _..t.__

___L I _ ___J__ __ L ._ __LI_ __...JIL___ _JI_
___J

0 20 40 60 80 100 120 140 160

Downstream distance (em)
 lll

(a)

'
.c: 1.0

95

' 61
1.4
1.2
1.0
0.8
0.6 0. 7 ':-_J_--l._..L._--L---l_ _L__L_L.--.L-l
90 100 110 120 130
Downstream distance (em) 85 95 lOS liS 125 IJS
(c) I .6 r--r----r-;r--,--,----,--,--,----,-----,
(C)

21 JL: I
I.J
' 61
1.4
1.2
.c
.....
.c.
0

I .0
1.0
0.8
0.6
120 130 140 150 160 140 ISO 160 170
Downstream distance (em)
Downstream Distance <em)

Figure 6.21: G = 520.32, S = 676.65, j3 = 6.4deg and corresponds to the

experimental conditions of Liu & Gollub(1994)(Fig.3). A sinusoidal per-
turbation of the form h(O, t) = 1 + 0.05 * sin(2.4546t) is imposed at the inlet.
The wave profiles shown on the left are obtained computationally and
the ones shown at the right are those measured by Liu & Gollub(1994)
(Fig.3).
 112

The experimental conditions simulated next are same as the above except that
the inlet forcing frequency is increased to 4.5 Hz and the amplitude of the excitation

is taken to be 8 = 0.01. The free surface profiles in intervals of 0.104s are shown

in Figs.6.22(a,b,c). The inlet forcing frequency controls the spacing between the

solitary humps. In this case the frequency is quite high and the solitary humps are

too close to each other and never settle down to achieve steady traveling waveforms.

This is analogous to the quasi-periodic behavior we observed for some wave numbers

in the study of temporal instability. The experimental and numerical wave profiles

are compared in Fig.6.23. Just as in the previous case, the agreement is not very

good in the first part of the domain (Fig.6.23a). This either due to a difference

in the forcing amplitude at the inlet or due to a difference in the way the flow-is

excited at the inlet. Liu & Gollub(1994) impose fluctuations in the pressure, where

as in the numerical simulations the film thickness is excited. In the later part of the

domain the agreement is very good (Fig.6.23b,c). Since the solitary waves are too

closely packed we do not observe spatial saturation.

Lastly, keeping the same experimental conditions as above, but increasing the

inlet forcing frequency to 7 Hz, we obtain nearly sinusoidal film profiles as shown
in Fig.6.24. The waves are closely packed, nearly sinusoidal and symmetrical. The

forcing frequency in this case is close to the cut-off frequency, hence the almost

sinusoidal wave profile. The comparison with experiments is very good in this case

also (Fig.6.24(b,c)).
In §6.5, the experimental conditions of Kapitza & Kapitza (1949) were simu-

lated numerically in a periodic domain. For one of the experimental conditions

(Figs.6.8,6.9), our full-scale numerical simulations predict a quasi-periodic wave-

forms in contrast to the periodic waveforms reported by Kapitza & Kapitza (1949).
 113

(a)
I I I I I I I I

~~~
'-"'"'
2.5 ~~

V r'~V
"-" \r v
~
v
"-"
'{':
2.0 '{
A
'-"'"' A v
~

v
'A
"~ v
oA

..-.. 1.5
0
"'" A -
(]) ~

(/)

---
(])
v'l,~

-'\.,.
-
E
:;::: 1.0 "-" -

0.5 t-

0.0

I I I I I I
' '
0 20 40 60 80 100 120 140 160
Downstream distance (em)

Figure 6.22: Spatio-temporal evolution with f = 4.5H z periodic forcing

at the inlet. G = 520.32, S = 676.65, f3 = 6.4deg and corresponds to the
experimental conditions of Liu & Gollub(1994) (Fig. 7). A sinusoidal per-
turbation of the form h(O, t) = 1 + 0.01 * sin(7.3638t) is imposed at the inlet.
The free surface profiles are shown in intervals of 0.104s from: (a) t = O.Os
to t = 2.6s; (b) t = 2.6s to t = 5.2s; and (c) t = 6.656s to t = 9.4s.
 115

(c)

uQ)
(/)
;a.o
E ..--....~~ '"" ,
:;:;

60 80 100 120 160

Downstream distance (em)
 116

(a)

1.21 0 1.0
1.0 .c.
0.8 '.c. 0.8
0.6L-------~~--~----------~
50 60 70 80 90
Downstream distance (em)
(b)
(b)

1.21 0
1.0
1.0 .c.

0.8
'.c. 0.8

0.6L-------------~~~~~~~
90 1 00 11 0 1 20 1 30
Downstream distance (em) 135
(c)
(C)

0
1.21 .c.
1.0
'.c. 0.8
0.8
0
·Y2o 130 140 1so 160 150 160 170
Downstream distance (em)
Downstream· Distance <cml

Figure 6.23: G = 520.32, S = 676.65, j3 = 6.4deg and corresponds to the

experimental conditions of Liu & Gollub (1994)(Fig.7). A sinusoidal per-
turbation of the form h(O, t) = 1 + 0.05 * sin(7.3638t) is imposed at the inlet.
The wave profiles shown on the left are obtained computationally and
the ones shown at the right are those measured by Liu & Gollub (1994)
(Fig.7).
 117

_9.0
1Jl
~8.5
E
o+=> 8.0

7.5

0 eo ao 100 120
Downstream distance (em)

1.1
1.0 0 1.0
h 0.9 .c.
ho o.a
.....
..c: 0.8
0.7
0.6
eo 70 ao so 100
Downstream distance (em) 65 75 85 95 105
Downstream Distance <cml

Figure 6.24: G = 520.32, S = 676.65, j3 = 6.4deg and corresponds to the

experimental conditions of Liu & Gollub (1994) (Fig.8). A sinusoidal
perturbation of the form h(O, t) = 1 + 0.05 * sin(11.4548t) is imposed at the
inlet. (a) film thicknesses in intervals of 0.104s from t = 7.6s to t = 9.6s;
(b) numerically computed free surface profile; and (c) experimentally
reported film thickness profile (Fig.8 in Liu and Gollub (1994)).
 118

The same experimental conditions are now numerically simulated in a long non-

periodic domain and the spatial evolution of the disturbance in a manner similar

to the physical experiments is simulated and the results are shown in Fig.6.25. No

spatial saturation is observed in the stream-wise direction and the wave profile is

quasi-periodic consistent with the temporal stability predictions.

6.8 Concluding Remarks

A numerical procedure based on the Finite Element Method has been developed to

study the interfacial instabilities in isothermal thin film flows flowing down vertical

or inclined planes. The numerical code is not very expensive in terms of compu-

tational time. A mesh of 2001 X 11 was used in the simulation of Liu and Gollab

(1994)'s experimental results. A simulation time of 30 non-dimensional time units

which correspond to 10 real seconds were sufficient to observe spatial saturation.

This took about 30000 time steps and approximately 30 CPU hours on a SPARC-10.

All the simulations presented in this chapter have been performed on a SPARC-10

computing machine.

The work presented in this chapter is the first comprehensive study of thin film

instability based on the solution of the transient N a vier-Stokes equations ( Chippada,

Ramaswamy & Joo 1995). Salamon etal (1994) did a detailed study of thin film

instability, but a priori assume steady traveling waveforms and solve the steady state

N a vier-Stokes equations in a moving frame of reference. However, the stability of the

waveform can only be ascertained through the solution of the initial value problem.

Extensive numerical simulations based on the direct solution of the N a vier-

Stokes equations reveal some interesting features of the thin film instability. In

concurrence with what is known about thin film instability based on linear and ap-
 119

--
0
Q)
(/)
;o.7s
-
E

0.65

0.60 ' - - - - ' - - - - - ' - - - - - ' - - - - - ' - - _ _ _ , J ._ _

..__~_ _ j __ __J__ __...J

0 2 4 6 8 10 12 14 16 18 20
Downstream distance (em)

Figure 6.25: The experimental condition of Kapitza & Kapitza (1949)

simulated in a long non-periodic domain. G = 60.0, S = 4410, f3 = 1r /2
and w = 5.707. At the inlet the free surface is perturbed as follows:
h(O, t) = 1 + 8 sin (wt). The free surface profiles in steps of 0.033sec are
shown.
 120

proximate nonlinear theories, finite amplitude waveforms are the stable solution for

wavenumbers smaller than the linear cut-off wavenumber kc. For wavenumbers close

to kc, the waveforms are nearly sinusoidal. For wavenumbers much smaller than kc,

the waveforms are solitary-wave-like. This transition from nearly sinusoidal wave-

forms to solitary waveforms seem to pass through a quasi-steady regime, in which

the spatial harmonic coefficients are in a state of constant fluctuations. The phase

boundary de-lineating this regime has been obtained through extensive numerical
simulations.

The spatial stability analysis of the thin film instability has also been studied

by considering a very long domain with periodic forcing at the inlet and absorbing

boundary conditions at the exit. Excellent agreement with the experiments of Li-u

& Gollub (1993) has been obtained. Depending on the frequency of excitation,

the waves formed downstream are either nearly sinusoidal or solitary-like or quasi-

periodic (in streamwise direction).

Due to the amplitude dependence of the wave speed complex wave interactions

are likely to occur on the gas-liquid interface. Waves with larger amplitude travel

faster and coalesce with smaller waves on the wave. This wave interaction is found
to be completely inelastic and the resultant wave grows further in amplitude and

travels downstream leaving behind a nearly fiat interface. However, there appears

to be a natural wavelength the system tries to have in the solitary wave regime. The

resultant wavelength downstream is also weakly dependent on the initial condition.

A powerful numerical technique has been developed and used successfully to

under stand many aspects of the surface wave instability in the thin film flows. A

natural extension of the numerical procedure would be to do three-dimensional sim-

ulations. Very far away from the source the nonlinear waves are three-dimensional.
 121

A numerical procedure based on the direct solution of the three-dimensional N a vier-

Stokes equations would help us understand this nonlinear wave regime.

 122

Chapter 7

Turbulent Fluid Flow Over Streamwise Periodic

Artificial Bedforms

7.1 Introduction

An erodible bed composed of non-cohesive sediment particles with turbulent fluid

flow on top of it is rarely flat, e.g. river beds and desserts. More common is

the presence of bedform features such as ripples, dunes and antidunes collectively
referred to as bedforms. Ripples and dunes are found to exist in flows with Froude

number less than 0. 7. Ripples occur at length scales independent of the fluid depth,
and depend only on the local bed properties such as grain size. Dunes occur at

larger length scales that are dependent on the fluid depth. Antidunes are found to
exist in flows with Froude number greater than 1.0. Unlike ripples and dunes, in the
case of anti dunes, strong coupling exists between the free surface (air-water) waves
and the bed waves. A more detailed description and classification of bedforms can

be found in ASCE (1966) and Kennedy (1969).

The formation of the bedforms is due to a complex interplay between the fluid

and the bed. In particular three important features can be isolated, namely, the
turbulent fluid flow, the sediment transport and the bed morphology. These three
phenomena are coupled and collectively influence the formation and migration of
bedforms. For example, the local sediment concentration is dependent on the fluid
velocities, stresses and turbulence intensities. The sediment concentration field in
turn determines the erosion and deposition of sediment particles at the bed thus
 123

controlling the formation of bedforms. The shape of the bed (bed morphology) in
turn influences the flow and turbulence field.

Several theories have been proposed which try to explain the formation of bed-

forms as a problem in stability. The most important and the most difficult aspect

in these models is the determination of the turbulent fluid flow over the bedforms.

The earliest models assumed potential flow (Kennedy 1963). A more realistic model

for the fluid flow, namely, real viscous fluid flow with vorticity, was introduced

by Engelund (1970) and Fredsoe (1974). However, Engelund (1970) and Fredsoe

( 197 4) assume a constant turbulent eddy viscosity which is determined from the

undisturbed fluid flow. Richards ( 1980) extended the theory of Engelund ( 1970)

to incorporate a spatially varying eddy viscosity. They determine the turbulen.ce

eddy viscosity by solving the turbulence kinetic energy equation, but assume lin-

early varying length scale away from the bed. All these theories lose their validity

when the amplitude of the bedform becomes finite. The fluid flow is known to

undergo separation resulting in the formation of recirculation zones downstream of

the dune peak (see Fig.7.1). Widely varying length scales exist in different parts

of the fluid domain, and no complete analytical theory exists that can successfully

model this type of complex turbulent flow. McLean & Smith (1986) developed a

semi-analytical theory which uses different solution techniques in different parts of

the fluid domain. The internal boundary layer is solved using boundary layer as-

sumptions and asymptotic expansions and match it with the external wake solution
calculated based on the wake solution behind a circular cylinder. Nelson & Smith

(1989) extend the theory to include the effects of short wavelengths and recirculation

zone. These semi-analytical theories are not universal and involve some parameters

which need to be fine-tuned to fit the field conditions.

 124

~g
AIR

WATER

---x

Figure 7.1: Schematic description of fluid flow over periodic bedforms

 125

There remains a need for obtaining accurate mean flow and turbulence charac-

teristics of the fluid flow over bedforms. It is believed that a good understanding

of the fluid flow mechanisms are essential to improve the existing bedform stabil-

ity models. Therefore, the goal of this work is to develop a numerical procedure

capable of predicting the turbulent fluid flow over artificial fixed bedforms. In addi-

tion, an accurate computation of fluid flow over bedforms can be used for improved

prediction of effective roughness of channels which would be of use to hydrologists.

Due to the limitations in the existing turbulence theories. the only means of

obtaining detailed mean !low and turbulence characteristics are experimental and
numerical. Raudkivi (1966). Lyn (1993) and Nelson, McLean & Wolfe (1993) are

some of the experimental studies that have measured the mean flow and turbulence

characteristics in fluid flow over bedforms. The presence of sediments in suspension

complicate the measurements. In addition if the bedforms are erodible, the bedforms

migrate slowly downstream due to the continuous erosion and deposition of sedi-

ments particles at the bed and the fluid flow is not stationary in the strict sense. To

overcome these difficulties, the laboratory plume studies measure the mean flow tur-

bulence characteristics for fluid flow over artificial, fixed and non-erodible bedforms,

whose size and shape corresponds to naturally occurring bedforms. Numerical stud-

ies invariably solve the Reynolds Averaged Navier-Stokes (RANS) equations instead
of the original Navier-Stokes equations due to the limitations imposed by the exist-

ing computational architecture. Due to the presence of variety of flow regimes such
as internal boundary layer and recirculation zone, the mixing length closure mod-

els are inaccurate and higher order turbulence models need to be used. Richards

& Taylor (1981) used a one-equation turbulence kinetic energy closure model and
specified the length scale as linearly varying away from the wall. Mendoza & Shen
 126

(1990) use a two-equation k- E model along with an Algebraic Stress Model to

compute the turbulence field. The presence of recirculation zone requires that at

the minimum a two-equation closure model be used to compute the turbulence field
with reasonable accuracy.

To our knowledge, all the previous numerical studies make the rigid-lid approx-

imation, i.e., they replace the unknown air-water interface with a known rigid, free-

slip wall. This eliminates the free boundary simplifying the problem considerably.
The rigid lid approximation gives reasonably accurate results for flows with small

Froude numbers ( F « 1). For higher Froude number flows it is no longer a rea-

sonable assumption since the free surface deflection is not insignificant anymore.

In this work, we develop a numerical procedure capable of handling the air-wat~r

free surface in a simple manner. Thus the numerical algorithm developed in this

work is capable of simulating high Froude number flows as well. The turbulence is

modeled through a k- € closure model coupled with the Boussinesq Eddy Viscosity

hypothesis and the RANS equations are solved using the Galerkin Finite Element

Method. The emphasis of this work is on the accurate determination of the mean

flow and turbulence field. But, we also compute the sediment load using the existing

empirical models to understand the influence of the flow field on the sediment load

characteristics.

7.2 Problem Definition and Fluid Flow Modeling

The fluid flow geometry is shown schematically in Fig.7.1. The bottom boundary

consists of bedforms which are periodic in the stream-wise direction with wavelength
L. Even though, the bedforms are known to migrate downstream, the phase velocity

of the bedform migration is very small compared to the stream-wise fluid velocities.
 127

Hence. we assume that the bedform geometry does not change with time, i.e. the

bed is assumed to be non-erodible and determine the final steady-state equilibrium

velocity and turbulence field. The channel bed is inclined at a small angle () with

the horizontal. Thus. the flow is driven by gravity . .r and y axes are as shown in

the figure and b( x) represents the bed-liquid interface. h( x, t) is the height of the
liquid-air interface from the x-axis.

The Reynolds number based on the mean fluid depth is usually much greater

than the cut-off Reynolds number for transition to turbulence in open channel flows.

Hence, the fluid flow occurring in nature over bedforms is turbulent in nature.

This turbulent fluid flow is modeled through the Reynolds Averaged N a vier-Stokes

equations (RANS), written as:

au au au - - 1 op 2
- + u- + v - = gsmB- - - + (v + vr) \7 u (7.1)
at ax ay pox

-av
at
+ ov av
u- + v- =
ax 8y
1 8p
-g cos() - - -
p oy
+ ( 2
v + vy) \7 v (7.2)

(7.3)

The turbulent eddy viscosity vy is determined through the two-equation k-t closure

model, which has been explained in chapter 4.

The boundary conditions at the bed are the no-slip conditions given by:

u = 0 }
at y = b(x). (7.4)
v = 0

The boundary conditions at the free surface are given by:

 128

aun
-p+ 2f-t-
an
au_,. aun)
1-l ( - + - 0 at y = h(x, t). (7.5)
an aT
ah ah
-+u- 0
at ax
The first two stress conditions are imposed in the solution of the Navier-Stokes

equations and the third condition, namely, the kinematic free surface condition is

used to update the free surface location every time step. The moving boundary

problem is handled through the Arbitrary Lagrangian-Eulerian (ALE) procedure

which has been explained in chapter 3. The boundary conditions in the streamwise
direction are the periodic conditions given by:

~(X = 0) = ~(X = L) (7.6)

where, ~ stands for u, v, p, k, t and h. Even though we are interested only in

the steady state flow and turbulence field, since the procedure used is a transient

one, we also need to impose appropriate initial conditions. The initial velocity field

is specified based on the universal log-law and in such a way that it satisfies the

divergence free condition (Eq.7.3). The initial conditions for k and t are imposed

so that they vary linearly in the following manner (Alfrink & van Rijn 1983):

y- b(x) ) u;(x)
k(x,y,t = O) = ( 1 - h(x,t = 0)- b(x) jC:. (7.7)

y- b(x) ) u~(x) (7.8)

t(x, y, t = O) = ( 1 - h(x, t = 0)- b(x) --;;;;;-

As explained in chapter 5, the Galerkin Finite Element Method along with a Chorin-

type projection scheme is used to solve the RANS.

 129

/\ll the results presented in this chapter have been obtained usmg a 61 x 11
mesh on SPARC-10 computing devices. A typical Finite Element Mesh is shown in
Fig.7.2.

7.3 Sediment Thansport Modeling

Though, the primary focus of this work is the computation of the flow and turbu-

lence field, a numerical model capable of simulating the sediment transport is also

implemented in the code. and the details of this sediment transport model are briefly
explained in this section.

The transport of sediment particles, which could be sand, silt, gravel or even
boulders, downstream by the fluid flow is termed sediment transport (Henders()n

1966). The sediment materials can be broadly divided into two classes, namely,

cohesive and noncohesive materials. In cohesive sediment materials there exist ad-

hesive forces between the sediment particles, e.g., clay. In noncohesive sediment

materials, no attractive forces exist between the sediment particles, and only this

type of sediment transport is modeled in this work.

For the sake of argument, consider a sediment bed with stagnant fluid on top

of it. No type of sediment transport is expected in this state. Let the fluid be set

in motion through some means such as a pump, or by tilting the plume. It will be

found that for low fluid velocities, still there would be no movement of the sediment

particles from the bed. As the fluid motion is continuously increased, at some

point we would find a perceptible motion of the sediment particles. Thus, a certain
threshold fluid velocity or shear stress is required to set the sediment particles in

motion. The most commonly used expression for this thres~old condition is that

given by Shields (1936). When the sediment particles start to move, initially, they
 130

Figure 7.2: Typical Finite Element mesh

 131

roll and slide along the bed and are in continuous contact with the bed. With

further increase in the fluid velocity, the particles undergo occasional jumps called

saltations. This type of transport of the sediment particles where thev are either

rolling, sliding or saltating is called the bed load transport. A further increase in

the fluid velocities causes the sediment particles to go into suspension. The weight
of the sediment particles is supported by turbulence forces and the transport of the

sediment particles in suspension is termed the suspended load transport. The bed

load transport is influenced to a large extent only by the fluid friction at the bed,

whereas. the suspended load transport JS influenced by both, the bed friction and

the fluid velocities in bulk flow.

Many of the mechanisms of the sediment transport are not very well understood

and the theories proposed to compute the sediment transport in general tend to be

empirical in nature. In general the naturally occurring sediment beds are composed
of sediment particles that vary widely in sizes and shapes. The process of entrain-

ment of the sedim~nt particles are stochastic in nature, and no clear understanding
exists about the influence of the suspended sediment particle on the fluid flow and

turbulence characteristics. Thus, most of the models proposed to compute the bed

load transport are only heuristic in nature. We chose a sediment transport model

proposed by van Rijn (1984a,b,c) to compute the sediment transport. The hydrody-

namics code developed in chapters 3, 4 & 5 is coupled with the sediment transport

model and the influence of the bed morphology on the sediment load distribution

is studied. With the availability of better sediment transport models, they can be

linked to the hydrodynamics code in a straight forward manner.

The sediment transport model developed by van Rijn (1984a,b,c) is briefly ex-

plained next. In this sediment transport model, the bed load transport is computed
 132

using an empirical equation which is a function of the local bed shear stress and
the suspended sediment concentration is modeled through a continuum advection-
diffusion transport equation. First, the following three nondimensional parameters
are computed:
1/3
D = D (s- 1) g (7.9)
• 50 [ 2 ]
v

T= (7.10)

(7.11)

D. is called the particle parameter, T is called the transport parameter and Z js

called the suspension parameter. D 50 is 50th percentile sediment particle diameter,

s = Psi p is the specific density of the sediment material, g is the gravitational

acceleration, v is the kinematic viscosity of the fluid, u. = ;;:JP is the friction

velocity at the bed, Tw is the shear stress at the bed, u.,cr is the critical bed shear

velocity needed for the initiation of the sediment movement as proposed by Shields
(1936), and "' is the von Karman constant usually set to 0.40. The Shields critical

bed shear velocity u•,cr is computed as follows:

D. ::::; 4, 8cr = 0.024 (D.)- 1

4 <D. ::::; 10, 8cr =0.014(D.)- 0 "64

10 < D. ::::; 20, 8cr = 0.04 ( D.)- 0 ·10 (7.12)

20 <D. ::::; 150, 8cr = 0.013 (D.) 0 .29

150 <D., ecr = 0.055

where, 8cr = u;,cr/ ( ( s - 1) gDso).

 133

The bed load transport S; is computed from the following empirical equation
(van Rijn 1984a):

yz.1
S; = 0.053 [(s -l)y) 05 D~ 0 5
-
.
D0.3
(7.13)

The bed load transport is influenced by the local slope of the bed also. For instance,

on the downstream slope of the dune, the sediment particles have a tendency to roll

down. This must be taken into account in the computation of the bed load transport.

Let T9 represent the effective shear stress at the bed due to local slope of the bed

and it is calculated as:

_ 2 sin(o:)
Tg - u. CT • ~ (7.14)
· sm o:o

where o: is the local slope of the bed, and o: 0 is the angle of repose usually set to

30°. This effective shear stress due to the local slope of the bed is added to the
shear stress at the bed T w, and this total shear stress at the bed is used in the

determination of the transport parameter T to be used in the computation of the

total bed load transport. Note that, the local bed slope only influences the bed load

transport and the suspended sediment transport is independent of it.

Let c represent the suspended sediment concentration defined as volume of sed-

iment material per unit volume of the sediment-laden water. For dilute suspen-

sions ( c ~ 1) the suspended sediment concentration can be assumed to obey the

advection- diffusion equation of the form:

ac ac ac 2 a
ot +uax +vay =(v+vr)\7 c- oy(wsc) (7.15)

The above equation is solved through the ALE formulation and the numerical pro-

cedure outlined in chapters 3 and 5. Note that in the above equation, it is assumed
 134

that the turbulent mixing coefficient vrc of the sediment concentration is same as

the turbulent eddy viscosity vy. Various researchers have previously used mixing

coefficients of the form vrc = f3vr, where f3 is a constant with values ranging from

0.8- 1.3. We arbitrarily set the value of f3 to be 1.0, for lack of better information.

w. is the settling velocity of the sediment particles due to its weight and is computed
as follows:
1(s-l)gD;
D.< lOO~tm,
18 v

100~tm :<:; Ds < 1000~tm, W8

__
10 ~
s
[
D 1+
0.01 (s- 1)gm]o.s _
v 2
1
(7.16)

Ds > 1000~tm, w, = 1.1 [(s- 1)gD.]

0 5
·

D. is the effective particle size of the suspended sediment computed as follows: -

D. =1+0.011(<7.-1)(T-25), (7.17)
Dso

where <78 is the standard deviation of the the particle size distribution. The fall

velocity calculated from the above expression is corrected for the local sediment

concentration effect as follows:

Ws,m = (1 - c) 4 W 8• (7.18)

The boundary condition at the free surface is the no flux condition given by:

ac
vr Y=w.c at y=h(x,t) (7.19)
0

At the sediment bed we impose a reference concentration Ca based on the local shear

stress at a small reference level above the bed a.

c = Ca at y =a+ b(x) (7.20)

 135

The reference level a and the reference concentration Ca are computed as follows:

a= max (k" 0.01h), (7.21)

and

D yts
= 0.015 ~O DD. 3 .
Ca
. (7.22)

Thus, the reference level is set to either the equivalent sand roughness height ks or
1 %of the flow depth, which ever is more. Note that, in the k- E turbulence model,

the wall functions approach is used for the specification of the boundary conditions

at the bottom bed. The first numerical grid point is placed at a small distance away

from the bed such that it lies in the fully turbulent log layer. With the inclusion of
the suspended sediment model, there is an additional restriction on the placement

of the first numerical node near the wall. The first grid point is placed at a small

distance away from the wall so that the conditions set by the k - E model and the

suspended sediment model are both satisfied. Similar to the flow variables, periodic

boundary conditions are imposed in the streamwise direction. The initial condition

for the suspended sediment profile is imposed in the following manner:

_ (a(h-y))z
y < 0.5h, C-Ca y(h-a)
(7.23)
y 2: 0.5h,
_ (ay(h-a)
C-Ca
(h- Y)) z e-4Z(yfh-o.s)

Once the suspended sediment concentration is computed, the total suspended

sediment load is computed as follows:

y=h(x,t)
Ss(x, t) = J.
y=b(x)
u(x, y, t)c(x, y, t) dy (7.24)
 136

From the model outlined in the previous few paragraphs, it is apparent that the
sediment transport model is very complex with several empirical parameters in it.

van Rijn (1984a,b,c) has fine tuned his constants to obtain good agreement with a

wide variety of plume and field data. Several other assumptions are made in the

flow and sediment transport model. With the presence of sediments in suspension,

the fluid flow is a two-phase flow. However, we model the suspended sediment

concentration c as a passive scalar quantity which follows the fluid motion without

it in turn influencing the fluid flow. With all these assumptions it is obvious that

the sediment transport model is only of an approximate nature. But, unfortunately

no good rigorous model exists that can model the sediment transport accurately.

7.4 Results and Discussion

Lyn (1993) obtained experimental measurements of the mean flow and turbulence

characteristics for flow over two-dimensional periodic array of artificial bedforms as

shown in Fig.7.3. The Type I bedform is a periodic array of 45° triangular elements

with height 1.2cm and wave length 15cm. The Type II bedform differs from the

Type I bedform in the upstream geometry. The Type II bedform comprises of

triangular elements with gentle upstream slope and sharp downstream slope and

resembles more closely the naturally occurring bedforms. The above described two

types of bedforms are studied to understand the influence of the bedform shape on

the mean flow and turbulence characteristics.

In the first part of our work, we study the mean flow and turbulence charac-

teristics invoking the rigid-lid approximation. The unknown air-water interface is

replaced with a frictionless impermeable surface and the surface pressure is not

fixed and is computed as part of the solution. The computational mesh points are
 137

(a) g

y
H

(b) g

Parameter case 1 case 2

Bedform Type I Type II
Flow depth, H (em) 6.1 6.1
Wave length, L( em) 15 15
Bedform height, D( em) 1.2 1.2
Energy slope, (} 0.00525 0.00145

Figure 7.3: Bedform shapes: (a) Type I bedform with 45° upstream and
downstream faces; (b) Type II bedform with 45° downstream slope and
a gentle upstream slope.
 138

thus fixed in space and the ALE formulation reduces to pure Eulerian formulation.
The rigid-lid approximation is thought to be a good approximation for low Froude

number flows. In the second part of our work, the air-water interface is allowed to
deform and its equilibrium position is determined as part of the solution. The influ-

ence of the free surface on the quantities of interest such as shear stress at the bed
is examined. For flow over Type II bedform the sediment load distribution along

the channel length is also determined. Of special importance is the suspended load

distribution in the flow separation region as it is thought that the high levels of tur-

bulence associated with flow separation enhance the ability of the flow to transport
sediment particles in suspension.

7 .4.1 Mean flow and turbulence characteristics with the free surface
assumed to be rigid frictionless lid

The mean flow and turbulence characteristics with the free surface assumed as a
frictionless rigid lid are presented in this section.

The steady state mean flow velocities u and v at various cross-sections for flow
over Type I bedform with rigid-lid approximation are shown in Fig.7.4(a,b). The

Nikuradse equivalent sand roughness height k. was taken to be equal to the sand

grain diameter ( = 0.25mm) in this simulation. At a short distance downstream

of the dune crest, the shear layer separates, resulting in a region of recirculation,

which can be easily identified by negative u velocities. The separated shear layer

re-attaches to the bed some point downstream. The reattachment length Lr is

found to be 3.914cm = 3.26D, where D is the height of the dune. Lyn (1993)

measured velocity profiles at only four cross-sections and couldn't determine the
exact reattachment length. From the point of reattachment, an internal boundary
 139

(a)

-2 0 2 4 6 8 10 12 14 16 18
X
(b)

-2 0 2 4 6 8 10 12 14 16 18
X
(c)
10
----
Au 5
-;;r

0
10-1 10°
ylh(x)

Figure 7.4: Mean flow characteristics for flow over Type I bedform with
rigid-lid approximation: (a) horizontal velocity u (in em/ s); (b) vertical
velocity v (in cm/s); (c) Velocity defect profiles.
 140

layer begins to grow and the velocity profiles begin to relax towards the fiat bed

profiles. Following Lyn (1993), the velocity defect D.u = Umax- u, normalized with

the total friction velocity u~ = v'gH sin(} is plotted against y I H, where y is the

vertical distance from the mean bed level (Fig.7.4c). D.ulu~ collapse very nicely

onto a single curve for y I H > 0.4, consistent with the experimental findings of Lyn

( 1993). Thus, the flow is disturbed due to the presence of triangular bedform and

this disturbance persists up toy I H < 0.4 in the channeL The fluid above y I H > 0.4
is relatively undisturbed and the bedform can be thought to act as macro-roughness

element.

The most important feature of the velocity profiles shown in Fig. 7.4( a,b) is the

large vertical and horizontal velocities at the dune crest. To obtain a better insight

into the flow field, the area contour plots for u, v and dynamic pressure p= p-

g( h - y )cos(} are shown in Fig. 7.5. Flow near the dune crest is similar to flow
around a sharp corner and is associated with large pressure gradients and centrifugal

accelerations due to the streamline curvature. This local acceleration results in

large horizontal and vertical velocities at the crest of the dune. · Lyn(1993) also

reported the significant positive vertical velocity at the dune crest. The horizontal

and vertical velocities just upstream and downstream of the dune crest are shown

in Fig. 7.6( a,b ). The local acceleration is not limited to the dune crest alone and

is felt in a small region around the tip of the dune. Due to the local acceleration

associated with the streamline bending, the shear stress at the dune crest is very

high (Fig.7.6c). If the bed were erodible, the tip of the dune would be washed

away immediately and the resulting bedform would be smooth without any sharp

corners. Since we used the rigid-lid approximation, looking at the surface pressures

on the rigid-lid can give us a crude estimate of what the surface deformation is likely
 l4o1

(a)
.....
.....
.......".2T..,
- - ~- ,. - - - ~ . .....
''"'-r~\~; .. ~-·.·~~~~"':~:~~-~-~~~,.:·: :~£=:~ -~-~-,~~ -~ ......
21!1.08

111.-'CI
II.ID

...
10.80

'·""
-L ...
-<.2·
(b)

(c)
111.1
110..6
IOU

....
••u
.....
.....
.....
lll.J.
... 7
-aLe
-103.0
-1'74..3
-1<1.7
-117.0

Figure 7.5: Mean flow characteristics for flow over Type I bedform with
rigid-lid approximation: (a) horizontal velocity u (in cmjs); (b) vertical
velocity v (in cm/s); (c) dynamic pressure p/p = pjp- g(h- y)cosO (in
cm 2 /s 2 ).
 142

(a) (b)
8 8
-x=7.39~:
- X=7.499i
6 6 · -x=7.603(

>-4 >-4

2 2
--------· .-.---c.·-

0 0
0 20 40 60 -20 0 20 40
u u
(c) (d)
~~----~----~----~

15

10
:r,._ a.
p 5

-5L-----~----~----~ ~o~----~----~--~
0 5 10 15 0 5 10 15
X X

Figure 7.6: Turbulent flow over Type I bedform with rigid-lid approxi-
mation: (a) horizontal velocities in the vicinity of dune crest (in cm/s);
(b) vertical velocities in the vicinity of dune crest (in cmjs); (c) shear
stress (Tw/P) at the bed in (cm 2 /s 2 ); (d) surface pressure on the air-liquid
interface (in cm 2 j s 2 ).
 143

to be. The maximum deviation of the surface pressure form the mean pressure is
about 50cm 2 /s 2 as shown in Fig.7.6(d). If the free surface were allowed to deform

the hydrostatic pressure would compensate for this. Thus the maximum surface

deflection would be approximately .SO/ gcosB ~ 0.05cm. The surface deflection is

not very significant and the rigid-lid approximation is probably quite good in this

case. A better estimate of the surface deflection can be obtained if we let the free

surface to deform and these results are discussed in the latter part of this section.

The original motivation for choosing to simulate the experimental cases of Lyn

(1993) is to make quantitative comparisons. Lyn (1993) reports the flow rate for

the above described case as 154.90cm 2 / s. However, at steady state, the flow rate

predicted by the numerical simulation is 215cm 2 / s. We believe the discrepancy-is

due to the fact that, we can not obtain perfect periodic conditions in experiments

as we do in numerical simulations. The equilibrium is achieved when the gravity

force is exactly balanced by the fluid friction at the bed. The friction at the bed

retards the flow, and this effect diffuses slowly upward based on the molecular and

turbulent viscosity. Sufficiently large number of bedforms are required, so that

the momentum diffusion over the entire fluid layer can take place, and the fluid

flow reach equilibrium (steady state). This is in many respects, similar to the

spatial and temporal stability analyses explained in the previous chapter on thin-

film flows. What we are investigating numerically is the temporal stability, wherein

the domain is periodic in space and evolving in time. Whereas, what is being
determined experimentally by Lyn (1993) and other experimenters is the spatial

stability, wherein the disturbance is evolving in the downstream direction. In their

earlier studies (Nelson & Smith 1989; Nelson, McLean & Wolfe 1993) they use

only five periodic bedforms, and find that the flow starts reaching only around fifth
 144

bedform. In their subsequent study, they increase the number of periodic bedforms

to 20 (McLean, Nelson & Wolfe 1994). Due to this discrepancy in flow rates, our
velocities are higher than Lyn (1993) and it is not possible to make quantitative

comparison. Qualitatively however, the numerical predictions compare very well

with the experimental predictions of Lyn (1993).

With the flow rate of 215cm 2 / s and a mean flow depth of 6.lcm, the Froude

number of this flow Fr = Q/ ,f9H3 turns out to be 0.45. Note that, this is no longer

a small Froude number flow, and we should expect to observe significant free surface

deformation.

The turbulence characteristics, namely, the turbulence kinetic energy k, the tur-

bulence eddy viscosity VT and the Reynolds shear stress -u'v', at various cross-

sections are shown in Fig.7.7(a,b,c). Large values of k, VT and -u'v' are found

near the dune crest which are a direct result of the local acceleration associated

with sharp corner. Separation turbulence dominates the wall generated turbulence

which has also been observed experimentally by Nelson eta! (1993). The maximum

k occurs at about the center of the wake. Downstream of the reattachment point
the profiles relax towards the flat bed profiles. However, due to the relatively short

wavelength, the fluid quickly runs into another dune and the flow separation results

without fluid ever attaining flat bed turbulence characteristics. The turbulence eddy
viscosity VT has two maxima for some distance downstream of the dune crest. The

top one is due to the local acceleration at the dune tip and the bottom one is due

to the free stream turbulence associated with flow separation. The Reynolds shear

stress -u'v' has a similar qualitative behavior as k and VT. To obtain a better view
of the turbulence field, the area contour plots of turbulence parameters k, E and VT
 145

(a)

-2 0 2 4 6 8 10 12 14 16 18
X
(b)

-2 0 2 4 6 8 10 12 14 16 18
X
(c)

-2 0 2 4 6 8 10 12 14 16 18
X

Figure 7. 7: Turbulence characteristics for flow over Type I bedform with

rigid-lid approximation: (a) turbulence kinetic energy k (in cm 2 / s 2 ); (b)
eddy viscosity VT (cm 2 /s); (c) Reynolds shear stress -u 1v' (in cm 2 /s 2 ).
 146

are shown in Fig.7.8. To particularly note are the large eddy viscosities upstream
of the dune crest which is due to the local acceleration at the dune crest.

Flow over Type II bedform with energy slope 0.00145 is simulated next. The

equilibrium horizontal and vertical velocities at different cross-sections are shown

in Fig.7.9(a,b). The separation region can be identified by the negative u and v

velocities. The reattachment length Lr was found to be 2.36cm in this case. The

considerably smaller separation zone in this case compared to the Type I bedform is

due to the different upstream geometry. An internal boundary layer grows from the

point of reattachment and the flow profiles relax towards the flat bed profiles. Due

to the shorter reattachment length, the flow field has a longer redevelopment length

and shows greater relaxation towards the flat bed profiles than in the case of Typ~ I

bedform. Just as in the previous case, the flow undergoes local acceleration resulting

in very high vertical and horizontal velocities at the dune crest. The velocity defect

Ll.u = Umax - u, the normalized with the total shear ·velocity u"[ = VgH sinO at

four different cro~s-sections is shown in Fig.7.9(c). Except for the velocity defect

profile at the dune crest xl L = 0.50, at all other three cross-sections, the Ll.u profiles

collapse very nicely onto the same curve. The velocity defect profile at the dune

crest differs from the rest due to the local acceleration associated with stream line

bending.

Similar to the previous case, the equilibrium flow rate computed numerically

is 280cm 2 Is as opposed to the experimental flow rate of 164.90cm 2 Is reported by

Lyn (1993). Based on a flow rate of 280cm 2 Is and mean flow depth of 6.lcm, the

Froude number for this flow is Fr = 0.59. Thus in this case, the Froude number

is more than in the previous case, and we should observe even greater free surface
 147

(a)
.....
.....
117.3

ta.•
1117.2
150.5
133.3
117.1

........
100.4

....
au

ILl

(b)

(c)
....
L-
L-
UIN

,_
UN

1.1111

,....
U6<

.........
UIOII
,
.......
~

Figure 7.8: Turbulence flow characteristics for flow over Type I bedform
with rigid-lid approximation: (a) turbulence kinetic energy k (in cm 2 / s 2 );
(b) rate of dissipation of turbulence kinetic energy E (in cm 2 /s 3 ); (c) eddy
viscosity vr (in cm 2 /s).
 148

{a)

-2 0 2 4 6 8 10 12 14 16 18
X
{b)
8
f
' ' I I

0 15.0
6
>-4
2
0f- I I I
'
I I~
'
.J,
~( \i
-2 0 2 4 6 8 10 12 14 16 18
X
{c)
10
\ /
·~

Au 5 - x/L=O.O
~ -- xiL=0.25
- - x/L = 0.50
··· x/L = 0.75
0
10-1 10°
y/h(x)

Figure 7.9: Mean flow characteristics for flow over Type II bedform with
rigid-lid approximation. (a) horizontal velocity u (in cm/s); (b) vertical
velocity v (in cm/s); (c) velocity defect profiles.
 149

deformation. We will simulate this case without making the rigid-lid approximation
in the next subsection to determine the free surface deformation.

The turbulence kinetic energy k, eddy viscosity vy and the Reynolds shear stress

-u'v' at various cross-sections are shown in Fig.7.10(a,b,c). Very large values are

found at the dune crest due to local acceleration. Wake turbulence associated with
flow separation dominates the wall generated turbulence. Rest of the qualitative

details are similar to that observed in the previous case.

The shear stress at the bed Tw(x) is shown in Fig.7.11(a). Very large values are

observed at dune crest due to local acceleration. It would be interesting to compute

the sediment load distribution in the channel for this flow field. The turbulent flow

field picks up the sediment particles and transports it downstream in the form of

bed load and suspended sediment load. As explained in the previous section on

the mathematical formulation, the bed load is computed using an empirical formula

proposed by van Rijn(l984a) which is a function of sediment particle parameters

such as density, diameter and the shear stress at the bed. As would be expected, the

bed load transport(Sb(x)) is very high in regions with large shear stress (Fig.7.11b)

such as the dune crest. The suspended sediment concentration is computed using

an advection-diffusion transport equation for the suspended sediment concentration

(Eq.7.15). Integrating in they direction the suspended load (S,(x)) as a function x

is computed as follows:

h(x)
s.(x) = l
b(x)
ucdy (7.25)

The suspended load (S,(x)) distribution along the channel is shown in Fig.7.ll(b).
The suspended sediment distribution is dependent on both the flow field and the

shear stress at the bed and is thus only weakly dependent on the shear stress.
 150

(a)
I I I I
8 .
0 50.0
6 1-
>-4 -

2
0 I~ I~
I I
I)
I
I ) I)
I
I
I
L) L~ I~ I~
-2 0 2 4 6 8 10 12 14 16 18
X
(b)
I I I I I I

8 0 1.0
6
>-4
2
0 1- I
I
I
I
I
I
I 1--
I
I__.)
I
J
[
l ) L ) I
I I
.

-2 0 2 4 6 8 10 12 14 16 18
X
(c)

-2 0 2 4 6 8 10 12 14 16 18
X

Figure 7.10: Turbulence characteristics for flow over Type II bedform

with rigid-lid approximation: (a) turbulence kinetic energy k (in cm 2 ls 2 );
(b) eddy viscosity VT (in cm 2 Is); (c) Reynolds shear stress -u'v' (in cm 2 I s 2 ).
 151

(a)
15 I

10 -

~
p 5-

I I
-5
0 5 10 15
x(cm)
(b)
0_1 -I I

0.08- -
0.06-
Sb,Ss
0.04

0.02
.. I
0
0 5 10 15
x(cm)

Figure 7.11: Flow over Type II bedform with rigid-lid approximation:

(a) shear stress at the bed (in cm 2I s 2 ); and (b) sediment load in the form
of bed load (continuous line, in cm 2 Is) and suspended load ( dott~d line,
in cm 2 ls).
 152

Downstream of the dune crest in the flow separation region, the bed load transport

is negligible, where as there is significant amount of suspended load. Thus models

which neglect the suspended load and take into account only the bed load would be
inaccurate in the separation region. The significant suspended sediment concentra-

tion is due to the wake turbulence which maintains the sediment particles in suspen-

sion. In our suspended sediment modeling, the sediment concentration is modeled

as passive scalar and hence doesn't influence the flow field. Near the bed, a small
layer of sediment particles of the order of few particle diameters thickness moves

as bed load and is known to enhance the effective roughness of the bed. Following
van Rijn(1984c), the effective sand roughness height ks is taken to be equal to 3D 50

where D 50 is the median diameter of the sediment particle. The flow field has been

computed with ks = 0.25mm and the sediment transport calculations correspond to

sediment bed with 0.08mm sediment particle size. The sediment particle size is very

small and consequently, we find significant levels of suspended sediment transport.

However, the suspended sediment load is known to decrease with an increase in the

particle size. Thus, whether suspended sediment load also needs to be considered

along with the bed load is dependent on the sediment material.

7.4.2 Mean flow and turbulence characteristics without the rigid-lid as-
sumption

In this section we present the results for the simulations in which the air-liquid

interface is allowed to deform. The location of the free surface is computed along

with the fluid flow and turbulence field using the ALE formulation.

The horizontal velocity u at various cross-sections for the flow over Type I bed-

form are shown in Fig.7.12. The horizontal velocity profiles look very similar to
 153

(a)

0 2 4 6 8 10 12 14 16 18
X
(b)

0 2 4 6 8 10 12 14 16 18
X
(c)

····· - xll=O.O
··············· ....
------ --- --- ·· ..
-- xll= 0.25
6.u 5
;r - - xll= 0.50
..... xll=0.75

ylh(x)

Figure 7.12: Flow over Type I bedform with the free surface allowed to
deform: (a) horizontal velocity u (in em/ s), (b) vertical velocity v (in
cm/s), and (c) the velocity defect profiles.
 154

when the free surface was assumed to be rigid (Fig.7.4). The most important fea-

tures are the local acceleration at the tip of the dune, separation of the boundary

layer downstream of the crest, and the reattachment and subsequent growth of the

internal boundary. The recirculation length is found to be 3.25cm as opposed to a

3.9cm recirculation length observed when the free surface was assumed to be a rigid

lid. The vertical velocity profiles are shown in Fig.7.12(b). The vertical velocity

profiles especially near the free surface are markedly different from that obtained
previously with the rigid-lid assumption (Fig. 7.4b ). The velocity defect profiles are

plotted in Fig.7.12(c), and fall on a single curve for yjh(x) > 0.4 just as in the
previous simulations with rigid-lid assumption (Fig. 7.4c ).

The turbulence kinetic energy k, the eddy viscosity vr and the Reynolds shear
stress -u'v' are shown in Fig. 7.13. The qualitative behavior is seen to be very similar

to that observed previously with the rigid-lid assumption. The turbulence associated

with separation wake completely dominates the wall generated turbulence.

The equilibrium free surface profile is shown in Fig.7.14(a) and compared with

the rigid-lid assumption. The maximum free surface deflection is approximately

0.3cm. The minimum free surface height occurs near the point of reattachment and

the maximum free surface height occurs well upstream of the dune crest. The most

important quantity of interest is the shear stress at the bed, and a comparison is

made between that obtained with and without the rigid-lid approximation. The

shear stress is found to be lower in the case where the free surface was allowed
to deform, compared to that obtained with the rigid-lid assumption. A very big

difference in the shear stress is found near the tip of the dune, and this can have

profound influence on the sediment transport and bedform formation modeling.

 155

(a)

-2 0 2 4 6 8 10 12 14 16 18
X
(b)

-2 0 2 4 6 8 10 12 14 16 18
X
(c)
150.0

-2 0 2 4 6 8 10 12 14 16 18
X

Figure 7.13: Flow over Type I bedform with the free surface allowed to
deform: (a) turbulence kinetic energy k (in cm 2 /s 2 ), (b) eddy viscosity VT
(in cm 2 /s), and (c) Reynolds stress -u'v' (in cm 2 fs 2 ).
 156

(a)
6.5r---------,---------.------

I
I

Figure 7.14: Comparison between the solutions obtained with and with-
out the rigid-lid assumption: (a) The continuous line corresponds to the
equilibrium free surface profile (in em), the broken line corresponds to
the bed profile (shifted upwards), and the dotted line is the rigid-lid;
(b) the shear stress (in cm 2 / s 2 )shown in continuous line is obtained by
allowing the free surface to deform and the broken line is obtained by
making the rigid-lid approximation.
 157

Flow over Type II bedform (Case 2) with energy slope 0.00145 is simulated next.
The rigid-lid approximation is not made, and the free surface location is determined

through the ALE formulation. Unlike the previous case, the fluid flow does not go
towards steady state. Instead waves are found to occur in this case on the air-liquid

interface, leading to a time-periodic free surface height as shown in Fig. 7.15. Note
that, oscillatory behaviour at x = 0 is different form that at x = L/2. At x = 0,
the free surface height oscillates between a minimum of 5.9cm and a maximum of
7.1cm. Where as, at x = L/2, the free surface height oscillates between 5.3cm and
6.0cm. This tells us that this is not a traveling wave. It is probably a standing
wave. This can be better understood through the spatia-temporal evolution of the
free surface profile shown in Fig. 7.16. The wave is found to propagate both upstream
and downstream, alternately. This is unlike that seen in thin-film flows (chapter 6),

where the surface waves were found to travel downstream. This is due to the fact

that the fluid flow is in a subcritical state (Fr < 1), and the disturbance can travel
both upstream and downstream. The time history of the shear stress at the bed

is shown in Fig. 7.17. The shear stress at the bed is also seen to be time-periodic.
This type of periodic behaviour can have profound influence on the sediment load

and consequently, the bedform formation. Also, the deformation of the free surface
is quite large, and the simulations made with the rigid-lid approximation would be

expected to perform poorly.

As mentioned briefly in chapter 2, gravity surface waves are found to exist on
open-channel flows with Froude numbers greater than a cut-off value usually in the
range 1 - 2. However, these estimates are for flow over a smooth wall. In our case,
the bed is filled with bedforms, which probably, reduce the cut-off value for initiation
 158

(a)

5.5L__L__L__L__L__L__L__L__.L___.L______j
0 5 10 15 20 25 30 35 40 45 50 -
time( sec)
(b)

5 10 15 20 25 30 35 40 45 50
time(sec)

Figure 7.15: Flow over Type II bedforms. Time history of the free surface
height at (a) x = 0 and (b) x = 7.5cm (the crest of the bedform).
 159

40 I I

-- ~
35

30
t:-::: ~
--.

~
(.)
Q)
~
1/)
._,
Q)
E
:;::;
-
f.---:: ~
::::::::--==::: ~ ~
25

20
~

\_ I
15
0 5 10 15
x(cm)

Figure 7.16: Spatia-temporal evolution of the free surface profiles in the

case of flow over Type II bedform without making the rigid-lid approxi-
mation.
 160

(a)
4

:c... 2
p

0
0 5 10 15 20 25 30 35 40 45 50
time( sec)
(b)
20

5 10 15 20 25 30 35 40 45 50
time( sec)

Figure 7.17: Flow over Type II bedform without making the rigid-lid
approximation. Time history of the shear stress (rw/P) at the bed at (a)
x = 0; and (b) x = 7.5cm (the crest of the bedform).
 161

of gravity waves on the free surface. Thus, even with a Froude number flow of only
Fr = O.S9, we are observing surface waves on the free surface.

vVith the same bedform geometry as the previous case, that is Type II bedform,
but with energy slope 0.00400 is simulated next. Due to the greater energy slope,

the flow rate increases resulting in higher Froude number flow, which in this case

is close to Fr = 1.0. Steady traveling waves on the free surface are found to exist

in this case (Fig./.18\. The waves are traveling downstream and the flow is in a
supercritical state.

7.5 Concluding remarks

A numerical method has been developed to study the turbulent fluid flow over

stream-wise periodic bedforms. The procedure can handle the free boundary in a

straight forward manner through the ALE procedure, and gives us the capability

to simulate high Froude number flows. Numerical procedures with the rigid-lid

approximation developed previously (Mendoza & Shen 1990) are strictly valid only

for very small Froude number flows.

Study of the turbulent flow over artificial bedforms has revealed certain inter-

esting features not previously reported. Due to the local acceleration associated

with streamline bending, very large velocities and shear stresses exist at the tip of

the dune. If the bed were erodible, the tip of the bed would be the first area that

would be washed out, and the resulting bedform would be smooth. The influence
of the bedform is restricted to the bottom half of the fluid layer, and velocity defect

profiles for the top half of the fluid layer are found to be self similar with respect to

the total friction velocity u~ = v'g H cos(). The separation wake turbulence is found

to dominate the wall generated turbulence completely, and the maximum eddy vis-
 162

7.5 I

7
- -
6.5"" -
-
......
.-
(J
Ql
en
"-'
Ql
61-
E
:;:

1-
5.5 1-

=-
-
5

-
4.5'-----------'L--------l----------l
I I
-
0 5 10 15
x(cm)

Figure 7.18: Spatio-temporal evolution of the free surface profiles in the

case of flow over Type II bedform with energy slope 0.004, and without
making the rigid-lid approximation.
 163

cosities and Reynolds stresses occur at a distance, approximately equal to the dune

height, away from the bed.

In the case of flow over Type I bedform with energy slope 0.00525. a maximum

free surface deflection of 0.3cm was found. However, the biggest difference between

the two simulations, namely. with and without making the rigid-lid assumption, is in

the shear stress values. The simulations without making the rigid-lid approximation

predict a much smaller \'alue of shear stress at the peak.

In the case of Type II bedform with energy slope 0.00145, no steady state solution

was found. Instead standing waves were found to occur on the free surface and

the waves propagate both upstream and downstream since t'he fluid flow is in a

subcritical state. This periodic free surface waves, result in periodic shear stresses

at the bed and this can have profound influence on the sediment load characteristics

and bedform formation.

 164

Chapter 8

Hydraulic Jump

8.1 Introduction

Hydraulic jump is one of the most interesting phenomena occurring in open channel

flows and is analogous to shock wave in compressible flows. The free surface rises

sharply across the hydraulic jump and the flow character changes across the jump,

going from a supercritical state to a subcritical state. There is intense turbulent

mixing in the hydraulic jump region and a large amount of energy dissipation occurs

in the jump region. Because of its energy-dissipating capacity the hydraulic jump

is widely used as an energy dissipator below dams, spillways and other outlet works

(Chow 1959; Rajaratnam 1967). The hydraulic jump has many other practical

applications as described by Chow(1959). The hydraulic jump may occur when

there is supercritical flow in a channel having an obstruction or a .rapid change

in cross- sectional area (Shames 1992). Usually a stationary jump is referred to

as a hydraulic jump and a jump propagating into quiescent fluid is referred to as

hydraulic bore. A hydraulic bore is very similar to a hydraulic jump and can be

visualized as a hydraulic jump in a moving frame of reference. The important

difference between a bore and a jump is that the bottom boundary layer is weaker

in the bore compared to hydraulic jump (Madsen and Svendsen 1983).

Rajaratnam (1967) defines the classical jump as a jump formed in a smooth,

wide, and horizontal rectangular channel. A typical hydraulic jump is shown schemat-
ically in Fig.8.l. h 1 and U1 are respectively, the depth and mean horizontal velocity
 165

L·
J
~
t
Lr
1:

air u2
h2
water
Ut
hl

Figure 8.1: Schematic sketch of a hydraulic jump

 166

of the supercritical flow just upstream of the jump. h2 and u2 are respectively, the
flow depth and mean horizontal velocity of the subcritical flow just downstream of

the jump. A hydraulic jump is characterized by the formation of the surface roller
and air-entrainment into the roller. Lr is the length of the surface roller and L j is

the length of the jump. The most important non-dimensional number in the open
channel flows is the Froude number defined as:

u
Fr = y'gl! (8.1)

where U is the mean velocity and h is the flow depth. There are three distinct
regimes in an open channel flow:

• Fr > 1, supercritical flow;

• Fr = 1, critical flow; and

• Fr < 1, subcritical flow

Froude number is the ratio of the mean velocity to the gravity wave celerity. Thus in

a supercritical flow, the wave celerity is less than the mean velocity and the distur-

bance can only travel downstream. In a subcritical flow, however, the disturbance
can travel both upstream and downstream. The Froude number is thus analogous

to Mach number in compressible flows.

The quantities of primary interest in hydraulic jump are the downstream quan-

tities such as h 2 and U2 , the profile of the jump and the energy loss across the jump.

Knowing upstream quantities U1 and hr, the downstream quantities h 2 and U2 can
be determined analytically if the following assumptions are made:

• uniform velocity profile;

 167

• pressure is hydrostatic:

• bottom friction is negligible and can be neglected: and

• surface tension can be neglected

Using the integral conservation of mass and conservation of momentum, a relation

relating h1 and h 2 , can be derived and is also called as the Belanger momentum
equation (Rajaratnam 1967).

~: = ~ ( j1 +SF( - 1) (S.2)

Even though there is conservation of momentum, there is a loss of energy across the

hydraulic jump. Define the specific energy as E = U 2 /2 + gh. The loss of specific

energy across the jump, EL = E2 - E 1 , can be expressed as (Rajaratnam 1967):

EL 1 (J1+SF(-3r
(S.3)
£ 1 = S(2 + F() ( ,j1 +SF(- 1)

The specific energy loss across the jump depends on F 1 , with a larger F 1 resulting

in a greater energy loss across the jump.

The nature of the hydraulic jump is strongly dependent on the supercritical

Froude number F 1 • Chow (1959) classified the hydraulic jump into the following

categories based on F 1 :

• F1 = 1- 1.7, undular jump;

• F 1 = 1. 7- 2.5, weak jump;

• F1 = 2.5 - 4.0, oscillating jump;

• F 1 = 4.5 - 9.0, steady jump

 168

• F1 ~ 9.0, strong jump

The above listed various types of jumps are shown schematically in Fig.8.2. The
various jumps as listed above are not completely distinct and overlap to certain

extent depending on the local conditions.

8.2 Literature Review

Hydraulic jump was first described by Leonardo da Vinci (Rouse & Ince 1957). The

first systematic analysis of the hydraulic jump was undertaken by Bidone (1820) and

Belanger (1828). The hydraulic jump has been extensively studied since then and

are listed by Chow (1959). A good review of all the work done on hydraulic jump

prior to 1967 can be found in Rajaratnam(1967). Most of the earliest investigations

have been experimental with the aim of finding the external characteristics of the

jump such as length of the jump and the roller, surface profile of the jump and the

sequent depth ratio h2/ h1 of the jump.

The earliest numerical studies neglected the variations in the physical quantities

along the vertical (and cross-sectional) direction and assumed the pressure to be
hydrostatic. This results in one dimensional St. Venant equations and the shear

stress at the bed was modeled using uniform flow relations such as Manning's or
Chezy's equation (Lai 1986). Basco (1983) has developed one dimensional open

channel flow equations which take into account the non-hydrostatic effects of pres-

sure and these are called the Boussinesq equations. Most of the numerical study of
open channel flows has been to solve these one dimensional St. Venant equations

or Boussinesq equations (Lai 1986; Basco 1983). The method of characteristics, the

Finite Difference Method and the Finite Element Method have been used in the
solution of these one-dimensional equations. The hydraulic jump is modeled as a
 169

-
'/ffi./J//////1.7/I//I/////7///7////7//7////ff)///)/
F, •1-t7 Undulcr jump

~=:; s~- - -
"""''l!!"".....

-- -
'l//7/////J/)///////1'///7////ff//////////1////d/
F1 • t1- 2.5 Weak JUmP

,_$}'/?")
--............. _
-::..
--=---//'/ _..,
/////U/7////7/////H/////////////////1/ffd///h
F1 • 4.5-9.0 Steady jump

W/////7//7/7/H)//,,,Y////H///////ff///////////h
'• >9.0 S!rMq jump

Figure 8.2: Various types of hydraulic jump (adapted from Chow (1959))
 170

discontinuity and the surrounding fluid flow is computed using the one-dimensional
equations. These one-dimensional equations however, do not describe the internal

characteristics of the hydraulic jump and depend for their accuracy on the empirical

relations employed to model the bed friction. A some what more accurate one-

dimensional model based on the depth-averaged flow and k- c turbulence equations

and hydrostatic pressure has been developed by Madsen and Svendsen ( 1983) for

the study of hydraulic jump.

For a description of the internal flow field and the turbulence characteristics of the

hydraulic jump, the two-dimensional or three-dimensional Navier- Stokes equations

without any depth averaging will need to be solved. Two- dimensional modeling of

hydraulic jump is however computationally very intensive and only recently with the

increased computing power it is becoming feasible. Rahman ( 1991) solved the steady

two- dimensional flow equations to simulate a hydraulic jump. Finite- differences

in curvilinear boundary-fitted coordinate system was used and the free surface was

determined iteratiydy by minimizing the error in pressure at the free surface. Lemos

(1992a,b) solved the unsteady flow equations and updated the free surface in time

using the Volume-Of-Fluid (VOF) method. Both of them modeled turbulence using

the k - c closure equations.

8.3 Problem Definition

A typical 2-D hydraulic jump flow configuration is shown in Fig.8.3. The flow is
bounded below by a non-erodible plane wall and at the top by a moving free surface

whose height from the bed is taken to be h(x, t). x = 0 is the inlet cross-section and
is designated as station 1. x = Lis the exit cross-section and is designated as station
2. section 1' corresponds to the cross-section just upstream of the jump when the
 171

lj
I·
lr
I· ·I FREE SURFACE
..

(!)
CD :
y :
: h(x,t)

4ep
1------l
--~~~i'77TJ7777/77fTT~At~777-

L--------------------~

Figure 8.3: Typical Geometry for a hydraulic jump

 172

free surface begins to rise abruptly and section 2' corresponds to the cross-section

just downstream of the hydraulic jump. The subscripts 1,1',2 and 2' refer to the
cross-section at which the quantity in question is evaluated.

The governing equations for the fluid flow are the Reynolds Averaged Navier-
Stokes (RANS) equations and the continuity equation.

U,x + V,y = 0 (8.4)

U,t + UU,x + VU,y = -gh,x- !P,x

p
+ (v +liT) 'i7 2 u (8.5)

(8.6)

u and v are the velocities in x and y directions respectively. t is the time and v is

the laminar viscosity. liT is the kinematic eddy viscosity and is determined using

the k - c model described in chapter 4. p is the non-hydrostatic pressure and is

related to the pressure p as follows:

p=pg[h(x)-y]+p (8.7)

The boundary conditions at the wall are the usual no-slip, impermeability con-

ditions:

u = 0 }
at y = 0. (8.8)
v = 0

Since, the turbulence model used is only valid in the fully turbulent outer layer, the
wall function approach is used to impose the boundary conditions at the wall and

these are described in chapter 4.

At the free surface we have the stress-continuity conditions and the kinematic
free surface condition. Assuming the datum atmospheric pressure to be zero, and
 173

neglecting the viscous stresses at the free surface. since in general they are small

compared to the inertial forces. we have the following simplified boundary conditions
at the free surface:

p 0

au 0
an at y = h(:r, t).
av 0
(8.9)
an
oh oh
at + u ax == u

The first condition is imposed in the pressure Poisson equation, and the next two

are imposed in the solution of the momentum equations. The last condition is used

to update the free surface location.

The boundary condition for u at the inlet depends on the state of the flow

entering the fluid domain. If the flow at the inlet is assumed to be undeveloped

then the boundary conditions for velocity at the inlet is given by:

u = ul ) at x = 0. (8.10)
v = 0

On the other hand if the flow at the inlet is assumed to be fully developed then

velocity field satisfying the logarithmic law is to be specified at the inlet.

u = u.lnEyu.)
"' v at x = 0. (8.11)
v = 0

u. = ;;::JP is the shear velocity, Tw is the shear stress at the wall, E is a constant

dependent on the roughness of the wall and "' is the von Karman constant taken to

be equal to 0.4. In a very long fluid column, if the hydraulic jump occurs very close

the inlet, then the fluid flow is most probably in an undeveloped state. Thus, the
 174

fluid flow entering the fluid domain is assumed to be potential flow, and without

turbulence. If the hydraulic jump occurs far downstream from the inlet, there is

sufficient time for the flow to come to equilibrium, and the fluid flow just upstream
of the hydraulic jump is a fully developed turbulent flow.

At the exit a fully developed velocity profile boundary condition is imposed.

au
ax = 0 }
av at X= L. (8.12)
- = 0
ax
The finite element mesh used was 401 X 41 is the case of inlet supercritical Froude

number F1 = 2.0, whereas a numerical mesh of 801 X 11 was used for F 1 = 4.0. In the
case of F1 = 2.0, uniform mesh spacing was used. However, in the case of F 1 = 4-0

the following mesh spacing was used in the x direction:

L 1 L
for 0 <X<-
- - 3
3 200
L 1 L 2L
~x= for - < x < - (8.13)
3 400 3 - - 3
L 1 2L
for -<x<L
3 - -
3 200

This typeof a mesh spacing is used for F 1 = 4.0 to obtain fine mesh spacing in the
jump region. A typical finite element mesh is shown in Fig.8.4. Note that, the mesh

shown is only for part of the physical domain.

8.4 Results and Discussion

Let y 1 and y2 be the flow depths upstream and downstream of the hydraulic jump

respectively and U1 the mean horizontal velocity of flow at the inlet. A hydraulic

jump on a straight 2-D horizontal channel is influenced by the following parameters:

1. Inlet Reynolds number Re1 = U1ydv.

 175

Figure 8.4: Finite Element Mesh for part of the domain

 176

2. Inlet Froude number F 1 = UJ/ .,fffYl.

3. Wall roughness height k 5 •

4. Type of flow at the inlet; whether fully developed (FD) or undeveloped (UD).

The governing equations have been solved in a non-dimensional form with the inlet
flow height Y1 and inlet mean velocity U1 used as reference length and reference
velocity respectively. We simulated the hydraulic jump with four different inlet
conditions which are listed in Table 8.1.

For all our simulations inlet Reynolds number Re 1 = 1.0 x 106 and and the wall
is hydraulically smooth (ks = 0.0).
Mathematically, the jump will occur when y 1 , y 2 and F 1 satisfy the Belanger
momentum equation (Chow 1959):

y2 = J8F{ + 1 - 1
(8.14)
Y1 2

In the numerical model the downstream subcritical flow height y 2 is determined from

the above equation and is specified as an essential boundary condition in the solution
of free surface kinematic equation. In nature however, hydraulic jump occurs when

supercritical flow encounters an obstacle or sudden change of cross-section. It is

being assumed that the properties of the hydraulic jump numerically simulated by

case FI Re1 k. inflow state YI Y2 L

2(UD) 2 1.0 X 10 6 0.0 Undeveloped 1.0 2.37 237.0
2(FD) 2 1.0 X 106 0.0 Fully developed 1.0 2.37 237.0
4(UD) 4 1.0 X 10 6 0.0 Undeveloped 1.0 5.17 517.0
4(FD) 4 1.0 X 106 0.0 Fully developed 1.0 5.17 517.0

Table 8.1: List of simulation parameters

 17i

specifying the downstream subcritical flow height y 2 are nearly the same as the

naturally occurring ones. The nature of the hydraulic jump is mostly dependent

on the upstream supercritical flow conditions and the etfect of the downstream

boundary condition is mostly on the location of the hydraulic jump.

The flow after it enters the domain slows down due to friction at the wall and the

free surface height rises. As shown in Fig.8.3, 1' is the cross-section just upstream of

the jump and 2' is the cross-section just downstream of the jump. The flow heights

and Froude numbers for sections 1' and 2' are listed in Table 8.4. For comparison

with theoretical and experimental results we will need to use supercritical Froude

number Fl' and not inlet Froude number F 1 .

The velocity vectors in the hydraulic jump region for the four different cases

that have been simulated is shown in Fig.8.5. The plot is drawn to scale to draw

attention to the vastly differing length scales in horizontal and vertical directions.

The subsequent figures are not drawn to scale and we caution the reader not to

be misled by what appears to be rapid changes in free surface height. The x and

y coordinates of the four corner points are shown on all the plots so that it will

be clear what part of the fluid domain we are looking at. Due to reasons already

explained in the section on turbulence the first numerical grid point is not placed

on the wall but is placed at a small distance away from the wall in such a way that

it lies in the completely turbulent wall layer. Thus the bottom surface of all the

plots do not correspond to the wall and a small layer of fluid close to the wall is not

represented in our plots.

In Fig.8.6 the variationuj ul, vI ul and and PI pUf for the complete flow domain,
i.e. 0.0 ~ x ~ L, is shown for the case F 1 = 2(UD). According to Chow (1959)

an undular jump is formed for 1.0 < F 1 < 1.7 and a weak jump is formed for 1.7 <
 case XI' g, F2' Yl 1 Y2' E1' E2' &
E, Lr/Y2 L1/Y2 L,ep/Y2
2(UD) 100.2 1.78 0.53 1.07 2.36 2.79 2.71 0.028(0.056) 0.00 15 0.25
2(FD) 103.2 1.82 0.54 1.06 2.38 2.83 2.74 0.031(0.062) 0.25 15 0.25
4(UD) 245.5 3.27 0.37 1.13 5.52 7.27 5.52 0.240(0.296) 1.83 15 0.00
4(FD) 204.5 3.42 0.35 1.10 5.63 7.63 5.63 0.262(0.318) 2.33 15 0.00

Table 8.2: External characteristics of hydraulic jump

,_.
-l
C/)
 179

(100.2. 1.08) rJt1lTfJnmmmTr (135.8. 2.37)

llllllllll~lllllliliililllillilllliilllli
ooo.2. o.oo) (a) (135.8. o.ooJ

(138.8. 2.311)

(l0 ~~i~~••••• rtfiiJ]JIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIi

3

(103..2, 0.00) (b) (138.8. 0.00)

(24S.S. 1.1~)

(2~S.S. 0.00) (c) (323.2. 0.00)

(204.5. 0.00) (d) (282.2. 0.00)

Figure 8.5: Velocity vectors in the hydraulic jump region: (a) F 1 = 2(U D),
(b) F 1 = 2(FD), (c) F 1 = 4(UD), (d) F1 = 4(FD).
 180

L003

(237.2, 2.37)
....
·""'
,.,.,
....
·''"
·'"
.475
.400

....
·""
. 17-'

·""'
.023
-.OIU
( 0.0, 0.00) (237.2, 0.00)
(a)
.....
.....
...,.
..,..
.3148

.2201
.lMl
.!...
•1<102
.00?3
.ot.a
_,...
-.0'716
-.11+4
-.115'1'4
( 0.0, 0.00) (237 .2, 0.00)
(b)
.!...
......
.....
·"""
.....
.....
---.-
-.0153

__ .,.,.
-......
-.1188
-.18113
-.18011
-.1801
( 0.0, 0.00) (237.2, 0.00)
(c)

Figure 8.6: Mean flow field for F 1 = 2(UD) (a)u/U1 ; (b)v/U1 ; and (c)p/pU{.
 181

F1 < 2.S. In our case F1' = l./8 and the jump appears to exhibit characteristics

of both undular jump and weak jump. The free surface in the jump region is wavy

and the flow downstream of the jump is smooth and steady. The vertical velocity v

and non-hydrostatic pressure correction pare significant only in the jump region.

To get a better understanding of the nature of flow in the jump, u/U1 , vjU1

and fJ/ pUf are plotted for the jump region only, i.e. x 1 , :S: x :S: x 2 ,, in Fig.8. 7 for

F1 = 2(UD). A small recirculation region exists under the jump. Boundary layer
separation in hydraulic jumps has been suspected for long and Leutheusser & Alemu

( 1979) presented experimental evidence confirming the presence of flow separation

under hydraulic jump. Even though the flow appears to be be slowing down and

stagnating near the free surface, there is no surface roller i.e, there is no backwa~d

flow on the free surface. Departure of pressure from the hydrostatic value is quite

large with fJ/ pUf being as high as 0.18.

The variation -::-rpk ,

Pt
~';
PtYl
and u"x
tYt
is shown in Fig.8.8 for F 1 = 2(UD). The

maximum value of turbulent kinetic energy k is only 1.2% of the inlet kinetic energy

and occurs near the separation zone. The dissipation rate of turbulence kinetic

energy t takes very large values near the wall. In fact the boundary condition for t

at the wall was tp = u~ + . In the rest of the flow domain t is not as high
~<yp(l.O-e-•P 126 )
as it is at the wall and hence the linear plot shows an essentially constant field for t.

The region of maximum turbulence kinetic energy and maximum turbulent viscosity

are not the same. Maximum turbulence kinetic energy occurs further downstream of

the maximum turbulence kinetic energy cross-section. This is due to the important

role played by convection which transports k downstream.

The production( G), dissipation( t) and convection of turbulent kinetic energy for

F1 = 2(UD) are plotted in Fig.8.9 for the flow domain obtained after eliminating
 182

·-....
"""
8 32

"'
.1111
.537

(100.2. 0.00) (a) (135.6. 0.00)

·""'
.4008
...,,.
.3145
2720

..
.2201

,,.
.11181

. 1002

·""
-.-
.0143

-.071&
-.11+4
-.115'1'4
(100.2, 0.00) (135.8, 0.00)
(b)
·"""
·-...............
·"""'

-.....
-.0183

- .....
-.am
__ ..,.
-.111111
-.lSICJ
-.11100
-.1108
( 100.2, 0.00) (135.8, 0.00)
(c)

Figure 8.7: Mean flow field for F 1 = 2(UD) (a)u/Ut; (b)v/Ut; and (c)P/ pUf.
 183

.........
·-..............
,OOI!i101

-
.006111

.......
.......
.0021Ut

.oourrr

.......
.0014.06

( 100.2. 0.00)
(a) (135.8, 0.00)
·-,_
·-
.1401

.l.la?
.UUI

·-.........
.100<

,
·-·
....
.....
.....
.0111

(135.8, 0.00)
(b)
.......
·--..............
....,
,

·-·-........
·-... ...
.0Dm14
,
.......,
(1oo.2. o.oo) (c) (135.8, 0.00) '

Figure 8.8: Turbulence characteristics for F 1 2(UD) (a) k; (b) t; and

(c)vT.
 184

a small layer of fluid near the wall ( 4 finite element cells) where production and
dissipation of turbulence kinetic energy is very high. This makes it possible to study

the variation of G and E in the rest of the flow domain. Production and dissipation
are both large in the separated flow region. The difference between production

and dissipation equals convection. Convection of turbulent kinetic energy is quite

significant in some regions of the flow. It is as large as E itself in some areas. The

zero equation turbulence models like Prandtl's mixing length models assume that

that production and dissipation are in balance every where and neglect transport

of turbulent kinetic energy. But in separation flows convection of kinetic energy is

quite significant and justifies our use of higher order turbulence model.

The variation of u/U1 , v/U1 and P/ pU{ when the flow at the inlet is fully devel-

oped and the F 1 = 2 (designated as 2(FD))is shown in Fig.8.10. The jump is less

undular compared to the case when the inflow was undeveloped. This could be due

to higher supercritical Froude number F 1, = 1.82 compared to the previous case.

The biggest difference from the previous case (F1 = 2(UD)) is the presence of a
surface roller.

~·
In Fig.8.11 the variation of turbulence parameters ::u'Ik ,
P 1 1
P 1 Yl
and ·u"x
IYI
is shown
for F 1 = 2(FD). The maximum turbulence kinetic energy k is about 5.0% of the

inlet kinetic energy and occurs in the surface roller region. The higher values of k

and vr compared to the previous case when inflow was undeveloped is solely due

to the presence of surface roller. Surface roller from our simulations appears to be
the most important turbulence generating mechanism and is consistent with the

findings of Rouse, Siao & Nagaratnam (1958).

The production, dissipation and convection of turbulent kinetic energy are shown

in Fig.8.12 for F 1 = 2(FD). As in the previous turbulence budget plot a very small
 185

.......
_..,..,
........
..,..,.
_.....,
.001?111
.001611:1

--
.001:!34
.OOUM

........
JIOO&'Pl
_ _,.,
......
.00111110

(100.2, 0.01)
(a) (135.8, 0.02)
..... ,
.......
_..,.,.
.
. 001174
.OOl?M
.00111115

...,...
.OCH4315

.00101"1
........
.......
--
·"""""
........
(100.2. 0.01) (b) (135.6. 0.02)

--.,.
.......
·--
--.......
.-.

-·-
- -ll
.........
.......
. 0000181
-
-.OOCUcrnl
-.0001-
(100.2, 0.01) (135.6. 0.02)
(c)

Figure 8.9: Turbulence kinetic energy budget for F 1 = 2(UD) (a) produc-
tion ; (b) dissipation; and ( c )convection.
 l86

....
·"''
.....
"'
'"'
.,.
...
=
,.,.
.30C

...
""'
·""
_.,.
-002

(103.2, 0.00) (138.8, 0.00)

(a)
.....
·"""
-"'1
~·?0

.....,
-""'
....,
.1?0<

·'""
.00'10

-"""
......
-.0121

(103.2, 0.00) (136.6, 0.00)

- .....
-.04il!l

(b)
.L213
_,...,
.....
.....
.....
......
.....
-.-.....
-__ ....,
_..,....
- .....
-.102'1
-.l.a14
-.1.01

(103.2, 0.00) (138.8, 0.00)

(c)

Figure 8.10: Mean flow for F 1 = 2(FD) (a)u/U1 ; (b)v/U1 vertical velocity;
and ( c )P/ pU{.
 1~7

.-..
......
""''
. 01124
.017UI
.01674
.013VIiil
.012211

_,,....
.010!10

.OCI'fOO
.QQ080

-""""
"""'
(I 03.2, 0.00)
(a) (138.8, 0.00) ""'
·'""
_,..,.
_, ...
. 1611'
......
...
. 1111

----
_,

.......

--
.OllilA

.01111
.QQQQ

(I 03.2, 0.00) (138.8, 0.00)

(b)
......
.014211

......
.oum

......
.01018

.oc:.te
.QQTll
.001111

·""""'
·""""
......
·""""'
.001. .

(103.2, 0.00) (138.8, 0.00)

-""""
(c)

Figure 8.11: Turbulence characteristics for F 1 2(FD) (a) k; (b) t:; and
(c)vT.
 188

·""""'
.0074111

·"'"""'
·""'""
.Oilfo'ro?
.006131!
.oo4Me
.0038its
.003<24
.002803

·"""""'
.0017111
.001141
.0006?0

(103.2, 0.01) (138.8, 0.02) ·""""'

(a)
........
........
.OOIS614

·""'"''
·""""'
·""""'
·""'"
.001813
.001883
.001SS3
.OOl!Hr.l
000812

·"""""
..,.,..,.
·"""""
·""""
.0041"
.OOS?liO
.........
·"""'""
""""'
.0018oU
.001600

"'""
.OOCMill
.0001111!1

-""""
-.000718
-.ooua
-.ootear
(103.2, 0.01) (138.8, 0.02)
(c)

Figure 8.12: Turbulence kinetic energy budget for F 1 2(FD) (a) pro-
duct ion ; (b) dissipation; and ( c )convection.
 189

layer ( 4 finite element cells) near the wall is not represented in this plot. The

production( G) and dissipation( E) of turbulent kinetic energy are very large in the

surface roller region. Dissipation is only half of production in the surface roller

region and the rest of the produced kinetic energy is convected away by the mean

flow and once again justifies the use of k - E closure model.

For the cases F1 = 2(UD) and F 1 = 2(FD) it was possible to obtain steady

state solutions. The code took around 4000 time steps to reach steady state. For
F1 = 4 however it was not possible to obtain steady state solution. The code did
not converge to steady state even after 10,000 time steps. Chow (1959) classifies

hydraulic jumps with 2.5 < F1 < 4.5 as oscillating jump. We are not certain if our
inability to achieve steady state solutions for F1 = 4.0 is due to the physical nature

of the oscillating jump. Thus the results presented by us for F1 = 4(UD) and

F 1 = 4(FD) are not steady state solutions but solutions obtained after simulating
the jump for long times (after about 10,000 time steps and 2200yJ/U1 time units).

For F 1 = 4(UD) and F 1 = 4(FD) no separation of boundary layer was observed.

The jump did not exhibit undular behavior and a strong surface roller was found

to occur. In Fig.8.13 ujU1 , -k

pu;
and vr/U1 y1 are shown for F 1 = 4(UD). Just as
in the previous cases the k has large values in the roller region whereas vr achieves

maximum values downstream of the roller. Based on Leutheusser & Alemu (1979)'s

experiments we would expect a larger separation zone for F 1 = 4. Leutheusser &

Alemu (1979) found the length of recirculation to increase with increasing F1 and

decrease with increasing inlet Reynolds number Re 1 • The reasons for the absence of

recirculation zone for F 1 = 4 are not clear. For F 1 = 4(FD) the maximum values of
k and vr were found to be higher than for F 1 = 4(UD). The rest of the qualitative

behavior is the same.

 190

(323.2, 5.17)
,.,.
.....
_,

..........
·""'
_,..
,..
·""
.080

.
-0011
_.,..
_,
_.,.
(245.5, 0.00)
(a) (323.2, o.ool
(323.2, 5.17
.......
.....
·""""'
......
......
.0:.1~

......
.08132

......
.01867
.OlMI!J
.OUT-4.

·""'""
·"""'
(245.5, 0.00) (323.2, o.ool '
(b) (323.2, 5.17
.......
.063'111

......
.04.810

......
·"""'
......
•1130?3

......
.011121
.0168'7
.01163
.......
......
(245.5, 0.00) (323.2, 0.00)
(c)

Figure 8.13: Mean flow and turbulence characteristics for F 1 = 4(UD) (a)
uj U1 ; (b) k/ pUf; and (c )vr /UtYt·
 191

The various parameters such as length of the roller L"' length of the recirculation
zone Lsep and specific energy loss across the jump EL/ £ 1 , for all the four cases we

simulated are listed in Table 8.4. The length of the hydraulic jump Lj is defined to

be the horizontal distance from the toe of the jump to a point downstream where the

influence of the jump has become negligible and flow conditions are determined by

the characteristics of the channel alone. This definition for LJ has been proposed by

Leutheusser & Alemu ( 1979) and found the length of the jump to be approximately
12 to 15 times the downstream subcritical flow height y2 . Determination of the

cross-section downstream where the influence of the jump is negligible is not very
easy and consequently we arbitrarily fixed the jump length to be 15y2 • The actual

jump length is lesser than 15Y2·

The specific energy E at a cross-section is defined as follows

1 lay 2
E = y +-
2g 0
u dy
I
(8.15)

The loss of specific energy EL = £ 1, - £ 2, assummg uniform flow, hydrostatic

pressure and negligible wall friction can be shown to be to be dependent on the

supercritical Froude number as (Rajaratnam 1967)

1 (J1+8Ff,-3?
(8.16)
8 (2 + Ff,)( J1 + 8Ff,- 1)
The values shown in the brackets for ELf £ 1, in Table 8.4 are those obtained from
Eq.8.16. The specific energy loss across a hydraulic jump is more if the flow at the

inlet is fully developed and is consistent with Leutheusser & Alemu (1979)'s finding.

Experimental measurements indicate that actual EL is smaller than that predicted

by Eq.8.16 for F 1 < 2 (Chow 1959). This explains the difference in EL/ E1' obtained

bv us and that given by Eq.8.16 for F1 = 2.

 192

The variation of skin friction coefficient C1 = T w/ ~pU[, is shown in Fig.8.14. The

value of cf upstream of the jump is approximately 0.0026 and decreases to about

0.0005 downstream of the jump. From experimental data, Rajaratnam ( 1967) found

the cf to decrease from about 0.0037 at the beginning of the jump to about 0.0001

near the end of the jump for the whole range of Froude numbers from 3.20 to 9.78.

The values for C1 obtained by us are lower than those presented by Rajaratnam.

This could be due to the assumption of hydraulically smooth wall by us.

8.5 Concluding remarks

A numerical model capable of simulating two-dimensional open-channel flows with-

out any depth averaging has been developed. The a priori unknown air-water ip-

terface is handled through a mixed Lagrangian-Eulerian procedure called the ALE

formulation which handles fairly large free surface deformation without appreciable

mesh distortion. The turbulence is modeled using a two equation k- E model which

gives us the ability to simulate complex flows involving boundary layer separation.

A very interesting but difficult problem in open-channel flows, namely, the hy-
draulic jump has been simulated using the numerical model developed in this work.

The free surface rises abruptly across the jump region and intense turbulent mixing

is known to occur in this region. Successful simulation of this very difficult problem

shows that the numerical procedure developed in this work is robust, and can handle

high gradient regions.

In the case of hydraulic jump with inlet supercritical Froude number F1 = 2.0,

small recirculation zone was found to exist at the bed in the jump region. Due to

the abrupt rise in the free surface, the flow tends to slow down near the free surface,

with the formation of a surface roller. The recirculation zone near the bed and
 193

..§
0
.
,;
•:
! •••:
,.1\
..8 I•
1\
••
1\
0 I•
~ 1\
0
~~
,,
80

I
0

0
0

Figure 8.14: Variation of skin friction coefficient C1 along the jump

 194

the surface roller at the free surface influence the turbulence field significantly and

dominate the wall generated turbulence. Significant departure of the pressure field

from the hydrostatic pressure is observed in the jump region, showing the necessity

for a full-scale numerical model without depth averaging in the simulation of rapidly

varying flows such as hydraulic jumps and bores.

For hydraulic jumps with higher supercritical inlet Froude numbers, breaking

of the free surface is known to occur with the entrainment of air into the water

phase. Thus, the fluid flow is no longer a single-phase flow, rather it is a two-phase

flow containing both air and water. Thus for successful simulation of high Froude

number flows, two-phase flow models may have to used.

 195

Chapter 9

General Conclusions

9.1 Conclusion

The following conclusions may be arrived at based on the work presented so far:

l. .-\n accurate and efficient numerical procedure has been developed which can

simulate many types of gravity-driven flows. The computational efficiency de-

rives in part from the ALE procedure, which allows us to handle the large free
surface deformations in a straight forward manner without any mesh entan-

glement. The computational efficiency is also in part due to the Chorin-type

projection scheme which decouples the pressure field from the velocity field,

thus allowif!-g us to compute them sequentially. All the results presented in the
thin-film flows and turbulent flows over bedforms, have been obtained using a

SPARC-10 machine. Only the hydraulic jump results have been obtained on

a CRAY-YMP. Quantitative comparisons with experiments in thin-film flows,

establishes the accuracy of the numerical procedure. The comparisons in the

case of turbulent flow simulations have been only qualitative, and some more

rigorous validation may be required to establish the accuracy in the case of

complex turbulent flows.

2. An extensive numerical study of the surface wave instability in isothermal

thin-film flows has been done. Both spatially periodic and non-periodic do-

mains have been considered. Extensive parametric revealed a phenomena in

 196

the finite amplitude wave regimes, namely, the quasi-periodic wave evolutions.
Even though, this type of a behaviour was alluded to previously, in the con-

text of other simplified systems, no systematic study, such as done in this

dissertation has been reported previously. Complex wave interactions such as

wave-splitting and wave-mergers have been simulated. The wave interactions

observed in our numerical simulations are in excellent qualitative agreement

with the experimental observations of Liu & Gollub (1994). Evolution of a

time-periodic disturbance in a a long non-periodic spatial domain has been

simulated and are in excellent qualitative and quantitative agreement with

the experiments of Liu & Gollub (1994). This is especially gratifying, since it

establishes that the numerical code developed in this dissertation can be us~d

to perform "computational experiments."

3. A numerical procedure capable of simulating turbulent fluid flows over arti-

ficial stream-wise periodic bedforms has been developed. In the context of

turbulent flows over bedforms, the work presented in this dissertation is the

first study which does not invoke the rigid-lid approximation, but computes

the equilibrium free surface shape along with the mean and turbulence field.

Thus, the numerical procedure can be used for studying high Froude num-

ber flows also. Our numerical simulations reveal the existence of very large

shear stresses at the tip of the dune due to local accelerations associated with

stream-line bending. This local acceleration is seen to be more when we make

the rigid-lid approximation. The turbulence associated with separation wake

is found to dominate the wall generated turbulence. Consequently, the maxi-

mum turbulence intensities occur at a distance approximately of the order of

dune height away from the bed. The suspended sediment load downstream of
 197

the dune peak is found to be significant, and needs to be taken into account
in the bedform formation theories, at least in some cases.

4. The celebratPd hydraulic jump problem in open-channel flows has been nu-

merically simulated. The numerical simulation of hydraulic jump is very chal-

lenging due to the presence of very high spatial gradients in the jump region,

and is in many respects similar to shock wave in compressible flows. In some

cases a small recirculation zone is found to be present at the foot of the jump.

The mixing layer turbulence associated with the recirculation zone and surface

roller are found to dominate the wall generated turbulence.

9.2 Future Work

Some of the possible future research directions are listed below:

1. In reality, the fluid flows existing in nature are three-dimensional. The numer-

ical code developed in this dissertation can be extended to three-dimensions

in a straight forward manner. However, the limitations arise not from imple-

mentation, but from the excessive computational time and memory required

for full-scale 3-D simulations.

2. In the thin-film flow study, it has been mentioned that the presence of surface

waves on the interface is known in increase heat and mass transfer in both the

liquid and gaseous phases. Thus extending the isothermal numerical procedure

to include heat transfer would help us study this aspect quantitatively.

3. Another area which can be improved upon is the turbulence modeling. The

k- c model used in this work, appears to give reasonably good results without
excessive computational cost. Other closure models that are supposedly more
 198

accurate, are very expensive computational requirements wise. Thus, there

exists a need to develop closure models that model the physics accurately

without being very expensive computationally. The Large Eddy Simulation

(LES) appears to be promising but lot more work needs to be done, before

they could be used for solving real-life fluid flow problems.

 199

Bibliography

[1] Alekseenko, S. V., Nakoryakov, V. Ye. and Pokusaev, B.C., "Wave formation
on a vertical falling liquid film," rUChE J., 31, pp.1446-1460 (1985).

[2] Alfrink, B. J. and van Rijn, L. C., "Two- equation turbulence model for flow
in trenches," J. Hydraulic Engg., 109, pp.941-958 (1983).

[3] Amsden, A. A. and Hirt, C. W., "YAQUI: an arbitrary Lagrangian-Eulerian

computer program for fluid at all speeds." Los Alamos Scientific Laboratory
Report, LA-5100 (1973).

[4] Anderson, D. A., Tannehill, J. C. and Pletcher, R. H., Computational Fluid

Mechanics and Heat Transfer, p.247, McGraw-Hill Book Co., (1984).

[5] Anshus, B. E. and Goren, S. L., "A method of getting approximate solutions
to the Orr-Sommerfeld equation for flow on a vertical wall," AIChE J., 12,
pp.1004-1008 (1964).

[6] A.S.C.E. Task Force on Bed Forms in Alluvial Channels, 1966: "Nomenclature
for bedforms in alluvial channels," Proc. A.S.C.E. J. Hyd. Div., 92, pp.51-64
(1966).

[7] ASCE Task Committee on Turbulence Models in Hydraulic Computations,

"Turbulence Modeling of Surface Water Flow and Transport," J. of Hyd. Engg.,
114, pp.970-1073, part:1-5, (1988).
 200

[8] Atherton, R. W. and Homsy, G. M., ''On the derivation of evolution equations
for interfacial waves," Chem. Engng Commun., bf 2, pp.57-77, (1976).

[9] Bach, P. and Villadsen, J., "Simulation of the vertical flow of a thin wavy film
' .
using a finite element method," Int. J. Heat Mass Transfer, 27, pp.815-827,
( 1984).

[10] Basco, D. R., "Introduction to the rapidly- varied unsteady free- surface flow

computation," USGS, Water Resour. Invest. Report N. 83-4284, U.S. Geological

Survey, Reston, Va. (1983).

[11] Behr, M. A., Franca, L. P. and Tezduyar, T. E., "Stabilized finite element

methods for the velocity-pressure-stress formulation of incompressible flows~"

University of Minnesota Supercomputer Institute Research Report, UMSI 92/8,

January (1992).

[12] Belanger, J. B., "Essai sur la solution numerique de quelques problemes relatifs
au mouvement permanent des eaux courantes (Essay on the Numerical Solution

of Some Problems Realtive to the Steady Flow of Water)," Carilian-Goeury,

Paris, (1828).

[13] Bell, J. B., Colella, P. and Glaz, H. M., "A second-order projection method
for the incompressible Navier-Stokes equations," J. Comp. Phys., 85, 257-283

(1989).

[14] Belytschko, T. B. & Flanagan, D. P., "Finite element methods with user-
controlled meshes for fluid-structure interaction," Comput. Meth. Appl. Mech.

Eng., 33, 669-688 (1982).

 201

[15] Benim, A. C. and Zinser, W., "Investigation into the finite element analysis of

confined turbulent flows using a k - t model of turbulence," Comput. Meth.

A.ppl. Mech. Engg., 51 .507-523, (1985).

[16] Benjamin, T. B., "Wave formation in laminar flow down an inclined plane," J.
Fluid l'vfech., 2, pp.554-574, (1957).

[17] Benney, D. J., "Long waves on liquid films,'' J. Math. fj Phys., 45, pp.150-155,

( 1966).

[18] G. Bidone, "Experiences sur le remous et Ia propagation des ondes (Exper-

iments on backwater and the propagation of waves)," l'vfemorie della Reale

Academia delle Scienze di Torino, Turin, 25, pp.21-112 (1820).

[19] Chan, R. K. -C., "A generalized arbitrary Lagrangian Eulerian method for

incompressible flows with sharp interfaces," J. Comp. Phys., 17, pp.311-331

(1975).

[20] Chang, H. -C., Demekhin, E. A. and Kopelevich, D. I., "Nonlinear evolution

of waves on a vertically falling film," ]. Fluid Mech., 250, pp.433-480 ( 1993).

[21] Chang, H. -C., "Wave evolution on a falling film," Ann. Rev. Fluid Mech., 26,

pp.103-136 (1994).

[22] Cheng, M. and Chang, H. -C., "Competition between subharmonic and side-

band secondary instabilities," Phys. Fluids, 7, pp.34-54 (1995).

[23] Chin, R. W., Abernathy, F. F. and Bertschy, J. R., "Gravity and shear wave

stability of free surface flows. Part 1. Numerical calculations," J. Fluid Mech.,

168, pp.501-513 (1986).

 202

[24] Chorin, A. J., "A Numerical Method for Solving Viscous Incompressible Flow
Problems," J. Comput. Phys., 2, pp.12-26, (1967).

[25] Chippada, S., Ramaswamy, B. and Wheeler, M. F., "Numerical simulation of

hydraulic jump," Int. J. Num. Meth. Engg., 37, pp.1381-1397 (1994).

[26] Chippada, S., Ramaswamy, B. and Joo, S. W., "A full-scale numerical study
of interfacial instabilities in thin film flows," to be submitted to J. Fluid Mech.
(1995).

[27] Chow, V. T., Open-Channel Hydraulics, McGraw Hill Book Co., (1959).

[28] Churchill, S. W., "Viscous Flows: The Practical Use of Theory," Buuterworth
Series in Chemical Engineering, (1988).

[29] Donea, J., Giuliani, S. and Halleux, J. P., " An arbitrary Lagrangian Eule-
rian finite element method for transient dynamic fluid-structure interactions,"
Comp. Meth. Appl. Mech. Eng., 33, 689-723 (1982).

[30] Donea, J., Giuliani, S. and Laval, H., "Finite element solution of the unsteady
Navier-Stokes equations by fractional step method," Comp. Meth. Appl. Mech.
Eng., 30, 53-73 (1982).

[31] Donea, J., "Arbitrary Lagrangian Eulerian finite element methods," m

Computational Methods for Transient Analysis (ed. Belytschko, T. B. &

Hughes, T. J. R.), 474-516 (1983).

[32] Dressler, R. F., "Mathematical solution of the problem of roll-waves in inclined

open channels," Commun. Pure Appl. Math., 2, pp.149-194 (1949).
 203

[:33] Dukler, A. E., ··Characterization, elfects and modeling of the wavy gas-liquid
interface,"' In Progress zn Heat and Mass Transfer ( ed. G.Hetstroni. S.sideman
& J.P. Hartnet), 6, pp.207-23~, Pergamon, (1972).

[34] Dukowicz, J. K. and Ramshaw, J.D., "Tensor viscosity method for convection
in numerical fluid dynamics, .. J. Comput. Physics, 32, pp.71-79, (1979).

[:35] Engelund, F,. "Instability of erodible beds." J. Fluid Mech., 42, pp.225-244
(1970).

[36] Finlayson. B. A.. Numerical methods for problems with moving fronts, Ravenna
Park Publishing Inc., Seattle. Washington, USA, 434-456 (1992).

[37] Floryan, J. M. and Rasmussen, H., "Numerical methods for viscous flows with
moving boundaries," Appl Mech Rev, 42, pp.323-341, (1989).

[38] Fredsoe, J., "On the development of dunes in erodible channels," J. Fluid
Mech., 64, pp.1-16 (1974).

[39] Fulford, G. D., "The Flow of Liquids in Thin Films," Adv. in Chemical Engg.,
5, Academic Press, N.Y., (1964).

(40] Gjevik, B., "Occurrence of finite-amplitude surface waves on flailing liquid

films," Phys. Fluids, 13, pp.1918-1925 (1970).

[41] Glowinski, R., "Splitting methods for the numerical solution of the incompress-

ible Navier-Stokes equations," in Vistas in Applied Mathematics, Optimization

and Software (ed. Balakrishna, A. V., etal), 57-95 (1986).

[42] Gresho, P. M., "On the theory of semi-implicit projection methods for viscous
incompressible flow and its implementation via a finite element method that
 204

also introduces a nearly consistent mass matrix. Part 1: Theory," Int. J. Num.
l\!!eth. Fluids, 11, pp.587-620, (1990).

[43] Gresho, P. M. and Chan, S. T., "On the theory of semi-implicit projection

methods for viscous incompressible flow and its implementation via a finite
element method that also introduces a nearly-consistent mass matrix Part 2:

Implementation," Int. J. Num. Meth. Fluids, 11, 621-660 (1990).

[44] Gresho, P. M., Chan, S. T., Christon, M. A. and Hindmarsh, A. C., "A little
more on stabilized Q1 Q1 for transient viscous incompressible flow." in Pro-
ceedings of the International Conference on Monsoon Validity and Prediction,

Trieste, Italy ( 1994).

[45] Harlow, F. H. and Welch, J. E., "Numerical study of large-amplitude free-

surface motions," Phys. Fluids, 9, pp.842-851 (1966).

[46] Henderson, F. M., Open Channel Flow, Macmillan, (1966).

[47] Hinze, J. 0., Turbulence, McGraw Hill Book Co., (1975).

[48] Hirt, C. W., Cook, J. L. and Butler, T. D., "A Lagrangian Method for Calcu-
lating the Dynamics of an Incompressible Fluid with Free Surface," J. Comput.

Phys., 5, pp.103-124, (1970).

[49] Hirt, C. W., Amsden, A. A. and Cook, J. L., "An arbitrary Lagrangian-
Eulerian Computing Method for all flow speeds." J. Comp. Phys, 14, pp.227-

253, (1974).

[50] Hirt, C. W. and Nichols, B. D., "Volume of fluid (VOF) method for the dy-
namics of free boundaries," J. Comput. Phys., 39, pp.201-255, (1981).
 205

[.51] Ho, L. -W. and Patera. A. T., .. A Legendre Spectral Element Method for
simulation of unsteady incompressible viscous free-surface flows," Camp. Meth.
Appl. ;vlech. Eng., 80, pp.J.S.S-366 ( 1990).

[52] Hooper, A. P. and Grimshaw. R., ·'Nonlinear instability at the interface be-
tween two viscous fluids," Phys. Fluids, 28, pp.37-45 (1985).

[53] Hopf, L., "Turbulenz bei einem Flusse,'' Ann. Phys., Ser. 4, 32, pp.777 (1910).

[.54] Huerta, A. and Liu, W. K., ·'Viscous flow with large free surface motion,"
Comput. :\feth. Appl. Mech. Eng., 69, pp.277-324 (1988).

[.55] Hughes, T. J. R., Liu, W. K. and Zimmerman, T., "The Lagrangian-Eulerian

finite element formulation for incompressible viscous flow," Comput. Meth.

Appl. Mech. Eng., 29, pp.329-349 (1981).

[56] Hughes, T. J. R., Franca, L. P. and Balestra, M., "A new finite element
formulation for computational fluid mechanics: V. Circumventing the Babuska-

Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem

accommodating equal order interpolation," Camp. Meth. Appl. Mech. Eng.,

59, pp.85-99 (1986).

[57] Joo, S. W., Davis, S. H. and Bankoff, S. G., "On falling film instabilities and
wave breaking," Phys. Fluids A, 3(1), pp.231-232 (199la).

[58] Joo, S. W., Davis, S. H. and Bankoff, S. G., "Long-wave instabilities of heated
falling films: two-dimensional theory of uniform layers," J. Fluid Mech., 230,

pp.117-146, (199lb).

[59] Joo, S. W. & Davis, S. H., "Irregular waves on viscous falling films," Chem.

Eng. Commun., 118, pp.lll-123 (1992).

 206

[60] Kawahara, M., Hirano, H. and Tsubota, K., "Selective lumping finite element
method for shallow water flow," Int. J. Num. Meth. Fluids, 2, pp.89-112,
(1982).

[61] Kapitza, P. L. and Kapitza, S. P., "Wave flow of thin layers of a viscous
fluid:III. Experimental study of undulatory flow conditions. Zh. Exp. Teor.
Fiz. 19, p.105, (1949). Also in Collected papers of P. L. Kapitza (ed. D. Ter
Haar), vol.2, pp.690-709, Pergamon, (1965)

[62] Kennedy, J. F., "The mechanics of dunes and antidunes in erodible channels,"
J. Fluid Mech., 16, pp.521-544 (1963).

[63] Kennedy, J. F., "The formation of sediment ripples, dunes, and antidunes,"
Ann. Rev. Fluid Mech., 1, pp.147-168 (1969).

[64] Kheshgi, H. S. and Scriven, L. E., "Disturbed film flow on a vertical plate,"

Phys. Fluids, 30(4), pp.990-997 (1987).

[65] Kim, J. and Moin, P., "Application of a fractional-step method to incompress-

ible Navier-Stokes equations," J. Comp. Phys., 59, 308-323 (1985).

[66] Krantz, W. B. and Goren, S. L., "Stability of thin liquid films flowing down a
plane," Ind. Eng. Chem. Fund., 10, pp.91-101 (1971).

[67] Lacroix, M. and Garon, A., "Numerical solution of phase change problems: an
Eulerian-Lagrangian approach," Num. Heat Transfer, 19, pp.57-78 (1992).

[68] Lacy, C. E., Sheintuck, M. and Dukler, A. E., "Methods of deterministic chaos
applied to the flow of thin wavy films," AIChE J., 37, pp.481-489 (1991).
 207

[69] Lai, C., "Numerical modeling of unsteady open-channel flow," Adv. in Hydro-
sczence, 14, 161-333 (1986).

[70] Launder, B. E. and Spalding. D. B., "The Numerical Compuation of Turbulent

Flows," Comput. Meth. Appl. Mech. Engg., 3, 269-289, (1974).

[71] Le, H. and Moin, P., " An improvement of fractional step method for the
incompressible Navier-Stokes equations," J. Comp. Phys., 92, 369-379 (1991).

[72] Lee, J.-J. and Mei, C. C., "Stationary waves on an inclined sheet of viscous
fluid at high Reynolds and moderate Weber numbers," Private Communication
(1995).

[73] Lemos, C. M., "A simple numerical technique for turbulent flows with free
surfaces," Int. J. Num. Meth. Fluids, 15, pp.127-146, (1992a).

[74] Lemos, C. M., Wave breaking : A numerical study, Lecture Notes in Engineer-
ing, ed. C.A. Brebbia and S.A. Orszag, Springer-Verlag, (1992b ).

[75] Leutheusser, H. J. and Alemu, S., "Flow separation under hydraulic jump," J.
Hyd. Research, 17, pp.193-206 (1979).

[76] Lin, S. P., "Finite-amplitude stability of a parallel flow with a free surface," J.
Fluid Mech., 36, pp.l13-126 (1969).

[77] Lin, S. P. and Wang, C. Y., "Modeling wavy film flows," in Encyclopedia of
Fluid Mechanics (ed. N. P. Chermemisinoff) 1, pp.931-951, Gulf, (1985).

[78] Liu, W. K., Chang, H., Chen, J. & Belytschko, T., "Arbitrary Lagrangian
Eulerain Petrov-Galerkin finite elements for for nonlinear continua," Comp.
Meth. Appl. Mech. Eng., 68, 259-310 (1988).
 208

[79] Liu, J., Paul, J. D. and Gollub, J. P., "Measurements of the primary instabilities
of flim flows," J. Fluid Mech., 250, pp.69-101, (1993).

[80] Liu, J. and Gollub, J. P., "Onset of spatially chaotic waves on flowing films,"
Physical Review Letters, 70, pp.2289-2292, (1993).

[81] Liu, J. and Gollub, J.P., "Solitary wave dynamics of film flows," Phys. Fluids,
6(5), pp.1702-1712 (1994).

[82] Lyn, D. A., "Turbulence measurements in open-channel flows over artificial

bedforms," J. Hyd. Engg., 119, pp.306-326 (1993).

[83] Madsen, P. A., and Svendsen, I. A., "Turbulent bores and hydraulic jumps,"
J. Fluid Mech., 129, pp.1-25 (1983).

[84] Malamataris, N. T. and Papanastasiou, T. C., "Unsteady free surface flows on

truncated domains," Ind. Eng. Chem. Res., 30, pp.2211-2219 (1991).

[85] McLean, S. R. and Smith, J.D., "A model for flow two-dimensional bedforms,"

112, pp.300-317 (1986).

[86] McLean, S. R., Nelson, J. M. and Wolfe, S. R., "Turbulence structure over
two-dimensional bedforms: Implications for sediment transprt," J. Geophysical

Research, 99, No. C6, pp.12,729-12,747 (1994).

[87] Mendoza, C. and Shen, H. W., "Investigation of flow over dunes," J. Hyd.

Engg., 116, pp.459-477 (1990).

[88] Mizukami, A. & Tsuchiya, M., "A finite element method for the three-
dimensional non-steady Navier-Stokes equations," Int. J. Num. Meth. Fluids,

4, 349-357 (1984).
 209

[89] Nakaya, C., "Long waves on a thin fluid layer flowing down an inclined plane,"
Phys. Fluids, 18, pp.1407-1412 (1975).

[90] Nelson, J. M. and Smith, J.D., ":\1echanics of Flow Over Ripples and Dunes,"
J. Geophysical Research, 94, No.C6, pp.8146-8162 (1989).

[91] Nelson, J. M., McLean, S. R. and Wolfe, S. R., "Mean flow and turbulence
fields over two-dimensional bedforms,'' Water Resources Research, 29, No.ll,
pp.3935-3953 ( 1993).

[92] Nusselt, W., "Die oberflachenkondensation des Wasserdampfes," Z. Ver. Deut.

lng., 60, pp.541-569 (1916).

[93] Oden, J. T., "Theory and implementation of high-order adaptive hp meth-

ods for analysis of incompressible viscous flows," in Computational Nonlinear
Mechanics in Aerospace Engineering, (ed. Atluri, S. N.), AIAA Inc., 321-363

(1992).

[94] Orlanski, I., J. Comp. Phys., 21, pp.251 (1976).

[95] Patankar, S. V., Numerical Heat Transfer f3 Fluid Flow, Hemisphere, Wash-
ington, D.C. (1980).

[96] Pierson, F. W. and Whitaker, S., "Some theoretical and experimental observa-
tions of the wave structure of falling iiquid films," Ind. Engng Chem. Fundam.,
16, pp.401-408 (1977).

[97] Pironneau, 0., "On the transport-diffusion algorithm and its applications to
the Navier-Stokes equations," Nuemrische Mathematic, 38, 309-332 (1982).
 210

[98] Portalski, S. and Clegg, A. J., "An experimental study of wave inception on
falling liquid films," Chem. Eng. Sci., 27, pp.1257-1265 (1972).

[99] Pracht, W. E., "Calculating three-dimensional fluid flows at all speeds with an

Eulerian-Lagrangian computing mesh," J. Comp. Phys., 17, 132-159 (1975).

[100] Prokopiou, Th., Cheng, M. and Chang, H. -C., "Long waves on inclined films
at high Reynolds number," J. Fluid Mech., 222, pp.665-691 (1991).

[101] Pumir, A., Manneville, P. and Pomeau, Y., "On solitary waves running down
an inclined plane," J. Fluid Mech., 135, pp.27-50 (1983).

[102] Rahman, M. M., Faghri, A. and Hankey, W. L., "Computation of turbulent

flow in a thin layer of fluid involving a hydraulic jump," J. Fluids Engg., 113,

pp.411-418, (1991).

[103] Rajaratnam, N., Adv. in Hydroscience, 4, 197-280, (1967).

[104] Ramaswamy, B. and Kawahara, M., "Arbitrary Lagrangian-Eulerian Finite

Element Method for unsteady, convective, incompressible viscous free surface
fluid flow," Int. J. Num. Meth. Fluids, 7, pp.1053-1075 (1987a).

[105] Ramaswamy, B. and Kawahara, M., "Lagrangian Finite Element Analysis

Applied to Viscous Free Surface Fluid Flow," Int. J. Num. Methods Fluids, 7,

pp.953-984, (1987b ).

(106] Ramaswamy, B. and Kawahara, M., "An equal order velocity pressure finite el-
ement formulation for solving the time-dependent incompressible Navier-Stokes

equations," Bull. Facul. Sci. f3 Eng., Chuo Univ., pp.63-104 (1987c).

 211

[107] Ramaswamy, B., Jue, T. C. and Akin, J. E., "Semi-implicit and explicit finite
element schemes for coupled fluid/thermal problems," Int. J. Numer. Meth.
Eng., 34, 67.5-696 ( 1992).

[108] Rannacher, R., ·'On Chorin's projection method for the incompressible
Navier-Stokes equations," in The Navier-Stokes Equations II- Theory and
Numerical Methods, (ed. J. G. Heywood eta!), Lecture Notes in Mathematics,
Vol.l-530, Springer, Berlin, ( 1992).

[109] Raudkivi, A. J., ''Bedforms in alluvial channels," J. Fluid Mech., 26, pp.-507-
514 (1966).

[110] Rice, J. G. and Schnipke, R. J., "An equal-order velocity-pressure formulation

that does not exhibit spurious pressure modes," Comp. Meth. Appl. Mech. Eng.,
58, pp.135-149 (1986).

[111] Richards, K. J., "The formation of ripples and dunes on an erodible bed,"
J. Fluid Mech., 99, pp.597-618 (1980).

[112] Richards, K. J. and Taylor, P. A., "A numerical model of flow over sand waves
in water of finite depth," Geophys. J. R. Astron. Soc., 65, pp.103-128 (1981).

[113] Rosenau, P. and Oron, A., "Evolution and breaking of liquid film flowing on
a vertical cylinder," Phys. Fluids A, 1(11), pp.1763-1766 (1989).

[114] Roskes, G. J., "Three-dimensional long waves on a liquid film," Phys. Fluids,
13, pp.1440-1445 (1970).

[115] Rosten, H. I. and Worrell, J. K., "Generalized wall functions for turbulent
flow,' PHOENICS J. Comp. Fluid Dynamics and Applications, 1, pp.81-109,

(1988).
 212

[116] Rouse, H., and Ince, S., History of Hydraulics, Iowa Institute of Hvdraulic
Research, State Univ. of Iowa, Iowa City, Iowa (1957).

[117] Rouse, H., Siao, T. T. and Nagaratnam, S., "Turbulence characteristics of the
hydraulic jump," J. Hyd. Div. ASCE, (1958).

[118] Salamon, T. R., Armstrong, R. C. and Brown, R. A., "Traveling waves on

vertical films: Numerical analysis using the Finite Element Method," Phys.
Fluids, 6, pp.2202-2220, ( 1994 ).

[119] Schneider, G. E., Raithby, G. D. and Yovanovich, M. M., "Finite-element

solution procedures for solving the incompressible, Navier-Stokes equations us-

ing the equal order variable interpolation," Num. Heat Transfer, 1, pp.443-451
(1978).

[120] Shames, I. H., Mechanics of Fluids, 3rd ed. McGraw-Hill, Inc., (1992).

[121] Shaw, C. T., "Using a segregated finite element scheme to solve the incom-
pressible Navier-Stokes equations," Int. J. Num. Meth. Fluids, 12, pp.81-92
(1991).

[122] Shen, J., "On Error Estimates of Projection Methods for Navier-Stokes Equa-
tions: First-Order Schemes," SIAM J. Numer. Anal., 29, pp.57-77, (1992).

[123] Shields, A., "Anwendung der Ahnlichkeits-Mechanik und der Turbulenz-

Forschung auf die Geschiebe-bewegung," Mitt. Preuss. Versuchsanst. Wasser

Schiffsbau, 26, pp.1-26 (1936).

[124] Silliman, W. J. and Scriven, L. E., J. Comput. Phys, 34, 287 (1980).
 213

[125] Sivashinsky, G. I. and Michelson, D. M., "On the irregular wavy flow of liquid
film down a vertical plane," Prog. Theor. Phys., 63, pp.2112 ( 1980).

[126] Soulaimani, A., Fortin, M., Dhatt, G. and Ouellet, Y., "Finite element simu-
lation of two- and three-dimensional free surface flows," Computer Methods in
Applied Mechanics and Engineering,' 86, pp.265-296 (1991 ).

[127] Stainthorp, F. P. and Allen, J. M., "The development of ripples on the surface
of liquid film flowing inside a vertical tube," Trans. Inst. Chem. Engrs. 43,

pp.85-91 (1965).

[128] Temam, R., On the theory and numerical analysis of the Navier-Stokes equa-
tions, North-Holland, Amsterdam (1971).

[129] Van Kan, J., "A second-order accurate pressure correction scheme for viscous
incompressible flow," SIAM J. on Scientific & Statistical Computing, 7, 870-

891 (1986).

[130] van Rijn, L. C., "Sediment Transport, Part 1: Bed Load Transport," J. of
Hyd. Engg., 110, pp.l431-1456, (1984a).

[131] van Rijn, L. C., "Sediment Transport, Part 2: Suspended Load Transport,"
J. of Hyd. Engg., 110, pp.1613-1642, (1984b).

[132] van Rijn, L. C., "Sediment Transport, Part 3: Bed Forms and Alluvial Rough-
ness," J. of Hyd. Engg., 110, pp.1733-1753, (1984c).

[133] Whitaker, S., "Effect of surface active agents on stability of falling liquid

films," Ind. Engng Chem. Fundam., 3, pp.132-142 (1964).

[134] White, F. M., Viscous Fluid Flow, McGraw-Hill, (1991).

 214

[135] Yih, C. -S., "Stability of parallel laminar flow with a free surface," Proc. 2nd
US Congr. Appl. Mech., pp.623-628, ASME, (1955).

[136] Yih, C. -S., "Stability of liquid flow down an inclined plane," Phys. Fluids,
6, pp.321-324, (1963).

[137] Zienkiewicz, 0. C. and Wu, J., "Incompressibility without tears- how to avoid
restrictions of mixed formulation," Int. J. Num. Meth. Eng., 34, pp.ll89-1203

(1991).

[138] Zienkiewicz, 0. C. and Wu, J., "A general explicit or semi-explicit algo-
rithm for compressible and incompressible flows," Int. J. Num. Meth. Eng.,

35, pp.457-479 (1992).