Long waves generated by ships moving in shallow water

(1)

Long waves generated by ships moving in shallow water

Farnaz Rezvany Hesary

Master’s Thesis, Spring 2019

(2)

of Mathematics, University of Oslo. The scope of the thesis is 60 credits.

The front page depicts a section of the root system of the exceptional

Lie groupE8, projected into the plane. Lie groups were invented by the Norwe- gian mathematician Sophus Lie (1842–1899) to express symmetries in differential equations and today they play a central role in various parts of mathematics.

(3)

Abstract

Motivated by the relatively new phenomenon observed in Oslofjord where ships moving with subcritical speeds across substantial depth changes generated long waves with heights up to 1.4m, the linear generation mechanism for these upstream waves is inves- tigated using the linear shallow water equations. The simulations are performed both with one and two horizontal dimensions where the average depth at the location where the bottom variation happens is twice the change in depth, ∆h/¯h= 0.5. Analytical cal- culations on the amplitude of the generated free waves as the source moves over a step in bottom topography shows good agreement with the numerical results. For ships moving from deep to shallow water, the maximum elevation of the generated waves grow accroding toF rⁿ, where n is in the range 3.6−4.6. The simulations with two horizontal dimensions are compared to a dispersive model to give a measure of how well the linear shallow water model captures this phenomenon.

iii

(4)

(5)

Acknowledgements

I wish to express my very deep gratitude and appreciation to my two supervisors, Geir K. Pedersen and John Grue, for their patient guidance, valuable and constructive sug- gestions, enthusiastic encouragement and useful critiques of this thesis. I wish to show special gratitude towards their willingness to give of their valued time so generously.

I also wish to express my appreciation and gratefulness towards my beloved husband, Faraz Ferdosian, for his unwavering support, positive reinforcements and enduring pa- tience during all my years of study and specially these last months.

v

(6)

(7)

List of Figures

3.1 Bottom elevation by a step . . . 8 3.2 Pressure distribution, where x₀ = 10,α = 6,pm = 6.8/H₀ and Lp =

224/H0. Pressure parameters are chosen according to (Col) . . . 10 3.3 A sketch of the geometry. (Constant depth) . . . 11 3.4 Staggered grid for the solution of equations: (3.42) and (3.43) (Pedersen

[2015]) . . . 12 3.5 Convergence of the numerical solution to the analytical one by grid refine-

ment. There have been used fewer markers than datapoints to make the figure more readable. . . 13 3.6 A sketch of the geometry. (Variable depth with linear slope) . . . 14 3.7 Wave propagation at constant depth . . . 15 4.1 Comparison of the analytical and numerical solutions for the evolution of

the forced wave moving over a step . . . 18 4.2 Evolution of the forced wave moving over a step . . . 19 4.3 The transformation of the forced wave and generation of free waves,

∆F r= 0.19 . . . 22 4.4 The transformation of the forced wave and generation of free waves,

∆F r= 0.08 . . . 24 4.5 The slopem = 0.1 and the depth before the slope h₁ = 1 is kept constant,

while the slope length and h₂ vary. ∆F r = 0.19,0.13 and 0.08 for h₂ =

0.6,0.7 and 0.8 respectively. t=200. . . 25 4.6 The effect of slope length on the maximum elevation and the wavelength

os the free upstream wave. ∆F r= 0.19,m = 0.1 (solid line),0.02 (dashed line) and 0.01 (dotted line) . . . 26 4.7 Stationary response due to the applied steady surface pressure moving

with the velocity U . . . 27 4.8 h1 = 1, h₂ = 0.6⇒¯h= 0.8, ^∆h_¯

h = 0.5 . . . 27 5.1 Pressure distribution, where x₀ = y₀ = 0, α = 6, p_m = 6.8/H₀,L_p =

224/H0 andW_p = 35/H0. Pressure parameters are chosen according to

(Col) . . . 30 5.2 Computational domain for the solution of (5.4) at t=0, shows the order-

ing of the unknowns . . . 32 5.3 Comparison of the solution with different tolerances that the solver needs

to achieve . . . 32 5.4 Verification of the solver by MMS . . . 33 5.5 Stationary response, ¯η at t=0 . . . 34

ix

(10)

5.6 Arakawa C-grid for space differences . . . 36 5.7 Initial velocities . . . 36 5.8 The solution along a transect at y = 0 There are a fewer markers than

there are data-points, to make the plots easier to study . . . 37 5.9 . . . 38 6.1 Free upstream wave, Fr in a subcritical range of 0.65−0.84, Ly = 40 . . . 40 6.2 Surface elevation in front of the pressure field moving in an horizontally

unbounded fluid . . . 40 6.3 A comparison of the amplitudes and wavelengths of the upstream free

waves generated by ships moving with the velocity U = 0.65, in channels of different widths. . . 41 6.4 Upstream wave elevation generated by a disturbance moving at subcriti-

cal velocities,L_y = 100 . . . 42 6.5 Maximum elevation’s growth by increasing Fr,L_y = 100 . . . 42 6.6 . . . 43 6.7 Surface elevation of the upstream waves when the ship moves from deep

to shallow water. With reference to the width of the fjord at Huk,L_y =

40 and the ship’s velocity is set to 0.6 corresponding to 18.3 knots. . . 43 6.8 η_max and η_min of the upstream waves when the ship moves from deep to

shallow water and then again to deeper water, U = 0.6 . . . 44 6.9 Surface elevation of the upstream waves along a transect at y= 0 . . . 45 6.10 The effect of the ship’s velocity on maximum elevation and length of the

upstream waves . . . 45 6.11 The initial response for ships with different displaced volume . . . 46 6.12 Upstream waves caused by ships of different displaced volume. . . 47 7.1 Generated waves both ahead and behind the pressure disturbance as it

