Van-Sang Ngo

DAMPING EFFECTS IN BOUNDARY LAYERS FOR ROTATING FLUIDS WITH SMALL VISCOSITY

Keywords: 1991 Mathematics Subject Classification. 76D03, 76D05, 76U05 Damping Euler, Navier-Stokes, Rotating fluids, Ekman boundary layer

The goal of this paper is to study the system of rotating fluids between two infinite parallel plates with Dirichlet boundary conditions and with small viscosity which vanishes when the Rossby number goes to zero. We want to improve the convergence result of [16] and show the global in time convergence of the weak solution of the system of rotating fluids towards the solution of a twodimensional damped Euler system with three components, using the decay in time of the H s -norm (s > 2) of the limiting solution.

Introduction

In this paper, we consider the system of fast rotating incompressible fluids (1.1)

       ∂ t u ε -ν h (ε)∆ h u ε -βε∂ 2 3 u ε + u ε • ∇u ε + e 3 ∧ u ε ε + ∇p ε = 0 div u ε = 0
u ε |t=0 = u ε 0 , for t > 0 and for x ∈ Ω = R 2 ×]0, 1[, coupled with the Dirichlet boundary conditions

u ε (t, x 1 , x 2 , 0) = u ε (t, x 1 , x 2 , 1) = 0.
Here, the fluid rotates around e 3 -axis and ε > 0 stands for the Rossby number, which will be small and goes to zero. In this system, the vertical viscosity is ν v = βε, with β > 0 and the horizontal viscosity ν h = ν h (ε) > 0 depends on ε and goes to zero as ε goes to zero. The classical isotropic diffusion -ν∆ is replaced by an anisotropic diffusion of the form

-ν h ∆ h -ν v ∂ 2 3 = -ν h (ε) ∂ 2 1 + ∂ 2 2 -βε∂ 2 3
, where ∂ i is the partial derivative with respect to the variable x i . This anisotropy is often considered in meteorology and oceanography, where vertical viscosity ν v is taken to be much smaller than the horizontal viscosity ν h . We refer to Pedlosky [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF] and Greenspan [START_REF] Greenspan | The theory of rotating fluids[END_REF] for more details about physical considerations.

It is clear that the Coriolis force term has no contribution in the energy estimates for the system (1.1). In the case where ν h > 0 and ν v > 0, independent of ε, the global existence of weak solutions of (1.1) in the sense of Leray can be proved by the same method as in the case of the Navier-Stokes equations (see J. Leray, [START_REF] Leray | Essai sur le mouvement d'un liquide visqueux emplissant l'espace[END_REF]). By using the same arguments as in H. Fujika and T. Kato [START_REF] Fujita | On the Navier-Stokes initial value problem I[END_REF] for the Navier-Stokes equations, one can show local existence and uniqueness of a strong solution of (1.1), the lifespan of which is bounded from below, uniformly with respect to ε. We also want to refer to [START_REF] Chemin | Anisotropy and dispersion in rotating fluids, Nonlinear Partial Differential Equations and their application[END_REF], [START_REF] Chemin | Fluids with anisotropic viscosity, Special issue for R. Temam's 60th birthday[END_REF], [START_REF] Iftimie | The resolution of the Navier-Stokes equations in anisotropic spaces[END_REF], [START_REF] Paicu | Étude asymptotique pour les fluides anisotropes en rotation rapide dans le cas périodique[END_REF], [START_REF] Paicu | Équation anisotrope de Navier-Stokes dans des espaces critiques[END_REF], . . . , and the references therein for the anisotropic cases where the vertical viscosity is zero. In the case where ν h and ν v vanish or depend on ε, using the same method as for the Euler system, we can prove that the system (1.1) is locally well-posed when the initial data u ε 0 belong to a smooth enough Sobolev space, independent of ν h and ν v .

We have not mentioned the role of the Coriolis force yet. Actually in fast rotating systems with small Rossby number, the Coriolis force plays an important role. In the experiment of G.I. Taylor (see [START_REF] Cushman-Roisin | Introduction to geophysical fluid dynamics[END_REF]), drops of dye injected into a rapidly rotating, homogeneous fluid, within a few rotations, formed perfectly vertical sheets of dyed fluid, known as Taylor curtains. In large-scale atmospheric and oceanic flows, the fluid motions also have a tendency towards columnar behaviors (Taylor columns). For example, currents in the western North Atlantic have been observed to extend vertically over several thousands meters without significant change in amplitude and direction ([26]).

Fast rotating flows have first been considered in periodic domains. In the work of Babin, Mahalov and Nicolaenko [START_REF] Babin | Global Splitting, Integrability and Regularity of 3D Euler and Navier-Stokes Equations for Uniformly Rotating Fluids[END_REF], [START_REF] Babin | Global regularity of 3D rotating Navier-Stokes equations for resonant domains[END_REF], the authors studied the isotropic case, where the diffusion term is -ν∆, with ν > 0 independent of ε. For any initial data in appropriate Sobolev spaces, they proved the global existence and uniqueness of a smooth solution when the rotation is fast enough (that is, the Rossby number ε is small enough). The same result was proven by Gallagher in [START_REF] Gallagher | Applications of Schochet's Methods to Parabolic Equation[END_REF], using the method of Schochet [START_REF] Schochet | Fast singular limits of hyperbolic PDEs[END_REF], for fast rotating fluid systems and also for more general parabolic systems. It was also shown in [START_REF] Babin | Global Splitting, Integrability and Regularity of 3D Euler and Navier-Stokes Equations for Uniformly Rotating Fluids[END_REF], [START_REF] Babin | Global regularity of 3D rotating Navier-Stokes equations for resonant domains[END_REF] and [START_REF] Gallagher | Applications of Schochet's Methods to Parabolic Equation[END_REF] that, when ε goes to zero, the fast rotating system converges to a two-dimensional Navier-Stokes system with three components, if ν h > 0 is fixed, or to a two-dimensional Euler system with three components, if ν h is zero or goes to zero as ε goes to zero. The convergence in the more general case, where the rotation axis is fixed but not the amplitude, was proven in [START_REF] Gallagher | Weak convergence results for inhomogeneous rotating fluid equations[END_REF] by Gallagher and Saint-Raymond. The anisotropic case, with zero vertical viscosity (ν v = 0), in periodic domains, that is the system (1.1), with ν h > 0 independent of ε, in the torus T 3 , was studied later by Paicu in [START_REF] Paicu | Étude asymptotique pour les fluides anisotropes en rotation rapide dans le cas périodique[END_REF].

In the case of rotating fluids in the whole space R 3 , it was proved in [START_REF] Chemin | Anisotropy and dispersion in rotating fluids, Nonlinear Partial Differential Equations and their application[END_REF], [START_REF] Chemin | Fluids with anisotropic viscosity, Special issue for R. Temam's 60th birthday[END_REF] and [START_REF] Gallagher | Weak convergence results for inhomogeneous rotating fluid equations[END_REF] that, if the initial data are divergence-free vector fields belonging to H 0,s (R 3 ) then, when ε goes to zero, the limiting system corresponding to (1.1) is zero. We remark that in this paper, in order to simplify the notations, we use the bold character X to indicate the space of vector fields, each component of which belongs to the space X. The main idea consists in proving and using Strichartz estimates for the associated linear system, which shows that, in the case of R 3 , the free oscillation waves of the fluid propagate to infinity and the energy of the fluid decreases. Using this property, it was also proved in [START_REF] Chemin | Fluids with anisotropic viscosity, Special issue for R. Temam's 60th birthday[END_REF] that the system (1.1), with ν h > 0 fixed, has a unique, global solution when the rotation is fast enough. The same result in the case where ν h slowly goes to zero as ε → 0, say ν h = ε α , with α ∈]0, α 0 ], for a certain α 0 > 0, was proved in [START_REF] Ngo | Rotating Fluids with small viscosity[END_REF]. If the initial data is not in H 0,s (R 3 ) but is the sum of a 2D part in L 2 (R 2 h ) 3 and a 3D part in H 0,s (R 3 ), s > 1 2 , then the limit system is not zero but a 2D Navier-Stokes system (ν h > 0 fixe) or a 2D Euler system (ν h → 0) with three components, as ε goes to zero. Using the global wellposedness of the limit system, the global existence of a strong solution when ε is small enough can be proved in the case of fixed ν h > 0 (see [START_REF] Chemin | Fluids with anisotropic viscosity, Special issue for R. Temam's 60th birthday[END_REF]).

The case of a domain between two parallel plates (with Dirichlet boundary conditions) is very different from the case of a domain without boundary (R 3 or T 3 ) mentioned above. Indeed, when the rotation goes to infinity (ε → 0), the fluid has the tendency to have a columnar behavior (the Taylor-Proudman theorem). However, the Taylor columns are only formed in the interior of the domain. Near the boundary, the Taylor columns are destroyed and very thin boundary layers are formed (Ekman boundary layers). Inside the boundary layers, the behavior of the fluid becomes very complex and the friction slows the fluid down in a way that the velocity is zero on the boundary. As a consequence, in the limiting system, an additional damping term of the form γu, γ > 0, appears. The coefficient γ was proved to be √ 2β in [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF] (we will explain the calculation of this coefficient in the appendix). This phenomenon is well known in fluid mechanics as the Ekman pumping.

For any ε > 0, and for any u ε 0 ∈ L 2 (Ω), the system (1.1) possesses a weak Leray solution

u ε ∈ L ∞ (R + , L 2 (Ω)) ∩ L 2 (R + , Ḣ1 (Ω)).
If ν h > 0 is fixed, then taking into account the Ekman pumping, the limiting system is the following 2D damped Navier-Stokes sytem

(1.2)        ∂ t u h -ν h ∆ h u h + u h • ∇ h u h + 2β u h = -∇ h p ∂ 3 u = 0, div h u h = 0, u 3 = 0 u |t=0 = u 0 = (u 1 0 , u 2 0 , 0). When ε → 0, u ε was proved to converge to the solution of the limiting system in L ∞ (R + , L 2 (Ω)), for Ω = T 2 h × [0, 1]
, by Grenier and Masmoudi [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF] and by Masmoudi [20] in the case where the initial data are well prepared (i.e. lim ε→0

u ε 0 = u 0 = (u 1 0 (x 1 , x 2 ), u 2 0 (x 1 , x 2 ), 0) in L 2 (T 2 h × [0, 1]
) and by Masmoudi [START_REF] Masmoudi | Ekman layers of rotating fluids: the case of general initial data[END_REF] in the case of general initial data. The same result for Ω = R 2 h × [0, 1] was prove by Chemin et al. in [START_REF] Chemin | Ekman boundary layers in rotating fluids, ESAIM Controle optimal et calcul des variations[END_REF], using Strichartz estimates.

If ν h = ν h (ε) → 0 when ε → 0, the limiting system is now the following 2D damped Euler system

(1.3)        ∂ t u h + u h • ∇ h u h + 2β u h = -∇ h p ∂ 3 u = 0, div h u h = 0, u 3 = 0 u |t=0 = u 0 = (u 1 0 , u 2 0 , 0).
In [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF], [START_REF] Masmoudi | The Euler Limit of the Navier-Stokes Equations, and Rotating Fluids with Boundary[END_REF], [START_REF] Masmoudi | Ekman layers of rotating fluids: the case of general initial data[END_REF] (for Ω = T 2 h × [0, 1]) and in [START_REF] Chemin | Ekman boundary layers in rotating fluids, ESAIM Controle optimal et calcul des variations[END_REF] (for Ω = R 2 h × [0, 1]), the solution u ε can still be proved to converge to the solution of (1.3) in L ∞ loc (R + , L 2 (Ω)). The local convergence with respect to the time variable is due to the local in time boundedness of the H σ (R 2 h )-norm (σ > 2) of the solution of (1.3) used in the previous results.

The goal of this paper is to improve the above results and to prove that, in the case where ν h = ν h (ε) goes to zero as ε goes to zero, the convergence of u ε towards u is global in time. We consider the system (1.1) in the domain Ω = R 2 h × [0, 1] (the same result can be proved for the domain Ω = T 2 h × [0, 1] with minor changes) with Dirichlet boundary conditions and with the well-prepared data as in [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF] or [START_REF] Masmoudi | The Euler Limit of the Navier-Stokes Equations, and Rotating Fluids with Boundary[END_REF]. The main result of the paper is the following theorem.

Theorem 1.1. We suppose that ν h (ε) → 0 as ε → 0. Let u ε 0 ∈ L 2 (R 2 h × [0, 1]
) be a family of initial data such that lim

ε→0 u ε 0 = u 0 = (u 1 0 (x 1 , x 2 ), u 2 0 (x 1 , x 2 ), 0) in L 2 (R 2 h × [0, 1]),
where u 0 is a divergence-free two-dimensional vector field in

H σ (R 2 h ), σ > 2.
Let u be the solution of the limiting system (1.3) with initial data u 0 and, for each ε > 0, let u ε be a weak solution of (1.1) with initial data u ε 0 . There exists a positive constant

C = C(u 0 ) > 0 such that, if ν h (ε) ≥ C(u 0 ) ε, then lim ε→0 u ε -u L ∞ (R+,L 2 (R 2 h ×[0,1])) = 0.
The key of the proof of Theorem 1.1 consists in the careful study of the limiting system (1.3). We prove that the damping term √ 2β u in the system (1.3) plays a very important role which allow us to prove the exponential decay in time of the H σ (R 2 h )-norm, σ > 2, of the solution. More precisely, in Section 3, we prove the following theorem.

Theorem 1.2. Let u 0 = (u 1 0 (x 1 , x 2 ), u 2 0 (x 1 , x 2 ), 0) ∈ L 2 (R 2
h ) be a divergence-free vector field, the horizontal vorticity of which

w 0 = ∂ 1 u 2 0 -∂ 2 u 1 0 ∈ L 2 (R 2 h ) ∩ L ∞ (R 2 h
). Then, the system (1.3), with initial data u 0 , has a unique, global solution

u ∈ C(R + , L 2 (R 2 h )) ∩ L ∞ (R + , L 2 (R 2 h )). Moreover, if u 0 belongs to H σ (R 2
h ), σ > 2, then there exists a positive constant C depending on β and u 0 H σ (R 2 h ) such that, for any t > 0,

u(t) H σ (R 2 h ) ≤ Ce -t √ 2β .
We note that a unique global solution of (1.3) exists thanks to the Yudovitch theorem [START_REF] Yudovitch | Non stationnary flows of an ideal incompressible fluid[END_REF] and the exponential decay is not true in the case of the two-dimensional non-damped Euler system, where the L ∞ -norm of the gradient of the vorticity may have a double exponential growth in time (see Bahouri & Chemin [START_REF] Bahouri | Équations de transport relatives à des champs de vecteurs non-lipschitziens et mécanique des fluides [Transport equations for non-Lipschitz vector fields and fluid mechanics[END_REF], Kiselev & Šverák [START_REF] Kiselev | Small scale creation for solutions of the incompressible two-dimensional Euler equation[END_REF] and the references therein). The goal of Section 3 is to prove the exponential decay properties of this solution in the case of damped Euler system.

The paper is organized as follows. In the next section, we briefly recall the dyadic decompositions and the Littlewood-Paley theory which will be used in Section 3 to prove Theorem 1.2. In the last section, we prove the main theorem (Theorem 1.1) on the convergence of weak solutions of the system (1.1) towards the solution of the system (1.3). Finally, for the convenience of the reader, we will very briefly recall the construction of Ekman boundary layers and the formation of the limiting system in the appendix. We will also explain why the damping term appears in the limiting system (the Ekman pumping).

Brief recall of dyadic decompositions

In this section, we briefly recall the properties of dyadic decompositions in the Fourier space and give some elements of the Littlewood-Paley theory. The complete details can be found in [START_REF] Bony | Calcul symbolique et propagation des singularités pour les équations aux dérivées partielles non linéaires[END_REF], [START_REF] Chemin | Fluides parfaits incompressibles[END_REF] and [START_REF] Chemin | Flot de champs de vecteurs non Lipschitziens et équations de Navier-Stokes[END_REF]. Let F and F -1 be the Fourier transform and its inverse, and that we also write u = Fu. For any d ∈ N * and 0 < r < R, we denote

B d (0, R) = ξ ∈ R d | |ξ| ≤ R and C d (r, R) = ξ ∈ R d | r ≤ |ξ| ≤ R .
The following Bernstein lemma gives important properties of a distribution u when its Fourier transform is well localized. We refer the reader to [START_REF] Chemin | Fluides parfaits incompressibles[END_REF] for the proof of this lemma.

Lemma 2.1. Let k ∈ N, d ∈ N * and r 1 , r 2 ∈ R satisfy 0 < r 1 < r 2 .
There exists a constant C > 0 such that, for any a, b ∈ R, 1 ≤ a ≤ b ≤ +∞, for any λ > 0 and for any u ∈ L a (R d ), we have

(2.1) supp ( u) ⊂ B d (0, r 1 λ) =⇒ sup |α|=k ∂ α u L b ≤ C k λ k+d( 1 a -1 b ) u L a , and
(2.2) supp ( u) ⊂ C d (r 1 λ, r 2 λ) =⇒ C -k λ k u L a ≤ sup |α|=k ∂ α u L a ≤ C k λ k u L a .
Let ψ be an even smooth function in C ∞ 0 (R), whose support is contained in the ball B 1 (0, 4 3 ), such that ψ is equal to 1 on a neighborhood of the ball B 1 (0, 3 4 ). Let

ϕ(z) = ψ z 2 -ψ(z).
Then, the support of ϕ is contained in the ring

C 1 ( 3 4 , 8
3 ), and ϕ is identically equal to 1 on the ring

C 1 ( 4 3 , 3
2 ). The functions ψ and ϕ allow us to define a dyadic partition of R d , d ∈ N * , as follows

∀z ∈ R, ψ(z) + j∈N ϕ(2 -j z) = 1.
Moreover, this decomposition is almost orthogonal, in the sense that, if |j -

j | ≥ 2, then supp ϕ(2 -j (•)) ∩ supp ϕ(2 -j (•)) = ∅.
We introduce the following dyadic frequency cut-off operators. We refer to [START_REF] Bony | Calcul symbolique et propagation des singularités pour les équations aux dérivées partielles non linéaires[END_REF] and [START_REF] Chemin | Fluides parfaits incompressibles[END_REF] for more details.

Definition 2.2. For any d ∈ N * and for any tempered distribution u ∈ S (R d ), we set

∆ q u = F -1 ϕ(2 -q |ξ|) u(ξ) , ∀q ∈ N, ∆ -1 u = F -1 (ψ(|ξ|) u(ξ)) , ∆ q u = 0, ∀q ≤ -2, S q u = q ≤q-1 ∆ q u, ∀q ≥ 1.
Using the properties of ψ and ϕ, one can prove that for any tempered distribution u ∈ S (R d ), we have

u = q≥-1 ∆ q u in S (R d ),
and the (isotropic) nonhomogeneous Sobolev spaces H s (R d ), with s ∈ R, can be characterized as follows

Proposition 2.3. Let d ∈ N * , s ∈ R and u ∈ H s (R d ). Then, u H s := R d (1 + |ξ| 2 ) s | u(ξ)| 2 dξ 1 2 ∼   q≥-1 2 2qs ∆ q u 2 L 2   1 2
Moreover, there exists a square-summable sequence of positive numbers {c q (u)} with q c q (u)

2 = 1, such that ∆ q u L 2 ≤ c q (u)2 -qs u H s .
In what follows, we also use separate dyadic decompositions in the horizontal and vertical directions.

Definition 2.4. For any tempered distribution u ∈ S (R 3 ), we set

∆ h q u = F -1 ϕ(2 -q |ξ h |) u(ξ) , ∆ v q u = F -1 ϕ(2 -q |ξ 3 |) u(ξ) , ∀q ∈ N, ∆ h -1 u = F -1 (ψ(|ξ h |) u(ξ)) , ∆ v -1 u = F -1 (ψ(|ξ 3 |) u(ξ)) , ∆ h q u = 0, ∆ v q u = 0, ∀q ≤ -2, S h q u = q ≤q-1 ∆ h q u, S v q u = q ≤q-1 ∆ v q u, ∀q ≥ 1.

Exponential decay for the 2D limiting system

In this section, we prove the exponential decay in time of the solution of the 2D limiting system (1.3). We recall that throughout this paper, for any vector field u = (u 1 , u 2 , u 3 ) independent of the vertical variable x 3 , we denote by w the associated horizontal vorticity, w = ∂ 1 u 2 -∂ 2 u 1 . We also set γ = √ 2β to simplify the notation.

Lemma 3.1. Let u 0 = (u 1 0 (x 1 , x 2 ), u 2 0 (x 1 , x 2 ), u 3 0 (x 1 , x 2 )) ∈ L 2 (R 2
h ) be a divergence-free vector field, the horizontal vorticity of which

w 0 = ∂ 1 u 2 0 -∂ 2 u 1 0 ∈ L 2 (R 2 h ) ∩ L ∞ (R 2 h
). Then, the system (1.3), with initial data u 0 , has a unique, global solution

u ∈ C(R + , L 2 (R 2 h )) ∩ L ∞ (R + , L 2 (R 2 h )).

Moreover,

(i) There exists a constant C > 0 such that, for any p ≥ 2 and for any t > 0, we have

∇ h u h (t) L p (R 2 h ) ≤ CM p e -γt , (3.1)
u h (t) L p (R 2 h ) ≤ CM e -γt , (3.2)
where

M = max u h 0 L 2 (R 2 h ) , w 0 L 2 (R 2 h ) , w 0 L ∞ (R 2 h ) . (ii) For any p ≥ 2, if u 3 0 ∈ L p (R 2 h ), then, (3.3) u 3 (t) L p (R 2 h ) ≤ u 3 0 L p (R 2
h ) e -γt .

Proof

In (1.3), the first two components of u satisfy a two-dimensional damped Euler system. Then, thanks to the Yudovitch theorem [START_REF] Yudovitch | Non stationnary flows of an ideal incompressible fluid[END_REF] (see also [START_REF] Chemin | Fluides parfaits incompressibles[END_REF]), this system has a unique solution

u h ∈ C R + , L 2 (R 2 h ) ∩ L ∞ R + , L 2 (R 2 h ) such that the horizontal vorticity w ∈ L ∞ (R + , L 2 (R 2 h )) ∩ L ∞ (R + , L ∞ (R 2 h )).
Since the third component u 3 satisfies a linear transport-type equation, we can deduce the existence and uniqueness of the solution u of the limiting system (1.3).

By definition, the horizontal vorticity w satisfies the following equation

(3.4) ∂ t w + γw + u h • ∇ h w = 0
Taking the L 2 scalar product of (3.4) with |w| p-2 w, we get

1 p d dt w p L p (R 2 h ) + γ w p L p (R 2 h ) = 0. By using an interpolation between L 2 (R 2 h ) and L ∞ (R 2
h ), we deduce from the above equation that

(3.5) w(t) L p (R 2 h ) ≤ w 0 L p (R 2 h ) e -γt ≤ C w 0 2 p L 2 (R 2 h ) w 0 1-2 p L ∞ (R 2
h ) e -γt ≤ CM e -γt . We recall that ∇ h u h = R w where R is an homogeneous Caldéron-Zygmund operator of order 0. According to [ [START_REF] Chemin | Fluides parfaits incompressibles[END_REF], Theorem 3.1.1], for any p > 1 and for any t > 0,

∇ h u h (t) L p (R 2 h ) ≤ C p 2 p -1 w(t) L p (R 2 h
) , which implies (3.1), for any p ≥ 2. In particular, we have

(3.6) u h (t) Ḣ1 (R 2 h ) = ∇ h u h (t) L 2 (R 2 h ) = w(t) L 2 (R 2 h ) ≤ CM e -γt .
We remark that Inequality (3.2), in the case where p = 2, directly comes from the energy estimate for the damped Euler systems. So, using the Sobolev embedding

Ḣ p-2 p (R 2 h ) → L p (R 2 h ) and the fact that Ḣ p-2 p (R 2 h
) is an interpolated space between L 2 (R 2 h ) and Ḣ1 (R 2 h ), we deduce (3.2) from Estimate (3.6), when p ≥ 2. Since the vertical component u 3 satisfies the same linear transport equation as w, Inequality (3.3) can be proved in the same way as (3.5).

Lemma 3.2. Let r > 1.

There is a constant C > 0 such that, for any divergence-free vector fields u in

L 2 (R 2 h ), the horizontal gradient ∇ h u h ∈ L ∞ (R 2 h ) and the vorticity w = ∂ 1 u 2 -∂ 2 u 1 ∈ H r (R 2 h ), we have (3.7) ∆ h q (u h • ∇ h w) ∆ h q w ≤ Cd q 2 -2qr ∇ h u h L ∞ + w L ∞ w 2 H r .
where (d q ) q≥-1 is a summable sequence of positive numbers, the 1 -norm of which is independent of u.

In order to prove Lemma 3.2, we need the following estimate (the proof of which can be found in [START_REF] Paicu | Équation anisotrope de Navier-Stokes dans des espaces critiques[END_REF]). Let [., .] denote the usual commutator. Lemma 3.3. Let d ∈ N * . There exists a constant C > 0 such that, for any tempered distributions u, v in S (R d ), we have

[∆ q , u] v L 2 := ∆ q (uv) -u∆ q v L 2 ≤ C2 -q ∇u L ∞ v L 2 .

Proof of Lemma 3.2

First of all, using the anisotropic Bony decomposition into paraproducts and remainders as in Definition 2.4, we can write

(3.8) ∆ h q (u h • ∇ h w) ∆ h q w ≤ I 1 q (w) + I 2 q (w)
, where

I 1 q (w) = ∆ h q q -q≥N S h q +2 (∇ h w)∆ h q u h ∆ h q w , and
I 2 q (w) = ∆ h q |q -q|≤N S h q -1 u h ∆ h q (∇ h w) ∆ h q w ,
and where N > 0 is a fixed large enough integer.

Using Cauchy-Schwarz and Hölder inequalities, we get (3.9)

I 1 q (w) ≤ q -q≥N S h q +2 (∇ h w) L ∞ ∆ h q u h L 2 ∆ h q w L 2 .
The Bernstein lemma 2.1 implies that

S h q +2 (∇ h w) L ∞ ≤ C2 q S h q +2 w L ∞ ≤ C2 q w L ∞ , and ∆ h q u h L 2 ≤ C2 -q ∆ h q ∇ h u h L 2 .
Taking into account the fact that w = ∂ 1 u 2 -∂ 2 u 1 and div h u h = 0, using an integration by parts, we easily obtain

∆ h q w 2 L 2 = ∆ h q ∂ 1 u 2 2 L 2 + ∆ h q ∂ 2 u 1 2 L 2 -2 ∆ h q ∂ 1 u 2 | ∆ h q ∂ 2 u 1 = ∆ h q ∂ 1 u 2 2 L 2 + ∆ h q ∂ 2 u 1 2 L 2 + ∆ h q ∂ 1 u 1 2 L 2 + ∆ h q ∂ 2 u 2 2 L 2 = ∆ h q ∇ h u h 2 L 2 .
Therefore, using Proposition 2.3, we deduce from (3.9) that

I 1 q (w) ≤ C q -q≥N w L ∞ ∆ h q w L 2 ∆ h q w L 2 ≤ C c q (w) q -q≥N c q (w)2 -(q -q)r 2 -2qr w L ∞ w 2 H r ,
where (c q (w)) q≥-1 is a square-summable sequence of positive numbers such that q c q (w) 2 = 1. We remark that (2 -qr ) q≥-N is a summable sequence. So, as a convolution product of a 2 -sequence and a 1 -sequence,   q -q≥N c q (w)2 -(q -q)r   q is square-summable. Let d q = c q (w)

q -q≥N c q (w)2 -(q -q)r .

It is clear that (d q ) q≥-1 is a summable sequence of positive numbers. Thus, we have (3.10)

I 1 q (w) ≤ Cd q 2 -2qr w L ∞ w 2 H r .
In order to estimate I 2 q (w), we write (3.11) I 2 q (w) ≤ I 2A q (w) + I 2B q (w) + I 2C q (w), where

I 2A q (w) = S h q u h • ∇ h ∆ h q w ∆ h q w , I 2B q (w) = |q -q|≤N S h q -S h q -1 u h • ∇ h ∆ h q w ∆ h q w , and
I 2C q (w) = |q -q|≤N ∆ h q , S h q -1 u h ∇ h ∆ h q w ∆ h q w .
Using integrations by parts, we get (3.12)

I 2A q (w) = 1 2 S h q (div h u h )∆ h q w ∆ h q w = 0.
We remark that when q, q ≥ 1, S h q -S h q -1 does not contain the low Fourier frequencies. So, I 2B q (w) can be bounded in the same way as I 1 q (w), and we have (3.13)

I 2B q (w) ≤ Cd q 2 -2qr w L ∞ w 2 H r .
Finally, using Lemma 3.3 and the same calculations as for I 1 q (w) (since r > 1), we obtain There exists a constant C r > 0 such that, for any divergence-free vector fields

I 2C q (w) ≤ C2 -q |q -q|≤N S h q -1 ∇ h u h L ∞ ∇ h ∆ h q w L 2 ∆ h q w L 2 (3.14) ≤ C2 -2qr ∇ h u h L ∞ w 2 H r   c q (w) |q -q|≤N c q (w)2 -(q -q)(r-1)   ≤ Cd q 2 -2qr ∇ h u h L ∞ w 2
u in L 2 (R 2 h ), the horizontal gradient ∇ h u h ∈ L ∞ (R 2 h ) and the vorticity w = ∂ 1 u 2 -∂ 2 u 1 ∈ H r (R 2
h ), we have

(3.15) ∇ h u h L ∞ (R 2 h ) ≤ C r w L ∞ (R 2 h ) ln e + w H r (R 2 h ) w L ∞ (R 2 h )
.

The lemma 3.4 is only a special case of [ [START_REF] Chemin | Fluides parfaits incompressibles[END_REF], Theorem 3.3.2]. We refer to [START_REF] Chemin | Fluides parfaits incompressibles[END_REF] for a proof of this lemma.

Lemma 3.5. Let r > 1. Under the hypotheses of Lemma 3.1 and the additional hypothesis that w 0 belongs to H r (R 2 h ), there exist positive constants C 1 and C 2 depending on γ and

w 0 H r (R 2 ) such that (3.16) w(t) H r (R 2 h ) ≤ C 1 e -γt , and
(3.17) ∇ h u h (t) L ∞ (R 2 h ) ≤ C 2 e -γt .

Proof

In what follows, we use C to refer to a generic positive constant that may change from line to line. For any r > 1, using Lemma 3.2, we get the following energy estimate in the Sobolev H r -norm:

(3.18) 1 2 d dt w(t) 2 H r + γ w(t) 2 H r ≤ C w(t) L ∞ + ∇ h u h (t) L ∞ w(t)
M = max u h 0 L 2 (R 2 h ) , w 0 L 2 (R 2 h ) , w 0 L ∞ (R 2 h
) . Since C and M do not depend on p, we have

w(t) L ∞ ≤ CM e -γt .
We remark that x ln e + α

x , α > 0, is an increasing function. We deduce from ( Choosing

C 1 = CM (e + w 0 H r ) 2e CM γ , we get w(t) H r ≤ C 1 e -γt .
Combining the above estimate with (3.15) and using again the fact that x ln e + α x is an increasing function, we obtain the existence of a positive constant C 2 , depending on γ and w 0 H r , such that

∇ h u h (t) L ∞ ≤ C 2 e -γt .
In order to prove Theorem 1.2, we finally need the following lemma. Lemma 3.6. Let σ > 2. For any divergence-free vector fields u = (u 1 (x 1 , x 2 ), u 2 (x 1 , x 2 ), u 3 (x 1 , x 2 )) in H σ (R 2 h ), there exist a constant C > 0 independent of u, a summable sequence d q (u h ) q≥-1 of positive numbers such that

(3.22) ∆ h q (u h • ∇ h u h ) ∆ h q u h ≤ Cd q (u h )2 -2qσ ∇ h u h L ∞ u h 2
H σ , and a summable sequence (d q (u)) q≥-1 of positive numbers such that

(3.23) ∆ h q (u h • ∇ h u 3 ) ∆ h q u 3 ≤ Cd q (u)2 -2qσ u h H σ u 3 2 H σ .

Proof

We use the similar Bony decomposition as in the proof of Lemma 3.2, with w replaced by u h .

(3.24) ∆ h q (u h • ∇ h u h ) ∆ h q u h ≤ I 1 q (u h ) + I 2A q (u h ) + I 2B q (u h ) + I 2C q (u h
). Then, using Young's inequality, Bernstein lemma 2.1 and Proposition 2.3, we have

I 1 q (u h ) = ∆ h q q -q≥N S h q +2 (∇ h u h )∆ h q u h ∆ h q u h (3.25) ≤ q -q≥N S h q +2 (∇ h u h ) L ∞ ∆ h q u h L 2 ∆ h q u h L 2 ≤ C q -q≥N ∇ h u h L ∞ c q (u h )2 -q σ u h H σ c q (u h )2 -qσ u h H σ = C   c q (u h ) q -q≥N c q (u h )2 -(q -q)σ   2 -2qσ ∇ h u h L ∞ u h 2 H σ = Cd q (u h )2 -2qσ ∇ h u h L ∞ u h 2 H σ .
As in the proof of Lemma 3.2, integrations by parts imply (3.26)

I 2A q (u h ) = S h q u h • ∇ h ∆ h q u h ∆ h q u h = 0.
Since for any q, q ≥ 1, S h q -S h q -1 does not contain the low Fourier frequencies, I 2B q can be bounded in the same way as I 1 q , and we have

(3.27) I 2B q (u h ) = |q -q|≤N S h q -S h q -1 u h • ∇ h ∆ h q u h ∆ h q u h ≤ Cd q (u h )2 -2qσ ∇ h u h L ∞ u h 2 H σ .
The last term on the right hand side of (3.24) can be bounded exactly in the same way as in (3.14)

I 2C q (u h ) = |q -q|≤N ∆ h q , S h q -1 u h ∇ h ∆ h q u h ∆ h q u h (3.28) ≤ C2 -q |q -q|≤N S h q -1 ∇ h u h L ∞ ∇ h ∆ h q u h L 2 ∆ h q u h L 2 ≤ C2 -2qσ ∇ h u h L ∞ u h 2 H σ   c q (u h ) |q -q|≤N c q (u h )2 -(q -q)(σ-1)   ≤ Cd q (u h )2 -2qσ ∇ h u h L ∞ u h 2 H σ .
Combining Estimates (3.24) to (3.28), we obtain (3.22). Now, to prove (3.23), as in (3.24), with u h replaced by u 3 , we also write

(3.29) ∆ h q (u h • ∇ h u 3 ) ∆ h q u 3 ≤ I 1 q (u 3 ) + I 2A q (u 3 ) + I 2B q (u 3 ) + I 2C q (u 3 ).
We will estimate I 2A q (u 3 ) and I 2C q (u 3 ) exactly in the same way as I 2A q (u h ) and I 2C q (u h ). We have (3.30)

I 2A q (u 3 ) = S h q u h • ∇ h ∆ h q u 3 ∆ h q u 3 = 0,
and, using the inclusion

H σ-1 (R 2 h ) → L ∞ (R 2 h ) for σ > 2, I 2C q (u 3 ) = |q -q|≤N ∆ h q , S h q -1 u h ∇ h ∆ h q u 3 ∆ h q u 3 (3.31) ≤ C2 -q |q -q|≤N S h q -1 ∇ h u h L ∞ ∇ h ∆ h q u 3 L 2 ∆ h q u 3 L 2 ≤ C2 -2qσ ∇ h u h L ∞ u 3 2 H σ   c q (u 3 ) |q -q|≤N c q (u 3 )2 -(q -q)(σ-1)   ≤ Cd q (u 3 )2 -2qσ u h H σ u 3 2 H σ .
To estimate I 1 q (u 3 ), we write

I 1 q (u 3 ) = ∆ h q q -q≥N S h q +2 (∇ h u 3 )∆ h q u h ∆ h q u 3 (3.32) ≤ q -q≥N S h q +2 (∇ h u 3 ) L ∞ ∆ h q u h L 2 ∆ h q u 3 L 2 ≤ C q -q≥N ∇ h u 3 L ∞ c q (u h )2 -q σ u h H σ c q (u 3 )2 -qσ u 3 H σ = C   c q (u 3 ) q -q≥N c q (u h )2 -(q -q)σ   2 -2qσ ∇ h u 3 L ∞ u h H σ u 3 H σ = Cd q (u)2 -2qσ u h H σ u 3 2 H σ .
Finally, the term I 2B q (u 3 ) can be bounded in the same way as I 1 q (u 3 ), (3.33)

I 2B q (u 3 ) = |q -q|≤N S h q -S h q -1 u h • ∇ h ∆ h q u 3 ∆ h q u 3 ≤ Cd q (u)2 -2qσ u h H σ u 3 2 H σ .
Combining Estimates (3.29) to (3.33), we obtain (3.23).

Proof of Theorem 1.2

We recall that C to refer to a generic positive constant that may change from line to line. Since u h satisfies a two-dimensional damped Euler system, using Lemma 3.6, we deduce the following energy estimate in the Sobolev H σ -norm for u h 1 2

d dt u h (t) 2 H σ + γ u h (t) 2 H σ ≤ C ∇ h u h (t) L ∞ u h (t) 2 H σ . Since u h 0 ∈ H σ (R 2 h ), σ > 2, we have w 0 = ∂ 1 u 2 0 -∂ 2 u 1 0 ∈ H σ-1 , with σ -1 > 1.
Using Estimate (3.17) of Lemma 3.5, we get

d dt u h (t) H σ + γ u h (t) H σ ≤ CC 2 e -γt u h (t) H σ , which leads to d dt ln u h (t) H σ e γt ≤ CC 2 e -γt .
Integrating this inequality, we obtain

(3.34) u h (t) H σ ≤ u h 0 H σ + CC 2 γ e -γt .
Now, using Lemma 3.6, we get the following energy estimate in the Sobolev H σ -norm for u 3 1 2

d dt u 3 (t) 2 H σ + γ u 3 (t) 2 H σ ≤ C u h (t) H σ u 3 (t) 2 H σ .
Using (3.34), we can write

d dt u 3 (t) H σ + γ u 3 (t) H σ ≤ C u h 0 H σ + C 2 γ e -γt u 3 (t) H σ ,
which finally implies that

(3.35) u 3 (t) H σ ≤ C u 3 0 H σ + u h 0 H σ γ + C 2 γ 2 e -γt .
Theorem 1.2 is then proved.

Remark 3.7. If σ ∈ N, we do not need Lemma 3.6. Indeed, Theorem 1.2 can be immediately deduced from Lemma 3.5 using the identity

∇ α h ∇ h u h L 2 = ∇ α h w L 2 ,
which holds for all α ∈ N 2 .

Proof of the main result

In this paragraph, we provide the needed a priori estimates and a sketch of the proof of Theorem 1.1. These a priori estimates can be justified by a classical approximation by smooth fonctions (see for instance [START_REF] Chemin | Mathematical Geophysics: An introduction to rotating fluids and to the Navier-Stokes equations[END_REF]). We recall that it is already known that convergence in L ∞ loc R + , L 2 holds, as proved in [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF] and [START_REF] Chemin | Ekman boundary layers in rotating fluids, ESAIM Controle optimal et calcul des variations[END_REF]. The important point here is to prove global convergence in time.

4.1.

The corrective boundary layers. Let u be the solution of the limiting system (1.3), the third component u 3 = 0 and we set curl u h = ∂ 1 u 2 -∂ 2 u 1 . Then, as in [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF], we can write the boundary layer part in the following form (4.1)

B = B 1 + B 2 + B 3 + B 4 ,
where B i , i ∈ {1, 2, 3, 4}, are defined as described below:

1. The term B 1 is the main boundary layer term (see (A.7), (A.8) and (A.9) in the appendix)

(4.2) B 1 =     -e - x 3 ε √ 2β u 1 cos x3 ε √ 2β + u 2 sin x3 ε √ 2β -e - 1-x 3 ε √ 2β u 1 cos 1-x3 ε √ 2β + u 2 sin 1-x3 ε √ 2β -e - x 3 ε √ 2β u 2 cos x3 ε √ 2β -u 1 sin x3 ε √ 2β -e - 1-x 3 ε √ 2β u 2 cos 1-x3 ε √ 2β -u 1 sin 1-x3 ε √ 2β ε √ β G(x 3 ) curl u h    
where

G(x 3 ) = -e - x 3 ε √ 2β sin x 3 ε √ 2β + π 4 + e - 1-x 3 ε √ 2β sin 1 -x 3 ε √ 2β + π 4 .
2. The corrective terms B 2 and B 3 are added to B 1 to ensure the Dirichlet boundary conditions at {x 3 = 0} and {x 3 = 1}, and we have

B 2 =   ε √ 2β u 2 -ε √ 2β u 1 ε √ 2β 1 2 -x 3 curl u h   (4.3) B 3 = e -1 ε √ 2β cos 1 ε √ 2β   u 1 u 2 0   . (4.4)
3. Finally, the term B 4 is added to the sum B 1 + B 2 + B 3 to ensure the divergence-free property of B,

(4.5) B 4 = f (x 3 )   u 2 -u 1 0   + g(x 3 )   0 0 curl u h   , where f (x 3 ) = a e - x 3 ε √ 2β + e - 1-x 3 ε √ 2β + b, (4.6) g(x 3 ) = -ε β e -1 ε √ 2β sin 1 ε √ 2β + π 4 - x3 0 f (s)ds,
and where (a, b) is the solution of the linear system (4.7)

       1 + e -1 ε √ 2β a + b = -ε 2β + e -1 ε √ 2β sin 1 ε √ 2β 2ε 2β 1 -e -1 ε √ 2β a + b = 2ε βe -1 ε √ 2β sin 1 ε √ 2β + π 4 .
We remark that the determinant of the system (4.7) is

D = 1 + e -1 ε √ 2β -2ε 2β 1 -e -1 ε √ 2β
.

Thus, for ε > 0 small enough, we have D > 1 2 and (4.7) always has the following solution

(4.8) a = J ε -K ε D and b = K ε 1 + e -1 ε √ 2β -2J ε ε √ 2β 1 -e -1 ε √ 2β D ,
where

J ε = -ε 2β + e -1 ε √ 2β sin 1 ε √ 2β , K ε = 2ε βe -1 ε √ 2β sin 1 ε √ 2β + π 4 .
It is easy to prove that if ε > 0 is small enough, then

|a| < 4(β + β)ε and |b| < 32βε 2 .
With the previously defined boundary layer term B, we can verify that div (u + B) = 0 and (u + B)

| {x 3 =0} = (u + B) | {x 3 =1} = 0. Now, let B 0 (x 3 ) =     -e - x 3 ε √ 2β cos x3 ε √ 2β -e - 1-x 3 ε √ 2β cos 1-x3 ε √ 2β -e - x 3 ε √ 2β sin x3 ε √ 2β -e - 1-x 3 ε √ 2β sin 1-x3 ε √ 2β e - x 3 ε √ 2β sin x3 ε √ 2β + e - 1-x 3 ε √ 2β sin 1-x3 ε √ 2β -e - x 3 ε √ 2β cos x3 ε √ 2β -e - 1-x 3 ε √ 2β cos 1-x3 ε √ 2β    
Then, we can write B in the following form (4.9)

B = M(x 3 )   u 1 u 2 curl u h   .
The matrix M(x 3 ) is defined by

(4.10) M(x 3 ) = M (x 3 ) 0 0 m(x 3 )
where

M (x 3 ) = B 0 (x 3 ) + ε 2β + f (x 3 ) 0 1 -1 0 + e -1 ε √ 2β cos 1 ε √ 2β 1 0 0 1 , m(x 3 ) = ε β G(x 3 ) + ε 2β 1 2 -x 3 + g(x 3 ).
We can also prove the existence of a constant C > 0 such that, for any p ≥ 1, we have (

           M(•) L p x 3 ≤ Cε 1 p , M(•) L ∞ x 3 ≤ C, M (•) L p x 3 ≤ Cε 1 p -1 , m(•) L ∞ x 3 ≤ Cε, m(•) L p x 3 ≤ Cε sup x3∈[0, 1 2 ] 4.11)
x 2 3 M (x 3 ) ≤ Cε and sup

x3∈[ 1 2 ,1] (1 -x 3 ) 2 M (x 3 ) ≤ Cε.

4.2.

Convergence to the limiting system. It is easy to prove that B is small, i.e., B goes to 0 in

L ∞ R + , L 2 (R 2 h × [0, 1]
) as ε goes to 0. Then, our goal is to prove that

v ε = u ε -u -B converges to 0 in L ∞ R + , L 2 (R 2
h × [0, 1]) as ε goes to 0. Since a two-dimensional divergence-free vector field (independant of x 3 ) belongs to the kernel of the operator P(e 3 ∧ •), where P is the Leray projection of L 2 (R 3 ) onto the subspace of divergence-free vector fields, e 3 ∧ u is a gradient term. Replacing u ε by v ε + u + B in the first equation of the system (1.1), and using the fact that u is the solution of the system (1.3) and u 3 ≡ 0, we deduce that v ε satisfies the following equation

(4.12) ∂ t v ε -ν h (ε)∆ h v ε -βε∂ 2 3 v ε + u ε • ∇v ε + v ε • ∇B + v ε,h • ∇ h u + ∂ t B + B • ∇B + B h • ∇ h u + u h • ∇ h B -L 1 -L 2 + e 3 ∧ v ε ε = -∇ p ε ,
where

L 1 = ν h (ε)∆ h u + ν h (ε)∆ h B, L 2 = βε∂ 2 3 B - e 3 ∧ B ε + 2β u, ∇ p ε = ∇p ε -∇ h p + e 3 ∧ u ε .
Taking the L 2 scalar product of (4.12) with v ε , then integrating by parts the obtained equation and taking into account the fact that v ε satisfies the Dirichlet boundary condition, we get (4.13) 1 2

d dt v ε 2 L 2 + ν h (ε) ∇ h v ε 2 L 2 + βε ∂ 3 v ε 2 L 2 = -u ε • ∇v ε , v ε -v ε • ∇B, v ε -v ε,h • ∇ h u, v ε -∂ t B, v ε -B • ∇B, v ε -B h • ∇ h u, v ε -u h • ∇ h B, v ε + L 1 , v ε + L 2 , v ε .
In what follows, we will separately estimate the nine terms on the right-hand side of (4.13). We recall that we will always denote B 1 , B 2 and B 2 (v ε,1 , v ε,2 and v ε,3 respectively) the three components of B (v ε respectively) and we also write B = (B h , B 3 ) and v ε = (v ε,h , v ε,3 ).

1. The estimate of the first term is classical, using the divergence-free property of u ε :

(4.14) u ε • ∇v ε , v ε = 0.
2. The second term is the most difficult to treat. First, we decompose this term as follows (4.15)

v ε • ∇B, v ε = v ε,h • ∇ h B h , v ε,h + v ε,h • ∇ h B 3 , v ε,3 + v ε,3 ∂ 3 B 3 , v ε,3 + v ε,3 ∂ 3 B h , v ε,h .
For the first term on the right-hand side of (4.15), Hölder's inequality implies that

v ε,h • ∇ h B h , v ε,h ≤ C v ε,h L 2 ∇ h B h L ∞ v ε,h L 2 ≤ C M(•) L ∞ x 3 ∇ h u h L ∞ x h v ε 2 L 2 .
Then, using Estimates (4.11) and Lemma 3.5, we obtain

(4.16) v ε,h • ∇ h B h , v ε,h ≤ C(u 0 ) e -t √ 2β v ε 2 L 2 .
Next, using integrations by parts and Hölder's inequality, we deduce that

v ε,h • ∇ h B 3 , v ε,3 ≤ C ∇ h v ε L 2 B 3 L ∞ v ε L 2 .
So, Estimates (4.11), Lemmas 3.1 and 3.5 and Young's inequality imply

v ε,h • ∇ h B 3 , v ε,3 ≤ C m(•) L ∞ x 3 curl u h L ∞ x h ∇ h v ε L 2 v ε L 2 (4.17) ≤ C(u 0 ) ε e -t √ 2β v ε L 2 ∇ h v ε L 2 ≤ C(u 0 ) ε e -t √ 2β v ε 2 L 2 + C(u 0 ) ε 8 ∇ h v ε 2 L 2 .
Performing an integration by parts, we can control the third term on the right-hand side of (4.15) in the same way as the second one:

v ε,3 ∂ 3 B 3 , v ε,3 = 2 B 3 v ε,3 , ∂ 3 v ε,3 = 2 B 3 v ε,3 , div h v ε,h (4.18) ≤ C(u 0 ) ε e -t √ 2β v ε 2 L 2 + C(u 0 ) ε 8 ∇ h v ε 2 L 2 .
In order to estimate the last term of the right-hand side of (4.15), we decompose it into two parts, the first part corresponding to the boundary layer near {x 3 = 0} and the other part corresponding to the one near {x 3 = 1}:

v ε,3 ∂ 3 B h , v ε,h = R 2 h ×[0, 1 2 ] (v ε,3 ∂ 3 B h ) • v ε,h dx + R 2 h ×[ 1 2 ,1] (v ε,3 ∂ 3 B h ) • v ε,h dx
For the first part, since v ε vanishes on {x 3 = 0}, using Hölder's inequality and Hardy's inequality, we get

R 2 h ×[0, 1 2 ] (v ε,3 ∂ 3 B h ) • v ε,h dx ≤ sup x3∈[0, 1 2 ] x 2 3 M (x 3 ) u L ∞ x h v ε,3 x 3 L 2 v ε,h x 3 L 2 ≤ sup x3∈[0, 1 2 ] x 2 3 M (x 3 ) u L ∞ x h ∂ 3 v ε,3 L 2 ∂ 3 v ε,h L 2 .
We recall that ∂ 3 v ε,3 = -div h v ε,h . Then, Lemmas 3.1 and 3.5, Estimates (4.11) and Young's inequality imply

R 2 h ×[0, 1 2 ] (v ε,3 ∂ 3 B h ) • v ε,h dx ≤ C(u 0 ) ε e -t √ 2β div h v ε,3 L 2 ∂ 3 v ε,h L 2 (4.19) ≤ C(u 0 ) ε 8 ∇ h v ε 2 L 2 + βε 4 ∂ 3 v ε 2 L 2 .
For the second part concerning the boundary layer near {x 3 = 1}, since v ε = (v ε,h , v ε,3 ) vanishes on {x 3 = 1}, Hardy's inequality implies that

I v = R 2 h 1 1 2 v ε,3 (x h , x 3 ) 1 -x 3 2 dx 3 dx h = R 2 h 1 2 0 v ε,3 (x h , 1 -x 3 ) x 3 2 dx 3 dx h ≤ C R 2 h 1 2 0 ∂ 3 v ε,3 (x h , 1 -x 3 ) 2 dx 3 dx h ≤ C ∂ 3 v ε,3 2 L 2 = C div h v ε,h 2 L 2 . Likewise, I h = R 2 h 1 1 2 v ε,h (x h , x 3 ) 1 -x 3 2 dx 3 dx h ≤ C ∂ 3 v ε,h 2 L 2 .
Thus, using Hölder's inequality, we get

R 2 h ×[ 1 2 ,1] (v ε,3 ∂ 3 B h ) • v ε,h dx ≤ sup x3∈[ 1 2 ,1] (1 -x 3 ) 2 M (x 3 ) u L ∞ x h I v I h (4.20) ≤ C(u 0 ) ε e -t √ 2β div h v ε,3 L 2 ∂ 3 v ε,h L 2 ≤ C(u 0 ) ε 8 ∇ h v ε 2 L 2 + βε 4 ∂ 3 v ε 2 L 2 .
3. The third term on the right-hand side of (4.13) can be treated using Hölder's inequality and Lemma 3.5. We have

(4.21) v ε,h • ∇ h u, v ε ≤ C ∇ h u h L ∞ (R 2 h ) v ε 2 L 2 ≤ C(u 0 ) e -t √ 2β v ε 2 L 2 .
4. We recall that

B = M(x 3 ) t A = M (x 3 ) 0 0 m(x 3 )   u 1 u 2 curl u h   .
Then, we can write the fourth term as

∂ t B, v ε = M (x 3 )∂ t u h , v ε,h + m(x 3 )∂ t curl u h , v ε,3 . Let P R 2
h be the Leray projection of L 2 (R 2 h ) onto the subspace of two dimensional divergence-free vector fields. Since u 3 = 0, and ∂ 3 p = 0, using Estimate (4.11) and Lemmas 3.1 and 3.5, we have

M (x 3 )∂ t u h , v ε,h ≤ M (x 3 )P R 2 h u h • ∇ h u h , v ε,h + 2β M (x 3 )u h , v ε,h (4.22) ≤ M (•) L 2 x 3 P R 2 h u h • ∇ h u h L 2 + 2β u h L 2 v ε,h L 2 ≤ Cε 1 2 u h L 2 1 + ∇ h u h L ∞ v ε L 2 ≤ C(u 0 ) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 .
Applying the operator curl to the first two equations of the system (1.3) (we recall that in this paper, curl only acts on the horizontal components and we already defined curl

u h = ∂ 1 u 2 -∂ 2 u 1 ), we obtain ∂ t (curl u h ) + u h • ∇ h (curl u h ) + 2β (curl u h ) = 0. Then, m(x 3 )∂ t curl u h , v ε,3 ≤ m(x 3 )u h • ∇ h curl u h , v ε,3 + 2β m(x 3 )curl u h , v ε,3 (4.23) ≤ m(•) L 2 x 3 u h • ∇ h curl u h L 2 + 2β curl u h L 2 v ε,3 L 2 ≤ Cε u h L ∞ ∇ h curl u h L 2 + curl u h L 2 v ε L 2 ≤ C(u 0 ) ε e -t √ 2β 1 + v ε 2 L 2 .
5. We decompose the fifth term into two parts:

B • ∇B, v ε = B h • ∇ h B, v ε + B 3 ∂ 3 B, v ε .
The "horizontal" part can be bounded using Hölder's inequality, Estimates (4.11) and Theorem 1.2 as follows

B h • ∇ h B, v ε ≤ C B h L ∞ ∇ h B L 2 v ε L 2 (4.24) ≤ C M (•) L ∞ x 3 u h L ∞ M(•) L 2 x 3 ∇ h u h L 2 + ∇ h curl u h L 2 v ε L 2 ≤ Cε 1 2 u h 2 H σ v ε L 2 ≤ C(u 0 ) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 .
Likewise, we have the following estimate for the "vertical" part:

B 3 ∂ 3 B, v ε ≤ B 3 L ∞ ∂ 3 B L 2 v ε L 2 (4.25) ≤ C(u 0 ) e -t √ 2β m(•) L ∞ x 3 M (•) L 2 x 3 v ε L 2 ≤ C(u 0 ) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 .
6. For the sixth term, taking into account the fact that u is independent of x 3 , Estimates (4.11) and Lemma 3.1 imply

B h • ∇ h u, v ε ≤ B h L 4 x h L 2 x 3 ∇ h u L 4 x h v ε L 2 (4.26) ≤ M (•) L 2 x 3 u 2 L 4 + ∇ h u 2 L 4 v ε L 2 ≤ C(u 0 ) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 .
7. The seventh term can be bounded in a similar way as the sixth one.

u h • ∇ h B, v ε ≤ u h L ∞ ∇ h B L 2 v ε L 2 (4.27) ≤ C u h L ∞ M(•) L 2 x 3 ∇ h u h L 2 + ∇ h curl u h L 2 v ε L 2 ≤ C(u 0 ) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 .
8. We write the eighth term as

L 1 , v ε = ν h (ε) ∆ h u h , v ε,h + ν h (ε) ∆ h B, v ε .
Then, we have

ν h (ε) ∆ h u h , v ε,h ≤ ν h (ε) u h H σ v ε L 2 (4.28) ≤ C(u 0 ) ν h (ε) e -t √ 2β 1 + v ε 2 L 2 and |ν h (ε) ∆ h B, v ε | = ν h (ε) | ∇ h B, ∇ h v ε | ≤ ν h (ε) ∇ h B L 2 ∇ h v ε L 2 (4.29) ≤ ν h (ε) M(•) L 2 x 3 ∇ h u h L 2 + ∇ h curl u h L 2 v ε L 2 ≤ C(u 0 ) ν h (ε) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 .
9. We will estimate the ninth term as in [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF] or [START_REF] Masmoudi | Ekman layers of rotating fluids: the case of general initial data[END_REF]. We have

L 2 , v ε = βε∂ 2 3 B, v ε - e 3 ∧ B ε , v ε + 2β u, v ε .
We recall the definition (4.1) of B B = B 1 + B 2 + B 3 + B 4 , and for any i ∈ {1, 2, 3, 4}, we set B i = (B h i , B 3 i ), where B h i and B 3 i denote the horizontal and vertical components of B i respectively. We also recall that B = (B h , B 3 ). Then, the following identities are immediate ∂ 2 3 B h 3 = 0,

βε∂ 2 3 B h 1 - e 3 ∧ B 1 ε = 0, βε∂ 2 3 B h 2 - e 3 ∧ B 2 ε + 2β u = 0.
For the remaining terms, we have

βε ∂ 2 3 B 3 , v ε,3 ≤ βε ∂ 3 B 3 L 2 ∂ 3 v ε,3 (4.30) ≤ βε ∇ h B h L 2 ∇ h v ε L 2 ≤ βε M (x 3 ) L 2 x 3 ∇ h u h L 2 x h ∇ h v ε L 2 ≤ C(u 0 ) ε 2 e -t √ 2β + C(u 0 ) ε 8 ∇ h v ε 2 L 2 ; e 3 ∧ B 3 ε , v ε,h ≤ Cε -1 e -1 ε u L 2 v ε L 2 ≤ C(u 0 ) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 . (4.31)
We recall that

B h 4 = f (x 3 ) u 2 -u 1 ,
where f is defined in (4.6),

f (x 3 ) = a exp - x 3 ε √ 2β + exp - 1 -x 3 ε √ 2β + b,
and that, if ε > 0 is small enough, we have |a| < 4(β + β)ε and |b| < 32βε 2 .

Then, End of the proof of Theorem 1.1: Summing all the inequalities from (4.14) to (4.33), we deduce from (4.13) that, for ε > 0 small enough

βε ∂ 2 3 B h 4 , v ε,h ≤ βε f (•) L 2 x 3 u h L 2 v ε,h L 2 (4.32) ≤ Cε 1 2 u h L 2 v ε L 2 ≤ C(u 0 ) ε 1 2 e -t √ 2β 1 + v ε 2 L 2 .
d dt v ε 2 L 2 + ν h (ε) ∇ h v ε 2 L 2 + βε ∂ 3 v ε 2 L 2 ≤ C(u 0 ) ε 1 2 + ν h (ε) e -t √ 2β + C(u 0 ) e -t √ 2β v ε 2 L 2 + C(u 0 ) ε ∇ h v ε 2 L 2 .
If C(u 0 ) ε ≤ ν h (ε), for any ε > 0 small enough, by integrating the above inequality with respect to the time variable, we get

(4.34) v ε (t) 2 L 2 ≤ v ε (0) 2 L 2 + C(u 0 ) ε 1 2 + ν h (ε) √ 2β + C(u 0 ) t 0 e -s √ 2β v ε (s) 2 L 2 ds.
We recall that v ε = u ε -u -B. Thus,

v ε (0) 2 L 2 ≤ u ε (0) -u(0) 2 L 2 + B(0) 2 L 2 ≤ u ε 0 -u 0 2 L 2 + M(•) 2 L 2 x 3 u 0 2 L 2 ≤ u ε 0 -u 0 2 L 2 + Cε 1 2 u 0 2 L 2 .
According to the Gronwall lemma, it follows from (4.34) that

v ε (t) 2 L 2 ≤ u ε 0 -u 0 2 L 2 + Cε

  H r . From (3.10)-(3.14), we easily deduce (3.7) and Lemma 3.2 is proved. Lemma 3.4. Let r > 1.

Finally, we have e 3

 3

1 2 u 0 2 L 2 + 1 2

 221 C(u 0 ) ε + ν h (ε) √ 2β exp C(u 0 ) √ 2β , which leads to lim ε→0 v ε L ∞ (R+,L 2 )= 0, and Theorem 1.1 is proved.

E-mail address: van-sang.ngo@univ-rouen.fr

Appendix A. Contruction of the Ekman layers

The construction of Ekman boundary layers is classical and can be found in [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF] or [START_REF] Chemin | Mathematical Geophysics: An introduction to rotating fluids and to the Navier-Stokes equations[END_REF] for instance. In this paper, we will only give very brief recall of this construction. The typical approach consists in looking for the appropriate solutions of the system (1.1) in the following form (A.1) u ε = u I,0 + u B,0 + ũB,0 + ε u I,1 + u B,1 + ũB,1 + ε 2 u I,2 + u B,2 + ũB,2 + . . . p ε = 1 ε p I,-1 + p B,-1 + pB,-1 + p I,0 + p B,0 + pB,0 + ε p I,1 + p B,1 + pB,1 + . . . .

The "interior part" functions with the index "I" depend on (t, x 1 , x 2 , x 3 ) and the "boundary layer part" with the index "B" consists in smooth functions of the form

and which rapidly (exponentially) decrease when the third space variable goes to infinity. The justification of the size of the boundary layers being ε and the discussion about other sizes of the layers can be found in [START_REF] Grenier | Ekman layers of rotating fluid, the case of well prepared initial data[END_REF], [START_REF] Masmoudi | Ekman layers of rotating fluids: the case of general initial data[END_REF], [START_REF] Chemin | Ekman boundary layers in rotating fluids, ESAIM Controle optimal et calcul des variations[END_REF] and [START_REF] Chemin | Mathematical Geophysics: An introduction to rotating fluids and to the Navier-Stokes equations[END_REF].

A.1. In the interior part of the domain. Away from the boundary, all the boundary terms in (A.1) rapidly decrease to zero. Thus, at the leading order ε -1 , the divergence-free conditions and the fast rotation yield

The limiting velocity u I,0 is then a two-dimensional divergence-free vector field with three components, which means that the fluid has tendency to move in columns when the rotation is fast, as predicts the Taylor-Proudman theorem. Now, looking at the zeroth order (ε 0 ), one obtains the governing equation of the fluid at the limit as ε → 0

Here, we remark that the fast rotation implies the action of the first order term, i.e. the added vector field -u I,1 2 , u I,1 1 , 0 , on the limiting behavior of the fluid. In the next paragraph, we will show that the interaction with the boundary allows to explicitly calculate this vector field, which turns out to be a dissipative term (the Ekman pumping).

A.2. In the Ekman boundary layers. We will focus on the boundary layer near {x 3 = 0}. The other layer near {x 3 = 1} can be obtained in a similar way. Performing the change of variable y =

x 3 ε , we deduce the following "divergence-free properties" of u B,j , for any j ≥ 0,

We remark that the first equation of (A.4) implies that u B,0

3

≡ 0 because it goes to zero as y → +∞.

The Dirichlet boundary condition u ε (t, x 1 , x 2 , 0) = 0 is rewritten as follows

which implies u I,0 3 ≡ 0, since u I,0 3 is independent of x 3 and u I,0 (t, x 1 , x 2 , 0) = -u B,0 (t, x 1 , x 2 , 0) = 0. Now, putting the Ansatz into the first equation of (1.1) and looking at the leading order ε -2 , we simply get ∂ y p B,-1 = 0, and so, p B,-1 = 0 since p B,-1 → 0 as y → +∞. This is a classical principle of fluid mechanics which claims that the pressure does not vary in a boundary layer. At the next order ε -1 , combining the boundary condition (A.5) and the fact that u I,0 = u I,0 (t, x 1 , x 2 ) is independent of x 3 , we recover the following system

This system can be explicitly solved and that gives (A.7)

Using (A.4), (A.7) and the fact that ∂ 1 u I,0 1 + ∂ 2 u I,0 2 = 0, we have

Integrating the above equation with respect to y and recalling that u B,1 3 → 0 as y → +∞, we obtain

Remark A.1. We remark that unlike the irrotational case where the boundary layers are discribed by the Prandtl equations, in this case, the main boundary layer term u B,0 1 , u B,0 2 , u B,1 3 can be explicitly computed from the limiting velocity in the interior domain, which allows to directly treat these layers without using the Prandtl equations.

The same arguments near the boundary {x 3 = 1} lead to pB,-1 = 0 and (A.9)

Since u I,0 is independent of x 3 , by integrating (A.10) with respect to x 3 over [0, 1], one gets


Saeid Nezamabadi
email: saeid.nezamabadi@umontpellier.fr

Xavier Frank

Jean-Yves Delenne

Julien Averseng

Farhang Radjai

Parallel implicit contact algorithm for soft particle systems

Keywords:

This paper presents a numerical technique to model soft particle materials in which the particles can undergo large deformations. It combines an implicit finite strain formalism of the Material Point Method and the Contact Dynamics method. In this framework, the large deformations of individual particles as well as their collective interactions are treated consistently. In order to reduce the computational cost, this method is parallelised using the Message Passing Interface (MPI) strategy. Using this approach, we investigate the uniaxial compaction of 2D packings composed of particles governed by a Neo-Hookean material behaviour. We consider compressibility rates ranging from fully compressible to incompressible particles. The packing deformation mechanism is a combination of both particle rearrangements and large deformations, and leads to high packing fractions beyond the jamming state. We show that the packing strength declines when the particle compressibility decreases, and the packing can deform considerably. We also discuss the evolution of the connectivity of the particles and particle deformation distributions in the packing.

Introduction

The macroscopic behaviour of particulate materials is controlled by the microscopic mechanisms in terms of the interactions between individual particles as well as interactions with a surrounding fluid or confining walls. Understanding these mechanisms can be effectively achieved via particle scale simulation techniques based on microdynamic information. The Discrete Element Method (DEM) [1,2] and Contact Dynamics (CD) method [3][START_REF] Jean | Proceeding of The International Deep Drawing Research Group[END_REF][START_REF] Radjai | [END_REF] are recognised as efficient research tools for the investigation of the micromechanics of particulate materials. These methods are capable of dealing with different loading conditions, particle size distributions and physical properties of the particles. Such discrete simulations can provide detailed local information such as the trajectories of individual particles and transient forces acting on them that can be difficult to obtain by physical experimentation.

In the context of DEM methods, the particles are assumed to be hard or weakly deformable through different contact theories such as the Hertz contact theory, which is only valid up to about 10% of strain. However, this assumption is too crude in the application to highly soft particles such as metallic powders, many pharmaceutical and food products, and colloidal suspensions [6][7][8][9][10]. Soft particles may undergo large deformations without rupture. * Corresponding author.

Hence, as the classical DEM techniques are intrinsically unable to account for realistic constitutive models for individual particles and large particle deformations, soft particle materials require a methodology capable of treating the contact interactions between particles as well as individual particle deformations. We previously proposed a numerical procedure based on an implicit material point method (MPM) coupled with the CD method [11,12]. In the MPM, each particle is discretised by a set of material points carrying all state variables such as stress and velocity field. The MPM algorithm also uses a background grid for solving the momentum equations. The material points are assigned fixed masses during computation so that the conservation of mass is satisfied implicitly. The momentum changes are interpolated from the grid to the material points so that the total momentum is conserved. The implicit formulation allows for efficient coupling with implicit modelling of unilateral contacts and friction between the particles as in the CD method [3,[START_REF] Acary | Numerical Methods for Nonsmooth Dynamics[END_REF].

In the present paper, we propose a parallel implicit MPM procedure for the simulation of deformable particles in the context of the finite strain theory as an extension of our previous model based on the infinitesimal strain hypothesis [11,12]. This novel formulation allows for applying a large class of material behaviours like hyperelasticity [START_REF] Nezamabadi | [END_REF]. Furthermore, a parallel algorithm based on MPI (Message Passing Interface) is proposed in the context of the MPM. It permits to improve considerably the computational performance of our MPM framework. We apply this method to study the compaction of a packing of soft particles. The soft-particle packings may undergo volume change as a consequence of particle rearrangements as in hard-particle materials. But, their property of volume change by particle shape and size change under moderate external loads, leads to enhanced space filling. It allows the packing fraction to exceed the random close packing (RCP) limit [15][16][17]. The compaction and other rheological properties of soft-particle systems beyond this 'jamming' point are still poorly understood. Our results show the capability of the MPM coupled with CD for the investigation of soft particle packings beyond the RCP limit. We focus on the evolution of the packing and effects of particle shape change. As we shall see, the particle material behaviour affects the stress level and its evolution during compaction.

The paper is organised as follows. In Section 2, the new MPM formulation based on the finite strain theory and our contact algorithm are introduced. Section 3 is devoted to the presentation of the implicit MPM resolution. Then, in Section 4 we describe the parallelisation procedure of our MPM-CD method. In Section 5, we focus first on the behaviour of a single particle subjected to axial strain. Then, we analyse the compaction process of a packing of soft circular particles. We conclude with a brief summary and perspectives of this work.

Material point method formulation

In this section, we describe the basic formulation of the material point method in the context of finite strain theory. Similar formulations have been presented in our previous papers [11,12] in which the infinitesimal strain theory has been considered, where for modelling soft particles, the MPM has been coupled with the CD method for the treatment of frictional contacts between particles.

Let Ω t be a domain in R D , D being the domain dimension, associated with a continuum body, in its actual configuration at time t. Its conservation of mass is described by this continuity equation:

∂ρ( t x, t) ∂t + ∇t x ( ρ( t x, t) • v( t x, t) ) = 0 in Ω t , (1)
where ρ( t x, t) indicates the material density and v( t x, t) denotes the velocity field at position t x (the prefix superscript 't' indicates the time) in the actual configuration Ω t at time t. The conservation of linear momentum for this continuum body is defined by

∇t x • σ( t x, t) + b( t x, t) = ρ( t x, t) a( t x, t) in Ω t , (2)
where σ(

{ u( t x, t) = û(t) on ∂Ω u t , σ( t x, t) • n = f(t) on ∂Ω f t , (3)
where u( t x, t) and û(t) are the displacement field and the prescribed displacement, respectively. Here, n denotes the outward unit normal vector to ∂Ω t and f(t) is a prescribed load.

In the MPM, the continuum body domain is divided into N p infinitesimal constant mass elements called material points. Because of this assumption (constant material point mass), the mass conservation relation (1) is self-satisfied. Furthermore, the material points serve as integration points to compute the FEM integrals. The MPM then discretises these integrals through a Dirac delta function by considering a fixed material point mass. Hence, the weak form of the equation of motion (2) in its discretised version can be written as follows by considering the contact interactions between several bodies [11]:

M a node (t) = f int (t) + f ext (t) + f c (t) , (4)
where a node is the nodal acceleration, f c denotes the contact force, which will be illustrated below, and

M = Np ∑ p=1 m p N p lumped mass matrix, f int (t) = - Np ∑ p=1 G p σ p (t) V p (t) internal force vector, f ext (t) = Np ∑ p=1 N p b p (t) + f s (t)
sum of body forces and surface tractions f s .

In the above relations, V p denotes the material point volume and N p is the interpolation matrix or the shape function matrix at a material point p. It relates the quantities associated with the material points (displacement, position• • •) to nodal variables of the element to which the material point belongs. G p denotes the gradient of the shape function N p . Since there are generally more material points than grid nodes, a weighted squares approach is used to determine nodal velocities v node from the material point velocities v p . Hence, the nodal velocities are obtained by solving the relation

P node (t) = M v node (t) = ∑ Np p=1 m p N p v p (t) , (5)
where P node is the nodal momentum.

It is also important to note that, as we deal with deformable particle systems, the contact forces f c between particles need to be computed using a contact algorithm that accounts for the condition of impenetrability of matter as well as the Coulomb friction law. This contact algorithm combines the MPM and CD methods that was presented in detail in our previous paper [11]. For clarity, in the following, we briefly describe this algorithm.

Let us consider two deformable particles (α and β); see Fig. 1.

In the context of the multi-mesh algorithm, a proper background mesh is attributed to each particle. A contact point at the interface between the two particles may be treated by introducing a common background mesh with the same type of grids for the transfer of nodal quantities from the proper meshes to the common mesh.

The contact points between the particles α and β are treated at the neighbouring nodes belonging to the common background mesh. Their nodal values involve contributions from the two particles. At a potential contact node i, a normal unit vector n i , oriented from particle β to particle α, and a tangential unit vector t i are defined [18]. As long as the normal velocity remains positive, the normal force f n is identically zero. But when v n = 0, a non-negative (repulsive) normal force f n is mobilised at the contact node. These conditions define the velocity-Signorini complementary condition as shown in Fig. 2(a) [START_REF] Jean | Mechanics of Geomaterial Interfaces[END_REF][START_REF] Brogliato | Nonsmooth Mechanics[END_REF]. On the other hand, by combining the equations of motion

v n (v n = (v α i -v β i ) • n i )
P α node = M α v α
f n = 1 ∆t m α i m β i m α i +m β i v n + k n , (6)
where m α i and m β i are the nodal masses of bodies of α and β, respectively, ∆t denotes the incremental time, and k n is an offset force which depends on other contact forces exerted by the neighbouring bodies of α and β. The normal force at all contact nodes are obtained through an iterative process by intersecting the above linear relation with the Signorini graph, as shown in Fig. 2(a).

In a similar vein, the Coulomb law of dry friction is a complementarity relation between the friction force f t and the tangential

velocity v t (v t = (v α i -v β i ) • t i ) at the contact node; see Fig. 2(b).
Like the Signorini graph, the Coulomb law is a complementarity relation in the sense that it cannot be reduced to a single-valued function. The equations of motion at the common node i yield

f t = 1 ∆t m α i m β i m α i +m β i v t + k t , (7)
which is intersected with the Coulomb graph to calculate the friction force f t simultaneously at all contact nodes in the same iterative process used to calculate the normal forces. The convergence to the solution both for contact forces and internal stresses is smooth, and a high precision may be achieved through the convergence criterion.

It is worth noting that in the presented algorithm, a contact may occur between the particles even if they are not physically in contact. Indeed, since the contact is computed on the nodes of the background mesh (not on the material points), the distance between the particles in contact can vary within one element size. The contact force accuracy depends hence on the particle discretisation as well as time discretisation (time step). This issue exists in all contact algorithms and depending on the necessary solution accuracy required for a specified problem, one can adjust the time and/or space resolution. In our case, the proposed contact algorithm allows us to treat rapidly and accurately enough the contact between deformable particles with any arbitrary shape although some local parameters such as contact surface may not be accurately defined.

A finite strain formulation for MPM

To complement the continuity equation ( 1) and the momentum equation ( 2), we consider a constitutive relationship in the context of the finite strain theory:

t 0 Π( 0 x, t) = F (r) ( t 0 F( 0 x, t)) , (8)
where t 0 Π( 0 x, t) is the first Piola-Kirchhoff stress tensor at position 0 x in the initial configuration and at time t. Let t 0 F( 0 x, t) = ∇0 x u( t x, t) + I be the deformation gradient tensor, where I is the second-order identity tensor. Note that t 0 Π and t 0 F are defined at the actual configuration 't' with respect to the initial configuration

at time t = 0. t σ is also related to t 0 Π through t σ( t x, t) = 1 0 J t 0 Π( 0 x, t) ( t 0 F( 0 x, t)) T , (9)
with 0 J = det( t 0 F( 0 x, t) ). Note that, by virtue of the definition of the constitutive relation (8), this framework corresponds to the finite strain theory, and it is not specifically designed for a particular constitutive law. Hence, the material behaviour can cover various nonlinear and complex physical and geometrical evolutions of the continuum body.

Finite strain MPM: an implicit-type formalism

In our previous paper [11], a MPM algorithm with an implicit time integration was introduced. In this section, we adopt this approach to our new formulation in the framework of the finite strain theory. Note that the implicit resolution concerns only the nodal parameters whereas those related to the material points are determined explicitly.

Let us advance the solution of (4) from 't' to 't + ∆t' in the context of the implicit resolution. We consider that f ext (t + ∆t) is known, and the grid kinematics is advanced in time as follows:

u node (t + ∆t) = ∆t v node (t + ∆t) , (10)
v node (t + ∆t) = v node (t) + ∆t a node (t + ∆t) . (11)
Note that, in Eq. ( 10), we have u node (t) = 0 since u node (t + ∆t) is in fact the grid displacement from 't' to 't + ∆t'. From Eqs. (10) and (11), the nodal acceleration at time t + ∆t is given by

a node (t + ∆t) = 1 ∆t 2 u node (t + ∆t) - 1 ∆t v node (t) . (12)
In the context of the finite strain theory, the evaluation of the material point volume V p changes from t to t + ∆t according to

V p (t + ∆t) = t J p V p (t) (13)
with t J = det( t+∆t t F ). In an incremental-iterative resolution algorithm, a new estimation of the nodal displacement u k node (t + ∆t) at iteration k is obtained by adding the incremental displacement ∆u k node to the previous estimated displacement:

u k node (t + ∆t) = u k-1 node (t + ∆t) + ∆u k node . (14)
To obtain ∆u k node at iteration k, we solve

K k-1 ∆u k node = R k , ( 15
)
where K is the stiffness matrix and R refers to the residual term. This equation is the incremental form of relation [START_REF] Jean | Proceeding of The International Deep Drawing Research Group[END_REF]. The terms K and R are defined in Appendix A.

This incremental algorithm finds a nodal displacement u node (t + ∆t) that minimises the residual term, R. So, as in [START_REF] Guilkey | [END_REF], we introduce two convergence criteria:

C 1 = ∥∆u k node ∥ ∥∆u max node ∥ < ϵ 1 and C 2 = ∥∆u k node R k ∥ ∥∆u 1 node R 1 ∥ < ϵ 2 , ( 16
)
where ϵ 1 and ϵ 2 are tolerance parameters on velocities and energy, respectively, ∥ • ∥ is the norm operator, ∥∆u max node ∥ denotes the maximum value of the norm of the incremental displacement, and

∥∆u 1 node R 1
∥ indicates the initial value of the inner product of the incremental displacement and residual term.

MPM and parallel computation

The MPM simulations involving a modest number of particles and material points can be performed in a reasonable time on a single-processor workstation. As the number of particles and physical complexity of the numerical model increase, so does the computational resources required. The simulation of 1000 particles, for example, with a sufficient number of material points for their discretisation is not thus suited to a single processor. Hopefully, the required level of computational power can be obtained by parallel MPM simulations. Herein, the goal is to produce a portable parallel implementation that would exhibit good performance scaling up to several hundred processors for large-scale simulations. The code is required to run on parallel computer systems having either shared or distributed memory.

Parallelisation procedure

The proposed MPM parallel algorithm is of SPMD (Single Program, Multiple Data) style, based on MPI. The particles in the simulation domain are divided into spatial sub-domains. Choosing the number of sub-domains equal to the number of available processors, the data associated with particles in the same sub-domain are stored in the memory of a single processor. The computational effort depends on the resolution of Eq. (4) for each particle, which is closely related to the particle volume. Assuming that the processors have equal performances, to achieve load balance (needed to minimise synchronisation delays), the sub-domains should be chosen such that the sum of volumes of particles is nearly the same in all sub-domains. Hence, the particles in a sub-domain are not necessarily neighbours.

In the proposed MPM algorithm, the resolution of Eq. ( 4) imposes the computation of the contact forces f c for each particle, which depend on the interactions between particles. Since the contact forces are computed on the master process, the only data exchange between the particles (and thus the processors) are the contact forces f c ; see Fig. 3. It is also worth noting that the computation of contact forces on the master process is constant independently of the number of processors and represents less than about 0.5% of the total computational time.

Load balance

Each MPI thread manages MPM computations of a set of particles, one particle being attributed to only one process. So, the elementary MPM computational weight unit is the weight of a particle. However, particles involved in the global simulation do not exhibit exactly the same volume, and thus the same computational weight. To reach correct performances, the particles have to be distributed in such a way that the weight is well balanced between processes. To do so, we propose the following algorithm for P processes:

1. The N particles are sorted by decreasing radius. 2. The P largest sorted particles are attributed to the P processes.

3. The P + 1 particle is attributed to the process with the lightest workload 4. Step 3 is iterated until all particles are attributed. Fig. 4 displays our load balance algorithm for an example of circular particles in 2D. As the weight here is proportional to the particle surface (in 2D), the total surface of particles ∑ R i 2 must be well balanced. This approach allows reaching very low load imbalances. For instance, it is less than 1% for N = 300 particles (radii of particles ranging from 0.7 mm to 1.2 mm) and P = 60 processes.

Scalability

The scalability of the code was studied with the help of the cluster of the Genotoul Bioinformatics Platform (Toulouse, France). Each compute node embeds 2 Ivy-Bridge 10 cores hyper-threaded microprocessors (2.5 GHz) and the nodes are interconnected through a QDR Infini-Band network for both MPI communications and IO. In our study, the number of processors P varies from 1 to 60, and according to this number at most 20 cores were used per compute node.

As shown in Fig. 5, the scaling is not linear and the efficiency decreases quite fast when P increases. Clearly, such a behaviour can be attributed to the communication bottleneck in rank 0 process. Collective communications involved in data exchange between MPM and Contact Dynamics lead to a large amount of data to be sent or received by rank 0 MPI thread. Such a point should be addressed in future improvements of our code.

Numerical examples

The accuracy and efficiency of the proposed algorithm within the finite strain theory are studied through several mechanical compaction tests. In our previous works [11,12], the performance of a similar approach in the framework of the infinitesimal strain hypothesis was shown. We propose two main applications. The first one deals with the uni-axial deformation of a single soft particle. The second example concerns the compaction of a packing of 300 soft particles. To avoid stress gradients in these examples, the gravitational acceleration is set to be zero.

In the MPM, two-dimensional simulations in plane strain conditions were performed. The computation domain was meshed with four-node quadrangular elements. For the applications that we target in the present work, we consider two types of material behaviours [START_REF] Nezamabadi | [END_REF]: a linear Saint-Venant Kirchhoff constitutive relation S = λ Tr(γ) I + µ γ , (17) and a nonlinear Neo-Hookean constitutive law

S = (λ ln(J) -µ) C -1 + µ I , ( 18
)
where S is the second Piola-Kirchhoff stress tensor and related to Π through Π = F S. Let C = F T F be the right symmetric Cauchy-Green tensor. γ denotes the Green-Lagrange strain tensor (γ = 1 2 (C -I)) and J = det (F). λ and µ represent the Lamé coefficients. Fig. 6. Contact geometry between a single particle and a rigid plate.

Axial compaction of a single particle

We consider here the case of a single cylindrical particle subjected to axial compression. The particle has a diameter of D = 20 mm and is compressed between two rigid walls as shown in Fig. 6. The bottom wall is fixed and the top wall moves downwards at a constant velocity of 0.5 m/s. The time step is set to ∆t = Fig. 7. Normal contact force applied on a single particle as a function of the displacement of the particle centre for two different behaviours.

0.1 µs. To compare the infinitesimal and finite strain formulations in the context of the MPM, the linear Hookean and Saint-Venant Kirchhoff elastic behaviours (see Eq. ( 17)) are considered. In the two cases, the Lamé coefficients and density of the particles were set to λ = 100 MPa, µ = 1.5 MPa and ρ = 990 kg/m 3 , respectively. Fig. 7 presents the normal contact force F as a function of displacement d of the centre of the particle. In the two cases, a quasi-linear evolution of force with displacement is observed but with a small deviation for d D > 0.05 towards a lower level of force.

This behaviour corresponds to the prediction of the Hertz analysis for a cylinder of unit length [START_REF] Johnson | Contact Mechanics[END_REF]:

F = π 4 E * d , (19)
where E * is the effective elastic modulus defined as E * = E/(1-ν 2 ) with E being Young's modulus and ν Poisson's ratio. Furthermore, the predicted values of forces by the infinitesimal and finite strain formulations are not very different. This means that, as a result of the small value of the time step, the second order terms in the Green-Lagrange strain rate have little effect on the total strain.

We carried out the same test by considering a Neo-Hookean particle (see Eq. ( 18)). We set µ = 1.5 MPa, ρ = 990 kg/m 3 , and three values of λ = 0, 3 and 100 MPa. These different values of λ define the compressibility of the particle, i.e. for λ = 0 the particle is fully compressible whereas for λ = 100 MPa the particle is quasi-incompressible. The deformed particles at vertical strain γ yy = 20% for these values of λ are shown in Fig. 8.

We also note the decrease of the lateral extension as the particle compressibility increases. This extension is negligible for the fully compressible particle (λ = 0 MPa) (see Fig. 8). It can be explained by the fact that for the incompressible particle, the volume of the compressed portion can migrate to the non-contact portion more efficiently as a result of its dense structure [START_REF] Lin | [END_REF]. Fig. 9 displays the second Piola-Kirchhoff stress S yy as a function of the Green-Lagrange strain γ yy . We see that the stress increases as λ increases.

In other words, deforming less compressible particles requires a larger force.

Compaction of a packing of elastic particles

In this section, we investigate the compaction of a packing of elastic particles using the Neo-Hookean material behaviour. By means of MPM simulations, we study the evolution of different packing properties (packing fraction, connectivity etc.). We consider a packing of 300 particles confined inside a rectangular box. The initial configuration is prepared by means of DEM simulations. A uniform distribution of the particle diameters by volume fractions in the range [2,[START_REF] Jean | Proceeding of The International Deep Drawing Research Group[END_REF] mm is introduced. This polydispersity allows avoiding long-range ordering. We simulate the compaction process by moving the top wall downwards at constant velocity of 2 m/s and with a time step of ∆t = 0.1 µs. We consider the Neo-Hookean particles with the same material parameters as in the previous section. The gravitational acceleration is set to be zero in order to avoid stress gradients. There is no friction between the particles, and between the particles and the walls. is the volume of particles and V the total volume, increases by particle shape change and at the end of the compaction nearly the whole space is filled by the particles. The shapes of the particles gradually change from circular to nearly polygonal as shown in Fig. 10. Note that the gaps observed between particles are related to the meshing resolution, which may be increased for a finer discretisation of the contact zone. Moreover, as mentioned before, since the less compressible particles can elongate more, the pores between these particles are more rapidly filled even for a low global deformation. The above feature is more clearly highlighted in Fig. 11. It shows the cumulative volume deformation of the particles defined by ln(V S /V S i ), where V S i is the initial volume of the particles, and the cumulative vertical strain ε as a function of the packing fraction Φ.

The latter is expected to vary due to the elastic volume change of the particles as a result of elastic compressibility of the particles as well as the variation of the total volume V due to particle rearrangements and shape change. Since the width of the box is constant, we have ln

( V S V S i ) = ln ( Φ Φ i ) + ε (20)
with ε = ln(h/h i ), where h i is the initial height of the sample. In Fig. 11, the data for three values of λ coincide up to Φ ≃ 0.8. Beyond this packing fraction, ε and ln(V S /V S i ) vary at different nearly linear rates for each λ value. As expected, this rate increases as λ (or compressibility) decreases. Note that for λ = 100 MPa (quasi-incompressible particles), the volume variation of particles is negligibly small as shown in Fig. 11. Moreover, one can consider Φ ≃ 0.8 as the jamming point above which no particle rearrangements occur anymore, and thus the packing evolution is only due to particle shape change.

The mean coordination number Z is an important parameter, which evolves with the evolution of the packings. Fig. 12 displays Z as a function of Φ. We see that by normalising the excess coordination number Z -Z 0 , where Z 0 is the mean coordination number at the jamming point, by Z 1 -Z 0 , where Z 1 is the coordination number at Φ 1 ≃ 1, and Φ -Φ 0 by Φ 1 -Φ 0 , all data points collapse on a single plot that is well fitted by a power-law function:

Z -Z 0 Z 1 -Z 0 = ( Φ-Φ 0 Φ 1 -Φ 0 ) 0.5 . (21)
A similar power-law behaviour was observed by several authors specially in the case of emulsions and foams [17,24,25] in the following form: These two last equations coincide by setting β = 0.5 and z 0 = (Z 1 -Z 0 )/(Φ 1 -Φ 0 ) 0.5 . As in our study, for all cases, we have Z 0 ≃ 4 for Φ 0 and Z 1 ≃ 5.6 for Φ 1 , one gets z 0 ≃ 3.6. This value is fully consistent with O'Haren et al. ( 2003) [24] who predicted that z 0 is equal to 3.5 ± 0.3 and β ≃ 0.5± 0.03 in 2D. It is also interesting to note that the material behaviour of the particles has almost no effect on these results (see Fig. 12), in agreement with the references [17,24,25], which observed that the evolution of the coordination number with packing fraction is independent of space dimension, interaction potential and polydispersity.

Z -Z 0 = z 0 (Φ -Φ 0 ) β . (22)
The evolution of the applied stress σ beyond the jamming point allows for a macroscopic analysis of the packing evolution. Fig. 13 shows σ , computed from the contact forces acting on the bottom wall and normalised by the particle P-wave modulus M p (M p = λ+2µ) as a function of Φ. We note a nonlinear behaviour for three cases with different rates. As in the case of one particle simulations, we observe also that the required force to compress the packing increases with the particle compressibility. These observations may be explained by the fact that beyond the jamming point the packing behaves almost like a continuum medium as there are no more particle rearrangements. This assumption leads to a logarithmic relation between σ and Φ (see Appendix B):

σ M P = -Z Φ Z + Mp c 1 Kp (ln(Φ) + c 2 ) , (23)
where K p is the particle bulk modulus (K p = λ + µ in 2D), c 1 is a parameter depending on the particle material behaviour and c 2 is a constant term. The coordination number Z is also defined as a function of Φ by Eq. ( 21). The predictions of this model [START_REF] Lin | [END_REF] are in good agreement with our MPM simulations shown in Fig. 13 with c 2 ≃ 0.27, and c 1 ≃ 0.01 for λ = 100 MPa, c 1 ≃ 0.23 for λ = 3 MPa and c 1 ≃ 1.3 for λ = 0 MPa. Note that, although this model seems to predict well the applied stress σ as a function of Φ, at high packing fractions, the Neo-Hookean particles can only overfill the remained little pores for much higher stresses (involving smaller and smaller radii of curvature) but when the packing fraction tends to 1, the corresponding applied stress should tend to infinity. Hence, this model does not hold at high values of the packing fraction. In our simulations, to resolve correctly the small radii of curvature at the contact zones between particles, one should refine in the same proportion the discretisation.

In order to analyse the particle deformations and volume change during the compaction of the packings, we also consider the evolution of the distribution of the equivalent von Mises strain γ eq and the Jacobian of deformation gradient 0 J (see Eq. ( 9)), respectively.

Here, γ eq is defined as

γ eq = √ 2 3 γ d : γ d , (24)
with γ d = γ -1 3 Trace(γ) I. To characterise the shapes of these distributions, we consider here their standard deviation and excess kurtosis. Fig. 14 displays the standard deviation and excess kurtosis of γ eq as a function of Φ. The standard deviation of γ eq for the three values of λ coincide up to Φ ≃ 0.8, but beyond this value they vary at different rates. In the quasi-incompressible case (λ = 100 MPa), the standard deviation is larger than in compressible cases. This trend can be explained by the occurrence of more important stress chains between particles when their compressibility decreases.

The values of the excess kurtosis of γ eq distributions for several values of λ coincide, and they tend to a positive value about 2 beyond the jamming point. This value is compatible with a Leptokurtic distribution, which shows heavier tails than the normal distribution.

Finally, the standard deviation and the excess kurtosis of the Jacobian of deformation gradient 0 J as a function of Φ are shown in Fig. 15. As expected, the standard deviation for λ = 100 MPa (quasi-incompressible particles) is almost zero since there is no particle volume change. The standard deviation is larger for the fully compressible particles (λ = 0 MPa). It is due to the larger possible deformation of the particles in this case. However, the kurtosis is nearly zero, meaning that the 0 J distributions are nearly normal.

Conclusion

In this paper, we improved our approach for modelling softparticle systems developed in [11]. In this novel approach, the finite strain formulation is used in the context of the implicit Material Point Method (MPM). The MPM allows one to take into account the realistic mechanical behaviour of individual particles. Coupling the MPM with the Contact Dynamics (CD) method makes it possible to deal correctly with frictional contacts between particles.

It was shown that two MPM formulations (infinitesimal and finite strain) are similar. The finite-strain formulations can host more complex constitutive behaviours such as hyperelasticity for the particles. Furthermore, to improve computational performance, a parallelisation procedure was proposed in the framework of this algorithm. Although the efficiency of this procedure declines with increasing number of processors, it is still useful for decreasing the computational cost.

The uni-axial compaction of a packing of soft particles was simulated using MPM by considering several values of particle compressibility (from quasi-incompressible to fully compressible particles). The packing with more compressible particles can undergo larger deformations under the action of lower compressive stress due to considerable particle volume changes and occurrence of weaker stress chains between particles. It was shown that this stress beyond the jamming state varies logarithmically with packing fraction. This behaviour was explained by introducing a simple model. Another interesting result of this work concerns the evolution of the coordination number, which can be related to the packing fraction by a power-law function beyond jamming transition.

Fig. 1 .

 1 Fig. 1. Geometry of a contact between two soft particles discretised in multi-mesh MPM algorithm. The solid points represent the potential contact nodes; see text.

Fig. 2 .

 2 Fig. 2. Contact conditions: (a) Velocity-Signorini complementarity condition as a graph relating the normal relative velocity v n and normal force f n ; (b) Coulomb friction law as a graph relating the tangential velocity v t and friction force f t ; µ is the coefficient of friction. The dashed lines represent the linear relations obtained from a linear combination of the equations of dynamics; see text.

  node = M β v β node at the common node i, we get the following linear relation:

Fig. 3 .

 3 Fig.3. Scheme of the parallelisation procedure of the MPM-CD framework; N is the number of the particles and, f c , m, n and t represent the contact force, mass, normal and tangential vectors, respectively, for the potential contact nodes of each particle.

Fig. 4 .

 4 Fig. 4. Load balance algorithm. A set of circular particles (in 2D) is shown, R i 2 being written in each particle, and distributed in 3 processes. The total workload ∑ R i 2 is

Fig. 5 .

 5 Fig. 5. Measured speedup for the MPM-CD simulation of 300 particles. The line represents perfect scaling.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Geometry of a single Neo-Hookean particle and it's deformed configurations at vertical strain γ yy = 20% for different values of λ. Although the deformed particles seem to touch the two walls only over a short segment at the centre, the material points belonging to the boundary elements between the particle and the bottom and top walls are actually within the contact zone. The observed gap is due to the background mesh element thickness.

Fig. 10 .

 10 Fig. 10. A snapshot of the initial configuration (a) and three snapshots of the compaction of a packing of soft particles with Neo-Hookean behaviour for packing fraction of 0.97 and several values of the λ (b-d). Note that, despite the same value of the packing fraction, the packing volumes are different due to the different compressibilities of the particles. The black points represent the material points.

Fig. 10

 10 represents the snapshots of the compaction test for different values of λ . The packing fraction Φ = V S /V , where V S

Fig. 11 .

 11 Fig. 11. Evolution of the cumulative vertical strain ε and, the total particles volume change ln(V S /V S i ) as a function of the packing fraction Φ for several values of λ.

Fig. 12 .

 12 Fig.12. Evolution of the excess mean coordination number Z -Z 0 as a function of excess packing fraction Φ -Φ 0 for several values of λ. The solid line is power-law fit (Φ -Φ 0 ) 0.5 ; see Eq. (21).

Fig. 13 .

 13 Fig. 13. The normalised applied stress as a function of packing fraction normalised by the particle P-wave modulus M p for several values of λ. The lines represent the predicted behaviour by the model of compaction introduced in this paper; see Eq. (23).

Fig. 14 .

 14 Fig. 14. Evolution of the standard deviation SD and excess Kurtosis exKurt of the equivalent von Mises strain, γ eq as a function of the packing fraction Φ for several values of λ.

Fig. 15 .

 15 Fig. 15. Evolution of the standard deviation SD and the excess Kurtosis exKurt of 0 J as a function of the packing fraction for several values of λ.

  This continuum body is subjected to prescribed displacements and forces on the disjoint complementary parts of the boundary

	∂Ω u t (the Dirichlet boundaries) and ∂Ω f t (the Neumann bound-
	aries), both in the actual configuration, respectively. The boundary
	conditions are then defined by

t x, t) is the Cauchy stress tensor, b( t x, t) represents the body force and a( t x, t) denotes the acceleration at position t x and time t.

Acknowledgements

This work/project (ID 1502-607) was publicly funded through ANR (the French National Research Agency) under the ''Investissements d'avenir'' programme with the reference ANR-10-LABX-001-01 Labex Agro and coordinated by Agropolis Fondation, France under the frame of I-SITE MUSE (ANR-16-IDEX-0006). We are also grateful to the genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing computing resources.

Conflict of interest

Saeid Nezamabadi, Xavier Frank, Jean-Yves Delenne, Julien Averseng and Farhang Radjai state that there are no conflicts of interest.

Appendix A. Definitions of K and R

The implicit integration in the context of MPM takes into account the discretised equation of the motion:

By considering that the external force at time t + ∆t is known, f ext (t + ∆t), and by assuming an incremental-iterative Newton solution strategy, the linearised equation of motion at iteration k is

where

Note that in the last relation, t+∆t H p can be obtained using Eq. ( 9) and the constitutive relation (8).

Appendix B. Relation between the applied stress, σ, and the packing fraction, Φ, for a packing under uniaxial compression

We assume that the packing of particles behaves almost as a continuum medium beyond the jamming point under uniaxial compression. Hence, in this range the applied stress σ may be related to the cumulative vertical strain ε through an effective P-wave modulus M:

Here, the particle and pore volume changes can be assumed to be the same, implying that the effective P-wave modulus is proportional to the packing fraction: M = ΦM p with M p the particle P-wave modulus. One may further assume that the particle bulk modulus K p relates the volume increment dV S of particles to the effective stress increment dσ S in particles:

σ can be related to σ S as follows:

where c 1 is a material constant to determine. Given that dε = dV S /V S -dΦ/Φ and using Eqs. (B.1), (B.2) and (B.3), the following differential equation to solve is obtained:

By knowing that there is a relationship between Z and Φ (see Eq. ( 21)), the integration of the differential equation (B.4) is given:

where c 2 is the integral constant.