Numerical Methods for the Discretization of Model Equations for Surface Water Waves: Theory and Applications

(1)

Numerical Methods for the Discretization of Model Equations for Surface Water Waves:

Theory and Applications

Philosophiae Doctor Thesis

Magnar Bjørkav˚ag

Department of Mathematics University of Bergen

Norway

June 2012

(2)

(3)

Preface

This thesis is submitted for the degree of Philosophical Doctor in applied mathematics at the Department of Mathematics, University of Bergen. The work on the thesis started in January 2007, and it is a part of the project Wavemaker funded by The Research Council of Norway through Frinat grant 171267/V30. The working environment has been the Department of Mathematics in Bergen. The main supervisor has been Prof. Henrik Kalisch, and Prof. Jarle Berntsen was co-supervisor, both at the Department of Mathematics at the University of Bergen.

The thesis consists of two parts. General background information, motivation and a summary of the papers are found in Part I of the thesis. Part II consists of 5 papers written during the work with the thesis. The titles, authors and publishers of these papers are:

Paper A: Exponential convergence of a spectral projection of the KdV equation. M. Bjørkav˚ag and H. Kalisch. Published in Physics Letters A, 2007.

Paper B: Radius of analyticity and exponential convergence for spectral projections of the generalized KdV equation. M. Bjørkav˚ag and H. Kalisch.

Published in Communications in Nonlinear Science and Numerical Simulation, 2010.

Paper C:Energy budget in a dispersive model for undular bores. H. Kalisch and M. Bjørkav˚ag. Published in Proceedings of the Estonian Academy of Science, 2010.

Paper D: Wave breaking in Boussinesq models for undular bores. M.

Bjørkav˚ag and H. Kalisch. Published in Physics Letters A, 2011.

Paper E:Benchmarking of numerical discretizations of a Boussinesq system M. Bjørkav˚ag and H. Kalisch. To be submitted.

The papers have been written in collaboration with Prof. Henrik Kalisch.

(4)

couragement during the work with this thesis. I am really grateful for the patience, the motivation and the enthusiasm you have provided in the course of study. And I have also appreciated some of the off-topic talks and discussions we have shared these couple of years.

I would also like to thank the colleagues and friends at the Department of Mathematics, for the good, and informal, which suits me well, working environment at the department. But also for the fruitful discussions during lunch and the small chats when meeting during the day.

Most of all I would like to thank my family and friends. I could not have done this without your support. Especially I would like to thank you, Gunn Elisa- beth, for your love, patience and support during these years. And to my daughter, Mathilde, you are my everyday sunshine and inspiration.

In memory of my mother, Margrethe Bjørkav˚ag

Magnar Bjørkav˚ag Bergen, June, 2012.

(5)

Part I

General Background

(8)

(9)

Chapter 1 Introduction

The main focus in this thesis is the study of properties of propagating bores and the derivation and application of spectral methods. In some sense then, the thesis could be viewed as consisting of two parts; one part covering the modelling of a physical phenomenon and the second part covering theoretical work and analysis in connection with spectral methods. Both parts, however, are developed upon the theories and equations governing surface waves, which is a common denominator for the parts and papers of the thesis. Also the numerical strategies uniﬁes the two parts of the thesis.

One of the main objectives of the thesis have been on using higher-order dispersive systems in the modelling of surface wave phenomenons. The focus has been, as mentioned, on the modelling of bores, and to develop methods to gain further quantitive information on this speciﬁc phenomenon. By comparing data from bore simulations with measurements, indication on the applicability and the limitations of the underlying Boussinesq scaling regime could be revealed. Specif- ically it is the energy budget of the undular bore and the onset of breaking in these weak bores, which is studied in the thesis.

Spectral methods are powerful tools when solving differential equations modelling different physical phenomenons. Such methods are often limited by the geometry of the physical problem, but if they are applicable they usually perform better then other available discretization techniques. In this thesis, spectral methods are derived for both a Boussinesq system of equations and the unidirectional Korteweg and deVries equation. The former method is competing in a benchmark process in order to reveal the efﬁciency of the method compared with several ﬁ- nite difference methods. For the spectral method applied on the KdV equation, a convergence estimate is supplied. Moreover, based on this convergence rate, a numerical method is constructed to determine the radius of analyticity of the solution of the equation in question.

(10)

1.1 The undular bore

In this chapter, the undular bore is introduced. Some background matter on the bore, like frequency of observations and geographical occurrences, are provided.

This is found in Section 1.1. In this ﬁrst section we also discuss the energy transfers and the energy loss in bores, which is the topic of Paper C of the thesis. Then in Section 1.2 the phenomenon of breaking waves are introduced and discussed.

Breaking waves are a common and well-known phenomenon. In the thesis, the study is limited to quantifying the onset of breaking waves in bores. It is Paper D of the thesis which covers this study.

A bore is a phenomenon which arises in rivers and open-channel ﬂows. Tidal bores form typical in spring tide conditions in river estuaries for which the geometry is favorable. Such river bores are thus generally tidal phenomenons. Among the impacts of the bore, is a signiﬁcant sediment transport in the rivers. The bore causes bed erosion beneath the bore front, and the suspended matter are carried away by the wave motion. As the river bore can travel far, the effect of the bore may be felt along considerable distances upstream of the river estuaries. So in effect the bore carries the tidal energy far upstream in those rivers in which the bore occurs.

Tidal bores are a well known phenomenon, and are observed in many rivers around the world. Among the most famous are the Severn river bore in England, the bore on the Qiantang river in China and the bore on the Dordogne river in France. The estuary of the Severn river has the second largest tidal range, 14,5 metres, in the world. The rising water is funnelled up the Severn river during the highest tides. The largest bores occur during the spring tide, but smaller bores can be seen throughout the year. The largest bore, up to 9 metres high and traveling with a speed up to 11 metres per second, is the one on Qiantang river. The usual height of this bore is between 1.5 - 4.6 metres. Other rivers exhibiting bores are Batang Lupar (benak) in Malaysia, Styx River and Daly River in Australia, River Dee and River Trent in the UK, the Savannah River in US and Petitcodiak River in Canada, to name a few.