moves from deep to shallow water . . . 49 7.2 Time evolution of the elevation in front of the pressure disturbance mov-

ing with velocity U = 0.6 . . . 50 7.3 Comparison of the elevation for the upstream wave, generated by a dis-

turbance moving at different sub-critical velocities . . . 50 7.4 The growth of maximum elevation by increasing Froude number. The

pressure disturbance moves from deep to shallow water . . . 51 7.5 The growth in maximum height of the upstream wave by an increasing

Froude number. The disturbance has passed another depth change from

shallow to deep water . . . 51

(11)

List of Tables

4.1 Simulation parameters . . . 17 4.2 The difference between the analytical and numerical results . . . 19 4.3 Maximum elevation, area between the curve and the x-axis and the wave-

length for the forced wave before and after the step and for the free waves 20 4.4 Simulation parameters . . . 21 4.5 Simulation parameters . . . 23 4.6 Maximum elevation and the wavelength for the three components od the

solution after the depth change . . . 25 4.7 Maximum elevation and the wavelength of the upstream free wave . . . 26 4.8 The effect of pressure field’s velocity on the maximum elevation and the

wavelength of the upstream free wave . . . 28 5.1 Simulation parameters . . . 34 6.1 Maximum elevation and the wavelength of the upstream free waves in

shallow channels of different widths. . . 41 6.2 Maximum elevation and the wavelength (λ = V /ηmax), of the upstream

free waves generated by a ship as it moves from deep to shallow water . . . 46

xi

(12)

(13)

Chapter 1 Introduction

The generation of very long waves with relatively high amplitudes ahead of the largest cruise ferries operating in Oslofjord has been a curious topic, both for the people who have observed it and the scientists who have tried to explain it. Color Line’s ships operating between Oslo and Kiel, with a displaced volume of almost 37000m³ are the largest cruiseferries in the fjord that can also travel with speeds up to 22 knots. The displaced volume and the speed of these ships are directly connected to the upstream long waves generated by their movement across significant depth changes at the bottom of the fjord, since these waves have not been observed for the smaller ships which also move at smaller velocities.

The height of the upstream waves has been measured to reach a maximum of 1.4m, in about 30 to 60 seconds, which can cause a significant amount of erosion when the water retreats. The bath houses and piers along the shore of the fjord has been observed to lean towards the sea which could be caused by the erosion caused by these waves. This is probably the reason for the interest of the media in this topic.

John Grue, professor at the mathematical institute of the University of Oslo who used the term mini-tsunami for these waves, was contacted to solve the mystery and try to explain the unfortunate phenomenon happening in Oslofjord.

Since this long-wave phenomenon is relatively new, there is not much coverage in the scientific literature about it, except for the work of John Grue published in the Journal of Fluid mechanics (Grue [2017]). Using a fully dispersive and linear analysis, he inves- tigates the upstream generation mechanism at positions where the bottom of the fluid domain changes significantly, such that the change in depth ∆h, is almost comparable to the average depth ¯h, at the position.

The generation of upstream waves caused by moving pressure fields at close to critical speeds has been a very interesting topic for a long time and there is a vast number of scientific litterature addressing the generation of a stationary pattern ahead of a moving disturbance. Seung-Joon Lee et al., using generalized Boussinesq equations and the Korteweg-de Vries equation found that a moving pressure field at transcritical speeds in shallow water generates a succession of solitary waves ahead of the disturbance (Lee et al. [1989]). C.C.Mei in (Mei [1986]) developed a theory for soliton radiation by slen- der bodies moving near the critical speed in an unbounded shallow channel. Tomas Torsvik in ((Torsvik et al.)) studied the generation of solitary waves with large amplitudes ahead of a pressure field. He concluded that the pattern of the generated wave

1

(14)

differs significantly for solitons caused by an increasing or a decreasing Froude number.

Geir K. Pedersen studied wave generation by moving disturbances using a Boussinesq type numerical model and could conclude that upstream radiation must be expected to occur for channels of any width(Pedersen [1988]).

All the studies mentioned above consider moving disturbances with a velocity close to the shallow water speed, where the non-linear effects are of major importance since linearized theory becomes resonant at critical Froude number. However the phenomenon observed in Oslofjord differs completely from soliton generation, as it is linear and the pressure fields move at velocities that are well below the shallow water speed.

The main objective of this thesis is to model and explain the generation of the upstream free waves at Huk in Oslofjord, where ∆h/h¯≈0.5, using the linear shallow water(LSW) equations. The ships will be modelled by a pressure field of similar length, width and draught. The results obtained by the LSW model will than be compared to a linear and dispersive model, to get a better idea of how well the LSW model can capture this phenomenon.

This thesis is organized as follows: chapter two presents the LSW equations. Chapter three gives a formulation of the two dimensional shallow water wave model where the LSW equations with one horizontal dimension and a added pressure term are presenred.

These equations will be solved first analytically for the case of constant depth and then a bottom elevation by a step is considered. Further in this chapter the equations are scaled and solved numerically. Chapter four shows the results of the numerical simulations, and the effect of different parameters on the results are discussed. Chapter 5 attends to the LSW equation with two horizontal dimensions with an additional surface pressure term. In this chapter a stationary response due to the presence of a surface pressure field is found numerically and then used to simulate the solution to the LSW equations in variable depth. The results concerning the three dimensional simulations of the developement of the stationary response in variabel depth is shown in chapter 6, where the effect of the lateral length of the domain, velocity of the pressure distribution and displaced volume of the ship, on the results are discussed. Chapter 7 gives a limited comparison of the results obtained by using a linear non-dispersive model to a linear but dispersive one, where we make some simulations with a Boussinesq solver called Globouss. Chapter 8 includes a conclusion.

(15)

Chapter 2 Shallow water equations

