R Chabreyrie

D Vainchtein

C Chandre

P Singh

N Aubry

Tailored mixing inside a translating droplet

Keywords: numbers: 47.51.+a, 47.61.Ne, 47.52.+j

Tailored mixing inside individual droplets could be useful to ensure that reactions within microscopic discrete fluid volumes, which are used as microreactors in "digital microfluidic" applications, take place in a controlled fashion. In this article we consider a translating spherical liquid drop to which we impose a time periodic rigid-body rotation. Such a rotation not only induces mixing via chaotic advection, which operates through the stretching and folding of material lines, but also offers the possibility of tuning the mixing by controlling the location and size of the mixing region. Tuned mixing is achieved by judiciously adjusting the amplitude and frequency of the rotation, which are determined by using a resonance condition and following the evolution of adiabatic invariants. As the size of the mixing region is increased, complete mixing within the drop is obtained.

I. INTRODUCTION

Droplets have been proposed as an alternative to standard fluid-stream microfluidics for lab-on-a-chip applications. This microfluidics approach, also referred to as "digital" because it uses "discrete" fluid volumes (droplets) rather than continuous streams, holds great promise due to the possibility of using single droplets as microreactors [1]. Efficient mixing, however, is needed for reactions to occur, but remains difficult to achieve because the Reynolds number (Re) is usually very small and so the flow is laminar. This issue has recently attracted much attention in the literature. For flows in microchannels, while there are many strategies based on altering the channel geometry, the use of forcing alone (see, e.g., [2,3,4,5,6]) has also proved to be efficient, especially in the case of low Re [7]. The combination of both geometry alteration and forcing has been explored as well [7,8,9,10]. For droplet-based microfluidics, the forcing alone is the preferred strategy as the deformation of the droplet is difficult to control. In almost all cases, the enhancement of mixing in miniature geometries is based on chaotic advection, the stirring phenomenon that stretches and folds fluid elements thus increasing the interfacial area between the two fluids to be mixed. Chaotic advection inside a liquid drop subjected to a forcing (at low Re) has been studied extensively [11,12,13,14,15,16,17] and was obtained experimentally by means of oscillatory flows [18,19]. In this letter, we focus on unsteady -yet periodic -forcing. From a dynamical systems viewpoint, the introduction of a time-dependent perturbation or forcing breaks the invariants (related to the symmetries of the unperturbed system), thus introducing resonances between the natural frequencies of the unperturbed problem and the frequency(-ies) of the forcing. Although such resonances create chaotic regions where mixing occurs, in general, chaotic and regular regions co-exist and unexpected regular sizable pockets persist. In many situations where it is indeed possible to create chaos, controlling the mixing region(s) remains a challenge. Such a control, however, should be possible since a chaotic system is sensitive to changes in parameter values (as it is to changes in initial conditions). These changes should generically modify the resonances, and thus the location and size of the chaotic regions. Our general approach along these lines is to consider a bounded three-dimensional (3D) flow, which is the superposition of an integrable flow v 0 with at least one invariant and a small time-dependent perturbation εv 1 (x, t), 0 ≤ ε ≪ 1. If v 0 has only one invariant, the phase space contains two-dimensional tori. In this case, the perturbed flow, ẋ = v 0 (x) + εv 1 (x, t), has poor mixing properties if the amplitude of the perturbation ε is small, since twodimensional (2D) tori act as barriers to chaotic diffusion (e.g., [20]). If, on the other hand, v 0 has two invariants, trajectories of this integrable flow are all periodic. Most of these periodic orbits are expected to be broken by a generic perturbation v 1 with an arbitrarily small amplitude ε. Efficient mixing properties might then be obtained with such perturbed flows. In this work, we consider an axisymmetric integrable flow possessing two invariants, thus possibly offering efficient mixing properties after being perturbed. While many previous works [11,12,21,22] have shown the existence of chaotic behavior in 3D bounded steady flows, we turn our attention to unsteady flows; the added unsteadiness targets the control of the chaotic behavior through resonance phenomena [17,23,24]. Specifically, we seek to create a mixing zone of tunable size which remains localized within a well-defined region of the drop. This should also provide a rationale for the route to complete mixing as the perturbation increases.

II. MODEL A. Flow equation and assumptions

We consider a spherical Newtonian drop immersed in an incompressible Newtonian flow in the case where the linear external field is characterized by translational velocity and vorticity vectors, similarly to [12]. As in the latter reference, we assume that the local Re is much smaller than one and that the interfacial tension is sufficiently large for the drop to remain spherical. The internal velocity field is obtained by solving the Stokes flow problem for both the internal and external flows satisfying the continuity of velocity and tangential stress conditions across the drop surface. In addition, we introduce unsteadiness in the problem by making the vorticity time dependent. In a Cartesian coordinate system translating with the center-of-mass velocity of the drop, and with the z axis in the direction of the translation, the paths of passive marker particles are given by the solution of the non-autonomous dynamical system:

u = ẋ = zx -a(t)ω z y, v = ẏ = zy + a(t) (ω z x -ω x z) , (1)
w = ż = (1 -2x 2 -2y 2 -z 2 ) + a(t)ω x y,
where all lengths and velocities have been nondimensionalized by the drop radius and the magnitude of the translational velocity. Here, the vorticity is de-

fined by ω = (ω x , ω y , ω z ) = (1/ √ 2, 0, 1/ √ 2)
, the unitary vector corresponding to the axis of rotation, and a(t) = ε/2 (1 + cos ωt), characterized by the frequency ω and the amplitude ε. In this letter, we consider only small amplitudes, i.e. for 0 ≤ ε ≪ 1. Note that the former equations are identical to those in [12] except that the constant vorticity vector has been replaced by a(t)ω. This can be done by either assuming unsteady vorticity in the external flow field, or by applying a time dependent body force. In practice, this could be realized, e.g., by creating a time dependent swirl motion in the external flow or by applying an electric field that exerts a torque on the drop (e.g., [25] it or work on electrorotation). This flow is the superposition of a Hill's vortex and an unsteady rigid body rotation, and the surface of the drop,