A bore can be described as a transition between two uniform free-surface flows with different flow depths. It may, without loss of generality, be assumed that one of these flow depths are the undisturbed water level, denoted byh0. See Figure 1.1. The incident water level is then given byh₀+a₀, and we denote byα=a₀/h₀ the strength of the bore. The bore, and more generally long-wave phenomena, have interested researchers for many decades. Among the early developments on waves are the work of Airy [1] in 1841 and Stokes [49] in 1847. Also Boussinesq [14, 15] and Lord Rayleigh [47, 48] contributed to the developments. Korteweg and deVries derived in their paper [39] the famous KdV equation in 1895. Later Lord Rayleigh in [48], Benjamin and Lighthill in [5] and Peregrine [46] all pub-

(11)

1.1 The undular bore 5

bottom h₀

a₀

Figure 1.1: A schematic picture of the development of an undular bore. The bore attains its maximum amplitude at the bore front, and trailing waves appears behind the front.

lished papers more directly related on the bore phenomenon. Early experimental studies include the study of Favre [31] and Binnie and Orkney in [6]. More recent experimental works could be found in [23, 44, 22], and the survey [21].

Bores are often observed in nature, but there are few ﬁeld measurements available, and laboratory experiments are important in order to gather any quantitative information on the phenomenon. Among the classical experimental investigations, are the study of Favre [31]. In his study, Favre focused on bores propagating into undisturbed water. He found that, depending on the strength of the bore,α=a₀/h₀, the bore fall into one of three categories. A purely undular bore, which is a bore for which none of the trailing oscillations are breaking, was pro- duced when α was below 0.28. With a ratio ofα between 0.28 and 0.75, there where still oscillations behind the bore front, but one or more of the waves were observed breaking. Forα>0.75, the bore was completely turbulent. The onset of breaking in the bore, which happened at the sharply deﬁned ratio ofα=0.28 are under review in Paper D of the thesis.

In Paper C, the focus is on the energy budget for undular bores. Lord Rayleigh [48] used, in his classical theory on bores, conservation of mass and momentum in the shallow-water equations, to predict an energy loss at the bore front. Lemoine tried to match the energy loss due to the shallow-water approximation with the wave energy contained in the trailing oscillations [41], but obtained only moderate agreement. Instead of using harmonic waves to match the oscillations, Benjamin and Lighthill [5] resorted to cnoidal functions, which are traveling-wave solutions of the Korteweg-de Vries (KdV) equation. They concluded that energy is indeed disseminated by the trailing oscillations, but an additional energy loss is required, and they argued that viscosity should be taken into account. Also Sturtevant [50]

and Byatt-Smith [18] concluded that some of the energy at the bore front are absorbed by viscous action, rather then being disseminated by oscillations. And surely viscosity will play a role in a physical channel ﬂow, but viscosity are not

(12)

modeled in the shallow-water equations, as they are based on the Euler equations and not the Navier-Stokes equations. Therefore the energy loss predicted by them should be explained without resorting to viscosity. In Paper C a numerical study is conducted to give quantitative evidence to the common view that the energy loss in weak bores are in fact disseminated by the trailing oscillations.

1.2 Breaking waves

Waves in general, and breaking waves in particular, are of importance in near- coastal areas. Predicting wave heights and the occurrences of breaking waves along the shore-line is important when planning structures and environmental is- sues off and near the coastline. When waves propagate from deep-water towards the shoreline, the waves will enter the surf zone. The waves break in the surf zone, and the water continues to run up the beach (the swash). In order to develop or extend models which covers the surf zone and the swash, it is important to be able to describe breaking waves.

When waves approach a sloping beach they will start to feel the effect of the bottom. When they enter shallow water their speed will diminish while they will grow in height. This process is called shoaling, and the shoaling waves eventually break in some form. Some waves leads to so-called spilling breakers, while others end up as the more violent plunging breakers. For the spilling breaker the wave will steepen until the crest becomes unstable, and turbulent white water are spilling down the face of the wave. See Figure 1.2 and page 114 in [27]. The plunging breaker on the other hand, steepens until the front of the wave becomes vertical and the wave curls over it’s top and into the trough in front of it. Waves can also break on open bodies of water. And in fact there is a vast literature concerning the development of breaking waves in a number of situations, including shoaling, wave breaking in open bodies of water, and breaking induced by a wavemaker. Many references may be found in the review [30] and textbooks on the subject, such as [52].

In Paper D the focus is on the onset of breaking. The study can be viewed as describing the very early stage of the so-called spilling breaker (Figure 1.2). It is shown in Paper D, that it is possible to formulate a breaking criterion, which, when applied to the bore, quantiﬁes the onset of breaking reasonable well, when compared to the experimental α =0.28 found by Favre. It should be empha- sized that the aim of the study is not to describe breaking waves, however, but to quantify the onset of breaking. The model is not valid after the ﬁrst wave has started to break. When describing breaking waves there are many physical phenomenons which come into play, which our equations do not model. Among these are viscosity, which may prevent the waves from breaking at all, as was shown by

(13)

1.2 Breaking waves 7

bottom

1 2 3 4 5

Figure 1.2: A schematic picture of the spilling breaker. We will be interested in the onset of breaking in the very early stages (1 & 2) of the breaking process.

numerical approximation in [18]. Capillary effects and shear and transverse ﬂows are also disregarded. The directional spreading of the two-dimensional wave ﬁeld is important in the study of breaking waves in general circumstances [17, 45].

Near the onset of breaking, capillary and shear effects also become important [30, 42, 43]. It is possible to derive systems of Boussinesq type both for gravity- capillary waves, and for waves on shear ﬂows [26, 54]. However, in breaking waves, these effects are localized near the crest of the breaking wave, and it appears that we obtain fair agreement with the experiments in [31] without including those additional physical effects.

In this thesis, Boussinesq equations have been used to model the bore and thereby gain quantitative information on phenomenons like the energy budget in bores and the onset of breaking waves in bores. In Chapter 2 of the thesis, we give additional background and theoretical aspects on the equations utilized in the investigations conducted in the thesis. The Korteweg-deVries equation are introduced in Chapter 3. Then in Chapter 4, the numerical methods used in the thesis are presented. This chapter could also be regarded as a summary of the papers A, B and E of the thesis. Finally, in Chapter 5, the quantitative study of the bore, with respect to energy transfers and the onset of breaking, are presented.

This chapter contains a summary of the papers C and D of the thesis.

(14)

(15)

Chapter 2 Modelling of surface waves