When we consider long waves with small amplitude propagating in shallow water, where the wavelength is large in comparison to the depth such that ^h_λ 1, and kh 1, the horizontal fluid particle velocity is much larger than the vertical one so the fluid particles move in elliptical paths and we assume that the vertical accelerations of the fluid particles have a negligible effect on the pressure in the fluid so that the pressure is almost only affected by the acceleration of gravity, therefore assumed hydrostatic, such that:

p=ρg(η−z) +p₀ ⇒ ∇_hp=ρg∇_hη (2.1)

Where p₀ is the pressure above the fluid,g is the gravitational acceleration,ρ is the fluid density and ∇h denotes the horizontal component of the gradient such that:

∇=∇_h+k ∂

∂z =i ∂

∂x +j ∂

∂y +k ∂

∂z (2.2)

So as we can see the pressure gradient is independent of the vertical coordinate z.

We also assume that the horizontal components of the acceleration and hence the velocity are also independent of z at any time t. Further we may define a horizontal volume flux integrated from the bottom to the surface of the fluid that allow us to employ a depth integrated continuity equation which together with the horizontal part of Euler’s equation of motion give us a a set of nonlinear equations for continuity and momentum that describes the propagation of long waves in shallow water. A detailed derivation of these equations is given in (Gjevik et al. [2015]).

The equations describing the motion of plane and nonlinear waves in shallow water are:

∂v_h

∂t +v_h· ∇_hv_h =−g∇_hη (2.3)

∂η

∂t =−∇_h ·((h+η)v_h) (2.4)

These are the nonlinear shallow water equations (NLSW) and as we can see these equations have one less spatial free variable compared to the standard continuity and momentum equation, which allows us to model three dimensional physical phenomena using only two dimensional mathematics.

Since our main interest in this thesis is to simulate and investigate waves that are generated by a moving pressure distribution, we need to include an additional pressure term

1

ρ∇_hpˆin the momentum equation (2.3).

3

(16)

Further assuming that both the particle velocity, v_h and the free surface elevation,η and their derivatives are all small quantities which is an assumption that yields for waves of small steepness, then the nonlinear terms will be much smaller than the linear ones, so the NLSW equations will simplify to the linear shallow water equations (LSW) which in the case of constant depth will give us the linear wave equation.The LSW equations for 3D motion are:

∂v_h

∂t =−g∇_hη− 1

ρ∇_hpˆ (2.5)

∂η

∂t =−∇_h ·(hv_h) (2.6)

The LSW equations are the simplest set of equations that describe the propagation of non-dispersive linear long waves where the horizontal scale of the fluid motion is much bigger than the depth of the fluid and therefore these equations are widely used to model tsunami propagation.

A linear model is not always a very good one, since even very small nonlinear terms can cause very large effects in time. But in this thesis we are more interested in the physical phenomenon itself so we assume that a linear model will be a good enough approximation.

(17)

Chapter 3 Two dimensional shallow water wave model

3.1 Mathematical formulation

Since we are investigating a linear long-wave phenomenon in this thesis, we will get a good representation using the LSW equations even though they give no information on the dispersive nor nonlinear effects, so the short waves are excluded. Considering this, we do not actually get any short waves as long as we are considering a 2D motion where the width of the ship is not represented.

For 2D motion, the LSW equations with a non-hydrostatic pressure term reduce to:

∂u

∂t =−g∂η

∂x − 1 ρ

∂p¯

∂x (3.1)

∂η

∂t =− ∂

∂x(hu) (3.2)

A complete derivation of the LSW equations are given by J.J. Stoker in (Stoker [1957]).

5

(18)

3.2 Analytical solution

3.2.1 Constant depth

The linearized shallow water equations with one horizontal dimension and a pressure term on constant depth are as follows:

∂η

∂t =−H∂u

∂x (3.3)

∂u

∂t =−g∂η

∂x − 1 ρ

∂p

∂x (3.4)

Where, p= ¯p(x−U t) =p_mT(x−U t), u= ¯u(x−U t),η = ¯η(x−U t) and p_m denotes the maximum amplitude for the pressure distribution function.

Substitutingx−U t with ξ results in an ODE as follows:

d¯η dξ = H

U d¯u

dξ (3.5)

d¯u dξ = 1

U(gdη¯ dξ + 1

ρ d¯p

dξ) (3.6)

Substitution of (3.5) in (3.6) and integration wrt. ξ gives:

u¯= gH

U²u¯+ 1

ρUp¯= c²

U²u¯+ 1

ρUp¯= u¯

F r² + 1

ρUp¯ (3.7)

¯ η= H

Uu¯ (3.8)

⇒ u¯= p¯

ρU(1− _{F r}¹2) (3.9)

⇒η¯= H ρU²

¯ p 1−_{F r}¹2

(3.10) where F r = ^U_c is the Froude number and c = √

gH is the shallow water wave celerity and we have assumed that ¯η,u¯ → 0 as x → ∞.

Knowing the particular solution and assuming an initial state as follows:

η(x,0) = ¯η(x,0) +η^(H)(x,0) = 0 (3.11)

u(x,0) = ¯u(x,0) +u^(H)(x,0) = 0 (3.12)

We can find a homogeneous solution to the LSW equations where p= 0.

The homogeneous solutions together with the particular solutions, will give us the complete analytical solutions to the LSW equations with an additional pressure term.

We assume:

η^(H) =F(x−ct) +G(x+ct) (3.13)

u^(H) = c H

F(x−ct)−G(x+ct) (3.14)

(3.15)

(19)

3.2. ANALYTICAL SOLUTION 7 At t= 0, we have:

η(x,0) = H ρU²

p_mT(x) 1−_{F r}¹2

+F(x) +G(x) = 0 (3.16)

