Continuation of a single layer flow over an obstacle

Continuation of a single

layer flow over an obstacle

Jan Willem Greidanus

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

Bachelor Thesis in Applied Mathematics

Supervisor(s): F.W. Wubs October 2008

Continuation of a single

layer flow over an obstacle

Jan Willem Greidanus

Supervisor(s):

F.W. Wubs

Institute of Mathematics and Computing Science P.O. Box 407

9700 AK Groningen

The Netherlands

Contents

1 Introduction 1

2 Basic Equations 3

2.1 Derivation of the Shallow Water Equations . . . 3

2.2 The Conservative Form . . . 4

2.3 Froude Number . . . 5

2.4 Solution to the Linear Hydrostatic Flow Problem . . . 6

2.5 Eigenvalue Problem . . . 7

3 Numerical Model 9 3.1 Space Discretization . . . 9

3.2 Solution Method . . . 10

3.3 Matlab Program . . . 11

3.4 Generalized Eigenvalue Problem . . . 11

4 Numerical Results 13 4.1 Numerical Eigenvalue Spectra for the Linear Hydrostatic Flow Problem . . . 13

4.2 Solutions Of The Model . . . 14

4.3 Subcritical To Supercritical . . . 15

5 Conclusions 17

iii

iv CONTENTS

Contents

0

Chapter 1

Introduction

In this text I will investigate steady single layer flow over an obstacle of which the height will be increased slowly to see what this does to the solutions and eigenvalues of the model. If the obstacle is not too high we will have completely subcritical flow. However, at a certain height the flow may become supercritical.

First I will introduce the basic equations for this problem. These will be derived from the Navier-Stokes equations. Using some assumptions and simplifications we arrive at the equations for our specific problem. After that I will derive the analytical solution to the linearized equations and look at what this solution can learn us about the transition from sub- to supercritical flow. This means we will look at the eigenvalue problem for a perturbation of the linear equations.

Furthermore we take a look at the model that we need for the numerical analysis and what methods are used. We can come up with an eigenvalue problem related to the discretized equations and this is important for the numerical research. After that I will investigate numerical solutions of the equations and look what happens with the eigenvalues when we go from sub- to supercritical flow

1

2 CHAPTER 1. INTRODUCTION

Chapter 2

Basic Equations

2.1 Derivation of the Shallow Water Equations

For research in fluid dynamics one usually starts with the Navier-Stokes equations. For a Newtonian fluid these equations are given by: The momentum equation

ρDu

Dt = −∇p + ρF + µ∆u (2.1)

and the continuity equation:

∂ρ

∂t + ∇ · (ρu) (2.2)

Here u = (u, v, w) is the velocity in Cartesian directions (x, y, z). Throughout this text incompressible flow will be considered, hence ^∂ρ_∂t = 0 and the continuity equation can be simplified to

∇u = 0

When working with an incompressible homogeneous inviscid fluid we get:

Du Dt = −1

ρ∇p + F where F will be given by gravity, so

Du Dt = −1

ρ∇p − gz (2.3)

Here ^Du_Dt is the material derivative:

Du Dt = ∂u

∂t + (u · ∇)u

z is in the vertical upward direction, furthermore a bottom topography is considered to be given by z = h(x, y). Because this topography or bottom profile has a considerable influence on the flow properties we will make different choices for the maximum of h. The boundary condition on this bottom for inviscid fluid, is that the velocity normal to h is zero, this can be written as

u ·



 1 h_x h_y



 = 0 (2.4)

3

4 CHAPTER 2. BASIC EQUATIONS on z = h(x, y). On the free surface η = d−d₀+h which is the displacement from the mean level z = d₀ of the free surface, we need two boundary conditions. Because the mass of any fluid above this surface is negligible, we can take the pressure at the surface to be constant, so p = p₀(hydrostatic pressure). We can now use Bernoulli’s Law and state the dynamic free surface boundary condition

· u v

¸

t

+1

2(u²+ v²) + gη = 0 (2.5)

The second boundary condition on the free surface is the kinematic free surface condition, it guarantees that there will be no fluid above the free surface. In mathematical terms this is

η_t+ uη_x+ vη_y = w (2.6)

In any flow it is possible to write the equations in terms of vorticity rather then velocity, this is done by taking the curl of the equation. (may have to apply the 2nd Fundamental Theorem of Calculus to do this rigorously ) in our case this gives

