Dimitrios N Papadopoulos
email: dimitris.papadopoulos@institutoptique.fr

Alain Pellegrina

Lourdes Patricia Ramirez

Patrick Georges

Frédéric Druon

Broadband high energy diode pumped Yb:KYW multipass amplifier

We report a diode pumped Yb:KYW amplifier delivering up to 27 mJ pulses at a repetition rate of 100 Hz and a spectral bandwidth of 5.5 nm, centered at around 1030 nm. The system is based on a double-head multipass amplifier configuration where the pump thermal load is distributed over two relatively thin laser crystals which permits a sufficiently large number of passes with minimal passive losses thus maximizing the energy extraction efficiency. The amplified pulse bandwidth theoretically supports 340 fs pulses and as a demonstration, a small fraction of the amplified pulses has been compressed to 560 fs.

Ultrafast diode pumped laser systems are of great interest either as direct high intensity systems or as pump sources for OPCPA based systems operating in the extremely short pulse duration regime. Yb-doped materials, due to the excellent Yb 3+ ion's spectroscopic properties as well as their thermal management advantages have been extensively investigated in this general context [1]. The most widely studied material is Yb:YAG as it combines both high gain and high thermal conductivity that has recently resulted in the amplification of picoseconds pulses to the Joule level [2][3][4]. The main limitation however of this material is the available amplification bandwidth, which in the case of high energy systems, restricts its application in the ps regime. Yb:CaF2 on the other hand is a very promising material mainly because of its broad gain bandwidth and high energy storage capacity. Regenerative amplification to several mJ in the kHz regime, both at room temperature and for cryogenically cooled crystals [5,6], as well as Joule level amplification in the femtosecond regime at 10 Hz [7] have been recently demonstrated. The main limitation in this case is the very high saturation fluence of Yb:CaF2, resulting in low gain amplifiers or reduced energy extraction efficiencies. Ideally Yb:CaF2 could be used in the final energetic stages of an amplification chain seeded by a sufficiently high energy, broadband signal to fully exploit its advantages. On the other hand, Yb-doped tungstates (Yb:KGW, Yb:KYW, Yb:KLuW) are well known both for their high gain and broad bandwidth and seem to be the right alternative towards a medium energy level amplification of broad bandwidth pulses. But due to their relatively limited energy storage capacity (with an upper state lifetime of 300 µs) they have been so far restricted to relatively low energy applications at high repetition rates.

In the framework of the Institut de la Lumière Extrême (ILE) project and the 10 PW Apollon source front end development, we study the development of a 100 Hz broadband Yb:KYW multi-mJ level (∆λ>3.5 nm, E>20 mJ), multipass preamplifier for Yb:CaF2 based main amplifiers to reach the required final Joule level output.

Several works have been reported on Yb:KYW or Yb:KGW in both regenerative and multipass configurations. The highest energy achieved so far from a diode pumped system based on a cryogenically cooled Yb:KYW regenerative amplifier, has been limited to 5.5 mJ [8], while pulses as short as 190 fs have been reported for sub-mJ setups [9]. 25 mJ pulses have been also reported but based on a low efficiency and complex Ti:Sapphire pumped Yb:KGW regenerative amplifier operating at 1 Hz [10]. In this study, we present a multipass amplifier based on a diode pumped Yb:KYW crystal operating at room temperature and delivering pulses with energies up to 27 mJ at 100 Hz (2.7 W average power). This is the highest energy ever reported for a Yb:KYW based amplifier. The overall optical-to-optical efficiency of the system is around 15% with an extraction efficiency estimated to be greater than 70%. The achieved bandwidth at the maximum energy is 5.5 nm, corresponding to 340 fs Fourier transform limited (FTL) pulses which, when compared to previously reported results, is an improvement by a factor of at least 2 for this type of amplifier and energy level.

In Figure 1, the general schematic of the complete chirped pulse amplification (CPA) chain studied here is presented. The starting point is a very broadband Ti:Sapphire (Rainbow from Femtolasers) oscillator delivering sub-7 fs pulses at 800 nm with a spectral content that is well extended in the 1030 nm regime. About 4 pJ pulses around 1030 nm are filtered and separated at the output of the oscillator. These pulses are first amplified in a Yd-doped fiber amplifier (Keopsys) to 4 nJ and then stretched in a classical double pass grating stretcher (stretching factor: 0.5 ns/nm). The following amplifier consists of a diode-pumped Yb:KYW regenerative amplifier (Amplitude Systemes) delivering pulses up to 2 mJ. The repetition rate of this amplifier is adjustable up to 1 kHz while the output spectrum can be finely tuned by an intracavity polarizer. In these experiments, the repetition rate was fixed at 100 Hz and the output spectrum was centered at 1030 nm with a bandwidth of 4.5 nm.

In the upper part of Fig. 1, we present a detailed design of the amplifier under study. The amplifier is very compact, covering an area of only 1.15x0.3 m 2 . The whole CPA setup could fit in less than half of an optical table. It is composed of a simple folded 4f imaging cavity where two laser crystals are used as end mirrors. The laser crystals are two identical 0.5º wedged, 5x5mm Yb:KYW disks of 700 µm thickness at b-cut (E//a) and 5% doping level. The crystals are soldered on a metal mount that is cooled by circulating water at about 18 ºC.

The specific cavity design permits both a large number of passes with minimal losses, but most importantly, the combination of two relatively thin amplifying disks that results in the pump thermal load distribution. As shown in Fig. 1, the crystals are pumped in series by a unique fiber coupled diode with a 400 µm core diameter, delivering up to 400 W at 981 nm (LIMO). The diode operates in QCW regime with variable repetition rate and pulse duration, allowing the fine temperature adjustment of the emission wavelength exactly to 981 nm (to match the absorption peak of Yb:KYW), with ~3.7 nm bandwidth. The repetition rate in these experiments is set at 100 Hz and the pulse duration at 450 µs, corresponding to 180 mJ pulse energy (18 W). At maximum pump power, about 135 mJ (75%) are absorbed in both crystals (~86 mJ in the first crystal and 49 mJ in the second one).

Only about 37 mJ is estimated to be the actual total stored energy. This is not surprising since for the relatively high incident pump intensity of more than 22 kW/cm 2 (5 times the saturation intensity of the Yb:KYW for unpolarized pump), the effective storage time is reduced to τsto≈50 µsec (τsto= τfluo/(1+Ip/Isat)). The 450 µsec long pump pulses were chosen however in order to increase the stored energy (with an upper duration limit set by the crystals thermal load) but at the cost of decreasing the global energy storage efficiency.

In the optimal cavity configuration, the beam passes 11 and 10 times through the first and the second crystal respectively (a complete round trip inside the crystal is counted as one pass). Both pump and signal diameter on the crystals planes is 1.5 mm. The cavity is optimized by adjustment of the folding mirror (Mf) position to compensate the induced lensing by the laser disks estimated (at maximum pump power) around 10-12 m. The cavity passive losses have been estimated at ~15% by replacing the two crystals with high reflection mirrors. The amplifier operates in air environment and no vacuum is required. Fig. 2. a) Output energy as a function of the input energy for maximum pump energy (180 mJ). The dashed line curve corresponds to the simulated output (Esto,tot: the total stored energy, Jsto,1/2: the stored fluence for each crystal, g0,1/2: the small signal gain coefficient for each crystal). b) Output energy as a function of the incident pump energy for different input pulse energy at 100 Hz.