In this chapter equations governing free-surface flow are introduced. In this derivation it is assumed that the the fluid is incompressible and inviscid. Further- more, it is supposed that the fluid-flow is irrotational and two-dimensional, and that the free surface can be described by a single-valued function η(x,t). With these assumptions, it is customary to introduce a velocity potentialφ(x,z,t). In the first section below, it is shown that on the background just mentioned, the surface-wave problem are given in terms of the Laplace equation forφ, along with a set of boundary conditions on the free surface and at the bottom.

The resulting equations from Section 2.1 are difﬁcult to solve, due to the nonlinear boundary conditions at the free surface, and the fact that these conditions should be evaluated on an a-priori unknown surface. Simpliﬁed sets of equations, which are supposed to be approximations of the full water-wave problem found in Section 2.1, have therefore been derived. Among these are the Boussinesq equations, which play a major role in this thesis. The derivation of these are found in Section 2.4. Before this, a linear solution of the surface-wave problem are given in Section 2.2. In Section 2.2, we also focus on the dispersion and the modelling of pressure, as these quantities have some relevance when describing the differ- ences between the different approximate models. In Section 2.3 the shallow-water equations are derived.

2.1 Equations

For an incompressible ﬂuid, the continuity equation is given as

∇1·�u=0. (2.1)

where ∇1 = (∂x,∂y,∂z) and�u= (u,v,w) is the fluid velocity. This equation ex- presses conservation of mass within the fluid. The motion of the fluid is governed

(16)

z=−h(x)

g

z=η(x,t) x z

❄

✲

✻

❆❆❯

Figure 2.1: Frame of reference. The description of ﬂuid ﬂow is restricted to two spatial dimensions. Still, phenomenons like long-crested waves, and waves in a canal, for which the variations in the y-direction (into the paper) are small, are well accounted for by such models.

by the Euler equations. If the densityρis considered constant, the Euler equation could be written

�u_t+ (�u·∇₁)�u=−∇p1

ρ −g∇₁z. (2.2)

Here p₁ denotes the pressure and g is the acceleration due to gravity. Equation (2.2) is basically Newton’s second law of motion. In addition to the equations (2.1) and (2.2) conditions on the boundary must be given. If the free surface is denoted by η(x,y,t) and the bottom is found at z =−h(x,y), the boundary conditions are given as

η_t+uη_x+vη_y−w=0, at z=η(x,z,t), p₁=p_a, at z=η(x,z,t), w+uh_x+vh_y=0, at z=−h(x,y),

(2.3)

where p_ais the constant atmospheric pressure.

In this thesis, and in the following, the description of the fluid flow will be restricted to two spatial dimensions. The vertical coordinate will be labeledz, while xwill be used in the horizontal direction. Phenomenons, like long-crested waves and flow in a canal, for which the variations in the second horizontal direction are small, would be well accounted for in such a set-up. See Figure 2.1, where the y-direction points into the paper. Also, the depth is taken to be constant in this thesis. Furthermore, as it is assumed that the flow is irrotational, the fluid velocity field can be represented by a velocity potential,φ, as

�u(x,z,t) =∇φ(x,z,t). (2.4)

(17)

2.2 Linear solutions 11

This leads to a reduction in the number of unknowns, and in terms of the velocity potential the equations (2.1), (2.2) and (2.3) could be expressed as

φxx+φzz=0, at −h₀<z<η(x,t) ηt+φxηx−φz=0, at z=η(x,t)

φt+¹₂φ²_x+¹₂φ²_z+gη=0, at z=η(x,t) φz=0, at z=−h₀,

(2.5)

where h₀ is the constant depth. The Euler equation and the dynamic boundary condition for the pressure have been combined to form the Bernoulli equation which is the third one from the top in (2.5).

The system of equation (2.5) is challenging to solve. The boundary conditions are non-linear, and should be evaluated on the unknown surfacez=η(x,t). There- fore simpliﬁed sets of equations have been developed. Next, linear solutions of the equations (2.5) are presented, before the simpliﬁed models of (2.5) are derived.

2.2 Linear solutions

Analysis of linear models and their solutions, give valuable physical insight in various problems. Linearization of the system (2.5) is done by assuming thatφ,η and their derivatives are small, so that products of these quantities are negligible compared with the remaining terms in the equations.

For the free surface problem, in order to apply the boundary conditions, an assumption of the form of η(x,t) is needed [40]. Taking η to have sinusoidal form, in turn dictates the form of the velocity potential. Assuming that

η(x,t) =a₀cos(kx−ωt),

φ(x,z,t) =φ₀(z)sin(kx−ωt), (2.6) where k and ω is the wavenumber and frequency of the wave, respectively, the linear versions of equations (2.5) gives the vertical dependenceφ₀(z)as

φ₀(z) = ωa₀ k

cosh(k(z+h₀))

sinh(kh0) . (2.7)

In order for this solution to be consistent with the assumptions inherent in the linearization process, the ratioa₀/h₀must be small [33].

For shallow water, the ratio between depth and wavelength should be small.

This is equivalent to the restriction that kh0 �1 in shallow water. For linear shallow water waves, this implies, from (2.6) and (2.7), that the horizontal and

(18)

vertical ﬂuid velocity to a good approximation are given as u(x,z,t) = ωa₀

kh₀ cos(kx−ωt), w(x,z,t) =ωa₀�

1+ z h₀

�sin(kx−ωt), (2.8)

respectively. In (2.8) it is seen that the horizontal velocity is independent of the vertical coordinate, and that the ratio of vertical to horizontal velocity goes as w/u∼kh₀, and sincekh₀�1 in shallow water, the ﬂuid motion is thus shown to be mainly horizontal for linear shallow-water waves.

The relation between ω and k, is called the dispersion relation. For the linearized versions of the equations (2.5), the dispersion relation is found as

ω²=k²gh₀tanh(kh0)

kh₀ . (2.9)

To classify the accuracy of different formulations of the Boussinesq equations, which will be derived in Section 2.4, it is common to analyze how the modelling of linear dispersion compares with equation (2.9). See for instance [11]. Since shallow-water is deﬁned for small values ofkh₀, a Taylor expansion of the dispersion relation (2.9) aroundkh0=0, is of interest. This is given as

c²=ω²

k² =gh₀� 1−1

3(kh₀)²+ 2

15(kh₀)⁴�

+

O

