Pascal Poullet
email: pascal.poullet@univ-antilles.fr

Priscilla Ramsamy
email: priscilla.ramsamy@univ-antilles.fr

Mario Ricchiuto
email: mario.ricchiuto@inria.fr

Residual based method for sediment transport

come

Introduction

This work is incorporated within the framework of the study of a sediment transport modelling. One aim of this contribution, is to provide first, an useful simulation tool, in the context of a 1D space-time problem. Considering the geophysical aspect, the sediment transport can be divided into several categories, in this paper we will be interested in bedload transport. The governing equations consist of a system of three equations, modelling the interaction between fluid and sediment in a river. The hydrodynamical component is described by the shallow-water equations (SWE) and the morphodynamical component, is given by a solid transport discharge due to the Exner law with the Grass formula. In this way, the model can be depicted by a non-conservative and non-linear hyperbolic system, and our main objective is to seek numerical solutions in accordance with these specific aspects. In particular, a good scheme for that model, must comply with the well-balanced property and the path-conservative character. Moreover, as the fluid interacts very weakly with the sediment and characteristic velocities being such a different magnitude, long time simulation and high order accuracy (at least second order) are needed.

To treat hyperbolic problem, finite volumes are the most popular methods as the Godunov scheme. For conservation laws, first attempts to propose approximate solver for hyperbolic systems in non-conservative form, were due to Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors and difference schemes[END_REF]. After that, several approaches have been introduced like approximate finite volume Roe with characteristic flux scheme [START_REF] Ghidaglia | Une méthode de volumes finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF], schemes based on exact or incomplete Riemann solvers [START_REF] Benkhaldoun | Solution on the sediment transport equations using a finite volume method based on sign matrix[END_REF], WENO schemes [START_REF] Xing | High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF], generalized Roe methods with or without WENO reconstruction [START_REF] Castro | High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. applications to shallow water systems[END_REF][START_REF] Castro Diaz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF][START_REF] Castro | On some fast well-balanced first order solvers for nonconservative systems[END_REF], or kinetic schemes [START_REF] Perthame | A kinetic scheme for the Saint-Venant system with a source term[END_REF]. In parallel to finite volume, another family of methods called Residual Distribution (RD) methods that emerged from Roe's works in the 80s, is used in this paper [START_REF] Roe | Approximate Riemann solvers, parameter vectors and difference schemes[END_REF][START_REF] Deconinck | Residual distribution schemes: foundations and analysis[END_REF]. Combining advantages from finite volume and finite element methods, their construction allows them nowadays to be monotone, conservative, well-balanced and easily high order accurate [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Warzynski | Runge-Kutta residual distribution schemes[END_REF][START_REF] Arpaia | An ALE formulation for explicit Runge-Kutta residual distribution[END_REF][START_REF] Abgrall | High-Order Methods for CFD[END_REF][START_REF] Mazaheri | Improved second-order hyperbolic residual-distribution scheme and its extension to third-orer onn arbitrary triangular grids[END_REF][START_REF] Lin | High order residual distribution for steady state problems for hyperbolic conservation laws[END_REF]. In this contribution, a RD scheme, viewed as a recast of the approximate Riemann solver, called FV-Roe approximation scheme is proposed with simulations. Often used to solve the shallow-water problem, the residual based method is adapted here to solve the coupling problem with sediment transport. And for that the use of a TVD Runge Kutta procedure, but also upwinding and a flux limiter procedures, have been added, to compute weak entropic solution, considering Lax entropy. To introduce our scheme, the present paper is organized according to the following outline: In section 2, the governing equations are introduced. In section 3, our scheme is proposed. And finally, numerical tests are given, in order to show how accurate the scheme is, but also its well-balanced preserving aspect, for the lake at rest.

A sediment transport in a shallow water

In the context of the study of sediment transport in shallow water, several morphodynamical models can be found in the litterature depending on the way of considering the displacement mode. In the case of bedload transport, among several models [START_REF] Meyer-Peter | Formulas for Bed-Load transport[END_REF][START_REF] Van Rijn | Sediment transport (i): bed load transport[END_REF][START_REF] Kalinske | Movement of sediment as bed load in rivers[END_REF][START_REF] Einstein | Formulas for the transportation of bed load[END_REF], the discharge is written by the Grass formula [START_REF] Grass | Sediments transport by waves and currents[END_REF] for simplicity. It involves that the critical shear stress is neglected, then the sediment is viewed as starting its own movement as soon as the fluid starts to move. About the hydrodynamical part, shallow-water equations are considered. The result is an hyperbolic system of three equations. The governing system for the 1D space-time problem, is as follows

               ∂ h ∂t + ∂q ∂ x = 0 ∂q ∂t + ∂ ∂ x q 2 h + 1 2 gh 2 = gh ∂S ∂ x -ghS f ∂S ∂t -ζ ∂q b ∂ x = 0, (1)