r 2 = x 2 + y 2 + z 2 = 1, is invariant under flow (1).

B. Integrable case

We now discuss some features of the unperturbed axisymmetric (2D) flow (ε = 0). The flow possesses two independent integrals of motion, e.g., the streamfunction ψ and the azimuthal angle φ:

ψ = 1/2ρ 2 1 -r 2 , φ = arctan y/x, (2)
where ρ 2 = x 2 + y 2 and ψ ∈ [0, 1/8]. The streamlines of the unperturbed system are lines of constant ψ and φ, denoted by Γ ψ,φ , and defined as (1 -2ρ 2 ) 2 + (2ρz) 2 = 1 -8ψ (see Fig. 1). Almost all streamlines are closed curves surrounding a circle of degenerate elliptic fixed points (ρ = 1/ √ 2, z = 0). In addition, there are two hyperbolic fixed points located at the poles of the sphere which are connected by heteroclinic orbits. The frequency of the motion on Γ ψ,φ is given by 2π

Ω(ψ) = π/2 -π/2 √ 2 dα 1 + γ(ψ) sin α = 2 √ 2 √ 1 + γ K 2γ 1 + γ , (3)
where γ(ψ) = √ 1 -8ψ and K is the complete elliptic function of the first kind. The frequency Ω is bounded by two limits, Ω(0) = 0 and Ω(1/8) = √ 2 (see Fig. 1). On every streamline Γ ψ,φ , we introduce a uniform phase χ mod(2π) such that χ = 0 on the x -y plane (with ρ ≤ 1/ √ 2) and χ = Ω (ψ). The unperturbed system, which can be rewritten in terms of (ψ, φ, χ) as ψ = 0, φ = 0, χ = Ω(ψ), belongs to the class of action-action-angle flows.

C. Perturbed case

In the perturbed case 0 < ε ≪ 1, the time evolution of the two invariants of the unperturbed system is given by ψ = -2a(t)ω x ψ sin φG (ψ, χ) ,

φ = a(t)ω z -a(t)ω x cos φG (ψ, χ) , (4)
where G(ψ, χ) = z/ρ is 2π periodic in χ and has zero average in χ. The time evolution equation for χ is

χ = Ω(ψ) + a(t)H(ψ, φ, χ),
where H is 2π periodic in χ. The dynamics possesses two time scales, a fast one (of order one) associated with χ, and a slow one (of order 1/ε) associated with ψ and φ.

If Ω and ω are incommensurate, then the averaging over Ω and over ω can be performed independently. In this case, the time-periodic terms in Eq. ( 4) average out, and the averaged system reduces to ψ = 0, φ = -ε/2. Thus in the averaged system the value of ψ is conserved as it was in the unperturbed system; in other words, ψ is an invariant of the averaged system. Each trajectory of the averaged system evolves on two-dimensional nested tori T ψ . In the perturbed system, ψ is an adiabatic invariant and the motion follows adiabatically the tori T ψ .

III. METHODS AND RESULTS

A. Mixing generation via resonance phenomena

We now turn to the generation of a 3D chaotic mixing region inside the drop, for which we seek to control both the location and the size. The strategy used for this purpose is to bring a chosen family of unperturbed tori T ψ into resonance with the perturbation a(t) by adjusting the frequency ω to satisfy the resonance condition

nΩ(ψ) -ω = 0, ( 5
)
for some n ∈ N (see Fig. 1). For any fixed ω we denote by T (n) (ω) |n ∈ N the set of resonant tori T ψ satisfying (5). Hereafter, we denote the chaotic mixing region generated around T (1) (ω) by CMR.

B. Control of the mixing

Figures 2 and 3 present Liouvillian sections of the perturbed system, which consist of 2D projections of timeperiodic 3D flows by a combination of a stroboscopic map and a Poincaré section (here, the y = 0 plane). Figure 2 shows that a perturbation a(t) creates a 3D CMR around T (1) (ω) and its location is controlled by varying ω according to Eq. (5). In what follows, we analyze the location and the size of the CMR as ω and ε vary.

For small values of ω, all resonances are located near the pole-to-pole heteroclinic connections (at ψ = 0, near the z axis and near the boundaries of the drop, see Fig. 2a). As ω is increased, the CMR penetrates deeper into the drop (Fig. 2b). In the interval 0 < ω < √ 2, the CMR is the largest chaotic region (compared to chaotic regions corresponding to higher order resonances), with all the other chaotic regions localized close to the z axis and near the drop boundaries (around the heteroclinic orbits); this is due to the shape of Ω (ψ). As ω is increased further, the CMR moves toward the location of the elliptic fixed points of the unperturbed system, closely following the location of the resonant torus T (1) (ω) (Fig. 2c). As the value of ω approaches √ 2, the CMR shrinks to the circle of elliptic fixed points (Fig. 5).

Whereas the frequency ω of the rigid body rotation is mostly responsible for the location of the CMR, it is its amplitude ε which mostly determines its size. Figure 3 shows that the size of the chaotic mixing regions created by the n = 1 resonance and by higher order resonances (mostly the n = 2 resonance) increases as the amplitude of the perturbation increases. Around ε ≈ 0.20, the chaotic regions around the heteroclinic orbits and the CMR join together to cover the entire drop volume.