u(x,0) = 1 ρU

p_mT(x) 1− _{F r}¹2

+ c H

F(x)−G(x)= 0 (3.17)

By adding (3.16) and (3.17), we can eliminate G(x), so we find an expression forF(x).

Then by subtraction of (3.17) from (3.16) and substitution of F(x) into the resulting expression, we findG(x).

We then have:

⇒F(x−ct) = (c+U)Hp_mT(x−ct)

2ρc(c²−U²) (3.18)

⇒G(x+ct) = (c−U)Hp_mT(x+ct)

2ρc(c²−U²) (3.19)

The final solutions are then:

η(x, t) = Hp_m 2ρc(c²−U²)

(c+U)T(x−ct) + (c−U)T(x+ct)−2cT(x−U t)

(3.20) u(x, t) = p_m

2ρ(c²−U²)

(c+U)T(x−ct) + (c−U)T(x+ct)−2U T(x−U t)

(3.21) The initial conditions are fulfilled since both η and u become identically zero when t= 0.

The first two terms in the solutions above represent free waves that propagate in opposite directions with the phase velocityc, also called the wave celerity. The last term is the forced solution that we get due to the presence of a pressure field that moves with the constant velocity U.

As we can see from (3.20) and (3.21), the surface elevation and the fluid particle velocity for the wave moving to the right is higher than the one moving to the left.

Further we may notice that as long asU < c such that we have a subcritical Froude number, then the forced solution is negative, but it will change sign whenF r >1.

We may also notice that these solutions are not valid for critical flow where F r = 1.

In that case we have to assume a solution of the formη = tY(x−U t), that we can substitute into:

∂²η

∂t² =c²∂²η

∂x² + c² ρg

∂²p

∂x² (3.22)

Which we obtain by first taking a temporal differentiation of (3.3) and then substituting (3.4) into the resulting expression.

Further an integration of (3.22) wrt. ξ gives:

η(x, t) = −tcp_m

2ρgT⁰(x−U t) (3.23)

As we can see when F r = 1 we get a solution that grows with time which means that the linear model is resonant and theoretically the linear response is infinite so in this case we will need to consider nonlinearity which will provide a way for energy transfer as opposed to a linear model.

(20)

3.2.2 Bottom elevation by a step

While the analytical solutions given in the previous section yield for waves propagating in shallow water over a uniform depth, we may want to investigate the solution for variable depth. Finding simple analytical solutions will give us more insight into the problem and is always a good starting point for further investigation. So let’s start by an analysis of what happens as a surface wave propagates over a step in the bottom topography.

A bottom profile is given by:

h(x) =







h₁, x <0

h₂, x >0 (3.24)

x H(x)

−h

₁

−h

₂

0 x<0 x=0 x>0

Figure 3.1: Bottom elevation by a step

Here we may choose the forced solution ¯η_I(x−U t) as the incoming wave propagating with the velocity U towards the step. After ¯η_I hits the step, there will be a wave going backward with the velocityc₁ = √

gh₁ which we may represent by η_B(x−c₁t) while we also get another forward moving wave propagating with the velocity c₂ = √

gh₂ represented by η_F(x−c₂t) and of course we will also have the forced solution after the step which we represent by ¯η_T(x−U t). To sum up we have:

x <0 : η₁ = ¯η_I(t− x

U) +η_B(t+ x

c₁) (3.25)

x >0 : η₂ = ¯η_T(t− x

U) +η_F(t− x

c₂) (3.26)

Where we have rewritten all the functions of the form Y(x±αt) to the form Y(t± _α^x) to simplify the problem. Rescaling the argument in these functions just make them easier to handle at the step wherex= 0 since the functions then do not on any velocity.

Requiring continuity at the step and according to (3.13), at x= 0 we have:

η₁ =η₂

⇒ η¯_I(t) +η_B(t) = ¯η_T(t) +η_F(t) (3.27)

− h₁p_m

ρ(c²₁ −U²)T(t) +η_B(t) = − h₂p_m

ρ(c²₂−U²)T(t) +η_F(t) (3.28)

(21)

3.2. ANALYTICAL SOLUTION 9 where we have used the forced term in the solution we found in (3.20) for ¯η_I(t) and

¯ η_T(t).

Further we use the same notation as before, but now for the horizontal fluid particle velocity functions.

x <0 : u₁ = ¯u_I(t− x

U) +u_B(t+ x

c₁) (3.29)

x >0 : u₂ = ¯u_T(t− x

U) +u_F(t− x

c₂) (3.30)

Requiring continuity of the volume flux at the step and according to (3.14), at x= 0 we have:

h₁u₁ =h₂u₂

⇒ h₁u¯_I(t) +h₁u_B(t) =h₂u¯_T(t) +h₂u_F(t) (3.31) h₁U p_m

ρ(c²₁−U²)T(t) +c1ηB(t) = h₂U p_m

ρ(c²₂−U²)T(t)−c2ηF(t) (3.32) Multiplying every term in (3.32) with −_U¹ and then adding (3.28) and (3.32) gives:

η_B(t) = 1−^c_U²

1 + ^c_U¹η_F(t) (3.33)

Plugging (3.33) into (3.28) eliminatesη_R(t) and gives:

η_F(t) = pm(U +c₁)(−h₂c²₁+h₂U²+h₁c²₂−h₁U²)

ρ(−c₁−c₂)(c²₂−U²)(c²₁−U²) T(t) (3.34)

⇒η_B(t) = p_m(U −c₂)(−h₂c²₁+h₂U²+h₁c²₂−h₁U²)

ρ(−c₁−c₂)(c²₂−U²)(c²₁−U²) T(t) (3.35) (3.36) Sincec₁ =√

gh₁ and c₂ =√