where x denote the horizontal variable at the axis of the channel and t the time variable (see [START_REF] Castro Diaz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF]). By h(x, t) we denote the height of the water column, q(x, t) is the discharge, g is the gravity constant, S f models the friction term and ζ a parameter linked to the sediment porosity (ζ = 1/(1ρ 0 ) with ρ 0 the porosity). The third equation of the model describes the sediment transport by the expression of the sediment volume equation, q b being the solid transport discharge obtained by Grass formula (here, q b = A g (q/h) 3 with A g related to the interaction between the fluid and the sediment). The variable S(x, t) is the distance from a given reference level to the bottom layer. A schematic description is provided (see Fig. 1), η denoting the extra height of the water column.

Neglecting S f in this study, a classical approach is to treat the system (1) as a hyperbolic system with a non-conservative term B:

∂W ∂t + A(W) ∂W ∂ x = 0, (2)
with W = W(x, t) and (x, t) ∈ I R × I R + . In fact, the vector of unknowns is

W = (h, q, S) T ,
F is the flux function and, A(W) equals to the difference between the Jacobian matrix of F and the non-conservative part:

A(W) = ∂F ∂W (W) -B(W).
More precisely,

F(W) = q q 2 h + 1 2 gh 2 -ζ q b , B = 0 0 0 0 0 gh 0 0 0 , A(W) = 0 1 0 - q 2 h 2 + gh 2 q h 0 -ζ ∂q b ∂ h -ζ ∂q b ∂q 0 , and A(W) = 0 1 0 - q 2 h 2 + gh 2 q h -gh -ζ ∂q b ∂ h -ζ ∂q b ∂q 0 .

A residual based predictor-corrector upwind discretization for 1D space-time sediment transport

Following the introduction, the scheme built in this work, can be viewed as an high order recast of the approximate FV-Roe as introduced by Ghidaglia et al. in [START_REF] Ghidaglia | Une méthode de volumes finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF]. Instead of using a high order reconstruction, our approach exploits the residual based approach discussed in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF].

To present the scheme, one neglects the friction term, and one considers the hyperbolic system with a non-conservative term. Also, to simplify writing in this section, in integrals the term dx will be omitted (all are space integrals). One will consider the intervals (computing cells) defined by

I i = [x i-1 2 , x i+ 1 2 ], i ∈ Z Z, but also the intervals I i+ 1 2 = [x i , x i+1 ]. Let the step ∆x i = x i+ 1 2 -x i-1 2
and that

x i-1 2 = i-1 k=1
∆x k is the intercell located at the middle of I i-1 ∪ I i . ∆t is the time step and t n = n∆t. As usual, we denote by W n i the approximate mean value of W in node x i and at time t n . The RD procedure consists of making the computation of the residual, called global, on a cell I i and then distributing fractions of this quantity to each of its verteces. Under these assumptions, one gets for a linear approximation of W(x, t n ),

W n i := 1 ∆x i ∫ I i W(x, t n ) ≈ W(x i , t n ).
Then, the residual-based predictor corrector method developed in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] can be written as follows

W n+1 i = W n i - ∆t ∆x 1 2 Θ i (W n ) + 1 2 Θ i (W * ) + Ψ i-1 2 + Ψ i+ 1 2 , (3)
with Θ i (W) corresponds to the approximate FV-Roe scheme discrete evolution operator (without the additional term)

Θ i (W) = P + (A i-1 2 )φ i-1 2 + P -(A i+ 1 2 )φ i+ 1 2 .
Concerning the non-conservative terms, one proceeds by linearization along the path joining W i and W i±1 , to compute an approximate value of

W i± 1 2 = 1 2 (W i + W i±1 ). We then denote by A i± 1 2 = A(W i± 1 2 ), B i± 1 2 = B(W i± 1 2
). Also, the fluctuations are:

φ i-1 2 = F i -F i-1 -B i-1 2 (W i -W i-1 ) and φ i+ 1 2 = F i+1 -F i -B i+ 1 2 (W i+1 -W i ),
and the projectors,

P ± (A i∓ 1 2 ) = 1 2 (I ± sign A i∓ 1 2 )),
with the sign of a matrix computed by eigen-decomposition

