Pierre Rognon

Cyprien Gay
email: cyprien.gay@univ-paris-diderot.fr

Soft Dynamics simulation: normal approach of deformable particles in a viscous fluid

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

Discrete simulation methods are efficient tools to investigate the complex behaviors of complex fluids made of either dry granular materials or dilute suspensions. By contrast, materials made of soft and/or concentrated units (emulsions, foams, vesicles, dense suspensions) can exhibit both significant elastic particle deflections (Hertz-like response) and strong viscous forces (squeezed liquid). We point out that the gap between two particles is then not determined solely by the positions of their centers, but rather exhibits its own dynamics. We provide a new discrete numerical method, named Soft Dynamics, to simulate the combined dynamics of particles and contacts. As an illustration, we present the results for the approach of two particles. We recover the scaling behaviors expected in three limits: the Stokes limit for very large gaps, the Poiseuille-lubricated limit for small gaps and even smaller surface deflections, and the Hertz limit for significant surface deflections. Larger scale simulations with this new method are a promising tool for investigating the collective behaviors of many complex materials.

Foams, emulsions and granular materials are made of interacting particles, respectively bubbles, droplets and grains, in a surrounding fluid. A great deal of research revealed their 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]. How are these behaviors related to the particle interactions? This issue is central in the current debate about constitutive laws, dilatancy [START_REF] Weaire | [END_REF]4], and shear banding [5,6,7,8,9]. Discrete simulation methods are efficient tools for such an investigation, since the motion of each particle is derived from the interactions directly. Molecular Dynamics (MD) for dry granular [10], Stokesian Dynamics (SD) for viscous suspensions [11] and Bubble Model (BM) for foams [12] implement various adjustable particle interactions and provide the corresponding macroscopic behaviors, with good agreement with experiments [13,14,15]. Hence, implementing relevant interaction models is crucial. Accurate models are already available for either deformable elastic-like grains without surrounding fluid (MD for dry granular), or non-deformable particles in a viscous fluid (SD for suspension). By contrast, the combination of both particle deformation and viscous flow has not been fully described yet, although it is central in such materials as foams and emulsions. BM provided a first step in that direction. Like with MD for grains, BM particles can slightly overlap with an associated normal elastic-like repulsion. Besides, a tangential viscous force arises from the shearing of the contact film. However, BM does not include any viscous effects for normal displacement, and is therefore not fully suitable for the normal approach or separation of particles, which are essential in T1 and other topological processes.

In this paper, we address the description of normal interactions between two elastic-like particles in a viscous fluid, as depicted in Fig. 1. We show that such a system exhibits two distinct dynamics for the center-to-center and surface-to-surface distances. As this feature was ignored so far in existing simulations, we here introduce a new discrete numerical method, named Soft Dynamics, and implement the simple case of two similar particles for illustration purposes. We recover the scaling behaviors expected in Stokes, Poiseuille-lubricated and Hertz limits.

R η E h W E R X δ FIG.
1: Sketch of normal particle interaction. A force between elastic surfaces (modulus E) is transmitted partly through the fluid (visosity η) and partly through possible remote interaction (labeled W ). It can deform them (deflection δ ), which affects the flow. The consequence is that the center-to-center distance X and the gap h between both surfaces exhibit two distinct dynamics.

Model of interaction.

Let us describe a minimal model of interaction between two identical, elastic spheres with Young modulus E and radius R within a Newtonian fluid of viscosity η. We neglect both particle and fluid inertia, thereby assuming the system to be in a Stokes regime. Furthermore, we focus on normal motion. An external force F ext induces a relative motion of particles, so that both their inter-centre distance X and the gap h evolve in time. Particle surfaces are then submitted to the viscous force F vis due to the motion of the fluid in the gap, and also to a possible remote interaction F adh that is not transmitted through the fluid. Submitted to these two forces, particle deform elastically (surface normal deflection δ , counted as positive while in compression). The force balance (with no inertia) and the geometry yield the following equations of the system (Fig. 1):

F ext = F ela = F vis + F adh , (1) X = 2R + h -2δ . (2)
According to Hertz [16], the elastic deformation of a sphere due to a force acting on a region of radius a is mainly located within a volume of size a, and scales like δ /a. Hence, the total elastic energy stored in the contact scales like Eaδ 2 and the associated force reads:

