Pierre Rognon

Cyprien Gay

Soft Dynamics simulation: 2. Elastic spheres undergoing a T1 process in a viscous fluid

Keywords: numbers: 02.70.Ns, 82.70.-y, 83.80.Iz

Robust empirical constitutive laws for granular materials in air or in a viscous fluid have been expressed in terms of timescales based on the dynamics of a single particle. However, some behaviours such as viscosity bifurcation or shear localization, observed also in foams, emulsions, and block copolymer cubic phases, seem to involve other micro-timescales which may be related to the dynamics of local particle reorganizations. In the present work, we consider a T1 process as an example of a rearrangement. Using the Soft dynamics simulation method introduced in the first paper of this series, we describe theoretically and numerically the motion of four elastic spheres in a viscous fluid. Hydrodynamic interactions are described at the level of lubrication (Poiseuille squeezing and Couette shear flow) and the elastic deflection of the particle surface is modeled as Hertzian. The duration of the simulated T 1 process can vary substantially as a consequence of minute changes in the initial separations, consistently with predictions. For the first time, a collective behaviour is thus found to depend on another parameter than the typical volume fraction in particles.

I. INTRODUCTION

Many materials are made of particles in a surrounding fluid. Among them foams, emulsions, granular matter, colloidal suspensions and micro gels are of daily use. A great deal of research revealed their complex behaviors including elastic, plastic and viscous characters [START_REF] Weaire | The Physics of Foams[END_REF][START_REF] Coussot | Rheometry of pastes, suspensions, and granular materials[END_REF][START_REF] Stickel | [END_REF]4]. This complexity results from the wide range of particle properties and particle interactions involved. Great hints to comprehensive rheological models were obtained by considering the dynamics of a single particle. Thus emerged the time m/RP for a single grain of mass m, accelerated by the normal stress P (force P R 2 ), to move over a distance comparable to its own size R [5,6,7], the time η/P for a grain immersed in a fluid of viscosity η subjected to the same normal stress [8], and the relaxation time η R 2 /σ for a bubble or a droplet with surface tension σ in a viscous fluid [9,10,11]. The effective viscosity was expressed as an empirical function of these microscopic timescales [8], thereby providing robust scaling expressions for various properties of grains [12,13,14] and bubbles [15,16,17,18,19,20].

Nevertheless, particulate materials exhibit some uncommon rheological properties which seem to involve other timescales. Oscillatory shear experiments [4], and more generally the delayed adaptation of the shear rate to a sudden change in the applied stress [21,22,23,24,25,26,27] reveal long internal relaxation processes. Other observations such as a critical shear rate below which no homogeneous flow exists [5,6,8,22,28], or the coexistence of liquid and solid regions (shear localization, shear banding, cracks) in emulsions [29,30], foams [31,32,33,34], wormlike micelles [35,36] and granular materials [6,14,37,38] also point to a complex internal dynamics. Usually, this internal * Electronic address: cyprien.gay@univ-paris-diderot.fr dynamics is qualitatively understood as the competition between external solicitations that the particles experience and their ability to move within their neighborhood [22]. Such a mechanism is the core of the definition of the jamming transition in glassy systems, which is a subject of intense debate [START_REF] Isa | Journal of Physics: Conference Series[END_REF][START_REF] Hecke | [END_REF]41,42,43,44]. Getting new insights into reorganization micro-timescale should therefore clarify the origin of such properties and should also provide useful hints to refine and generalize existing models of the material response.

A common reorganization process is the separation of particles while other particles approach and fill the void. When they involve four particles, these events, usually refered to as T 1 processes when dealing with foams and emulsions, occur in deformed regions at a frequency proportional to the deformation rate (see for example [31,32,45,46,47,48,49,50,51,52]). They relax some stress and dissipate some energy. The relation between the duration of a T 1 and the local stress is thus expected to affect the rheological behaviour of the material. For dry foams, the T 1 dynamics has been shown to result from the surface tension and the surface viscosity [53]. The stretching ability of the interfaces avoids the need to squeeze violently the fluid between the approaching bubbles, and its viscous dissipation is thus negligible. By contrast, the T 1 dynamics is less well described in wet foams or other less concentrated systems. A comprehensive description of their dynamics requires a careful description of the particle interaction. For instance, visco-elastic and even adhesive properties of particles were shown to be important [54,55,56].

In this paper, we show that squeezing the liquid between close particles (here in three dimensions) can give rise to long relaxation times. To this aim, we do not focus on a specific material which would include the interaction between solid grains, bubbles, droplets or colloidal particles. Rather, we address the ubiquitous situation of elastic-like particles in a Newtonian fluid. We consider a simple system of in-plane spheres undergoing a T 1 process, as depicted on Fig. 1. We discuss under which circumstances a T 1 process should in-deed occur and (if it does occur) the relative contribution to the dynamic of the normal approach and separation versus the tangential sliding of particles [63]. In three dimensions, for a dry foam or a concentrated emulsion, such a T 1 process with four topologically active particles can in fact be decomposed into two topologically simpler processes involving five particles. However, whatever the exact process, the dynamics will still involve normal motions and tangential sliding (as well as rotation in general). Because normal motions are stronger, as we show below, we believe that no essential new phenomenon will emerge from other reorganization processes as compared to the time scale evidenced in the present work. This paper is the second of a series which presents the physics of materials made of close-packed elastic-like particles immersed in a viscous fluid. In the first paper [57], we focused on the normal separation of two particles, and we showed that the flow between their close surfaces interplays with the particle deformation in a non-trivial manner. As this feature was ignored so far in existing discrete element simulations such as Molecular Dynamics [58] (for elastic grains without a surrounding fluid) for Stokesian Dynamics [59] (for non-deformable grains in a viscous fluid), we are introducing a new simulation method, named Soft-Dynamics, to account for it. In this paper, we include the tangential interaction, and we provide the main steps of the implementation of the Soft Dynamics method for the present context. This will serve as an introduction to the principle of larger scale simulations with this new method, which will include both particle rotation and boundary conditions, and which should constitute a promising tool for investigating the collective behaviors of many complex materials.

As we shall see, the geometry addressed in the present paper, although rather symmetric (the centers of the three dimensional particles are arranged within a plane, at the vertices of a losange), proves sufficiently rich to reveal how minute changes in the system configuration have an essential influence upon its dynamics.

(a) (b) (c) θ θ ≃ π/6 θ ≃ π/3
FIG. 1: Schematic representation of a T 1 process with four in-plane spheres. Due to the applied forces, the group of four particles swap neighbours. Two particles separate while two other particles establish contact. Meanwhile, the other particle pairs reorient, as shown by the evolution of angle θ from about π 3 to about π 6 .

II. MODELLING PARTICLE INTERACTIONS

When addressing the question of a T 1 process between elastic spheres in a viscous fluid, see Fig. 1, most of the interactions have already been described in the first paper [57]. The only new feature is particle sliding and, correspondingly, tangential forces. Hence, quantities such as viscous friction coefficients or spring constants are now tensorial. We express these interactions in the present section. Let us recall that we deal with three dimensional particles.

A. Pairwise interactions

As discussed in detail in the first paper [57], because we consider rather dense systems where each particle is close to several other particles (surface-to-surface gap much smaller than the particle size), we simply discard long-range, manybody interactions [64]. Furthermore, under such thin gap conditions, the interacting region between particles is much smaller than the particle size and the interactions between a particle and its neighbours are mostly independent from each other and can therefore be treated as a sum of pairwise interactions.

Particle i is subjected to some force -→ F i j by its neighbouring particles j, which can be decomposted into (i) the local pressure field in the fluid that results from a viscous lubrication interaction and (ii) a remote interaction (such as a damped electrostatic interaction, steric repulsion, van der Waals interactions, disjoining pressure, etc):

-→ F i j = -→ F vis i j + - → F rem i j (1)

B. Particle surface deflection

A fraction of the above force -→ F i j exerted by particle j transits through a small portion of the surface of particle i and deflects it elastically.

In practice, the viscous component -→ F vis of the force is entirely transmitted by the surface of particle i. The effect of the remote component -→ F rem is more subtle. Electrostatic forces between surface charges act upon the surface and contribute entirely to the elastic deflection. By contrast, Van der Waals interactions also act directly within particle i. However, most of such interactions occur within a depth comparable with the inter-particle gap, which is always much smaller than the depth of the region that is deformed elastically (see paragraph II F below).

Hence, for simplicity, it is reasonable to assume that the total force -→ F i j between both particles entirely contributes to the elastic surface deflection:

-→

F ela i j ≃ -→ F i j (2)
The expression of -→ F ela i j in terms of the corresponding surface deflection is discussed in paragraph II F below.

C. Force balance for each particle

The sum of all forces applied to particle i, both the external force -→ F ext i and the pairwise forces -→ F i j is equal to the mass m i times the acceleration. This is assumed to vanish due to the dominant effect of the fluid viscosity over inertia:

-→ F ext i + ∑ j -→ F i j = m i - → X i = - → 0 (3)
where -→ F ext i is an external force acting on grain i (such as gravity) and where the sum runs over the neighbours of particle i. In principle, there is another equation, similar to Eq. ( 3), for the torques applied to particle i. But as mentioned earlier [63], this is not needed for the present symmetric T 1 configuration such as that of Fig. 1.

The Soft Dynamics method [57] simulates the dynamics of such a system, determined by the system of Eqs. (??) for all interactions and Eqs (3) for all particles i. In the present work, for simplicity, we omit the remote interactions in Eq. (??) as we did before [57].

In order to specify the elastic and viscous forces, let us now describe the geometry and the kinematics of the interacting region between a pair of neighbouring particles.

D. Contact geometry and kinematics

Let i and j denote two interacting particles, as depicted on Fig. 2. As compared to the first paper, the positions of the particles centers are now vectors, labeled -→ OX i and --→ OX j , and -→ X i j = --→ OX j --→ OX i is the center-to-center vector. The deflections of the particle surfaces are also vectors, labeled -→ δ j i and -→ δ i j . Since all particles are identical and since the (lubrication) forces are pairwise and act locally, facing deflections are symmetric:

-→ δ j i + -→ δ i j = 0. Thus, for simplicity, we shall use the total deflection -→ δ i j = -→ δ j i --→ δ i j for each pair of interacting particles. The unit vector normal to the contact can be expressed as

-→ n i j = -→ X i j - - → δ i j | -→ X i j - - → δ i j | (4)
The gap h i j between both particle surfaces depends on both the center-to-center vector -→ X i j and the total deflection -→ δ i j :

h i j = | -→ X i j - - → δ i j | -2R, (5)
Similarly, the relative velocity of the material points that constitute each particle surface, -→ v s , involves the translation velocity of the particles (as already mentioned [63], the particles do not rotate in the present situation) and the evolution of the surface deflection:

-→ v s = - → X i j - - → δ i j . (6)
In order to specify viscous and elastic interactions, we will need to deal with projectors and tensors. We will use the symbol "•" for the tensor product (contraction of one coordinate index), and -→ u T will denote the transposed of vector -→ u .

Hence, -→ u T • -→ v = -→ v T • -→ u will be the scalar product of -→ u and -→ v , and -→ u • -→ v T their outer product, which is a tensor. In particular, we will make use of tensor α defined as the projector onto the normal direction:

α = -→ n • - → n T ( 7
)
0000000000000 0000000000000 0000000000000 1111111111111 1111111111111 1111111111111

h (a) δ n a (c) (b) -→ n -→ δ i j -→ δ j i --→ OX i --→ OX j -→ v s • a -→ v s • (i -a)
FIG. 2

E. Viscous force

For a pair of close spheres, as discussed earlier [57], the fluid region that mediates most of the force between both particles has a large aspect ratio, and the flow is essentially parallel to the solid surfaces: the lubrication approximation can be used (see for example [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]). As before, the fluid inertia is negligible (low Reynolds numbers) and the viscous force -→ F vis acting on the surfaces depends linearly on their relative velocity -→ v s :

-→

F vis = Z • - → v s (8) Z = ζ α + λ (1 -α) (9) ζ = 3πηa 4 2h 3 (10) λ = πηa 2 h ( 11
)
where the interparticle friction tensor Z has two components (normal and in-plane), expressed in terms of the unity tensor 1 and the projector α defined by Eq. ( 7). The normal viscous friction ζ is related to the Poiseuille flow induced by squeezing or pulling [57] (see Fig. 2b), while the in-plane friction coefficient λ reflects the tangential motion (sliding) between both particles, which generates a Couette (shear) flow in the gap (see Fig. 2c).

F. Elastic force

Let us assume that the size a (discussed in the next paragraph) of the interacting region between particles i and j is known. Then, as before [57], the force depends linearly on the surface deflection, but this time the relation is tensorial:

-→ F ela = aE • - → δ , (12)
E = E c n α + c t (1 -α) , (13)
where the tensorial proportionality constant E is essentially the (scalar) Young modulus E, but incorporates geometrical constants on the order of unity c n and c t for the normal and tangential responses, respectively. The elastic response of bubbles and droplets were found to deviate from such a Hertz elasticity [54,55]. Although they have no bulk elasticity, the surface tension σ confers them some elastic-like properties, and the elastic-like force mainly depends on the deflection δ , the size a of the interacting region and an effective Young modulus which scales like σ /R.

G. Size of the interacting region

The size of the interacting region, again [57], depends either on the gap thickness h (Poiseuille regime) when the particle surface is weakly deflected, or on the normal force (Hertz regime) when the particle surface can be considered planar. In the first case, it can be expressed as a ≈ √ 2 R h. In the second case, it is essentially independent of the tangential force [START_REF] Johnson | Contact Mechanics[END_REF] and can thus be expressed in terms of the normal deflection:

a ≈ R|δ n | = R | -→ n T • - → δ |.
As explained earlier [57], for the purpose of the Soft Dynamics method, we interpolate between both behaviours of a in a simple manner:

a(h, δ n ) = R(2h + |δ n |). (14)
The choice of this interpolation is not physically supported, but it does not affect assymptotic behavior in both limits.

III. METHOD OF THE SOFT-DYNAMICS SIMULATION

The Soft-Dynamics method aims at simulating the time evolution of a system of elastic particles and in a viscous fluid, such as depicted in previous sections. Like usual discrete simulation methods, the motion of each particle center results from the force balance, Eq. ( 3). The specificity is that the interaction evolution results from the decomposition of the center-to-center distance given by Eq. ( 5). As illustrated previously [57], this generates a Maxwellian contact dynamics through the combination of the elastic surface deflection and the viscous response of the fluid in the gap: it is possible to move the center-to-center distance -→ X i j while keeping constant the deflection -→ δ i j , and vice-versa. But as compared to a classical Maxwell behaviour, the elastic element does always behave linearly (Hertzian contact in the strong deflection regime), and the viscous element does not have a constant value, as it depend on the geometry of the gap, see Eqs. (8)(9)(10)(11).

The Soft-Dynamics method consists in calculating the rate of change of all center positions -→ OX i and all gap deflections -→ δ i j as a function of their current values, and integrating them over a small time step.

A. Equations of motion

The system satisfies one equation per interaction, namely Eq. (??), and one equation per particle, namely Eq. ( 3). We shall now see how it is possible to derive equations of motion.

For this, we need to express the unknowns velocities -→ δ i j and -→ X i j in terms of the current state of the system.

From Eqs. ( 4), ( 5) and ( 14), it appears that the size a of the interacting region can be expressed as a function of -→ X i j and -→ δ i j . It then follows from Eqs. ( 4) and ( 12) that the elastic force -→ F ela can also be expressed as a function of -→ X i j and -→ δ i j :

-→

F ela = -→ F ela ( -→ X i j , - → δ i j ) (15)
As a result, its time-derivative -→ F ela i j can be expressed as a sum two terms: one of them is linear in -→ X i j while the other is linear in -→ δ i j . The (tensorial) coefficient of each of these two terms is a function of the current system configuration, i.e., of all particle and gap variables -→ OX i and -→ δ i j . Now, it follows from

Eqs. ( 6), ( 8) and (??) that -→ δ i j is an affine function of -→ X i j :

-

→ δ i j = - → X i j + Z -1 i j • -→ F ela i j - - → F rem i j ( 16
)
where Z i j , -→ F ela i j and -→ F rem i j depend on the current system configuration. Hence, -→ F ela i j can be expressed as an affine function of -→ X i j :

- → F ela i j = G i j • - → X i j - - → b i j ( 17
)
where the coefficients G i j and -→ b i j depend only on the current system configuration. The detailed calculation of these coefficients is provided in Appendix A.

From this, the time derivative of Eq. (3) yields a system of equations for the particle center velocities. The equation that corresponds to particle i reads:

∑ j G i j • ( - → OX j - - → OX i ) = ∑ j -→ b i j - - → F ext i . ( 18
)
where the sums run over all neighbours of particle i.

Note that because G ji = -G i j and -→ b ji = --→ b i j , and if we assume that the sum of all external forces vanishes,

∑ i -→ F ext i = 0, (19)
then the sum of Eqs. ( 18) for all particles i vanishes. In other words, these vector equations are not independent: one of them must be replaced, for instance, by the condition that the average particle velocity is zero:

∑ i - → OX i = 0 ( 20
)
Let us consider the system of Eq. ( 20) (or a similar one) together with Eqs. (18), taken for all particles i except one. This system of equations can be inverted to obtain the particle center velocities -→ OX i . The gap velocities -→ F ela i j are then calculated from Eq. ( 17).

B. Choice of a numerical step

Gaining the center velocity -→ OX i requires to solve the linear system (18). Standard and efficient procedures are available to inverse it. We used a second order Newtonian scheme for the numerical integration of particle position and as well as deflections. A typical time in the problem is the Stokes time τ taken by a single particle submitted to a typical force F to move over a distance R in a fluid with viscosity η, see Eq. (29) below. The numerical time step is set to 10 -3 in units of τ for all simulations. Other numerical schemes, such as Runge Kutta method, should make simulations faster. Furthermore, a study of the optimal required time step will be necessary when dealing with significantly more than only four particles.

IV. T 1 DYNAMICS

Let us now use the Soft-Dynamics method to simulate a single T 1 process. The system is depicted on Fig. 3: initially, particles B and C are aligned horizontally, with a small gap h 0 , while particles A and D are aligned vertically. The diagonal gaps (between A and B, etc) have thickness h 0 too.

A horizontal force N x is applied on particles B and C while a vertical force N y is applied on A and D. Various evolutions are possible depending on these two forces, which may or may not give rise to a T 1 process (see Fig. 3). Basically, a T 1 occurs only if the interaction between particles B and C is tensile. The criterion for the occurrence of a T 1 process will be derived below, as well as a scaling for its dynamics. The duration of a T 1 will then be measured from the simulation.

A. Theoretical predictions

The T 1 process, which consists in a separation of the horizontal pair of particles (BC) and an approach of the vertical pair of particles (AD), implies some sliding of the diagonal pairs (see Fig. 1).

At the early stages of the process, when θ ≈ π/3, the external forces N x and N y can be expressed in terms of the normal forces in the horizontal (N h ) and diagonal (N d ) pairs of particles, and in terms of the sliding force S d in the diagonal pairs:

N x = N h + 2 1 2 N d -2 √ 3 2 S d ( 21
)
N y = 2 √ 3 2 N d + 2 1 2 S d ( 22
)
In fact, as we shall now see, the tangential force S d is much smaller than the normal forces. To show this, let us first notice that the tangential velocity is related to the angle θ defined on Fig. 1: v t ≃ -R θ . In the Poiseuille regime, the particles surfaces are weakly deflected and the horizontal and diagonal gaps are related to angle θ through R+ 1 2 h h = (2R+ h d ) cos θ . Hence, the gap variations obey 1 2 ḣh ≈ ḣd cos θ -2R θ sin θ , FIG. 3: Schematic evolution of four particles subjected to external forces. Force N x is horizontal and acts on particles B and D. Force N y is vertical and acts upon particles A and C. Both N x and N y can be either compressive (> 0) or tensile (< 0). Regimes (I) and (II) correspond to compressive forces. In regime (I), the configuration remains mostly unaltered. By contrast, a topological rearrangement (T 1 process) occurs when N y √ 3N x , which corresponds to region (II). When N x or N y is tensile, the four beads do not remain together, as can be seen from the time evolutions sketched for regimes (III)-(VI). On the whole, a T 1 process always occurs when N y √ 3N x (regimes II, III and IV, white region). It is followed by particle separation when N x is tensile (regimes III and IV). By contrast, no T 1 process occurs when N y √ 3N x (regimes I, V and VI, light grey region).

i.e.:

1 2 ḣh ≈ ḣd cos θ + 2 v t sin θ (23)
Let us now transform each term of the above equation by expressing it as a function of the corresponding normal or tangential force by using the appropriate friction coefficient as defined by Eq. ( 9):

- 1 2 N h ζ ≈ - N d ζ 1 2 + 2 S d λ √ 3 2 (24)
The relative magnitude of friction coefficients ζ and λ can be derived from Eqs. ( 10) and ( 11):

ζ λ = 3 2 a h 2 . ( 25
)
where the size a of the interaction region is given by Eq. ( 14).

We thus have ζ /λ ≈ R/h in the Poiseuille regime and ζ /λ ≈ Rδ n /h 2 in the Hertz regime. Hence, except for very large gaps h comparable to the particle size R, the normal friction is much larger than the sliding friction: ζ ≫ λ . It follows that

S d ≃ 1 2 √ 3 λ ζ (N d -N h ) (26)
can be neglected in Eqs. (21)(22). Hence, the interaction force within the horizontal pair BC depends only on the applied forces:

N h ≈ N x - 1 √ 3 N y (27)
This implies that, as pictured on Fig. 3, the gap will open and the T 1 will proceed whenever N h is tensile, i.e., when N y N x √ 3 (white region of the diagram). By contrast, the particles will not swap neighbours when N y N x √ 3 (light grey region).

When N h is indeed tensile, we now wish to determine how long it takes for the horizontal pair of particles to separate.

The dynamics of such a normal motion was detailed in Ref. [57]. Let us define the reduced force

κ = |N h | ER 2 (28)
and the Stokes time

τ 6πηR 2 |N h | (29)
With the force N h acting within the horizontal pair BC, the initial configuration (gap h 0 ) corresponds to the Poiseuille regime if h 0 h HP and to the Hertz regime if h 0 h HP , where

h HP = R κ 2/3 (30)
The corresponding rate of change of the gap [57] can be expressed as:

ḣ = h τ (Poiseuille, h > h HP ) (31)
ḣ = h 3 τ R 2 κ -4 3 (Hertz, h < h HP ). (32)
Integrating these equations yields the typical time ∆ required to achieve the separation of the horizontal pair BC of particles from an initial gap h 0 to a much larger gap h f ≈ R:

∆ ≃ τ ln h f h 0 ≈ τ (Poiseuille, h 0 > h HP ) (33)
∆ ≃ τ κ 4 3 R 2 h 2 0 - R 2 h 2 HP + τ ln h f h HP ≈ τ κ 4 3 R 2 h 2 0 ≫ τ (Hertz, h 0 < h HP ) (34)
Once the gap h h of the horizontal pair BC becomes comparable to R, the diagonal pairs such as AB slide rather quickly (since their λ ≪ ζ ), and soon the gap h v of the vertical pair AD becomes significantly smaller than R. The time it then takes to reach the same value h 0 is again comparable to ∆.

Hence, the order of magnitude given by Eqs. (33)(34) for the time ∆ is typically the expected order of magnitude for the duration of the entire T 1 process. We will now test this prediction by comparing it with the simulation results.

B. Result from simulations

We implement the Soft Dynamics method to simulate a T 1 process such as that depicted on Fig. 3, varying the two control parameters we pointed out above: the initial gap 10 -3 < h 0 /R < 0.8 and the reduced force 10 -4 < κ < 0.1. For simplicity, there is no horizontal force (N x = 0). The reduced force given by Eq. ( 28) is then equal to κ = |N y |/ER 2 √ 3 and the Stokes time is τ = 6π √ 3ηR 2 /|N y |. Figure 4 displays the variations of several quantities in the course of a T 1 process with a given set of parameters (h 0 = 10 -2 R, κ = 3.10 -3 ). In order to avoid discontinuities in the simulation, the force N y is increased from zero to its nominal value within a time τ, and remains constant thereafter. From a macroscopic point of view, for instance through the variation of the angle θ , the system seems to be almost blocked (θ ≈ π/3) for a significant amount of time (t 100τ). It then starts moving to reach its final configuration (θ ≈ π/6) where it remains thereafter (t 250τ). During the "blocked" phase, the applied force N y is transmitted through the diagonal interaction such as AB, thereby inducing a tensile force N h ≈ -N y / √ 3 in the horizontal pair BC. Hence, despite the overall "blocked" appearance of the system, the horizontal gap h h between particles B and C slowly increases from its initial value h 0 . Correspondingly, the horizontal friction decreases.

The fast moving period starts as soon as this friction is low enough. Particles B and C then separate quickly while particles A and D in the vertical pair approach each other, thereby giving rise to sliding friction on the diagonal interactions. As particles A and D approach, the corresponding gap h v decreases and the friction increases. This approach then slows down. Thus, although the system keeps moving, it appears to reach a new "blocked" configuration, with no more sliding or horizontal traction, but only a vertical compression.

As the system subjected to a constant force keeps moving, we need to arbitrarily define the end of the T 1 process. Among various possible choices, we shall here consider that the T 1 process is completed when the vertical interaction transmits most of the applied force (N v = 0.99 N y ). The resulting duration of the T 1 process is plotted on Fig. 5 (other criteria would yield similar results). The first observation is that, for the range of initial gaps and particle stiffnesses we consider, the duration of the T 1 is distributed over a wide range of time scales, roughly between 3τ and 10 3 τ. Next, we observe that these results match our theoretical predictions reasonably: if the horizontal pair BC is in the Poiseuille regime, the T 1 duration ∆ scales like τ ln(R/h 0 ). It thus depends on particle radius, on the applied force and on the fluid viscosity through τ, as can be seen from Eq. ( 29), and slightly on the initial gap through the logarithmic factor. The T 1 duration is then just a few times larger than the Stokes time τ;