sign(A i∓ 1 2 ) = K i∓ 1 2 sign(L i∓ 1 2 )K -1 i∓ 1 2 .
The matrix K gathering the eigenvectors of the Roe matrix (along the path), sign(L) is the diagonal matrix whose coefficients are the sign of the eigenvalues (see [START_REF] Parés | On the well-balance property of roe's method for nonconservative hyperbolic systems. applications to shallow-water systems[END_REF] for example). We denote by W * , a predicted value of the solution that has been obtained, from the upwind scheme,

W * i =W n i - ∆t ∆x i Θ i (W n ).
Actually, if we introduce a parallel approach using Galerkin method by seeking a piecewise linear solution, in a domain

Ω =]0, L[, W h (x, t) = i ϕ i (x)W i (x i , t), and
F h = F(W h ),
with ϕ i representing the standard Lagrange basis fonctions associated to the node x i . One replaces the unknown solution by its approximate solution by finite element from the variational form

∫ Ω ϕ i ∂ t W h + ∫ Ω ϕ i ∂ x F h - ∫ Ω ϕ i B(W h )∂ x W h = 0. (4)
And one can notice that by linear approximation,

∫ Ω ϕ i ∂ x F h = ∫ I i-1 2 ϕ i ∂ x F h + ∫ I i+ 1 2 ϕ i ∂ x F h ≈ F i -F i-1 2 + F i+1 -F i 2 . ( 5
)
Hence, the resulting non-stabilized method reads

∫ Ω ϕ i ∂ t W h + 1 2 φ i+ 1 2 + 1 2 φ i-1 2 = 0, (6)
with the approximation of the first term ∆x i dW i dt , one recovers the mass lumping process.

But as it is known that Galerkin method suffers from lack of stability, one adds a residual based stabilization in the spirit of the streamline upwind method or Streamline Upwind Petrov Galerkin (SUPG) [START_REF] Hughes | Finite element formulation for computational fluid dynamics : I symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics[END_REF][START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF]. For a node x i , the stabilization operator S i reads

S i = ∫ Ω A(W h )∂ x ϕ i T r with r = ∂ t W h + ∂ x F h -B(W h )∂ x W h ,
the matrix T being a scaling factor guaranteeing the uniform boundedness of the stabilization w.r.t. the residual. As before, explicit computable expressions are obtained when introducing the linear finite element approximation and introducing appropriate mean value linearizations of the matrices that appear. Recalling that for a linear approximation

∂ x ϕ i I i± 1 2 = ∓1/∆x i± 1 2
, the stabilization term can be evaluated as

S i = A i+ 1 2 T i+ 1 2 ∫ I i+ 1 2 ∂ t W h + φ i+ 1 2 -A i-1 2 T i-1 2 ∫ I i-1 2 ∂ t W h + φ i-1 2 (7)
In one dimension, a typical definition of the scaling matrix T , also used here, being the following

T i± 1 2 = 1 2 ∆x i± 1 2 sign(A i± 1 2 )A -1 i± 1 2 = ∆x i± 1 2 2 |A i± 1 2 | -1 , (8)
the complete semi-discrete (in space) equations read ∫

I i-1 2 ϕ i + sign(A i-1 2 ) ∂ t W h + ∫ I i+ 1 2 ϕ i -sign(A i+ 1 2 ) ∂ t W h = -P + (A i-1 2 )φ i-1 2 + P -(A i+ 1 2 )φ i+ 1 2 . (9)
As it has been explained in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF], a simplified stabilization step can be added in the Galerkin process with the Runge-Kutta scheme without any loss of accuracy (here, an explicit RK2 of second order accuracy is used). Thus, the residual

r n+1 := W n+1 h -W n h ∆t + 1 2 (∂ x F h -B h ∂ x W h ) * + 1 2 (∂ x F h -B h ∂ x W h ) n (10)
can be replaced with a simplified residual

r * := W * h -W n h ∆t + 1 2 (∂ x F h -B h ∂ x W h ) * + 1 2 (∂ x F h -B h ∂ x W h ) n , (11)
to design the scheme by computing

∫ Ω ϕ i r n+1 + ∫ I i-1 2 ∪I i+ 1 2 A(W h )∂ x ϕ i T r * = 0. ( 12
)
Using an explicit scheme by performing mass lumping in the predictor step [START_REF] Benkhaldoun | A flux-limiter method for dam-break flows over erodible sediment beds[END_REF], and with midpoint rule to evaluate the integrals, one obtains after few calculations and recast equation ( 3), that one recalls

W n+1 i = W n i - ∆t ∆x 1 2 Θ i (W n ) + 1 2 Θ i (W * ) + Ψ i-1 2 + Ψ i+ 1 2 ,
with the aim to define Ψ i± 1 2 as:

ψ i± 1 2 = ∆t ∆x ∫ I i± 1 2 ϕ i - 1 2 (I ∓ sign(A i± 1 2 )) W * h -W n h ∆t ≈ 1 2 W * i -W n i -(I ∓ sign(A i± 1 2 ))(W * i± 1 2 -W n i± 1 2
) .

However, comparing to the FV-Roe scheme, the additional terms Ψ i± 1 2 that derive from residual stabilization are of second order, that cannot always match with non regular solution. To avoid their effects accross discontinuous features, cell based limiters have been introduced and finally,

Ψ i± 1 2 = δ i± 1 2 ψ i± 1 2 with δ i± 1 2
computed by means of a standard finite volume limiters using values of the MUSCL MinMod flux limiter function [START_REF] Bilanceri | Linearized implicit time advancing and defect correction applied to sediment transport simulations[END_REF].

Numerical results

We present the numerical results of several reference problems. In this section, we aim at validating our numerical scheme and highlighting its characteristic properties against classical tests. The first test consists in proving an approximate well-balanced property of our scheme. We then underline the ability of our scheme to simulate a parabolic sediment transport until a discontinuous solution is obtained. We also prove its high order accuracy by means of an order test problem that has been discussed in [START_REF] Castro Diaz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF]. Finally, we prove that our scheme is capable of faithfully reproducing a dam-break problem over a wet bottom topography [START_REF] Audusse | A simple three-wave approximate riemann solver for the saint-venant-exner equations[END_REF].