Recall that in the averaged system the adiabatic invariant ψ is constant. In the exact system, however, along a given trajectory starting at ψ = ψ 0 it varies between ψ -(ψ 0 ; ω, ε) and ψ + (ψ 0 ; ω, ε). The width ∆ψ = ψ + (ψ 0 ; ω, ε) -ψ -(ψ 0 ; ω, ε) is small away from the resonance, and increases significantly closer to the resonance. The projection of three characteristic trajectories onto the (ψ, φ)-plane (called the slow plane in dynamical sys- ε for ω = 0.55, 0.93, 1.28, 1.41 (a-d).

tems) is presented in Fig. 4. The narrow regions on the sides are off-resonance trajectories that stay quite close to the corresponding tori T ψ . In between, the middle trajectory deviates much further from its T ψ = T (1) (ω) and fills the entire CMR. The quantity ∆ψ is probably the most convenient quantity to estimate the size of the CMR (around T (1) (ω) for ω < √ 2). The volume between the tori T ψ -(ω) and T ψ + (ω) gives the CMR size in 3D. The dependence of the size of the CMR (in terms of ∆ψ) on ε and ω is illustrated in Fig. 5. The curves (a)-(d) in the upper and lower panels correspond to the Figs. 2a-d and 3a-d, respectively. For a given ω value (i.e. for a given T (1) (ω)), the size can be controlled by adjusting the value of ε; for example in the range of frequencies 1.181 ≤ ω ≤ 1.357, the entire droplet exhibits chaotic mixing for ε ≥ 0.175. For each smaller value of ε the size reaches a maximum for a certain value ω m (ε) of the frequency. On the one hand, this property can be used as an optimization technique to obtain the maxi-mal CMR size one can reach for a given amplitude ε of the rotation. On the other hand, ∆ψ versus ε increases quite monotonically for all values of ω. The derivation of the maxima locations and estimates of ∆ψ as functions of the parameters and the order of resonance, will be addressed elsewhere. The structure of the CMR in our case is rather different from that obtained in other problems that possess resonance-induced chaotic advection. Namely, here the size of the CMR vanishes as ε goes to 0 and the CMR is localized near the resonance. In contrast, in the flow considered in, e.g., [17], the mixing is caused by resonances, but the CMR occupies a volume on the scale of the whole system. The difference comes from the fact that the averaged change of the frequency of the fast system vanishes in the current system, thus preventing the trajectories starting away from the resonance from approaching it. This property makes the kind of flows investigated here useful as it may be advantageous to localize the mixing in certain parts of the system only.

IV. CONCLUSION

In summary, we have shown that by applying a judicious oscillatory rotation to a translating drop (an integrable system), one can create a chaotic mixing zone with a prescribed location and size. The appropriate values of the parameters of the perturbation (here, a rotation of a given frequency and amplitude) are determined by quantitative features of the integrable case. For any amplitude of the rotation, the frequency optimizing the CMR size has been obtained. Such an optimization could be useful in guiding the design of practical mixing devices aiming at the best possible mixing rate within individual drops.

FIG. 1 :

 1 FIG.1: Streamlines inside the drop (without rotation) and their frequencies Ω (ψ) as given by Eq. (3).

FIG. 2 :

 2 FIG.2: Liouvillian sections for the amplitude ε = 0.03 and the frequencies ω = 0.55, 0.93, 1.28, 1.41 (a-d). The (red) dashed line inside the CMR is the torus T(1) .

FIG. 4 :FIG. 5 :

 45 FIG. 4: Projection of three characteristic trajectories on the slow phase plane, with φ0 = 0 and ψ0 = 0.010, 0.073, 0.125.

Acknowledgments

This article is based upon work partially supported by the NSF (grants CTS-0626070 (N.A.), CTS-0626123 (P.S.) and 0400370 (D.V.)). D.V. is grateful to the RBRF (grant 06-01-00117) and to the Donors of the ACS Petroleum Research Fund. C.C. acknowledges support from Euratom-CEA (contract EUR 344-88-1 FUA F).


G Kermouche

E Barthel

D Vandembroucq

Ph Dubujet

Mechanical modelling of indentation-induced densification in amorphous silica

Keywords: Pressure dependent plasticity, densification maps, Nanoindentation, Finite element method

