Benoît Noetinger
email: benoit.noetinger@ifpen.fr

Laurène Hume
email: laurene.hume@univ-pau.fr

Robin Chatelin
email: robin.chatelin@enise.fr

Philippe Poncet
email: philippe.poncet@univ-pau.fr

Effective viscosity of a random mixture of fluids

We propose an estimation of the effective viscosity of a random mixture of Newtonian fluids that ignores capillary effects. The local viscosity of the mixture is assumed to be a random function of the position. Using perturbation expansions up to the second order, the resulting formula can be recast under the form of a simple power averaging mixing low. Numerical tests are used to assess the validity of the formula and the range of its applicability.

I. INTRODUCTION

Computing the effective viscosity of complex fluids such as suspensions and emulsions is an old problem studied by Einstein and Taylor [START_REF] Einstein | Investigations on the Theory of the Brownian Movement[END_REF][START_REF] Taylor | The viscosity of a fluid containing small drops of another fluid[END_REF] among many others. Such calculations have applications in rheology and related areas. The general issue is to find the connection between the large-scale rheological behavior, mainly the viscosity of the mixture, and information regarding the fluid microstructure. In the case of suspensions, these data can be the shapes, sizes, volume fractions, and pair correlation functions of the suspended particles [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF][START_REF] Brady | Stokesian dynamics[END_REF][START_REF] Brady | The rheological behavior of concentrated colloidal dispersions[END_REF]. Similar approaches have been used for emulsions and immiscible inclusions of one fluid into another [START_REF] Taylor | The viscosity of a fluid containing small drops of another fluid[END_REF][START_REF] Danov | On the viscosity of dilute emulsions[END_REF]. In another context, closer to thermodynamics, the viscosity of gas or liquid mixtures can be estimated from statistical mechanics principles using the properties of each fluid and microscopic interaction parameters. Many formulas have been proposed to estimate the viscosity (see Ref. [START_REF] Williams | The prediction of viscosity for mixtures using a modified square well intermolecular potential model[END_REF] and references therein). In the oil industry, glycerol-water mixtures have been studied experimentally in the context of oil refining [START_REF]A method for predicting the performance of unstable miscible displacement in heterogeneous media[END_REF][START_REF] Blackwell | Factors influencing the efficiency of miscible displacement[END_REF], leading to the so-called "quarter power mixing rule":

η eff = η -1/4 -4 = cη -1/4 Glycerol + (1 -c)η -1/4 Water -4 .
In this equation, η eff is the effective viscosity of the mixture, η Glycerol and η Water are the respective viscosities of pure fluids, and c is the glycerol concentration; • is the arithmetic averaging. To the best of our knowledge, these empirical observations do not have any theoretical underpinning.

In order to improve our understanding of viscosity homogenization, we will follow here the so-called stochastic approach [START_REF] Landau | Electrodynamics of Continuous Media[END_REF][START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF][START_REF] Dagan | Flow and Transport in Porous Formations[END_REF][START_REF] Noetinger | Use of the Fourier-Laplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs[END_REF][START_REF] Stepanyants | Effective hydraulic conductivity of a randomly heterogeneous porous medium[END_REF][START_REF] Willot | Elastic behavior of composites containing boolean random sets of inhomogeneities[END_REF], which is popular in the fields of random and composite materials and in hydrogeology. The key idea is to consider that the mixture can be represented by a single fluid characterized by a local viscosity modeled by a stationary random function which depends on the position. This random function is described by its two first moments: its local average and its two-point correlation function. Such approaches are still the subject of active investigations in the study of random materials when determining both the effective Young modulus and the electrical conductivity [START_REF] Noetinger | Use of the Fourier-Laplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs[END_REF][START_REF] Stepanyants | Effective hydraulic conductivity of a randomly heterogeneous porous medium[END_REF][START_REF] Willot | Elastic behavior of composites containing boolean random sets of inhomogeneities[END_REF][START_REF] Jikov | Homogenization of Differential Operators and Integral Functionals[END_REF] and references therein. Related mathematical techniques range from homogenization and stochastic perturbation techniques to diagrammatic techniques of theoretical physics. Due to the formal analogies in the underlying equations, such approaches were also followed to study the permeability of heterogeneous rocks [START_REF] Landau | Electrodynamics of Continuous Media[END_REF][START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF][START_REF] Dagan | Flow and Transport in Porous Formations[END_REF][START_REF] Stepanyants | Effective hydraulic conductivity of a randomly heterogeneous porous medium[END_REF][START_REF] Noetinger | The effective permeability of a heterogeneous porous medium[END_REF][START_REF] Wit | Correlation structure dependence of the effective permeability of heterogeneous porous media[END_REF][START_REF] Renard | Calculating equivalent permeability: A review[END_REF][START_REF] Quintard | Transport in chemically and mechanically heterogeneous porous media. III. Large-scale mechanical equilibrium and the regional form of Darcy's law[END_REF]. As far as we know, stochastic homogenization has never been carried out for viscosity. This is probably due to the emphasis on the study of suspensions and emulsions which require specific techniques, because most of the analytical difficulty comes from the boundary condition at the interface separating both fluids. In these systems, there is a strong focus on the relationship between the microstructure of the suspension and its time evolution which can lead to quite subtle organizational effects [START_REF] Petit | Shear-induced structures in macroscopic dispersions[END_REF][START_REF] Scirocco | Effect of the viscoelasticity of the suspending fluid on structure formation in suspensions[END_REF][START_REF] Malkin | Self-organization in the flow of complex fluids (colloid and polymer systems): Part 1: Experimental evidence[END_REF].

It seems rather difficult to study and set up real systems with random viscosity. First, emulsions with low surface tension seem to be good prospects. Indeed, the low surface tension hypothesis ensures that the pressure jump at the bubbles' interfaces which maintains their sphericity does not change the analysis. Otherwise a Taylor-like calculation accounting for the details of the flow inside and outside a single bubble [START_REF] Taylor | The viscosity of a fluid containing small drops of another fluid[END_REF][START_REF] Danov | On the viscosity of dilute emulsions[END_REF] is better suited. Moreover, a variable viscosity Stokes model is typical of biomicrofluidics with heterogeneities. A second example is shear-thinning fluids exhibiting space-variable concentrations [START_REF] Chatelin | Numerical and experimental investigation of mucociliary clearance breakdown in cystic fibrosis[END_REF] C in a globally constant shear fluid for which the Stokes equations are a good model as a first approximation (such as mucus or blood plasma). This leads to a viscosity μ(C,D) whose fluctuation correlation length remains statistically isotropic over time (D denotes the strain rate).