F ela (a, δ ) = αEaδ . ( 3
)
where the geometrical constant α is equal to 4/9 for elastic solids. The Young modulus E ranges typically from 100 GPa for metals to 100 MPa for polymers or 1 MPa for biomaterials. By contrast, Eq. ( 3) is not so accurate for bubbles and droplets since their bulk is not elastic. Nevertheless, due to the surface tension, σ , their surface deflects in an elastic-like manner. The effective Young modulus is much smaller than for solids since it scales like σ /R and thus ranges from 1 Pa to 10 4 Pa (σ ≈ 10 -2 Nm -1 , R ≈ 1 cm to 1 µm). Although other parameters can play a role in the elastic-like properties of bubbles and droplets, such as the number of contacts [17,18], the elasticlike force mainly depends on the deflection δ and size a of the interacting region.

The viscous force exerted by the fluid on the particle surfaces, due to their normal motion ḣ, is well established in three limits. First, the Stokes limit holds for very large gaps (h ≫ R): the viscous force 6πηR ḣ is transmitted through a region of size a ≈ R. Secondly, for smaller gaps (h R), and negligible surface deflections (|δ | ≪ h), the interaction can be regarded as a Poiseuille lubrication flow between two spheres. Hence, the viscous force 6πηR 2 h ḣ is mainly transmitted through a region of size a where the surface-to-surface distance lies between h and 2h, yielding a ≈ √ 2Rh (Fig. 2 a). Thirdly, the Hertz limit corresponds to elastic deflections larger than the gap (|δ | ≫ h). In this case, the force is transmitted mainly through a region of size a determined by Hertz theory, a ≈ 2R|δ |, and the contact can be regarded as two parallel disks, so that the viscous force reads 3πηa 4 2h 3 ḣ. In order to interpolate between these three limits, we express the viscous force as:

F vis (a, h, ḣ) = -ζ ḣ, (4)
where

ζ (a, h) = 3πηa 4 2h 3
is the effective normal friction coefficient of the contact, with h = Rh/(R + h). Other choices are possible for this coefficient, but it will at least depend on the contact geometry through h and a. We choose an expression for the size of the interacting region that interpolates between the asymptotic expressions valid in all three regimes (Stokes, Poiseuille and Hertz):

a(h, δ ) = 2R(h + |δ |) (5)
Other choices for a are possible, but a will depend at least on h and δ .

R -δ a R R a (a) (b) 2h
FIG. 2: Scheme of the region through which both particles interact (bold lines). (a) In the Poiseuille lubricated limit for sligthly deflected surfaces (δ ≪ h), the size a of the interacting region depends mainly on h and marginally on δ . (b) By Contrast, in the Hertz limit for significant deflections (δ ≫ h), it depends mainly on δ and marginally on h.

A remote interaction between both particle surfaces can also be present. It may include steric repulsion (with a closest approach distance related to the size of the surface asperities in solid grains), surfactant repulsion or, more generally, disjoining pressure in foams or emulsions [19]. It may also include some adhesion. Various models can be developped for such remote forces, which will at least depend on the gap h and on the size a of the interacting region: F adh (a, h). In the present work, for simplicity, we omit such an interaction.

Let us now consider two particles with elastic deflection and viscous interactions. This situation evokes a Mawxell fluid whose elastic and viscous coefficients would depend on the gap h and on the size a of the interacting region. We can derive the expected scaling behaviors of δ , a and ḣ with respect to h in several asymptotic regimes, and express them simply (see Fig. 3 a,b) in terms of a typical Stokes velocity V and a renormalized applied force κ defined by:

V = F ext /(ηR), κ = F ext /(ER 2 ). (6)
We do not address large elastic deformation (δ ≃ R) that corresponds to κ > 1. In the Stokes regime (h ≫ R) the deflection δ can be neglected in a first approximation, and we find ḣ ≃ V and a ≃ R. According to Eq. ( 3), we then obtain the deflection δ ≃ κ. The crossover to the Poiseuillelubricated regime occurs when the gap h becomes comparable to R, while remaining much larger than δ . Again, in the Poiseuille regime, δ can be neglected at first, and we have a ≃ √ Rh. Thus, from Eqs. ( 4) and ( 6), ḣ ≃ V h/R. Then, from Eqs. (3) and ( 6), we get δ ≃ (R 3 /h)κ. The crossover between Poiseuille and Hertz regimes occurs when the elastic deflection becomes comparable to the gap δ ≃ h, which corresponds to κ ≃ (h/R) 3/2 . In the Hertz regime, a ≃ √ Rδ . According to Eqs. ( 3) and ( 6), this implies a ≃ Rκ 1/3 . Finally, using Eqs. ( 4) and ( 6), we obtain ḣ ≃ V κ -4/3 (h/R) 3 . 10 -4 10 -6