Test of well-balanced property

To check the property, the following numerical test is used [START_REF] Qian | Well-balanced central weno schemes for the sediment transport model in shallow water[END_REF]. It deals with the ability of the scheme to reproduce the behaviour of the steady state. Thus, if the numerical scheme is well-balanced, a small difference should be observed between the initial solution and the solution obtained at the final instant. For this, the interval [0, 10] is assumed as the physical domain and the simulation is performed up to T = 0.5s with 100 and 200 cells. A discontinuity in the bed is assumed as the initial condition, so the thickness of the sediment layer is considered,

z b (x, 0) = 4 if 4 ≤ x ≤ 8 0 elsewhere ( 13
)
and q(x, 0) = 0, h(x, 0) + z b (x, 0) = 10. The results in Table 1, show that the scheme preserves the approximate wellbalanced property [START_REF] Munoz-Ruiz | On the convergence and well-balanced property of pathconservative numerical schemes for systems of balance laws[END_REF]. The differences between the initial solution and the solution at final time are very small. More precisely, the ratio is of around 1.5 between the doubled gridpoints and the coarser one. The accuracy seems to be of order 1.5.

A 1D space-time dune test case

To verify the shock capturing property, the classical transport of parabolic sediment layer has been taken. For this case proposed in [START_REF] Hudson | Numerical techniques for morphodynamic modelling[END_REF][START_REF] Bilanceri | Linearized implicit time advancing and defect correction applied to sediment transport simulations[END_REF][START_REF] Castro Diaz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF], the interval [0, 1000] is assumed as the physical domain and a strong interaction between the fluid and the sediment is taken (A g = 1). The initial conditions are given as follows: the bathymetry is of 10 m, q(x, t) = 10 for the discharge of the fluid, h(x, t) = 10z b (x, t) for the water column height with the sediment layer thickness,

z b (x, 0) = 0.1 + sin 2 π(x-300) 200 if 300 ≤ x ≤ 500 0.1 elsewhere . (14)
Numerical solutions are generated up to T = 700s. As shown in Fig. 3, the solution of our scheme (solid line) is compared to the solution of the FV-Roe approximation Fig. 2 The dune at initial time. (dotted line), and we note that the RD scheme is less diffusive than the first and the shock is more visible (thanks to the limiters). Therefore, the results show that our numerical scheme seems to satisfy the shock capturing property for this test. It can be noted that it is also stable, more precisely the RD scheme does not generate oscillations even though the interaction between the fluid and the sediment is strong (somewhat artificial with the initial conditions chosen [START_REF] Hudson | Numerical techniques for morphodynamic modelling[END_REF]). The comparison between the evolution of the dune between the RD scheme and the SRNH scheme (a secondorder scheme extracted from [START_REF] Benkhaldoun | A flux-limiter method for dam-break flows over erodible sediment beds[END_REF], see fig. 4 and5) seems to confirm that our scheme is of the second order and is less diffusive than the other.

A test of order

To check the accuracy of the numerical scheme, let us introduce the following one-dimensional problem for which initial conditions are: This test problem has been considered in a previous work by Castro Diaz et. al. [START_REF] Castro Diaz | Sediment transport models in shallow water equations and numerical approach by high order finite volume methods[END_REF], and as the exact solution is unknown, we use a reference solution obtained by a fine mesh of 5120 volumes (as it has been done in their work) with a medium interaction (A g = 0.3). Numerically, the results for h, q and the sediment layer thickness (z b ) and different error norms, show that the second order of accuracy seems to be obtained (see Tables 2, 3 and4).

q(0, x) = 0, h(0, x) = 2 -0.1 exp(-x 2 ), z b (0, x) = 0.1 -0.01 exp(-x 2 ). (15)

A dam break test over a wet bottom topography

In this classical test case ( [START_REF] Audusse | A simple three-wave approximate riemann solver for the saint-venant-exner equations[END_REF]) a dam break is considered over a flat wet bottom, in a channel of 10m long. A low interaction between the fluid and the sediment is taken (A g = 0.005), and the initial conditions are,

h(x, 0) = 2 if x ≤ 5m 0.125 if x > 5m , (16)
q(x, 0) = 0m/s and the bottom topography z b = 0m. The numerical test is performed until T = 1s. The results confirm that our scheme keep the stability as attempted for the coupled approaches (see Fig. 6,7), in comparison the approach by splitting [START_REF] Audusse | Sediment transport modelling: relaxation schemes for saint-venant-exner and three layer models[END_REF]. The accuracy of our RD scheme is of course better than those obtained by the FV-Roe scheme (8). More precisely, as expected and considering the bottom firstly (Fig. 6) the solution computed with the RD scheme forms two plates without any oscillation (in a critical region includes the interval [START_REF] Audusse | A simple three-wave approximate riemann solver for the saint-venant-exner equations[END_REF][START_REF] Benkhaldoun | A flux-limiter method for dam-break flows over erodible sediment beds[END_REF] with -0.1 < z b < -0.05, and in a region where x > 6.5 with z b > 0.5). Then, for the free surface (Fig. 7), the solution of the RD scheme (dashed line) is decreasing along the time, and a plate is reached without oscillation in the critical region [START_REF] Audusse | A simple three-wave approximate riemann solver for the saint-venant-exner equations[END_REF][START_REF] Benkhaldoun | A flux-limiter method for dam-break flows over erodible sediment beds[END_REF]. For both unknowns, a comparison with the FV-Roe approximation (dotted line) is proposed, underlying that the results are quite similar. However the RD scheme is more accurate for the bottom (see Fig. 8 for which the L 2 errors are produced from a referent solution computed with 2000 elements grid).