Though brittle, silicate glasses undergo plastic deformation at the microscopic scale. Mechanical contact and indentation are the most common situations of interest. The plasticity of glasses is characterized not only by shear flow but also by a permanent densification process. After indentation, densification can locally reach 20% in a pure silica glass. In this paper, a new constitutive model, derived from experimental observations is presented to account for the plasticity of fused silica. The use of nanoindentation tests to identify the plastic behavior of amorphous silica is discussed. A set of material properties is determined by comparing experimental load/displacement indentation curves to the results of finite element simulations. The numerical results show good agreement with recent experimental indentationinduced densification maps obtained by Perriot et al. [

Introduction

Though primarily brittle [START_REF] Zarzycki | Glasses and the vitreous state[END_REF], silicate glasses can undergo plastic deformation at the micron scale. his can occur under contact loading [START_REF] Ernsberger | [END_REF] or in the close vicinity of a crack tip. In that respect glasses do not differ from standard materials, except that the characteristic length scale below which volume deformation (plasticity) is energetically more favorable than surface creation (crack propagation) is smaller [3]. Although known for half a century [4] and of primary interest for many micron-scale applications, the plasticity of silicate glasses has up to now remained little studied [5][6][7][8].

Standard dislocation-based mechanisms of crystalline solids is of no direct use to explain the plasticity of silicate glasses because they are amorphous.

Instead, according to the dominant scenario, plastic deformation of glasses seems to result from a succession of localized structural rearrangements [9,10].

On the plasticity of amorphous materials and more specifically of metallic glasses, see for instance the recent review of Schuh [11].

Moreover, plasticity of silicate glasses is characterized by one more specific feature. Indeed, it appears that a more or less significant phenomenon of permanent densification [12] takes place in addition to the isochoric plastic strain classically observed in metallic materials. This characteristic leads to a distinction between normal glasses and anomalous glasses: normal glasses exhibit little densification while anomalous glasses exhibit large densification as well as a different indentation fracture pattern [13]. Fused silica, the most anomalous silicate glass of industrial interest, can achieve a significant 20% densification [5].

The lack of studies on this subject is directly related to the brittle nature of the material and the necessity to resort to mechanical measurements at a micron scale. Progress has been obtained thanks to the recent developments in instrumented nanoindentation techniques [START_REF] Loubet | Microindentation techniques in Material Science and Engineering[END_REF][START_REF] Oliver | [END_REF]. Nanoindentation tests are micromechanical tests where the load applied to a tip is continuously measured as a function of the penetration depth. Recently it has been successfully used to extract material properties such as hardness or elastic modulus of very thin devices [16][17][18]. Nevertheless, a single load-penetration curve does not provide a unique set of material properties for even moderately complex materials such as von Mises elastoplastic solids with power-law isotropic strain hardening [19].

The plasticity of amorphous silica is much more involved and an attempt to describe its behavior requires additional data to identify a constitutive law with some confidence.

Recently, Perriot et al. [7] provided a map of the local indentation-induced densification in a fused silica sample. The latter was obtained from Raman microspectroscopic mappings [20] on cross-sections of indents. For that purpose, they used the work of Sugiura et al. [21] which correlates the position of the D 2 line [22] of the Raman Spectrum of silica with the ratio of the sample density to its initial density. They also provided clear evidence for a densification-induced hardening process with diamond-anvil cell (DAC) measurement [23]. Ji et al. [8] have also evidenced such a densification-induced hardening in window glass.

These experimental results provide a rich and complex set of data from which to identify a constitutive equation for silica. In this paper a new plastic yield criterion is proposed and its implementation in a finite element software is explained in details. Then, the identification of the material parameters is carried out using instrumented indentation results. Finally the relevance of the new criterion is assessed by comparison with the experimental density maps [7,[START_REF] Perriot | Nanoindentation de couches minces deposées sur substrat de verre de silice[END_REF].

Constitutive modelling of amorphous silica

In standard metallic materials, plasticity is pressure-independent [START_REF] Ashby | Materials[END_REF]. Hence, for such materials, the yield function is often written :

f (σ ij ) = q -Y (1)
where σ ij denotes the Cauchy stress tensor, q is the von Mises equivalent stress and Y is the yield stress. The macroscopic von Mises stress is proportional to the shear stress [START_REF] Holzapfel | Non linear solid mechanics -a continuum approach for engineering[END_REF] and is given by :

q = 3 2 S ij S ij (2)
In equation ( 2), S ij is the stress deviator defined by

S ij = σ ij + pδ ij
where p = -σ kk /3 is the hydrostatic pressure and δ ij is the Kronecker tensor.

In associative plasticity the plastic rate of deformation is normal to the yield surface [START_REF] Hill | The mathematical theory of plasticity[END_REF]. This leads to :

εp ij = λ ∂f ∂σ ij ( 3
)
where λ is the plastic multiplier which can be determined from the consistency condition [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF], f (σ ij ) = 0 during plastic yielding.

For obvious reasons, such a criterion cannot account for the silica glass densification under pure hydrostatic loading and a relation between f (σ ij ) and the hydrostatic pressure p has to be taken into account instead. Recently, Xin et al. [START_REF] Shorey | Deformation of fused silica : Nanoindentation and densification[END_REF][START_REF] Xin | Densification of fused silica : Effects on nanoindentation[END_REF] proposed the following criterion for amorphous silica :

f (σ ij ) = αq + (1 -α)p -β (4)
This pressure-dependent yield criterion resembles the Mohr-Coulomb criterion initially developed for frictional materials [START_REF] Khan | Continuum Theory of Plasticity[END_REF] and involves two physical characteristics: i) the yield stress β; ii) the densification factor α which balances the respective influence of applied hydrostatic pressure, p, and von Mises stress, q. Under the assumption of an associated perfect-plastic behavior, Xin and Lambropoulos [START_REF] Xin | Densification of fused silica : Effects on nanoindentation[END_REF] identified these two parameters by comparing experimental load/displacement indentation curves to finite element simulation results.

With the values α= 0.6 and β = 9.4 GPa, they adequately reproduced their experimental data.

Recently, Perriot et al. [7] used Raman micro-spectroscopy to estimate indentationinduced densification in silica. The micrometric spatial resolution of this technique allowed them to perform a mapping of the permanent densification below an indentation imprint. It was found that the densification increases gradually from 0 to 20% when approaching the contact area while the Kim model, however, predicts a sharp transition between a fully densified hemispherical inclusion and a non densified bulk. Hence the simple criterion (4) needs to be refined to account for these additional experimental results.

Two comments regarding the limitations of such a Mohr-Coulomb like criterion should be made at this stage. First, the linearity of (4) together with the hypothesis of associated flow rule, implies that the direction of the plastic strain rate does not depend on pressure p and shear stress q. Moreover, no hardening is taken into account in this constitutive equation. Even if the relevance of hardening in metallic glasses is still a matter of discussion [11],

the existence of a hardening effect has been clearly established in the case of silica [7]. To overcome these limitations, we propose in the following a new constitutive equation. We focus first on the yield criterion and then discuss how the hardening can be incorporated.

A new yield criterion for fused silica