A third candidate could be a system with random temperature fields inducing in turn viscosity variations appearing as random. In cases in which the thermal diffusion coefficient is smaller than the dynamic viscosity associated with the momentum diffusion, the predictions of present work could be tested experimentally. Finally, in the widely studied case of suspensions, at a given scale the fluctuations of averaged volume fraction of the suspended particles also induces fluctuations of the local effective viscosity of the suspension. These fluctuations have an overall effect that can be studied following our approach. So in turn, these fluctuations can be related to the overall effective viscosity of the mixture. The present work could furthermore contribute to our understanding of fluctuation effects in suspensions, which are known to be of great importance (see Ref. [START_REF] Mucha | A model for velocity fluctuations in sedimentation[END_REF] and references therein). Another goal is to justify the emergence of the quarter power mixing rule. In short the present study provides a quantification of any fluctuating effect on viscosity as long as the one-phase Newtonian model is acceptable. In addition, our result can provide an analytic benchmark to test numerical Stokes solvers, as is done in random porous media [START_REF] Romeu | Calculation of internodal transmissivities in finite difference models of flow in heterogeneous porous media[END_REF].

The paper is organized as follows: in Sec. II we present the Stokes equation and notations. In Sec. III we describe the random viscosity model. The definition of the effective viscosity is given in Sec. IV. In Sec. V A the perturbation method is carried out by means of a fluctuation expansion of the Stokeslet, as a power series of the viscosity fluctuations. Averaged results up to the second order are presented in Sec. V B. In Sec. VI we present the numerical methodology and results, followed by concluding remarks in Sec. VII.

II. BASIC HYPOTHESIS AND EQUATIONS

A. Basic equations and notations

We consider the flow of an incompressible Newtonian fluid in an infinite three-dimensional (3D) domain in the context of low Reynolds number hydrodynamics. The fluid's motion is governed by Stokes equations describing momentum and mass conservation:

∇ • τ (r) + f(r) -∇p(r) = 0, τ (r) = η(r){∇v(r) + [∇v(r)] t }, (1)
∇ • v(r) = 0, 004100-2
where r, v(r), η(r), and τ (r) denote spatial position vector, the local fluid velocity, the local viscosity, and the local viscous stress tensor, respectively. The quantities f(r) and p(r) denote an arbitrary body force field and the associated pressure, respectively. The boundary conditions at infinity are implicitly assumed to be v(r) = 0, p(r) = 0, if the force field f(r) decays sufficiently fast.

B. The homogeneous case

In this section we consider that the viscosity is uniform: η(r) = η 0 . Due to the underlying linearity of the Stokes problem we can write

v(r) = d 3 r G(r -r ) • f(r ), (2)
p(r) = d 3 r P(r -r ) • f(r ). (3)
The second order tensor G(r) is known as the Oseen tensor and is given by G(r) = 1 8πη 0 ( 1 |r| + rr |r| 3 ).

Standard bold notations are used for vectors. The first order tensor P(r) = r 4π|r| 3 is the pressure tensor.

The convolution form suggests using Fourier transforms defined by

h(q) = d 3 re iq•r h(r). ( 4
)
This yields simpler linear relations between the Fourier transforms:

v(q) = G(q) • f(q), (5)
p(q) = P(q) • f(q), (6)
where

G(q) = (1 -q q) η 0 q 2 , ( 7
)
P(q) = iq q 2 . ( 8
)
Here q = q q is the unit vector built with q. Using components we get more explicit expressions:

G αβ (q) = δ αβ -qα qβ η 0 q 2 , ( 9
)
P α (q) = iq α q 2 . ( 10
)
When there is no ambiguity r and q will denote the real and Fourier variables, respectively.

III. THE RANDOM MIXTURE

In this section, we use the stochastic approach to introduce the random mixture model. The idea is to consider the local viscosity as being a random variable given by η(r) = η 0 + δη(r), [START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF] where η 0 is the arithmetic mean viscosity. The quantity δη(r) is assumed to be a stationary random variable, depending on the location r, with the following properties:

η(r) = 0, ( 12
)
δη(r)δη(r ) = C(r -r). (13)
The brackets • represent the averaging over all possible realizations. The pair correlation function C(rr) of the viscosity fluctuation describes the spatial correlation between fluctuations at two points. For rr = 0, C(0) is the local variance of the viscosity. In order to have a well-defined Fourier transform C(q) we assume that C(rr) decays sufficiently fast at infinity. In the isotropic case, the correlation function depends only on the modulus of |r -r|, so we have C(rr) = C(|r -r|). We also have C(q) = C(q) for the associated Fourier transform. In the simplest model, δη(r) is assumed to be a multi-Gaussian variable. There are standard and efficient algorithms [START_REF] Le Ravalec | The FFT moving average (FFT-MA) generator: An efficient numerical method for generating and conditioning Gaussian simulations[END_REF] to generate 3D maps of δη(r) for an arbitrary covariance function C(r) of the spatial range l c . In the sequel we will consider an isotropic two-point correlation function depending only on the modulus of the lag vector (rr).

If we are not careful the viscosity can become negative in Eq. ( 11) because the noise δη(r) has a multi-Gaussian distribution. In order to avoid this problem, we consider that the logarithm log η of the local viscosity follows a normal distribution, i.e., η follows a log-normal distribution, which ensures positiveness. Therefore η(r) = η g exp[σ Z(r)] where the quantities η g and σ 2 are the geometric average of the viscosity and the variance of the logarithmic viscosity, respectively. The quantity Z(r) is a random function of position assumed to be multi-Gaussian and having a unit variance. Moreover σ 2 is dimensionless. This gives a nonambiguous mathematical meaning to a small variance.

IV. DEFINITION OF THE EFFECTIVE VISCOSITY

In this section, we define the effective viscosity η eff of the random mixture. We follow a technique that was proposed in the study of suspensions [START_REF] Beenakker | The effective viscosity of a concentrated suspension of spheres (and its relation to diffusion)[END_REF][START_REF] Noetinger | A two fluid model for sedimentation phenomena[END_REF] or random porous media [START_REF] Noetinger | The effective permeability of a heterogeneous porous medium[END_REF].

Let f(r) be an arbitrary body force field acting on the fluid. In the case of a homogeneous fluid of viscosity η eff , one gets, using Fourier transforms:

v(q) = (1 -q q) η eff q 2 • f(q), (14)
p(q) = iq q 2 • f(q). (15)
In the case of random mixture we expect after averaging a linear response of the form

v(q) = (1 -q q) η eff (q)q 2 • f(q), (16)
p(q) = p(q) iq q 2 • f(q).
Such a form can be expected due to the linearity of the Stokes equations and to the statistical translational invariance of the system that leads to a convolution form after averaging. The second order calculation, presented in Sec. V, confirms this point. The so-called effective or equivalent fluid viscosity corresponds to the limit of these equations when q tends to 0. An algebraic definition of η eff can thus be proposed, if the limit exists:

η eff = lim q→0 η eff (q). ( 17
)
By analogy with random porous media [START_REF] Noetinger | The effective permeability of a heterogeneous porous medium[END_REF] and using the isotropy of the system, we can rigorously predict the equality p(q) . = 1. This can be justified using arbitrary force fields derived from a potential: we consider a force field f(r) = -∇φ(r), where the force potential φ(r) is arbitrary with Fourier transform φ(q). In this case one obtains v(q) = 0 and p(q) = φ(q). As this equality is valid for any φ(q), we deduce that p(q) = 1. This means that the average pressure tensor remains equal to its homogeneous counterpart after averaging. In other words, the role of pressure is to project the equation on divergence-free fields. This can be better understood by writing Eqs. ( 16) in the physical 004100-4 space:

-

d 3 r η eff (|r -r|)∇ 2 v(r ) = f(r) -∇ p(r) , ∇ • v(r) = 0. ( 18
)
This equation once again means that p(q) = 1. Using similar methods [START_REF] Noetinger | Use of the Fourier-Laplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs[END_REF], we can show that the viscosity kernel η eff (r) has a spatial range equal to the underlying correlation length l c of the input random viscosity. Given the low frequency limit, the convolution in Eq. ( 18) can be approximated by a local product

d 3 r η eff (|(r -r)|)∇ 2 v(r ) ∼ d 3 r η eff (|(r )|)∇ 2 v(r)
. The large-scale effective viscosity can be defined as η eff = d 3 r η eff (|(r )|) = η eff (q = 0). Alternatively, using formal perturbation expansion techniques like in Ref. [START_REF] Noetinger | Use of the Fourier-Laplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs[END_REF], η eff can also be derived using Feynman graphs and Dyson's equation, which involves many resummations of higher order terms.

This provides directly the form of the effective equation driving the average velocity v(r) , like

Eq. ( 18), driving the mean velocity rather than directly the mean velocity under the form [START_REF] Jikov | Homogenization of Differential Operators and Integral Functionals[END_REF]. In other words, with this resummation we can obtain directly an expansion of η eff as a power series of the viscosity fluctuation instead of 1/η eff .

In the next section, we present the basis of this perturbation expansion, without presenting the complete resummation techniques, which are not useful at the second order.

V. SETTING UP A PERTURBATION EXPANSION

A. General approach

In this section, we set up the perturbation expansion by writing formally η(r) = η 0 + εδη(r) and expanding the solution of Eqs. (1) as a formal power series of ε, called a Neumann series expansion.

Here ε plays the role of a scaling parameter of the viscosity mean-square deviation. Strictly speaking, in order to avoid convergence issues, it would be more accurate to work directly with the logarithm of the viscosity. However, in the present paper the calculation of the Stokeslet solution will be restricted to second order terms of ε. Hence it is more convenient and it simplifies the presentation to work with this representation of the mixture. The log-normal transformation will be carried out only at the end of the section, keeping in mind that the log viscosity variance σ 2 is assumed to be small. Analogous situations and choices of variables arise when evaluating the effective conductivity of composites or of random conductivity materials [START_REF] Landau | Electrodynamics of Continuous Media[END_REF][START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF][START_REF] Dagan | Flow and Transport in Porous Formations[END_REF][START_REF] Noetinger | The effective permeability of a heterogeneous porous medium[END_REF][START_REF] Wit | Correlation structure dependence of the effective permeability of heterogeneous porous media[END_REF][START_REF] Neuman | Prediction of steady state flow in nonuniform geologic media by conditional moments: Exact nonlocal formalism, effective conductivities, and weak approximation[END_REF].

We next introduce the sequence of velocities v 0 (r), v 1 (r), . . . , v n (r), . . . , and p 0 (r), . . . , p n (r), . . . defined recursively by the solution to the following set of Stokes equations:

-η 0 ∇ 2 v 0 (r) = f(r) -∇p 0 (r) = 0, (19)
∇ • v 0 (r) = 0. ( 20
)
We recognize the unperturbed Stokes equations. The recursion may be written as

-η 0 ∇ 2 v n+1 (r) -∇ • {εδη(r)[∇v n (r) + ∇v n (r) t ]} = -∇p n+1 (r), (21)
∇ • v n+1 (r) = 0. ( 22
)
Ignoring convergence issues we can check that v(r) = ∞ n=0 v n (r) and p(r) = ∞ n=0 p n (r) are solutions of Eqs. [START_REF] Einstein | Investigations on the Theory of the Brownian Movement[END_REF]. It can be shown by direct recursion using Eq. ( 21) that each term v n (r) of this expansion is of order ε n . Using the formal solution of the homogeneous case Eq. ( 2), one obtains

v n+1 (r) = -ε d 3 r G(r -r )∇ • {δη(r )[∇v n (r ) + ∇v n (r ) t ]}, (23)
which may be rewritten using the tensorial components:

v n+1 α (r) = -ε d 3 r G αβ (r -r )∂ ν 1 δη(r ) ∂ ν 1 v n β (r ) + ∂ β v n ν 1 (r ) . ( 24
)
004100-5

Here ∂ x, y, or z indicate spatial derivations with respect to the space variable x,y,z as well as the variable if there is some ambiguity. The Einstein summation convention over repeated indices is fully satisfied. Finally, using integration by parts, one obtains

v n+1 α (r) = -ε d 3 r ∂ ν 1 r G αβ (r -r ) δη(r ) ∂ ν 1 v n β (r ) + ∂ β v n ν 1 (r ) . ( 25
)
This formula shows that the nth order term involves rather complex integrations of products of n viscosity fluctuations evaluated at n different points. After averaging, the first order term will vanish.

Combined with cumulant expansions of the nth order correlation functions of the viscosity, this formula is the starting point of diagrammatic expansions that are beyond the scope of this paper. The rest of the paper is limited to second order term (n = 2).

B. Second order results