Concluding remarks

This contribution proposed a new predictor-corrector scheme, based on the residual, to simulate a sediment transport problem. Numerical tests have highlighted its high order accuracy, its approximate well-balancedness property and its stability for some test problems. In particular, the solutions obtained for the problem of dam-break with wet topography are sharp. Work is in progress to take into account dry bottoms, for example. The extension of this model to take into account a coastal configuration (2D physical domain), by parallel programming, will also be done in a future contribution.

Fig. 1 A

 1 Fig. 1 A sediment layer in a shallow water

Fig. 3 z

 3 Fig. 3 z b at 700 s for RD scheme and FV-Roe.

Fig. 4 zFig. 5 z

 45 Fig. 4 z b at 700 s for RD scheme and SRNH.

Fig. 6

 6 Fig. 6 Dam break over a wet bottom: h at t = 1s for RD scheme with 1000 cells.

Fig. 7

 7 Fig. 7 Dam break over a wet bottom: z b at t = 1s for RD scheme with 1000 cells.

Fig. 8

 8 Fig. 8 Dam break over a wet bottom: error analysis for RD scheme.

Table 1

 1 Accuracy of the scheme for the well-balanced test property

	precision	L 2 -error h	ratio	L 2 -error q	ratio
	100	8.6052 × 10 -16	-	7.1712 × 10 -15	-
	200	5.2050 × 10 -16	1.65	5.8743 × 10 -15	1.22

Table 2

 2 Accuracy of the RD-scheme with distributed residual scheme for h

	# gridpts	L 1 -error	Order L 2 -error	Order L ∞ -error	Order
	20	5.628 × 10 -3	-	8.600 × 10 -3	-	2.088 × 10 -2	-
	40	2.421 × 10 -3	1.22	4.159 × 10 -3	1.05	1.086 × 10 -2	0.94
	80	7.918 × 10 -4	1.61	1.408 × 10 -3	1.56	4.036 × 10 -3	1.43
	160	2.072 × 10 -4	1.93	3.815 × 10 -4	1.88	1.211 × 10 -3	1.74
	320	5.257 × 10 -5	1.98	9.693 × 10 -5	1.98	3.159 × 10 -4	1.94
	640	1.314 × 10 -5	2.00	2.423 × 10 -5	2.00	7.927 × 10 -5	1.99

Table 3

 3 Accuracy of the RD-scheme with distributed residual scheme for q, the discharge

	# gridpts	L 1 -error	Order L 2 -error	Order L ∞ -error	Order
	20	2.350 × 10 -2	-	3.771 × 10 -2	-	9.248 × 10 -2	-
	40	1.042 × 10 -2	1.17	1.835 × 10 -2	1.04	4.783 × 10 -2	0.95
	80	3.435 × 10 -3	1.60	6.205 × 10 -3	1.56	1.783 × 10 -2	1.42
	160	8.995 × 10 -4	1.93	1.681 × 10 -3	1.88	5.349 × 10 -3	1.73
	320	2.281 × 10 -4	1.98	4.269 × 10 -4	1.98	1.397 × 10 -3	1.94
	640	5.697 × 10 -5	2.00	1.067 × 10 -4	2.00	3.510 × 10 -4	1.99

Table 4

 4 Accuracy of the RD-scheme with distributed residual scheme for z b , the height of the

	sediment layer
	# gridpts	L 1 -error	Order L 2 -error	Order L ∞ -error	Order
	20	2.969 × 10 -5	-	5.443 × 10 -5	-	1.423 × 10 -4	-
	40	1.265 × 10 -5	1.23	2.739 × 10 -5	0.99	8.168 × 10 -5	0.80
	80	4.469 × 10 -6	1.50	9.014 × 10 -6	1.60	2.861 × 10 -5	1.51
	160	1.200 × 10 -6	1.90	2.536 × 10 -6	1.83	8.578 × 10 -6	1.74
	320	3.019 × 10 -7	1.99	6.409 × 10 -7	1.98	2.217 × 10 -6	1.95
	640	7.516 × 10 -8	2.01	1.596 × 10 -7	2.00	5.571 × 10 -7	1.99

Acknowledgements

Computational tests have been performed using Wahoo, the server of the Centre Commun de Calcul Intensif of the Université des Antilles.