gh₂ are both bigger than one and c₁ > c₂, we can see from (U +c₁) and (U −c₂) thatηF(t) has a higher amplitude than ηB(t) and as long as the

flow is subcritical all the way such thatU < c₂, η_B(t) will be negative.

(22)

3.3 Pressure distribution

To be able to make a simulation of the two-dimensional motion of the waves modelled by (3.1-3.2) we may choose a one dimensional surface pressure function of the form:

¯

p(x−U t) =p_me

−

(

^2(x−x_L⁰^{)−U t}

p

)

^α

(3.37) where p_m is the maximum amplitude of the pressure distribution which could be inter- preted as a ship’s draught while L_p gives the horizontal length at half the maximum of the function which we can think of as the length between perpendiculars of a ship or the length of the ship at waterline. x₀ is the center of the maximum of the pressure distribution so the distribution is symmetric around x₀ and α which gives the slope of the function is also the shape parameter of the pressure distribution which we choose to set to 6 for a wide and more ship-like pressure distribution function. The pressure field moves to the right with a constant velocity U.

Figure 3.2: Pressure distribution, where x₀ = 10, α= 6, p_m = 6.8/H₀ and L_p = 224/H0.

Pressure parameters are chosen according to (Col)

(23)

3.4. SCALING 11

3.4 Scaling

x H(x)

x

₀

−H 0

L L

_p

H

₀

U

Figure 3.3: A sketch of the geometry. (Constant depth)

By applying a characteristic length scale H₀ which we can set to be the maximum depth, we can make the following non-dimensional variables:

x^∗ = x

H₀ , t^∗ = t H₀

q

gH₀, h^∗ = h H₀ η^∗ = η

H₀ , u^∗ = u

√gH₀ , p^∗ = p ρgH₀

where the non-dimensional quantities are marked by a star. The scaled non-dimensional equations are:

∂η^∗

∂t^∗ =− ∂

∂x^∗(h^∗u^∗) (3.38)

∂u^∗

∂t^∗ =−∂η^∗

∂x^∗ − ∂p^∗

∂x^∗ (3.39)

Also both lengths L_p and p_m in the pressure distribution function (3.37) will also be scaled by the reference depth H₀.

All equations in the rest of this chapter are dimensionless, but we drop the stars(*) to simplify the notation.

(24)

3.5 Numerical solution

3.5.1 Discretization of the governing equations

The governing equations given by (3.3) and (3.4) are discretized on a staggered grid shown in Figure:(3.4) as follows:

η_i−ⁿ 1 2

−η_i−ⁿ⁻¹1 2

∆t =−h_iuⁿ⁻

1 2

i −hi−1uⁿ⁻

1 2

i−1

∆x (3.40)

uⁿ⁺

1 2

i −uⁿ⁻

1 2

i

∆t =−

η_i+ⁿ 1 2

−ηⁿ_i−1 2

∆x

−

pⁿ_i+1 2

−pⁿ_i−1 2

∆x

(3.41)

⇒η_i−ⁿ 1

2 =−∆t

∆x

h_iuⁿ⁻

1 2

i −hi−1uⁿ⁻

1

i−12

+ηⁿ⁻¹_i−1 2

(3.42)

⇒uⁿ⁺

1 2

i =−∆t

∆x

ηⁿ_i+1 2

−η_i−ⁿ 1 2

+pⁿ_i+1 2

−pⁿ_i−1 2

+uⁿ⁻

1 2

i (3.43)

where n= 0,1, ..., N and i= 0,1, ..., M. ∆t and ∆x denote the temporal and the spatial spacing respectively, such that t=n∆t and x=m∆x.

x t

O

∆t

∆x

u η

Figure 3.4: Staggered grid for the solution of equations: (3.42) and (3.43) (Pedersen [2015])

As we can see in Figure:(3.4), we use centered two-point spatial differences so we store the elevation η and the velocity uin staggered locations. This way we can calculate the gradient of η accurately where u is and we can calculate the gradient of u whereη is stored.

(25)

3.5. NUMERICAL SOLUTION 13

3.5.2 Implementation

We may start the simulation by setting the forced solution as the initial state, such that we only follow this solution, and the free waves are excluded. But for now we will start the simulation from rest. Also an initial state where the time-array index ˆn = 0 and the fluid is at rest so bothη and uare zero and the pressure field is defined by

p = p

m

e

⁽^2x^L⁾⁶. We may notice that the array index ˆn differs from the index n used in the discretization such that at ˆn= 0, η is at n = 0 while u is at n= ¹₂.

After initialization, we use the discretized continuity equation (3.42) to find η in the whole domain at a new time-step given by n, and then we update the pressure pat the new time-step n∆t. We then use the updated η and pto find velocity u in the whole spatial domain.

There’s assumed a wall at both ends of the spatial domain and so the no slip boundary condition is applied there such that u[0] =u[M] = 0.

To limit dispersion caused only by discretization parameters as described in (Langtan- gen [2013]), we test the numerical model with different spatial step sizes to see if a grid refinement would affect the solution, see Figure:(3.5).

Figure 3.5: Convergence of the numerical solution to the analytical one by grid refinement. There have been used fewer markers than datapoints to make the figure more readable.

Figure:(3.5) shows that for a coarse grid where dx = 1, we can clearly see numerical dispersion. A grid refinement by setting dx= 0.5, reduces dispersion but we can still see

(26)

some points ending up well outside the analytical solution. A further refinement of the grid by choosingdx = 0.2 makes the numerical solution coincide with the analytical one almost completely except atx≈155, where we can see that the numerical solution falls a bit under the analytical one. So all simulations of the two dimensional fluid motion in this thesis is done by choosing a fine grid where dx= 0.1, where we can see that the two solutions coincide completely.

Further to make sure that the solution is stable, we choose ∆t according to the CFL criterion as described in (Langtangen and Linge [2017]):