((kh₀)⁶). (2.10) Herecis the wave phase velocity. For linear shallow water waves it is seen from (2.10) thatcapproachesc=c₀=√gh₀whenkh₀→0. This is the linear limiting long- wave speed for gravity waves.

Finally the pressure distribution beneath the waves is examined. For linear gravity waves the pressure-ﬁeld could be found from Bernoulli’s equation on the

form p₁=−ρgz−ρφt. (2.11)

The ﬁrst term on the right hand side is the hydrostatic pressure. Inserting p= p₁+ρgzin (2.11), an expression for the pressure change due to the wave-motion appears. For shallow water, it is found from equations (2.6), (2.7) and (2.9) that

p=ρgη. (2.12)

The pressure change due to the wave-motion for linear shallow-water gravity waves is seen to be independent of depth and equal to the increase in pressure because of the free-surface elevation. The pressure is thus found to be completely hydrostatic, which means that the pressure at a point is caused by the weight of the water above this point.

(19)

2.3 Shallow-water equations 13

2.3 Shallow-water equations

The shallow-water equations are derived from the equations (2.1) and (2.2) under the assumption that the pressure is hydrostatic. As is seen from (2.12) this is a good approximation as long as the water is shallow enough. With a purely hydrostatic pressure, the vertical momentum equation becomes

p=ρg(η−z). (2.13)

Using the pressure term (2.13) in the two-dimensional horizontal momentum equation gives

ut+uux+wuz=−gηx. (2.14) The right hand side in (2.14) is independent ofz. Ifuinitially was independent of z, it will continue to be so, and the last term on the left hand side equates to zero.

Even though the vertical acceleration could be neglected compared with the remaining terms in the vertical momentum equation, it must be kept in the equation of continuity [56]. Integrating this vertically from the bottom atz=−h₀ to the surface gives

0=

� _η

−h₀(u_x+w_z)dz

= (η+h₀)u_x+w|z=η−w|z=−h₀

= (η+h0)ux+ηt+uηx,

(2.15)

after applying the boundary conditions at the bottom and at the free surface, and sinceuis taken to be independent ofzin this approximation. The equations (2.14) and (2.15) are the shallow-water equations, which are often written as

ηt+ ((η+h₀)u)_x=0,

u_t+gηx+uu_x=0. (2.16)

In Paper C, the shallow-water equations are used in connection with energy calcu- lations beneath a propagating undular bore. The energy loss in the bore predicted by the shallow-water equations, are compared with the predictions of the energy transfers given by a Boussinesq system of equations.

Given the nonlinear hyperbolic character of equation (2.16), all solutions of this set of equations will develop discontinuities and eventually break. As the surface steepens, the assumption that the vertical acceleration can be ignored compared with the other terms in the vertical momentum equation is no longer valid [46], and corrections to the hydrostatic pressure ﬁeld is needed. In the next section a model which incorporates such pressure correction terms is derived.

(20)

2.4 Boussinesq equations

In this section a class of equations, which fits in between the linear approximation found in Section 2.2 and the shallow-water system (2.16), are presented. The linear solutions do not include any finite amplitude or non-linear effects, while the shallow-water system assumes a hydrostatic pressure and so do not incorporate any dispersive effects. The Boussinesq equations on the other hand, models both non-linear and dispersive effects. The Boussinesq scaling regime was introduced in [15], and different types of Boussinesq equations have been used in the study of fluids. For instance, certain aspects of an undular bore have been investigated in [46].

The Boussinesq equations derived here, are the one which have been used in the Papers C, D and E of the thesis. The equations are given as,

ηt+h₀u^θ_x+ (ηu^θ)_x−h²₀ 2

�θ²−1 3

�ηxxt =0,

u_t^θ+gη_x+u^θu^θ_x−h²₀ 2

�1−θ²�

u^θ_xxt =0.

(2.17)

In the system (2.17), u^θ(x,t) represents the horizontal ﬂuid velocity at a height 0<θh0<h₀, andη(x,t)describes the surface displacement from the rest position.

As above, h₀ is the undisturbed depth of water and g denotes the gravitational acceleration. The system (2.17) represents one of the models put forward in [11].

The derivation of the system (2.17) as a model for surface waves is now re- called. This derivation is well known, but it is included here as details from this derivation are used in papers C and D. The system (2.17) is derived from the surface-wave problem (2.5) using a formal asymptotic expansion. In the derivation a polynomial expression is used to approximate the vertical distribution of the ﬂow ﬁeld. The effect of this is that the vertical coordinate is eliminated from the problem as is seen below.

Non-dimensional variables are used to make explicit the different dominant length scales in the horizontal and vertical direction, and thus expedite the derivation of the dispersive system of equations. In Part I of this thesis, such non- dimensional variables are labeled with the tilde-symbol whenever there is room for confusion. Supposing that l₀ represents a dominant wavelength, a₀ denotes a typical wave amplitude, and c₀=√gh₀ is the limiting long-wave speed, the non-dimensional variables are deﬁned by

l₀x˜=x, h₀z˜=z+h₀, a₀η˜ =η, l₀

c₀t˜=t, gl₀a₀

c₀ φ˜ =φ. (2.18) Using these non-dimensional variables, the surface-wave problem (2.5) trans-

(21)

2.4 Boussinesq equations 15

forms to,

βφ˜_x˜_˜_x+φ˜_z˜_˜_z=0, at 0<z˜<1+αη(˜ x,˜ t)˜ , φ˜_z_˜=0, at ˜z=0,

η˜t˜+αφ˜x˜η˜x˜−_β¹φ˜z˜=0, at ˜z=1+αη(˜ x,˜ t)˜ , η˜+φ˜t˜+¹₂αφ˜²_x_˜+¹₂^α_βφ˜²_z_˜ =0, at ˜z=1+αη(˜ x,˜ t)˜ .

(2.19)

The non-dimensional quantitiesα=a₀/h₀andβ=h²₀/l₀²governs non-linear and dispersive effects, respectively. These parameters are used to measure the amplitude and wavelength of a typical wave in comparison to the undisturbed depth of the ﬂuid. The Boussinesq scaling regime applies whenαandβare small and of the same order of magnitude, and this assumption will be the basis for the following development.

This derivation follows the analysis given in [56]. An ansatz on the form of the velocity potential is made by assuming that ˜φcan be written as

φ(˜ x,˜ z,˜ t) =˜

∑

^∞

