Contents

Miko laj Szydlarski

CHIARA SIMEONI

List of Figures

Simply stated, open-channel flow is the flow of a liquid in a conduit with a free surface. There are many practical examples,both artificial (flumes, spillways, canals, weirs, drainage ditches, culverts) and natural (streams, rivers, estuaries, floodplains). This section introduces the elementary analysis of such flows, which are dominated by the effects of gravity.

The presence of the free surface, which is essentially at atmospheric pressure, both helps and hurts the analysis. It helps because the pressure can be taken constant along the free surface, which therefore is equivalent to the hydraulic grade line of the flow. Unlike flow in closed ducts, the pressure gradient is not a direct factor in open-channel flow, where the balance of forces is confined to gravity and friction. But the free surface complicates the analysis because its shape is a priori unknown, the depth profile changes with conditions and must be computed as part of the problem, especially in unsteady problems involving wave motion (in our case, long waves in shallow water phenomena).

An open channel always has two sides and a bottom, where the flow satisfies the no-slip condition. Therefore even a straight channel has a threedimensional velocity distribution. Some measurements of straight-channel (From Ref. [START_REF] Frank | Fluid Mechanics Mc-Graw-Hill[END_REF] ) velocity contours are shown in Fig (1.1). The profiles are quite complex,with maximum velocity typically occurring in the mid-plane about 20 percent below the surface. In very broad shallow channels the maximum velocity is near the surface, and the velocity profile is nearly logarithmic from the bottom to the free surface. In noncircular channels there are also secondary motions. If the channel curves or meanders, the secondary motion intensifies due to centrifugal effects, with high velocity occurring near the outer radius of the bend. Curved natural channels are subject to strong bottom erosion and deposition effects. With the advent of the supercomputer, it is possible to make numerical simulations of complex flow patterns such as in Fig. (1.1).

Miko laj Szydlarski

However, the practical engineering approach, that is used here, is to make a one-dimensional-flow approximation, as shown in Fig. (1.2). Since

Miko laj Szydlarski

Flow classification by depth variation

The most common method of classifying open-channel flows is by the rate of change of the free-surface depth. The simplest and most widely analyzed case is uniform flow, where the depth (hence the velocity in steady flow) remains constant. Uniform-flow conditions are approximated by long straight runs of constant-slope and constant-area channel.

If the channel slope or cross section changes or there is an obstruction in the flow, then the depth changes and the flow is said to be varied.Theflowis gradually varying if the one-dimensional approximation is valid and rapid ly varying not. Some examples of this method of classification are shown in Fig. (1.3). The classes can be summarized as follows:

1. Uniform flow (constant depth and slope)

Varied flow

• Gradually varied (one-dimensional)

• Rapidly varied (multidimensional)

Typically uniform flow is separated from rapidly varying flow by a region of gradually varied flow. Gradually varied flow can be analyzed by a first-order differential equation ,but rapidly varying flow usually requires experimentation or three-dimensional potential theory. The Froude-number denominator (gy) 1/2 is the speed of an infinitesimal shallow-water surface wave. We can derive this with reference to Fig. (1.4a), which shows a wave of height δy propagating at speed c into still liquid.

To achieve a steady-flow inertial frame of reference, we fix the coordinates on the wave as in Fig. (1.4b), so that the still water moves to the right at velocity c.

For the control volume of (1.5)

The "stronger" the wave height δy, the faster the wave speed c.I nt h el i m i t of an infinitesimal wave height δy → 0, the speed becomes:

c 2 0 = gy. (1.6)
This is the surface-wave equivalent of fluid sound speed, and thus the Froude number in channel flow Fr = V /c 0 is the analog of the Mach number.

As in gas dynamics, a channel flow can accelerate from subcritical to critical to supercritical flow and then return to subcritical flow through a sort of normal shock called a hydraulic jump. This is illustrated in Fig. (1.5) The flow upstream of the sluice gate is subcritical. It then accelerates to critical and supercritical flow as it passes under the gate,which serves as (From Ref. [START_REF] Frank | Fluid Mechanics Mc-Graw-Hill[END_REF] ) a sort of "nozzle".Further downstream the flow "shocks" back to subcritical flow because the downstream "receiver" height is too high to maintain supercritical flow.

Miko laj Szydlarski

As suggested by Bakhmeteff 1 in 1911, the specific energy E is a useful parameter in channel flow, it is defined as

E = y + V 2 2g . (1.7)
Where y is the water depth. It is seen from Fig. (1.6), that E is the height of the energy grade line (EGL) above the channel bottom. For a given flow rate, there are usually two states possible for the same specific energy. Consider the possible states at a given location. Let q = Vy be the discharge per unit width of a rectangular channel. Then, with q constant, Eq. (1.7) becomes:

E = y + q 2 2gy 2 .
(1.8)

Figure (1.6b) is a plot of y versus E for constant q from Eq. (1.8). There setting dE/dy = 0 at constant q, we find that E min occurs at

y = y c = q 2 g 1 3
.

(1.9)

The critical depth y c is sketched as a dashed line in Fig.

E min = E(y c )= 3 2 y c (1.10)
The critical depth y c corresponds to some critical channel velocity V c , that turns out to be equal to the shallow-water wave propagation speed c 0 from Eq. (1.6). To see this, we rewrite Eq. (1.9) as

q 2 = gy 3 c =(gy c )y 2 c = V 2 c y 2 c (1.11)
By comparasion it follows that the critical channel velocity is Miko laj Szydlarski