CF L≡ |c|∆t

∆x ≤1

Where cis the wave celerity which in the scaled form simplifies to √ H.

In case of a variable depth, we will need a function for the depth change h(x). For a shoaling case with linear slope,h(x) is set to h₁ for x < L₁ and to h₂ forx > L₁+L_s and a linear slope from h₁ to h₂ is given by a linear functiony =mx+b where the slope m is given by ^∆H_L

s . All the parameters needed to run the algorithm for a shoaling case on linear slope are shown below in Figure:(3.6).

x H(x)

x

₀

h

₂

h

₁

0 L

₁

L

_s

L

₂

L L

_p

x

₁

x

₂

U

Figure 3.6: A sketch of the geometry. (Variable depth with linear slope)

3.5.3 Verification

In order to verificate the numerical method for solving the linear shallow water equations on constant depth and our implementation of it, we can adjust the depth function given in Section 3.5.2 to yield for a constant depth by giving the same value to h₁ and h₂ and then we can compare the numerical results with the analytical solution derived

(27)

3.5. NUMERICAL SOLUTION 15 in Section 3.2. We find the three components to the solutionT(x−U t), T(x−ct) and T(x+ct) from the pressure distribution function defined in (3.37).

The center of the pressure distribution is at x₀ = 100. All parameters necessary to make the simulation are scaled with a reference depth H₀ = 25 according to Section 3.4.

The pressure parameters pm, Lp and α discussed in Section 3.3 are set to 0.27, 9.00 and 6.00 respectively. The pressure velocity U is set to 0.65 which corresponds to a Froude number of equal value whenH = 1.

Figure:(3.7) below shows the results for running the simulation with the parameters chosen as described above and we can see that the two solutions coincide completely.

(a) Initial state (b)

(c) (d)

Figure 3.7: Wave propagation at constant depth

(28)

We can clearly see the generation and the propagation of the three components of the solutions (3.20) and (3.21) from the figures above. We have the positive free wave solutions moving with the shallow water wave celerity in opposite directions, where the one moving to the right has a higher amplitude then the one moving to the left. We can also see the negative forced solution, which is generated as a response due to the presence of a moving surface pressure field.

The velocity of the free waves are higher than the forced one c > U so the flow is subcritical. We can see from Figure:(3.7c) and Figure:(3.7d) that the positive wave travel- ling to the right which propagates with a higher velocity than the negative forced wave gets further and further apart from the forced solution.

We may also notice that all three components of the solution have the exact same length. This is in agreement with (3.20) and (3.21) where we can see that all three components of the solution have the same length scale.

(29)

Chapter 4 Results and discussion: Two dimensional case

4.1 Bottom elevation by a step

We have found an analytical solution for the case where the stationary response comes to a step-type bottom topography. We now want to investigate this case numerically by running the simulation using the following parameters.

H1=1 pm=0.27 H₂=0.6 U=0.65

L₁=40 α=6 Ls=dx=0.1 Lp=9

Table 4.1: Simulation parameters And a bottom profile as in Figure:(3.1) is given by:

h(x) =







h₁ = 1 , x <40

h₂ = 0.6 , x >40 (4.1)

We set the forced solution as initial function for the surface elevation, the initial free waves will no longer exist and we will only follow the forced wave so that we can see what happens to the stationary response as it comes to a sudden change in depth.

We should also remark that the analytical solution for the stationary response that we found in Section 3.2.1 yield both before and after the step since the depth is constant at both sides of the step, and the analytical solution for the forward and backward propagating free waves given in Section 3.2.2 will occur instantly after the step, while the numerical solution will take some time to adjust itself to the change in depth.

17

(30)

(a) Transformation of the analytical solution for the stationary response before and after the step

(b) Numerical and analytical solutions after the step

Figure 4.1: Comparison of the analytical and numerical solutions for the evolution of the forced wave moving over a step

Figure:(4.1a) shows the analytical solution for the stationary response and we can see that right after the step, where the depth decreases, the response immediately gets deeper.

Figure:(4.1b) shows the analytical solutions we found for the forward and backward propagating free waves, in addition to the numerical and the analytical solution for the evolution of the forced wave. As we can see here because the change in depth happens so abruptly, the generation of the free waves happens right after the step and then it takes a while for all three components of the solution to get adjusted to this change and become steady. In the next figure we can see how the numerical solution evolves

(31)

4.1. BOTTOM ELEVATION BY A STEP 19 and eventually becomes almost equal to the analytical solutions we have found for the stationary response and the forward and backward free waves. Table:(4.2) shows there is a small difference between the numerical and analytical results.

Analytical solution Numerical solution |analytical sol. - numerical sol.|

¯

ηT,max 0.919437 0.919443 6.00E−6

η_F,max 0.416955 0.416888 6.70E−5

η_B,max 0.031486 0.031820 3.34E−4

Table 4.2: The difference between the analytical and numerical results

(a)

(b)

Figure 4.2: Evolution of the forced wave moving over a step

(32)

Figure:(4.2) shows the evolution of the forced wave coming to a step. As we discussed earlier and we can see in the figure above, the generation of the free waves happens after the step. As the response grows deeper, an upstream free positive wave and a negative one with small amplitude moving in the other direction are generated. By taking the integral A=^R_x^x₀^M ηdx for each of the functions for the surface elevation of the response before and after the step and also for the forward and backward wave, we find:

Z xM

x0

¯

η_T,t=200−

Z xM

x0

¯ η_I,t=0 =

Z xM

x0

η_F,t=200−

Z xM

x0

η_B,t=200 (4.2)

The wavelength is defined by:

λ= 1 a

Z ∞

−∞ηdx (4.3)

where a=|η_max|, (Yamashita et al. [2016])