A strong dependence on pressure is the main characteristic of the plastic behavior of fused silica. For that reason silica may appear as closer to powders or geomaterials than to metals or crystalline materials. Yield functions designed to account for the compaction of metal powders [START_REF] Brown | [END_REF]33] thus exhibit similar dependences with hydrostatic pressure. In the present case, however, there is no evidence of frictional effect and no reason to depart from associative plasticity. Moreover, in the absence of indications of a pressure effect in tension we will ignore it for negative hydrostatic pressure.

To account for the plasticity of fused silica we propose, for positive pressures, an elliptic criterion (see Fig. 1) similar to the one used by Shima and Oyane [34] to reproduce experiments on copper. For negative pressure we assume a simple von Mises criterion. In compression (p > 0), the proposed yield criterion is:

f (σ ij ) = q q c 2 + p p c 2 -1 (5)
while in tension (p < 0), one has:

f (σ ij ) = q -q c ( 6
)
where p c and q c represent the hydrostatic plastic limit in pure hydrostatic state and the shear limit in pure deviatoric state respectively. From equations

(3) and ( 5), the plastic strain rate (for positive pressure) is then given by :

εp ij = λ 3S ij - 2q 2 c 3p 2 c pδ ij
One observes from (5) that the direction of the plastic strain rate now strongly depends on the hydrostatic pressure p. For a pure deviatoric stress state, no plastic densification appearsthe plastic strain rate is isochoricwhereas for a pure hydrostatic stress only plastic densification can occur.

Finally, in the absence of experimental data on a potential shear hardening effect, we shall only consider density hardening: as sketched in Fig. 1, p c increases with plastic densification whereas q c remains constant. The density hardening behavior of silica is addressed in the next section.

Densification-induced hydrostatic hardening

The densification-induced hardening has been measured by Perriot et al. [START_REF] Perriot | Nanoindentation de couches minces deposées sur substrat de verre de silice[END_REF] with the help of a Diamond-Anvil Cell apparatus [23,35]. Amorphous silica samples have been submitted to successive loading/unloading cycles of growing maximum pressure P max . At the end of each cycle, a Raman spectrum was measured and the densification was quantified from the shift of the D2 peak

[21] [7]. The hardening curve obtained by this procedure is presented in Figure 2.

Note that densification, i.e. the relative volume change can be related to the first invariant of the small strain tensor ε m which is called mean strain in the present study:

∆V V = ε kk = 3ε m ,
where V denotes the initial volume. It should also be noted that the small strain assumption has been implicitly applied to describe the densification in fused silica. This hypothesis is expected to be valid for silica since saturation of the densification is reached around 20 % [5].

Strain hardening phenomena occur when materials exhibit plastic deformation. Therefore, the densification-induced hardening should depend on the plastic densification, which leads us to write :

p c = g(ε p m ) ( 7
)
where g is a function which describes the evolution of p c with plastic densification. In this paper, it has been assumed that the hardening is isotropic and that the hydrostatic hardening does not depend on the deviatoric plastic strain. The mean plastic strain is computed as follows :

ε p m = t 0 | εp m |dτ
In view of its implementation in a finite element software, the function g has been fitted with a linear relation as displayed in Figure 2. Hence, Eq. ( 7) reduces to :

p c = ζε p m + p c0 (8)
where p c0 = 11.5 GPa is the hydrostatic pressure for which the first plastic densification occurs and ζ = 100 GPa is the hardening slope. Let us note that this value is only slightly smaller than the bulk elastic modulus, which relates the mean strain to the mean stress in linear elasticity, and is about K = 120 GPa for amorphous silica.

For classic metallic solids, where hardening mainly results from dislocation pinning or entanglement, the shear hardening slope is commonly one hundred times lower than the elastic modulus. In the case of amorphous materials, where plasticity is believed to originate from successive local reorganizations [11], there is no equivalent mechanism. A statistical mechanism based on the progressive exhaustion of the most unstable regions has been recently

proposed [36]. In the present case of hardening, the progressive approach to the limit density for amorphous silica, equal to the density of the most dense crystalline analog must be taken into account. In a recent study of the plasticity of window glass which is known to exhibit low densification, Ji et al. [8] used a multi-anvil apparatus to measure the bulk plastic modulus for densification of window glass. They obtained a value close to 500 GPa, i.e. four times higher than the bulk elastic modulus. To be consistent within the elastic-plastic formalism, this result requires a very large increase in the bulk elastic modulus with plastic densification. This also points to some of the difficulties which underlie the measurement and the analysis of the hydrostatic hardening of glasses.

Note finally that Perriot et al. [7] have observed that the zone where the densification saturates is small compared to the contact area. Hence the saturation process will not be really taken into account in this paper.

3 Finite element modelling of indentation-induced densification in fused silica

Numerical implementation

The yield criterion developed in the previous section has been implemented in the finite element software Systus R [37], using a fully implicit backward

Euler scheme [START_REF] Belytschko | Non linear Finite Elements for Continua and Structures[END_REF] to perform the stress update.

A stress update algorithm consists of determining at a time t + ∆t, the stress tensor, from the knowledge of the various physical quantities at time t and the strain increment ∆ε ij . In a fully implicit scheme, the yield condition is enforced at the end of the step, contrary to fully explicit schemes. Hence, a fully implicit scheme is more accurate than an explicit scheme but the implementation is far more complex.

In the constitutive model, the additive decomposition of the strain rate tensor into elastic and plastic parts is assumed. The elastic response is thus hypoelastic [START_REF] Holzapfel | Non linear solid mechanics -a continuum approach for engineering[END_REF]. In incremental form the constitutive equation is taken to be :

σij = C ijkl ( εkl -εp kl ) (9)
where σij is a suitable objective stress rate and C ijkl is the fourth rank elastic tensor. Over one step, Eq. ( 9) can be approximated by:

∆σ ij = C ijkl (∆ε kl -∆ε p kl ) (10)
Eq. ( 3) can also be approximated by :

∆ε p ij = ∆λ ∂f ∂σ ij (11)
Substituting Eq. ( 11) into Eq. ( 10), and using the consistency condition and Eq. ( 8) leads to a set of non-linear equations where the unknowns are ∆λ and ∆ε p m . This set of equations is solved using a Newton scheme [START_REF] Belytschko | Non linear Finite Elements for Continua and Structures[END_REF].

Finite element model

Calculations have been performed using axisymmetric elements to model cone indentation and using a large displacement / large strain option (updated Lagrangian formulation, logarithmic strain). The objective rate of the stress tensor used is the Green-Naghdi rate [START_REF] Belytschko | Non linear Finite Elements for Continua and Structures[END_REF]. The mesh is specially refined near the contact zone. Its spatial extent is large enough to make the results unsensitive to the location of the outer boundaries [START_REF] Kermouche | [END_REF]. The size of the global mesh is one hundred times smaller than the contact radius. For a good representation of the contact geometry, the width of the elements has been determined so that at least 40 nodes are in contact for the deepest penetration. The whole mesh contains about 9000 elements and 9400 nodes (see Fig. 3). The contact between the indenter and the workpiece is assumed to be frictionless and loading is achieved by imposing a quasi-static displacement of the indenter which is pushed vertically into the workpiece. The indenter is assumed to be perfectly rigid.

Identification and numerical results

The constitutive model developed in this paper depends on two main physical quantities p c and q c which represent respectively the plastic limit in a pure hydrostatic state and the shear limit in a pure deviatoric state. The first quantity p c depends on the densification and has been identified from DAC experiments whereas the second (q c ) is constant and unknown. Before being able to confront the results of the finite element simulations with the experimental densification map obtained by Raman micro-spectroscopy, the value of q c has to be identified. Following Xin et al. [START_REF] Xin | Densification of fused silica : Effects on nanoindentation[END_REF], we used an inverse analysis based on the comparison of experimental load/displacement indentation curves with finite element simulations. For that purpose, five nanoindentation tests on fused silica samples were performed with the help of a Nano-indenter XP (MTS). The maximum penetration depth (about 2 µm) has been chosen large enough so that the influence of indentation length scales [40] on the measured load-penetration curve could be neglected. It has also been verified a posteriori by checking the hardness-penetration curve.

Notice that a Berkovich indenter has been used in this study whereas a Vickers indenter was used by Perriot et al. [7]. As both Vickers and Berkovich indenters are usually, as here, replaced by an equivalent cone of semi-apical angle of 70 degrees for load-penetration calculations, this substitution is not prejudicial.

One important result of Perriot et al. [7] is the densification gradient in the bulk of the material induced by indentation testing. The densification gradient induced by the equivalent cone indentation will be approximately the same as obtained with a 3D analysis, except just under the edge of the Vickers indenter. It seems also not wise to try fitting perfectly the densification map of Perriot et al according to the high measurement uncertainities. Therefore an axi-symmetric approximation is sufficient to test this constitutive model.

Inverse analysis

When elastoplastic solids are indented by rigid self similar indenters, the principle of geometric similarity can be applied [17,[START_REF] Kermouche | [END_REF]41]. This principle [START_REF] Tabor | The hardness of solids[END_REF] states that if two indentations are made by the same geometric shapes, then, whatever their size, the strain ε ij and stress σ ij distributions around the indentation will be geometrically similar. Hence, neglecting indenter tip defect and size effects [40], the loading curve is given by

F = Ch 2 (12)
where F is the load applied on the indenter, h is the penetration depth and C is a function of the mechanical properties [START_REF] Loubet | Microindentation techniques in Material Science and Engineering[END_REF]. One consequence of this principle is that there is one unique value of q c for which the numerical loading curve matches the experimental loading curve, the other mechanical properties being already known. Nevertheless, C is often difficult to quantify accurately by nanoindentation due to tip deffects and size effects. In this paper we pro-pose to base the inverse method on the evaluation of the plastically dissipated energy W pl and the elastically restored energy W el (see Figure 4) during one indentation cycle [17]. The maximal load used for the Raman mapping [7] being much higher than the load used for the nanoindentation testings, the ratios W pl /W tot and W el /W tot where W tot = W el + W pl have been chosen to perform the inverse analysis. These ratios represent the percentage of plastically dissipated energy and elastically restored energy over one indentation cycle. They are theoretically independent on the maximal load applied on the indenter and are thus well suited for this study.

Therefore, finite element simulations have been performed with different values of q c . The evolution of W el /W tot and W pl /W tot with q c are displayed in Figure 5. W el /W tot increases (respectively W pl /W tot decreases) linearly with q c . This result is not surprising since the higher q c , the lower the plastic strain. The experimental values of W pl /W tot and W el /W tot are respectively 0.33 ±0.01 and 0.67 ±0.01. The best fit have been obtained for q c = 6.5 GPa. The values of the different parameters for this model are summarized in table 1.

The experimental and numerical load-penetration curves are plotted in Fig. 6.

Both loading and unloading curves give a satisfactory agreement with the experimental nanoindentation curve.

Densification map

The densification map obtained with this constitutive model is plotted in Fig. 7. As observed with Raman microspectroscopy [7], the iso-densification lines are concentric. Densification is maximum just under the tip and decreases when the distance to the tip increases. Let us note that the maximum densification resulting from the finite element analysis is higher than the maximum Table 1 Mechanical properties of the constitutive model developed in this paper (equations ( 5) and ( 8))