V c =(gy c ) 1 2 = c 0 ,F r =1. ( 1

Physical and Mathematical Modelling

Upstream is a zone of silence, and a small obstruction in the flow will create a wedge-shaped wave. The angle of these waves must be

µ =sin -1 c 0 V =sin -1 (gy 1 2 ) V (1.13)
The wave angle and the depth can thus be used as a simple measurement of supercritical-flow velocity.

V 1 y 1 = V 2 y 2 V 2 1 2g + y 1 = V 2 2 2g + y 2 + ∆h (1.14)
Eliminating V 2 between these two gives a cubic polynomial equation for the water depth y 2 over the bump:

y 3 2 -E 2 y 2 2 + V 2 1 y 2 1 2g =0 where E 2 = V 2 1 2g + y 1 -∆h (1.15)
This equation has one negative and two positive solutions if ∆h is not too large. Its behavior is illustrated in Fig. (1.7b) and depends upon whether condition 1 is on the upper or lower leg of the energy curve. The specific energy E 2 is exactly ∆h less than the approach energy E 1 ,and point 2 will lie on the same leg of the curve as E 1 . A subcritical approach, Fr 1 < 1.0, will cause the water level to decrease at the bump. Supercritical approach flow, Fr 1 > 1.0, causes a water-level increase over the bump.

If the bump height reaches

∆h max = E 1 -E c , as illustrated in Fig. (1.7b),
the flow at the crest will be exactly critical (Fr > 1.0). If ∆h>∆h max ,there are no smooth physically correct solutions to Eq. (1.15). That is, a bump too large will "choke" the channel and cause frictional effects, typically a discontinuous solution as a hydraulic jump (see Subsection 1.1.3).

Miko laj Szydlarski 15

CHIARA SIMEONI Being extremely turbulent and agitated, the hydraulic jump is a very effective energy dissipator and is a feature of stilling-basin and spillway applications. Figure (1.8) shows the naturally occurring hydraulic jump formed at the bottom of a river.

The principal parameter affecting hydraulic-jump performance is the upstream Froude number

Fr 1 = V 1 /(gy 1 ) 1 2
. The Reynolds number and channel geometry for real flows, have only a secondary effect. The following ranges of operation can be outlined, as illustrated in Fig. (1.9).

Miko laj Szydlarski 17

CHIARA SIMEONI

Physical and Mathematical Modelling

Miko laj Szydlarski 19

CHIARA SIMEONI

Physical and Mathematical Modelling

Fr 1 < 1.0: Jump impossible, violates second law of thermodynamics.

Fr 1 =1.0 to 1. A jump which occurs on a steep channel slope can be affected by the difference in water-weight components along the flow. The effect is small, however, so that the classic theory assumes that the jump occurs on a horizontal bottom. We have already analyzed this problem in Sec.

V 2 1 = 1 2 gy 1 η(η + 1) (1.16)
where η = y 2 /y 1 . Introducing the Froude number

Fr 1 = V 1 /(gy 1 ) 1
2 and solving this quadratic equation for η, we obtain:

2y 2 y 1 = -1 + (1 + 8Fr 2 1 ) 1 2
(1.17)

Miko laj Szydlarski 20

CHIARA SIMEONI

Physical and Mathematical Modelling

With y 2 thus known, V 2 following form the wide-channel continuity relation:

V 2 = V 1 y 1 y 2 (1.18)
Finally, we can evaluate the dissipation loss across the jump from the steadyflow energy equation:

∆E = E 1 -E 2 = y 1 + V 2 1 2g y 2 + V 2 2 2g (1.19)
Introducing y 2 and V 2 from Eqs. (1.17) and (1.18), we find after considerable algebraic manipulation that:

∆E = (y 2 -y 1 ) 3 4y 1 y 2 (1.20)
Equation (1.20) shows that the dissipation loss is positive only if y 2 >y 1 .

which is a requirement of the second law of thermodynamic. Equation (1.17) then requires that Fr 1 > 1.0; that is, the upstream flow must be supercritical. Finally, Eq. (1.18) shows that V 2 <V 1 and the downstream flow is subcritical.

Miko laj Szydlarski 21

CHIARA SIMEONI

Physical and Mathematical Modelling

Systems of conservation laws

Before we pose a model of shallow water, we will focus on systems of conservation laws, as the shallow water equations are a specific case of such systems.

To see how conservation laws arise from physical principles, we will consider the fundamental case of deriving the equation for conservation of mass in a one-dimensional gas dynamic problem, for example the flow in a tube, where properties of the gas such as density and velocity are assumed to be constant across each section of the tube. Then we will show how this gas-model is related with shallow water system.

Let x represents the distance along the tube and let ρ(x, t)bethedensity of the gas at point x and time t. This density is defined in such way that the total mass of gas in any given section between x 1 and x 2 , is given by the integral of density:

mass in [x 1 ,x 2 ] at time t = x 2 x 1 ρ(x, t)dx (1.21)
If we assume that the walls of the tube are impermeable and the mass is neither created nor destroyed, then the mass in one reference section can change only because of gas flowing across the endpoints x 1 or x 2 . Now let v(x, t) be the velocity of the gas at point x and time t,t h e nt h e rate of flow, (or flux) of gas, past trough this point is given by mass flux at (x, t)=ρ(x, t)v(x, t).

(1.22) By our comments above, the rate of change of mass in [x 1 ,x 2 ] is given by the difference in fluxes at x 1 and x 2 , namely d dt Note that these equations involve another quantity, the pressure p,w h i c h must be specified as a given function of ρ, ρv, and E in order for the fluxes to be well defined functions of the conserved quantities alone. This additional, specific equation is called equation of state and depends on physical properties of the gas under study.

x 2 x 1 ρ(x, t)dx = ρ(x 2 ,t)v(x 2 ,t) -ρ(x 1 ,t)v(x 1 ,t) (1.
If we introduce the vector u ∈ R 3 as

u(x, t)=     ρ(x, t) ρ(x, t)v(x, t) E(x, t)     ≡     u 1 u 2 u 3     (1.32)
then the system of equations (1.28), (1.30), (1.31) can be written simply as

u t + f (u) x = 0 (1.33)
where

f (u)=     ρv ρv 2 + p v(E + p)     =     u 2 u 2 2 u 1 + p(u) u 2 (u 3 +p(u)) u 1     .
(1.34)

Again, the form (1.33)i st h edifferential form of the conservation laws, which holds in the usual sense only when u is smooth. More generally, the integral form for a system of m equations reads d dt

x 2 x 1 u(x, t)dx = f (u(x 2 ,t)) -f (u(x 1 ,t)) (1.35)
for all x 1 , x 2 , t. Equivalently, integrating from t 1 to t 2 gives

x 2 x 1 u(x, t 2 )dx = x 2 x 1 u(x, t 1 )dx + t 2 t 1 f (u(x 2 ,t))dt - t 2 t 1 f (u(x 1 ,t))dt, (1.36)
for all x 1 , x 2 , t 1 and t 2 . These integral forms of the conservation law will be fundamental in later analysis.

Miko laj Szydlarski

Euler equations and Isentropic flows

The Euler equations of gas dynamics are a particularly important example because it is closely connected with the shallow water model. The continuity equation (conservation of mass) has been derived in the previous section. Here, we will consider the momentum and energy equation in more detail, as well as the equation of state and a few other quantities of physical (and mathematical) significance, such as the entropy. We will also look at some simplifications namely, the isentropic case, which moves directly to the shallow water model. We "sketch" the derivation here, with an emphasis on the main ideas. A more thorough introduction can be found in Whitham [START_REF] Whitham | Linear and Nonlinear Waves[END_REF].

Recall that ρ is the density, v the velocity, E the total energy, and p the pressure of the gas. The continuity equation is reads

ρ t +(ρv) x =0, (1.37)
where the mass flux is given by ρv. More generally, for any quantity z that is advected with the flow there will be a contribution to the flux for z of the form zv. Thus, the momentum equation has a contribution Ev.

In addition to advection, there are forces on the fluid that cause acceleration due to Netwon's laws, and hence changes in momentum. If there are no outside forces, then the only force is due to variations in the fluid itself, and is proportional to the pressure gradient which is simply p x in one dimension.

Combining this with the advective flux gives the momentum equation

(ρv) t +(ρv 2 + p) x =0. (1.38)
The total energy E is often decomposed as

E = 1 2 ρv 2 + ρe, (1.39)
where the first term is the kinetic energy, while ρe is the internal energy.

The variable e, internal energy per unit mass, is called the specific internal energy (In general "specific" means "per unit mass"). Internal energy includes rational and vibrational energy and possibly other form of energy in more complicated situations. In Euler equations we assume that the gas is

Miko laj Szydlarski 25

CHIARA SIMEONI

Physical and Mathematical Modelling

in chemical and thermodynamic equilibrium and that the internal energy is a known function of pressure and density, that is e = e(p, ρ).

(1.40) This is the "equation of state" for the gas, which depends on the particular gas under study.

The total energy advects with the flow, but is also modified due to work done on the system. In the absence of outside forces, work is done only by the pressure forces and is proportional to the gradient in one dimension, of vp. The conservation law for total energy thus takes the following form

E t +[v(E + p)] x =0. (1.41)
Putting these equations together gives the system of Euler equations

    ρ ρv E     t +     ρ ρv 2 + p v(E + p)     x =0. (1.42)
We still need to specify the equation of state relating the internal energy to pressure and density. For an ideal gas, internal energy is a function of temperature alone, e = e(T ), and T is related to p and ρ by the ideal gas law,

p = RρT, (1.43)
where R is a constant. To good approximation, the internal energy is simply proportional to the temperature, Note that by the ideal gas law,

e = c v T, (1.44
c p -c v = R. (1.50)
The equation of state for a polytropic gas turns out to depend only on the ratio of specific heats, usually denoted by

γ = c p c v . (1.51)
Internal energy in a molecule is typically split up between various degrees of freedom (translational, rotational, vibrational, etc.). How many degrees of freedom exist depends on the nature of gas. The general principle of equipartition of energy says that the avarage energy in each of these is the same. Each degree of freedom contributes an average energy of 1 2 kT per molecule, where k is Boltzmann's constant. This gives a total contribution of α 2 kT per molecule if there are α degrees of freedom. Multiplying this by n, the number of molecules per unit mass (which depends on the gas), gives e = α 2 nkT.

(1.52)

The product nk is precisely the gas constant R, so comparing this to (1.44) gives

c v = α 2 R (1.53)

Miko laj Szydlarski

From (1.50) we obtain

c p = 1+ α 2 R , (1.54)
and so

γ = c p c v = α +2 α (1.55)
For monatomic gas the only degrees of freedom are the three translational degrees so α = 3 and γ = 5 3 . For diatomic gas (such as air, which is composed primarily on N 2 and O 2 ), there are also two rotational degrees of freedom and α = 5, so that γ = 7 5 . Note that T = p/Rρ so that by (1.50) and (1.51),

e = c v T = c v R p ρ = p (γ -1)ρ
, equation of state for a polytropic gas.

(1.56)

Using this in (1.39) gives the common form of the equation of state for a polytropic gas:

E = p γ -1 + 1 2 ρv 2 (1.57)
Another important thermodynamic quantity is the entropy. Roughly speaking, this measures the disorder in the system. The entropy S is defined up to an additive constant by

S = c v log( p ρ γ ) + constant. (1.58)
This can be solved for p to give

p = κe S cv ρ γ , (1.59)
where κ is a constant.

From the Euler equations we can derive the following relation

S t + vS x =0, (1.60)
which says that entropy is constant along particle path in region of smooth flow. In fact (1.60) can be derived from fundamental principles and this equation, together with the conservation of mass and momentum equation, gives an alternative formulation of the Euler equations (though not in conservation form):

ρ t +(ρv) x =0 (ρv) t +(ρv 2 + p) x =0 S t + vS x =0
(1.61)

Miko laj Szydlarski

It turns out that equation of state then gives p as a function of ρ and S alone, e.g. (1.59) for polytropic gas. In this form the partial derivative of p with respect to ρ (holding S fixed) plays a fundamental role: its square root c is the local speed of sound in the gas. For a polytropic gas we have

c 2 = ∂p ∂ρ S=constant = γκe S cv ρ γ-1 = γ p ρ (1.62)
and so

c = γp ρ (1.63)
From our standpoint the most important property of entropy is the fact that is smooth flow entropy remains constant on each particle patch, while if a particle cross a shock then the entropy may jump, but must increase. This is the physical entropy condition for shocks. Note that, along a particle path of smooth flow, since S is constant we find by (1.59) that

p =κρ γ , (1.64)
where κ = κe S/cv is a constant which depends only on the initial entropy of the particles. This explicit relation between density and pressure along particle paths is sometimes useful. Of course, if the initial entropy varies in space then κ will be different along different particle paths.

If the entropy is constant everywhere then (1.64) holds with the same value of κ and the Euler equations simplify. This is the case, for example, if we consider fluid flows that start at uniform rest state (so S is constant) and remains smooth (so S remains constant). Then using (1.64), the equation of state (1.57) reduces to en explicit expression for E in terms of ρ and ρv.The energy equation then becomes redundant and the Euler equations reduce to a system of two equations, the equations of isentropic gas dynamics,

ρ ρv t + ρv ρv 2 +κρ γ x =0.
(1.65)

Shallow water equations

To derive the one-dimensional shallow water equations, we consider fluid in a channel of unit width and assume that the vertical velocity of the fluid is

Miko laj Szydlarski 29

negligible and the horizontal velocity v(x, t) is roughly constant throughout any cross section of the channel. This is true if we consider small-amplitude waves in a fluid that is shallow relative to the wavelength. We now assume the fluid is incompressible, so the density ρ is constant. Instead we allow the depth of the fluid to vary, and it is this depth, or height h(x, t), that we wish to determine. In analogy with (1.21), the total mass in [x 1 ,x 2 ] at time t is given by

x 2
x 1 ρh(x, t)dx.

(1.66)

The density of momentum at each point is ρv(x, t), and integrating this vertically gives the mass flux to be ρv(x, t)h(x, t). The constant ρ drops out of the conservation of mass equation (1.28), so that we get

h t +(vh) x = 0 (1.67)
The quantity hv is often called discharge in shallow water theory, since it measures the rate of water past a point.

The conservation of momentum also takes form as in Euler equations

(1.38), namely

(ρhv) t +(ρhv 2 + p) x =0, (1.68)
but now p is determined from a hydrostatic law, stating that the pressure at distance h -y below the surface is ρg(h -y), where g is the gravitational constant. This pressure arises simply from the weight of the fluid above.

Integrating this vertically from y =0toy = h(x, t) gives the total pressure felt at (x, t), the proper pressure term in the momentum flux:

p = 1 2 ρgh 2 . (1.69)
Using this in (1.68) and canceling ρ out gives

(hv) t + hv 2 + 1 2 gh 2 x =0. ( 1

Miko laj Szydlarski

Note that this is equivalent to the equations of isentropic gas dynamics (discussed in Section 1.2.1) with the value γ = 2, since setting p(ρ)= 1 2 gh 2 in (1.61) gives the same system.

If we assume that h and v are smooth, then equation (1.70) can be simplified by expanding the derivatives and using (1.67) to replace the h t term. Then several terms drop out, and (1.70)i sr e d u c e dt o

v t + 1 2 v 2 + gh x =0. (1.72)
This equation is equivalent to the previous set (1.71) for smooth solutions, but it is important to note that the manipulations performed above depend on smoothness. For problems with shock waves, the two sets of conservation laws are not equivalent, and we know that it is crucial that we use the correct set in calculating shock waves. The form (1.71), which is derived directly from the original integral equations, is the correct set to use.

Since we will be interested in studying shock waves, we use the form (1.71) and take

q(x, t)= h hv = q 1 q 2 , f (q)= hv hv 2 + 1 2 gh 2 = q 2 (q 2 ) 2 q 1 + 1 2 g(q 1 ) 2 .
(1.73)

For smooth solutions, these equations can equivalently be rewritten in the quasi-linear form

q t + f (q)q x =0, (1.74)
where the Jacobian matrix f (q) is given by

f (q)=   01 - (q 2 ) q 1 2 + gq 1 2 q 2 q 1   = 01 -v 2 + gh 2v . (1.75)
The eigenvalues of f (q) are

λ 1 = v - gh, λ 2 = v + gh, (1.76)
with corresponding eigenvectors

r 1 = 1 v - √ gh ,r 2 = 1 v + √ gh .
(1.77)

Note that the eigenvalues and eigenvectors are functions of q for this nonlinear system.

Remarks:T h es y s t e m( 1.74) is called hyperbolic if the matrix f (q)i s diagonalizable with real eigenvalues and it is called strictly hyperbolic if the eigenvalues are distinct.

If we wish to study waves with very small amplitude, then we can linearize these equations to obtain a linear system. Suppose the fluid is essentially at constant depth h 0 > 0 and moving at constant velocity v 0 (which may be zero), and let q now represent the perturbations from this constant state, so

q = h -h 0 hu -h 0 v 0 and q 0 = h 0 h 0 v 0 .
Then expanding the flux function and dropping terms of O(q 2 ) gives the linear system q t + Aq x =0w h e r eA = f q 0 . Hence small-amplitude waves move at the characteristic velocities λ 1 0 = v 0 -c 0 and λ 2 0 = v 0 + c 0 ,w h e r e c 0 = √ gh 0 . These waves propagate at speed ±c 0 relative to the fluid. These shallow water waves should not be confused with sound waves, however.

Sound does propagate in water, due to its slight compressibility, but in the shallow water equations we are ignoring this compressibility and hence ignoring sound waves. The waves we are modeling are often called gravity waves, since they are driven by the hydrostatic pressure resulting from gravity.

They typically propagate at a speed √ gh that is much less than the speed of sound in water.

Note that λ 1 and λ 2 can be of either sign, depending on the magnitude of v relative to c. In shallow water theory the ratio

Fr = |v| c (1.78) is called the Froude number (see subsection 1.1.2)
The wave speed √ gh 0 depends on the depth of the fluid; waves in deeper water move faster. Note that within a wave the depth of the fluid varies (it is deeper at a crest than in a trough), and so we should expect the crest of a wave to propagate slightly faster than a trough. If the amplitude of the wave is very small compared to h 0 , then we can safely ignore this slight variation in speed, which is what we do in linearizing the equations. Then

Miko laj Szydlarski

Physical and Mathematical Modelling

all parts of the wave travel at the same speed based on the background depth h 0 , and the wave propagates with its shape unchanged. For waves with larger amplitude, however, the deformation of the wave due to differing wave speeds may be quite noticeable. In this case the linearized equations will not be an adequate model and the full nonlinear equations must be solved.

The nonlinear distortion of a wave leads to a steepening of the wave in the region where the fast-moving crest is catching up with the slower trough ahead of it (a compression wave), and a flattening of the wave (an expansion or rarefaction) in the region where the crest is pulling away from the following trough.

This behavior is familiar from watching waves break on the beach. Far from shore the waves we normally observe have a wavelength that is very small compared to the water depth, and hence they are governed by surfacewave theory rather than shallow water theory. Near the beach, however, the water depth is small enough that nonlinear shallow water theory applies.

In this shallow water, the difference in h between crests and troughs is significant and the waves steepen. In fact the crest is often observed to move beyond the position of the preceding trough. At this point the assumptions of shallow water theory no longer hold, and a more complicated set of equations would have to be used to model breakers. Beyond the breaking time the depth h is triple-valued, a situation that obviously can't occur with other systems of conservation laws (such as gas dynamics) where the corresponding variable is a density that must be single-valued.

This extreme behavior of breaking waves results from the additional complication of a sloping beach. This leads to a continuous decrease in the fluid depth seen by the wave and a severe accentuation of the nonlinear effects.

(The sloping beach, or more generally any variation in the bottom topography, also leads to additional source terms in the shallow water equations.)

Shallow water waves in a domain with a flat bottom will typically not exhibit this type of breakers. Instead the gradual steepening of the wave due to nonlinearity would be counterbalanced by other effects such as surface tension (and also the vertical velocity, which is ignored in the one-dimensional model). Modeling these other effects would lead to higher-order derivatives

Miko laj Szydlarski 33

CHIARA SIMEONI

Physical and Mathematical Modelling

in the equations (with small coefficients) and consequently the equations would have smooth solutions for all time. When these coefficients are small, the wave can become nearly discontinuous, and the shock-wave solution to the hyperbolic system gives a good approximation to such solutions. In shallow water flow, a shock wave is often called a hydraulic jump A number of fundamental assumptions are inherent within the model and these can be summarized as follows:

The Saint-Venant System

• the flow is one-dimensional, the mean velocity is constant over a cross section and the water level is horizontal;

• the vertical component of the acceleration of the fluid is negligible so that the pressure variation with depth is hydrostatic;

Miko laj Szydlarski

• friction and turbulence can be represented using the same empirical laws that govern steady state flows.

Summarizing, the Saint-Venant equations are a particular case of shallow water equations described in Section (1.2.2), they are commonly used to describe physical situations such as flows in rivers or costal areas, specially the one-dimensional version is well adapted for ideal rectangular rivers. The system allows to describes the flow, at time t ≥ 0 and at point x ∈ R, through the height of water h(t, x) ≥ 0 and the mean velocity v(t, x) ∈ R, by the hyperbolic system

h t +(hv) x =0, (1.79) (hv) t + hv 2 + gh 2 2 x + ghZ (x)=0, (1.80)
where g denotes the gravity intensity and Z(x) is the bottom height, therefore h + Z is the level of the water surface (in what follows, we also denote the discharge by q = hv). We recall some properties of the Saint-Venant system. We take them into account later in order to develop a numerical method so as to be coherent with the physical model.

First of all, the system in naturally posed for h ≥ 0 and it is strictly hyperbolic for h>0. The water height h can indeed vanish (flooding zones,

Miko laj Szydlarski 35

CHIARA SIMEONI

Physical and Mathematical Modelling

dry soils, tidal flats); this facts leads to a theoretical and numerical difficulty, because the system loses hyperbolicity at h = 0. (See Section 1.2.2).

The Saint-Venant system admits a mathematical entropy, which is also the physical energy,

E(h, v, Z)= hv 2 2 + gh 2 2 + gZh, (1.81)
which satisfies the "entropy inequality"

E t + v E + gh 2 2 x ≤ 0. (1.82)
Another fundamental property is related to the stationary solutions, that are characterized by the relations

hv = C 1 , v 2 2 + g(h + Z)=C 2 .
(1.83)

Where C 1 and C 2 are two arbitrary constant. The equations (1.83) leads to several cases of special interest:

1. In the case of stationary state there is no motion, v = 0, and the water surface is flat, h + Z = C st , for some constant C st . This situation is refereed to as the lake at rest.

2.

A quasi-stationary state is generated by slightly perturbing the height of the stationary state. In the one-dimensional case, assume a bottom topography in the shape of a bump and a perturbation of the form

h = C st -Z + ,w h e
h(x, 0) = h l if x<0, h r if x>0, u(x, 0) = 0, (2.1)
where h l >h r ≥ 0. This is a special case of Riemann problem for which v l = v r = 0, and it is called dam-break problem because it models what happens if a dam separating two levels of water bursts at time t = 0. This is the shallow water equivalent of the shock-tube problem of gas dynamics.

We assume h r > 0. in each wedge where q = hv is constant (see Section 1.2.2).

Miko laj Szydlarski

In Figure (2.3a) we see that the 1-characteristics behave near the 1rarefaction wave just as we would expect from the nonlinear scalar case.

Miko laj Szydlarski 39

CHIARA SIMEONI again as we would expect from the scalar theory of hyperbolic conservation laws (refer to [START_REF] Lax | Hyperbolic systems of conservation laws and the mathematical theory of shock laws[END_REF]). These characteristics cross the 1-rarefaction with a smooth change in velocity. This is the standard situation for many nonlinear systems of equations.

Riemann Problem for Shallow Water Equations

For a system of m equations, there will be m characteristic families and m waves in the solution to the Riemann problem. If the p th wave is a shock, then characteristics of families 1 through p -1 will cross the shock from left to right, characteristics of family p+1 through m will cross from right to left, and characteristics of family p will impinge on the shock from both sides. This classical situation is observed in many physical problems, including the Euler equations of gas dynamics, and is the case that is best understood mathematically.

The shallow water equations are strictly hyperbolic and genuinely non-

Miko laj Szydlarski

linear (see Section 1.2.2) and so the Riemann problem always consists of two waves, each of which is a shock or rarefaction. The following example shows that other combinations are possible.

Consider the Riemann data To understand the characteristic structure of the shallow water equations, it is useful to consider what happens in the solution to the Riemann problems discussed above in the case where the initial jump is so small that the linearized equation gives a good model.

h(x, 0) ≡ h 0 ,v (x, 0) = v l if x<0, -v l if x>0, (2.2
Consider the data (2.1), for example, with h l = h 0 + and h r = h 0for some h 0 . Then if we solve the Riemann problem for the linearized equation with v 0 =0i n( 1.75), we find that the solution consists of two acoustic waves with speeds ± √ gh 0 , separated by a state (h m ,v m )w i t h

h m = h 0 ,v m = √ gh 0 (2.3)
The solution consists of two discontinuities. If we solved the nonlinear equations with this same data, the solution would look quite similar, but the

Miko laj Szydlarski 42

CHIARA SIMEONI

Riemann Problem for Shallow Water Equations

Shock Waves and Hugoniot Loci

We consider a shallow water 2-shock such as the right-going shock in section (2.1). This shock connects some state q m =(h m ,h m v m ) to the right state q r =(h r ,h r v r ) from the Riemann data. We will view q r as being fixed and determine all possible states q =(h, hv) that can be connected to q r by a 2-shock. We will find that there is a one-parameter family of such states, which trace out a curve in state space as shown in Figure (2.6b). Here the state space (phase plane) is the h -hv plane. This set of states is called a Hugoniot locus. Which one of these possible states corresponds to q m in the solution to the Riemann problem depends not only on q r but also on q l . The state q m must lie on the curve shown in Figure (2.6b), but it must also lie on an analogous curve of all states that can be connected to q l by a 1-wave, as determined below.

We now consider the problem of determining all states q that can be connected to some fixed state q * =(h * ,h * v * )(representingeitherq l or q r )by a shock. Across any shock the Rankine-Hugoniot condition must be satisfied, so for the shallow water equation, this gives a system of two equations that must simultaneously be satisfied:

s(h * -h)=h * v * -hv s(h * v * -hv)=h * v 2 * -hv 2 + 1 2 g(h 2 * -h 2 ), (2.4)
where s is the speed of the shock. Recall that q * is fixed and we wish to find all states q and corresponding speeds s satisfying these relation. We thus have two equations with three unknowns, so we expect to find a oneparameter family of solutions. In fact there are two distinct families of solutions, corresponding to 1-shocks and 2-shocks. For the time being we use the term "shock" to refer to a discontinuous (weak) solution satisfying the

Miko laj Szydlarski 44

CHIARA SIMEONI Hugoniot locus of points in shallow water state space that can be connected to a given state q r by a 2-shock satisfying the Rankine-Hugoniot conditions (2.4)( From Ref. [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] )

Rankine-Hugoniot condition. Later we will consider the additional admissibility condition that is required to ensure that such a solution is truly a physical shock wave. There are many different ways to parameterize these families. Fairly simple formulas result from using h as the parameter. For each value of h we will determine the corresponding v and s, and plotting hv against h will give the curves shown in Figure (2.6). We first determine v by eliminating s from the system (2.4). The first equation gives

s = h * v * -hv h * -h, (2.5)
and substituting this into the second equation gives an equation relating v to h. This is a quadratic equation in v that, after simplifying somewhat, becomes

v 2 -2v * v + v 2 * - g 2 h * h - h h * (h * -h) with roots v(h)=v * ± g 2 h * h - h h * (h * -h). (2.6)
Note that when h = h * this reduces to v = v * , as we expect, since the curves we seek must pass through the point q * . For each h = h * there are two different values of v, corresponding to the two families of shocks. In the case of a very weak shock (q ≈ q * ) we expect the linearized theory to hold, and so we expect one of these curves to be tangent to the eigenvector r 1 (q * ) at q * and the other to be tangent to r 2 (q * ). This allows us to distinguish which curve corresponds to the 1-shocks and which to 2-shocks. To see this more clearly, we multiply (2.6)b yh and reparameterize by a value α,w i t h

h = h * + α so that h = h * at α = 0, to obtain hv = h * v * + α v * ± gh * 1+ α h * 1+ α 2h * . (2.7)
Hence we have

q = q * + α 1 v * ± gh * + O(α) as α → 0 (2.8)
For α very small (as q approaches q * ), we can ignore the (α) term and we see that these curves approach the point q * tangent to the vectors 1

v * ± √ gh *
which are simply the eigenvectors of the Jacobian matrix (1.75) at q * . From this we see that choosing the -sign in (2.7) gives the locus of 1-shocks, while the + sign gives the locus of 2-shocks.

The Entropy condition

We consider a general Riemann problem with data q l and q r , and suppose we know that the solution consist of two shocks. We can then solve the Riemann problem by finding the state q m that can be connected to q l by 1-shock and also to q r by a 2-shock.

We found in the previous section that though the point q r there is a curve of points q that can be connected to q r by a 2-shock. For a shallow water

Miko laj Szydlarski

equations, these points satisfy (2.7) with the plus sign and with q * = q r .

Since q m must lie on this curve, we have

h m v m = h r v r +(h m -h r ) v r + gh r 1+ h m -h r h r 1+ h m -h r 2h r ,
which can be simplified to give

v m = v r +(h m -h r ) g 2 1 h m + 1 h r . (2.9)
Similarly, there is a curve through q l of states that can be connected to q l by a 1-shock, obtained by setting q * = q l and taking the minus sign in (2.7).

Since q m must lie on this curve, we find that

v m = v r -(h m -h r ) g 2 1 h m + 1 h r . (2.10)
We thus have a system of two equations (2.9) and (2.10) for the two unknowns h m and v m . Solving this system gives the desired intermediate state in the Riemann solution. We can easily eliminate v m from this system by noting that this appears only on the left of each equation, and the left-hand sides are equal, so equating the right-hand sides gives a single equation involving only the one unknown h m . This can be solved by an iterative method for nonlinear equations, such as the Newton's method.

What happens if we apply this same procedure to a Riemann problem where the physical solution should not consist of two shocks? For example, consider the dam-break Riemann problem from section (2.1), where the solution should consist of a 1-rarefaction and a 2-shock. We can still solve the problem in terms of two "shock waves" that satisfy the Rankine-Hugoniot jump conditions, as illustrated in Figure (2.7). This gives a weak solution of the conservation laws, but one that does not satisfy the proper entropy condition for this system, as discussed in the next section. The procedure for finding the physically correct solution with a rarefaction wave is given in Section (2.5). shows the structure for the weak solution for dam-break problem described at the end of previous subsection. Solution consists of two discontinuities.

The 1-characteristics are not impinging on the 1-shock as they should, an indication that this structure is not stable to small perturbations and that this shock should be replaced by a rarefaction wave. This suggests the following criterion for judging whether a given weak solution is in fact the physically correct solution, a generalization of the Lax Entropy Condition to systems of equations (refer to [START_REF] Lax | Hyperbolic systems of conservation laws and the mathematical theory of shock laws[END_REF]).

Lax Entropy Condition. A discontinuity separating states u

l = h l h l v l and u r = h r h r v r
,propagating at speed s, satisfies the Lax entropy condition if there is an index p such that: λ p (u l ) >s>λ p (u r ), so that p-characteristics are impinging on the discontinuity, while the other characteristics are crossing the discontinuity, λ j (u l ) <s and λ j (u r ) <s for j<p , λ j (u l ) >s and λ j (u r ) >s for j>p , (2.11)

In this definition we assume the eigenvalues are ordered so that λ 1 < λ 2 <...<λ m in each state.

Miko laj Szydlarski 48

CHIARA SIMEONI (From Ref. [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] )

For the shallow water equations there is a simple criterion that can be applied to determine which parts of each Hugoniot locus give physically correct shock waves satisfying the Lax entropy condition. Across a 1-shock connecting q l to a state q m , we require that the characteristic velocity λ 1 = v-√ gh must decrease. In conjunction with the Rankine-Hugoniot condition, it can be shown that this implies that h must increase, so we require h m >h l .

Similarly, a 2-shock connecting q m to q r satisfies the Lax entropy condition if h m >h r . Note from Figure (2.2) and Figure (2.4) that this also means that fluid particles experience an increase in depth as they pass through a shock. This is similar to the physical entropy condition for gas dynamics, that gas particles must experience an increase in physical entropy as they pass through a shock wave. Figure (2.7) shows the portions of each Hugoniot locus along which the entropy condition is satisfied as solid lines. These are simply the portions along which h is increasing. The portions indicated by dashed lines are states that can be connected by a discontinuity that satisfies the Rankine-Hugoniot condition, but not the entropy condition. We see from Figure (2.7b) that the solution to the dam-break Riemann problem consisting of two shocks fails to satisfy the entropy condition. Instead we must find a solution to the Riemann problem that consists of a 1-rarefaction and a 2-shock. In the next section we investigate rarefaction waves and

Miko laj Szydlarski

will see that the Hugoniot locus through q l must be replaced by a different curve, the integral curve of r 1 . The intersection of this curve with the 1shock Hugoniot locus will give the correct intermediate state (h m ,v m ) for the solution.

Simple Waves and Rarefaction

In this section we will investigate solutions that are smoothly varying but which also have the property that they are associated with only one characteristic family of the system. Such waves are called simple waves.

In particular, the centered rarefaction waves that arise in the solution to

Riemann problems for nonlinear systems are simple waves, but these are just one special case. They are special in that they also have the property that they are similarity solutions of the equations and are constant along every ray x t = constant. They arise naturally from Riemann problems because of the special data used, which varies only at a single point x = 0, and hence all variation in the solution flows out from the point x = t = 0. Let q(ξ) be a smooth curve through state space parameterized by a scalar parameter ξ. We say that this curve is an integral curve of the vector field r p if at each point q(ξ) the tangent vector to the curve, q (ξ), is an eigenvector of f (q(ξ)) corresponding to the eigenvalue λ p (q(ξ)). If we have chosen some particular set of eigenvectors, that we call r p (q), then q (ξ) must be some scalar multiple of the particular eigenvector r p (q(ξ)), q (ξ)=α(ξ)r p (q(ξ)),

(2.12)

The value of α(ξ) depends on the particular parameterization of the curve and on the normalization of r p , but the crucial idea is that the tangent to the curve is always in the direction of the appropriate eigenvector r p evaluated at the point on the curve. Figure (2.9) shows integral curves of r 1 and r 2

for the shallow water equations, for which the eigenvectors are given by (1.77). As an example of how these curves can be determined, consider r 1 and set α(ξ) ≡ 1, which selects one particular parameterization for which

Miko laj Szydlarski 50

CHIARA SIMEONI the formulas are relatively simple. Then (2.12)r e d u c e st o q (ξ)=r 1 (q(ξ)) = 1 q2 q1 -√ g q1 (2.13) by using (1.77). This gives two ordinary differential equations for the two components of q(ξ);

Riemann Problem for Shallow Water Equations

(q 1 ) = In terms of the velocity instead of momentum, we can rewrite this as

v = v * +2 gh * - gh .
(2.20)

Similarly, the integral curve of r 2 passing though the point q * can be shown to have the form

v = v * -2 gh * - gh .
(2.21)

The expression (2.20) describes an integral curve r 1 ,w h e r eq * is an arbitrary point on the curve. This can be rewritten as

u +2 gh = v * +2 gh * .
Since q * and q are any two points on the curve, we see that the function

w 1 (q)=v +2 gh (2.22)
has the same value at all points on this curve. This function is called Riemann invariant for the 1-family, or simply a 1-Riemann invariant. It is a function of q whose value is invariant along any integral curve of r 1 , though it will take a different value on a different integral curve.

Similarly, from (2.22) we see that

w 2 (q)=v -2 gh, (2.23)
is a 2-Riemann invariant, a function whose value is constant along any itegral curve of r 2 .

A simple wave is a special solution to the conservation law in which q(x, t)=q(ξ(x, t)), (2.24)

Miko laj Szydlarski 52

CHIARA SIMEONI

where q(ξ) traces out an integral curve of some family of eigenvectors r p and ξ(x, t) is a smooth mapping from (x, t) to the parameter ξ. This means that all states q(x, t) appearing in the simple wave lie on the same integral curve.

Note that any p-Riemann invariant is constant throughout the simple wave.

But not every function of the form (2.24) will satisfy the conservation law.

The function ξ(x, t) must be chosen appropriately. We compute q t =q (ξ(x, t))ξ t and q x =q (ξ(x, t))ξ x , so to satisfy q t + f (q) x =0w em ustha v e ξ t q (ξ)+ξ x f (q(ξ))q (ξ)=0.

Since q (ξ) is always an eigenvector of f (q(ξ)), this yields

[ξ t + ξ x λ p (q(ξ))]q (ξ)=0,
and hence the function ξ(x, t) must satisfy

ξ t + λ p (q(ξ))ξ x =0. (2.25)
Note that this is a scalar quasilinear hyperbolic equation for ξ.

In particular, if we use initial data q(x, 0) that is restricted entirely to this integral curve, so that q(x, 0) = q( ξ(x)) for some smooth choice of ξ(x), then (2.24) will be a solution to the conservation law for t>0 provided that ξ(x, t) solves (2.25) with initial data ξ(x, 0) = ξ(x), at least for as long as the function ξ(x, t) remains smooth. In a simple wave the nonlinear system of equations reduces to the scalar nonlinear equation (2.25) for ξ(x, t).

A centered rarefaction wave is a special case of a simple wave, in which ξ(x, t)= x t , so that the solution is constant on rays through the origin. A centered rarefaction wave has the form

q(x, t)=        q l if x t ≤ ξ 1 , q x t if ξ 1 ≤ x t ≤ ξ 2 , q r if x t ≥ ξ 2 ,
(2.26)

where q l and q r are two points on a single integral curve with λ p (q l ) < λ p q r . This condition is required so that characteristics spread out as time advances

Miko laj Szydlarski 53

CHIARA SIMEONI

Riemann Problem for Shallow Water Equations

and the rarefaction wave makes physical sense. For a centered rarefaction wave a particular parameterization of the integral curve is forced by the fact that we set ξ = x t . Rewriting this as x = ξt, we see that the value q(ξ) observed along the ray x t is propagating at speed ξ , which suggests that ξ at each point on the integral curve must be equal to the characteristic speed λ p (q(ξ)) at this point. This is confirmed by noting that (2.25) in this case becomes

x t = λ p q( x t ) .
(2.27)

In particular, the left edge of the rarefaction fan should be the ray x t = λ p (q l ) so that ξ 1 = λ p (q l )i n( 2.26), while the right edge should be the ray x t = λ p (q r ) so that ξ 2 = λ p (q r ). We thus have

ξ 1 = λ p (q l ), q(ξ 1 )=q l , ξ 2 = λ p (q r ), q(ξ 2 )=q r ,
(2.28)

To determine how q(ξ) varies for ξ 1 < ξ < ξ 2 through the rarefaction wave

(2.26), rewrite (2.27) as

ξ = λ p (q(ξ)) (2.29)
and differentiate this with respect to ξ to obtain 1=∇λ p (q(ξ)) • q (ξ) (2.30) Using (2.12)i n( 2.30) gives 1=α(ξ)∇λ p (q(ξ)) • r p (q(ξ))

and hence

α(ξ)= 1 ∇λ p (q(ξ)) • r p (q(ξ)) . (2.31)
Using this in (2.12) gives a system of ODEs for q(ξ):

q (ξ)= r p (q(ξ)) ∇λ p (q(ξ)) • r p (q(ξ)) (2.32)
This system must be solved over the interval ξ 1 ≤ ξ ≤ ξ 2 using either of the conditions in (2.28) as an initial condition. Note that the denominator is nonzero provided that λ p is monotonically varying. A rarefaction wave would not make sense past a point where the denominator vanishes.

Finally, for the shallow water equations (1.73)w eh a v e

λ 1 = v - √ gh = q 2 q 1 - √ gq 1 , ∇λ 1 =   -q 2 (q 1 ) 2 -1 2 g q 1 1 q 1   , r 1 = 1 q 2 q 1 - √ gq 1 ,
(2.33) and hence

∇λ 1 • r 1 = - 3 2 g q 1 ,
(2.34) so that the equation of this system is

q (ξ)=- 2 3 h(ξ) g . (2.35) The general solution is h = 1 9 g(A -ξ) 2 , (2.36)
for some constant A. This constant must be chosen so that (2.28) is satisfied, i.e., so that h = h l at ξ = v l -√ gh l and also h = h r at ξ = v r -√ gh r .

Provided that q l and q r both lie on an integral curve of r 1 , as they must if they can be joined by a centered rarefaction wave, we can satisfy both of these condition by taking

A = v l +2 gh l = v r +2 gh r .
(2.37)

Recall that v +2 √ gh is a 1-Riemann invariant, which has the same value at all points on the integral curve. We see that h varies quadratically with ξ = x t through a rarefaction wave (2.26).

Once we know h as a function of ξ , we can use the formula (2.21), which holds through any simple wave, to determine how v varies through the rarefaction wave. (i.e., we use the fact that the Riemann invariant is constant). Note that since we know the relation between h and v from having previously found the Riemann invariants, we do not need to solve both the ODEs in the system (2.35). We have chosen the simpler one to solve. This trick is often useful for other systems as well.

Miko laj Szydlarski

Solving the Dam-Break Problem

We start by solving the Riemann problem for which we know the solution consist of two rarefaction waves. So lets consider the shallow water equations with data (2.2), but now take v l < 0. This corresponds to two streams of water that are moving apart from one another. Again the solution will be symmetric but will consist of two rarefaction waves. To solve this specific Riemann problem, we can proceed in a manner similar to what we did in Section (2.3) for the all-shock solution. There is an integral curve of r 1 through q l consisting of all states that can be connected to q l by a 1rarefaction, and an integral curve of r 2 through q r consisting of all states that can be connected to q r by a 2-rarefaction. These are illustrated in Figure The intermediate state q m in the Riemann solution must lie on both of these curves, and hence is at the intersection as shown in Figure (2.10a). For this particular example q m lies on the h-axis due to symmetry. In general we can find the intersection by using the fact that q m must lie on the curve Miko laj Szydlarski 56 described by (2.20)w i t hq * = q l and on the curve described by (2.21)w i t h

q * = q r ,s o v m = v l +2 √ gh l - √ gh m , v m = v r -2 √ gh r - √ gh m , (2.39)
This is s system of two nonlinear equations for h m and v m . Equating the right-hand sides gives a single equation for h m , which can be explicitly solved to obtain

h m = 1 16 g v l -v r +2 gh l + gh r (2.40)
This is valid provided that the expression being squared is nonnegative.

When it reaches zero. the outflow is sufficiently great that the depth h m goes to zero. For a given state q l only some points on the integral curve of r 1 can be connected to q l by a rarefaction wave that makes physical sense, since we are assuming q l is the state on the left of the rarefaction wave. We must have λ 1 (q l ) < λ 1 (q) for all states q in the rarefaction wave, and hence q must lie on the portion of the integral curve. Similarly, if q r is the state to the right of a 2-rarefaction, then states q in the rarefaction must satisfy λ 2 (q) < λ 2 (q r ) and must lie on the solid portion of the integral curve for r 2 sketched through q r in Figure (2.10a). For the data shown in this figure, there is a state q m that can be connected to both q l and q r by physically correct rarefaction waves, and the Riemann solution consists of two rarefactions.

For other data this might not be the case. Figure (2.10b) shows the data for the dam-break Riemann problem, with h l = 3, h r = 1, and v l = v r = 0.

We can still use (2.40) to compute an intermediate state q m lying at the intersection of the integral curves, as illustrated in Figure (2.10b), but the resulting 2-wave does not make physical sense as a rarefaction wave, since this problem that consists of two shock waves,one of which does not satisfy the Lax entropy condition. To find the correct solution we must determine an intermediate state q m that is connected to q l by a 1-rarefaction wave and simultaneously is connected to q r by a 2-shock wave. The state q m must lie on an integral curve of r 1 passing through q l ,s ob y( 2.20)w em u s th a v e

λ 2 (q m ) > λ 2 (q r ). Compare Figure (2.
v m = v l +2 gh l - gh m . (2.41)
It must also lie on the Hugoniot locus of 2-shocks passing though q r ,s ob y (2.9) it must satisfy

v m = v r +(h m -h r ) g 2 1 h m + 1 h r . (2.42)
We can easily eliminate v m from these two equation and obtain a single connecting q l to q m can then be determined using theory of Centered Rarefaction Waves. Note that the intermediate state q m resulting from this procedure will be slightly different from that obtained in either Figure

The General Riemann Solver for Shallow Water equations

For the dam-break problem we know that the 1-wave is a rarefaction while the 2-wave is a shock, leading to the system of equations (2.41) and

(2.42) to solve for h m and v m . For general values of q l and q r we might have any combination of shocks and rarefactions in the two families, depending on the specific data. In general, to find the state q m we can define two functions φ l and φ r by

φ l (h)=      v l +2 √ gh l - √ gh if h<h l , v l -(h -h l ) g 2 - 1 h -1 h l if h>h l , and
φ r (h)=      v r -2 √ gh r - √ gh if h<h r , v r +(h -h r ) g 2 - 1 h + 1 hr if h>h r ,
For a given state h, the function φ l (h) returns the value of v such that (h, hu) can be connected to q l through a physically correct 1-wave, while φ r (h) returns the value such that (h, hv) can be connected to q r through a physically-correct 2-wave. We want to determine h m so that φ l (h m )= φ r (h m ). This can be accomplished by applying a nonlinear root finder to the function φ(h) ≡ φ l (h) -φ r (h), such as the Newton method.

In a scalar equation, such a Burgers equation, when two shock waves collide, they simply merge into a single shock wave with a larger jump.

Miko laj Szydlarski 59

CHIARA SIMEONI Figure 2.12: (a) Initial states q 1 and q 3 are each connected to q 2 by 2-shock.

After collision, solving the Riemann problem between q 3 and q 1 gives a new state q 4 and a reflected 1-wave. (From Ref. [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] )

For a system of equations, the result of a shock collision is not so simple, even if the two shocks are in the same characteristic family. The result will include a stronger shock of this same family, but the collision will typically also introduce waves of the other families. We consider initial data for the shallow water equations (1.71) consisting of the three states shown in Figure (2.12a). These three states all lie on the Hugoniot locus S 2 (q 2 ), so there is a 2-shock connecting q 1 to q 2 and a slower 2-shock connecting q 2 to q 3 .I fw e solve the shallow water equations with data

q(x, 0) =        q 1 if x<x 1 , q 2 if x 1 ≤ x ≤ x 2 , q 3 if x>x 2 , (2.43)
for some initial shock locations x 1 <x 2 , then the solution consists of the these two shocks, which eventually collide at some point x c .A tt h et i m et c when they collide, the state q 2 disappears and the solution has the form

q(x, t c )= q 1 if x<x c , q 3 if x>x c , (2.44)
To determine the solution beyond this time, note that this has the form of Riemann problem data with left state q 1 and right state q 3 . The Riemann solution is not a single 2-shock, because q 1 will not lie on the Hugoniot locus Miko laj Szydlarski 60 S 2 (q 3 ). Instead, a 1-wave must be introduced to connect q 1 to a new state q 4 that lies on this Hugoniot locus, as illustrated in Figure (2.12b). We see that the 1-wave must be a rarefaction wave, since h 4 <h 1 and h 4 is determined by the intersection of the integral curve R 1 (q 1 ) with the Hugoniot locus S 2 (q 3 ).

Miko laj Szydlarski 61

h = R f (h, ξ -v)dξ, (3.3
)

hv = R ξf (h, ξ -v)dξ, (3.4
)

hv 2 + gh 2 2 = R ξ 2 f (h, ξ -v)dξ. (3.5)
The two Saint-Venant equations are obtained by taking the moments of the kinetic equation (3.1)i ndξ against 1, ξ and ξ 2 repectively: the righthand side vanishes by (3.2) and left-hand sides coincide exactly thanks to hypothesis (3.3)- (3.5). Therefore, the non-linear shallow water system can be viewed as a single linear equation on a non-linear quantity f , for which it is easier to find simple numerical schemes with good theoretical properties. In particular, any numerical method known for the simple advection equations can directly be applied to (3.1), provided that an explicit form of the kinetic density function f is introduced.

For the numerical scheme used here, the density of particles f (t, x, ξ)i s given by a so-called Gibbs equilibrium

f (t, x, ξ)=f (h, ξ -v)= h(t, x)χ ξ -v(x, t) h(t, x) , (3.6)
with χ defined by

χ(ω)= √ 2 π √ g 1 - ω 2 2g 1 2 + . (3.7)
We also recover the energy by

E(h, v, Z)= R ξ 2 2 f (ξ)+ π 2 g 2 6 f 3 (ξ)+gZf(ξ) dξ. (3.8)

Miko laj Szydlarski 63

CHIARA SIMEONI

Numerical Simulation of the Shallow Water equations

We finally point out that that, choice of χ(x)i n( 3.7) is the only function such that the density (3.6) satisfies the equation

ξ ∂f ∂x -gZ (x) ∂f ∂ξ = 0 (3.9)
on the steady state of the lake at rest

v(t, x)=0,h (t, x)+Z(x)=H, ∀t ≥ 0.
This allows to construct a consistent numerical scheme for the Saint-Venant equations with the good properties of preserving the stationary solutions as well as being endowed with a discrete entropy inequality (see Section 3.3)

Finite Volume Methods for Conservation Laws

Before we present the numerical scheme for the one-dimensional Saint-Venant system, based on the kinetic approach described in Section 3.1,i ti s necessary to introduce some basic concept of numerical methods for conservation laws. For the most of practical interest, it is not possible to compute an exact solutions of the shallow water equations, by using analytical techniques such as the method of characteristics. This has lead to develop to numerical methods where the continuous problem, i.e. the governing equations, is transformed into a discrete form which then results in a series of algebraic equations that can be solved on a computer. The solution to the discrete problem represents an approximation to the solution of the continuous problem and various concepts have been developed in an attempt to quantify how good the calculated numerical solution compared to the analytic solution. The following general description of basic concepts for numerical analysis of conservation laws are taken from LeVeque [START_REF] Leveque | Numerical methods for conservation laws[END_REF] and [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF].

In one space dimension, a finite volume method is based on subdividing the spatial domain into intervals (the "finite volumes", also called grid cells) and keeping track of an approximation to the integral of the conserved quantities u of he model, over each of these volumes. In each time step we update these values using approximations to the flux through the endpoints of the intervals. Denote the i th grid cell by

C i =(x i-1 2 ,x i+ 1 2 ),
Miko laj Szydlarski 64 ). The value U n i will be an approximation of the average value of u over i th interval at time t n :

U n i ≈ 1 ∆x x i + 1 2 x i -1 2 u(x, t n )dx ≡ C i u(x, t n )dx, (3.10)
where ∆x = x i+ 1 2

-x i-1 2
is the length of the cell. For simplicity we will generally assume a uniform grid, but is not required.

If u(x, t) is a smooth function, by applying Taylor's theory, the integral in (3.10) agrees with the value of u at the midpoint of the interval to O(∆x 2 ).

By working with cell averages, however, it is easier to use important properties of the conservation law in deriving numerical methods. In particular, we can insure that the numerical method is conservative in a way that mimics the true solution, and this is extremely important in accurately calculating shock waves. This is because

N i=1
U n i ∆x approximates the integral of u over the entire domain of computation, and if we use a method that is in conservation form (as described below), then this discrete sum will change only due to fluxes at the boundaries of that domain. The total mass within the computational domain will be preserved, or at least will vary correctly provided the boundary conditions are properly imposed. The integral form

Miko laj Szydlarski

of the conservation law (1.23) gives:

d dt C i u(x, t)dx = f (u(x 1-1 2 ,t)) -f (u(x 1+ 1 2 ,t)).
(3.11)

We can use this expression to develop an explicit time-marching algorithm.

Given U n i , the cell averages at time t n , we want to approximate U n+1 i ,t h e cell averages at the next time t n+1 after a time step of length ∆t = t n+1 -t n .

Integrating (3.11) in time from t n to t n+1 yields

C i u(x, t n+1 )dx - C i u(x, t n )dx = t n+1 tn f (u(x 1-1 2 ,t))dt - t n+1 tn f (u(x 1+ 1 2 ,t))dt.
Rearranging this and dividing by ∆x gives

1 ∆x C i u(x, t n+1 )dx = 1 ∆x C i u(x, t n )dx -1 ∆x t n+1 tn f (u(x 1+ 1 2 ,t))dt - t n+1 tn f (u(x 1-1 2 ,t))dt . (3.12)
This tells us exactly how the cell average of u from (3.10) should be updated in one time step. In general, however, we cannot evaluate the time integrals on the right-hand side of (3.12) exactly, since u(x i± 1 2 ,t) varies with time along each edge of the cell, and we do not have the exact solution to work with. This does suggest that we should study numerical methods of the form

U n+1 i = U n i - ∆t ∆x (F n i+ 1 2 -F n i-1 2 ), (3.13)
where

F n i-1 2
is some approximation to average flux along x = x i-1 2 : For a hyperbolic problem, information propagates with finite speed, so it is reasonable to first suppose that we can obtain

F n i-1 2 ≈ 1 ∆t t n+1 tn f (u(x 1-1 2 ,t))dt. ( 3
F n i-1 2
based only on the

Miko laj Szydlarski

values U n i-1 and U n i , the cell averages on either side of this interface. Then we might use a formula of the form

F n i-1 2 = F(U n i-1 ,U n i ) (3.15)
where F is some numerical flux function.T h em e t h od( 3.13) then becomes

U n+1 i = U n i - ∆t ∆x [F(U n i ,U n i+1 ) -F(U n i-1 ,U n i )] (3.16)
The specific method obtained depends on how we choose the formula F, but in general any method of this type is an explicit method with a threepoint stencil, meaning that the value U n+1 i will depend on the three values

U n i-1 ,U n i and U n i+1
at the previous time level. Moreover, it is said to be in conservation form, since it mimics the property (3.12) of the exact solution.

Note that if we sum ∆xU n+1 i from (3.13) over any set of cells, we obtain

∆x N i=1 U n+1 i = ∆x N i=I U n i - ∆t ∆x ( F n N + 1 2 -F n 1- 1 2
).

(3.17 The method (3.16) can be viewed as a direct finite difference approximation to the conservation law u t + f (u) x = 0, since it gives

U n+1 i -U n i ∆x + F n i+ 1 2 -F n i-1 2 ∆x =0. ( 3

.18)

Many methods can be equally viewed as finite difference approximations to this equation or as finite volume methods.

There are several considerations that go into judging how good a particular flux function is for numerical computation. One essential requirement is that the resulting method should be convergent, i.e., the numerical solution should converge to the true solution of the differential equation as the grid is refined (as ∆x, ∆t → 0). This generally requires two conditions:

• the method must be consistent with the differential equation, meaning that it approximates it well locally;

Miko laj Szydlarski

• the method must be stable in some appropriate sense, meaning that the small errors made in each time step do not grow too fast in later time steps.

The numerical flux should approximate the integral in (3.14). In particular, if the function u(x, t) ≡ u is constant in x,thenu will not change with time and the integral in (3.14)s i m p l yr e d u c e st of (u). As a result, if U n i-1 = U n i = u, then we expect the numerical flux function F of (3.15)toreducetof (u), so we require

F(u, u)=f (u) (3.19)
for any value u. This is part of basic consistency condition. We generally also expect continuity in this function as U i-1 and U i vary, so that

F(U i-1 ,U i ) → f (u) as U i-1 ,U i → u.
Typically some requirement of Lipschitz continuity is made, e.g., there exist a constant L so that

|F(U i-1 ,U i ) -f (u)| ≤ L max(|U i -u|, |U i-1 -u|) (3.20)
The so called CFL condition, must be satisfied by any finite volume or finite difference method if we expect it to be stable and converge to the solution of the differential equation as the grid is refined. It simply states that the method must be used in such a way that information has a chance to propagate at the correct physical speeds, as determined by the eigenvalues of the flux Jacobian f (u).

With the explicit method (3.16) the value U n i+1 depends only on three values U n i-1 ,U n i , and U n i+1 at the previous time step. Suppose we apply such a method to the advection equation u t + au x =0w i t ha>0 so that the exact solution simply translates at speed a and propagates a distance a∆t over one time step.

Friedrichs and Lewy [2]:

A numerical method can be convergent only if its numerical domain of dependence contains the true domain of dependence of the differential equation, at least in the limit as ∆t and ∆x go to zero.

It is very important to note that the CFL condition is only a necessary condition for stability (and hence convergence). It is not always sufficient to guarantee stability.

The numerical domain of dependence of a method is defined as the set of points where the initial data can possibly affect the numerical solution at any point (x, t).

In order for the CFL condition to be satisfied, the domain of dependence of the true solution must lie within computational domain. This leads to a relation between the time step ∆t and space step ∆x through the Curantnumber, to be checked in the implementation of the scheme. For a three-point method the CFL condition leads to a necessary condition ν ≤ 1. Note that if the method has a wider stencil, then the CFL condition will lead to a more lenient condition on the time step.

For hyperbolic equations we typically use explicit methods and grids for which the Courant number is somewhat smaller than 1. This allows keeping ∆t ∆x fixed as the grid is refined, which is sensible in that generally we wish to add more resolution at the same rate in both space and in time in order to improve the solution.

Kinetic scheme for the Saint-Venant System

We introduce a finite volume scheme for the one-dimensional Saint-Venant system, based on the kinetic approach described in Section 3.1, that is taken from Perthame, Simeoni [START_REF] Perthame | Kinetic scheme for the Saint-Venant system with a source term[END_REF].

If Z(x) is the function describing the bottom height, the cell averages are given by Z i = 1 we perform a discretization directly on the density of particles

f n+1 i (ξ) -f n i (ξ)+ ∆t ∆x ξ f - i+ 1 2 (ξ) -f + i-1 2 (ξ) =0, (3.22)
where the interface densities f ± i+ 1 2 are defined later. As usual, the "collision term" Q(x, t, ξ) in the kinetic representation, which relaxes the kinetic density f to a Gibbs equilibrium (3.6), is neglected in the numerical scheme; at each time-step we project f n i (ξ) on the equilibrium, which is way to performe all collisions at once and to recover the Gibbs equilibrium without computing it.

Note that the fluxes can also be written as

f - i+ 1 2 (ξ)=f i+ 1 2 (ξ)+ f - i+ 1 2 (ξ) -f i+ 1 2 (ξ) (3.23)
and the quantity δf - (ξ) = 0 for ξ ≥ 0 in the scheme below. This is the principle of interfracial upwind sources: the source is not treated as a volumic term but at the interfaces and it is upwinded. Now, we integrate Eq. (3.22)i ndξ against 1 and ξ, with notation

i+ 1 2 (ξ)=f - i+ 1 2 (ξ) -f i+ 1 2 (
U n+1 i =(h n+1 i , (hu) n+1 i
),

h n+1 i = R f n+1 i (ξ)dξ, (hu) n+1 i = R ξf n+1 i (ξ)dξ
and obtain the macroscopic scheme

U n+1 i -U n i + ∆t ∆x F - i+ 1 2 -F + i-1 2 =0. (3.24)
The numerical fluxes are thus given by the kinetic fluxes

F - i+ 1 2 = R ξ 1 ξ f - i+ 1 2 (ξ)dξ, (3.25)
F + i-1 2 = R ξ 1 ξ f + i-1 2 (ξ)dξ. (3.26)
In order to take the neighboring cells into account in a natural interpretation of the microscopic features of the system, we formulate a peculiar discretization for the fluxes in (3.22), computed by the upwind formulas

f - i+ 1 2 (ξ)=f n i (ξ)I ξ≥0 + f n i+ 1 2 (ξ)I ξ≤0 , (3.27)
f + i-1 2 (ξ)=f n i-1 2 (ξ)I ξ≥0 + f n i (ξ)I ξ≤0 , (3.28)
where we define The conservation of water height and momentum is also obvious for the system with flat bottom: the continuous system (1.79)-(1.80) becomes homogeneous (Z (x) = 0) and we obtain a conservative scheme, with the fluxsplitting form of the standard kinetic scheme.

f n i+ 1 2 (ξ)=f n i (-ξ)I |ξ| 2 ≤2g∆Z i+ 1 2 + f n i+ 1 2 (- |ξ| 2 -2g∆Z i+ 1
Remark: We emphasize that demonstrating the numerical scheme (3.24)-(3.26) to be consistent can not be achieved in the classical manner. Because of the presence of the source term and the decision to process it implicitly.

We state here theoretical properties of the numerical scheme introduced above, which represent the discrete analogues of the main properties of the Saint-Venant system stated in Section 1.2.3.W er e f e rt o [START_REF] Perthame | Kinetic scheme for the Saint-Venant system with a source term[END_REF] for the proof.

We assume the CFL condition, and the discrete energy

∆t
E n i = h n i |v n i | 2 2 + g(h n i ) 2 2 + gZ i h n i ;
3. the scheme (3.24)-(3.26) preserves the steady states of the system given by a lake at rest,

v n i =0,h n i + Z i = H, ∀i ∈ Z, ∀n ∈ N.
To proceed to the actual implementation of the scheme (3.24)-(3.26), we have to compute the numerical fluxes explicitly. Since their expressions are not always immediate to calculate, it is necessary to use an approximation technique for some of them. We refer to to [START_REF] Perthame | Kinetic scheme for the Saint-Venant system with a source term[END_REF] for the details.

Numerical Tests

We conclude this chapter with numerical examples that illustrate the results stated in the previous sections and we check the properties of the scheme on different test cases for which analytical solutions of the equations are available (see Chapter 2). For each test, the channel length is L = 200m and the computational domain is chosen to be symmetric around the point x = 0, the mesh size is ∆s = L/100 and the time-step ∆t is computed according to the CFL condition (3.29).

Miko laj Szydlarski 74

Small perturbation

We now consider a test case concerning the small perturbation on water surface, over flat bottom. This situation is related to discussion about linearization for shallow water equations form Section (1.2.2). Figures (3.

Introduction 1

 1 Physical and Mathematical Modelling 1.1 Physical basis: Open -channel flows .............. 1.1.1 Flow classification by depth variation ......... 1.1.2 Flow classification by Froude Number, Surface wave speed, Specific Energy and Critical Depth ....... 1.1.3 Frictionless Flow over a Bump and Hydraulic Jumps .

1. 2

 2 Systems of conservation laws .................. 1.2.1 Euler equations and Isentropic flows .......... 1.2.2 Shallow water equations ................. 1.2.3 The Saint-Venant System ................ 2 Riemann Problem for Shallow Water Equations 2.1 Dam-Break and Riemann Problems ............... 2.2 Shock Waves and Hugoniot Loci ................ 2.3 The Entropy condition ...................... 2.4 Simple Waves and Rarefaction ................. 2.5 Solving the Dam-Break Problem ................ 2.6 The General Riemann Solver for Shallow Water equations .. 3 Numerical Simulation of the Shallow Water equations 3.1 Kinetic approach for the Saint-Venant System ......... 3.2 Finite Volume Methods for Conservation Laws ........ 3.3 Kinetic scheme for the Saint-Venant System .......... 3.4 Numerical Tests ..........................

1. 1 1

 11 Isovelocity contours in typical straight open channel flows .. 1.2 Geometry and notation for open-channel flow ......... 1.3 Open-channel flow classified by regions ............. 1.4 Analysis of small surface wave propagating into still shallow water ................................ 1.5 Flow under a sluice gate ..................... 1.6 Specific-energy considerations .................. 1.7 Frictionless two-dimensional flow ower a bump ........ 1.8 Naturally occurring hydraulic jump ............... 1.9 Classification of hydraulic jumps ................ 1.10 The level of water surface .................... 2.1 Solution of the dam-break Riemann problem for the shallow water equation .......................... 2.2 Structure of the similarity solution of the dam-break Riemann problem for the shallow water equation ............. 2.3 Solution of the dam-break Riemann problem for the shallow water equation show in xt plane ............... 2.4 Structure of the similarity solution of the two-shock Riemann problem for the shallow water equations ............ 2.5 Solution of the two-shock Riemann problem for the shallow water equations .......................... 2.6 Hugoniot locus of points q =( h, hv) in shallow water state space ................................ 2.7 All-shock solutions to the shallow water Riemann problem .. Physical basis: Open -channel flows

Figure 1 . 1 :Figure 1 . 2 :

 1112 Figure 1.1: Measured isovelocity contours in typical straight open-channel flows (From Ref. [9] )

  the liquid density is nearly constant, the steady-flow continuity equation reduces to constant-volume flow Q along the channelQ = V (x)A(x)=const (1.1)where V is average velocity and A the local cross-sectional area, as sketched in Fig.(1.2).

1. 1 . 2 Figure 1 . 3 :Fr < 1

 12131 Figure 1.3: Open-channel flow classified by regions of rapidly varyfing flow (RVF), gradually varyfing flow (GVF), and uniform-flow depth profiles. (From Ref. [9] )

Fig. ( 1 1 . 3 )Figure 1 . 4 : 4 ) 2

 1131442 Figure 1.4: Analysis of small surface wave propagating into still shallow water: (a) moving wave, nonsteady frame; (b) fixed wave, inertial frame of reference. (From Ref. [9] )

Figure 1 . 5 :

 15 Figure 1.5: Flow under a sluice gate accelerates from subcritical to critical to supercritical flow and then jumps back to subcritical flow.

Figure 1 . 6 :

 16 Figure 1.6: Specific-energy considerations: (a) illustration sketch; (b) depth versus E from Eq. (1.8), showing minimum specific energy occurring at critical depth. (From Ref. [9] )

( 1 . 5 )

 15 for reference. It is an important parameter in characterizing open-channel flow. The associated energy is:

  .12) For E<E min solution does not exists in Fig.(1.6), and thus such a flow is impossible physically. For E>E min solutions are possible: (1) large depth with V<V c , called subcritical, and (2) small depth with V>V c , called supercritical. In subcritical flow, disturbances can propagate upstream because wave speed c 0 >V. In supercritical flow, waves are swept downstream:

1. 1 . 3 Frictionless

 13 Flow over a Bump and Hydraulic Jumps Now, let us consider a open-channel flow over a bump, as in Fig. (1.7). The behavior of the free surface is sharply different according to whether the approach flow is subcritical or supercritical. The height of the bump also can change the character of the results. For frictionless steady one-dimensional flow, sections 1 and 2 in Fig. (1.7a) are related by continuity and momentum balance:

Figure 1 . 7 :

 17 Figure 1.7: Frictionless two-dimensional flow ower a bump: (a) definition sketch showing Froude-number dependence; (b) specific-energy plot showing bump size and water depths. (From Ref. [9] )

Figure 1 . 8 :Figure 1 . 9 :

 1819 Figure 1.8: Naturally occurring hydraulic jump observed on the Upper Spokane Falls north channel (USA). (Photography from Wikipedia)

  (1.1.2) A hydraulic jump is exactly equivalent to the strong fixed wave in Fig. (1.4), where the change in depth δy is not neglected. If V 1 and y 1 upstream are known, V 2 and y 2 are computed by applying continuity and momentum across the wave, as in Eqs. (1.3) and (1.4). Equation (1.5) is therefore the correct solution for a jump if we interpret C and y in Fig. (1.4) as upstream conditions V 1 and y 1 ,w i t hC -δV and y + δy being the downstream conditions V 2 and y 2 , as in Fig. (1.9) Equation (1.5) becomes:

  .70) We can combine equation (1.67), (1.70) into the system of one-

A

  lot of mathematical models have been developed to describe fluid flows, the most general is the Navier-Stokes equations, that are used to predict the behavior of viscous compressible/incompressible fluids in three dimensions. In practice when building a mathematical model, many assumptions are made to simplify the problem under consideration, and the most basic equations for describing the required phenomena are used. In open channel flows the most commonly used model is the shallow water system of equations, in which it is assumed that the flow is shallow with respect to the dimension of the considered framework. The basis of shallow water model is the continuity equation, corresponding to conservation of mass, and an equation of motion, as we have seen in Section (1.2.2). Additional terms may be incorporated to include other effects such as friction, geometry variation, viscosity etc. and these are referred to as the source terms which generally correspond to some form of loss or gain terms in the equations of the system. In the case of modeling predominantly one-dimensional flows, the Saint-Venant equations are the most commonly used system for solving open channel flow problems, and these describe the gradually varied flow (see Section 1.1.1) of an incompressible inviscid fluid. The equations consist of a continuity or mass equation, and an equation of motion which is formed by applying Newton's Second law of motion along the channel.

Figure 1 . 10 :

 110 Figure 1.10: The level of water surface. (From Ref. [1] )

Figure ( 2 . 1 )

 21 shows the evolution of the depth and fluid velocity for the dam-break problem with data h l = 3 and h r = 1. Figure (2.2) shows the structure of this solution in the x-t plane. Water flows from left to right in a wedge that expands from the dam location x = 0. At the right edge of this wedge, moving water with some intermediate depth h m and velocity v m > 0 slams into the stationary water with h = h r , accelerating it instantaneously through a shock wave. The water is accelerated away from the deeper stationary water through the structure of a centered rarefaction wave.

Figure 2 . 1 :Figure ( 2 . 3 )

 2123 Figure 2.1: Solution of the dam-break Riemann problem for the shallow water equation with v l = v r = 0. On the left is the depth h and on the right is the momentum hv.( From Ref. [5] )

Figure 2 . 2 :Figure 2 . 3 :

 2223 Figure 2.2: Structure of the similarity solution of the dam-break Riemann problem for the shallow water equation with v l = v r = 0. The depth h, velocity v, and vertically integrated pressure are displayed as function of xt . The structure in the x-t plane is also shown with particle path indicated for a set of particles with the spacing between particles inversely proportional to the depth. (From Ref.[START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] )

)

  If v l > 0, then this corresponds to two streams of water slamming into each other, with the resulting solution shown in Figure(2.4) for the case h 0 = 1 and v l = 1. The solution is symmetric in x with h(-x, t)=h(x, t) and v(-x, t)=-v(x, t) at all times. A shock wave moves in each direction, bringing the fluid to rest, since the middle state must have v m =0b y symmetry. The characteristic structure of this solution is shown in Figure(2.5). Note again that 1-characteristics impinge on the 1-shock while crossing the 2-shock, whereas 2-characteristics impinge on the 2-shock. Note that if we look at only half of the domain, say x<0, then we obtain the solution to the problem of shallow water flowing into a wall located at x =0w i t h velocity v l . A shock wave moves out from the wall, behind which the fluid is at rest. The numerical tests performed in Chapter (3) refer to the physical situations described in the present chapter.

Figure 2 . 4 :Figure 2 . 5 :

 2425 Figure 2.4: Structure of the similarity solution of the two-shock Riemann problem for the shallow water equations with v l = -v r .Thedepthh,velocity v, and vertically integrated pressure are displayed as functions of x t .T h e structure in the x -t plane is also shown with particle paths indicated for a set of particles with the spacing between particles inversely proportional to the depth. (From Ref.[START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] )

Figure 2 . 6 :

 26 Figure 2.6: (a) Hugoniot locus of points q =( h, hv) in shallow water state space that can be connected to a given state q l by a 1-shock satisfying the Rankine-Hugoniot conditions. Only some of these states (on the solid portion of the curves) satisfy the entropy condition; see Section (2.3). (b)

Figure ( 2 . 3 ) 8 )

 238 Figure(2.3) shows the characteristic structure for the physically correct solution to the dam-break Riemann problem from section (2.1). Figure(2.8)

Figure 2 . 7 :

 27 Figure 2.7: All-shock solutions to the shallow water Riemann problem can be constructed by finding the intersection of the appropriate Hugoniot loci. (a) For Riemann problem with h l = h r = 1, v l =0.5 and v r = -0.5. (b) An entropy-violating Riemann solution for dam-break. (From Ref. [5] )

Figure 2 . 8 :

 28 Figure 2.8: Entropy-violating solution of the dam-break Riemann problem for the shallow water equations, shown in the x-t plane. The dark lines show the shocks. The lighter lines show 1-characteristics and 2-characteristics.

Figure 2 . 9 :

 29 Figure 2.9: (a)Integral curves of the eigenvector r 1 for the shallow water equations. The eigenvector r 1 (q) evaluated at any point on a curve is tangent to the curve at that point. (b)Integrals curves for r 2 (From Ref. [5] )

Figure 2 .

 2 Figure 2.10: (a) Construction for a all-rarefaction Riemann solution for a problem form section (2.1). (b) The Physically incorrect all-rarefaction Riemann solution for a dam-break problem. (From Ref. [5] )

( 2 .

 2 10a) for the Riemann data v l = -0.5,v r =0.5, and h l = h r =1.(2.38)

For

  the symmetric data (2.38) used in Figure (2.10a), the expression (2.40) gives h m =( 4 √ g -1) 2 /16g = 9 16 ,s i n c eu s i n gg = 1. Then either equation from (2.39) gives v m =0 The integral curves in Figure (2.10a) are shown partly as dashed lines.

  Figure (2.2). In Figure (2.7) we saw how to construct a weak solution to

Figure 2 .

 2 Figure 2.11: (a) The Hugoniot loci from Figure (2.7) together with the integral curves. (b) Close-up of the region where the curves intersect. S 1 : Entropy-violating 1-shocks; R 1 : 1-rarefactions; S 2 : 2-shocks: R 2 : unphysical 2-rarefactions. (From Ref. [5] )

Figure ( 2 .

 2 Figure (2.10), since the Hugoniot loci are different from the integral curves. This is illustrated in Figure (2.11), where the curves from Figures (2.7) and

7 ) and ( 2 . 10 )

 7210 are plotted together. Figure (2.10b) shows a close-up near the points of intersection of these curves. The correct solution to the dam-break Riemann problem has the intermediate state at the point where the two solid lines cross.

1 1 ) 2 )

 112 Kinetic approach for the Saint-Venant SystemBy analogy with the Euler equations of gas dynamics, we link the Saint-Venant system to a microscopic description of the fluid and we explain how a kinetic model can be used (in general) to approximate solutions to hyperbolic conservation laws.In the continuum description (1.80) of the fluid, the conserved quantities depth/mass, momentum, and energy are h, hv and E,r e s p e c t i v e l y .T h i s continuum description can be refined by taking the particle structure of the fluid into account. As for the case of rarefied gases, this can be done within the theory of Boltzmann equation. The basic quantity is a distribution function f (t, x, v) which describes the density of particles with velocity v at position x and time t. The fluid particles move freely in space unless they undergo collisions. The corresponding evolution of f is given by a Bolzmann-The left hand side of (3.1) describes free flow of particles whereas collisions are described by the operator Q, which satisfies, for a.e. (t, x),Miko laj Szydlarski 62A connection between the two descriptions (1.79)-(1.80) and (3.1)i s (formally) obtained in a limit where particle collisions are dominant. In this asymptotic case, the state variables f and (h, v, E) as well as the evolutions (1.79)-(1.80) and (3.1) are equivalent. More precisely, the particle distribution function f corresponding to the macroscopic conserved quantities satisfies the following relations,

Figure 3 . 1 :

 31 Figure 3.1: Illustration of a finite volume method for updating the cell average U n i by fluxes at the cell edges. Shown in x -t space. (From Ref. [5] )

  .14) If we can approximate this average flux based on values U n ,t h e nw ew i l l have a fully discrete method. See Figure (3.1) of a schematic for this process.

) 1 2 and F n 1 -1 2 ,

 112 In fact, the sum of the flux differences cancels out except for the fluxes at the extreme edges. Over the full domain we have exact conservation except for fluxes at the boundaries, denoted by F n N + w eh a v et od e fi n e according to the model under study.

Figure ( 3 . 1 2 1 2Figure 3 . 2 :

 31132 Figure 3.2: Characteristic for the advection equation, showing the information that flows into cell C i during a single time step. (a) For a small enough time step, the flux x i-1/2 , depends only on the values in neighboring cellsonly on U n i-1 in this case where a>0. (b) For a larger time step, the flux should depend on values father away. (From Ref. [5] )

For a hyperbolic system

  of equations there are generally a set of m wave speeds λ 1 ,...,λ m where m ∈ N + (see Section 1.2.2 for example of water Miko laj Szydlarski 69 CHIARA SIMEONI 3. Numerical Simulation of the Shallow Water equations small waves). In this case we define the Courant number by

  dx. Starting from the microscopic equation (3.1),

  ξ) holds for the discrete contribution of the source term hZ (x) in the system, for negative velocities;

Figure 3 . 3 : 2 ) 2 of 2 , 2 =(F h ) + i+ 1 2 ,

 3322222 Figure 3.3: Illustration of the typical sytuation occuring in a celll C i of the mesh, centered at the point x i ∈ R.( From Ref. [7] )