Dω

Dt = ω · ∇u

In our case the flow is initially rotation free (ω = 0) and that means it will be so throughout.

If everywhere the vertical accelerations ^D_Dt^2w2 are much less then gravity the flow can be called hydrostatic, so we have −_ρ¹

0

∂p

∂z−g = 0. We call the scale of the horizontal fluid motions L, this means that the above approximation holds if (^d_L⁰)² << 1. The pressure within the fluid will then be given by p = p₀+ ρ₀g(d₀+ η − z) and so the horizontal pressure gradient is independent of z, so if u and v are initially independent of z, this will be the case throughout.

We can now state that the following set of equations describe the motion of the flow:

u_t+ uu_x+ vu_y = −gη_x (2.7)

v_t+ uv_x+ vv_y = −gη_y (2.8)

η_t+ (du)_x+ (dv)_y = 0 (2.9)

Where d = d₀+ η − h is the thickness of the fluid layer.

2.2 The Conservative Form

In this text I will consider a flow of a layer of fluid of thickness d(x, t) with velocity u(x, t).

The bottom will be given by h(x) in which the flow will encounter an (long) obstacle, in the middle of the region of interest, with given height and shape. The Width of the channel is considered equal at every point in the flow. η Is the displacement from the equilibrium depth d₀. These assumptions lead to the following differential equations:

u_t+ uu_x= −gη_x (2.10)

d_t+ (du)_x= 0 (2.11)

2.3. FROUDE NUMBER 5 The second is already in conservative form. To get the first equations into conservative form, we multiply (2.10), (2.11) by d and u respectively, so we have

du_t+ (du)u_x+ gdη_x = 0 (2.12)

d_tu + u(du)_x= 0 (2.13)

The reason to use this conservative form [WIKI] is that they can deal with phenomena that normally break down the assumptions essential for the usage of the shallow water equations, like hydraulic jumps which are possible in this setting [BAI03]. when we add (2.12) and (2.13) and use (du)_t= du_t+ ud_tand (du²)_x = (du)u_x+ u(du)_x= (du)u_x+ udu_x+ u²d_x. We end up with

(du)_t+ (du²)_x+ gdη_x= 0 (2.14)

or using η = d − d₀+ h and making the equation time independent (du²)_x+ (1

2gd²)_x+ gdh_x = 0 (2.15)

For the time independent case equation 2.11 is simply equal to

∂du

∂x = 0 (2.16)

or

u_rd_r= u_ld_l

Where u_r and d_r are the velocity and depth at a point right and u_l and d_l are the velocity and depth at a point to the left. Those points can be any points in the x direction, because of the steady flow we consider. This means we can easily choose the left and right points exactly at our mesh points. This equation guarantees that the same amount of fluid passes through every point all the time.

2.3 Froude Number

The so-called Froude Number F , named after William Froude (1810-1879), is an important dimensionless parameter, for shallow flow over an obstacle. It has the form F = ^√U_gl with U the fluid speed, g the gravitational acceleration and l is a length. There are different choices for the length l. It may be the length or the height of the obstacle. Another possibility is to take the depth of the stream. These different choices for l lead to different dynamical interpretations of the Froude Number. This is in contradiction to the Reynolds number in which also different lengths can be used, but the interpretations are the same, namely viscous flow for small Reynolds number and turbulent for large Reynolds numbers analogous to sub- and supercritical flow for the Froude Number.

In this study I will use the local depth of the stream d for the length in the Froud Number, however it was never used in this way by William Froude, it is the most common choice in hydraulics of streams. In our single-layer flow it then represents the local fluid speed divided by the local significant wave speed. Hence the Froude number depends on the position in the computational domain. The Froude number thus gives a measure to what extend long linear gravity waves will propagate against the stream(First derived by Lagrange).

In general flow with a Froude number F below one is called subcritical, F equal to one critical and F larger then one supercritical.

6 CHAPTER 2. BASIC EQUATIONS

2.4 Solution to the Linear Hydrostatic Flow Problem

We consider a uniform stream of magnitude U in the x-direction. In this flow we insert a long obstacle with small height, with topography given by z = h(x, y). We linearize the equations about the initial state (u = U + u⁰) and get:

u⁰_t+ U u⁰_x = −gη_x v_t+ U v_x = −gη_y η_t+ U η_x+ d₀(u⁰_x+ v_y) = U h_x