In Figure 2a, we plot the the output energy as a function of the signal energy at maximum pump power. We also plot the predicted (dashed line) output energy based on a Frantz-Nodvik model, taking into account the measured pump absorption and passive losses as well as the calculated total stored energy and its relative distribution in each crystal (model parameters are indicated in the figure). For a 2 mJ input, we obtained 27 mJ at the output, well into the saturated, energy extraction regime of operation of our amplifier. In Figure 2b, we plot the output energy as a function of the total incident pump energy for various input pulse energy levels (0.2-2 mJ) at 100 Hz. For 140 mJ pump power, we observe the reduction of the efficiency slope due to the pump absorption saturation and energy storage efficiency reduction. It is also interesting to notice that the output energy barrier of 20 mJ is still surpassed for ~0.5 mJ input (40x gain). At maximum output energy the global efficiency is 15% and the extraction efficiency ~72%.

In Figure 3, the input and output beam profiles of the amplifier are given for increasing pump energy. The M 2 of the beam was optimized for maximum pump power and 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1. 4

E in =2mJ E in =1.5mJ E in =1mJ E in =0.6mJ E in =0.2mJ
measured to be better than 1.2 in both directions (1.1 for the input beam), indicating a high quality output beam.

The main difficulty at this point was the rising astigmatism induced by the laser crystal with increasing pump power [11]. For moderate pump power level (<20W), proper orientation of the crystals (common for both crystals) and respecting the E//a polarization direction, could be beneficially used to compensate the small input beam astigmatism (comparison of beam profiles at Ep=90 and 180 mJ in Fig. 3). Further increase of the repetition rate was possible up to 200 Hz but with a moderate output energy reduction of 10-15% due to the progressive thermal degradation of the beam quality and increasing crystal temperature. Above 200 Hz however the induced asymmetric thermal lens resulted in strong beam diffraction at the edges of the cavity mirrors and the fast decrease of the output energy. The stability of the system has been also examined. At 27 mJ output and over 20 minutes, the pulse-to-pulse rms variation was measured to be 0.8% (0.4% at the input). The beam pointing stability was found <50 µrad.

The evolution of the output spectrum as a function of the output energy is also examined. In Figure 4, we plot the input/output spectra for 1 mJ input energy and for 8 to 22 mJ output energy. We observe the spectrum progressive blue shift towards 1025 nm, resulting in a rectangular-like spectrum of 5.5 nm bandwidth at 22 mJ. This behavior is the combined result of (i) the Yb:KYW (E//a) gain peak around 1025 nm, (ii) the gain saturation in combination to the CPA that counteracts the blue shift and (iii) the specific spectral form of the input pulses centered at 1030 nm, extending to almost 1025 nm. This spectral bandwidth corresponds to about 340 fs FTL pulses. As a demonstration of the compressibility of these pulses, based on a matched grating compressor, a small fraction of the output (few µJ) has been compressed around 560 fs (Fig. 4 inset), about 1.6 times the FTL duration This deviation is most probably due to the complicated higher order spectral phase generated in the Ti:Sapphire oscillator in the 1030 nm regime.

In conclusion, we presented the development of an Yb:KYW multipass amplifier that delivers broadband pulses of 5.5 nm on a record energy level of 27 mJ. The system operates at a high repetition rate of 100 Hz, providing 2.7 Watts of average power with an excellent beam quality. Scaling of the amplifier performance in terms of output energy seems to be possible based on the use of slightly thicker crystals. According to our model, replacing both crystals with 1 mm thick crystals would result in more than 50 mJ output energy. The repetition rate in this case, could be maintained at 100 Hz with a carefully designed cavity and proper cut of the crystal to compensate the astigmatic thermal lens [11].

This work has been funded by the ILE 07-CPER 017-01 contract.

Fig. 1 .

 1 Fig. 1. Experimental setup. Schematic of the CPA injection system (bottom part). Double head Yb:KYW multipass amplifier (upper part).

Fig. 3 .

 3 Fig. 3. Input and output beam profiles (near field) for 90 mJ, 120 mJ and 180 mJ pump energy at 100Hz (Ein=1mJ).

Fig. 4 .

 4 Fig.4. Output spectrum as a function of the output energy for Ein=1mJ (dashed line curve corresponds to the input spectrum). Inset : autocorrelation trace of the recompressed pulses at 22 mJ.


A Jullien

X Chen

A Ricci

J P Rousseau

R Lopez-Martens

L P Ramirez

D Papadopoulos

A Pellegrina

F Druon

P Georges

High-contrast Ultrabroadband Frontend Source for High Intensity Few-Cycle Lasers

Keywords:

An ultrabroadband seed source for high-power, high-contrast OPCPA systems at 800 nm is presented. The source is based on post compression in a hollow-core fiber followed by crossed polarized waves (XPW) filtering and is capable of delivering 80µJ, 5fs, CEP-stable (0.3rad RMS) pulses with excellent spectral and temporal quality.

High temporal quality few-cycle pulses are needed as seeds for high-power OPCPA type systems to enhance the final contrast [START_REF] Tavella | Investigation of the superfluorescence and signal amplification in an ultrabroadband multiterawatt optical parametric chirped pulse amplifier system[END_REF]. Such ultrashort pulses must feature both high incoherent contrast as well as high spectral quality to minimize the relative intensity of the coherent pedestal and satellite pulses. The preservation of the carrier-envelope phase (CEP) value is also required. XPW generation in χ (3) -anisotropic crystals (BaF2) is a well-known technique leading to significant improvement of the contrast of femtosecond pulses as well as remarkable spectral smoothing [START_REF] Jullien | 10-10 Temporal contrast for femtosecond ultraintense lasers by cross-polarized wave generation[END_REF][START_REF] Jullien | Highly-efficient nonlinear filter for femtosecond pulse contrast enhancement and pulse shortening[END_REF]. It has been recently demonstrated that XPW generation preserves the CEP [START_REF] Osvay | Preservation of the carrier envelope phase during cross-polarized wave generation[END_REF]. Furthermore, in our previous work, we have shown that this achromatic third-order nonlinear effect is well suited for processes involving ultra-broadband spectra [START_REF] Jullien | Nonlinear spectral cleaning of sub-10 fs pulses via cross-polarized wave (XPW) generation[END_REF]. In this paper, we demonstrate the efficient temporal cleaning of high energy (550 µJ), 5 fs pulses by the XPW filter. On a day-to-day basis 80 µJ pulses with a duration of 5.5 fs and excellent spectral quality (quasi-Gaussian XPW spectrum over 350 nm) are generated. The temporal characterization of these ultrashort XPW pulses shows, for the first time, such an impressive improvement of their coherent and incoherent contrast. The experimental setup (Fig. 1(a)) includes a commercial 1 kHz CEP-stabilized laser delivering up to 1.5 mJ, 25 fs pulses (Femtolasers, GmbH). To reach the few-cycle regime, about 1.2 mJ, are focused into a hollow-core fiber (HCF) filled with neon and compressed by a set of broadband chirped mirrors and a pair of fused-silica wedges. The obtained 550 µJ, 5 fs compressed pulses are then seeded into the XPW filter. The beam is focused by a f=2 m mirror into a vacuum chamber to avoid nonlinear pulse distortion in air. The 1 mm thick nonlinear crystal (BaF2, [011] crystallographic orientation) is placed out of focus to reach the adequate peak intensity for XPW generation. The XPW signal is then discriminated by a broadband thin film polarizer with low dispersion and good transmission over a broad spectral bandwidth, allowing the compression of the XPW pulse. However, its polarizing

(a) (b)

efficiency is limited so that, with an input polarizer before the HCF, the setup extinction ratio is 10 -3 -10 -2 :1. The XPW pulse is sent to a second chirped mirror compressor for temporal characterization. According to our previous investigations [START_REF] Jullien | Nonlinear spectral cleaning of sub-10 fs pulses via cross-polarized wave (XPW) generation[END_REF], accurate tuning of the input pulse dispersion (± 4 fs 2 ) allows fine optimization of the XPW process. The filter transmission efficiency is then 15% and 80 µJ pulses are generated. Under the specific compression conditions, the fast electronic response of the XPW process enables efficient spectral cleaning: the XPW output spectrum exhibits a smooth, near-Gaussian shape over 350 nm (Fig. 1(b)). This excellent spectral property is a consequence of the temporal quality enhancement occurring during the process. Fig. 2 (a) shows a comparison of the FROG traces measured before and after XPW. For the XPW pulses, temporal side lobes and sharp features typical of few-cycle pulses are decreased and the spectro-temporal energy distribution is balanced and homogeneous. The measured duration of the XPW pulse is 5.5 fs and the coherent background is clearly reduced on the femtosecond time-scale. After calculated additional propagation in glass (~100 µm) to remove the low residual quadratic phase, we found that the main peak contains 91% of the energy after XPW, for only 79% after the HCF. A rough estimation of the improvement of the temporal contrast on a picosecond time-scale is provided by a third-order high-dynamical cross-correlator measurement (Fig. 2(b)). The device is not designed for ultra-broadband spectra and consequently over-estimates the background intensity relative to the main peak. Furthermore, due to the low energy of the XPW signal, the dynamic range of the device is limited to 8 orders of magnitude in the XPW case. However, the measurement indicates a contrast improvement of at least 2 orders of magnitude (consistent with the extinction ratio of the setup) and for the first time, a 5 fs pulse is shown to present a steep rising edge (<10 -7 at t = -2 ps) and a temporal contrast better than 10 -8 (detection limited). Finally, using a home-made f-to-2f interferometer, we verified the CEP stability after these two consecutive nonlinear stages (HCF+XPW) measured to about 0.3 rad RMS. The CEP stability could be further improved by simple protection of the whole experimental setup to avoid air fluctuations.

Fig. 2

 2 Fig. 2 : a) Experimental setup, b) Spectra after HCF (grey area), XPW (red line) and its Gaussian fit (red dotted line).

Fig. 2

 2 Fig. 2 : a) FROG measurements of the compressed few-cycle pulses after HCF and after XPW, b) Temporal contrast of the few-cycle pulses measured by a third-order high-dynamical cross-correlator before and after the XPW filter.

Ackowledgements: The authors gratefully acknowledge P. Monot for the high-dynamical contrast measurements.


Julien Chauchat
email: julien.chauchat@polytech.univ-mrs.fr

Mickael Pailha

Pascale Aussillous

Marc Médale

Elisabeth Guazzelli
email: elisabeth.guazzelli@polytech.univ-mrs.fr

Bedload Transport. Part 1: Two-Phase Model and 3D Numerical Implementation

Keywords: Sediment transport, two-phase model, granular rheology, numerical model

 where the authors investigate experimentally the behaviour of the flowing granular layer. The presented model and numerical results will be used for a detailed comparison in a near future.

I Introduction

The transport of sediment or more generally the transport of particles by a fluid flow is a problem of major importance in geophysical flows such as coastal or river morphodynamic or in industrial flows with the hydrate or sand issues in oil production and granular transport in food or pharmaceutical industries. This problem has been extensively studied in the literature since the middle of the twentieth century but poorly understood actually (Einstein, 1942;[START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF][START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF][START_REF] Yalin | An expression for bed-load transportation[END_REF].

Recently, Ouriemi et al. (2009a) have proposed a two-phase model describing the bed-load transport in laminar flows that allows to incorporate more physics than in previous modelling based on particle flux or erosion deposition approaches. This two-phase model is based on a Newtonian rheology for the fluid phase and a frictional rheology for the particulate phase [START_REF] Forterre | Flows of dense granular media[END_REF] while the fluid-particle interaction is assumed to follow a Darcy law. This approach allows to predict the threshold of motion for the particle phase and give a description of the flow inside the mobile granular layer. Away from the threshold of motion, a simple analytical model for the particle flux is obtained which gives a quite satisfactory description of experimental observations of bed-load transport in pipe flows (Ouriemi et al., 2009a;[START_REF] Pailha | Bed-load Transport. Part 2: The Mobile Granular Layer[END_REF].

Based on this theoretical model we have developed a 3D Finite Element numerical model that allows to simulate bed-load transport in 2D or 3D configurations (Chauchat and Médale, 2010). It is restricted to the cases where the granular bed does not change its shape in the course of time. The viscoplastic behaviour of the granular phase characterized by the existence of a yield stress is implemented using a regularisation technique. This numerical model is used here to investigate the influence of the granular rheology, Coulomb or µ(I) [START_REF] Forterre | Flows of dense granular media[END_REF], on the mobile granular layer characteristics (velocity profiles, thickness of the layer and particle flux) for a wide range of fluid flow rate.

II Two-phase model and Numerical Implementation

The present model is based on [START_REF] Jackson | The dynamics of fluidized particles[END_REF] averaged equations using the closures developed by Ouriemi et al. (2009a). These equations are summarized hereafter in dimensionless form using the following scaling: the length are scaled by H, the channel height, and the stresses are scaled by ∆ρgH, and therefore the time is scaled by η/∆ρgH where ∆ρ = ρ p -ρ f . The general problem is expressed in terms of the solid volume fraction φ, the mixture velocity -→ u m and the particulate velocity -→ u p . In this paper we only present the mixed-fluid model equation ( 1) which is the sum of the fluid and particulate phases momentum equations assuming no-slip velocity between the fluid and the particles ( -→ u p = -→ u m ) (Chauchat and Médale, 2010).

∇.

-→ u m = 0

H 3 d 3 (1 + R ρ ) D -→ u m Dt = -∇p f -∇p p + ρ m -→ g ∆ρ -→ g + ∇. η e + η p η γm (1)
In these equations,

γm = ∇ -→ u m + (∇ -→ u m ) T , R ρ = ρ f /ρ p
represents the density ratio and Ga = d 3 ρ f ∆ρg/η 2 is the Galileo number where d is the particle diameter. The Galileo number is a Reynolds number based on the settling velocity of particles. p f and p p represent the fluid and particulate pressure respectively. The latter is assumed to be hydrostatic and then its vertical gradient is equal to the apparent weight of the particles. To close these equations we need to prescribe the effective fluid viscosity η e = η(1 + 5/2φ) [START_REF] Einstein | Eine neue bestimmung der molekuldimensionen[END_REF]. Following [START_REF] Jop | A constitutive law for dense granular flows[END_REF], the frictional stress can be written as τ p = η p γp , with η p = µ(I)p p / γp where the friction coefficient is expressed as µ(I) = µ s + (µ 2 -µ s ) / (I 0 /I + 1) as a function of the inertial number I = γp η f /(α p p ) [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF][START_REF] Forterre | Flows of dense granular media[END_REF]. The parameter α is the linked to the permeability K which is defined following Kozeny-Carman [START_REF] Happel | Low Reynolds number hydrodynamics[END_REF]. The coefficient µ s corresponds to the static friction and µ 2 to the dynamical friction whereas I 0 is an empirical parameter that have been estimated to I 0 = 1 by [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF]. This formulation of the frictional viscosity is based on the µ(I) rheology but we can assume a simpler rheology by considering a constant friction coefficient µ s corresponding to a Coulomb rheology therefore

K = (1-φ) 3 k CK φ 2 d 2 = α d 2 with k CK ≈ 180
η p = µ s p p / γp .
The main issue for the numerical solution of these equations comes from the divergence of the particulate viscosity as the particulate shear rate tends toward zero (i.e. in the static zone). The basic idea to overcome this issue consists in regularizing the viscosity by adding a small quantity (λ) to the denominator of the particulate viscosity η p = µ s p p /( γp +λ) then the divergence is controlled by this parameter and the viscosity is kept finite. In other words, the static zone in the frictional rheology is replaced by a very viscous one. We have analysed the influence of the parameter λ on the model solution using a Coulomb rheology in Chauchat and Médale (2010), the reader is referred to this article for more details. It is possible to regularize the µ(I) rheology using the same technique:

η p = µ s p p γp +λ + (µ 2 -µ s )p p
λ+ γp +I 0 α p p /η f The present numerical model results from the discretisation by the standard Finite Element Method of the variational formulation associated with equation ( 1). The non-linearities in the governing equations are solved by a Newton-Raphson algorithm by taking the first variation of the variational formulation. In our implementation, we use piecewise quadratic polynomial approximation for the velocity and piecewise linear discontinuous approximation for the pressure. In the computations, we have employed a 27-nodes hexahedra element (H27) for the velocities.

The incompressibility constraint is solved by a penalisation method (penalty parameter set to 10 9 ).

III Results and discussion

We present the results of the previous two-phase numerical model applied to the flow of a Newtonian fluid over a granular bed. The geometry used for the simulations is the same as in [START_REF] Pailha | Bed-load Transport. Part 2: The Mobile Granular Layer[END_REF]. This is a rectangular duct of aspect ratio W/H = 0.538 and fulfilled with a fluid-particles mixture at φ = 0.55 in the lower part (H p = 7/8) and with pure fluid (φ = 0) in the upper part of the duct. The regularisation parameter is set to λ = 10 -6 s -1 and we solved by FEM the mixed-fluid formulation of the two-phase flow model for a 3x40x80 mesh with a requested absolute residual lower than 10 -11 per degree of freedom. Since the vertical plane at the center of the duct (y = 0) is a symmetry plane of the problem we only solve for one half of the domain with appropriate symmetry conditions on the plane (y = 0). Figure 1 shows the mixture velocity profile (u m ) in a cross section of the duct. The red thick solid line represents the position of the "static granular bed". The fluid and the mixture are sheared in both z and y directions inducing an increase in the friction compared with the twodimensional case. Due to this complex shear the use of a three-dimensional model is needed. This figure also shows the influence of the granular rheology (Coulomb or µ(I)) on the shape of the mobile granular layer for a given fluid flow rate (Q f ) . In particular, we observe that a greater longitudinal pressure gradient is needed to reach the same flow rate for the µ(I) rheology than for the Coulomb one. The thickness of the flowing granular layer is thinner with the µ(I) rheology compared to the one predicted with a Coulomb rheology. This is due to the increased frictional viscosity in the µ(I) rheology therefore the driving fluid shear stress at the bed interface is dissipated on a thinner layer in this latter case.

We have conducted a parametric study to analyse the influence of the driving pressure drop on the particulate flow rate thanks to an Arc Length Continuation algorithm. We have validated this algorithm by comparison with the analytical solution of Ouriemi et al. (2009a) for the thickness of the mobile granular layer H p -H c and the particle flux Q p in the 2D case (See Chauchat et al. (2010b) for details). Figure (2a) shows the velocity profiles for a given flow rate (Q f = 1e -4). In the µ(I) case the vertical velocity gradient is very small in a rather thick layer above the fixed bed compared with the Coulomb case. As a consequence, for a given flow rate the particle flux predicted using a Coulomb rheology is greater than the one predicted using the µ(I) rheology. Figure 2b illustrates the influence of the granular rheology on the particle flux Q 2D

p /φ 0 . At a non-dimensionalised flow rate of Q f = 10 -4 a decrease of approximately 65% of the the 2D particle flux is observed with the µ(I) rheology. This sensitivity of the particle flux to the rheology is one of the important question that we would like to answer in the joint theoretical, numerical, and experimental investigation presented in [START_REF] Pailha | Bed-load Transport. Part 2: The Mobile Granular Layer[END_REF].

Figure 1 :

 1 Figure 1: Velocity profile obtained by numerical simulations for a fluid flow rate Q f = 1.79 10 -5 m 3 /s with a) the Coulomb friction rheology (dp/dx = -0.163) and b) the µ(I) rheology (dp/dx = -0.2).

Figure 2 :

 2 Figure 2: Comparison of the numerical results for the experimental 3D configuration using Coulomb friction and µ(I) friction rheology for a) the velocity profiles for both rheology at a flow rate Q f = 1e -4 and b) the particle flux Q 2Dp /φ 0 versus the fluid flow rate. Figure1shows the mixture velocity profile (u m ) in a cross section of the duct. The red thick solid line represents the position of the "static granular bed". The fluid and the mixture are sheared in both z and y directions inducing an increase in the friction compared with the twodimensional case. Due to this complex shear the use of a three-dimensional model is needed. This figure also shows the influence of the granular rheology (Coulomb or µ(I)) on the shape of the mobile granular layer for a given fluid flow rate (Q f ) . In particular, we observe that a greater longitudinal pressure gradient is needed to reach the same flow rate for the µ(I) rheology

Acknowledgement

Funding from the Institut Français du Pétrole and Agence Nationale de la Recherche (Project Dunes ANR-07-3_18-3892) are gratefully acknowledged. This work was performed using HPC resources from GENCI-IDRIS (Grant 2010-96212).


Mickael Pailha
email: mickael.pailha@manchester.ac.uk

Julien Chauchat

Pascale Aussillous

Marc Médale

Elisabeth Guazzelli
email: elisabeth.guazzelli@polytech.univ-mrs.fr

Bedload Transport. Part 2: The mobile granular Layer

Keywords: Sediment transport, two-phase model, granular rheology

come

I Introduction

Recently, a two-phase model describing the bed-load transport in laminar flows have been proposed by Ouriemi et al. (2009a), with a Newtonian rheology for the fluid phase and a frictional rheology for the particulate phase [START_REF] Forterre | Flows of dense granular media[END_REF]). This model can be solve analytically away from the threshold of motion and predict a quadratic evolution for the velocity field inside the flowing granular layer. This model gives a quite satisfactory description of indirect experimental observations of bed-load transport in pipe flows Ouriemi et al. (2009a).

In this part we present a direct experimental investigation of the velocity field inside the flowing granular layer in a rectangular duct partially filled with transparent spherical particles driven by an index-matched fluid.

II Experimental setup

The central part of the experimental setup is a rectangular translucent pipe made of glass, 1 m long, 3.5 cm wide and 6.5 cm height as shown in Fig. 1. The pipe is partially filled with a layer of beads and full of fluid. A gear pump generates the flow with a flow rate up to 7 L/min which allows a recirculation of the fluid with a constant differential pressure. A reservoir contains the fluid which is transfered to the pipe with the pump. At the entrance of the pipe, the fluid pass through a porous media made of 1cm glass beads in order to laminarise. The fluid reaches a poiseuille profile on top of the layer of beads and also flows in the layer. Finally it goes back to the reservoir through a sieve to retain the particules which go out.

In order to measure the velocity profile inside the layer, an index matching method is used which make paricles invisible and allows to see through the layer. The optical index of 1.49 is determined by particles made of pmma, 1 and 2 mm with a polydispercity <5%. The fluid, Triton X-100 is choosen as it has the advantage of matching with the particules. The density of First, the granular media is initially dilated in order to prevent from crystallization by putting the pipe upside down once and letting the beads sediment. The pipe is tilted to bring back the granular media to the entrance of the pipe. Then a slow flow rate is applied to create a flat interface. Once the fluid area equal to 4 mm, a choosen flow rate is applied. The interface decreases until a height at 7/8 filling where we finally record a movie of the enlightened cross section.

III Experimental results and discussion