• if pair BC is in the Hertz regime, the T 1 duration ∆ scales essentially like τ κ

4 3 R h 0 2
, which implies a much stronger dependence on h 0 , and longer durations since the particles are soft. In this case, ∆ can be much longer than the Stokes time τ.

Note that in the latter case, the separating pair of particles leave the Hertz regime and enter the Poiseuille regime in the late stages of separation (h > h HP ). However, because the evolution is much slower in the Hertz regime, see Eqs. (31)(32), these late Poiseuille stages contribute very weakly to the overall T 1 duration ∆.

In summary, the numerical result for the duration of a T 1 process presented on Fig. 5 are compatible with Eqs. (33)(34). They demonstrate that the duration of a T 1 is hardly larger than the Stokes time τ given by Eq. ( 29) as long as the surface deflection is small compared to the inter-particle gap (Poiseuille regime) and thus depend mainly on the applied force, on the fluid viscosity and on the particle size. Remarkably, in the opposite regime where the deflection is larger than the gap (Hertz regime), the T 1 duration depends strongly on the interparticle gap and can reach very large values, as illustrated by Fig. 5.

V. CONCLUSION: BEYOND VOLUME FRACTION

In this paper, we studied one of the simplest reorganization processes for immersed, closed-packed, elastically deformable particles in a simple geometry. We showed that the time needed for this process results principally from the viscous flow of the fluid into or out of the gap between pairs of almost contacting particles: it is always mostly driven by the normal approach or separation, while the role of tangential sliding is negligible.

We also showed that the time needed can be very long when particles are close or soft (more explicitely, when the gap is much thinner than the particle surface deflection). This is the central result of the present study and, as we show below, it pleads towards going beyond the sole usual volume fraction to describe the state of a particulate material.

A. Volume fraction and interparticle gap

Let us consider four particles in a compact configuration such as that on Fig. 1a (angle θ ≈ π/3). More precisely, let us consider two variants of this configuration, with two different values of the interparticle gap h 0 , say h 0 ≃ 10 -2 R and h 0 ≃ 10 -3 R. Let us now apply weak forces (say κ ≃ 10 -2 ). In both situations, because the force is weak and the gaps are small, the center-to-center distances are almost identical. Hence, both situations cannot be distinguished at first sight.

Yet it can be seen from Fig. 5 that the duration of the T 1 process will then differ substantially.

Similarly, with a large, disordered assembly of grains, it is anticipated that there will exist different situations where the volume fraction is almost identical but where a change in the typical value of the interparticle gaps causes a dramatic alteration of the delay ∆ after application of the stress for the system to set into appreciable motion. This conjecture will be tested in a future work, using simulations with a large number of particles.

B. Dilatancy and permeation

In even larger samples of granular materials in a compact state, it is anticipated that the need for some additional fluid to enable reorganization processes (a phenomenon called dilatancy and illustrated on Fig. 1) will become the main source of delay: fluid from loose or particle-free regions needs to permeate through the granular material which behaves as a porous medium [27]. Describing such a phenomenon requires introducing the liquid pressure, and is not included in our simulation so far.

C. Towards other materials

As such, the present work applies to soft, plain, elastic particles (such as elastomer beads or latex particles) immersed in a very viscous fluid. We showed that a possible physical origin for a delay in the system response is the viscous flow in the thin gap between neighbouring particle surfaces.

In other materials, however, other ingredients may also influence this delay or even become dominant. For instance, for objects enclosed in fluid interfaces (vesicles, onions, bubbles, droplets, etc), phenomena such as Marangoni effects, surface viscosity, the dynamics of surfactant adsorption and Gibb's elasticity should play a role. Finer phenomena should also be considered, such as the hydrodynamics involved either near a moving "contact" line between two such objects or within Plateau borders. By contrast, for solid grains, very different phenomena may come into play, including solid friction. For each of these phenomena, a simplified yet realistic pair-wise interaction law will need to be expressed and can then be included in the simulation rather easily.

D. Perspectives

The present study suggests that further investigations using the Soft-Dynamics method with larger systems (including particle rotation as well as boundary conditions) should provide interesting results, not only with the present system of plain, elastic beads in a viscous fluid, but also with different types of particle interactions. By testing ideas such as the influence of the typical interparticle gap (or other quantities if the interactions are different), it should also provide hints for analytical modelling beyond the role of the particle volume fraction.

  FIG. 2: Model of interaction for two elastic spheres in a viscous fluid. (a) elastic deflection of the surfaces (in traction); (b) normal dissipation due to Poiseuille flow in the gap; (c) tangential dissipation due to the Couette flow. The force is transmitted from a particle to another through the fluid and through a possible remote force. Such a system behaves like a Maxwell fluid (a dashpot and a spring in series). The effective friction is a function of the gap h and of the size a of the surface through which the force is transmitted.

FIG. 4 : 1 .

 41 FIG. 4: Time evolution of various quantities in the course of a T 1 process, with parameters κ = 310 -3 , h 0 = 10 -2 R and N x = 0. The gap h (a), the normal and tangential forces (c) and the ratio |δ n |/h (d) are plotted for the vertical (solid black lines, AD), the horizontal (dotted black lines, BC) and the diagonal (solid grey lines, AB etc) pairs of particles. On graph (c), the dotted grey line represents the tangential force S d of a diagonal pair such as AB (which is zero for the vertical pair AD and horizontal pair BC). Graph (b) shows the angle θ such as defined on Fig. 1. On graph (d), a pair of particles for which |δ n |/h > 1 is in the Hertz regime. If |δ n |/h < 1, it is in the Poiseuille regime.

FIG. 5 :

 5 FIG.5:(Color online) Typical duration ∆ of a T 1 process as a function of the initial gap h 0 and of the dimensionless applied force κ. The data points were obtained through the Soft-Dynamics simulation presented here. Blue open circles correspond to T 1 s where the horizontal pair has remained in the Poiseuille regime during the entire process. Full red squares correspond to T 1 s, such as that represented on Fig.4, whose horizontal (separating) pair has been in the Hertz regime for part of the time. The surface is that defined by the theoretical model for both regimes (Eqs. 33-34 with h f = 2.5R).

Acknowledgments

We gratefully acknowledge fruitful discussions with Franc ¸ois Molino and with participants of the GDR 2983 Mousses (CNRS). This work was supported by the Centre National de la Recherche Scientifique (CNRS), by the Université Paris-7 Paris Diderot, and by the Agence Nationale de la Recherche (ANR-05-BLAN-0105-01).

APPENDIX A: DYNAMICS OF PARTICLES

In this Appendix, we deduce the particle dynamics, given by Eqs. (18), from the physical model of interactions and the mechanical equilibria described in Sec. III. We start from the time derivative of the particle force balance given by Eq. ( 3):

Let us express the above as a sum of (i) -→ F ext i which is supposed to be known, (ii) terms that are linear in the particle center velocities, and (iii) another term that is explicitly known from the current state of the system, i.e., from -→ X i j and -→ δ . To this aim, using Eq. ( 12), let us express -→ F ela i j in terms of the partial derivative of -→ F ela (a, E, -→ δ ):

Combination

We can now easily express each terms of Eq. (A2) as a function of -→ X i j . The first term is directly given by ( 6):

The second term involves the time derivative Ė = E(c nc t ) α of the contact stiffness expressed by Eq. (13). Using (A13), it can be expressed as:

Note that this term vanishes for c n = c t . The third term involves the time derivative of a, which we express through its partial derivatives: ȧ(δ n , h) = a δ n δ n + a h ḣ (we used the notation a δ n = ∂ a ∂ δ n and a h = ∂ a ∂ h ). Replacing δ n and ḣ by their respective expressions in terms of -→ X i j , Eqs. (A14) and (A15) lead to:

Finally, substituting the three results of Eqs. (A3), (A6) and (A8) into Eq. (A2) yields:

where

-→ δ and -→ X i j Then, for each particle i, summing on its interacting particles j and using the force balance (3) yields the system of equation (18).

Preliminary differentiations

According to the definition of the normal vector,

, and to that of the associated projector, α = -→ n • -→ n T , we obtain their time derivatives:

as two an explicit functions of -→ δ and -→ X i j . Indeed, according to Eq. ( 16), -→ v s can be expressed as a function of the elastic force:

The evolution of the normal deflection δ n = -→ δ T

• -→ n can be expressed as a function of -→ X i j by using Eqs. (A12) and ( 6):

Finally, from Eq. ( 6), we deduce the expression of the gap evolution ḣ as a function of -→ X i j :