0 z˜ⁿf˜n(x,˜ t)˜ . (2.20) Substituting this expansion into the Laplace equation in (2.19) and applying the boundary condition on the bottom, yields

φ(˜ x,˜ z,˜ t) =˜

∑

^∞

0

(−1)^m z˜^2m (2m)!

∂^2mf˜(x,˜ t)˜

∂x˜^2m β^m, (2.21) where ˜f = f˜₀ is the non-dimensional velocity potential on the bottom. Using the expression (2.21) in the free surface conditions in (2.19) the system

η˜t˜+v˜_x_˜+α(η˜v)˜ _x_˜−¹₆βv˜_x_˜_x_˜_x_˜=

O

^�αβ,β²� ,

˜

vt˜+η˜_x_˜+αv˜v˜_x_˜−¹₂βv˜_x˜_˜x˜t =

O

^�αβ,β²�

, (2.22)

is obtained. The velocity ˜vin (2.22) is the horizontal velocity at the bottom. This velocity is thus given as the ﬁrst term in the expansion

φ˜x˜(x,˜ z,˜ t) =˜ v(˜ x,˜ t)˜ −1

2β˜z²v˜_x_˜_x_˜(x,˜ t) +˜

O

^�β²�

. (2.23)

Now, the expression (2.23) can be used to derive systems of equations for which the horizontal velocity is modeled at a non-dimensional depth θ rather then on the bottom. Denoting by ˜u^θ(x,˜ t) =˜ φ˜_x_˜(x,˜ z,˜ t)˜|z=θ the non-dimensional horizontal velocity at ˜z=θ, it is clear from (2.23) that ˜u^θis found as

˜

u^θ(x,˜ t) =˜ v(˜ x,˜ t)˜ −1

2βθ²v˜_x_˜_x_˜(x,˜ t) +˜

O

^�β²�

. (2.24)

(22)

Then, by inverting (2.24), the relation

˜

v=u˜^θ+1

2βθ²u˜^θ_x_˜_x_˜+

O

^�β²�

. (2.25)

is found. Replacing ˜v in (2.22) with the the expression from (2.25), the non- dimensional equations

η˜t˜+u˜^θ_x_˜+α� η˜u˜^θ�

˜ x+¹₂β�

θ²−¹₃�

˜

u^θ_x_˜_x_˜_x_˜=

O

^�αβ,β²� ,

˜

u^θ_t_˜+η˜_x_˜+αu˜^θu˜^θ_x_˜−¹₂β�

1−θ²�

˜

u^θ_x_˜_x˜_˜_t=

O

^�αβ,β²�

, (2.26)

appears. This is the non-dimensional version of equations (2.17), except for a difference in the dispersive term in the ﬁrst equation. By using lower order relations (see [11, 56]), in higher order terms, it is possible to further manipulate the system (2.26). Since ˜η˜t=−u˜_x+

O

(α,β), it is seen that it is possible, in this approximation, to replace the dispersive ˜u^θ_x˜_˜_x˜_xterm, with−η˜_x˜_˜x˜t in (2.26), and thus producing the equivalent non-dimensional system

η˜t˜+u˜^θ_x_˜+α� η˜u˜^θ�

˜

x−¹₂β�θ²−¹₃�

η˜_x_˜_x_˜t˜=

O

^�αβ,β²� ,

˜

u^θ_t_˜+η˜_x_˜+αu˜^θu˜^θ_x_˜−¹₂β�

1−θ²�

˜

u^θ_x_˜_x˜_˜_t=

O

^�αβ,β²�

. (2.27)

This is the equations (2.17) on non-dimensional form.

In Paper D, the expression in (2.23), giving the vertical distribution of the horizontal velocity in the ﬂuid column, is exploited when constructing an effective breaking criterion. When deriving expressions for the energy in Paper C, it is advantageous to have gained some insight in the procedure above of deriving the Boussinesq equations.

The dispersion relation corresponding to the system (2.17) is given as ω²(k;θ)

k² =gh₀ 1

�1+¹₂(θ²−¹₃)(kh0)²��

1+¹₂(1−θ²)(kh0)²�. (2.28) The equation (2.28) for the phase speedc=ω/kare derived by inserting waves on the form e^ikx⁻^iωt into the linearized version of the system (2.17) above. It is seen that the dispersion relation depends onθ, and is immediate from the formula that the dispersion relation becomes singular for θ²<1/3. Thus, the relevant values ofθ are given by ¹₃≤θ²≤1. Closer inspection of the dispersion relation above reveals that if θ²=1 or θ²=1/3, the phase speed has quadratic decay as k→∞. As the dispersion relation (2.9) for the full water-wave problem also exhibits quadratic decay whenk→∞, it is clear that values ofθclose to 1/3 and 1, gives a dispersion relation which most closely resemble the target formula (2.9) with which it should be compared.

(23)

2.4 Boussinesq equations 17

−4 −3 −2 −1 0 1 2 3 4

0 0.2 0.4 0.6 0.8 1 1.2

�tanh(kh0) kh0

θ= 0.95 θ=�

7/9

CC₀

kh₀

Figure 2.2: The linear dispersion relation for the systems(2.17)withθ=0.95(dashed line) and withθ=�

7/9(dashed-dotted line), compared with the linear dispersion char- acteristics for the full water wave problem (solid line).

In Figure 2.2, the phase speeds corresponding to the system (2.17) for two different values ofθare compared to the dispersion relation for the full water wave problem (solid curve). It can be seen that the curves corresponding to both values ofθeventually separate from the solid curve representing the target formula. This is expected as the equations (2.17) are derived under the assumptions that the values ofkh₀are small. In Figure 2.2 it is also evident that the curve forθ=0.95 is closer to the solid curve then is the line corresponding to the dispersion relation of (2.17) withθ²=7/9. These values ofθare displayed as they are used in Paper D of the thesis.

(24)

(25)

Chapter 3 The KdV equation

The KdV equation has been useful as a model equation in a variety of contexts, including the study of water waves, particle physics and ﬂow in blood vessels [15, 20, 29, 36, 39], just to name a few. The discovery by Zabusky and Kruskal of the elastic interaction of solitary waves [57], and the subsequent formulation of a solution algorithm by way of solving an inverse-scattering problem [4, 29], excited interest in the equation from both the mathematical and physical point of view. Along with the nonlinear Schr¨odinger equation, the KdV equation has become a paradigm for nonlinear wave equations featuring competing nonlinear and dispersive effects.