In this section, we compute the second order correction obtained after two applications of the recursion equation Eq. ( 25):

v 2 α (r) = ε 2 d 3 r 1 d 3 r 2 ∂ ν 1 r 1 G αβ (r -r 1 ) δη(r 1 ) ∂ ν 1 r 1 ∂ ν 2 r 2 G βγ (r 1 -r 2 ) δη(r 2 ) × ∂ ν 2 v 0 γ (r 2 ) + ∂ γ v 0 ν 2 (r 2 ) + ∂ βr 1 ∂ ν 2 r 2 G ν 1 γ (r 1 -r 2 ) δη(r 2 ) ∂ ν 2 v 0 γ (r 2 ) + ∂ γ v 0 ν 2 (r 2 ) . ( 26
)
Since C(r 2r 1 ) = δη(r 1 )δη(r 2 ) , one obtains after averaging over the viscosity fluctuations:

v 2 α (r) = ε 2 d 3 r 1 d 3 r 2 ∂ ν 1 r 1 G αβ (r -r 1 ) ∂ ν 1 r 1 ∂ ν 2 r 2 G βγ (r 1 -r 2 ) C(r 2 -r 1 ) × ∂ ν 2 v 0 γ (r 2 ) + ∂ γ v 0 ν 2 (r 2 ) + ∂ βr 1 ∂ ν 2 r 2 G ν 1 γ (r 1 -r 2 ) C(r 2 -r 1 ) ∂ ν 2 v 0 γ (r 2 ) + ∂ γ v 0 ν 2 (r 2 ) . ( 27
)
Up to the second order, this equation describes the average flow modification due to pair correlations between two viscosity fluctuations. Due to the statistical homogeneity, the final result appears under the form of several convolution products which are written under the form of simple products in the Fourier domain:

v 2 α (q) = ε 2 iq ν 1 (1 -q q) αβ η 0 q 2 H ν 1 ν 2 βγ (q) iq ν 2 (1 -q q) γ δ η 0 q 2 + iq γ (1 -q q) ν 2 δ η 0 q 2 f δ (q) + H βν 2 ν 1 γ (q) iq ν 2 (1 -q q) γ δ η 0 q 2 + iq γ (1 -q q) ν 2 δ η 0 q 2 f δ (q) . ( 28
)
In this equation the tensor H ij kl (r) is defined by

H ij kl (r) = ∂ i ∂ j G kl (r)C(r), (29)
and its Fourier transform is given by

H ij kl (q) = d 3 re iq•r H ij kl (r).
Using our definition of the effective viscosity Eq. ( 17), we are interested in the low q (q → 0) behavior of Eq. ( 28) that is expected to be proportional to 1/q 2 . This can be verified by direct inspection of factors involving iq (1-q q) η 0 q 2 inside Eq. ( 28), assuming that H ij kl (q = 0) has a finite value. Assuming that the correlation function decays sufficiently fast at infinity [meaning lim q→0

H ij kl (q) = H ij kl (q = 0) 004100-6
exists], we check that the low-q behavior of Eq. ( 26) is the same as

v 2 α (q) = ε 2 iq ν 1 (1 -q q) αβ η 0 q 2 H ν 1 ν 2 βγ (q = 0) iq ν 2 (1 -q q) γ δ η 0 q 2 + iq γ (1 -q q) ν 2 δ η 0 q 2 f δ (q) + H βν 2 ν 1 γ (q = 0) iq ν 2 (1 -q q) γ δ η 0 q 2 + iq γ (1 -q q) ν 2 δ η 0 q 2 f δ (q) . ( 30
)
In order to simplify notations, we can now set ε = 1.

In the last equation, the fourth rank tensor H(q = 0) is defined by

H ij kl (q = 0) = d 3 r∂ i ∂ j G kl (r)C(r). ( 31
)
This equality can be transformed using Parseval's identity:

H ij kl (q = 0) = 1 (2π ) 3 d 3 qiq i iq j (1 -q q) kl η 0 q 2 C(q) = 1 (2π ) 3 d 3 qi qi i qj (1 -q q) kl η 0 C(q). ( 32
)
We show in the Appendix that in the isotropic case, the following equality holds:

H ij kl (q = 0) = - 4 15 δ ij δ kl + 1 15 {δ ik δ jl + δ il δ jk } C(r = 0) η 0 . ( 33
)
The tensor contractions involved when inserting Eq. ( 32) in Eq. ( 30) are carried out in the Appendix.

The final second order correction to the Stokeslet is equal to 2 5

(1

-q q) η 2 0 q 2 × C(r = 0) η 0 . ( 34
)
As expected, this is proportional to the Oseen tensor 1-q q q 2 . This result is remarkable because it shows that up to the second order, in the Fourier domain, the q → 0 limit of the second order correction of η eff depends only on the local variance of the viscosity fluctuations, and not on the particulars of the whole correlation function. With this form we can use the equivalent viscosity η eff , Eq. ( 16) and Eq. ( 17) to get by direct identification the formula

η eff = η 0 1 - 2 5 C(0) η 2 0 , (35)
which is valid up to the second order. It is now more convenient to return to logarithms using the formula [START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF], which is exact for log-normal variables and which is derived in the Appendix:

η ω 1 ω = η g (e (ωσ ) 2 /2 ) 1 ω = η g exp ωσ 2 2 . ( 36
)
In the general case, this formula relating first and second moments of a random variable and its logarithm is still valid up to the second order of interest. Using the same formula, η 0 = η = η g exp σ 2 2 . Thus one obtains at the second order:

η eff = η 0 1 - 2 5 C(0) η 2 0 = η g 1 + 1 2 C(0) η 2 0 - 2 5
C(0) η 2 0 = η g 1 + 1 10
C(0) η 2 0 . 004100-7
The result in Eq. ( 35) may then be recast under the more compact and robust form of a power averaging formula, sharing the same second order expansion with respect to the log-viscosity variance σ 2

C(0) η 0 2
for small σ 2 :

η eff = η ω 1 ω (37)
with ω = 1 5 .

We recall that strictly speaking, this result is valid only at the second order. In practice it can be computed even for a large viscosity variance. Its theoretical validity domain is controlled by the validity of the second order expansion, and by the robustness of the proposed compact expression.