5 Conclusion 3 FE mesh used for the modelling of the cone indentation test 23 4

 5323 densification measured by Perriot et al. It could be explained by the saturation of the densification process which is badly taken into account in the present constitutive model. In spite of that, the densification gradient resulting from the finite element analysis is in very good agreement with the experimental data. It can be considered that the constitutive model developed in this paper allows to reproduce succesfully the indentation-induced densification process of silica. We have developed a new constitutive law to model the plastic deformation of amorphous silica. The model takes into account the densification-induced hardening observed under purely hydrostatic loading with DAC experiments.The set of material properties has been identified with an inverse analysis based on the load-penetration curve on a fused silica sample. The resulting indentation-induced densification map obtained with this model is in good agreement with the experimental map obtained by Raman microscopy measurements. Therefore, this constitutive law for the plastic behaviour of amorphous silica accurately accounts for the existing indentation data. Still, one noticeable drawback of this model is that shear hardening is not considered. This shortcoming points out that the densification map and the load-penetration curve are not sufficient to identify the full plastic behavior of amorphous silica. For a more complete constitutive law valid for a wider class of loadings, it is necessary to define new experimental tests capable of providing additional data on the shear yielding of such materials. Hence, the shear behavior of silica should be more accurately identified and the role of shear on the densification process could be clarified. For this purpose, we are currently working on the finite element analyses of 3D Berkovich indentation, 3D Vickers indentation and Cube Corner indentation. These results will then be compared to experimental results in order to check the relevance of this constitutive model on other indentation data. Beyond the case of indentation, we are also working on the question of densification under scratch loading[43]. The deformation level in this test is much higher than in indentation experiments. It should thus be possible to observe new phenomena, and perhaps to get a better understanding of the scratch resistance of amorphous silica. 1 Mechanical properties of the constitutive model developed in this paper (equations (??) and (??)) 21 List of Figures 1 The proposed yield criterion for amorphous silica 21 2 Densification-induced hardening obtained from DAC experiments (extracted from Perriot [24]) and the linear law used in this paper 22 Definition of the plastically dissipated energy W pl and the elastically restored energy W el from a typical load-penetration curve 23 5 FE analysis of the dependence of W el /W tot and W pl /W tot on shear yield limit q c . The experimental values of W pl /W tot and W el /W tot are respectively 0.33 ±0.01 and 0.67 ±0.01.

Fig. 1 .Fig. 2 . 22 Fig. 3 .Fig. 4 .Fig. 5 .

 1222345 Fig. 1. The proposed yield criterion for amorphous silica

Fig. 6 .Fig. 7 .

 67 Fig. 6. Nanoindentation experimental and numerical load-penetration curves

Acknowledgements

The authors would like to thank A. Perriot, J.M. Bergheau and J.L. Loubet for helpful discussions.


N Coquillé

S Stachowski-Haberkorn

S Morin

É Parlanti

D Ménard

J Rouxel

L Haugarreau

V Dupraz

M Eon

J Vedrenne

S Boutry

J Rosebery

J Ezzedine

J C Grégoire

H Budzinski

N Tapie

P Pardon

L Chevance-Demars

INFLUENCED BY THE NATURAL DISSOLVED ORGANIC MATTER (DOM)?

As primary producers, microalgae are the basis of aquatic food webs. Thus, they can be directly affected by herbicides [1]. In their natural environment, these organisms evolve in complex and changing conditions, where various interactions take place. For instance, dissolved organic matter (DOM) may interact with pesticides. Indeed, DOM was shown to interact with pollutants and to affect their transport [2], fate, bioavailability, biodegradation and toxicity on organisms [3].

As primary producers, microalgae are the basis of aquatic food webs. Thus, they can be directly affected by herbicides [1]. In their natural environment, these organisms evolve in complex and changing conditions, where various interactions take place. For instance, dissolved organic matter (DOM) may interact with pesticides. Indeed, DOM was shown to interact with pollutants and to affect their transport [2], fate, bioavailability, biodegradation and toxicity on organisms [3].

) influence the toxicity of three herbicides (irgarol (I), diuron (D) and S-metolachlor (S)), singly and in mixtures (M1 and M2), to two marine phytoplankton species: Chaetoceros calcitrans (Cc) and Tetraselmis suecica (Ts).

) and the same 9 conditions with natural DOM added to the cultures (twice the environmental concentration).

) influence the toxicity of three herbicides (irgarol (I), diuron (D) and S-metolachlor (S)), singly and in mixtures (M1 and M2), to two marine phytoplankton species: Chaetoceros calcitrans (Cc) and Tetraselmis suecica (Ts).

) and the same 9 conditions with natural DOM added to the cultures (twice the environmental concentration).

Material & Methods Material & Methods

Results & Discussion Results & Discussion

Our results demonstrated the high toxicity of irgarol and herbicide mixtures for two marine phytoplankton species cultivated in absence of DOM. When DOM was added in the culture media, this study showed, in controlled conditions, that natural DOM seems to interact with herbicides and/or microalgae, leading to increased or decreased herbicide effects, depending on the species. Finally, this study demonstrates the importance to consider DOM as a major factor possibly involved in toxicity modulation in the environment.

Our results demonstrated the high toxicity of irgarol and herbicide mixtures for two marine phytoplankton species cultivated in absence of DOM. When DOM was added in the culture media, this study showed, in controlled conditions, that natural DOM seems to interact with herbicides and/or microalgae, leading to increased or decreased herbicide effects, depending on the species. Finally, this study demonstrates the importance to consider DOM as a major factor possibly involved in toxicity modulation in the environment.

Conclusions Conclusions

Control S-metolachlor 0,05 µg.L -1 0,5 µg.L -1 0,05 µg.L -1 0,5 µg.L -1 0,05 µg.L -1 0,5 µg.L -1