In this chapter we will ﬁrst introduce the solitary waves. The KdV equation possesses such solitary waves as solutions. Having an analytical solution at hand is valuable when for instance conﬁrming a numerical implementation of the equations or a convergence rate of a numerical scheme. In the next section, a derivation of the KdV equation is furnished, before, in the last section, solutions of some equivalent forms of the KdV equation are shown. The KdV equation and a generalized version of it, along with the solitary wave solutions, play a key role in the papers A and B of the thesis.

3.1 Solitary waves

The solitary wave was ﬁrst observed on the Edinburgh to Glasgow canal in 1834 by the Scottish engineer John Scott Russell [28]. He called it “the great wave of translation”, and he reported his observations to the British Association in 1834 as this: “[...]. I was observing the motion of a boat which was rapidly drawn along a narrow channel by a pair of horses, when the boat suddenly stopped - not so the mass of water in the channel which it had put in motion; it accumulated round the prow of the vessel in a state of violent agitation, then suddenly leaving

(26)

it behind, rolled forward with great velocity, assuming the form of a large solitary elevation, a rounded, smooth and well-deﬁned heap of water, which continued its course along the channel apparently without change of form or diminution of speed. [...]”. Russell also did laboratory experiments, in which he generated solitary waves by dropping a weight in one end of a channel with water. Based on these experiments, he deduced, among other things, that the velocity of the wave was given as

c²=g(h₀+a₀), (3.1)

wherea₀ was the amplitude of the wave. It is clear that the larger waves travels faster then waves with smaller amplitudes.

Later, both Boussinesq (1871) and Rayleigh (1876), deduced Russell’s empir- ical formula (3.1) forcfrom the equations of motion of an incompressible ﬂuid.

In the derivation they assumed that the ratio of waterdepth to wavelength was small. Moreover, they also showed that the wave could be expressed with the now familiar “sech”-proﬁle

ζ(x,t) =a₀sech²�x−ct b

�, (3.2)

whereb=4h²₀(h₀+a₀)/3a0.

Although Boussinesq and Rayleigh both came up with the surface proﬁle (3.2), they did not derive the equation for which (3.2) is a solution. This equation was derived by Korteweg and deVries in 1895. In the original derivation [39], the effect of surface tension was included, and the original equation was given as

ζt= 3 2

�g h₀

�ζζx+2

3αζx+1 3σζxxx

�. (3.3)

Hereαis a small constant andσ=h³₀/3−T h0/gρ, in whichT denotes the surface tension.

3.2 Derivation

The system of equation (2.17), which on non-dimensional form was found in (2.27) as

η˜t˜+u˜^θ_x_˜+α�η˜u˜^θ�

˜

x−¹₂β�θ²−¹₃�η˜_x_˜x˜_˜t=

O

^�αβ,β²� ,

˜

u^θ_t_˜+η˜x˜+αu˜^θu˜^θ_x_˜−¹₂β�

1−θ²�

˜

u^θ_x_˜_x_˜_t_˜=

O

^�αβ,β²�

, (2.27)

allows for counterpropagating waves. The KdV equation on the other hand, only models waves traveling in one spatial direction. For waves traveling in one direction, it can be shown that the the dependent variables in (2.27) must be related

as u˜^θ=η˜+

O

(α,β), (3.4)

(27)

3.2 Derivation 21

to lowest order in the parametersαandβ. Furthermore, for one-way propagating waves, the relation

η˜˜t+η˜x˜=O(α,β), (3.5) appears at lowest order. Such lower order relations are now employed in a sys- tematic manner to derive the well known KdV equation. This derivation follows along the lines given in Whitham [56], and the startingpoint for the derivation is the system (2.27) found above. Introduce some ﬁrst order, inαandβ, corrections of the dependent variables, as

˜

u^θ=η˜+αf˜+βg˜+

O

^�α²,αβ,β²�

. (3.6)

The function ˜f and ˜g can, because of the lower order relations (3.4) and (3.5) above, be considered functions of ˜ηand ˜x-derivatives of ˜η. Plugging this relation into the system (2.27), the equations

η˜t˜+η˜_x_˜+α�f˜_x_˜+2 ˜ηη˜_x_˜� +β�

˜

g_x_˜−¹₂�θ²−¹₃� η˜_x_˜_x_˜t˜�

=

O

^�α²,αβ,β²� , η˜t˜+η˜_x_˜+α�f˜t˜+η˜η˜_x_˜�

+β�

˜ g˜t−¹₂�

1−θ²� η˜_x_˜_x_˜t˜�

=

O

^�α²,αβ,β²�

, (3.7) appears. Since ˜ηt˜=−η˜_x_˜at lowest order, (3.5), the time derivatives in the higher order terms are replaced by−∂x˜+

O

(α,β)in (3.7) to give

η˜t˜+η˜x˜+α�f˜_x_˜+2 ˜ηη˜x˜� +β�

˜

g_x_˜+¹₂�θ²−¹₃�η˜x˜x˜x˜�

=

O

^�α²,αβ,β²� , η˜t˜+η˜x˜+α�

− f˜_x_˜+η˜η˜x˜� +β�

−g˜_x_˜+¹₂�

1−θ²�η˜x˜x˜x˜�

=

O

^�α²,αβ,β²�

. (3.8) Since both of these equations are correct to ﬁrst order, it is found that

f˜_x_˜=−¹₄� η˜²�

˜ x,

˜

g_x_˜= ¹₂�₂

3−θ²�η˜x˜x˜x˜. (3.9) Using these expressions in either of the equations in (3.8) above, the KdV equation

η˜˜t+η˜x˜+3

2αη˜η˜x˜+1

6βη˜x˜˜x˜x=

O

^�α²,αβ,β²�

. (3.10)

appears. In physical variables, the KdV equation (3.10), is given as ηt+c₀ηx+3

2 c₀

h₀ηηx+1

6c₀h²₀ηxxx=0, (3.11) withh₀being the undisturbed depth andc₀=√gh₀the limiting long-wave speed.

(28)

3.3 KdV equations and their solutions

By using different magniﬁcations and translations of the independent and dependent variables appearing in equation (3.11), the KdV equation can be written in many equivalent forms. If the transformations

X =x−c₀t, T =t , (3.12)

are used, this corresponds to using a frame of reference moving in the direction of the wave at speedc₀. In this frame of reference the equation (3.11) becomes