This power law expression (37) corresponds to assumptions regarding higher order terms that are convenient to give the compact form (37). Presently, the validity of such assumptions may be evaluated only by means of laboratory or numerical experiments, or by a more complex theory that is still missing. Some clues of its existence can be guessed in the case of random material conductivity, in which these power law averages are known to be exact in one or two dimensions [START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF]. In the general case, fourth order calculations or renormalization group arguments help in such procedures [START_REF] Stepanyants | Effective hydraulic conductivity of a randomly heterogeneous porous medium[END_REF][START_REF] Noetinger | The effective permeability of a heterogeneous porous medium[END_REF][START_REF] Indelman | A higher-order approximation to effective conductivity in media of anisotropic random structure[END_REF][START_REF] Noetinger | Computing the effective permeability of log-normal permeability fields using renormalization methods[END_REF]. Direct calculations done at the sixth order (corresponding to ε 6 carried out independently by Refs. [START_REF] Wit | Correlation structure dependence of the effective permeability of heterogeneous porous media[END_REF] and [START_REF] Abramovich | Effective permittivity of log-normal isotropic random media[END_REF]) show that the assumptions leading to such resummations are incorrect. Moreover the corresponding high-order terms depend explicitly on the particulars of the correlation function. We expect these approximations to be valid if the correlation length is smaller than the characteristic length of the samples [START_REF] Noetinger | The effective permeability of a heterogeneous porous medium[END_REF]. Even if in practice such formulas work very well, their theoretical foundation remains unknown. A mathematical study of the convergence of perturbed Neumann expansions would be useful for the community of disordered materials.

VI. NUMERICAL TESTS

In this section we present numerical evidence of the validity of Eq. (37) using a Stokes solver for viscous flows. This validation requires the generation of random viscosity maps and the solution to the Stokes equations with variable viscosity. These tests are difficult to carry out because our homogenization result is valid only for a small viscosity variance, and simulations in these conditions may yield noisy results. A compromise must then be found between an accurate signal-to-noise ratio (which means working with quite large σ 2 ) and a sufficiently small value of σ (in which case we remain in the theoretical validity domain).

A. Numerical method and implementation

We now set up a work flow to compare the effective viscosity derived from the computed velocity field and the value given by Eq. ( 37). We consider a cubic domain filled by a random mixture of viscosity η. Inside the domain, fluid flows periodically between the lateral faces {x = x min },{x = x max },{y = y min }, and {y = y max } under a uniform force in the x direction. On the horizontal faces {z = z min } and {z = z max }, no-slip conditions are imposed for the velocity field v:

v = 0.
The velocity field v satisfies the Stokes equations:

-div[2ηD(v)] = f -∇p in , ( 38
)
where D(v) = (∇v + ∇v T )/2 is the strain rate tensor, p is the pressure, and f is the external force.

We also set the incompressibility condition ∇ • v = 0 in and periodic boundary conditions on lateral faces. We note that this flow corresponds to the Poiseuille flow in the constant viscosity case. is the projection on divergence-free fields, with no slip-through condition (i.e., no normal velocity) 004100-8 at boundaries. It has been shown [START_REF] Chatelin | A hybrid grid-particle method for moving bodies in a 3D Stokes flow with variable viscosity[END_REF][START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF] that introducing the sequence {v * k } satisfying

-η∇ 2 v * k+1 = f + [2D(Pv * k ) + (∇ • v * k )I d]∇η in (39)
leads to the solution v of Eq. ( 38) by means of v = lim k→∞ Pv * k .

Boundary conditions for Eq. ( 39) at each iteration are v * k+1 = ∇ζ (v * k ) on {z = z min } and {z = z max }, and periodic conditions for other faces. Moreover, solving Eq. ( 39), a Poisson equation, requires only the use of straightforward finite-difference stencils and a fast Fourier transform (FFT) solver. This leads to a robust numerical method, well suited for variable viscosity flows even with fluid-structure interactions [START_REF] Chatelin | Numerical and experimental investigation of mucociliary clearance breakdown in cystic fibrosis[END_REF][START_REF] Chatelin | Hybrid grid-particle methods and penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT[END_REF][START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF].

Once this calculation is carried out and the velocity field is computed, the effective viscosity η eff of the mixture can be identified using the equivalent flow rate obtained solving the Poiseuille flow with constant viscosity η eff . If we solve analytically the Stokes equation with homogeneous viscosity:

η eff ∇ 2 v eff = f in , ( 40
)
with the same boundary conditions as in Eq. ( 38), we can then estimate η eff by direct identification.

The effective viscosity η eff is defined such that the velocity field v eff solution to Eq. ( 40) has the same mean flow rate as the field v solution to Eq. ( 38). Given that the external force f has only an

x component, then using the boundary conditions, we can write v eff = (v(z),0,0) T . The analytical solution of the Stokes equation may be found in textbooks and yields

η eff = (y max -y min )z 3 max 12Q num f, (41)
where the global flow rate Q num is computed from the numerical integration of the solution v. The resulting effective viscosity η eff can be compared to the result given by Eq. (37).

In our computations, the values x min ,y min ,z min are set to 0 and x max ,y max ,z max are set to 1, so that the domain is the unit cube. This cell is then discretized with a Cartesian grid having resolution N 3 (N points in each direction, N = 129,193, or 257).

B. Random field generation

Random viscosity fields are built using a 3D version of FFT-MA correlated random field generator [START_REF] Le Ravalec | The FFT moving average (FFT-MA) generator: An efficient numerical method for generating and conditioning Gaussian simulations[END_REF]. This algorithm produces random fields efficiently by taking advantage of the Fourier transform. In the algorithm's implementation, we first generate a symmetric normalized covariance field C N , following the so-called Gaussian model:

C N (h) = exp - h l c 2 , ( 42
)
where h is the space lag and l c is the correlation length. In what follows the length l c is always expressed with respect to the space step δx of a given discretization grid for the unit cube. This characteristic length is then written in the form n c δx, where n c is an integer. The associated variance is taken equal to unity. A Gaussian white noise z is also generated, and the final random field can be written as

Z(r) = Z 0 + F -1 [ F(C N )F(z)(r)], (43)
where F (respectively F -1 ) is the discrete Fourier transform (respectively the discrete inverse Fourier transform) in the 3D space, and Z 0 is the arithmetic average of Z.

The final viscosity field is then given by

η = η g exp{σ F -1 [ F(C N )F(z)]}. (44)
In the sequel, the quantity η g denotes the expected geometric average of η and will be set to 1. In that the histograms of viscosity logarithms are Gaussian. The small discrepancies may be explained by finite-size effects.