The toxic effects were measured on the doubling time (T D , hours), the photosynthetic efficiency (Y eff ), the intracellular ROS relative content using H 2 DCFDA and the intracellular lipid relative content using BODIPY (both by flow cytometry).

The toxic effects were measured on the doubling time (T D , hours), the photosynthetic efficiency (Y eff ), the intracellular ROS relative content using H 2 DCFDA and the intracellular lipid relative content using BODIPY (both by flow cytometry).

Influence of DOM on the pesticide effects regarding microalgae physiology

After exposure to herbicides together with DOM, significant effects were observed with I0.05, I0.5, M1 and M2 in Cc and with I0.5, D0.05, D0.5, S5, M1 and M2 in Ts (Fig. 4).

Influence of DOM on the pesticide effects regarding microalgae physiology

After exposure to herbicides together with DOM, significant effects were observed with I0.05, I0.5, M1 and M2 in Cc and with I0.5, D0.05, D0.5, S5, M1 and M2 in Ts (Fig. 4).

The exposure to DOM together with I0.5 or M2 caused quite similar effects on doubling time and photosynthesis efficiency in Cc. For Ts exposed to I0.5 and M2, results exhibited a much stronger toxicity of M2 compared to I0.5. An effect of D0.05, D0.5, S5 and M1 was also obtained on H 2 DCFDA fluorescence, meaning an increase of reactive oxygen species relative content when DOM was added.

When these results are compared to the exposures without DOM, it appears that for Cc, the presence of DOM seems to decrease the herbicide toxicity, while it is increased for Ts. The decrease or increase of toxicity can be explained by the complexation between DOM and herbicides. Indeed, this complexation can lead to a lower bioavailability of herbicides (because linked to DOM). The difference in herbicide toxicity between the two species in presence of DOM remains to be explained. Analysis of herbicide concentrations, dissolved organic carbon concentration and DOM (in progress) will provide valuable information to improve the understanding of interactions between microalgae, herbicides and DOM in these experiments.

The exposure to DOM together with I0.5 or M2 caused quite similar effects on doubling time and photosynthesis efficiency in Cc. For Ts exposed to I0.5 and M2, results exhibited a much stronger toxicity of M2 compared to I0.5. An effect of D0.05, D0.5, S5 and M1 was also obtained on H 2 DCFDA fluorescence, meaning an increase of reactive oxygen species relative content when DOM was added.

When these results are compared to the exposures without DOM, it appears that for Cc, the presence of DOM seems to decrease the herbicide toxicity, while it is increased for Ts. The decrease or increase of toxicity can be explained by the complexation between DOM and herbicides. Indeed, this complexation can lead to a lower bioavailability of herbicides (because linked to DOM). The difference in herbicide toxicity between the two species in presence of DOM remains to be explained. Analysis of herbicide concentrations, dissolved organic carbon concentration and DOM (in progress) will provide valuable information to improve the understanding of interactions between microalgae, herbicides and DOM in these experiments. Results obtained after exposure of Cc to M2 and I0.5 exhibited similar effects: the impacts induced by the mixture 2 might be due to irgarol. In Ts, M2 induced a greater growth and photosynthesis inhibition than I0.5 and a greater enhanced H 2 DCFDA fluorescence. This higher toxicity would suggest a possible synergistic effect when herbicides are mixed at the highest concentrations. Further experiments are needed to properly investigate this assumption.

Results obtained after exposure of Cc to M2 and I0.5 exhibited similar effects: the impacts induced by the mixture 2 might be due to irgarol. In Ts, M2 induced a greater growth and photosynthesis inhibition than I0.5 and a greater enhanced H 2 DCFDA fluorescence. This higher toxicity would suggest a possible synergistic effect when herbicides are mixed at the highest concentrations. Further experiments are needed to properly investigate this assumption.

Herbicide effects on microalgae

Significant effects were observed for almost all parameters after exposure to I0.5, and M2 for both species (Fig. 3): an increase of the doubling time, corresponding to a strong growth inhibition, a decrease of the photosynthetic yield and a decrease of the relative lipid content. In addition, Ts exhibited an increase of relative ROS content.

Herbicide effects on microalgae

Significant effects were observed for almost all parameters after exposure to I0.5, and M2 for both species (Fig. 3): an increase of the doubling time, corresponding to a strong growth inhibition, a decrease of the photosynthetic yield and a decrease of the relative lipid content. In addition, Ts exhibited an increase of relative ROS content.

Figure 2 -FigureFigure 3 -

 23 Figure 2 -Treatments applied to each species Figure 2 -Treatments applied to each species

Figure 3 -

 3 Figure3-Significant effects (ANOVA, p<0.05) detected in both species after exposure to S-metolachlor, irgarol and diuron, singly and in mixtures. Figure3-Significant effects (ANOVA, p<0.05) detected in both species after exposure to S-metolachlor, irgarol and diuron, singly and in mixtures.

Figure 4 -

 4 Figure 4 -Significant effects (ANOVA, p<0.05) detected in both species after exposure to S-metolachlor, irgarol and diuron, singly and in mixtures with and without DOM.Figure 4 -Significant effects (ANOVA, p<0.05) detected in both species after exposure to S-metolachlor, irgarol and diuron, singly and in mixtures with and without DOM.

Figure 4 -

 4 Figure 4 -Significant effects (ANOVA, p<0.05) detected in both species after exposure to S-metolachlor, irgarol and diuron, singly and in mixtures with and without DOM.Figure 4 -Significant effects (ANOVA, p<0.05) detected in both species after exposure to S-metolachlor, irgarol and diuron, singly and in mixtures with and without DOM.

Acknowledgement -This study has been carried out with financial support from the French National Research Agency (ANR) in the frame of the Investments for the future Programme, within the Cluster of Excellence COTE (ANR-10-LABX-45).