Table:(4.3) shows maximum elevation and wavelength for all four components of the solution.

|η_max| Area λ

¯

η_I 0.471 3.93 8.35

¯

η_T 0.919 7.68 8.35 η_F 0.417 4.15 9.97 η_B 0.032 0.40 12.71

Table 4.3: Maximum elevation, area between the curve and the x-axis and the wavelength for the forced wave before and after the step and for the free waves

4.2 Shoaling on linear slope

We now want to investigate what happens when the initial response comes to a linear slope and the depth decreases. A sketch of the geometry is given in Figure:(3.6). Fur- ther we will compare the results to what we had for constant depth. Since we know that a change in depth will cause the Froude number to vary, F r= ^U_c

=

^√^U

h, we will make the simulation for some different depth changes to see how ∆F r = F r_h₂ −F r_h₁ affects the evolution and the transformation of the forced solution.

4.2.1 Case 1: ∆ Fr = 0.19

Scaled with the maximum depth H₀ = 25m, the length of the fluid domain is L = 200.

The pressure field is symmetric around x₀ = 10. A linear slope starts at x = 40 such that atx = 44 the depth is decreased from h₁ = 1 toh₂ = 0.6 which means that the Froude number increases from 0.65 to 0.84, so the flow is subcritical all the way.

The forced solution is chosen as the initial condition as shown in Figure:(4.3a), so we follow only the forced wave which is a response to the presence of the moving pressure field and hence we exclude the free waves until x = 40, where the slope stars. This way we are able to study the transformation of the forced wave alone on the slope. Ta- ble:(4.4) gives an overview of the chosen values for the parameters necessary to make

(33)

4.2. SHOALING ON LINEAR SLOPE 21 the simulation. Section 3.3 and Figure:(3.6) explain each parameter involved in the simulation.

H₁=1 p_m=0.27 H₂=0.6 U=0.65

L₁=40 α=6 L_s=4 L_p=9

Table 4.4: Simulation parameters

(a)

(b)

(34)

(c)

Figure 4.3: The transformation of the forced wave and generation of free waves,

∆F r= 0.19

We can see that the generation of free waves starts where the slope starts. At time t= 58, we can see that the moving pressure field has passed the slope so the bottom is again flat, but the forced wave is still adapting to the depth change by an increase in the amplitude of the free wave while the forced wave is getting deeper.

The forced wave continues moving with the velocity U but the free wave that is generated with two components that propagate in opposite directions move with the shallow water celerity c_h1 andc_h2 and in agreement with the analytical solution we found for constant depth in Section (3.2.2) the negative backward going free wave has a much lower amplitude than the forward going one.

The results are quite similar to the case with the step, so L_s = 4. is not long enough to have any significant effect on the results. The is a 0.20% increase in maximum elevation and a 0.24% decrease in the wavelength. Figure:(4.3c) shows the last stages of the transformation of the forced wave and both the free wave and the forced wave come to a steady state. The free wave has then the same shape as the forced wave, but it is longer and moves faster so the free upstream wave will get further and further away from the forced wave behind it.

We can also see that the maximum elevation of the free forward moving wave, is connected to the change in the forced wave before and after the slope, see (4.2) in the previous section. This is only logical since the water has to go somewhere asη grows deeper and as result of this, the water gets passed on as this free forward moving wave.

To sum up our observations on Figure:(4.3), When the depth decreases, the pressure field has the same length and strength as before but the Froude number increases, resulting in a response that goes deeper. Since the bottom doesn’t allow for the water to pass through and the water has to go somewhere, it passes on as a free positive wave

(35)

4.2. SHOALING ON LINEAR SLOPE 23 going forward in the positive x-direction and a smaller negative wave going backward.

We would have a similar result if the pressure field was to accelerate, which would again result in a increasing Froude number.

4.2.2 Case 2: ∆ Fr = 0.08

To get a better grasp of the effect of increasing Froude number on the generation of the free waves after the forced wave reaches the slope, we can make the simulation using another value forh₂ such that ∆F r is different from the first case.

This time we choose h₂ closer toh₁ which will result in a less increase of the Froude number. For the comparison to be fair, the slope length L_s is set to 2.0 to make sure that the slope m, remains the same as before. All simulation parameters are given in Table:(4.5).

H₁=1 p_m=0.4 H₂=0.8 U=0.65 L₁=40 α=6

L_s=2 L_p=8

Table 4.5: Simulation parameters

Choosingh₁ and h₂ as shown in Table:(4.5) results in an increase in the Froude number from 0.65 to 0.73 which results in a ∆F r= 0.08 which is less than the first case. To be able to compare the results, we observe the solutions at the same timestep as what we had in the first case. Figure:(4.4) shows the solution at some of the same timesteps as in the previous case.

(a)

(36)

(b)

(c)

Figure 4.4: The transformation of the forced wave and generation of free waves,

∆F r= 0.08

We can see the same phenomenon here as we did before. Decreasing depth results in a deeper response which gets passed on as a free wave moving forward and a smaller negative wave moving backward. But in this case, the difference in the surface elevation η before and after the depth change is much smaller than the first case since ∆h is smaller, so the resulting free waves are much smaller in amplitude but a bit longer in wavelength.

(37)

4.2. SHOALING ON LINEAR SLOPE 25

|η_max| λ

¯

η_T 0.5764 8.35 η_F 0.092 11.55 η_B 0.014 12.72

Table 4.6: Maximum elevation and the wavelength for the three components od the solution after the depth change

4.2.3 The effect of the depth change ∆h and the slope length L

_s

on the upstream free wave

In the previous section we saw that a small change in the length of the slope did not really affect the length of the upstream free wave, but a smaller ∆h result in a longer wave with a smaller amplitude. Figure:(4.5) shows the solution for case1, case2 and also for ∆F r = 0.13, where the length of the slope is set to 4,3 and 2 respectively, so the slope m= 0.1 remains the same.