With this particular Fourier-based method, we can build periodic fields in the x and y directions in order to ensure consistency with boundary conditions of the reference Poiseuille problem and in order to avoid spurious discontinuities. In the z direction, the presence of the nonslip conditions does not impose any z periodicity. After the observation of results with periodic or nonslip conditions in this direction, the direct effect on the estimated effective viscosity is negligible. We recall that the different correlation lengths are expressed under the form n c δx, where n c varies from 1 to m, m being the greater integer such that m × dx 0.04. We also recall that the input viscosity geometric average is set to η g = 1.

Figure 3 shows expected and computed properties for 900 independent viscosity random fields generated by FFT-MA with resolution 257 3 . On the top plot, each point corresponds to a posteriori estimation of the geometric average using volume average of the corresponding viscosity map. We note that the computed average is then close to the input geometric mean (equal to 1). The dispersion 004100-10 of the data depends strongly on the correlation length. We can expect such behavior due to the central limit theorem: the variance of the volume average scales as σ 2 (l c 3 /V ), l c V 1/3 where V is the volume of the unit cell (V = 1 in present case). We also note the large scaling factor of the plot in Fig. 3. In order to conclude this section, the bottom graph in Fig. 3 presents the empirical log viscosity variance computed for each of the 900 fields versus the input imposed variance σ . We observe a perfect agreement as the correlation between observed and input variance is excellent.

C. Interpretation methodology

In the sequel we present the simulation results on the unit cell discretized with 129,193, and 257 grid points in each direction for both the log viscosity and the numerical solution of the Stokes equation. This discretization choice corresponds to space steps δx = δy = δz going from 7.8 × 10 -3 to 3.9 × 10 -3 .

In order to set up our comparisons, we still take advantage of the log-normality to rewrite the power average formula (37) with ω = 1/5 in the form ln(η eff ) = 1 5

σ 2 2 . ( 45
)
This result shows that by plotting ln(η eff ) with respect to σ 2 /2, a slope of 1/5 is expected at least for small variances. We recall that our theoretical result (37) was obtained using averaging of the apparent viscosity over all the detailed viscosity realizations. In order to test this theoretical result, we have to compute the effective viscosity of several independent realizations of the input viscosity maps, and then average these values. The associated linear regression should then exhibit the 0.2 slope. For each realization, the effective viscosity η eff is computed with the expression (41).

We set up a Monte Carlo study by computing the effective viscosity of several independent viscosity realizations sharing the same covariance structure. Each experiment is conducted as follows:

first, we generate a random viscosity field η with FFT-MA (44) with unit geometric mean and expected variance σ 2 ; then we solve Stokes equations with variable viscosity field η; finally we compute ln(η eff )

for this input viscosity field using Eq. (41). A linear regression is then performed using the resulting set of values ln(η eff ) i with respect to σ 2 i , where i denotes the ith experiment with the associated ith viscosity random field. For each grid size and correlation length, 10 sets of viscosity random fields 004100-11 For a given variance or correlation length, we generate additional sets and go on until the regression slope converges, with a maximum of 40 sets in total. We stop adding sets when the resulting slope is close to the initial one. This method leads to a massive data storage and the manipulation of many data files. For example the total size of all the viscosity fields generated for 257 3 resolution represents more than 800 GB.

D. Results

On Fig. This may be explained by a systematic bias introduced by the poor discretization. However, this observation is in good agreement with the results obtained by Ref. [START_REF] Romeu | Calculation of internodal transmissivities in finite difference models of flow in heterogeneous porous media[END_REF] in the case of conductivity and permeability. In the opposite case of large correlation lengths, a regression slope close to unity is expected. For the very large correlation lengths (compared to the overall system size), the viscosity map is in fact homogeneous. The viscosity η eff obtained by simulation will then be equal to the uniform common value, which remains a random log-normal variable (one single value per map).

For other in-between correlation lengths, the apparent averaging exponent is expected to increase as the correlation length increases, with a plateau at 0.2. Figures 4 and5 show the regression plots, which highlight the averaging exponent for different l c at resolution 257 3 .

In Fig. 6 the averaging exponents given by through the regressions are plotted with respect to the correlation length for three grid resolutions. The curve starts from small negative values for uncorrelated media, and then reaches a plateau close to the theoretical value of 1/5 before increasing to unity. We note that a very high resolution (257 3 ) is needed to validate the theoretical result.