CHIARA SIMEONI 3 .

 3 Numerical Simulation of the Shallow Water equations3.4.1 Dam BreakWe begin with a Dam-Break test in a rectangular channel with flat bottom (Z = 0). The initial conditions arev(0,x)=0 h(0,x)= h l for x ≤ 0 h r for x>0,where h l >h r in order to be consistent with the physical phenomenon of a Dam-Break from the left to the right. Note that this case corresponds to a Riemann problem described in Section (2.1) Figures (3.

  4) and (3.5)presen ts the results observed at time t =0s and t = 120s, for a Dam-Break on a wet bed (h l =4.0m, h r =1.25m).

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Numerical Test: Dam break on a wet bed -initial water level.

6 )

 6 presents the initial water level with small rectangular projection, symmetric around the point x = 0 and the Figure (3.7) presents the results observed at time t = 120s.

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Numerical Test: Small perturbation -initial water level.

Figure 3 . 8 :

 38 Figure 3.8: Numerical Test: Perturbation over bump -initial water level.

Figure 3 . 9 :

 39 Figure 3.9: Numerical Test: Perturbation over bump -final water level

Figure ( 3 .Figure 3 . 10 :Figure 3 .

 33103 Figure (3.12) shows the program in action with the windows open for each program module. In the case when the PLOT program is not available on the computer, the base program skips plot and animation process and it just writes the output data in series of *.dat files (more informations about using the program can be found in the PDF version of the user guide inside the CD jointed to the present manuscript). The scheme of interface that we

CHIARA SIMEONI