For the one-dimensional problem we can drop the second equation and take v_y = 0 thus for an even bottom we have

u_t= −U u_x− gη_x (2.17)

d_t= −(du)_x= −U d_x− d₀U_x (2.18) We can write this system of differential equations in matrix vector notation

· u d

¸

t

= −

· U g d₀ U

¸ · u d

¸

x

(2.19)

This system is diagonizable, but by symmetrizing the matrix we can easily find the resulting scalar equations.

"

1 0

0 qg

d0

# · U g d₀ U

¸ " 1 0 0

qd0

g

#

=

· U √

gd₀

√gd₀ U

¸

(2.20)

Thus the system can be written as

"

qu

g d0d

#

t

= −

· U √

gd₀

√gd₀ U

¸ " qu

g d0d

#

x

(2.21)

hence

(u + r g

d₀d)_t= −(U +p

gd₀)(u + r g

d₀d)_x (2.22)

(u − r g

d₀d)_t= −(U −p

gd₀)(u − r g

d₀d)_x (2.23)

so

(u + r g

d₀d) = f (x − (U +p

gd₀)t) (2.24)

(u − r g

d₀d) = g(x − (U −p

gd₀)t) (2.25)

2.5. EIGENVALUE PROBLEM 7

2.5 Eigenvalue Problem

Because we are interested in the perturbation of steady states we will now take a look at the eigenvalue problem related to (2.17) and (2.18). In matrix vector notation this is

· −λ + U_∂x^∂ g_∂x^∂ d_0∂x^∂ −λ + U_∂x^∂

¸ · u⁰ η⁰

¸

=

· 0 0

¸

(2.26) with u⁰ and η⁰ the perturbations to the steady state. We now have a system of differential equations for these perturbations, with boundary conditions:

u⁰(0) = 0 (2.27)

η⁰(L) = 0 (2.28)

Eliminating u⁰ we obtain

[(−λ + U ∂

∂x)²− gd₀( ∂

∂x)²]η = 0 (2.29)

and now we try solutions of the form η = e^γx

(−λ + U γ)²− gd₀γ²= 0 ⇔ (U γ − λ)²= gd₀γ² (2.30)

⇔ U γ − λ = ±p

gd₀γ ⇔ γ = λ U ±√

gd₀ ⇔ U ±p

gd₀ = λ

γ (2.31)

The latter equality we can use to write down a general solution for η, we have η(x) = Ae

λ U +√

gd0x

+ Be

λ U −√

gd0x

(2.32)

and ∂η(x)

∂x = A λ U +√

gd₀e

λ U +√

gd0x

+ B λ

U −√ gd₀e

λ U −√

gd0x

(2.33) from now on we will substitute

α₁ = λ U +√

gd₀ (2.34)

α₂ = λ U −√

gd₀ (2.35)

A and B are unknown up to here. Now we can use the second equation of the system (2.26) to get a system of equations for A and B. We have

d₀∂u⁰

∂x − λη⁰+ U∂η⁰

∂x = 0 (2.36)

or equivalent

d₀∂u⁰

∂x = λη⁰− U∂η⁰

∂x (2.37)

we can now use (2.32) and (2.33) in the last equation, hence

∂u⁰

∂x = λ

d₀(Ae^α1x+ Be^α2x) −λU

d₀ (Aα₁e^α1x+ Bα₂e^α2x) (2.38)

8 CHAPTER 2. BASIC EQUATIONS Since this is a first order differential equation we can solve by separation of variables:

Z du⁰=

Z λ

d₀(Ae^α1x+ Be^α2x) (2.39)

−λU

d₀ (Aα₁e^α1x+ Bα₂e^α2x)dx After integration we have two equations

u⁰ = λ d₀(A 1

α₁e^α1x+ B 1

α₂e^α2x) (2.40)

−λU

d₀ (Aα₁ 1

α₁e^α1x+ Bα₂ 1

α₂e^α2x) + constant u⁰=

r g

d₀(Ae^α1x− Be^α2x) + constant (2.41) Because u⁰(0) = 0 we have

A − B = 0 (2.42)

Also we fill in the boundary condition for η⁰(L) = 0 and we have Ae

λL U +√

gd0 + Be

λL U −√

gd0 = 0 (2.43)

