CFD Online Discussion Forums - Shallow Water Equations: Roe Approximate Solver First-Order with Wetting-Drying

Hi,

Firstly, I apologise for the length of this post, but it requires quite a bit of detail and explanation! Thanks in advance for any help you might provide...

I'm currently writing a structured 2D solver using a Godunov-type Finite Volume Method (first-order, i.e. no slope-reconstruction) for the shallow water equations with variable bed slope and friction, in C++.

The code uses the Roe approximate Riemann solver with incorporated bed slope flux, with friction solved in a point-wise manner.

To date the solver works with all test cases I have tried where the water height h (h being the height above the bed) exceeds the height of the bed.

I'm having some trouble with creating effective conditions for what is known as the "wetting-drying condition", which is essentially a moving boundary problem where the equations are solved at adjacent cells, one being "wet" and the other "dry". Dry is defined here as the water depth being less than a pre-specified tolerance value.

I've spent a good bit of time looking at some papers ([1], [2]) but I'm finding it tricky to tease out the algorithm as to how a wetting-drying front is actually implemented. There is mention of a situation where the total surface height H = h + z, where z is the bed height, must have H_L >= H_R, otherwise a modification to the bed slope, z, is needed.

My current algorithm for the Roe fluxes looks a bit like this:

// Only showing x-direction, y-direction is similar
Loop over i
----Loop over j
--------set "left cell" and "right cell"
--------if both cells are "dry" (i.e. h_L < tolerance && h_R < tolerance)
------------set all fluxes to zero
--------else // at least one is dry
------------if it's left one and H_L < H_R
----------------// use correction as in [1]
----------------z_R = z_L - (h_R - h_L);
----------------// set u, v velocities to zero
------------if it's right one and H_R < H_L
----------------// use correction as in [1]
----------------z_L = z_R - (h_L - h_R);
----------------// set u, v velocities to zero
--------// Do flux calculations (i.e. work out alpha, eigenvalues and eigenvectors)
--------// for fluxes and bed slop terms

Firstly, does anybody have any experience implementing an algorithm like this?

I've been trying for the last couple of days with various modifications to the above algorithm and nothing seems to be working. I've -nearly- managed to replicate the periodic tidal case in [1], but this is only with setting a rather high tolerance (>0.05m).

So my questions are:

1) Is my algorithm even correct? Do I need to set some flux terms to zero? Do I need to set more velocity terms to zero?

2) What is a good value to use as a tolerance? Most papers seem to ignore this, but it's an essential model parameter, so I'm left trying arbitrary/reasonable values!

Please let me know if you need more detail in answering the question. I can provide the source code, if necessary.

Thanks in advance for any help provided!

Cheers,

Mike.

---------

[1] Brufau et al. - A numerical model for the flooding and drying of irregular domains - International Journal For Numerical Methods in Fluids, May 2002.

[2] Brufau et al. - Zero mass error using unsteady wetting-drying conditions in shallow flows over dry irregular topography, International Journal For Numerical Methods in Fluids, August 2004.