L Batteux
email: lea.batteux@univ-antilles.fr

J Laminie

J C Latché

P Poullet
email: pascal.poullet@univ-antilles.fr

3-dimensional particulate flow modelling using a viscous penalty combined with a stable projection scheme

Keywords: 65M08, 76D05, 76T20

Introduction

This work is focused on the modelling of fluid-solid systems in a 3-dimensional domain. To reproduce faithfully the fluid-solid interactions is a problem of large interest due to numerous processes in industrial applications. There are several methods to attempt to model such a problem, but in this study, one considers rigid solid inclusions in an incompressible viscous fluid flow and one enforces a strong coupling between both phases. The motion of the solid domain may then be described using Newton laws for rigid bodies. As we are concerned with the efficiency and computational costs, we resort to an Eulerian formulation for the fluid flow and extend the fluid problem inside the solid domain in the manner of fictitious domain methods [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Hoffmann | On a motion of a solid body in a viscous fluid. twodimensional case[END_REF][START_REF] Carbou | Boundary layer for a penalization method for viscous incompressible flow[END_REF]. In our case we enforce some kind of Brinkmann law inside the solid domain by adopting the H1-penalty method [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Hoffmann | On a motion of a solid body in a viscous fluid. twodimensional case[END_REF]. In practice, it will come down to penalizing the tensor of deformation term where the particles are in the nonhomogeneous Navier-Stokes equations. Let us denote by Ω the domain containing the particulate flow. The continuous problem is given by:

                 ∂ ρ ∂t + ∇ • (ρv) = 0 in R + × Ω ∂ (ρv) ∂t + ∇ • (ρv ⊗ v) -2∇ • (µ(ρ)D(v)) + ∇p = f in R + × Ω ∇ • v = 0 in R + × Ω (1)
with the unknowns being the density ρ, velocity v and pressure p. The source term is denoted by f and contains the forces applied to the particles among other exterior forces. The advection of the rigid particles is carried out by the mass balance equation rather than Newton laws. The viscosity µ is continuously dependent on ρ, and taking µ → ∞ inside the solid domain allows for the penalization of the tensor of deformation rate D(v). We aim to tend towards D(v) L 2 (Ω s (t)) = 0 which is equivalent to a rigid motion velocity field inside the particles Ω s (t). The system of equations ( 1) is complemented with initial conditions for (ρ, v).

The flow being incompressible, one often resorts to a projection scheme [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Févrière | On the penalty-projection method for the Navier-Stokes equations with the mac mesh[END_REF] to solve problem [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]. Many variants can be found in the litterature especially when considering multiphase flows. However a lot of problems remain open regarding the stability or convergence of those schemes. Again, for efficiency reasons we will rely on the scheme introduced in [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure poisson equation[END_REF]. In this article, the authors circumvent the relatively expensive computational cost of the elliptic pressure problem of standard projection problem for incompressible variable density flows, by switching to an approximate and more efficient formulation of the latter. In this paper, we adapt the scheme presented in [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure poisson equation[END_REF] for finite elements to the MAC discretization, while keeping the same stability properties. Additionally, the discontinuity and the jumps of viscosity along the fluid/solid interface require an accurate tracking of the surface of the particles.

To this end, we replace the discrete mass balance equation with an anti-diffusive advection scheme introduced in [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] that is similar to the Ultra-Bee method. The scheme is adapted to the dimension d = 3 by considering an alternate direction variant.

In the following section we introduce the notations, meshes as well as the full discrete scheme. In a final section we carry out and comment on the simulation of the fall of a rigid sphere.

Numerical Method

Let 0 = t 0 < t 1 < • • • < t N = T be a uniform partition of the time interval [0, T ].
We note δt = T /N the time step so that t n = nδt, for n ∈ 0, N . The incremental projection scheme from [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure poisson equation[END_REF] reads,

ρ n+1 -ρ n δt + ∇ • (ρ n+1 v n ) = 0 in Ω (2a) ρ n+1 v n+1 -ρ n v n δt + ∇ • (ρ n+1 v n+1 ⊗ v n ) -∇ • ((µ n+1 D(v n+1 ))) +∇(2p n -p n-1 ) = f n+1 in Ω v n+1 = 0 on ∂ Ω (2b) -δt χ ∆ (p n+1 -p n ) = -∇ • v n+1 in Ω (2c) v n+1 = v n+1 -δt χ ∇(p n+1 -p n ) in Ω (2d)
given for any time increment t n+1 . We denote f(t n+1 ) = f n+1 and we take χ = min x∈Ω ρ 0 .

We aim to compute the sequence of discrete solution (v n+1 , v n+1 , p n+1 , ρ n+1 ) for n ≥ 0. For each time step, we obtain a tentative velocity v n+1 by solving equation (2b). Using the Helmholtz decomposition of L 2 (Ω ), the substep (2c) acts as the projection of

v n+1 on H = {v ∈ L 2 (Ω ), ∇ • v = 0, v • n | ∂ Ω = 0}
to get a corrected divergence-free velocity v n+1 in (2d). The advection of the density is carried out by (2a). The specificity (and interest) of the present scheme is to replace the actual density by the constant coefficient χ ; indeed, a direct adaptation of the classic incremental projection scheme [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] to nonhomogeneous fluids would result in a variable Poisson problem of the form ∇ • ((ρ n+1 ) -1 ∇Φ n+1 ) = 0 with Neumann boundary conditions. For discontinuous densities with a high ρ max /ρ min ratio, this problem can be expensive to solve due to its ill-conditioned status.

Notations, mesh and discrete projection scheme

The domain Ω is discretized according to a staggered MAC mesh D = (M , E ) so the scheme (2) can benefit from the infsup stability property. Let the primal grid M consists in a conforming structured partition of Ω using rectangular parallelepipeds elements. The parallelepipeds will be defined as primal cells and noted K. Therefore we have ∪ K∈M K = Ω . We may assume that the faces of the primal cells are normal to one vector of the standard basis of R 3 , denoted by (e 1 , . . . , e d ). A face of the primal cell K ∈ M will be noted σ ∈ E (K), E (K) referring to the set of all faces of K.

The staggered grid is completed by defining the dual grid E as the set of all edges of M : {σ ∈ E (K)|K ∈ M }. We note E = E int ∪ E ext , where E int (resp. E ext ) are the edges of E that lie in the interior (resp. on the boundary) of the domain. The set of faces that are orthogonal to the i th unit vector e i of the canonical basis of R d is denoted by E (i) , for i = 1, . . . , d. Correspondingly we introduce

E (i) int = E int ∩ E (i) , E (i) ext = E ext ∩ E (i) so E (i) = E (i) int ∪ E (i)
ext . For σ ∈ E int , we note σ = K|L for (K, L) ∈ M 2 such that ∂ K ∩ ∂ L = σ and we associate the dual cell D σ with D K,σ ∪ D L,σ , where D K,σ (resp. D L,σ ) is the halfpart of K (resp. L) adjacent to σ . If σ ∈ E ext is adjacent to the cell K, then D σ = D K,σ . We can define Ω from the dual mesh: Ω = ∪ σ ∈E i D σ , i = 1, . . . , d. A dual face separating two duals cells D σ and D σ is denoted by ε = σ |σ . In agreement with the staggered MAC scheme, we will define the unknowns (ρ, p) on the primal grid and the ith component of the velocity on E (i) in such a way that we deal with quantities (ρ K , p K ) K∈M and (u σ ) σ ∈E . The grids and notations for d = 2 are illustrated in Fig. 1.

D σ D σ σ = K|L K L σ σ = K|L σ u σ • e 2 × u σ • e 2 × ε = σ |σ D K,σ ρ K , p K × Fig. 1 2-dimensional representation of (M , E )
The discretization of problem (2) begins with the discrete approximation of the mass balance equation (2a). For K ∈ M , we resort in the classic way to the divergence formula for the computation of (2a) integrated over the primal cell K. This yields:

ρ n+1 K -ρ n K δt + 1 |K| ∑ σ ∈E (K) σ ∈E (i) int F n+1 K,σ = 0, (3)
where F n+1 K,σ refers to the mass flux across the primal face σ outward K. In section 2.2, we introduce the two techniques we adopt to compute F n+1 K,σ ; namely the classic upwind scheme and an antidiffusive scheme inspired by [START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF].

Let us focus on the discretization of the prediction step. For σ = K|L ∈ E (i) int , the approximation of (2b) is obtained by integrating the ith prediction equation over the associated dual cell D σ . In order for the scheme to be stable and provide the desired estimates for the unknowns, we must pay a particular attention to the discretization of the convective term in (2b):

1 |D σ | D σ (ρ n+1 u n+1 -ρ n u n )/δt + ∇ • (ρ n+1 u n+1 ⊗ u n ) • e i dx
so that the discretization of the prediction step is compatible with the the discrete mass balance. This yields an approximation of the form:

(ρ n+1 D σ u n+1 σ -ρ n D σ u n σ )/δt + 1 |D σ | ∑ ε∈E (D σ ) F n+1 σ ,ε u n+1 ε,i
where the values for the density on D σ , denoted ρ D σ (resp. the fluxes on the dual faces σ , noted as F n+1 σ ,ε ) are given as functions of the density on the primal cells K, L (resp. the fluxes on the primal faces). This is achieved by averaging the discrete mass balance over K and L to obtain a consistent mass balance equation on D σ . Therefore we define

|D σ |ρ n D σ = |D K,σ |ρ n K + |D L,σ |ρ n L .
The dual fluxes F n+1 σ ,ε are defined as the average of the fluxes on matching primal faces -the primal faces with coinciding outward normals [START_REF] Herbin | Kinetic energy control in the MAC discretization of compressible Navier-Stokes equations[END_REF][START_REF] Gallouët | An unconditionally stable pressure correction scheme for the compressible barotropic navier-stokes equations[END_REF]. The approximation u n+1 ε,i of the ith component of the velocity valued on the dual face ε is obtained by the upwind scheme. Finally the discrete ith prediction step is given by

(ρ n+1 D σ u n+1 σ -ρ n D σ u n σ )/δt + 1 |D σ | ∑ ε∈E (D σ ) F n+1 σ ,ε u n+1 ε,i -(∇ • (µ n+1 D(u n+1 ))) D σ + (∇(2p n -p n-1 )) D σ = f n+1 σ with |D σ | f n+1 σ = D σ f n+1 dx.
The remaining terms are discretized in a straightforward way, with (∇(2p np n-1 )) D σ approximated by

(|σ |/|D σ |)(ϕ n L -ϕ n K )n K,σ
• e i , if we note ϕ n = 2p np n-1 and define n K,σ as the normal to the face σ outward K. For the viscous term, precautions must be taken given the discontinuous nature of the viscosity. It comes down to the discretization of the following term:

- ∂ D σ µ n+1 D(u n+1 ) : (n ε,σ ⊗ e i )dx
Therefore involving the value of µ ε , the viscosity valued on the faces of D σ . We refer to [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF] for the complete approximation of the viscous term.

ρ * K = ρ K -ν (1) σ + (ρ K -ρ (1) σ + ) -ν (1) σ -(ρ (1) σ --ρ K ) ρ * * K = ρ * K -ν (2) σ + (ρ * K -ρ * σ + (2) ) -ν (2) σ -(ρ * σ - (2) -ρ * K ) ρ * * * K = ρ * * K -ν (3) σ + (ρ * * K -ρ * * σ + (3) ) -ν (3) σ -(ρ * * σ - (3) -ρ * * K )
3 Numerical Test -Dropping a ball in a viscous fluid

We drop a rigid heavy sphere in a viscous fluid and observe it reach its terminal velocity according to the principle that the drag force exerted on the particle by the fluid as well as buoyancy balances the gravity applied to the sphere. We define the fluid using the density and viscosity ρ f = µ f = 1. The sphere with radius r = 0.08 and density ρ s = 100 is falling in the rectangular domain [0, 1] × [0, 1] × [0, 3] to which we applied channel-flow boundary conditions. The gravity constant applied to the ball is g = 9.81. We take for the penalty viscosity µ s = 10 4 . For the time step we will be using δt = 0.001. The spatial step h is such that h = max i=x,y,z h i = 1/50. The initial position for the ball is (0.5, 0.5, 1). In this particular test we resort to the upwind scheme for the approximation of fluxes on the primal faces. However this technique produced a large diffusion of the discontinuous quantities and a deformability of the particle (that is ensured by a well-advected viscosity), thus resulting in an incorrect terminal velocity. The advection of the solid phase requires an antidiffusive scheme as introduced in Section 2.2, and complementing numerical tests are on the way.

Fig. 2

 2 Fig. 2 State of the Velocity component in the Z-direction at times t = 0.01, t = 0.02, t = 0.05, t = 0.07.

Antidiffusive transport scheme for the particles

For the advection of the density we resort on one hand to the classic upwind scheme, and on the other hand to an antidiffusive transport technique we introduce below. As stated in the previous section, the difference in those methods primarily involves the computation of F n+1 K,σ in (3). For the upwind scheme the latter is defined as

The updated density at the face σ = K|L, denoted ρ n+1 σ , is given by:

However in the context of the non-homogeneous Navier-Stokes equations penalized by the H1-penalty method, this approximation of the mass balance generates a large numerical diffusion around the solid phase (as observed in Section 3). We replace the diffusive upwind technique for the transport of the density with an antidiffusive scheme (AD-scheme) based on [START_REF] Bokanowski | Anti-dissipative schemes for advection and application to Hamilton-Jacobi-Bellmann equations[END_REF][START_REF] Després | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF] and adapted to the dimension d = 3 of the problem by considering an alternate directions variant.

Let us focus on the transport of the density by the AD-scheme in the direction e i . For K ∈ M we note ρ * K the updated value of the density on K computed from its previous value ρ K . Let (K -, K + ) ∈ M 2 so that the primal cells K -, K, K + are consecutive and such that σ -= K -|K and σ + = K|K + are in E (i) . We reorder the cells by imposing n K,σ + • e i ≥ 0 and n K,σ -• e i ≤ 0. For the time being, Let us assume that the velocities are positive. The transport of the density in ith direction is carried out by:

where u σ + (resp. u σ -) is the value of the ith component of the velocity on σ + ( resp. σ -). The density on the faces σ + and σ -, noted ρ σ + and ρ σ -, are to be determined. An equivalent formulation yields:

by defining ν σ + = δt|σ |u σ + /|K| and ν σ -= δt|σ |u σ -/|K| as local Courant numbers. We compute the value ρ σ + in such a way that

using the classic minmod formula. We carry out a similar process for any values of u on faces σ + , σ - [START_REF] Bokanowski | Anti-dissipative schemes for advection and application to Hamilton-Jacobi-Bellmann equations[END_REF]. We extend the advection to other directions to obtain the alternate direction variant of the AD-scheme: