M Rébillat

X Boutillon

Identification of honeycomb sandwich properties by high-resolution modal analysis

A method is proposed to identify the mechanical properties of the skin and core materials of honeycomb sandwich. All the elastic coefficients and loss-factors that matter in the dynamics of a panel in the thick-plate approximation are identified. To this end, experimental natural modes (i.e. eigenmodes of the damped system) are compared to the numerical modes of a large sandwich panel (l x,y /h ≃ 80). The chosen generic model for the visco-elastic behaviour of the materials is E(1 + jη). The numerical modes are computed by means of a Rayleigh-Ritz procedure and their dampings are predicted according to the visco-elastic model. The frequencies and dampings of the natural modes of the panel are estimated experimentally by means of a high-resolution modal analysis technique. An optimisation procedure yields the desired coefficients. A sensitivity analysis assess the reliability of the method.

Introduction

Because of their light weight and the easy adjustment of their mechanical properties, honeycomb sandwich panels are widely used nowadays. However, their structural mechanical properties are difficult to predict accurately on the basis of the material properties and identification procedures are often needed. Mixed numerical (Num) / experimental (XP) methods are used to identify the parameters of a model by comparing simulated and measured characteristics (for example modal dampings α Num n vs. α XP n and frequencies f Num n vs. f XP n of the first modes of the system). In order to obtain good identification results, the model parameters must be sensitive to the measured characteristics.

Several authors have addressed the problem of the identification of elastic and damping properties of sandwich panels [START_REF] Bastos | Numerical and experimental approach for identifying elastic parameters in sandwich plates[END_REF][START_REF] Matter | Numerical-experimental identification of the elastic and damping properties in composite plates[END_REF]. The honeycomb sandwich panels raise special difficulties. In order to consider the honeycomb core as homogeneous in the in-plane directions, up to a given frequency f , the corresponding wavelength λ must contain at least 50 cells [START_REF] Burton | Assessment of continuum models for sandwich panel honeycomb cores[END_REF]. For a typical cell side-length s Hexa and height h this implies that the dimensions of a panel must be such that l x,y > λ > 50 s Hexa . On the other hand, the panel must appear as a thick-plate (instead of thin-plate) if the out-of-plane elastic and damping properties are to be identified. For flexural waves this implies that high-enough frequencies are at stake: λ/h ≤ 6 [START_REF] Cremer | Structure-Borne Sound: Structural Vibrations and Sound Radiation at Audio Frequencies Third Edition[END_REF]. In other words, the panel must be large enough and the observed dynamics must include high-enough (but not too high. . . ) modes. Due to the intrinsic dissipations of the materials, the modal characteristics of high modes may be difficult to measure with the Fourier transform (FT) which is limited to modal overlaps of ≃ 30 % in most implementations. The high-resolution modal analysis (HRMA) technique [START_REF] Ege | High-resolution modal analysis[END_REF] is an alternative to the FT for the estimation of modal parameters up to a modal overlap of ≃ 70 %.

In the present work, the identification of most elastic and damping properties of a honeycomb sandwich structure is considered by means of the modal analysis of a large panel (l x,y /h ≃ 80). The HRMA technique is used to estimate modal frequencies and dampings of the first 45 modes of the panel. An optimisation procedure, based on a numerical thick-plate model is used afterwards to identify the corresponding elastic and damping properties.

Mechanical model of the honeycomb sandwich panels 2.1 Hypotheses

The sandwich panel is made of 3 layers: two identical skins and a core (Fig. 1). The thicknesses of the core and the skins are h c and h s respectively. The thickness of the panel is h = h c + 2h s . In the following, "panel" designs the physical structure whereas "plate" refers to the idealised structure made out of the equivalent homogeneous material. The following hypotheses are made on the dynamics of the panel and plate:

• Displacements are small so that the materials and structures behave linearly.

• Only flexural waves are considered.

• The plate is considered to follow the Reissner-Mindlin approximations (thick-plate model).

• The wavelengths include at least 50 cells; according to Burton et al. [START_REF] Burton | Assessment of continuum models for sandwich panel honeycomb cores[END_REF], this ensures that the error on the modal frequencies of the plate (with a homogeneous equivalent core) are less than 2% when compared to those of the panel as computed by various FE-models.

The skin and core materials are each considered as homogeneous, orthotropic in the x and y directions, and viscoelastic.

The formalism chosen for describing the viscoelastic behaviour is that of complex moduli E = E(1 + jη) which do not depend on the frequency (see the model of materials in section 2.2). The Young's and shear moduli and the Poisson coefficient of the core are E c

x , E c y , E c z , G c xy , G c xz , G c yz , ν c yx , ν c xz , ν c yz and ν c xy . The same parameters for the skins are denoted by the s index. The properties of the homogeneous material equivalent to the whole sandwich are denoted by the H index.

The following hypotheses are made on the sandwich panel:

• The sandwich panel is symmetric with respect to its mid-plane.

• Skins are very thin compared to the core so that shear stress in the skin can be ignored:

h s G s xz ≪ h c G c
xz (and the same in the y direction).

• The core is considered to be very soft (E c x ≪ E s x , E c y ≪ E s y and G c xy ≪ G s xy ). Given the generic expression of the moduli of the homogeneous equivalent material

E H = h c h 3 E c +        1 - h c h 3        E s
, this ensures that all in-plane stress in the plate are entirely due to those in the skins. Following these hypotheses, the influence of

E c,s,H z , ν c,s,H xz , ν c,s,H yz , G s xz , G s yz , E c x , E c y , G c xy , ν c xy , ν c yx is
neglected in what follows. These hypotheses are generally fulfilled in common honeycomb sandwich panels. The typical orders of magnitude for the considered parameters in this kind of sandwich panels are:

h s /h c ≃ 10 -1 E c x /E s x ≃ E c y /E s y ≃ G c xy /G s xy ≃ 10 -5
(1)

Model of the materials

The damping of plate vibrations has different origins. In the present study, it is assumed that panels vibrate below their coincidence acoustical frequencies. Consequently, the damping due to acoustical radiation in surrounding air is very small compared to the structural damping. Among the different structural damping models, we have retained the standard hysteretic model (which is frequency-independent, see for example [START_REF] Muravskii | On frequency independent damping[END_REF]). The relationship between the stress ǫ γ and the strain σ γ in each γ-material (γ = s, c, or H) involves 7 complex numbers and can be written, to first order in η as:

σ γ =                     E γ x (1 + jη γ x ) ν γ yx E γ x [1 + j(η γ ν yx + η γ x )] 0 0 0 ν γ xy E γ y [1 + j(η γ ν xy + η γ y )] E γ y (1 + jη γ y ) 0 0 0 0 0 G γ xz (1 + jη γ xz ) 0 0 0 0 0 G γ yz (1 + jη γ yz ) 0 0 0 0 0 G γ xy (1 + jη γ xy )                     ǫ γ (2)
The symmetry of the strain/stress relation leads to the additional relationships

ν γ xy E γ y = ν γ yx E γ x and η γ ν xy + η γ y = η γ ν yx + η γ
x which leaves 12 independent real parameters to be identified for each material (24 altogether). In order to keep a formal symmetry in the mathematical treatment, one defines 1 :

ν γ = ν γ xy ν γ yx η γ ν = η γ ν xy + η γ y (3)

Equivalent thick plate

Under the hypothesis and for the orders of magnitude given in section 2.1, the honeycomb sandwich panel behaves in the low frequency range like a homogeneous thick-plate. The thickness of the plate is chosen to be h. Its mechanical properties are given in Eq. ( 4) and ( 5) as functions of the mechanical and geometrical properties of the skins and the core.

                 E H x = E s x        1 - h c h 3        E H y = E s y        1 - h c h 3        ν H = ν s G H xy = G s xy        1 - h c h 3        G H xz = G c xz G H yz = G c yz (4)                  η H x = η c x E c x E s x h c h 3 + η s x        1 - h c h 3        η H y = η c y E c y E s y h c h 3 + η s y        1 - h c h 3        η H xy = η c xy G c xy G s xy h c h 3 + η s xy        1 - h c h 3        η H xz = η c xz η H yz = η c yz η H ν = η s ν (5)
The 12 independent real parameters

{E H x , η H x , E H y , η H y , G H xy , η H xy , G H xz , η H xz , G H yz , η H yz , ν H , η H ν }
are to be identified. Their knowledge yields the elastic and damping properties of each layer of the honeycomb sandwich panel provided that the 12-equation system formed by Eqs. ( 4) and ( 5) is invertible. A sufficient condition for that is: [START_REF] Muravskii | On frequency independent damping[END_REF] with already

η c x E c x E s x ≪ η s x η c y E c y E s y ≪ η s y η c xy G c xy G s xy ≪ η s xy
E c x E s x ≪ 1, E c y E s y ≪ 1, and G c xy G s xy ≪ 1 (see section 2.1
). This condition is not satisfied only if the η ccoefficients are several orders of magnitude larger than the η s -ones. Since this is not the case here and rarely the case in general 2 , the identification of the E H x , etc . . . yields a measurement of the mechanical properties of the skin and core materials.

Potential, kinetic and dissipated energies in the equivalent thick-plate

To model the dynamical behaviour of the equivalent homogeneous thick plate, the first order Reissner-Mindlin theory [START_REF] Ayorinde | Elastic constants of thick orthotropic composite plates[END_REF] has been chosen. The displacements {u, v, w} in the {x, y, z}-directions respectively are:

u(x, y, z, t) = -zΦ x (x, y, t) v(x, y, z, t) = -zΦ y (x, y, t) w(x, y, z, t) = w 0 (x, y, t) (7)
The potential energy of the plate is:

1 One must keep in mind that η γ ν is not the imaginary part of ν γ . 2 It can be the case when skins are made of metal.

U = 1 2 V (σ H ) T ǫ H dτ = 1 2 S D 1 ∂Φ x ∂x 2 + D 2 ∂Φ x ∂x ∂Φ y ∂y + D 3 ∂Φ y ∂y 2 + D 4        Φ 2 y -2Φ y ∂w 0 ∂y + ∂w 0 ∂y 2        + ... D 5        Φ 2 x -2Φ x ∂w 0 ∂x + ∂w 0 ∂x 2        + D 6        ∂Φ x ∂y 2 + 2 ∂Φ x ∂y ∂Φ y ∂x + ∂Φ y ∂x 2        dxdy (8)
with

D 1 = E H x h 3 12(1 -ν xy ν yx ) D 2 = ν xy E H y h 3 6(1 -ν xy ν yx ) D 3 = E H y h 3 12(1 -ν xy ν yx ) D 4 = 2κ 2 yz hG yz D 5 = 2κ 2 xz hG xz D 6 = G yz h 3 6 (9)
The shear correction factors κ 2 yz and κ 2 xz account for the fact that Eq. ( 7) is an approximation: the Φ x,y coefficients depend lightly on z and sections of the plate do not remain plane in the flexural deformation. The values κ yz = κ xz = 1 have been chosen according to the recommendations of [START_REF] Birman | On the choice of shear correction factor in sandwich structures[END_REF] for sandwich panels.

Based on the material model shown in (section 2.2), the lost fraction of energy during one cycle is given in Eq. ( 10) as:

∆U = -T V (σ H ) T ∂ǫ H ∂t dτ dt = -π S η H x D 1 ∂Φ x ∂x 2 + η H ν D 2 ∂Φ x ∂x ∂Φ y ∂y + η H y D 3 ∂Φ y ∂y 2 + η H yz D 4        Φ 2 y -2Φ y ∂w 0 ∂y + ∂w 0 ∂y 2        + ... η H xz D 5        Φ 2 x -2Φ x ∂w 0 ∂x + ∂w 0 ∂x 2        + η H xy D 6        ∂Φ x ∂y 2 + 2 ∂Φ x ∂y ∂Φ y ∂x + ∂Φ y ∂x 2        dxdy (10)
The kinetic energy T of the system is given in Eq. ( 11) as a function of Φ x , Φ y , and w 0 . In this expression, ρ H is the density of the equivalent homogeneous thick plate which is given by hρ H = h c ρ c + 2h s ρ s .

T = ρ H ω 2 2 (V) u 2 + v 2 + w 2 dτ = ρ H ω 2 2 (S) h 3 12 (Φ 2 x + Φ 2 y ) + hw 2 0 dxdy (11)

Numerical model of the thick plate

In order to compare experimental results to numerical simulations, it is necessary to evaluate the damping factors of numerical modes. The dynamics of the panel is given by the hypotheses listed in section 2.1, the Eqs. (2), and the boundary conditions. Instead of a direct time-integration of the motion, we model here the damping of the numerical modes of the associated conservative system, under the hyptothesis of light damping. The problem consists in evaluating the relationships between the α Num n damping factors and the η H loss-factors.

Modal representation of the system

The honeycomb sandwich panel is considered here as a non conservative system P NC with N degrees of freedom q = {q n }. The damping model presented in section 2.2 corresponds to viscous damping. Under this hypothesis, the equation of the free motion of P NC can be written as:

M q + C q + Kq = 0 (12)
where M, C and K are the mass, damping, and stiffness matrices. In what follows, the modes of P NC are called natural modes. We also refer to the associated conservative system P C corresponding to C = 0, whose modes will be called normal modes ξ n and the normal frequencies f n . If P NC is lightly damped, it can be shown [START_REF] Geradin | Mechanical Vibrations: Theory and Applications to Structural Dynamics[END_REF] that the natural modes are ξ n and the natural frequencies are f n + jα n to first order.

Let U n be the potential energy associated with the n th mode of P NC . It varies in time as exp(-2α n t) so that the energy lost by this mode during one cycle ∆U n is:

∆U n = -2 α n f n U n (13)
Once the mode shapes of P C are known, Eqs. ( 8), [START_REF] Bhat | Natural frequencies of rectangular-plates using characteristic orthogonal polynomials in rayleighritz method[END_REF], and (13) yield the modal dampings α n of P NC .

Rayleigh-Ritz procedure

A Rayleigh-Ritz procedure has been used to derive the mode shapes ξ Num n and the modal frequencies f Num n of P C . To this end, Φ x (x, y), Φ y (x, y), and w 0 (x, y) have been projected on an orthonormal polynomial basis of order Q satisfying the free-free boundary conditions [START_REF] Bhat | Natural frequencies of rectangular-plates using characteristic orthogonal polynomials in rayleighritz method[END_REF]:

Φ x (x, y) = i, j L i j p i (x)p j (y) Φ y (x, y) = i, j M i j p i (x)p j (y) w 0 (x, y) = i, j N i j p i (x)p j (y) (14)
The Hamilton principle is applied, leading to Eq. ( 15). The kinetic and potential energies T and U, defined in section 2.4 are calculated with Eqs. (14).

∀(i, j) ∈ [0, Q -1] 2 : ∂(T -U) ∂L i j = 0 ∂(T -U) ∂M i j = 0 ∂(T -U) ∂N i j = 0 (15)
The above system of 3Q 2 linear equations can be re-written as [K-4π 2 f 2 M]ξ = 0 where f is the eigenfrequency and ξ is the eigenvector of unknown coefficients L i j , M i j and N i j . The resolution of this eigenvalue problem gives a straightforward access to the modal frequencies f Num n and mode shapes ξ Num n of P C .

Derivation of α Num n

By introducing the numerical mode shapes ξ Num n and frequencies f Num n found in section 3.2 in the energies expressions of section 2.4, the relations Eqs. ( 16) are obtained. The coefficients t n and u nk depend only on the geometry and mass parameters of the plate and on the modal shape ξ Num n . For the subscripts of η, {x, ν, y, yz, xz, xy} have been replaced by {1, 2, 3, 4, 5, 6}.

∀n ∈ [1, N] : T n = 4π 2 f 2 n t n U n = 6 k=1 D k u nk ∆U = -π 6 k=1 η k D k u nk (16)
Using relations Eqs. ( 16), the expression Eq. ( 17) of the modal dampings α Num n can be deduced from Eq. (13) using the fact that T n = U n for P C . One can notice than α n is a linear combination of the η k .

α n = f n ∆U n 2T n = 6 k=1 η k D k u nk 4π f n t n (17)
4 Experimental study of a honeycomb sandwich panel

Experimental setup

A rectangular honeycomb sandwich panel (Nomex R honeycombs core and paper skins) is studied experimentally. The known parameters of the panel are given in table 1. Since l x,y /h ≃ 80, the panel is considered to satisfy the Reissner-Mindlin approximations. The panel is suspended by thin wires in order to ensure free-free boundary conditions. It is acoustically excited by an electro-dynamical loudspeaker placed in its vicinity and driven by a wide-band electrical signal. The panel response was measured with a laser vibrometer pointing in one corner. This ensures that all modes are present in the response. By means of a specially designed excitation signal, the impulse response of the panel can be reconstructed [11].

l x l y h s h c
s Hexa ρ c ρ s 39.15 cm 59.10 cm 0.2 mm 4.88 mm 4 mm 37.8 kg/m 3 713 kg/m 3 Table 1: Geometry and mass of the honeycomb sandwich panel. The side-length of the hexagonal core-cells is s Hexa .

High resolution modal analysis

The impulse response of the non-conservative system can be expressed as a summation over its natural modes:

h(t) = N n=1 ξ n exp( j2π f n t -α n t + jφ n ) (18)
In order to extract their experimental modal frequencies f XP n and dampings α XP n , a recently developed modal analysis method [START_REF] Ege | High-resolution modal analysis[END_REF] has been applied to velocity impulse responses of honeycomb sandwich panels as obtained in section 4.1. In the available noise conditions, the parameters of the 45 first modes could be extracted. The modal damping of the highest modes was ≃ 50 %, which is out of reach of traditional implementations of the Fourier transform, hence our use of the new method.

Using several bandpass filters associated with the ESPRIT and ESTER algorithms (see reference [START_REF] Ege | High-resolution modal analysis[END_REF] for details), it is shown that this method yields a precise estimation of f n and α n in presence of moderate noise: the modal frequencies f n can be estimated with a precision of ≃ 0.01 % and the modal dampings with a precision of ≃ 1 %. Moreover, this method allows for the identification of modal parameters of modes having a modal overlap up to 70 %. 19)) and an optimisation procedure based on the gradient-method has been implemented.

Optimization procedure

Estimation method

C E = N n=1 f XP n -f Num n f XP n 2 (19)
It has been shown in section 3.3 that the damping coefficients {α Num n } n∈ [1,N] can be expressed as linear combinations of the {η H

x , η H y , η H xy , η H xz , η H yz , η H ν } loss factors. Therefore, the latter can be obtained by a simple least mean square method, with the constraint that loss factors remain positive.

Results

The optimisation is performed on the 45 first modal frequencies and dampings. The numerical model used a Q = 14-order basis which proved to ensure the convergence of the higher modes values. The identified viscoelastic parameters of the equivalent homogeneous plate are given in Tab. 2. The relative errors in modal frequencies and dampings are shown in Fig. 2. It can be seen that the agreement is very good for modal frequencies (mean absolute error of 1.8 %). The predicted modal dampings fit well the mean measured ones, but the difference is more important (mean absolute error of 10.2 %).

Sensitivity analysis

The sensitivities of the frequency values f n to the coefficients {E H

x , E H y , G H xy , G H xz , G H yz } are defined as S f n (X) = ∂ f n ∂X f n X -1
. They reflect the information contained in a modal frequency relatively to the elastic parameter X.

Results are presented in Fig. 3. Since the modal frequencies are very little sensitive to the Poisson coefficients, their sensitivities to these parameters have not been represented. As expected, it can be seen in Fig. 3a that modes of the form (0, i) or ( j, 0) convey a lot of information relatively to E x and E y respectively. Since the thick-plate model differs from the thin-plate model for the higher frequencies, it is normal that there is almost 10 times more information relative to G xz and to G yz in the higher modes than in the lower ones (Fig. 3b). The lower sensitivity of G xz to the modal frequencies than that of G yz is simply due to the aspect ratio of the plate (l x < l y ).

Conclusion

This identification method yields all the mechanical parameters of the sandwich materials that matter dynamically, under only mild hypotheses. Compared to the method proposed by ref. [START_REF] Bastos | Numerical and experimental approach for identifying elastic parameters in sandwich plates[END_REF], it yields loss factors. Compared to the method proposed by ref. [START_REF] Matter | Numerical-experimental identification of the elastic and damping properties in composite plates[END_REF], it is considerably easier and faster to implement: only one vibrating point is measured, no FEM is needed; it also reaches frequency domains that are usually out of reach of the modal analyses based on the Fourier transform. Incidentally, the method presented here could be used to access the frequency-dependence of the loss factors by considering only modes in a given frequency range. Compared to the static investigations on each sandwich component, this dynamical method is non-destructive and the experimental test needs very little time. Avoiding heavy lab-equipment, it is a good candidate for industrial in-line process of quality control.

Figure 1 :

 1 Figure 1: Geometry of the sandwich plate.

Figure 2 :

 2 Figure 2: Comparisons between measured and predicted modal frequencies (left) and dampings (right).

Figure 3 :

 3 Figure 3: Sensitivities of the mechanical parameters to the modal frequencies. In the figures, each column is the sum of the sensitivities of the different involved parameters. On the top of the column the number of nodal lines in the x (top) and y (bottom) directions of this mode is specified .

  XP n . Since the modal frequencies of the conservative and the real systems are equal to first order (section 3.1), it is valid to find separately and successively the elastic constants and the loss factors.To first order, the modal frequencies depend only on the elastic constants of the homogeneous equivalent thickplate model {E H x , E H y , G H xy , G H xz , G H yz , ν H }. Since this dependence is non-linear, a cost function C E is defined (Eq. (

	This section describes how to derive {E H x , η H x , E H y , η H y , G H xy , η H xy , G H xz , η H xz , G H yz , η H yz , ν H , η H ν } from the experimental val-
	ues of the modal frequencies f XP n and dampings α

Table 2 :

 2 Identified parameters of the homogenised model. The coefficients in the two last columns are mutually related by the symmetry relationships (cf. section 2.2).

		E H x	E H y	G H xy	G H xz	G H yz	ν H xy	ν H yx
	Real part	1.0 GPa 1.4 GPa 0.46 GPa 12 MPa 26 MPa 0.23 0.33
	Loss factor (%)	1.3	1.4	1.1	4.4	8.1	0	0.1


Marc Rébillat

Romain Hennequin

Etienne Corteel

Brian F G Katz

Étienne Corteel

Prediction of harmonic distortion generated by electro-dynamic loudspeakers using cascade of Hammerstein models

Introduction

Even if great efforts have been done to reduce or control their nonlinear behavior, loudspeakers are still the least linear part of the audio chain. The study of their non-linearities is thus of great importance. In electrodynamic loudspeakers, a motor converts the electrical signal in motion and makes a cone vibrate. The piston-like movement of the cone generates the sound field. The motor induces non-linearities because of non-uniform magnetic field, Eddy currents and variations of the electrical inductance with displacements [START_REF] Klippel | Tutorial: Loudspeaker nonlinearities -causes, parameters, symptoms[END_REF][START_REF] Colloms | Distortion mechanisms of distributed mode (DM) panel loudspeakers[END_REF]. Significant excursion induces nonlinear bending in the cone and a nonlinear behavior of the suspensions [START_REF] Quaegebeur | Nonlinear vibrations of loudspeaker-like structures[END_REF].

Traditional nonlinear measurement methods [START_REF] Czerwinski | Multitone testing of sound system components some results and conclusions -part 2: Modeling and application[END_REF] give total harmonic distortion (THD), harmonic distortion of order n (HD n ) or inter-modulation products (IMP) generated by loudspeakers. These quantities are measured using pure tones at a given amplitude and frequency. They do not describe non-linearities themselves but only some of their effects for these arbitrary excitations. Moreover, experimental processes involved in those methods are very time-consuming if a wide range of frequencies and amplitudes is to be covered. There is thus a real need for rapid model based procedures to measure non-linearities.

In order to represent the nonlinear behavior of electro-dynamical loudspeakers, different physical models have already been built. Their formulation was either completely analytical [START_REF] Kaizer | Modeling of the nonlinear response of an electrodynamic loudspeaker by a volterra series expansion[END_REF][START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF] or based on the finite element method [START_REF] Tsuchiya | Finite element simulation of non-linear acoustic generation in a horn loudspeaker[END_REF]. Klippel [START_REF] Klippel | Nonlinear system identification for horn loudspeakers[END_REF] proposed to use a Volterra series expansion [START_REF] Hasler | Phénomènes non linéaires, Chapitre 3 : Séries de Volterra[END_REF] and identified it from measurements. However, all these physical models remains complex and heavy to manipulate. In Refs. [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Hawksford | System measurement and identification using pseudorandom filtered noise and music sequences[END_REF], it is thus suggested to use simpler nonlinear models, cascade of Hammerstein models, to represent and identify different audio systems, including acoustical transducers. Cascade of Hammerstein models [START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF] are a subclass of Volterra models and represents exactly systems having diagonal Volterra Kernels. These models are composed of N branches in parallel, made of one static polynomial non-linearity followed by one linear filter (see Fig. 1).

Some methods exist to identify the structural elements of a cascade of Hammerstein models. Gallman [START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF] and Hawksford [START_REF] Hawksford | System measurement and identification using pseudorandom filtered noise and music sequences[END_REF] proposed a method using Gaussian noise at different amplitudes as inputs. The estimation procedures are strongly based on the knowledge of the order of non-linearity of the polynomial expansion, which is unknown in practical cases. Farina [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF] proposed another method using sine sweeps with frequency varying exponentially with time. Only an upper bound of the order of non-linearity of the model has to be assumed. However this method permits the separation of the contributions of the different orders of non-linearity and not the complete identification of the Kernels of the system. Recently, Novák et al. [START_REF] Novák | Modeling of nonlinear audio systems using swept-sine signals: application to audio effects[END_REF] have identified Kernels from the contributions of the different orders of non-linearity using a least mean square minimization procedure. Of the examples given in Ref. [START_REF] Novák | Modeling of nonlinear audio systems using swept-sine signals: application to audio effects[END_REF], this method does not seem to be very accurate. An extension of the method of Farina [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF] is proposed here. It yields directly all the Kernels of a cascade of Hammerstein models from the contributions of the different orders of non-linearity obtained as in Ref. [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF].

In this paper, cascade of Hammerstein models corresponding to an electro-dynamical loudspeaker are identified using the previously introduced method from a single measurement. With the identified models, total harmonic distortion (THD) and n th -order harmonic distortion (HD n ) are predicted for any desired amplitude and frequency. The THD and HD n at different frequencies and amplitudes are also evaluated by classical means. The agreement between THD and HD n predicted by the identified cascade of Hammerstein models and measured classically is very good in the entire audio frequency range, for a wide range of amplitudes. This method is thus a very reliable manner to quickly characterize non-linearities generated by electro-dynamic loudspeakers.

Mathematical foundations of the method

Mathematical foundations of the method used for direct identification of the elements of a cascade of Hammerstein models are given in this section. This method is an extension of the procedure initially proposed by Farina [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF].

Cascade of Hammerstein models

A cascade of Hammerstein models is a nonlinear model made of N parallel branches. Each branch is composed of one nonlinear static polynomial element followed by a linear one h n (t) [START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF], as shown in Fig. 1. Mathematically, the relation between the input e(t) and the output s(t) of such a system is given by Eq. ( 1), where * denotes the convolution.

s(t) = N n=1 h n * e n (t) (1)
In this model, each impulse response h n (t) is convolved with the input signal elevated to its n th power and the output s(t) is the sum of these convolutions. The first impulse response h 1 (t) represents the linear response of the system. The other impulse responses {h n (t)} n∈{2...N } model the non-linearities. All impulse responses are assumed to be zero-mean. The family {h n (t)} n∈{1...N } will be referred to as the Kernels of the model. Any cascade of Hammerstein models is fully represented by its Kernels.

A cascade of Hammerstein models fed with sine sweeps

To experimentally cover the frequency range on which the system under study is to be identified, sines with time-varying frequencies are interesting signals:

∀t ∈ R e(t) = cos [Φ(t)] (2)
If e(t) is the input of the cascade of Hammerstein models, the output of the nonlinear block of the i th branch (see Fig. 1) is:

e i (t) = cos i [Φ(t)] (3)
Using Chebyshev polynomials, e i (t) is rewritten in Eq. ( 4) as a linear function of {cos[kΦ(t)]} k∈ [1,i] . The coefficient of the matrix C of order 8 are given by Eq. ( 15)

∀i ∈ {1 . . . N } cos i [Φ(t)] = i k=0 C(i, k) cos[kΦ(t)] (4)

Exponential sine sweeps

When the instantaneous frequency of e(t) is increasing exponentially from

f 1 to f 2 (f 1 , f 2 > 0) in a time
T , such a signal is referred to as an "exponential sine sweep" [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Norcross | A survey on the effects on nonlinearity on various types of transfertfunction measurements[END_REF] and its instantaneous phase is given by :

∀t ∈ R Φ(t) = 2π f 1 T ln f2 f1 (e t T ln f 2 f 1 -1) -π/2 (5)
The corresponding instantaneous frequency of e(t) is :

∀t ∈ R f (t) = Φ ′ (t) 2π = f 1 e t T ln f 2 f 1 (6) Thus f (0) = f 1 and f (T ) = f 2 . The frequency range [f 1 , f 2 ]
corresponds to the band of interest of the system under test.

Fundamental phase property

From Eq. ( 5), it can be shown that this type of signal exhibits the phase property given by Eq. [START_REF] Tsuchiya | Finite element simulation of non-linear acoustic generation in a horn loudspeaker[END_REF].

∀k ∈ N * , ∀t ∈ R : kΦ(t) = Φ(t + T ln k ln f2 f1 ) -(k -1)( π 2 + 2πf 1 T ln f2 f1 ) (7)
By choosing

T m = (2mπ - π 2 ) ln f 2 /f 1 2πf 1
with m ∈ N * , the second term in Eq. ( 4) becomes a multiple of 2π and one obtains Eq. ( 8) which is another way to express the k th term of the linearization given in Eq. ( 4).

∀k ∈ N * , cos [kΦ(t)] = cos [Φ(t + ∆t k )] (8a) with ∆t k = T m ln k ln f 2 /f 1 (8b)
For any T m -long logarithmic sweep, multiplying the phase by a factor k results in the same signal, but in advance in the time domain by ∆t k . As can be seen from Eqs. (8a) and (8b), this time advance depends only on the sweep parameters T m , f 1 , f 2 and on k. In Refs. [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Morise | Warped-TSP: An acoustic measurement signal robust to background noise and harmonic distortion[END_REF], similar time advances were obtained using different arguments.

The fact that T m must take only a discrete set of values to ensure the fundamental phase property (8a) has been first highlighted in Ref. [START_REF] Novák | Modeling of nonlinear audio systems using swept-sine signals: application to audio effects[END_REF] but is mathematically demonstrated here for the first time.

One should note that e(t) has been designed for all t with its instantaneous frequency increasing from f 1 to f 2 between t = 0 and t = T . In practice, signals are defined only on [0, T ]. Thus the phase property is not valid on the whole support of the function. The phase property becomes false when the instantaneous frequency of cos [kΦ(t)] is outside the frequency range of interest (i.e. [f 1 , f 2 ]).

Inverse convolution

Using the Eqs. (8a) and (1), one obtains:

s(t) = N n=1 g n * e(t + ∆t n ) with g n (t) = N k=n C(k, n)h k (t) (9)
In order to separately identify each Kernel of the cascade of Hammerstein models, a signal y(t) which looks like an inverse in the convolution sense of e(t) is needed. Unfortunately, such an inverse does not necessarily exist mathematically. A band-limited inverse y(t) can however easily be defined such that it satisfies:

y(t) * e(t) = sinc(2f 2 t) -sinc(2f 1 t) = d(t) (10)
with sinc(x) = sin(πx)/πx. d(t) can be seen as a band-limited Dirac Function, since its Fourier transform is

1 [-f2,-f1]∩[f1,f2] (f ).
Then, Y (f ), the Fourier transform of the inverse filter y(t) can be built in the frequency domain using Eq. [START_REF] Hawksford | System measurement and identification using pseudorandom filtered noise and music sequences[END_REF].

Y (f ) = 1 E(f ) 1 [-f2,-f1]∩[f1,f2] (f ) ≃ 1 E(f ) 1 1 + ǫ(f ) (11)
In practice, the filter Y (f ) should be build by replacing the discontinuous function

1 [-f2,-f1]∩[f1,f2] (f ) by a function of the type 1 1 + ǫ(f )
. This function ensures a smoother transition between the two frequency domains and thus generates less unwanted side effects in the time domain. ǫ(f ) is a frequency dependent real parameter chosen as equal to 0 in the bandwidth and as having a large value outside of it, with a continuous transition between the two domains.

In so doing, y(t) can be considered as an inverse of e(t) in the sense of convolution in the frequency range [f 1 , f 2 ].

Kernels Identification in the temporal domain

After convolving the output of the cascade of Hammerstein models s(t) given in Eq. ( 9) with y(t), one obtains :

y * s(t) = N i=1 g i (t + ∆t n ) * d(t) (12)
Since it is assumed that the system under study has no significant behavior outside [f 1 , f 2 ], g i (t + ∆ t n ) * d(t) reduces to g i (t + ∆t n ), and then :

y * s(t) = N i=1 g i (t + ∆t n ) ( 13
)
In Fig. 2, y * s(t) is represented. Because ∆t n ∝ ln(n) and f 2 > f 1 , the higher the order of linearity n is, the more in advance the corresponding g n (t) will be. Thus, if T m is long enough, the different g n (t) will not overlap. It is then easy to separate them by windowing in the time domain. The separation of the contribution of the different orders of non-linearity by using exponential sweeps, which is mathematically demonstrated here, is already experimentally well known in the audio community [START_REF] Norcross | A survey on the effects on nonlinearity on various types of transfertfunction measurements[END_REF][START_REF] Craven | Practical adaptative room and loudspeaker equaliser for hi-fi use[END_REF][START_REF] Griesinger | Beyond MLS -occupied hall measurerement with FFT techniques[END_REF].

After that step, using Eq. ( 14), the family {h n (t)} n∈ [1,N ] of the Kernels of the cascade of Hammerstein models can be fully extracted.

   h 1 (t) . . . h N (t)    = A T c    g 1 (t) . . . g N (t)    (14)
Figure 2: Separation of the different order of non-linearity after convolution with y(t).

The matrix A of order 8, which is sufficient for practical use, is given in Eq. ( 15). A T c is the matrix A T without the first column and the first row, and (.) T stands for matrix transposition. Then, as all the Kernels of the system are known, the system under study is completely identified.

A = C -1 =    

Practical implementation

In this section, the practical discrete-time implementation of the method presented in section 2 is described.

Overview of the method

In Fig. 3, a global overview of the procedure is given. It can be decomposed in the following steps:

1. Design of the input sweep e(t) using Eq. ( 5). The choice of f 1 , f 2 and T is discussed in section 3.2.

2. Playing e(t) and recording s(t). The sampling frequency f s must be chosen to avoid any aliasing effects caused by the digital to analog converter in the frequency range of interest [f 1 , f 2 ].

3. Generation of the inverse filter y(t) according to Eq. [START_REF] Hawksford | System measurement and identification using pseudorandom filtered noise and music sequences[END_REF]. A convenient way to implement this filter is described in section 2.5.

4. Convolution of the output signal s(t) with the inverse filter y(t) as in Eq. ( 13). This can be done in the frequency domain with a sufficient number of points to avoid temporal aliasing.

5.

Windowing in the temporal domain (cf. Fig. 2) to obtain the {g k (t + ∆t k )} k∈ [1,N ] . Rectangular windows of constant size τ = 1/2(∆t N -∆t N -1 ) can be chosen to separate the different orders of non-linearity. N is the highest desired order in the cascade Hammerstein model.

6. Temporal shift of the different orders of non-linearity to recover {g k (t)} k∈ [1,N ] . A shift of a noninteger number of samples can be performed with a phase shift in the frequency domain.

7.

Multiplication with A T c to access {h k (t)} k∈[1,N ] , according to Eq. ( 14). The matrix A of order 9 is given in Eq. (15).

Choice of the parameters (f 1 , f 2 , T and N )

For satisfactory measurements, the sweep parameters f 1 , f 2 , T and N must be well chosen. These choices must be made considering the following aspects: • The frequencies f 1 and f 2 must be chosen such that the interesting behavior of the system under study is in the frequency range [f 1 , f 2 ].

• The influence of noise on the identification results should also be minimized [START_REF] Morise | Warped-TSP: An acoustic measurement signal robust to background noise and harmonic distortion[END_REF]. By itself, the exponential sweep rejects correctly uncorrelated noise in quiet environment [START_REF] Stan | Comparison of different impulse response measurement techniques[END_REF]. Moreover its energy repartition in frequency is often adapted to the ambient noise [START_REF] Griesinger | Beyond MLS -occupied hall measurerement with FFT techniques[END_REF][START_REF] Muller | Transfer-function measurement with sweeps[END_REF]. The choice of this signal is thus interesting from this point of view. If an excellent signal to noise ratio (SNR) is needed, the longer T , the better the SNR after step 4 will be at a given amplitude of the input signal.

• The number of points to be convolved at step 4 is limited by the available computational power. Thus, as T increases, the computation time will increase. T must not be too high in order to avoid long calculation time.

• N should not be underestimated to guarantee good accuracy in identification. The optimal N is reached when it is impossible to extract the corresponding N th impulse response from the background noise.

• The different pseudo-impulse responses g k (t) which appear in the temporal domain after the convolution with the inverse signal (step 4, see Fig. 2) must not overlap each other. The global decay time of the system τ global is an upper bound of the decay times of each orders of non-linearity. Parameters f 1 , f 2 and T such that ∆t N -∆t N -1 > τ global will thus avoid overlapping of the different orders of non-linearity [START_REF] Muller | Transfer-function measurement with sweeps[END_REF][START_REF] Majdak | Multiple exponential sweep method for fast measurement of head-related transfer functions[END_REF]. Because ∆t N = T m ln N ln f 2 /f 1 , considering the chosen value for N , T must be chosen to be long enough and f 2 /f 1 not so large in order to respect the previous condition.

Modeling electro-dynamic loudspeakers with cascade of Hammerstein models

In this section, an electrodynamic loudspeaker is identified as cascade of Hammerstein models using the method presented in section 2.

Experimental setup

Experiments have been conducted on a standard electrodynamic loudspeaker. All measurements have been made on axis at one meter from the motor with a microphone. Measurements have been done in a semi-anechoic room. In these conditions, a signal to noise ratio of ≃ 50 dB below 50 Hz and of ≃ 70 dB above has been reached for pressure amplitudes of ≃ 110 dB.

Measured cascade of Hammerstein models Kernels

The Kernels corresponding to the system have been measured using the previously described experimental setup. As the cascade of Hammerstein models is a nonlinear model, its Kernels should be independent of the amplitude of measurement. To assess this, measurements of the Kernels have been done using parameters given in Tab. 1 for 10 different amplitudes ranging from 95 to 115 dB in pressure.

Parameter Value f 1 20 Hz f 2 20 kHz f s 192 kHz T 15 s N 5
Table 1: Parameters chosen to identify the simulated system.

Mean measured Kernels and their corresponding standard deviation are given in Fig. 4 for the electrodynamical loudspeaker. The amplitude of the different Kernels of order n ≥ 2 decreases with frequency, which is consistent with the physical analysis of [START_REF] Klippel | Modeling the nonlinearities in horn loudspeakers[END_REF]. The linear part, which is by definition independent of the amplitude, exhibits no variations among the different measurements. The non-linearities have thus been removed successfully from the linear part using the proposed method. The identified Kernels of order n ≥ 2 depends slightly on the amplitude at which they have been measured. As a consequence, the assumption that electro-dynamic loudspeakers can be represented by cascade of Hammerstein models is a correct approximation in the chosen range of amplitudes.

Prediction of the harmonic distortion generated by transducers

In this section, the previously identified cascade of Hammerstein models will be used to predict the harmonic distortion generated by the electro-dynamic loudspeaker.

Link between HD n , THD and cascade of Hammerstein models

To characterize distortion generated by an acoustic transducer, the following approach is classically adopted. The input of the system is assumed to be sinusoidal, and non-linearities generate harmonic components at frequencies higher than the input fundamental frequency. The amplitudes of these harmonics compared to the amplitude of the fundamental are considered as representative of the non-linearity of the transducer. Total harmonic distortion (THD) and harmonic distortion of order n (HD n ) are common tools to quantify this [START_REF] Czerwinski | Multitone testing of sound system components some results and conclusions -part 2: Modeling and application[END_REF]. The THD is the square root of the ratio of the power contained in the harmonics to the power contained in the fundamental. The HD n is the same but for the n th harmonic only.

For a sinusoidal input signal x(t) = X cos(2πf t) which enters a cascade of Hammerstein models identified at the amplitude X 0 , the output signal z(t) can be written as in Eq. ( 16) by using Eqs. ( 1) and [START_REF] Czerwinski | Multitone testing of sound system components some results and conclusions -part 2: Modeling and application[END_REF].

z(t) = N n=1 |Γ n (X, nf )| cos [2πnf t + ∠(Γ n (X, nf ))] with Γ n (X, f ) = N k=1 X X 0 k C(k, n)H k (f ) (16)
THD and HD n can thus be directly identified from Eq. ( 16) and expressed as:

HD n (X, f ) = Γ n (X, nf ) Γ Tot (X, f ) (17a) THD(X, f ) = N n=2 Γ n (X, nf ) Γ Tot (X, f ) 2 (17b) with Γ Tot (X, f ) = N n=1 [Γ n (X, nf )] 2 (17c)
The knowledge of the Kernels in the frequency range [f 1 , f 2 ] allows for the direct computation of the THD and HD n using Eqs. [START_REF] Griesinger | Beyond MLS -occupied hall measurerement with FFT techniques[END_REF]. This can be done for any value of input amplitude X and for any frequency f in [f 1 , f 2 ]. In the present study, the prediction of the inter-modulation product (IMP) has not been included for the sake of brevity. Nevertheless, this model is fully able also to predict IMP, or to compute the output corresponding to any desired signal.

Prediction of HD n and THD at a given amplitude

Using the different sets of Kernels measured in the previous section, the HD n and THD have been predicted from Eqs. [START_REF] Griesinger | Beyond MLS -occupied hall measurerement with FFT techniques[END_REF]. To compare with, traditional measurements using pure tones have been done using the experimental protocol depicted in section 4.1. HD n and THD have been measured this way for 50 frequencies between 50 Hz and 12 kHz and for 10 different amplitudes ranging from 95 to 115 dB in pressure. In Fig. 5 the predictions made using equations [START_REF] Griesinger | Beyond MLS -occupied hall measurerement with FFT techniques[END_REF] with identified Kernels and the traditional measurements are shown. It can be seen that the agreement between measured and predicted data is very good over the entire frequency range.

Prediction of HD n and THD for different amplitudes

To have an overview of the quality of the predictions depending on the amplitude at which Kernels have been measured (X m ) and on the amplitude at which prediction are made (X p ), a mean error has been introduced. This error is defined in Eq. 18. This error has been computed for the following frequency bands: [45, 355] Hz (octave bands centered at 63 Hz, 125 Hz, and 250 Hz), [355, 2800] Hz (octave bands centered at 500 Hz, 1 kHz and 2 kHz) and [2800, 11200] Hz (octave bands centered at 4 kHz and 8 kHz). The error in dB in each of these frequency bands is shown in Fig. 6.

∆ [fA,fB ] (X m , X p ) = ... 20 log 10 mean [fA,fB ] |THD m (X p , f ) -THD p (X m , X p , f )| (18)
In Fig. 6 (a), the error in the frequency band [45,355] Hz is shown. It can be seen that this error is low, around -11 dB. The minimum of -15 dB is reached when the THD is predicted for low values of X p .

Conclusion

In this article a method for the identification of a class of nonlinear systems, referred to as cascade of Hammerstein models, has been mathematically described and applied on an electro-dynamic loudspeakers. This method permits the complete identification of the constitutive elements of the model from only one single measurement. This method is very fast compared to other methods [START_REF] Hawksford | System measurement and identification using pseudorandom filtered noise and music sequences[END_REF][START_REF] Gallman | Iterative method for identification of nonlinear-systems using a Uryson model[END_REF] where several repeated measurements are necessary. The mathematical formalism developed by the authors allows moreover an extension beyond the simple separation of the different orders of non-linearities initially proposed by Farina [START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF].

Harmonic distortion of order n (HD n ) and the total harmonic distortion (THD) are afterwards predicted using the identified Kernels. The originality of the present approach is that HD n and THD are predicted at different amplitudes. This is an advance compared to the literature where HD n and THD are usually predicted only for a given amplitude [START_REF] Kaizer | Modeling of the nonlinear response of an electrodynamic loudspeaker by a volterra series expansion[END_REF][START_REF] Farina | Simultaneous measurement of impulse response and distortion with a swept-sine technique[END_REF][START_REF] Muller | Transfer-function measurement with sweeps[END_REF]. Results obtained when comparing this approach to the traditional one at different amplitudes for HD n and THD are very satisfying.

This method is also interesting for transducer quality assessment. Indeed, the present approach gives a straightforward access to a compact and very fine representation of the linear and nonlinear characteristics of real transducers. As has been shown in this paper, total harmonic distortion and n th -order harmonic distortion are precisely predicted by the identified models. However, they badly correlate with subjective experiments involving transducer's quality [START_REF] Voishvillo | Assessment of nonlinearity in transducers and sound systems -from thd to perceptual models[END_REF]. Nevertheless, as an input-independent model is available, simulations of the nonlinear responses of identified transducers for any desired input signal (noise, music, ...) can be easily performed. This can be very useful when searching pertinent cues to assess the decreases of quality caused by non-linearities in acoustical transducers.

In a future work, models of panel-based transducers (distributed mode loudspeakers [START_REF] Harris | Introduction to distributed mode loudspeakers (DML) with first-order behavioural modelling[END_REF] or multiactuator panels [START_REF] Kuster | Structural and acoustic analysis of multiactuator panels[END_REF]) will be identified with the present method. As their physical underlying principle is different it is not obvious that cascade of Hammerstein models represent them correctly. The influence of the order of non-linearity N and the sensitivity of the method to ambient noise will also be studied.

Figure 1 :

 1 Figure 1: Block diagram representation of a cascade of Hammerstein models.

Figure 3 :

 3 Figure 3: Overview of the method used to identify a Cascade of Hammerstein models.

Figure 4 :

 4 Figure 4: Mean measured Kernel of the cascade of Hammerstein models (solid line) and the corresponding standard deviations (dashed lines).

Figure 5 :

 5 Figure 5: Comparisons between measurements (circles) and predictions (solid lines) at 106 dB for the HD 2 , HD 3 and THD. Kernels identified at 107 dB have been used for predictions.

Figure 6 :

 6 Figure 6: Mean error in the frequency band [45, 355] (a), [355, 2800] (b) and [2800, 11200] Hz (c) for the electrodynamic loudspeaker. Amplitudes of measurement of the Kernels are given on the x-axis.


Julien Chauchat

Marc Médale

A three-dimensional numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow

Keywords: Two-phase flow model, Fluid-particle coupling, Viscoplastic flow, Regularisation technique, Quadratic Finite element

A numerical model for the simulation of incompressible two-phase flows of a flat granular bed submitted to a laminar shearing flow is presented, considering a two-fluid model and a mixed-fluid one. The governing equations are discretized by a finite element method and a penalisation method is introduced to cope with the incompressibility constraint. A regularisation technique is used to deal with the visco-plastic behaviour of the granular phase. Validations are carried out on three flow test cases: a Bingham fluid between two infinite parallel planes, a Bingham fluid in a square lid-driven cavity and a Newtonian fluid over a granular bed in a two-dimensional configuration, for which we compare our numerical results with existing analytical or numerical results. The accuracy and efficiency of the numerical models have been compared for the two formulations of the two-phase flow model. It turns out that the two-fluid model requires ten times more CPU time than the mixed-fluid one for a comparable accuracy, which can be achieved provided one takes a smaller regularisation parameter in the latter model. Finally, three-dimensional computations are presented for the flow of a Newtonian fluid over a granular bed in a square and circular cross-section ducts.

Introduction

Particles transport occurs in a variety of environmental and industrial flows such as sediment transport in rivers or at coasts, hydrate formation in pipelines (oil production) or granular transport in food or pharmaceutical industries.

There are mainly two modes for particles to be transported by a flow: suspended-load or bed-load. Suspended load is the part of the load where particles are carried without contact with the bed. On the other way, bed-load is the part of the load that is carried with intermittent contact with the bed, by rolling, sliding and bouncing [START_REF] Fredsoe | Mechanics of costal sediment transport[END_REF]. In this paper, we focus only on the bed-load transport and more precisely in laminar flow conditions.

The bed-load transport is by nature a two-phase problem (fluid-particles). The particles of the bed are moving depending on the value of the shear stress exerted by the fluid, the shear stress being usually made dimensionless by the apparent weight of a single particle the so-called Shields number θ [START_REF] Shields | Application of similarity principles and turbulence research to bed-load movement[END_REF]. If the Shields number based on the fluid bed shear stress is lower than a critical value θ c the particle flux is zero otherwise it evolves as a function of θ and θ c [START_REF] Einstein | The bed load function for sedimentation in open channel channel flows[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][START_REF] Charru | Erosion and deposition of particles on a bed sheared by a viscous flow[END_REF][START_REF] Ouriemi | Determination of the critical shields number for particle erosion in laminar flow[END_REF]. The bed-load is usually modelled as a two-layer problem: a pure fluid layer at the top and a fluid-particles mixture layer at the bottom, the two layers being separated by the bed upper surface. The fluid motion in the upper layer is solved assuming weak interactions between the two layers. From this calculation the fluid bed shear stress is known and the particle flux in the lower layer is deduced from an algebraic relationship. A different approach for the bed-load has been proposed by Ouriemi et al. [START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF]. The authors have proposed some closures of the two-phase model which are appropriate to a situation in which the sediment can be considered as a mobile granular medium where the particles are in contact: the interphase force is then Darcy drag and buoyancy, the fluid phase stress is of Newtonian form, and the particle phase stress is described by a granular rheology (Coulomb friction).

The granular rheology shares some properties with visco-plastic rheology, in particular it exhibits a threshold of motion due to the friction between grains. The archetype model for visco-plastic material is the Bingham model [START_REF] Bingham | Fluidity and Plasticity[END_REF]. It is possible to identify the Coulomb friction model for the granular media with the Bingham model in which the fluid viscosity vanishes. The simulation of Bingham fluid flows have been the subject of many papers in the literature (see [START_REF] Dean | On the numerical simulation of bingham visco-plastic flow: Old and new results[END_REF] for a recent review). There are mainly two approaches to deal with the yield stress in the Bingham model: the Augmented Lagrangian technique [START_REF] Fortin | Augmented Lagrangians : Application to the Numerica Solution of Boundary Value Problems[END_REF][START_REF] Glowinsky | Augmented Lagrangians and Operator-Splitting Methods in Nonlinear Mechanics[END_REF] and the regularisation technique [START_REF] Bercovier | A finite-element method for incompressible non-newtonian flows[END_REF][START_REF] Papanastasiou | Flows of materials with yield[END_REF][START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF]. In the Augmented Lagrangian technique, the discontinuity in the Bingham constitutive relationship is treated by introducing a new primal variable and a Lagrange multiplier that enforces it to be equal to the strain rate tensor. This method is particularly accurate to capture and predict the yielded regions of the flow. In the regularisation approach, the Bingham viscosity is "regularised" by adding a small quantity to the magnitude of the rate-of-strain tensor in the denominator. The solid regime is replaced by a very viscous one. But the cost overrun for the Augmented Lagrangian Methods compared with regularization one is obvious in terms of memory due to the introduction of an additional tensor variable (the "true" strain rate), whereas in terms of CPU time the comparison is not known a priori. Actually, the regularized problem is equivalent to the flow of a shear thinning material that induces additional non-linearity in the equations.

The computation of which can significantly increase the CPU time and do not allow to conclude on the most efficient method. However, let us mention that the implementation of the regularization method is easier than the Augmented Lagrangian one. Therefore we choose the regularisation technique to deal with the yield stress in our two-phase flow model. This method is advantageous for its simplicity but one must be careful of the induced creeping flow in the yielded regions that arises when using regularisation.

In this paper we present a three-dimensional Finite Element Method (FEM) model of the two-phase incompressible flow model for bed-load transport presented by Ouriemi et al. [START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF]. Our first concern is to propose a numerical model able to predict accurately the bed-load transport in laminar shearing flows. 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 have considered two formulations of the two-phase model. In the two-fluid model the unknowns are the velocities and pressure in each phase (fluid and particles) whereas in the mixed fluid model the fluid-particles mixture is only considered assuming that fluid and particles have the same velocity. The computational efficiency of the numerical models associated with both formulations is investigated in in terms of accuracy, CPU time and memory usage. The two-phase flow model equations and the numerical modelling are presented in Section 2. Section 3 is devoted to the validation of the model by comparison with analytical solutions or published numerical results. Firstly, we have validated the numerical model for Bingham fluid flows on two test cases. We have compared the numerical model with an analytical solution for the flow of a Bingham fluid between two infinite parallel planes.

We have also compared our model with the numerical results of Mitsoulis and Zisis [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF] for the flow of a Bingham fluid in a square lid-driven cavity. Then we have validated the numerical model with the analytical solution presented by Ouriemi et al. [START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF] for the flow of a Newtonian fluid over a granular bed in a two-dimensional configuration. After the validation, we present in Section 4 the application of the two-phase model to simulate the bed-load transport in threedimensional configurations, a square cross-section and a circular cross-section ducts. Finally, we give concluding remarks in section 5.

The two-phase flow model

Following Ouriemi et al. [START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF], we present here the formulation of the two-phase flow model for bed-load transport in laminar shearing flows.

Mathematical formulation

Governing equations

Given a cartesian coordinate system (O, x, y, z) where x represents the stream-wise direction, y the lateral direction and z the vertical upward direction, the velocity vector of the k phase and its cartesian components are respectively denoted by

- → u k = (u k , v k , w k ).
k is taken to be f for the fluid phase and p for the particulate one. We start from Jackson's equations [START_REF] Jackson | The dynamics of fluidized particles[END_REF] to get the set of governing equations for the two-phase problem.

For the fluid phase, the continuity equation reads:

∂ǫ ∂t + ∇. ǫ - → u f = 0, (1)
where ǫ designates the volume fraction of the fluid phase. The particulate phase continuity equation has the same form:

∂φ ∂t + ∇. φ - → u p = 0, ( 2
)
where φ is the particulate phase volume fraction. The global volume conservation imposes φ + ǫ = 1.

The momentum equations for the fluid and particulate phases are respectively:

ρ f          ∂ǫ - → u f ∂t + ∇. ǫ - → u f ⊗ - → u f          = ∇. σ f -n - → f + ǫρ f - → g , (3)
ρ p        ∂φ - → u p ∂t + ∇. φ - → u p ⊗ - → u p        = ∇. σ p + n - → f + φρ p - → g , (4)
where σ f and σ p represent the stress tensor associated with the fluid and particulate phases respectively. n -→ f represents the average force exerted by the fluid on the particles and -→ g is the gravity acceleration vector.

The set of partial differential equations ( 1)-( 4) introduces more unknowns than the number of equations then closure relationships are needed to solve the problem. These relations are of two types: fluid-particle interactions and stress tensor expressions.

Closures

Interaction term. Following Jackson [START_REF] Jackson | The dynamics of fluidized particles[END_REF] the average force exerted by the fluid on the particles can be decomposed in two contributions. The first one corresponds to the generalized buoyancy force and the second one gathers all the remaining contributions.

n - → f = φ∇. σ f + n -→ f 1 (5)
For a viscous fluid flow in a porous media, the remaining contributions reduce to the viscous drag force due to the relative motion between phases. Using the Darcy law, the term n -→ f 1 can be written:

n -→ f 1 = η ǫ 2 K - → u f - - → u p , (6)
where η is the dynamic viscosity of the pure fluid. The coefficient of permeability is empirically linked to ǫ and the particle diameter d by the Carman-Kozeny relationship:

K = ǫ 3 d 2 k CK (1 -ǫ) 2 (7)
A typical value for k CK ≈ 180 is proposed by Happel and Brenner [START_REF] Happel | Low Reynolds number hydrodynamics[END_REF] and Goharzadeh et al. [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF].

Stress tensors. The fluid phase has been assumed to be a Newtonian viscous liquid in which the Einstein dilute viscosity formula has been chosen to be applied to the concentrated situation :

σ f = -p f I + τ f = -p f I + η e ∇ -→ u m + (∇ -→ u m ) T , (8)
where η e = η(1 + 5φ/2) is the effective viscosity of the mixture and

-→ u m = ǫ - → u f + φ - → u p
is the velocity of the mixture.

The particle phase has been assumed to be described by Coulomb solid friction in which the extra stress is proportional to the particle pressure.

σ p = -p p I + τ p (9)
In the frame of the three-dimensional model we express the Coulomb friction model in tensorial form following the idea of Jop et al. [START_REF] Jop | A constitutive law for dense granular flows[END_REF]:

τ p = η p ( γp , p p ) γp , (10)
with

η p ( γp , p p ) = µ p p γp , (11)
where the rate of strain tensor γp is defined as γp = ∇ -→ u p + (∇ -→ u p ) T and its magnitude is given by the square root of its second invariant γp = 1 2 T r γ2 . We emphasise that this expression in the one dimensional case reduces to the classical Coulomb expression: τ p xz = µ p p .

Dimensionless equations

In the preceding subsections we have presented the ingredients of the two-phase model for bed-load transport, we summarise here the model equations to be solved. We consider two formulations of the model. The first one, called the two-fluid model, is based on the solution of mass and momentum conservation equations for each phase [START_REF] Dean | On the numerical simulation of bingham visco-plastic flow: Old and new results[END_REF]. The second one, called the mixed-fluid model, is based on the solution of mass and momentum equations for the mixture (13) (i.e.: a single effective phase is considered). In this latter formulation, the mass and momentum equations are simply obtained by summing the corresponding equations over each phase (fluid and particles).

In the following, we assume that the volume fractions are constant in space and time meaning that the interface between the fluid-particle mixture and the pure fluid region is fixed and no dilatation occurs. This assumption implies that fluid and particulate phase as well as the mixture are incompressible. We also express the equations for the fluid phase in terms of the mixture velocity.

Following Ouriemi et al. [START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF] we make all the values dimensionless by scaling the length by H, the height of the flow, and the stresses by ∆ρgH, and therefore the time by η/∆ρgH where ∆ρ = ρ pρ f . Using these scales one obtains the following dimensionless equations for the two formulations of the two-phase model.

Two-fluid model.

                                                         ∇. -→ u m = 0 ∇. - → u p = 0 Ga H 3 d 3 D -→ u m Dt = -∇p f + ∇. η e η ∇ -→ u m + ∇ -→ u m T - H 2 K -→ u m - - → u p + ρ f - → g ∆ρ - → g Ga H 3 d 3 R ρ φ D - → u p Dt = -∇p p + ∇. η p η ∇ - → u p + ∇ - → u p T -φ∇p f + φ∇. η e η ∇ -→ u m + ∇ -→ u m T + (1 -φ)H 2 K -→ u m - - → u p + φ - → g - → g ( 12
)
In these equations, 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 η ∇ -→ u m + ∇ -→ u m T + ∇. η p η ∇ -→ u m + ∇ -→ u m T (13)

Numerical model

The Finite Element Method (FEM) leads to the discretisation of the variational formulation of equations ( 12) for the two-fluid model and on equation ( 13) for the mixed-fluid model.

Weak formulations

Weak formulation for the two-fluid model. Let us define the physical domains associated with the pure fluid and the fluid-particles mixture by Ω f and Ω p respectively and their respective boundaries by Γ f and Γ p (See Figure 1). From the previous system of partial differential equations ( 12) one obtains the following weak formulation of the momentum equations for the two-fluid formulation [START_REF] Carey | Finite Elements[END_REF], [START_REF] Pironneau | Méthodes des éléments finis pour les fluides[END_REF]: 14), ( 15) and find

Find -→ u m h ∈ U h and p f ∈ Q h satisfying (
- → u p h ∈ U h and p p ∈ Q h satisfying (16), (17) ∀ -→ δu ∈ V h and ∀ δp ∈ Q h where U h = { - → u ∈ H 1 (Ω) | - → u = - → u Dirichlet on ∂Ω Dirichlet }, Q h = {δp ∈ L 2 (Ω)}, V h = { -→ δu ∈ H 1 0 (Ω) | -→ δu = - → 0 on ∂Ω Dirichlet }. 0 = Ω f ∪Ω p δp.∇ -→ u m h dΩ (14) 0 = - Ω f ∪Ω p Ga H 3 d 3 -→ δu. D -→ u m h Dt dΩ - Ω f ∪Ω p -→ δu.∇p f h dΩ + Ω p H 2 K -→ δu.( - → u p h - -→ u m h ) dΩ - Ω f ∪Ω p η e η ∇ -→ δu + ∇ -→ δu T . ∇ -→ u m h + ∇ -→ u m h T dΩ + Ω f ∪Ω p ρ f ∆ρ -→ δu. - → g - → g dΩ (15) 0 = Ω p δp.∇ - → u p h dΩ (16) 6 0 = - Ω p Ga H 3 d 3 R ρ φ -→ δu. D - → u p h Dt dΩ - Ω p -→ δu.(∇p p h -φ∇p f h ) dΩ - Ω p (1 -φ)H 2 K -→ δu.( - → u p h - -→ u m h ) dΩ - Ω p η p ∇ -→ δu + ∇ -→ δu T ∇ - → u p h + ∇ - → u p h T dΩ + Γ p η p -→ δu. ∇ - → u p h + ∇ - → u p h T . - → n dΓ - Ω p φ η e η ∇ -→ δu + ∇ -→ δu T ∇ -→ u m h + ∇ -→ u m h T dΩ + Γ p φ η e η -→ δu. ∇ -→ u m h + ∇ -→ u m h T . - → n dΓ + Ω p φ -→ δu. - → g - → g dΩ ( 17
)
Weak formulation for the mixed-fluid model. For the mixed fluid model, the weak formulation is given by: Find -→ u m h ∈ U h and p f and p p ∈ Q h satisfying ( 18) and ( 19)

∀ -→ δu ∈ V h and ∀ δp ∈ Q h where U h = { - → u ∈ H 1 (Ω) | - → u = - → u Dirichlet on ∂Ω Dirichlet } Q h = {δp ∈ L 2 (Ω)} V h = { -→ δu ∈ H 1 0 (Ω) | -→ δu = - → 0 on ∂Ω Dirichlet } 0 = Ω f ∪Ω p δp.∇u m h dΩ (18) 0 = Ω f ∪Ω p Ga H 3 d 3 (1 + R ρ ) -→ δu. Du m h Dt dΩ - Ω f ∪Ω p -→ δu.(∇p p h + ∇p f h ) dΩ - Ω f ∪Ω p η p + η e η ∇ -→ δu + ∇ -→ δu T ∇ -→ u m h + ∇ -→ u m h T dΩ + Ω f ∪Ω p ρ m ∆ρ -→ δu. - → g - → g dΩ (19)

Algorithms associated with the two formulations

The first algorithm dealing with the two-fluid formulation is based on the solution of the system ( 14)-( 17) in a weakly coupled way. For each Newton-Raphson iteration we solve for the fluid momentum equation ( 15) whereas the particulate phase velocity is taken at its previous iteration value. From the value of the mixture velocity obtained at this stage the particulate phase velocity is computed by solving equation [START_REF] Glowinsky | Augmented Lagrangians and Operator-Splitting Methods in Nonlinear Mechanics[END_REF]. On the other hand the algorithm devoted to the mixed-fluid formulation is based on the solution of the system [START_REF] Goharzadeh | Transition layer thickness at a fluid-porous interface[END_REF][START_REF] Guyon | Hydrodynamique physique[END_REF].

The non-linearities in the governing equations are solved by a Newton-Raphson algorithm by taking the first variation of the variational formulations ( 15)-( 17) and [START_REF] Guyon | Hydrodynamique physique[END_REF]. Only the advective terms and the particle stress term (visco-plastic) induce a contribution in the first variation of the variational formulations.

Regularisation technique

The stress tensor for the particulate phase [START_REF] Carey | Finite Elements[END_REF] reveals the presence of a yield stress that depends on the granular pressure. The mixed-fluid model ( 13) can be identified as a Bingham fluid in which the yield stress vary linearly with the granular pressure. This remark allows us to use some well-known methods derived for Bingham flows to deal with the visco-plastic behaviour of the fluid-particle mixture.

As with a Bingham model, the particulate viscosity diverges when the shear rate tends toward zero, rising evident numerical problems. The idea of a regularisation technique from a mathematical point of view is to smooth the divergence of the viscosity function. The consequence on the computational behaviour of the yield stress fluid results in a very viscous fluid at zero rate of strain instead of a pure rigid body behaviour. One of the easiest solutions to regularise the viscosity is the following [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF]:

η(γ) = µ + τ 0 γ +λ ,
where λ is the regularisation parameter (λ << 1) and τ 0 is the yield stress. We have chosen this expression in our implementation of the numerical model. A review of regularisation technique can be found in Frigaard and Nouar [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF].

Implementation

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][START_REF] Balay | PETSc users manual[END_REF][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 [2, 1, 3] with a penalty parameter set to 10 9 for all the simulations presented in this paper.

Validations

In this section we address a twofold goal for the validation of our numerical model. The first one concerns the implementation of viscoplastic flow using a regularisation technique (subsection 3.1) and the second one deals with the implementation of the two formulations associated with the two-fluid model and the mixed-fluid one (subsection 3.2). In each case we discuss the accuracy and numerical efficiency of the implemented numerical model.

Test cases for Bingham fluid flows

We have begun with test cases on Bingham fluid flows because of the close link that exists between the Coulomb friction for the particulate phase [START_REF] Carey | Finite Elements[END_REF] and the Bingham fluid model. These two models exhibit the presence of a yield stress. This is particularly important to validate the implementation of the regularisation on "well-known" configurations (Bingham fluid flows) in order to avoid questions on it when performing two-phase flows simulations.

We have done two test cases. The first one is the flow of a Bingham fluid between two infinite parallel planes. The second one is the flow of a Bingham fluid in a square lid-driven cavity.

Flow of a Bingham fluid between two infinite plane

In this first test case we study the flow of a Bingham fluid between two infinite parallel planes. This is a 2D problem where the longitudinal velocity component only depends on the vertical coordinate. This simple problem is particularly interesting for the validation since it possesses an analytical solution for the longitudinal velocity (See Guyon et al. [START_REF] Guyon | Hydrodynamique physique[END_REF] p. 219-220). Figure 2 shows the sketch of the problem and boundary conditions. The simulation are performed using a regularisation method with a regularisation parameter λ = 10 -4 and the Bingham number (Bn) is equal to 10.

Figure 3(a) shows the longitudinal velocity profiles for the numerical solution and the analytical one. In order to study the spatial convergence of the numerical model we have considered five meshes. The ratio of the mesh size in the three directions have been kept constant ∆x/∆z = ∆x/∆y = 10 to avoid errors induced by the distortion of the mesh. The five mesh sizes are reported in Table 1 as well as the CPU time requirements. All the simulations have been performed on a cluster composed of AMD 280 processors at 2.4 GHz and 8 Go of RAM by node. The CPU time evolution with the mesh size is more or less proportional to the number of degrees of freedom. Figure 3(b) shows the Root Mean Square (RMS) error between the numerical solution and the analytical one versus the mesh size (h, defined as an element height ∆z).

The longitudinal velocity profiles obtained by numerical simulation are in close agreement with the analytical solution for the five meshes. Figure 3(b) reveals that the spatial order of convergence is approximately of order three which corresponds to the optimal convergence rate of the method (second order approximation) [START_REF] Oden | Finite Elements[END_REF][START_REF] Carey | Finite Elements[END_REF][START_REF] Pironneau | Méthodes des éléments finis pour les fluides[END_REF][START_REF] Zienkiewicz | The Finite Element Method[END_REF]. This particularly good spatial convergence rate can be explained by the fact that the analytic solution is quadratic in the plastic region and constant in the plug located in the central region of the flow. Therefore our tri-quadratic finite element approximation is able to represent at convergence this analytical solution. So, this test case reveals that a correct implementation of the model to simulate Bingham fluid flows as the optimal asymptotic order of convergence has been reached, see figure 3(b).

Bingham flow in a lid-driven cavity

The second test case concerns the flow of a Bingham fluid in a square lid-driven cavity. Its main goal is to qualitatively validate the 3D implementation of the regularisation technique. The sketch of the lid-driven cavity is presented in Figure 4. A tangential velocity is imposed at the top of the cavity whereas the three other boundaries are considered as walls. The geometry is a unit square cavity and the velocity at the top is also chosen to be unity. This test case has been extensively studied and a great amount of literature on viscoplastic flows exists concerning this problem [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF][START_REF] Lalli | Fluid-particle flow simulation by averaged continuous model[END_REF][START_REF] Dean | On the numerical simulation of bingham visco-plastic flow: Old and new results[END_REF]. We have chosen the Mitsoulis and Zisis [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF] results as reference who also used a FEM model and a regularisation technique to deal with the yield stress. In the following simulations we have neglected inertial effects (i.e.: Re =0). The dimensionless number controlling the flow is the Bingham number (Bn = τ 0 H/ηU) where τ 0 is the yield stress.

The mesh is composed with 60x1x60 quadratic elements (H27). The regularisation parameter is fixed to λ = 10 -4 s -1 . This value is equivalent to the one chosen by Mitsoulis and Zisis [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF] but the regularisation technique is different. Mitsoulis and Zisis [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF] used the Papanastasiou [START_REF] Papanastasiou | Flows of materials with yield[END_REF] regularisation whereas we have used the simple regularisation (Cf. 2.3). In the simulations the Bingham number has been considered in the range [0 -20] to validate our model.

Figure 5 shows the velocity fields in the cavity for various Bingham number (Bn= 0, 2 and 20). We observe that a rigid region appears at the bottom of the cavity as the Bingham number is increased. Also the velocity decreases as well as the vortex intensity.

Figure 6 shows the yield zone position in the lid-driven cavity for various Bingham number. We compare our results with the results from Mitsoulis and Zisis [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF]. The position of the yield zone is defined as the position where the material flows (yields) i.e. : where the magnitude of the stress tensor τ exceeds the yield stress τ 0 ( τ ≥ τ 0 ).

The results obtained are in good agreement with Mitsoulis and Zisis [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF] and validate the implementation of the Bingham model in our FEM model.

As a conclusion on the two Bingham fluid flow test cases one can notice that the regularisation technique implemented give satisfactorily accurate results for Bingham number in the range 0-20 for a regularisation parameter λ = 10 -4 . Moreover with this numerical parameter a third order (optimal) asymptotic convergence rate has been reached for our tri-quadratic finite element approximation.

Two-phase simulation of bed-load transport in 2D

In this subsection we present results on the flow of a Newtonian fluid over a granular bed. The aim of this subsection is first to validate quantitatively the two formulations of the two-phase flow model by comparison with the analytical solution of Ouriemi et al. [START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF] for the bed-load transport in laminar shearing flows and secondly to assess computational efficiency of the numerical model associated with both formulations.

The sketch of the problem and boundary conditions are given in figure 7. 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. There are several choice for the definition of the Bingham number in the bed-load problem. Actually, in dense granular media the yield stress varies with the normal stress. Therefore one has to choose a pertinent value of the yield stress. Here we choose the yield stress corresponding to the first granular layer, this choice is natural since the relevant length scale for the estimation of the yield stress is the height of the moving granular bed that is of the order of few grain diameters. Assuming a hydrostatic pressure for the granular phase, it reads:

τ 0 = µ s ∆ρgd,
where d is the particle diameter. With this yield stress definition the value of the Bingham number in the bed-load problem considered here is Bn = 2 10 4 . This value is greatly higher than the one used in the validation tests for Bingham fluid flows (Bn # 10) and implies accordingly a smaller value of λ than the one used in previous Bingham flow test cases. So 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 (See subsection 2.1.3) corresponding to the coupled system ( 15)-( 17) for the two-fluid model and equation [START_REF] Guyon | Hydrodynamique physique[END_REF] for the mixed-fluid model. Numerical solutions for the two-fluid model and the mixed-fluid model are presented in figure 8 for a 6x1x40 mesh with a requested absolute residual lower than 10 -11 per degree of freedom.

Figures 8a and8b 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. [START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF]. In both figures the black solid line represents the analytical solution. A good agreement can be noticed between the numerical results and the analytical solution which validate the FEM model for the bed-load transport in a 2D configuration.

Moreover we have performed a spatial convergence analysis for both formulations (see figure 9). Concerning the two-fluid model solution the spatial convergence is of order one. Two hypotheses could explain the first order convergence, the presence of a predominant Darcy term in the momentum equations and the regularisation technique applied to the particulate phase rheology. Moreover one can observe that the error on the particulate phase velocity is greater than the fluid one. Therefore the difference could be attributed to the regularisation technique that can be reduced using a smaller regularisation parameter λ. For the mixed-fluid model solution a third order spatial convergence is reached for the coarsest meshes where the discretisation error dominates. On the other hand for finer meshes the order of convergence reduces and tends to a first order one in the region where it is dominated by the modelling error introduced by the regularisation technique. One can also notice that the RMS error is one order of magnitude greater than for the two-fluid model solution for the same regularisation parameter as in the two-fluid model. Reducing the regularisation parameter λ by one order of magnitude for the mixed-fluid model reduces the error accordingly (see figure 9). Therefore the same accuracy as the two-fluid model is recovered but with a tenfold CPU time reduction for the mixed-fluid model compared with the former formulation (see Tab. 2). We point out that the CPU time for the mixed-fluid model with a regularisation parameter λ = 10 -7 is only twenty percent higher than the one with λ = 10 -6

whereas the error is reduced by one order of magnitude. More surprisingly, it takes the same CPU time for two uniformly refined meshes for a comparable accuracy solution (see Tab. 2). Indeed the iteration count for the finest mesh is five times less than for the coarser meanwhile the CPU time per iteration is five times higher.

Consequently the two-fluid model is much more expensive (ten times in CPU time) than the mixed-fluid one for a comparable accuracy provided one takes a regularisation parameter sufficiently small (one order of magnitude smaller than in the two-fluid model). But one should keep in mind that the mixed-fluid model is based on the strong assumption of zero relative velocity between the fluid and particulate phases, which could be restrictive in actual problems.

It turns out from the results obtained in this section that the modelling error associated with the implemented regularisation technique is connected to the regularisation parameter value λ. Therefore in order to achieve compu-tations with controlled accuracy we suggest to link the value of the regularisation parameter to the Bingham number according to the following empirical relationship: λ = min 10 -4 , 10 -4 Bn .

Two-phase simulation of bed-load transport in 3D configuration

We have shown in the last section the convergence of the two-phase numerical model compared with analytical solution for a two-dimensional configuration. We now apply the model to three-dimensional configurations: a square and a circular cross-section ducts. As in the previous test case the values of the dimensionless numbers are: Re = 2 10 -2 , Ga = 11, R ρ = 0.4, d/H = 30 and Bn = 2 10 4 . However the characteristic length in the definition of the dimensionless numbers is the side length for the square duct whereas it is the diameter in the cylindrical one.

Square section duct

We have performed simulations for a square cross-section duct with two meshes: 6x20x40 and 6x40x80 elements uniformly distributed, only the half of the domain has been solved in the transverse direction for obvious symmetry reasons. Figure 10 shows the velocity profile in a cross section of the duct. 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 twodimensional case. Due to this shear increase the velocity is lower than in the two-dimensional case. In order to test the convergence of the solution with respect to the mesh we present in figure 11 the comparison of the velocity profiles on the plane of symmetry of the duct (y=0) for the two meshes. The numerical solutions superimpose themselves except in the neighbourhood of the yield zone where the finest mesh better resolves the yield zone. Figure 12 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 for the square cross-section duct (6x20x40). This illustrates the good behaviour of the numerical model in this three-dimensional configuration.

Circular section duct

As for the square cross-section duct we have performed two simulations on the circular cross-section duct. The mesh sizes are 6x896 and 6x3596 quadratic elements in the x direction and in the cross-section respectively. Figure 13 shows the velocity profile in a cross-section of the duct (6x896). As for the square cross-section duct, the hick solid line at z = 0 represents the position of the granular bed and the contour colors represent the mixture velocity (u m ). Here again the friction is increased in this geometry compared with the two-dimensional configuration. We have compared the velocity profiles on the vertical plane of symmetry to show the convergence of the solution with respect to the mesh size (See figure 14). Again figure 15 illustrates the good behaviour of the numerical model for three-dimensional flow configurations.

Concluding remarks

In conclusion, we have developed a numerical model to simulate incompressible two-phase flow of a Newtonian fluid over a granular bed. The model is based on a penalisation method for incompressible flows and a regularisation technique for the visco-plastic behaviour of the granular phase. Validations have been carried out on three test cases:

the flow of a Bingham fluid between two infinite parallel planes, the flow of a Bingham fluid in a square lid-driven cavity and the flow of a Newtonian fluid over a granular bed. One can notice the very good agreement between the computations and the existing analytical solutions [START_REF] Guyon | Hydrodynamique physique[END_REF][START_REF] Ouriemi | Sediment dynamics. part i: Bed-load transport by shearing flows[END_REF] or numerical results [START_REF] Mitsoulis | Flow of bingham plastics in a lid-driven square cavity[END_REF]. To get these results we have considered that the modelling error associated with the implemented regularisation technique should be correlated to the regularisation parameter value λ according to the empirical relationship: λ = min 10 -4 , 10 -4 /Bn .

Concerning the two-phase flow formulation of the bed-load problem we have shown that the mixed-fluid model is computationally more efficient than the two-fluid one (roughly ten times faster and requires 20% less memory).

Moreover, in order to achieve a comparable accuracy of the two models one has to choose a regularisation parameter ten times smaller for the mixed-fluid model than for the two-fluid one. But one should recall that the mixed-fluid formulation is based on the strong assumption that the fluid-particles relative velocity is negligible which could limits its validity range to the cases of small particulate Reynolds number.

Finally we have performed three-dimensional numerical simulation of the bed-load transport in a square and a cir- 2 where N is the number of nodes in the mesh.

u = u analytic z x u = v = w = 0 H u = v = w = 0 v = w = 0 v = w = 0 ∂u ∂x = 0
= v = w = 0 z u = v = w = 0 u = v = w = 0 H v = w = 0 u = U 0 x
p ∂x = 0 v f = w f = u p = v p = w p = 0 and ∂u f ∂z = 0 u f = v f = w f = 0 v f = w f = 0 and ∂u f ∂x = 0 v p = w p = 0 and ∂u p ∂x = 0
e = 1 N N i=1 (U i -U ana i )

  cular cross-section ducts illustrating the capability of the model to deal with arbitrary geometries where no analytical solution exists. Future developments will concern the implementation of a numerical technique to simulate the motion of the fluid-granular bed interface. Our final goal is to perform three-dimensional simulations for the formation of ripples and dunes that are observed in experiments [? ].

Figure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: Sketch of the problem and boundary conditions of a Bingham fluid between two infinite parallel plane in 2D.

  Error analysis: RMS error against analytical solution e = 1N N i=1 (U i -U ana i )2 where N is the number of nodes.

Figure 3 :

 3 Figure 3: Flow of a Bingham fluid between two infinite parallel plane: longitudinal velocity profiles a) and error analysis b).

  u

Figure 4 := 20 Figure 5 :X 2 0 5 0 20 Figure 6 :

 420525206 Figure 4: Sketch of the lid-driven cavity.

  w p = 0 and ∂u

Figure 7 :Figure 8 :

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

Figure 9 :

 9 Figure 9: 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:

Figure 10 :Figure 11 :

 1011 Figure 10: Velocity profile obtained by numerical simulations with the mixed-fluid model for the square cross-section duct (6x20x40).

5 YFigure 12 :

 512 Figure12: Velocity profile obtained by numerical simulations with the two-fluid model for the square cross-section duct (6x20x40). 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 13 :

 13 Figure 13: Velocity profile obtained by numerical simulations with the mixed-fluid model for the circular cross-section duct (6x896).

Figure 14 :

 14 Figure 14: Mesh sensitivity of the velocity profiles obtained by numerical simulations with the mixed-fluid model for the square cross-section duct (6x896 and 6x3596).

5 YFigure 15 :

 515 Figure15: Velocity profile obtained by numerical simulations with the two-fluid model for 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.

Table 1 :

 1 Mesh definition and CPU time per iteration

	Tables
	Mesh definition (NX x NY x NZ) 3x1x10 6x1x20 12x1x40 24x1x80 48x1x160
	Degrees of freedom	1323	4797	18225	71001	280233
	CPU time (s/it. on 8 proc.)	0.67	1.88	6.98	27.79	116.74

Table 2 :

 2 CPU time and number of iterations

Acknowledgements

We would like to thank P. Aussillous and É. Guazzelli for fruitful discussions regarding the two-phase flow model and Y. Forterre and O. Pouliquen for discussions regarding the granular rheology. The authors also acknowledge Y.

Jobic for technical support. Funding from Agence Nationale de la Recherche (Project Dunes ANR-07-3 18-3892) is gratefully acknowledged.