From the movie, we can measure the velocity profile and the thickness of the mobile particle layer. Then the the particle flux can be decduced by intergating the velocity profile over the whole layer. Fig. 2 show an image of the central section of the pipe for three flow rates. The velocity profiles are superimposed on the images of the flow and appears in black for the granular phase and in white for the fluid phase. The velocity profile of the granular layer doesn't seem to be quadratic as predicted by the simple two dimensional model. The application of a three dimensional numerical model based on this two-phase approach (Ouriemi et al., 2009a) to this experimental conditions is presented in companion paper [START_REF] Chauchat | Bedload Transport. Part 1: Two-Phase Model and 3D Numerical Implementation[END_REF]. The two major quantities that can be measured experimentally and compared with the three dimensional model are the thickness of the mobile particle layer and the particle flux. The comparison between experiments and numerical results will allow to better understand the sediment transport in this peculiar bed-load regime. In particular we would like to confirm the choice of a frictional rheology for the granular phase in this problem.

Figure 1 :

 1 Figure 1: Experimental setup

Acknowledgement

Funding from the Institut Français du Pétrole and Agence Nationale de la Recherche (Project Dunes ANR-07-3_18-3892) are gratefully acknowledged. Figure 2: Velocity profile for the granular phase(black) and the fluid phase (white) for PMMA particles, 2mm in diameter in Triton. The fluid flow rate increases from left to right a.Q = 3.10 -5 , b. Q = 5.10 -5 and c. Q = 9.10 -5 m3/s.


Malika Ouriemi

Julien Chauchat
email: julien.chauchat@polytech.univ-mrs.fr

Pascale Aussillous

Marc Médale

Elisabeth Guazzelli
email: elisabeth.guazzelli@polytech.univ-mrs.fr

Sediment transport and dunes in pipe flow

Keywords: Sediment transport, two-phase model, granular rheology

come

Introduction

Figure 1: Bed-load transport.

When particle beds are submitted to shearing flows, the particles at the surface of the bed can move as soon as hydrodynamic forces acting on them exceed a fraction of their apparent weight. Bed-load refers to the sediment in transport that is carried by intermittent contact with the stream-bed by rolling, sliding, and bouncing, as shown in Fig. 1. This situation occurs in a wide variety of natural phenomena, such as sediment transport in rivers or by air, and in industrial processes, such as hydrate or sand issues in oil production and granular transport in food or pharmaceutical industries. A very common feature that arises is the formation of ripples, i.e. small waves on the bed surface, or of dunes, i.e. larger mounds or ridges.

The widely recognized mechanism for dune or ripple formation is the fluid inertia or more precisely the phaselag between the bottom shear stress and the bed waviness generated by the fluid inertia, see e.g. [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid. Part 1. Steady flow[END_REF] and references therein. In that case, the shear stress, the maxima of which are slightly shifted upstream of the crests, drags the particles from the troughs up to the crests. However, a complete description of the bed instability is still lacking as the coupling between the granular media and the fluid is poorly understood. Empirical and/or phenomenological laws relating the particle flux to the bottom shear stress have been used in the literature.

The situation addressed here is that in which the bedload can be considered as a mobile granular medium where the particles mainly interact through contact forces. In the past decade, advances has been made in the understanding of granular flows. In particular, it has been shown that a simple rheological description in terms of a friction coefficient may be sufficient to capture the major properties of granular flows [START_REF] Gdr Midi | On dense granular flows[END_REF]. This description has been found to be also successful when an interstitial fluid is present [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF]. We have chosen to use this rheology to describe the granular phase.

We have used a two-phase model having a Newtonian rheology for the fluid phase and friction for the particulate phase. The model equations have been solved numerically in one-dimension and analytically in asymptotic cases. This continuum approach is able to provide an onset of motion for the particle phase and a description of the flow of the mobile granular layer. At some distance from threshold, we obtain the simpler analytical result of the particle flux varying cubically with the Shields number. This algebraic formulation seems quite satisfactory for describing experimental observations of bed-load transport in pipe flows (Ouriemi et al. 2009 Part 1). Based on this particle flux, a simple linear stability analysis has been performed to predict the threshold for dune formation. This basic analysis accounts reasonably well for the experimental observations for dune formation (Ouriemi et al. 2009 Part 2).

The two-phase model

We have chosen to use the standard two-phase equations, see e.g. [START_REF] Jackson | The dynamics of fluidized particles[END_REF], and to propose some closures appropriate to the present problem.

The equations of continuity for the fluid and the particle phases are respectively

∂ǫ ∂t + ∂(ǫu f i ) ∂x i = 0, (1)
∂φ ∂t + ∂(φu p i ) ∂x i = 0, (2)
where u f i is the local mean fluid velocity, u p i the local mean particle velocity, φ the particle volume fraction, and ǫ = 1 -φ the void fraction or fraction of space occupied by the fluid.

The momentum equations for the fluid and particle phases are respectively

ρ f D f (ǫu f i ) Dt = ∂σ f ij ∂x j -nf i + ǫρ f g i , (3)
ρ p D p (φu p i ) Dt = ∂σ p ij ∂x j + nf i + φρ p g i , (4)
where g i is the specific gravity force vector, ρ f the fluid density, ρ p the particle density, n the number density (number of particles per unit volume). The force f i represents the average value of the resultant force exerted by the fluid on a particle. The stress tensors σ f ij and σ p ij may be regarded as effective stress tensors associated with the fluid and particle phases, respectively. We need some closure for the interphase force and stresses of the two phases.

The interphase force can be decomposed into a generalized buoyancy force and a force which gathers all the remaining contributions [START_REF] Jackson | The dynamics of fluidized particles[END_REF]. In the present bed-load case, the remaining contribution reduces to the dominant viscous Darcy drag. The effective stress tensor associated with the fluid phase is supposed to be of Newtonian form with an effective viscosity η e (that for simplicity we can take as the Einstein viscosity) while the stress tensor of the particle phase comes only from direct particle-particle interactions and is described by a Coulomb friction model where:

• the tangential stress is proportional to the load (i.e. the particle pressure p p ) when the granular shear rate is positive (i.e. is equal to µp p with a friction coefficient µ which mostly depends upon the particle geometry and which is given by the tangent of the angle of repose),

• the tangential stress is indeterminate when the granular shear rate is zero.

Bed-load transport in pipe flow

In the calculation, we have considered a flat bed of thickness h p consisting of particles having a diameter d. This bed is submitted to a stationary and uniform Poiseuille flow in a two dimensional channel of thickness D, as depicted in Fig. 2. The two-phase equations presented above are shown to reduce to the momentum equation for the mixture (particles + fluid) and the Brinkman equation for the fluid velocity. Calculations of bed-load transport have been performed numerically but also analytically at some distance from threshold where the mobile layer is larger than a grain size. The details of these calculations can be found in Ouriemi et al. (2009 Part 1).

Here we present a much simpler calculation which gives the flavor of the physical mechanisms involved. The mo- mentum equation of the mixture can be written as

τ p (y) + τ f (y) = τ f (h p ) - ∂p f ∂x (h p -y), (5)
where ∂p f /∂x is the constant fluid pressure gradient driving the Poiseuille flow. This equation shows that the mixture shear stress τ f +τ p increases linearly with depth from the surface value due to the horizontal pressure gradient ∂p f /∂x. It also describes the exchange between the shear stress of the fluid phase τ f and that of the solid phase τ p . On the one hand, at the top of the granular bed h p , the particle shear stress is zero and increases inside the bed until it reaches µp p where the granular medium starts to be sheared. The particle pressure is proportional to the apparent weight of the solid phase

p p = φ 0 ∆ρg(h p -y), (6)
and increases inside the bed. The particle shear stress can keep the value µp p until it reaches y = h c inside the bed. On the other hand, at the top of the granular bed h p , the fluid shear stress is equal to τ f (h p ) and goes to zero at h c below which the granular medium is immobile and behaves as a porous medium. Moreover, the Darcy drag term is dominant in the Brinkman equation for the fluid velocity. Therefore, inside the bed, there is very little slip between the two phases and both particle and fluid phases move at the velocity of the mixture with the following parabolic profile, see Fig. 3, The solid lines represent the numerical integration of equation ( 11). The time-scale is (6Re/Ga) 1/2 (D/d)(η e /∆ρgd) where Re is the Reynolds number and Ga the Galileo number Ga = ρ f ∆ρgd 3 /η 2 . The bed-load thickness has a very simple linear variation in Shield number θ = τ f (h p )/∆ρgd (ratio of the shear stress at the top of the bed to the apparent weight of a single particle)

u p ≈ u f ≈ (µφ 0 ∆ρg + ∂p f ∂x ) η e (y -h c ) 2 2 . (7)
h p -h c ≈ θd µφ 0 + ∂p f ∂x /∆ρg . ( 8
)
Assuming that the critical Shields number for the onset of grain motion corresponds to the thickness of the mobile layer being half a particle, one has

θ c ≈ µφ 0 /2. (9)
Finally, the particle flux can be expressed as

q p / ∆ρgd 3 η e = φ 0 θ c 24 θ θ c 3 , ( 10
)
where η e is the effective viscosity of the fluid phase. These predictions have been tested against experimental measurements of bed profile evolution in a pipe flow. The principle of the experiment is (i) to fill the pipe with fluid and particles and to build an uniform flat bed, and (ii) to record the evolution of the bed height as a function of time using a laser sheet illumination for a given flow rate. If we choose to be in the condition of bed-load transport, i.e. above the critical Shields number, the bed height is always seen to decrease with increasing time as the test section is not fed in with particles. The indirect method for measuring the flux of particles is to use the equation for the conservation of particles

φ 0 ∂h p ∂t + dq p dh p ∂h p ∂x = 0, (11)
which can be written as a kinematic wave equation in dimensionless form as a function of the nondimensional Shields number, see details in Ouriemi et al. (2009 Part 1).

There is good agreement between the experimental observations and the theoretical prediction based on the flux of particles found earlier in the new two-phase model. First, the prediction θ c = µφ 0 /2 with φ 0 = 0.55 in the bulk of the bed and µ = 0.43 [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF] agrees well with the experimental value 0.12 ± 0.03 [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. Second, the cubic law for the particle flux seems satisfactory for describing the velocity of the kinematic waves which are triggered at the entrance of the tube, see Fig. 4, as well as the time-evolution of the bed height in conditions of bed-load transport for pipe flows, see Fig. 5.

Dune formation in pipe flow

Different dune patterns are observed as the flow rate is increased from the laminar to the turbulent regimes. Small To predict the small dune formation, we have performed a simple linear stability analysis where inertia in the fluid produces a phase-lag in the shear stress which is destabilizing, see Fig. 9, while the component of gravity down an incline stabilizes the perturbations, see Fig. 10. Following [START_REF] Charru | Phase diagram' of interfacial instabilities in a two-layer Couette flow and mechanism of the long-wave instability[END_REF], we first calculated the perturbed fluid flow over a wavy bottom considered as if fixed. Then, we used the particle flux found by Ouriemi et al. (2009 Part 1) to relate the bed height evolution to the shear stress at the top of the bed through the particle mass conservation. The threshold for dune formation is found to be controlled by the Reynolds number, see further details in Ouriemi et al. (2009 Part 2).

Figure 11 presents the phase diagram of the dune patterns in the plane Re pipe , Ga(h f /d) 2 . We choose this plane to exhibit both the threshold for incipient particle motion controlled by the Shields number and that for dune instability predicted to be controlled by the Reynolds number in the linear stability analysis. In this plane, the predicted threshold for particle motion (θ c = 0.12) is the dashed line. The predicted instability threshold is the horizontal solid line. The dotted line indicates the domain of validity of the algebraic law relating the dimensionless particle flux to the Shields number found by Ouriemi et al. (2009 Part 1) and thus indicates the domain of validity of the instability threshold prediction of Ouriemi et al. (2009 Part 2).

The three regimes of 'no motion', 'flat bed in motion', and 'small dunes' are well delineated by these boundaries in the given limit of validity. Clearly, the threshold for incipient particle motion and that for small dune instability are observed to differ as there is a large region of 'flat bed in motion' without any dune formation. The threshold prediction of the simple linear stability with a single adjustable parameter that we have taken to be realistic, Re c ∝ 1/µ ≈ 37.5 with µ = 0.43 [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF], provides a correct boundary for the 'small dune' instability. Furthermore, the regimes of 'vortex dunes' and 'sinuous dunes' seem separated and their thresholds also well described by the Reynolds number of the pipe as a control parameter.

Conclusions

In summary, we have proposed a two-phase model to describe bed-load transport in the laminar viscous regime, i.e. the flux of particles in a flat mobile bed submitted to laminar flows. The fluid phase has been assumed to be a Newtonian viscous liquid with Einstein dilute viscosity formula applied to the concentrated situation while the particle phase to have Coulomb solid friction with the shear stress proportional to the pressure.

We have applied this continuum model to bed-load transport in pipe flows. The relevant equations have been found to be the Brinkman equation for the fluid phase and the momentum balance equation for the mixture. They have been solved numerically but also analytically. A very simple analytical model, valid when the Poiseuille flow is not significantly perturbed, finds that the bed-load thickness varies linearly with the Shields number whereas the particle flux cubically with it.

We have compared these predictions with experimental observations in pipe flow. The cubic law for the particle flux seems satisfactory for describing the velocity of the kinematic waves which are triggered at the entrance of the tube as well as the time-evolution of the bed height in conditions of bed-load transport for pipe flows.

We have also examined the different dune patterns which are observed when a bed composed of spherical particles is submitted to a shearing flow in a pipe. 'Small dunes' present small amplitudes and only exist in laminar flow. 'Vortex dunes' are characterized by the existence of vortices at their front and are found either in laminar or turbulent flow. 'Sinuous dunes', showing a double periodicity, appear in turbulent flow. While the threshold for incipient motion is determined by the Shields number, that for dune formation seems to be described by the Reynolds number.

To predict the small dune formation, we have performed a simple linear stability analysis. This analysis containing the basic ingredient of the destabilizing fluid inertia and stabilizing gravity is found sufficient to provide realistic predictions.

The two-phase nature of the problem has been only accounted for in the particle conservation equation in the present study and thus a full two-phase analysis should be performed in the future, in particular through full three-dimensional computations as undertaken in Chauchat & Médale (2010) and in the accompanying paper (Chauchat et al. 2010).

Figure 2 :

 2 Figure 2: Sketch of a particle bed submitted to a Poiseuille flow in a two dimensional channel.

7Figure 3 :

 3 Figure 3: Typical parabolic profile predicted inside the bed.

Figure 4 :Figure 5 :

 45 Figure 4: Dimensionless kinematic velocity as a function of the initial Shields number for different combinations of particles and fluids. The solid line shows the prediction c * (θ) ∝ ∂qp ∂θ ∝ θ 7/2 . The vertical dotted line shows the limit of this prediction for large θ.

Figure 6 :

 6 Figure 6: Profile of the small dunes.

Figure 7 :

 7 Figure 7: Side and top views of the vortex dunes.

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Top and side views of the sinuous dunes.

Figure 11 :

 11 Figure 11: Phase diagram of the dune patterns.

th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Acknowledgement

Funding from the Institut Français du Pétrole and Agence Nationale de la Recherche (Project Dunes ANR-07-3_18-3892) are gratefully acknowledged.


Julien Chauchat
email: julien.chauchat@polytech.univ-mrs.fr

Malika Ouriemi

Pascale Aussillous

Marc Médale

Elisabeth Guazzelli
email: elisabeth.guazzelli@polytech.univ-mrs.fr

A 3D Two-Phase Numerical Model for Sediment Transport

Keywords: Sediment transport, two-phase model, granular rheology, numerical model

. We have considered two formulations of the model based on a two fluid or a single mixed-fluid description. The governing equations are discretized by a finite element method and a penalisation method is introduced to cope with the incompressibility constraint whereas a regularisation technique is used to deal with the visco-plastic behaviour of the granular phase. The accuracy and efficiency of the numerical models have been compared with the analytical solution of Ouriemi et al. (2009a). It turns out that one must takes a smaller regularisation parameter (one order of magnitude) in the mixed-model than in the two-fluid one for a comparable accuracy. Using an Arc Length Continuation algorithm coupled with this model we have investigated the evolution of the solution in terms of the height of the flowing granular layer and the particle flux against the longitudinal pressure gradient. The results of these simulations are in good agreement with the analytical solution in the 2D case and three-dimensional computations have been carried out in a square and circular cross-section ducts showing that the effect of the geometry is non-trivial.

Introduction

The transport of sediment or more generally the transport of particles by a fluid flow is a problem of major importance in geophysical flows such as coastal or river morphodynamic or in industrial flows with the hydrate or sand issues in oil production and granular transport in food or pharmaceutical industries. This problem has been extensively studied in the literature since the middle of the twentieth century but poorly understood actually [START_REF] Einstein | Formulas for the transportation of bed load[END_REF][START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF][START_REF] Bagnold | The flow of cohesionless grains in fluids[END_REF][START_REF] Yalin | An expression for bed-load transportation[END_REF].

Recently, Ouriemi et al. (2009a) have proposed a twophase model describing the bed-load transport in laminar flows that allows to incorporate more physics than in previous modelling based on particle flux or erosion deposition approaches. This two-phase model is based on a Newtonian rheology for the fluid phase and a frictional rheology for the particulate phase (µ(I) [START_REF] Forterre | Flows of dense granular media[END_REF] or Coulomb friction) while the fluidparticle interaction is assumed to follow a Darcy law. This approach allows to predict the threshold of motion for the particle phase and give a description of the flow of inside the flowing granular layer. Away from the threshold of motion, a simpler analytical model for the particle flux is obtained which gives a quite satisfactory description of experimental observations of bed-load transport in pipe flows Ouriemi et al. (2009a).

Based on this theoretical model we have developed a 3D numerical model that allows to simulate bed-load transport in 2D or 3D configurations [START_REF] Chauchat | A 3D numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF]. It is restricted to the cases where the granular bed does not change its shape in the course of time, consequently ripples and dunes formation are beyond the scope of this paper. We consider a mixed-fluid and a two-fluid formulation, the models equations are discretized by a finite element method and a penalisation method is used to impose the incompressibility constraint. The granular rheology is analoguous to a viscoplastic behaviour especially by the existence of a yield stress with the particularity that this yield stress depends on the depth inside the granular layer. A regularisation technique is used to deal with the viscoplastic rheology and is shown to give satisfactory results.

Two-phase model

The present model is based on [START_REF] Jackson | The dynamics of fluidized particles[END_REF][START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[END_REF] averaged equations using the closures developed by Ouriemi et al. (2009a). These equations are summarized hereafter in dimensionless form using the following scaling: the length is scaled by H, the channel height (see figure 1), and the stresses is scaled by ∆ρgH, and therefore the time is scaled by η/∆ρgH where ∆ρ = ρ p -ρ f . The problem is expressed in terms of the solid volume fraction φ, the mixture velocity -→ u m and the particulate velocity -→ u p .

Two-fluid model

∇. -→ u m = 0 ∇. -→ u p = 0 Ga H 3 d 3 D -→ u m Dt = -∇p f + ∇. η e η γm - H 2 K -→ u m - -→ u p + ρ f -→ g ∆ρ -→ g Ga H 3 d 3 R ρ φ D -→ u p Dt = -∇p p -φ∇p f +∇. η p η γp + φ∇. η e η γm + (1 -φ)H 2 K -→ u m - -→ u p + φ -→ g -→ g (1) In these equations, γk = ∇ -→ u k + (∇ -→ u k ) T with k = m or k = p, R ρ = ρ f /ρ p
represents the density ratio and Ga = d 3 ρ f ∆ρg/η 2 is the Galileo number where d is the particle diameter. The Galileo number is a Reynolds number based on the settling velocity of particles.

Mixed-fluid model

∇.

-→ u m = 0

Ga H 3 d 3 (1 + R ρ ) D -→ u m Dt = -∇p f -∇p p + ρ m -→ g ∆ρ -→ g +∇. η e + η p η γm
(2) To close these equations we need to prescribe the effective fluid viscosity η e = η(1 + 5/2φ) [START_REF] Einstein | Eine neue bestimmung der molekuldimensionen[END_REF], the particulate viscosity η p = µ s p p / γp [START_REF] Jop | A constitutive law for dense granular flows[END_REF] with µ s the internal friction coefficient and the permeability K = (1-φ) 3 d 2 k CK φ 2 with k CK ≈ 180 [START_REF] Happel | Low Reynolds number hydrodynamics[END_REF].

Numerical model

The Finite Element Method (FEM) is based on the discretisation of the variational formulation associated with equations ( 1) for the two-fluid model and on equation ( 2) for the mixed-fluid model.

The two-fluid formulation is based on the solution of the system (1) in a weakly coupled way. The algorithm associated with the mixed-fluid formulation is based on the solution of the system (2). The non-linearities in the governing equations are solved by a Newton-Raphson algorithm by taking the first variation of the variational formulations.

The specific issue raised by the previous two-phase model lie in the calculation of the frictional stress. Following [START_REF] Jop | A constitutive law for dense granular flows[END_REF]), the frictional stress can be written as

τ p = η p γp ,
with η p = µ s p p / γp . The particulate viscosity diverges as the particulate shear rate tends toward zero (i.e.: in the static zone) raises obvious numerical problems. The basic idea to overcome this issue consists in regularizing the viscosity by adding a small quantity (λ) to the denominator of the particulate viscosity η p = µ s p p /( γp +λ) then the divergence is controlled by this parameter and the viscosity is kept finite. In other words, the static zone in the frictional rheology is replaced by a very viscous zone. In the next section of this paper we will discuss the influence of the value of the parameter λ on the model solution.

In our implementation, we use piecewise quadratic polynomial approximation for the velocity and piecewise linear discontinuous approximation for the pressure. In the computations, we have employed a 27-nodes hexahedra element (H27) for the velocities. The incompressibility constraint is solved by a penalisation method. The code is developed with the PETSc library [START_REF] Balay | PETSc Web page[END_REF]Balay et al. ( , 2004[START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF]) which provides several parallel iterative and direct solvers. As we use a penalisation method to cope with the incompressibility constraint, all the algebraic systems have been solved by the MUMPS direct solver [START_REF] Amestoy | Multifrontal parallel distributed symmetricand unsymmetric solvers[END_REF][START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF] with a penalty parameter set to 10 9 for all the simulations presented in this paper.

Results

We present the results of the previous two-phase numerical model applied to the flow of a Newtonian fluid over a granular bed. First we compare the results of the numerical model with the analytical solution of Ouriemi et al. (2009a) for the bed-load transport in laminar shearing flows to validate quantitatively the two formulations of the twophase flow model. We have also looked at the computational efficiency of the numerical model associated with both formulations.

v f = w f = u p = v p = w p = 0 and ∂u f ∂z = 0 u f = v f = w f = 0 u f = w f = 0 ∂u f ∂x = 0 v f = w f = 0 ∂u f ∂x = 0 v p = w p = 0 ∂u p ∂x = 0 ∂u p ∂x = 0 v p = w p = 0
The sketch of the problem and boundary conditions are given in figure 1. The lower half of the domain is filled with particles at φ = 0.55 immersed in a fluid and the upper part is filled with pure fluid (φ = 0). Therefore in this problem the values of the dimensionless numbers are: Re = 2 10 -2 , Ga = 11, R ρ = 0.4 and d/H = 30. The regularisation parameter is set to λ = 10 -6 s -1 when its value is not mentioned in the figure captions. We solved by FEM the two formulations of the two-phase flow model for a 6x1x40 mesh with a requested absolute residual lower than 10 -11 per degree of freedom.

Figures 2a and2b present the comparison of the horizontal velocity profiles obtained for the two formulations of the two-phase flow model (two-fluid and mixed-fluid) by numerical simulations compared with the analytical solution proposed by Ouriemi et al. (2009a). In both figures the black solid line represents the analytical solution. The good agreement between the numerical solution and the analytical one gives a first qualitaitve validation of the FEM model for the bed-load transport.

The spatial convergence analysis for both formulations (see figure 3) shows that the regularisation of the granular rheology reduces the spatial order of convergence to order one whereas a second order is theoretically hopped. Also, to reach the same accuracy with both formulations we observe that the regularisation for the mixed-fluid model must be decreased by one order of magnitude compared with the two-fluid one.

We now apply the model to three-dimensional configurations: a square and a circular cross-section ducts, with the same values of the dimensionless numbers Re = 2 10 -2 , Ga = 11, R ρ = 0.4, d/H = 30 and Bn = 2 10 4 . We have performed simulations on the half of the domain in the transverse direction for obvious symmetry reasons.

Figure 4 shows the velocity profile in a cross section of the ducts. The contour colors represent the x-velocity of the mixture (u m ). The horizontal thick solid line at z = 0.5 represents the position of the granular bed. The fluid and the mixture are sheared in both z and y directions inducing an increase in the friction compared with the two-dimensional case. Due to this shear increase the velocity is lower than in the two-dimensional case. Figure 5 shows the velocity profiles of the fluid phase velocity in blue and the particulate phase velocity in red (an offset of 10 -3 has been added to make the particulate phase velocity visible) obtained with the two-fluid model. These results illustrate the good behaviour of the numerical model for three-dimensional flow configurations.

: e = 1 N ( N i=1 (U i -U ana i ) 2 ) 1/2
where N is the number of nodes in the mesh. ) the circular cross-section duct (6x896). The fluid phase velocity is in blue and the particulate phase velocity is in red. An offset of 10 -3 has been added to the velocity of the particulate phase (up) to make it visible.

Figure 6 shows the application of the mixed-fluid model coupled with an Arc Length Continuation algorithm that allows us to explore the parameter space in terms of the longitudinal pressure gradient ∂p/∂x. The agreement between the analytical solution of Ouriemi et al. (2009a) and the numerical model for the height of the flowing granular layer H p -H c and the particle flux Q p for a wide range of longitudinal pressure gradient gives another validation of our three dimensional numerical model. As for the velocity profiles we have explored the parameter space for three dimensional configurations that allows us to obtain the height evolution of the flowing granular layer H p -H c and the particle flux Q p in rectangular and circular cross section ducts. These results are presented in Figure 7 and illustrate the non-trivial effect of the geometry on the behaviour of the flowing granular layer. This method applied to complex geometry is a powerful tool to compare the prediction of the initial two-phase model with experiments in real geometry such as rectangular or cylindrical ducts. This is a strong argument for the development of a three dimensional numerical model.

Conclusions

In summary, we have developped a three dimensional numerical model to simulate incompressible particulate twophase flows for the flow of a Newtonian fluid over a granular bed. The incompresibility is imposed by a penalisation method and the viscoplastic characteristic of the granular rheology is dealt with a regularisation technique. We have studied the accuracy of our numerical solution by comparison with the analytical solution derived by Ouriemi et al. (2009a) in function of the spatial discretisation and the value of the regularisation parameter for the mixed-fluid and the two-fluid model. We have concluded that one must take a regularisation parameter one order of magnitude lower for the mixed-fluid model than for the twofluid one to reach the same accuracy.

The Finite Element method allows us to deal with arbitrary geometries and we have shown numerical results of three-dimensional simulation of the bed-load transport in a square and a circular cross-section ducts. This model had also been coupled with an Arc Length Continuation algorithm to describe the evolution of the height of the flowing granular layer and the particle flux with the longitudinal pressure gradient in both 2D and 3D configurations. This method is particularly interesting in the perspective of comparing the prediction of this model with experiments that are always carried out in 3D configurations. Ouriemi et al. (2009b[START_REF] Ouriemi | Sediment transport and dunes in pipe flow[END_REF] have performed a simple linear stability analysis which provides realistic predictions for the formation of small dunes. The aim of the present numerical model is to perform a full threedimensional stability analysis that accounts for the twophase nature of the problem.

7

  th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30 -June 4

Figure 1 :

 1 Figure 1: Sketch of the flow of a Newtonian fluid over a granular bed.

Figure 2 :

 2 Figure 2: Longitudinal velocity profiles for the flow of a Newtonian fluid over a granular bed between two infinite parallel planes obtained by numerical simulations for a) the two-fluid model and b) the mixed-fluid model compared with the analytical solution of Ouriemi et al. (2009a).

Figure 3 :

 3 Figure 3: RMS Error against analytical solution for the flow of a Newtonian fluid over a moving granular bed between two infinite parallel planes: e p and e f stands for the particulate phase and the fluid phase error respectively for the two-fluid model whereas e m designates the mixture velocity error for the mixed-fluid model. The value in brackets is the value of the regularisation parameter λ. The RMS error is defined as: e = 1 N (

Figure 4 :

 4 Figure 4: Velocity profile obtained by numerical simulations with the mixed-fluid model for a) the square cross-section duct (6x20x40) and b) the circular cross-section duct (6x896).

Figure 5 :

 5 Figure5: Velocity profile obtained by numerical simulations with the two-fluid model for a) the square cross-section duct (6x20x40) and b) the circular cross-section duct (6x896). The fluid phase velocity is in blue and the particulate phase velocity is in red. An offset of 10 -3 has been added to the velocity of the particulate phase (up) to make it visible.

Figure 6 :

 6 Figure 6: Comparison of the numerical results for the 2D (mesh size 4x1x640) configurations with the analytical solution of Ouriemi et al. (2009a) in terms of a) the height of the flowing granular layer H p -H c and b) the particle flux Q p versus the longitudinal pressure gradient.

Figure 7 :

 7 Figure 7: Comparison of the numerical results for the 2D and 3D configurations (rectangular and circular cross section ducts) in terms of a) the height of the flowing granular layer H p -H c and b) the particle flux Q p versus the longitudinal pressure gradient. The mesh sizes for the 2D case is 4x1x640, for the rectangular duct (aspect ratio of W/H = 0.2692) is 1x160x320 and for the cylindrical duct is 1x179998.

7

  th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30 -June 4, 2010

th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Acknowledgement

Funding from the Institut Français du Pétrole and Agence Nationale de la Recherche (Project Dunes ANR-07-3_18-3892) are gratefully acknowledged. This work was performed using HPC resources from GENCI-IDRIS (Grant 2010-96212).