Combining the last two equations we have a system of the equations for the unknowns A and

B "

1 −1

e

λL U +√

gd0 e

λL U −√

gd0

# · A B

¸

=

· 0 0

¸

(2.44) The determinant of this system is

e^α1L+ e^α2L (2.45)

and a non-trivial solution is found when the determinant is zero, so

e^α2L= −e^α1L (2.46)

or

e^(α2−α1)L= −1 = e^{−i(π+2kπ)} ∀k ∈ Z (2.47) Since

α₂− α₁ = 2λ√ gd₀

U²− gd₀ (2.48)

hence

2λ√ gd₀

U²− gd₀ = i(π + 2kπ)

L (2.49)

Now we rewrite the last equation for λ 2λp

gd₀ = i(π + 2kπ)(U²− gd₀)

L (2.50)

λ = i(π + 2kπ)(U²− gd₀) 2L√

gd₀ (2.51)

Hence the Eigenvalues are purely imaginary, for both sub- and supercritical flow. In the numerical analysis we will take a look at values for u and d, for which U²− gd₀ almost, since this means we have eigenvalues almost zero.

Chapter 3

Numerical Model

3.1 Space Discretization

We discretize equations 2.14 and 2.16 in space in the middle of two successive grid points, i and i + 1. The situation is visualised in 3.1. In the point i, u_i and d_i are given. The middle point i +¹₂ is where f_i and f_i+n are calculated and in i + 1, u_i+1 and d_i+1 are the variables.

We can write the equation as an integral from x_i to x_i+1.

Figure 3.1: the computational grid

Z _xi+1

xi

((du)_t+ (du²+1

2gd²)_x+ dh_x)dx = 0 (3.1) after integration this is

(du)_t∆x + (du²+1

2gd²)_i+1− (du²+ 1

2gd²)_i+ dh_x∆x = 0 (3.2) ((du)_t)_i+¹

2∆x + (du²+1

2gd²)_i+1− (du²+1

2gd²)_i+ dh_x∆x = 0 (3.3) Time independent we have

f_i = d_i+1· u²_i+1+1

2g_i+1d²_i+1− d_i· u²_i −1

2g_id²_i +1

2g(d_i+1+ d_i)h_x∆x (3.4) We can treat the second equation in a similar way:

0.5∆x∂

∂t(d_i+1+ d_i) = −u_i+1d_i+1+ u_id_i (3.5) or without the time derivative

f_i+n= u_id_i− u_i+1d_i+1 (3.6)

9

10 CHAPTER 3. NUMERICAL MODEL Of course we need to use boundary conditions in order to be able to work with this discretiza- tion. On the left side of the grid the velocity is prescribed and on the right the depth. This states that in the first, n-th (n + 1)-th and 2n-th equation we can simply substitute the given value for the either u or d.

The bottom profile h also appears in the equations and we have taken this to be zero everywhere except for the obstacle. The obstacle has a the form of a cosine. The obstacle is placed in the middle of the grid and covers halve of the grid length L. The mathematical shape is defined as:

h(x) = 1 − cos(2π(x − 0.25L)

0.5L ) (3.7)

3.2 Solution Method

At this point we have obtained a system of linear equations.

F (x) =





 f₁ f₂ ... f_2n





= 0

In order to solve the system of equations we are going to use a n-dimensional Newton Method [BUFA05]. This means we are looking for a fixed point of a multiple valued function G(x).

Where G is defined as:

G(x) = x − J(x)⁻¹F (x) (3.8)

The matrix J is the Jacobian matrix as defined in equation 3.9:

J(x) =









∂f1

∂x1

∂f1

∂x2 · · · _∂x^∂f_2n¹

∂f2

∂x1

∂f2

∂x2 · · · _∂x^∂f_2n² ... ... . .. ...

∂f2n

∂x1

∂f2n

∂x2 · · · ^∂f_∂x²ⁿ_2n







 (3.9)

Here x contains both the velocity in x-direction and the depth at every point on the grid, so x^T = [u, d] = [u₁, u₂, . . . u_n, d₁, d₂, . . . , d_n]

hence the Jacobian matrix in our case is











∂f1

∂u1

∂f1

∂u2 · · · _∂d^∂f_n¹

∂f2

∂u1

∂f2

∂u2 · · · _∂d^∂f_n² ... ... . .. ...

∂f2n

∂u1