We can also observe that, for each regression, the data deviation from the regression line is of the form ε (real σ ) 2 /2 (see, e.g., Fig.

VII. COMMENTS AND CONCLUDING REMARKS

Simulations provide results in good agreement with the power averaging formula (37), as the regression slopes present a plateau close to the theoretical value of 0.2. For the two coarser grids of 129 3 and 193 3 resolutions, an inflection in the slopes can be observed around this value. In order to obtain significant results, both a sufficient number of discretization points per correlation length and many correlation lengths are necessary to ensure the stabilization of the effective viscosity. There is actually no generalized result on the convergence with respect to the grid resolution, although we can expect such converging behavior because of the central limit theorem.

The range of σ 2 for viscosity random fields has been carefully chosen. Small values of σ 2 under 0.001 lead to meaningless results because small numerical errors dominate when viscosity fluctuations are so small. In the same way, values greater than 0.01 must be considered with caution 004100-15

For a given vector

  field u, we will denote by ζ (u) the solution of the Poisson equation -∇ 2 ζ = -∇ • u in with homogeneous Neumann boundary conditions, so that Pu = u -∇ζ

Fig. 1

 1 Fig.1three random fields are plotted with different correlation lengths (l c = 3δx, 10δx, and 25δx at resolution 2573 ). Plotting the log viscosity histograms of one particular realization in Fig.2 shows

FIG. 2 .

 2 FIG. 2. Histograms of the values of ln(η) (σ 2 = 0.3) for three correlation lengths.

FIG. 3 .

 3 FIG. 3. Empirical and expected properties of generated η fields.are generated where the ith field of a set corresponds to an input variance σ 2 = 0.001i. In order to get enough data to achieve the stabilization of linear regressions we generate between 10 and 40 independant sets. The number of generated sets for each grid resolution and l c value is displayed in TableI. For each value of σ 2 ∈ [0.001,0.01] we then obtain sets composed of 10 to 40 realizations.

FIG. 4 . 13 FIG. 5 .

 4135 FIG. 4. Regressions for l c = 6δx and 10δx on 257 3 grid.

14 FIG. 6 .

 146 FIG. 6. Regression slopes of ln(η eff ) vs (real σ ) 2 /2 for different values of l c and resolutions. Top : slope estimations; bottom: slope estimations with 95% confidence intervals.

TABLE I .

 I Number of realizations for each case.

	Resolution			129 3						193 3
	l c /δx	1	2	3	4	5	1	2	3	4	5	6	7
	No.sets	40	40	40	40	40	20	20	20	40		40	40
	Resolution					257 3
	l c /δx	1	2	3	4	5	6	7	8	9	1
	No.sets

because the theory is based on the series expansion of σ 2 . Furthermore, if we use low correlation lengths with σ 2 /2 > 0.22, the high values of ∇η can prevent the Stokes solver from converging.

We note that the variance of the observed effective viscosity depends strongly on the correlation length. This imposes severe constraints regarding the size of the computational grid.

In the case of the effective conductivity or permeability of random media, the power averaging formula (which is exact in one or two dimensions) is proved to be a very efficient approximation in three dimensions, up to log-conductivity variances of seven corresponding to variations of local conductivity of several orders of magnitude [START_REF] Neuman | Prediction of steady state flow in nonuniform geologic media by conditional moments: Exact nonlocal formalism, effective conductivities, and weak approximation[END_REF]. The results of this paper demonstrate that this approximation appears less robust in the context of the Stokes equation. A possible explanation is that due to its 1/q 2 factor, the Stokes equation propagator over-amplifies low frequency fluctuations.

Our observation is also reminiscent of the observed wide fluctuations of sedimentation velocity in suspensions and the appearance of large structures (see Ref. [START_REF] Mucha | A model for velocity fluctuations in sedimentation[END_REF] and references therein).

Finally, with reference to the known rheological literature, the power average formula with averaging power of 1/5 does not compare well with analogous results as provided by the quarter power mixing rule -1/4 [START_REF] Taylor | The viscosity of a fluid containing small drops of another fluid[END_REF][START_REF] Danov | On the viscosity of dilute emulsions[END_REF][START_REF]A method for predicting the performance of unstable miscible displacement in heterogeneous media[END_REF][START_REF] Blackwell | Factors influencing the efficiency of miscible displacement[END_REF]. We note that in the case of the effective viscosity of emulsions, the strong capillarity assumption which ensures overall sphericity of the bubbles changes the local boundary conditions of the flow at the bubbles' boundaries. As a result, an increased energy dissipation is expected, and thus a greater effective viscosity [START_REF] Taylor | The viscosity of a fluid containing small drops of another fluid[END_REF][START_REF] Danov | On the viscosity of dilute emulsions[END_REF] is obtained. This happens even if both fluids are assumed to share an almost common viscosity in order to be consistent with our calculations.

The effective viscosity averaging formula derived in this paper may be considered as a first estimator as far as minimal information is provided regarding the microstructure of a mixture.

Improvement can be obtained if additional information becomes available. From the point of view of computational fluid dynamics (CFD), the proposed formula can be used to test numerical methods aiming at solving Stokes equations with variable viscosity. These numerical issues are the source of major progress toward the development of micro-and nanofluidics modeling tools.

APPENDIX 1. Evaluation of tensor H in the isotropic case

In the case of isotropic correlations, the following quantity must be computed in order to obtain explicit second order results for the effective viscosity:

Due to the isotropy, the angular integral does not depend on q. The term C(r = 0) can then be recognized by computing the inverse Fourier transform. In order to compute the full tensor components, several integrals over the unit sphere must be computed:

By isotropy, the integral of 1 4π S 1 d 2 qi qi i qj is proportional to the unit tensor, so it is sufficient to evaluate its trace. We therefore obtain

The same type of reasoning can be followed for the evaluation of the other term 1 4π S 1 d 2 qi qi i qj qk ql , which is obviously symmetric over all its pairs of indices. The expression is thus proportional to the tensor δ ij δ kl + δ ik δ jl + δ il δ jk . The proportionality constant may be determined by evaluating its partial trace over any pair of indices. We finally obtain

Returning to tensor H, we then obtain

which is Eq. ( 33).

Tensor contractions

In this Appendix, we give some algebraic details about the simplification of Eq. ( 30), which leads to formula Eq. ( 35) when combined with Eq. ( 33). Using Eq. ( 30) combined with Eq. ( 33) (which is the sum of three elementary tensors), we observe that several tensor contractions must be carried out. The main steps are sketched below:

δ ν 1 β δ ν 2 γ iq ν 2 (1q q) γ δ η 0 q 2 +iq γ

(1q q) ν 2 δ η 0 q 2 = -1 5 iq ν 2 (1q q) αγ η 0 q 2 + iq γ (1q q) αν 2 η 0 q 2 iq ν 2 (1q q) γ δ η 0 q 2 + iq γ (1q q) ν 2 δ η 0 q 2 = 2 5

(1q q) αδ η 2 0 q 2 .

(A6)

Throughout this calculation, we have removed the contraction with f δ (q), which is of no interest.

Power average of log-normal variables

The general equality (36) may be derived using Gaussian integrals for log-normal variables:

We start from

The estimation of the average can be carried out, using the probability density associated with random variable ζ (with average 0 and unit variance):

where the normalization parameter N ensures that R 1 N e -x 2 /2 dx = 1.

004100-17

We then obtain = e (ωσ ) 2 /2 , which leads to Eq. [START_REF] Chatelin | A parametric study of mucociliary transport by numerical simulations of 3D non-homogeneous mucus[END_REF]. A direct application of previous calculations shows that η 2 = η 2 g exp(σ 2 ), η 2 = η 2 g exp(2σ 2 ).

The variance C(0) of η is then given by

For small log viscosity variance σ 2 , one obtains a direct relation between the variance of viscosity and the variance of the corresponding logarithm:

(A10)


Discrete element models for grain breakage

LGCIE, INSA Lyon, Villeurbanne, France, francois.nader@insa-lyon.fr b. LGCIE, INSA Lyon, Villeurbanne, France, claire.silvani@insa-lyon.fr c. LGCIE, INSA Lyon, Villeurbanne, France, irini.djeran-maigre@insa-lyon.fr d. LAMCOS, INSA Lyon, Villeurbanne, France, guilhem.mollon@insa-lyon.fr

Abstract:

The study of rockfill behaviour shows important experimental difficulties, due to the big grain size and the need of big-dimension apparatus to test and characterize these materials. The empirical methods used in building rockfill structures, the complexity of compaction procedures, and the effects of water and temperature, can result in particle breakage, which may lead to considerable, and sometimes dangerous, long term settlements. Numerical approaches are often adopted to describe and predict the behaviour of this material, offering the possibility to test various parameters and evaluate the effect of each parameter on the general behaviour of the material. The discrete element method is often used to represent the behaviour of such discontinuous materials, distinguished by grain breakage, taking in input physical parameters. In the work presented here, we adopted this type of simulation to model grain breakage in oedometric tests. Many models are used in literature to describe grain breakage in numerical discrete simulations. These models can be divided into two groups. In the first group, a breakage criterion is set consisting of a threshold value of a stress or a force, that may vary in accordance to the size of the grain following Weibull's distribution [START_REF] Weibull | A statistical distribution function of wide applicability[END_REF]. When a grainconsidered as one object-breaks, it is replaced by a certain number of smaller grains (Figure 1), either conserving the total mass, or losing part of it in the breakage. Most models are in 2D, with a few 3D models. Using a smooth contact method, in which overlapping generates a repulsion force, Åström and Herrmann [START_REF] Åström | Fragmentation of grains in a two-dimensional packing[END_REF] and Tsoungui [START_REF] Tsoungui | Numerical model of crushing of grains inside twodimensional granular materials[END_REF] tested two breakage criteria: the threshold value for the pressure on a grain, and the threshold value for the largest compressive contact force on a grain, replacing a broken disc with a number of smaller discs, conserving the total mass. Ben-Nun and Einav [START_REF] Ben-Nun | Force attractor in confined comminution of granular materials[END_REF] suggested a different technique, consisting of computing the mean force of the normal contact forces on a disc, and when the threshold is exceeded, the broken disc is replaced by discs small enough to avoid overlapping, then expanding their radius to conserve the mass. McDowell and De Bono [START_REF] Mcdowell | On the micro mechanics of one-dimensional normal compression[END_REF], dealt with 3D simulations, with a sphere model following the same breakage configuration as Åström and Herrmann [START_REF] Åström | Fragmentation of grains in a two-dimensional packing[END_REF], setting a threshold value for the induced tensile stress inside a sphere, and when broken, replacing it with two spheres overlapping, along the direction of the minor principal stress. Cantor et al. [START_REF] Cantor | New approach to grain fragmentation for discrete element methods[END_REF] simulated polygonal rigid grains that can split in the direction of the principal stress when the induced tensile stress inside the polygon 𝜎 𝑡 𝑛𝑒𝑡 exceeds a threshold value 𝜎 𝑐𝑟𝑖𝑡 determined experimentally. Using Abaqus simulations, they established a formula that allows estimating the induced tensile stress inside a polygon, taking its geometry into account. Lobo-Guerrero and Vallejo [START_REF] Lobo-Guerrero | Discrete element method evaluation of granular crushing under direct shear test conditions[END_REF] simulated discs that break when both the induced tensile strength 𝜎 𝑡 exceeds a threshold value 𝜎 𝑚𝑎𝑥 , and the coordination number of the disc is less or equal to 3. The replacing discs do not overlap, however, mass is lost in the procedure. In the second group of models, a grain is considered as an assembly of smaller particles, joined together by a cohesive bond (Figure 2). Breakage of these contact bonds leads to the breakage of the grain; no other breakage criterion is needed. Cheng et al. [START_REF] Cheng | Discrete element simulation of crushable soil[END_REF] joined spherical particles to create a spherical grain, and removed 20% of the particles to add defects to the grain. Silvani et al. [START_REF] Silvani | Discrete modelling of time-dependent rockfill behavior[END_REF] simulated two-dimensional oedometric tests using rigid discs bonded together by a time-dependent Mohr Coulomb interaction law, to reproduce long-term behaviour. Nguyen et al. [START_REF] Nguyen | Effects of shape and size polydispersity on strength properties of granular materials[END_REF] split a circular grain into smaller particles using Voronoï tessellation, and also applied a Mohr Coulomb type cohesive bond between particles.

Figure 1 -Grain breakage using the replacement method [START_REF] Cantor | New approach to grain fragmentation for discrete element methods[END_REF] Figure 2 -Grain breakage using the agglomerate method [START_REF] Cantor | New approach to grain fragmentation for discrete element methods[END_REF] To simulate grain breakage, we have also adopted the discrete element method (DEM), and more specifically the Non-Smooth Contact Dynamics (NSCD) method [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid bodies collections[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF] implemented in the software LMGC 90 . The goal of the NSCD method is to obtain the overall behaviour of a collection of objects by considering the dynamics of each element, taking into account the interactions between bodies. The problem is studied at two levels: movement equations of the particles expressed using global variables (kinematic variables) and the interactions described using local variables (contact variables). Today, there are many 2D models. In the work presented here, we enriched the modelling creating 3D grain models using the concept of particles assembly, and cohesive bonds. The first model consists of applying a cohesive bond between spherical particles, in order to create a spherical grain (Figure 3).

In order to find a model closer to reality, a second grain model was generated, with tetrahedral particles (Figure 4). Three input parameters are therefore set: the number of vertices and the minimal and maximal dimensions of the grain R min and R max . Using spherical coordinates, we place the vertices on semi-circles in space, at a distance between R min and R max from the centre (Figure 5). The grains are then divided into tetrahedral particles using an automatic mesh generator, and imported into LMGC 90 with a cohesive bond applied between the particles. A third model, inspired by the works of Mollon et al. [START_REF] Mollon | Generating realistic 3D sand particles using Fourier descriptors[END_REF], is also used: polyhedrons are generated and divided into Voronoï cells (Figure 6).

To validate the numerical simulations, oedometric tests were carried out in collaboration with IRSTEA of Aix-En-Provence, on limestone ballast of dimensions between 20 and 40mm (Figure 7). Correlation and stereo correlation techniques were used, to study displacement and breakage fields. Simulations of a single grain crushing test were performed with the different models. The tests showed that the grain strength depends on its dimensions, particles dimensions and cohesive strength of the bond joining the particles, as expected. Then, many samples were created, of various dimensions and shapes, in order to simulate the experimental oedometric tests. A parametric study is carried out to quantify the effect of factors influencing the grain breakage, under imposed displacement and imposed loading. The number of cohesive bonds is monitored throughout the simulations, in order to establish the grain breakage ratio and granular size distributions.

F.

  Nader a , C. Silvani b , I. Djeran-Maigre c , G. Mollon d a.

Figure 3 -Figure 5 -Figure 7 -

 357 Figure 3 -Grain composed of spherical particles Figure 4 -Grain composed of tetrahedral particles