1 0 κ ≃ ( h / R ) 3 / 2 κ ≃ 1 h/R ≃ 1 10 -2 10 -4 1 10 -2 10 -4 10 -6 1 10 2 10 -4 κ 1 κ 1/3 κ 2/3 R h κ h R h R 1 a/R δ /R ḣ/V h R 3 κ -4 3 (e) (f) (d) (b) (a) h R ḣ V κ h R κ h R h R κ λ 1 λ 1 λ 2 λ 2 Poiseuille (P) Hertz (H) Stokes (S) h R δ R (c) κ κ 10 -6
FIG. 3: Scaling behaviors of two elastic-like sphere in a viscous fluid. (a) Scaling expected from theory in the Stokes, the Poiseuille-lubricated and the Hertz limits; the elastic deflection δ /R, the contact size a/R and the gap rate of change ḣ/V are function of the current gap h/R and renoralized force κ. (b) Phase-diagram expected from theory for the crossover (large lines) between the three limits ; lines at constant κ come from Soft Dynamics simulations (see Fig. 4). (c-f) δ (h), ḣ(h) and their power law also from these simulations. Theoretical phase-diagram is replotted on the 3D plot basis.

Soft Dynamics simulation.

Let us now show that the gap h is not only dictated by the center-to-center position X, but exhibits its own dynamics. The system evolution is determined by Eqs. ( 1) and ( 2), and by the expressions for the forces F ela (a, δ ), F adh (a, h), F vis (a, h, ḣ) and a(h, δ ). The viscous force (Eq. 4) and force balance (Eq. 1) directly yield the gap evolution:

ḣ = - F ext -F adh ζ . ( 7
)
This is a purely viscous behavior with an effective friction ζ (a, h) and an applied force F ext -F adh which corresponds to the part of the total force that is transmitted through the fluid.

Obtaining the dynamics of X is also straightforward, although it requires more steps. To this aim, it proves helpful to express time derivative of a and F ela in terms of their partial derivatives (denoted by a subscript, for example F ela a = ∂ F ela /∂ a):

ȧ = a h ḣ + a δ δ , Ḟela = F ela a ȧ + F ela δ δ . ( 8
)
Substituting ȧ and 2 δ = ḣ -Ẋ (from Eq. 2) into the expression of Ḟela yields Ḟela = -C Ẋ + ḣ(C + F ela a a h ), where C = (F ela a a δ + F ela δ )/2. Then, using the force balance Ḟela = Ḟext and the expression for ḣ (Eq. 7) leads to the dynamics equation for X:

Ẋ = - Ḟext C - F ext -F adh D , (9)
with D = C ζ /(C + F ela a a h ). The above expression evokes a Maxwell behavior, involving the time derivative of the external force, with stiffness coefficient C and friction coefficient D. It thus differs from the purely viscous behavior of the gap by itself (Eq. 7). Both dynamics of X and h are here expressed in a generic and explicit manner in terms of the partial derivatives F ela a,δ and a h,δ and the effective friction coefficient ζ , and also a possible remote force F adh . The expressions for F ela (a, δ ), F adh (a, h) and a(h, δ ) can therefore be changed at will to describe various physical contact behaviours. The only restriction is that the viscous force must be linear in ḣ, consistently with the no-inertia assumption.

Knowing the distinct dynamics of h and X, we can now propose a method to simulate their time evolution in the presence of a known external force F ext (t) (see Fig. 4). Initial values of X and h yield the deflection δ according to Eq. ( 2), and then the size of the contact a(h, δ ). Quantities ζ , C and D are explicitely given by the expression of the forces and a. Finally, the external force F ext and its time derivative are needed to obtain both ḣ and Ẋ. By integration, new values for h and X are obtained.

We implemented the Soft Dynamics method in the simple case of two identical particles. They are initially separated by a large distance (h/R = 100), we then subject them to a constant, compressive force F ext , and the gap decreases continually. h may even tend to zero asymptotically since in this case, for simplicity, we have not included the remote force F adh . For various applied forces (10 -6 < κ < 0.3), Fig. 3b shows the range of gaps h reached throughout the approach. Plotting δ and ḣ as a function of the current value of the gap 4: Method for the Soft Dynamics simulation. From current center-to-center distance X and gap h, their respective time evolutions are explicitely calculated using the value of the external force F ext and its time derivative. The intermediate calculations involve analytical expressions that result from the physical description of the system (solid rectangular boxes). The dashed rectangular box (with C and D) represents intermediate calculation steps.

F ela F ela δ ,a F adh ζ Ẋ ḣ a a δ ,h C D F ext Ḟext δ X h FIG.
h allows us to extract the corresponding scaling behaviors. Figs. 3c,d represent function δ (h) and the associated exponent λ 1 (h) such that one locally has δ (h) ∝ h λ 1 . In both the Hertz and the Stokes regimes, the deflection δ depends only on the contact force, with the usual 2/3-power law in the Hertz regime and a linear law in the Stokes regime where the interacting region spans over one entire side of the particle. In the Poiseuille-lubricated regime, the width of the interacting region depends on the gap h: it causes δ to vary like (h/R) -1/2 . Function ḣ(h) is represented on Figs. 3e,f together with the associated local exponent: ḣ(h) ∝ h λ 2 . In the Stokes regime, the large gap h ≫ R implies that the interaction, mediated by the surrounding liquid, does not depend on distance: ḣ is then independent of h. By contrast, in both the Poiseuille and the Hertz regimes, the gap thickness affects the liquid flow, and ḣ varies respectively like ḣ ∝ h 1 and ḣ ∝ h 3 . Note that the range of gaps h where λ 1 ≈ -1/2 and λ 2 ≈ 1 shrinks and vanishes as the renormalized force κ is increased: strongly compressed particles cross over from Stokes to Hertz regimes directly as h deacreases. All these observations are consistent with the theoretical predictions summarized on Fig. 3a.

The conclusion of this work is that particles interacting in a viscous fluid can exhibit two distinct dynamics, which govern their overall motion on one hand and the gap evolution on the other hand. This behavior results from the interplay between the underlying Poiseuille and Hertz mechanisms: the viscous flow in the gap affects the transmitted force and hence the surface deflection, which in turn affects the flow in the gap. In Stokes and Poiseuille regimes, grain surfaces are remote and the deflection can be neglected: the gap dynamics is mainly dictated by the motions of the particle centers. However, as soon as the deflection becomes larger than the gap, both dynamics must be distinguished. This is easily encountered with soft particles, but can also occur with rather hard particles as soon as the distance h becomes small enough. It seems that this feature had not been implemented so far: we here provided a new simulation method, called Soft Dynamics, which expresses each of these two dynamics generically in terms of the forces involved (elastic, viscous as well as possible remote interaction). Here, we merely check that it describes the normal approach of two particles in a viscous fluid correctly. But Soft Dynamics can be implemented for many multicontacting particles [START_REF]Soft Dynamics" approach[END_REF]. Hence, like MD for dry grains and SD for diluted suspensions, Soft Dynamics should be able to provide new insights into the behavior of materials made of soft, closepacked units. Bubbly liquids, vesicules, foams, emulsions, micellar solutions and dense suspensions could thus be investigated, by including the specificities of their interactions.

We gratefully acknowledge fruitful discussions with Franc ¸ois Molino and with participants of the GDR 2352 Mousses (CNRS). This work was supported by the Agence Nationale de la Recherche (ANR05). works with several multi-contacting particles. Eqs. ( 7) and ( 9) now apply to each pair i j of contacting particles, with centerto-center vectors X i j = X j -X i , gaps h i j , and coefficients that vary from pair to pair. For each particle i, let us use the force balance equation F ext i = ∑ j F ela i j (see Eq. 1) and sum up the new Eqs. ( 9) that correspond to neighbouring particles j. We thus obtain the following system of vector equations (one equation per particle):

∑ j C i j Ẋj -Ẋi ∑ j C i j = -Ḟext i + ∑ j C i j D i j
F ela i j -F adh i j

Inverting this system of linear equations yields the particle center velocities Ẋi (in fact, these equations are not independent: one of them must be replaced, for instance, by the condition that the average particle velocity is zero). Using the Eqs. ( 7) then provides the gap evolutions. Thus, when several particles are present, the very last step of the method pictured on Fig. 4 now involves the inversion of a system of equations.