∂f2n

∂u2 · · · ^∂f_∂d²_nⁿ











If we now calculate all the partial derivatives in this Jacobian we can program a working Newton method. These derivatives for the equations will be:

∂f_i

∂u_j = (−1)^i+j+12d_ju_j (3.10)

3.3. MATLAB PROGRAM 11 fori = 1, 2, . . . , n, j = i, i + 1.

∂f_i

∂d_j = (−1)^i+j+1d_j (3.11)

for i = 1, 2, . . . , n, j = i, i + 1.

∂f_i

∂u_j = (−1)^i+j+1(u²_j + gd_j+ 1

2g∆xhx_j) (3.12)

for i = n + 1, n + 2, . . . , 2n, j = i, i + 1.

∂f_i

∂d_j = (−1)^i+j+1u_j (3.13)

for i = n + 1, n + 2, . . . , 2n, j = i, i + 1.

In most cases the derivatives are zero, because f_iwhich is defined between every two successive grid points only depends on the depth, fluid velocity and bottom profile in these points. To find the fixed point we iterate in x and we have the following iteration scheme

x^(k)= G(x^(k−1)) = x^(k−1)− J(x^(k−1))⁻¹F (x^(k−1))

Here k is the iteration number. When the starting values are chosen accurate enough and if J⁻¹ exists then this method gives quadratic convergence. We do not want to calculate the inverse of the Jacobian matrix explicitly, because this is a quite expansive operation. To do this we divide the method in two steps. First the system J(x^(k−1))dx = −F (x^(k−1)). Second this new solution is added to the existing one: x^(k) = x^(k+1)+ dx The stopping criterion is when the norm of dx is small enough.

3.3 Matlab Program

The program consists of five functions next to the main. One of those creates the bottom gradient. This is done in a way that guarantees we have obstacle with the shape of half a period of a sine. Then we have a function that creates the right hand side. Furthermore we have the a function for the Newton method that uses another to create the Jacobian matrix.

3.4 Generalized Eigenvalue Problem

When we take a look at the generalized eigenvalues of the problem we need to take into account the space discretization of the time dependent part. This is because when there is a perturbation we have to take into account the time dependent part of the equations again. We split d and u in a steady and an unsteady part. For the time derivative in the first equation (3.3) we have,

∂

∂t(M (du))_i+¹

2 (3.14)

M is the matrix that takes the mean value from the surrounding grid points, so we can rewrite

as ∂

∂t(M (du))_i+1

2 = ∂

∂t

((du)_i+1+ (du)_i)

2 (3.15)

12 CHAPTER 3. NUMERICAL MODEL Now we use the fact that

ab + cd = (a + c)(b + d)

2 + (a − c)(b − d) 2 so

∂

∂t

((du)_i+1+ (du)_i)

2 = ∂

∂t((d_i+1+ d_i)(u_i+1+ u_i)

4 +(d_i+1− d_i)(u_i+1− u_i)

4 ) (3.16)

and in matrix notation we have

∂

∂t(M dM u + DdDu) (3.17)

D is the matrix of the difference operator ((Dd)_i+1

2 = ^di+1₂^−di). The time derivative in de second equation 3.5 gives

∂

∂tM d (3.18)

After applying the product rule, the discretization of the time dependent parts will be re- spectively

M uM∂d

∂t + M dM∂u

∂t + DuD∂d

∂t + DdD∂u

∂t (3.19)

M∂d

∂t (3.20)

so the matrix A is defined as A = ∆x

· M dM + DdD M uM + DuD

0 M

¸

(3.21) Herewith we can write the discretized problem as

A(x)∂x

∂t = F (x) (3.22)

with x = [u, d] and A a matrix consisting of the difference and mean operator as mentioned before. F is given by equations 3.4 and 3.6. We take x = x_s+ v, with x_s the steady solution and v a small perturbation, because x_s is independent of t we have

A(x_s+ v)∂(x_s+ v)

∂t = (A(x_s) + A(v))∂v

∂t (3.23)

A(x_s+ v)∂v

∂t = F (x_s+ v) = F (x_s) + Jv = Jv (3.24) Because v is small, with respect to x_s we simplify to,

A(x_s)∂v

∂t = Jv (3.25)

We take

v = ˜ve^λt (3.26)

To find the eigenpair (˜v, λ), where ˜v = [˜u, ˜d]^T only depends on the space variables. Thus the generalized eigenvalue problem is

λA(x_s)˜v = J ˜v (3.27)

Chapter 4

Numerical Results

4.1 Numerical Eigenvalue Spectra for the Linear Hydrostatic Flow Problem

We will now try to find out what happens when parameters have critical values. First we look at the solutions and eigenvalues where U²− gd₀ is close to zero, because this is the case where the eigenvalue λ of the eigenvalue problem from section 2.5 is passing through zero. In this experiment we are not working with the obstacle. In our model we have taken g = 10 and for this experiment we keep the depth at the right side d₀ = 1, meaning we have a problem for U =√

10.

In the figure below one can see the eigenvalue spectra for U =√

10 − 0.0001 + 0.00011p for p ∈ 0, 1, 2 (4.1)

−1500 −1000 −500 0 500 1000 1500

−1500

−1000

−500 0 500 1000 1500

Re eigenvalues

√ 10− 0.0001

√ 10 + 0.00001

√10 + 0.00012

−0.01 −0.008 −0.006 −0.004 −0.002 0 0.002 0.004 0.006 0.008 0.01

−0.01

−0.008

−0.006

−0.004

−0.002 0 0.002 0.004 0.006 0.008 0.01

Re eigenvalues

√ 10− 0.0001

√ 10 + 0.00001

√10 + 0.00012

Figure 4.1: Eigenvalues of the generalized eigenvalue problem, the right picture gives a more detailed image of the siutation

In this figure there does not seem to be much of a difference for the different values of U , however if we take a closer look at the part around zero we get figure 4.1. What we in the right figure is that for values of U smaller then √

10 the eigenvalues in this part are in the shape of an arrow to the right and for U larger then√

10 this points in the opposite direction.

The eigenvalues go through the imaginary axis in this experiment, and this means the real parts go through zero, which indicates the solutions to be unstable.

In fact the continuous equation would give a whole bunch of eigenfunctions with eigenvalue zero. Here we see a clustering near zero because U is not precisely√

10. Moreover, in contrast 13

14 CHAPTER 4. NUMERICAL RESULTS to continuous equations there is some dissipation in the numerical model which becomes negative for some components if U >√

10.

4.2 Solutions Of The Model

If we run our program with some prescribed parameters we get solutions for u and d. These solutions are plotted and solutions for the depth d can be displayed as follows. In this figure

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

Figure 4.2: solution

the blue line represent the water level and the black one the bottom topography. For all experiments the length of the interval has length 10 and we choose an odd number of grid points, so the top of the obstacle is really one point and has the wanted height.

4.3. SUBCRITICAL TO SUPERCRITICAL 15

4.3 Subcritical To Supercritical

As stated in section 2.2 the Froude Number is an important quantity in our flow problem.

We will now look at the situation where we add an obstacle in the flow and slowly raise the height of this object, until we have a Froude Number larger than one for the first time on our interval of interest. The interval is divided in 31 grid points all at equal distances. When we start the solutions and local Froude number are as visible in the figures below.

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 0.2 0.4 0.6 0.8 1

F

Figure 4.3:

If we continue using our previous solution every time as the starting vector for the Newton method it slowly changes. The depth d gets higher in the first part of the interval and the Froude number F gets close to one at the end of the interval.

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 0.2 0.4 0.6 0.8 1

F

Figure 4.4:

From here on the steps with which the obstacle is raised are kept ten times smaller, so we can get a closer look at what’s the matter.

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 0.2 0.4 0.6 0.8 1

F

Figure 4.5:

16 CHAPTER 4. NUMERICAL RESULTS Now if we go on one more step in raising the obstacle height we have the following situation:

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0 1 2 3 4 5 6 7 8 9 10

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2

F

Figure 4.6:

In the figure below one can see eigenvalues for four solutions just before and after super- critical flow occurs. Blue are the eigenvalues corresponding to the solution in figure 4.5 and red for figure 4.6. We can see in (4.7) that for F > 1 some eigenvalues get a positive real part.

−3500 −3000 −2500 −2000 −1500 −1000 −500 0 500

−3000

−2000

−1000 0 1000 2000 3000

Re eigenvalues

subcritical subcritical subcritical supercritical

−500 −400 −300 −200 −100 0 100 200 300 400 500

−500

−400

−300

−200