η_T+3 2

c₀

h₀ηη_X+1

6c₀h²₀η_XXX =0. (3.13) To ﬁnd a non-dimensional equation, it would be natural to choose the length scale to beh0and the time scale ash0/c0. The dimensionless equation

η˜T˜+η˜η˜X˜+η˜X˜X˜X˜ =0, (3.14) appears from (3.13), when the transformations

X˜ =√ 6 X

h₀ , T˜ =√ 6c₀T

h₀ and ˜η= 3 2

η

h₀ , (3.15)

are used. Combined the two transformations above, (3.12) and (3.15), amounts to using

x= h₀

√6(X˜+T˜), t= h₀

√6c0

T˜ and η=2h0

3 η˜ , (3.16) as relations between the physical and dimensionless variables. The solutions of the equations transform accordingly. The solitary wave solution of equation (3.11) could be written (see [40]) as

η(x,t) = 4κ²h³₀

3 sech²�

κ(x−ct)�

, (3.17)

with

c=c0

�1+2 3h²₀κ²�

. (3.18)

The solitary wave solution of equation (3.14) is found from (3.17), via the relations (3.16), as

η(˜ X˜,T˜) = 3

2h0η� h₀

√6(X˜+T˜), h₀

√6c₀T˜�

= 3 2h0

4κ²h³₀

3 sech²� κ�h₀

√6(X˜+T˜)−c h₀

√6c0

T˜��

=12˜κ²sech²�

κ(˜ X˜−C˜T˜)� ,

(3.19)

(29)

3.3 KdV equations and their solutions 23

with

C˜= c

c₀−1 and ˜κ= h0

√6κ. (3.20)

From equations (3.18) and (3.20) it is seen that

C˜=4˜κ², (3.21)

showing that the speed depends on the amplitude for these waves. Once either the velocity or the amplitude of the wave are determined, the KdV solitary wave is given.

The KdV equation along with a generalized version of it has been used in the papers A and B. In those papers, solitary wave solutions, like (3.19), have been used to conﬁrm our numerical implemention, and especially to conﬁrm the rate of convergence of the numerical scheme. The solitary wave solutions are also used as a test case for our method to estimate the radius of analyticity of the solutions of the KdV equation and the generalized KdV equation. In for instance [12], solitary wave solutions of the latter equation are found.

(30)

(31)

Chapter 4 Numerical methods

In this chapter we will elaborate upon the numerical methods used in the thesis.

The content covers both different discretization techniques applied on the equations, as well as methods developed to investigate for instance the radius of analyticity of solutions of the KdV equation. The ﬁrst section will be concerned with Fourier spectral methods and the KdV equation, and can be regarded as a summary of the papers A and B of the thesis. In Section 4.2 numerical methods for solving a Boussinesq system of equations are presented along with results from a benchmarking process. This section summarizes the developments and results from Paper E.

4.1 Fourier spectral methods and the KdV equation

Since the discovery by Cooley and Tukey of a fast algorithm to compute the discrete Fourier transform [25], spectral methods based on the Fast Fourier Trans- form have become a popular choice for the spatial discretization of nonlinear par- tial differential equations. In particular, in wave propagation problems, spectral projections have been widely used in connection with the Fourier basis. In Papers A and B [7, 8] we have had occasions to apply this discretization technique on the KdV equation and a generalized version of it.

4.1.1 The Fourier spectral method

In this section, a derivation of a Fourier spectral method for the periodic KdV equation on the form

ηt+ηηx+ηxxx=0,

η(x,0) =η₀(x), (4.1)

(32)

is given. The numerical method is the one used in the numerical studies of Pa- per A. Without much effort the method could be adapted to the generalized KdV equation considered in Paper B. The variables appearing in Section 4.1 are non- dimensional.

The aim is to arrive at a semidiscrete set of ordinary differential equations. To achieve this, the ﬁnite dimensional space

S_N =�

e^ikx ��k∈Z, −N≤k≤N�

(4.2) of L²(0,2π) is used. The Galerkin approximation to the equation (4.1) is then given by a functionu_N from[0,T]toS_N satisfying the system of ordinary differential equations given as

� �∂tηN+¹₂∂x(η²_N) +∂³_xηN,φ�

= 0, t∈[0,T],

η_N(0) = P_Nη₀, (4.3) for all φ∈S_N. In (4.3), η₀ is the periodic, initial data. The operator P_N is the orthogonal projection

P_Nf(x) =

∑

−N≤k≤N

e^ikxfˆ(k), (4.4)

from L² onto SN. The coefficients ˆf(k) appearing in the equation (4.4), are the Fourier coefficients defined by

fˆ(k) = 1 2π

� _2π

0 e⁻^ikxf(x)dx, (4.5)

for 2π-periodic functions f. Also, as inner product on the space L²(0,2π), the expression

(f,g) =^� ^2π

0 f(x)g(x)dx. (4.6) is used.

In practice the Fourier coefﬁcients in (4.5) are evaluated discretely by the use of the discrete Fourier transform (dft). In the dft, the integrals are evaluated with the trapezoidal method, which for periodic functions corresponds to a Gauss quadrature rule (see [53]). The dft could be achieved in the order of

O

^(NlogN) operations through the use of thefast Fourier transform (fft) [25, 53]. Since the fft-algorithm is more effective when N is a product of small prime factors, it is common to letNbe even and some power of 2. In practice then, it is assumed that the solution of (4.3) is written as the sum

ηN(x,t) = ^N/2

∑

k=−N/2+1

ηˆ_k(t)e^ikx, (4.7)

(33)

4.1 Fourier spectral methods and the KdV equation 27

and the coefﬁcients in (4.7) is assumed evaluated discretely by the fft-algorithm.

The N equations deﬁning the N unknowns ˆηk are systematically derived from (4.3) by choosingφ=e^ilx forl=−N/2+1, . . . ,N/2. The equations become

0=

N/2

∑

k=−N/2+1

�∂_tηˆ_k(t)

� _2π

0 e^ikxe⁻^ilxdx+ηˆ_k(t)^� ^2π

0 (∂³_xe^ikx)e⁻^ilxdx +1

2ηˆ²_k(t)^� ^2π

0 (∂_xe^ikx)e⁻^ilxdx� .

(4.8)