Figure 4.5: The slope m = 0.1 and the depth before the slopeh1 = 1 is kept constant, while the slope length and h₂ vary.

∆F r= 0.19,0.13 and 0.08 for h₂ = 0.6,0.7 and 0.8 respectively. t=200.

It’s only logical to think that as long as we have not reached the maximum change from h₁ to h₂, the forced wave will continue to grow deeper and contribute to the generation of the upstream wave. Figure:(4.6) show the results of running the simulation for different values for the slope-length L_s while keeping all other parameters unaltered from case1 where h₁ = 1 andh₂ = 0.6 giving ∆F r= 0.19.

(38)

(a) t=116

(b) t=240

Figure 4.6: The effect of slope length on the maximum elevation and the wavelength os the free upstream wave. ∆F r = 0.19, m = 0.1 (solid line),0.02 (dashed line) and 0.01 (dotted line)

L_s = 4 L_s= 20 L_s = 40

|η_max| 0.4170 0.4132 0.3676 λ 9.9521 10.0153 11.2367

Table 4.7: Maximum elevation and the wavelength of the upstream free wave

As we can see the steeper the slope, the higher maximum elevation.

(39)

4.3. INVESTIGATION OF THE EFFECT OF PRESSURE DISTRIBUTION’S VELOCITY ON SURFACE ELEVATION27

4.3 Investigation of the effect of pressure distribu- tion’s velocity on surface elevation

As we can see in Figure:(4.7) and according to (3.20), in addition to a decrease in depth that results in a lower shallow water celerity, the pressure field’s velocityU will also result in a higher Froude number which will then affect the response, such that a higher velocity U gives a deeper response.

Figure 4.7: Stationary response due to the applied steady surface pressure moving with the velocity U

The higher the initial Froude number is, the higher increase we’ll have after the change in depth to a shallower region. Figure:(4.8) shows how the the velocity of the moving pressure distribution affects the maximum surface elevation of the generated upstream free wave. Here we have used the mean value for the depth to find the Froude number.

F r = ^√^U_¯

h, The method of least squares were used to find a fitted function for the data points. The generated free wave’s amplitude grows according to 1.32 F r^3.69

Figure 4.8: h₁ = 1, h₂ = 0.6⇒¯h= 0.8, ^∆h_¯_h = 0.5

(40)

Table:(4.8) shows how the maximum elevation increases with an increase in U, while the wavelength becomes shorter.

U = 0.4 U = 0.45 U = 0.5 U = 0.55 U = 0.6 U = 0.65

|η_max| 0.07 0.11 0.15 0.22 0.30 0.42

λ 16.05 14.46 12.92 11.49 10.70 9.95

Table 4.8: The effect of pressure field’s velocity on the maximum elevation and the wavelength of the upstream free wave

(41)

Chapter 5 Three dimensional shallow water wave model

5.1 Mathematical formulation

The linear shallow water equations with two horizontal dimensions and an additional pressure term which can be used to represent a ship are:

∂η

∂t =− ∂

∂x(hu)− ∂

∂y(hv) (5.1)

∂u

∂t =−g∂η

∂x − 1 ρ

∂p

∂x (5.2)

∂v

∂t =−g∂η

∂y −1 ρ

∂p

∂y (5.3)

We are going to solve these equations numerically to study the linear evolution of the response to a ship-like pressure disturbance as it moves across an appreciable depth change with a constant speed U.

Taking a temporal differentiation of (5.1) and then substituting (5.2) and (5.2) into the resulting equation, eliminates u and v and gives us the standard wave equation forη with forcing. We may then introduce a coordinate system that is fixed with the moving pressure disturbance, as we did for the 2D case, such that:

p= ¯p(ξ) = p_mT(ξ, y),u= ¯u(ξ, y), η= ¯η(ξ, y) and v = ¯v(ξ, y), where ξ =x−U t We will then have:

(U² −c²)∂²η¯

∂ξ² −c²∂²η¯

∂y² = c² ρg

∂²p¯

∂ξ² + ∂²p¯

∂y²

(5.4)

We are going to solve the equation above numerically on constant depth at t = 0, to find the initial state for the surface elevation which we will need for a further time integration to study evolution of the forced solution and the developement of waves upstream.

Further in this thesis we use only non-dimensionalized variables, unless the variables appear with a physical measurement unit. The scaling of the variables is done according to Section 3.4.

29

(42)

5.2 Pressure distribution

3D simulations of waves that are generated by a moving pressure disturbance are performed using a ship-like surface pressure field (Grue [2017]) of the following form:

¯

p(x−U t, y) =p_me

-(

^2(x−x_L⁰^{)−U t}

p

)

^α

- (

^2(y−y_W ⁰⁾

p

)

^α

(5.5) where y₀ denotes the center of the maximum pressure in the lateral direction, and W₀ gives the width of the pressure field. All other parameters are explaind in Section 3.3.

The pressure parameters L_p, W_p and p_m are connected to the ships displaced volume V_D such that V_D =L_pW_pp_m C_B, where C_B = 0.7 is the block coefficient giving the ratio of the displaced volume to the volume of a rectangular block with the same length, width and depth. (Newman [1977])

Figure 5.1: Pressure distribution, wherex₀ =y₀ = 0,α = 6,p_m = 6.8/H₀,L_p = 224/H0 and Wp = 35/H0. Pressure parameters are chosen according to (Col)

For a more accurate measure of the displaced volume we should have used

L_pp = 0.94L_p, which is a more accurate approximation for the length between perpendiculars of the ship. So the length used in this thesis to model the pressure distribution is a bit longer then what should have been used.

Long waves generated by ships moving in shallow water