−100 0 100 200 300 400 500

Re eigenvalues

subcritical subcritical subcritical supercritical

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2

x 10−3

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2x 10−3

Re eigenvalues

subcritical subcritical subcritical supercritical

Figure 4.7: Eigenvalues

First the tend to zero and when for F = 1 there is clustering near zero, and when the obstacle is raised again the spread out to the right of the imaginary axis. The number of eigenvalues that cluster to zero and get a positive real part is equal to the number of grid points where F > 1. If we repeat this research with 61 grid points (double the amount of before more or less) we can see from the figure 4.8 that the number of grid points where F > 1 is equal to (plus or minus one) the number of eigenvalues with positive real part. One can see that if we go on continuing after we already have F > 1 the eigenvalues to the right of the imaginary axis are spreading out, e.a. the are getting larger in norm.

−14000 −12000 −10000 −8000 −6000 −4000 −2000 0 2000

−1.5

−1

−0.5 0 0.5 1 1.5x 10⁴

Re eigenvalues

subcritical subcritical supercritical supercritical

−0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1

−0.1

−0.08

−0.06

−0.04

−0.02 0 0.02 0.04 0.06 0.08 0.1

Re eigenvalues

subcritical subcritical supercritical supercritical

Figure 4.8: Eigenvalues

Chapter 5

Conclusions

We started with the Navier-Stokes equations and derived the shallow water equations from them. This involves a number of assumptions, like incompressibility. Because it has usefulness in numerical computations we have written the equations for the one dimensional case in conservative form. The Froude Number is an important dimensionless quantity, because its value indicates whether we are dealing with a sub- or supercritical flow. In our interpretation the Froude number is defined locally, and its value is: the local fluid velocity divided by the local significant wave velocity. It is possible to come up with an analytical solution for the equations if we linearize them first.

We have also taken a look at the eigenvalue problem related to a perturbed system. We discovered the eigenvalues should be purely imaginary for sub- and supercritical flow. And the critical point is when U² = gd₀, because then we have eigenvalues zero, i.e. Froude number 1.

For the discretization we have used the equations in conservative form and we discretized in the middle of two successive grid points. From this we got the discretized equations. Because we want to take a look at the eigenvalues of our numerical model we came up with the generalized eigenvalue problem. Therefore we need to isolate the discretized time dependent parts. This eigenvalue problem is solved in the Matlab program.

Now we take a look at our numerical results. First we considered the situation without obstacle, and looked at the situation we came up with in section 2.5 where U²− gd₀ passes through zero. Some eigenvalues cluster to zero and after U²− gd₀ changes sign they spread out again. Furthermore we observe that the eigenvalues are not purely imaginary as expected from section 2.5, this is attributed to numerical dissipation.

If the obstacle is introduced we can do a same kind of analysis, by raising it slowly.

Initially all eigenvalues are in the left half plane. If we have F close to one somewhere on the interval then part of the eigenvalues of the generalized eigenvalue problem cluster and tend to zero. If we raise the obstacle further F becomes greater than one and the eigenvalues spread out again in the same shape, but now with positive real part. The number of eigenvalues for which this happens is equal to the number of grid points where we have F > 1 plus or minus one. It might be that also here the real part of the eigenvalues is caused by artificial diffusion. This is something left for future study. The main conclusion is that for both the even and uneven bottom case something similar occurs: part or all of the eigenvalues tend to zero. In conventional continuation we have that only 1 or a pair of eigenvalues pass the imaginary axis and standard techniques from bifurcation theory are known for handling that.

17

18 CHAPTER 5. CONCLUSIONS Here an infinite number of eigenvalues becomes zero and the textbooks do not say how to handle that. This is something for further research.

Bibliography

[BAI95] Baines, Peter G. Topographic effects in stratified flows. Cambridge University Press. 1995

[BAI03] Baines, Peter G./Whitehead J.A. On multiple states in single layer flows. in:

Physics Of Fluids volume 15, number 2, February 2003. American Institute of Physics

[BUFA05] Richard L. Burden, J. Douglas Faires Numerical Analysis 8th edition. Thomson Brooks/Cole. 2005

[HOO97] H.W. Hoogstraten Stromingsleer. Rijksuniversiteit Groningen, vakgroep Wiskunde, Oktober 1997

[WIKI] http://en.wikipedia.org/wiki/Shallow_water_equations

19