The linear part of the problem is comprised of the two ﬁrst terms on the right hand side of (4.8). These terms can be evaluated analytical. Due to the orthogonality of the complex exponentials, the integrals are zero unlessk=l. The quadratic term, η², is evaluated on the grid in physical space, before the product is transformed to the Fourier side by the fft. The Fourier coefﬁcients ofη² is denoted by ˆη²_k in (4.8). From (4.8) the system of ordinary differential equations,

ηˆ^�_k(t) =ik³ηˆ_k(t)−1

2ikηˆ²_k(t), k=−N/2+1,−N/2+2, . . . ,N/2, (4.9) emerges. Except for the use of the discrete Fourier transform, the spatial discretization of the linear part of equation (4.1) is exact. The non-linear term, however, might introduce aliasing effects. This can be remedied by the use of an de-aliasing method. The 3/2-rule from [19] has been implemented in (4.8) above.

With a large number of modes in the computations, we did not spot any difference between the aliased and de-aliased solutions of (4.3).

The system (4.9) can be solved by a number of methods. The explicit fourth- order Runge-Kutta method could for instance be used. This time marching technique has been employed on the Boussinesq systems in Paper E in the thesis.

In papers A and B, the Cranc-Nicholson method was used on the linear part of (4.9), while the non-linear part was solved by an iterative scheme. With a small time-step, the scheme needed only one iteration to converge.

4.1.2 Exponential convergence

In Paper A it is proved that the spectral projection (4.3) of the KdV equation, converges exponentially fast towards the exact solution of the problem. Such a result was also reported in [37]. The convergence result is stated in Theorem 4.1.1 below. The result in the theorem holds for solutions of (4.1) which are real- analytic functions of the spatial variablex. To quantify the domain of analyticity, the class of periodic analytic Gevrey spaces as introduced by Foias and Temam in [32] is used. Forσ>0, the Gevrey norm� · �Gσ is deﬁned by

�f�²G_σ =

∑

k∈Z

e^2σ√

1+|k|²|fˆ(k)|², (4.10)

(34)

where ˆf(k) are the Fourier coefﬁcients of the function f. Functions in the space G_σ are analytic in a strip of width 2σabout the real axis, and in the next section a method to estimate σbased on the convergence result in the Theorem 4.1.1 is given.

Theorem 4.1.1 Suppose η(x,t) is a continuous function of t with values in the space G_σ for some σ>0, and suppose that ηsolves(4.1). Given T >0 and an integer N>0, there exists a unique solutionηN of the ﬁnite-dimensional problem (4.3). Moreover, there is a constantΛ_T, such that

t∈sup[0,T]�η(·,t)−η_N(·,t)� ≤ Λ_TNe^−σN.

A similar convergence result also holds for spectral projections of the generalized KdV equation. The proof for the generalized KdV equation is contained in Paper B. One class of functions to which Theorem 4.1.1 applies directly is the family of periodic traveling waves of the KdV equation. Note however that in general, even if initial dataη(x,0)are given in the periodic analytic Gevrey spaceG_σ, it cannot be concluded that the solutionη(x,t)is inG_σ at any other timet >0. However, it can be shown thatη(x,t)is in G_σ_t, whereσ_t is a decreasing function oft. In the case of the real line, algebraic dependence ofσt ont has been recently proved [13]. In Section 4.1.3 below, we report on a method for determining the radius of analyticity (σ) of solutions of (4.1). One of the motivations for this work, was to give numerical evidence of the results in [13], where an algebraic dependence of σont was found.

4.1.3 Radius of analyticity

In this section a method for determining the radius of analyticity of solutions of equation (4.1) is given. The aim of the method is thus to estimateσ. The method is based on the convergence result given in Theorem 4.1.1. Basically σis found by computing theL²-norm of the difference between the numerical and the true solution of the equation (4.1), and then assuming that this difference should be given on the formΛTNe^−σN as indicated by the result in Theorem 4.1.1. A least square ﬁt between the computedL²-errors and the assumed functional form of the errors, gives an estimate of the radius of analyticity.

The centerpiece of the method is given by a robust way of computing the error committed by the Galerkin projection. Thus for ﬁxedN, an approximation of the integral

�_� _2π

0 |η(x,t)−ηN(x,t)|²dx

�1/2

, (4.11)

(35)

4.1 Fourier spectral methods and the KdV equation 29

is sought. The approximation of the integral (4.11) is often done on the grid used to obtain the numerical solution. That is, using a resolutionN corresponding to a grid spacingh_N =2π/N, theL²-error,E_D, is usually deﬁned as

E_D²(N) =h_N

∑

^N

j=1|η(xj,t)−ηN(x_j,t)|². (4.12) A problem with this deﬁnition is that both the error from the numerical method and the error from the approximation of the integral enter into the formula. How- ever, one really wants to remove the part ofEDwhich is not due to the numerical projection of the differential equation. In Paper B an alternative deﬁnition of the L²-errorE_F, namely

E_F²(N) =h_M

∑

^M

j=1|η(x_j,t)−η^M_N(x_j,t)|², (4.13) is introduced. Here h_M =2π/M and it is assumed that M>>N. The discrete solution on the ﬁner spacedM-grid is found by evaluating the trigonometric in- terpolant (see (4.7)) on the grid. While the difference between usingE_DandE_F appears unimportant when validating the convergence of a numerical implementation, the second deﬁnition is superior in determining the radius of analyticityσ from the convergence rate.

To illustrate the method of estimatingσ, and at the same time illuminate the difference between usingED orEF as deﬁnitions of the L²-errors, an example is now presented. In the example a solitary wave η(x,t) =ψ(x−ct) with transla- tional velocityc, is used. In order to analyze solitary wave solutions on a periodic domain, a parameterais introduced, and the equation (4.1) is written as

η_t+ηη_x+ 1

a³η_xxx=0, (4.14)

through the scaling η(x,t) =v(ax,t)/a. For equation (4.14), the solitary wave form is found whenψtakes the special form

ψ(x) =12κ²

a sech²(κax), (4.15) for anyκ>0, and with velocityc=4κ²/a. The Fourier transform ofψis given as

ψ(k) =ˆ 6πkcsch� πk 2κa

�

, (4.16)

revealing that ψ is in the spaceG_σ as long as σ<π/2κa. Thus the aim of the method is to ﬁnd an estimate ofσwhich is close toπ/2κa.

Numerical Methods for the Discretization of Model Equations for Surface Water Waves: Theory and Applications