P Fraundorf

Melanie Lipp

Taylor Savage

Analogs for Unlayered-Graphene Droplet-Formation in Stellar Atmospheres

The unlayered spherical graphene cores, found in a subset of micron-sized presolar graphite onions extracted from primitive meteorites, show graphene but not graphite inter-sheet spacings, sheet "coherence widths" of about 4nm, and (in electron-phase-contrast images) evidence of edge-on obtuse-angle graphene-sheet junctions. The literature on solidification of metallic liquids predicts supercooling by as much as 30% below the melt-temperature, suggesting that at low pressures carbon vapor will predictably condense first as a liquid. Density measurements, possible TEM evidence of faceted-pentacones, and preliminary studies of 5 and 6 atom loop formation in a solidifying melt, all suggest that these unlayered graphene cores (which bear isotopic signatures of formation in asymptotic-giant-branch star atmospheres after third dredge-up) may have resulted from dendritic solidification of "slow-cooled" carbon droplets. Preliminary notes on a strategy for slow-cooling carbon vapor-condensation in the laboratory (with help from an evaporating oven) are also provided. Micron-sized presolar carbon-spheres condensed before the origin of our solar system, in the atmosphere of asymptotic-giant-branch stars from carbon atoms dredged up following nucleosynthesis in the stellar interior, have been extracted by dissolution from primitive carbonaceous chondrite meteorites, identified by noblegas and major-element isotopic analysis, and available to electron microscopists on earth for over two decades [START_REF] Bernatowicz | Constraints on stellar grain formation from presolar graphite in the murchison meteorite[END_REF][START_REF] Fraundorf | The core structure of pre-solar graphite onions[END_REF][START_REF] Croat | Presolar graphite from AGB stars: Microstructure and s-process enrichment[END_REF][START_REF] Mandell | Electron beam characterization of carbon nanostructures[END_REF] . Some of these spheres appear to consist of graphite condensed on carbide grains, but others instead show graphite-onion rims that surround a spherical-core of unlayered-graphene i.e. carbon that in diffraction shows graphite (hk0) spacings only. The high spatial-frequency tail of these graphene lines 5 suggests a graphene-sheet coherence-width in the 4-nm size range 1 , consistent with electron-phase-contrast images of single (sometimes segmented) edge-on sheets [START_REF] Fraundorf | The core structure of pre-solar graphite onions[END_REF][START_REF] Mandell | Electron beam characterization of carbon nanostructures[END_REF] .

II. OPEN QUESTIONS

Three key questions are: What process in a lowpressure atmosphere can form condensed spheres of graphene with densities approaching that of graphite, but with graphitic-layering absent? How might we synthesize such material in the laboratory? What interesting properties might it have as a result?

Our group has been focusing on three aspects of this problem: (i) Measurements of core-material density (relative to the graphite rim) in ultra-microtomed spheres using zero-loss/deflection analysis 6 via energy-filtered transmission electron microscopy or TEM, (ii) modeling of graphene-sheet nucleation and growth from a carbon melt to get atom-positions for TEM image simulation and estimate the likelihood of pentagonal-loop nucleation as explanation of apparent edge-on facets 2 using the Tersoff potential [START_REF] Tersoff | Empirical interatomic potential for carbon, with applications to amorphous carbon[END_REF] , and (iii) the laboratory synthesis of unlayered-graphene via slow-cooling of carbon vapor.

III. DENSITY

Carbon plasmon-peak energies 1 suggest an unlayered graphene-density near the 2.3[g/cc] of graphitic carbon. Zero-loss/deflection analysis [6] shows promise for inelastic mean-free-path mapping in a wide range of disordered materials, but specular reflection by diffraction can complicate life (e.g. via Pendullosung or thickness-fringe contrast) by putting electrons into Bragg storage if the crystals are thick enough for multiple scattering within the same crystal. Although our first analyses suggest a density closer to 1.8[g/cc] for presolar unlayered-graphene, this may not be confirmed without further calibration studies on the presolar specimens themselves.

IV. DROPLETS AND SOLIDIFICATION

Condensation of particles in a low-pressure stellaratmosphere will require local density increases 1 e.g. asso-FIG. 1. On-line strong-phase-object simulation of 15 facetednanocones in an 8-nm sphere with 26, 846 "non-touching" C atoms. The projected-potential power-spectrum at right shows graphene {100} and {110} rings, as well as rel-rods associated with some of the edge-on graphene-sheets. ciated with jets of material around sunspots, and might take place one C atom at a time, by clustering of polycyclic aromatic hydrocarbons, or by condensation of droplets from a carbon vapor. The latter has the advantage that it will give us the spherical core-shape and density. Although thermodynamics predicts that below pressures of 10.8[MPa], carbon on warming sublimates to a vapor at around 3915K, containerless metals on cooling generally supercool in the liquid state until as much as 30% below the liquid condensation temperature 8 , de-pending on the amount of rearrangement needed between liquid and solid state configurations.

A liquid condensation temperature for carbon around 4600K, in this context, along with the major rearrangement between 12-neighbor coordination in an fcc liquid and the graphene-sheet structure, thus suggest that supercooling to temperatures as low as 3220K (2950C) might be expected in the absence of condensation seeds (like the carbide grains seen at the center of other particles). This synthesis path might also explain the possible TEM evidence 2,4 that the graphene sheets are flat, except where they are seen edge-on to connect at sharp bends whose geometry is suggestive of nucleation on pentagonal (rather than hexagonal) carbon-loops. So far, we are assembling graphene-sheet and graphene-pentacone structures for TEM image simulation (Fig. 1), but only beginning to look at simulated annealing models of pentagonal versus hexagonal loop nucleation.

V. SYNTHESIS IN THE LABORATORY

Finally, slow-cooled carbon-vapor synthesis in the inside of a resistively-evaporating carbon-cylinder now looks possible. However we are only beginning to examine carbon spheres created as a result.

Acknowledgments 2

References 2

ACKNOWLEDGMENTS

Thanks to Roy Lewis at U. Chicago and T. J. Bernatowicz at Washington U. for the presolar specimens used here.


José P Baltanás

José L Trueba

Miguel A F Sanjuán

Energy dissipation in a nonlinearly damped Duffing oscillator

Keywords: Energy dissipation, Duffing oscillator, Nonlinear friction

In this paper, we study the effect of including a nonlinear damping term proportional to the power of the velocity in the dynamics of a double-well Duffing oscillator. In particular, we focus our attention in understanding how the energy dissipates over a cycle and along the time, by the use of different tools of analysis. Analytical and numerical results for different damping terms are shown, and the presence of a discontinuity and an inversion of behavior depending on the initial energy are discussed. An averaged power loss in a period is defined, showing similar characteristics as the energy dissipation over a cycle, although no discontinuity is present. The discontinuity gap which appears for the energy dissipation at a certain value of the initial energy decreases as the power of the damping term increases and an associated scaling law is found.

Introduction

The phenomenon of dissipation is ubiquitous in physical and mechanical systems, always meaning a loss of energy due to friction forces. Friction has been a topic of technological attention since the early days of our civilization, and its interest in scientific and engineering research remains today [START_REF] Feeny | Friction as a nonlinearity in dynamics: a historical review[END_REF]. Most systems are subject to a different combination of terms of damping, each one of them with different effects. These dissipative forces may act on a system in the form of dry friction, such as Coulomb friction, which is typically a nonlinear damping mechanism, and which has deserved much attention recently in relation to its physical origins [START_REF] Persson | Sliding Friction: Physical Principles and Applications[END_REF][START_REF] Röder | Dry friction: modelling and energy flow[END_REF]. Other typical sources of friction is drag or aerodynamics friction, which is typically proportional to the square of the velocity, and viscous forces which are linear in the velocity. The main advantage of the viscous linear damping term comes from its simplicity, since the linearity of the governing differential equation, in the simple case of a harmonic oscillator, is not altered.

In ordinary applications of nonlinear dynamical systems such as the motion of an object through a fluid or gas, frictional and drag forces are usually very complex, and power-law dissipative terms in the velocity are used. However, one of the problems of considering these and other forms of damping, besides the common linear ones, is that they give rise to nonlinear terms in the governing differential equations even for the simplest system, without the possibility of obtaining exact solutions. The mathematical form of the friction forces is usually experimentally determined with the use of wind tunnels or water tanks. Among the possible models, one of the simplest empirical mathematical model of the dissipative force is taken as f (v) ∝ -v|v| p-1 , where v represents the velocity, with its positive or negative sign, and p ≥ 0. Notice that the absolute value of the velocity is needed to safeguard the fact that the dissipative force must be contrary to the motion of the system. Taking this fact into account, one can recover from this expression the important cases of Coulomb damping (p = 0), linear viscous damping (p = 1) and draglike friction (p = 2). A historical review of the applications of the Coulomb damping (p = 0) can be found in [START_REF] Feeny | Friction as a nonlinearity in dynamics: a historical review[END_REF], whereas the extensively used case p = 1 and its applications are considered, for example in [START_REF] Kelly | Fundamentals of Mechanical Vibrations[END_REF][START_REF] Virgin | Introduction to Experimental Nonlinear Dynamics[END_REF]. Phenomenological models accounting for nonlinear dissipation forces with higher values of p have been considered in some applied sciences, such as ship dynamics, vibration engineering, and motion of projectiles, among others [START_REF] Bikdash | Melnikov analysis for a ship with a general roll-damping model[END_REF][START_REF] Falzarano | Application of global methods for analyzing dynamical systems to ship rolling motion and capsizing[END_REF][START_REF] Chan | Estimation of nonlinear damping coefficients from large-amplitude ship rolling motions[END_REF][START_REF] Mestre | The Mathematics of Projectiles in Sport[END_REF][START_REF] Huilgol | The motion of a mass hanging from an overhead crane[END_REF][START_REF] Ravindra | Stability analysis of a non-linearly damped Duffing oscillator[END_REF]. A physical interpretation of these higher order damping terms may be given by assuming a power series expansion on the velocity of more complicated damping forces.

The effect of some physical parameters, such as wave amplitude and nonlinear damping, on the transient motions and the global system behavior of a capsizing ship was investigated in [START_REF] Falzarano | Application of global methods for analyzing dynamical systems to ship rolling motion and capsizing[END_REF]. In particular, they use a model of a double-well and a single-well Duffing oscillator with a linear plus a quadratic damping term. The stability of a nonlinearly damped hard Duffing oscillator was analyzed in [START_REF] Ravindra | Stability analysis of a non-linearly damped Duffing oscillator[END_REF], and the role of nonlinear dissipation in soft Duffing oscillators in [START_REF] Ravindra | Role of nonlinear dissipation in soft Duffing oscillators[END_REF]. Moreover, in [START_REF] Ravindra | Chaotic response of a harmonically excited mass on an isolator with non-linear stiffness and damping characteristics[END_REF] it was shown that nonlinear damping could be used as a mechanism to suppress chaos. The effect of nonlinear dissipation on some properties of the dynamics of the universal escape oscillator, such as the threshold of period-doubling bifurcations, fractal basins boundaries and the destruction of the basins of attraction was considered in [START_REF] Sanjuán | The effect of nonlinear damping on the universal escape oscillator[END_REF], and some analytical estimates of the nonlinear damping effects for the Duffing oscillator and for the pendulum, using Melnikov techniques, in [START_REF] Trueba | Analytical estimates of the effect of nonlinear damping in some nonlinear oscillators[END_REF]. A rather interesting example of application of this nonlinear damping appears in [START_REF] Pfenniger | Dissipation in barred galaxies: the growth of bulges and central mass concentrations[END_REF] in the context of dissipation in barred galaxies and the analysis of the mechanisms by which mass concentrations are produced, where a model of a quadratically damped pendulum is used. Dynamical friction of this sort is common in galactic dynamics [START_REF] Binney | Galactic Dynamics[END_REF].

A discussion of a general analysis of damping appears in [START_REF] Bikdash | Melnikov analysis for a ship with a general roll-damping model[END_REF] in the context of the naval architecture literature, including combinations of linear plus quadratic damping terms, linear plus cubic damping terms or even combinations of linear plus quadratic plus cubic terms. A general damping function of the form

D( ẋ) = N p=1 b p ẋ| ẋ| p-1 (1)
is used which includes a general combination of different damping terms (b p ≥ 0) that are powers in the velocity. In our analysis on the energy dissipation, we consider terms of the form b ẋ| ẋ| p-1 , and our results apply when a combination of terms of the above form is used.

From the nonlinear dynamics point of view, a dissipative system with a linear dissipative term has a rate of contraction of phase space which is constant. This is not the case when the dissipative force is nonlinear and is proportional to a power of the velocity, since it implies that the rate of contraction of phase space is velocity dependent. The idea of the volume contraction in phase space is somehow related to the logarithmic decrement of linear oscillators with linear damping and it might be also a generalization of the same idea for general nonlinear oscillators. It is common practice in physics and engineering to analyze the effects of damping by measuring the decay of free oscillations and the loss of energy in physical systems. A discussion of damping in nonlinear oscillators in the context of nonlinear dynamics appears in [START_REF] Virgin | Introduction to Experimental Nonlinear Dynamics[END_REF].

In the present work, we analyze the effect of using nonlinear damping terms on the loss of energy over a cycle and as a function of time, and the decay of the free oscillations for the nonlinear double-well Duffing oscillator. Commonly, this nonlinear oscillator has been analyzed with the inclusion of a simple viscous (linear) damping term, however, we have here a model of a nonlinear oscillator with a nonlinear damping term. This system plays a role of a paradigm in nonlinear dynamics and for this reason we aim at carrying out analytical and numerical computations which provide information on how the energy dissipates for a better understanding of its physical mechanisms.

The organization of the paper is as follows. A general overview of how the energy dissipates, including the analysis of the decay of the free oscillations, and a discussion on the results (analytical and numerical) concerning the loss of the energy over a cycle, are described in some detail in Section 2. In Section 3, we study how the energy dissipates as a function of time and analyze the logarithmic decrement associated to the decay of the amplitude of the oscillations. Since the energy dissipated over a cycle depends strongly on the period of the motion, we define the averaged power loss and devote Section 4 to its study and to analyze some scaling laws which appear in the analysis of the energy dissipated over a cycle for different values of the damping exponents. In Section 5 we give the conclusions. And finally, in Appendixes A and B, the analytical solutions for the undamped and unforced Duffing oscillator and a list of integrals involving Jacobian elliptic functions used throughout this work, are included.

Overview of the energy dissipation

Description of the model and decay of the oscillations

As we mentioned earlier, our goal is to analyze the effect of nonlinear damping terms on the behavior of a nonlinear oscillator. In particular, we focus our attention on the nonlinear Duffing oscillator. We consider a general nonlinear dissipative force proportional to the power of the velocity of the form

F p = -b ẋ| ẋ| p-1 , ( 2
)
where b (b 1) is the damping coefficient, p (p real) the damping exponent and ẋ stands for the velocity. Our purpose here is to estimate the energy dissipation of the Duffing oscillator when a dissipative term of the previous form is included. Thus, the general model equation for the nonlinearly damped Duffing oscillator reads

ẍ + b ẋ| ẋ| p-1 -x + x 3 = F cos ωt, ( 3
)
where F is the amplitude and ω the frequency of the external perturbation. For the objectives of this paper, our interest is focused in the decay of the free oscillations and how the energy dissipates over a cycle and as a function of time, so we are interested in the unforced oscillator, i.e., we fix F = 0 throughout the whole paper.

One of the first consequences of the damping on a free oscillator is the decay of the amplitude of the oscillations. In order to observe and analyze this phenomenon in our case, we have numerically computed the solutions of the free nonlinearly damped Duffing oscillator for different values of the initial energy. We have numerically integrated Eq. ( 3) using a fourth-order Runge-Kutta integrator. The temporal evolutions of the amplitude x(t) are plotted in Fig. 1, when b = 0.01, for different values of the damping exponent p and the energy E in order to visualize how the energy dissipates, and how it affects the amplitude of the oscillation. This figure shows six panels in three rows. The first row refers to the case of an orbit situated in the right well, and Fig. 1(a) shows the case of dry friction, where the energy dissipates rather quickly with an apparent linear decrement settling down at the bottom of the well. We can see in Fig. 1(b) for the same energy but with a cubic friction term that the energy takes longer time to decay. Note that one thing is how the energy dissipates as a function of time, and a different one is the dissipation over a cycle. Fig. 1(c) and(d) shows, for an orbit situated outside the well, that the amplitude decays faster for the case of dry friction than for the case of cubic friction. The last row, Fig. 1(e) and (f), corresponds to the case of an orbit of higher initial energy. For p = 0, the energy takes longer to decay than for p = 3, as in the previous cases. However, in this case, for p = 0 the energy dissipated over a cycle for p = 0 is lower than for p = 3, as we will see in the next section.

Energy dissipation over a cycle

Now we are interested in evaluating how the energy dissipates over a cycle. Following the argumentation given in [START_REF] Kelly | Fundamentals of Mechanical Vibrations[END_REF], the energy dissipated E over one cycle of period T due to the damping force F p is

E = T 0 F p ẋ dt. (4)
For simplicity during all our calculations, we define the quantity

Δ p ≡ E -b = T 0 ẋ2 | ẋ| p-1 dt, (5)
where this last integral can be estimated analytically in some cases (when p is an integer), whereas otherwise numerical techniques are used. For integer values of the damping exponent p, all we need to evaluate this last integral is the amplitude x(t) for the unforced and undamped oscillator and its period T . These expressions are given in Appendix A. All the necessary integrals involving Jacobian elliptic functions needed to compute these expressions have been included in Appendix B. We have carried out analytically the computations for the energy dissipated over one cycle Δ p for the damping exponents p = 0, 1, 2, 3, which correspond to the cases of Coulomb damping, linear viscous damping, quadratic and cubic damping, respectively. All the results are displayed in Table 1. In all these expressions, the parameter α is defined in terms of the energy E of the orbit as α = √ 1 + 4E, and K(m) and E(m) are the complete elliptic integrals of the first and second kind, respectively. Information concerning Jacobian elliptic functions is found in [START_REF] Lawden | Elliptic Functions and Applications[END_REF][START_REF] Milne-Thomson | Handbook of Mathematical Functions[END_REF]. The values of Δ p are computed for different orbits with initial energies below the local maximum and above the local maximum, which are referred in the table as regions I and II, respectively.

In Fig. 2, we have plotted the analytical results for Δ p for integer values of p increasing from p = 0 to 3, (solid line), and we have also included the corresponding numerical results (symbols). As is observed, the agreement between the exact solutions and the numerical computations is excellent.

+ b ẋ| ẋ| p-1 -x + x 3 = 0 a p Region Δ p 0 I Δ 0 = 2( √ 1 + α - √ 1 -α) II Δ 0 = 4 √ 1 + α 1 I Δ 1 = 4 3 1 + α 2 (α -1)K 2α 1 + α + E 2α 1 + α II Δ 1 = 4 3 √ α (α -1)K 1 + α 2α + 2E 1 + α 2α 2 I Δ 2 = 4 15 [(1 + α) 3/2 (3α -2) + (1 -α) 3/2 (3α + 2)] II Δ 2 = 8 15 (1 + α) 3/2 (3α -2) 3 I Δ 3 = 4 35 1 + α 2 (5α 2 -4)(α -1)K 2α 1 + α + 4(2α 2 -1)E 2α 1 + α II Δ 3 = 4 35 1 + α 2 (5α 2 -3α -4)(α -1)K 1 + α 2α + 8(2α 2 -1)E 1 + α 2α
a The region I indicates the region inside the double well where the initial energy of an orbit is -1 4 < E < 0, while the region II refers to an orbit with initial energy E > 0, i.e., outside the double well. The parameter α is related to the initial energy E of an orbit through the relationship α = √ 1 + 4E. K(m) and E(m) are the complete elliptic integrals of first and second kind, respectively, where 0 < m < 1 is the elliptic modulus. A plot of the energy dissipation over a cycle versus the initial energy for the different values of p considered here is represented in Fig. 2. After the observation of this figure, we can conclude that the energy dissipated over one cycle, Δ p , increases as α increases for a fixed p. This can be explained if we realize that, given the dependence of the damping term with the velocity of the system (see Eq. ( 2)), it is obvious that the faster the system moves, the stronger the damping effect and therefore the energy dissipation is. This indicates that the energy dissipation over a period for a fixed p should increase with α, i.e., with the total energy, since the velocity of the system also increases when the total energy increases.

On the other hand, it is worth noticing how Δ p behaves for a given α, when the damping exponent p varies. We observe that Δ p decreases when p increases for a certain range of values of α, but it inverts its behavior when α is above a value around 1.7. The explanation of this fact can be analyzed again in terms of the damping term. We may consider two values of p, say p 1 and p 2 , such that p 1 > p 2 . It can be shown that the velocity of the orbit ẋ is such that | ẋ| < 1 when α < √ 3 and that only when α ≥ √ 3 the velocity is greater than 1 for some instants of time. Taking into account that the nonlinear dissipative forcing term F p has a dependence on p given by Eq. ( 2), one would expect the damping effect to be weaker for p 2 than for p 1 when α is below a certain value close to √ 3 (| ẋ| < 1), whereas the opposite occurs when α is above this value (| ẋ| > 1). This implies that, roughly, Δ p is a decreasing function of p for α < √ 3 and an increasing function of p when α ≥ √ 3, which clearly corresponds to our results. Finally, the discontinuity at α = 1 is understood if we realize that for values of α < 1, the system is forced to move inside one of the two wells, whereas for values of α > 1 the system has enough energy to move outside the wells (see Fig. 3). So, it is straightforward to see that a discontinuity in both the period (see Fig. 4), and the velocity of the system, and therefore, in Δ p , must appear. It is interesting to note as well the gap of the curves at the discontinuity (α = 1), which is smaller as the damping exponent p increases. A scaling law associated to this phenomenon will be discussed in Section 4.

For completeness, we have carried out numerical computations involving noninteger values of p. Some curves have been depicted in Fig. 5(a) and (b). In these figures, we have displayed Δ p versus α for p ranging from 0 to 1 at intervals of 0.25 (Fig. 5(a)), and from p = 1 to 2 with the same interval (Fig. 5(b)). The dependence of Δ p on α for noninteger values of p is very similar to that exhibited by the system in Fig. 2. This means that the energy loss increases when α increases for a given p, whereas Δ p decreases when p increases for a given value of α, except from a certain value of α near √ 3. It is worthwhile to mention here the importance of the system on the initial conditions, since it has multiple equilibria. Numerical computations for higher values of the damping exponent p show similar behavior as the one described here.

Energy dissipation as a function of time

The results we have described previously are not easily obtained from the temporal evolution of the amplitude of the oscillation x(t) in Fig. 1. For example, at first glance, one would say that the energy dissipation for p = 0 and α = 0.5 (Fig. 1(a)) is stronger than for p = 0 and α = 5 (Fig. 1(e)), whereas our results show that it is just the opposite. To show that our computations make sense, we have displayed in Fig. 6 some numerical estimates of the total energy of the system E(t) for several values of p and α. It is clear that the behavior of E(t) as a function of p and α follows that encountered for Δ p , despite the fact that the solutions for the unforced and undamped Duffing oscillator used in the computation of Δ p are quite different from the true trajectories shown in Fig. 1. In any case, it should be noticed that Δ p represents in fact the loss of energy in the first cycle of the orbit, and that this dissipation of energy should be different in successive cycles. For example, in Fig. 6(c), which corresponds to α = 2, the strongest energy loss corresponds to p = 0 for long values of the time variable t, whereas Δ p is bigger for p = 3 than for p = 0, which coincides with the behavior of E(t) for small values of t.

Given a dynamical system, the divergence theorem relates the volume contraction in phase space with the divergence of the flow f in the following manner [START_REF] Virgin | Introduction to Experimental Nonlinear Dynamics[END_REF] V (t)

V (t) = div f. ( 6
)
In the linear damping case div f = -b and hence V (t) = V 0 e -bt , which shows an exponential contraction of the volume in phase space. In a general case, for our model of nonlinear damping, we have

V (t) V (t) = -bp| ẋ| p-1 , (7)
so the contraction is not exponential (note that this expression does not apply for the case p = 0, since the flow is discontinuous at ẋ = 0). The previous discussion is somehow related to the idea of the logarithmic decrement, which is a useful insight on the manner in which energy dissipates, defined as the natural logarithm of the ratio of the amplitudes of the oscillation on successive cycles [START_REF] Kelly | Fundamentals of Mechanical Vibrations[END_REF] by

δ n = ln x(t) + (n -1)T x(t + nT) . ( 8
)
In this expression, x(t) is the amplitude of the oscillation in the time t, T the period and the index n indicates the successive cycles. In the general case, δ n depends on the damping exponent p and the initial energy α. However, it is well known that in the case of a harmonic oscillator with viscous damping, where the above quantity usually appears, δ n depends only on the damping coefficient (b in Eq. ( 2)), but not on n or on the initial energy.

As a method to characterize the way in which the nonlinear oscillator dissipates energy we have numerically explored the behavior of δ n , since it gives us information about the decay of the oscillations for any general damping term. In Fig. 7, we have depicted δ n versus n for several values of p and α. This figure shows that the logarithmic decrement depends on n (i.e., it depends on the pair of successive cycles considered) and also on the initial energy, α, and the damping exponent p. We can analyze this dependence from different points of view. A result of our analysis is that, when p is fixed, δ 1 depends on α just in the way that one would expect from the behavior of Δ p , and the same can be said when α is fixed and p changes. Moreover, by the observation of the figures it can be deduced that this analysis can be extended to δ n when n is small enough (typically lower than 5), despite Δ p only concerns to the first cycle (n = 1). Another interesting result is that δ n is an increasing function of n for all α when p = 0, and a decreasing function of n for all α when p = 2, 3, but it changes from a decreasing to an increasing function of n if p = 1 when α is varied. In fact, this implies that Δ p passes through a maximum when considered as a function of the logarithmic decrement in the case of viscous damping, exactly as in the case of a harmonic oscillator [START_REF] Kelly | Fundamentals of Mechanical Vibrations[END_REF], whereas this is not true when p = 1.

Power dissipation and scaling law

The definition of the energy dissipation over a cycle that we have used (Eq. ( 5)) shows clearly that Δ p explicitly depends on the period T in such a way that the energy dissipation over a cycle increases as far as the period increases. This suggests that it would be convenient to make use of a quantity independent of the duration of the cycle. For this reason, we define the averaged power loss in a period, Δ p /T .

We have computed this magnitude for several values of the damping exponent p, and the results of our computations are plotted in Fig. 5(c) and (d). Analogous to Fig. 5 at α = 1 present in Fig. 2 for Δ p seems to disappear here, in other words, what we have termed averaged power loss looks like a continuous function of α. We have studied the behavior of Δ p /T for different values of the damping exponent p (analytically for integer values of p and numerically for noninteger ones), having found that Δ p /T is in fact a continuous function of α. It should be noticed as well the presence of local maxima in the curves for both cases α < 1 and α > 1. In the first case, the maximum location does not change when p varies, but it does in the case α > 1 (see Fig. 5(c), where two maxima are shown for α > 1).

As we mentioned in Section 2, the observation of Fig. 2 related to the energy dissipation Δ p versus the energy α implies that the gap that appears at the discontinuity α = 1 varies for the different damping exponents p. This suggest that a certain kind of scaling law should exist for λ(p), where λ(p) represents the gap in the discontinuity at α = 1 for the Δ p curves as a function of p. We have computed λ(p) numerically for different values of p varying from p = 0 to 15 (p being an integer, see Fig. 8). With the sequence we have obtained, we have applied the Aitken's method of quick convergence [START_REF] Press | Numerical Recipes in C[END_REF] in order to get evidence of such a scaling law. In order to apply this method to our case, we build a new sequence d p = λ(p -1)λ(p) from the sequence of λ(p), and afterwards define the quantity μ p = d p /d p+1 . Then we analyze how μ p behaves for increasing values of p. If μ p tends to a constant, then a scaling law exists for λ(p). We have found numerically that such a constant value exists and it is close to 1.48. This scaling law can also be calculated using a semi-analytical approach. If the points (p, λ(p)) are depicted in a semi-logarithmic scale, we see that the last ones (from p = 8 to 15) can be fitted to a straight line, i.e., they obey a decaying exponential law of the type λ(p) = λ 0 exp(-kp), where in our case, λ 0 = 1.24117 and k = 0.39102. Using this expression to calculate μ p , one obtains μ p = e k = 1.478 for all p, which is very close to the result obtained when the discrete sequence is used.

Conclusions

We have thoroughly analyzed the behavior of a nonlinearly damped Duffing oscillator with a dissipative term proportional to the power of the velocity. The energy dissipation over a cycle for different damping terms, such as dry, linear, quadratic and cubic friction, has been analytically calculated, providing explicit formulae for these cases. These formulae explicitly show the dependence of the energy dissipation on the initial energy α and on the damping exponent p. Numerical computations have also been performed for noninteger powers in the damping term, resulting in a dependence very similar to the case of integer exponent.

The presence of a discontinuity at the value of α = 1 and an inversion of the behavior of Δ p for values of the energy above α √ 3 have been analyzed and explained in terms of the dynamical behavior of the model. These results are contrasted with the evolution of the energy dissipation as a function of time for different values of the damping exponent p and initial energy α, finding that the dependence on these parameters for small values of time t is coherent with the results found for the dissipation of the energy Δ p . In order to avoid the dependence of the energy dissi-pation on the period, we have defined an averaged power loss in a period. The information provided by this new magnitude is rather similar to the one provided by the energy dissipation over a cycle, although no discontinuity is present in this case.

As a usual tool in the study of the decay of the oscillations, we have also computed the logarithmic decrement for our system for different values of p and α, having obtained consistence with our previous results and some resemblances with the harmonic oscillator when p = 1.

Finally, we have analyzed the discontinuity gap which appears for the energy dissipation at α = 1, and we have obtained numerically that this gap decreases as the power of the damping term increases. The analysis of the sequence of gaps for different damping exponents p has shown the presence of a scaling law associated with this dependence.

Fig. 1 .

 1 Fig. 1. Decay of the amplitude of the free oscillations of the nonlinearly damped Duffing oscillator ẍ + b ẋ| ẋ| p-1x + x 3 = 0 for different damping exponents p and different values of the initial energy (α = √ 1 + 4E): (a) p = 0 and α = 0.5, (b) p = 3 and α = 0.5, (c) p = 0 and α = 1.25, (d) p = 3 and α = 1.25, (e) p = 0 and α = 5, (f) p = 3 and α = 5.

Fig. 2 .

 2 Fig. 2. Dependence of the loss of the energy over a cycle Δ p versus the initial energy (α = √ 1 + 4E), for several integer values of the damping exponent p. The symbols refer to the numerical computations, whereas solid lines represent the analytical calculations. The discontinuity at α = 1 is a manifestation of the transition from the double well to the single well above the local maximum. Notice the inversion of the behavior around the value α = √ 3.

Fig. 3 .

 3 Fig. 3. Plot of the double well potential V (x) = 1 4 x 4 -1 2 x 2 , associated to the Duffing oscillator. Motions inside the double well corresponds to initial energies -1 4 < E < 0 and consequently 0 < α < 1. For values of the initial energy E > 0, α > 1, the motions are outside the double well.

Fig. 4 .

 4 Fig. 4. Period of the free Duffing oscillator as a function of α = √ 1 + 4E.Note that T (α) → ∞ both when α → 1 + and α → 1 -, but simultaneously when α → 1 + the value of T (α) is twice the value of T (α) when α → 1 -. This divergence is explained by the existence of a homoclinic orbit due to the presence of a hyperbolic fixed point in the phase space for α = 1.

Fig. 5 .

 5 Fig. 5. (a) and (b) Energy dissipation losses over a cycle for the nonlinearly damped Duffing oscillator as a function of α = √ 1 + 4E for different real values of the damping exponent p. (c) and (d) Averaged power losses over a cycle for the nonlinearly damped Duffing oscillator as a function of α for different real values of the damping exponent p. In the left panels (a) and (c), circles stand for p = 0, triangles for p = 0.25, crosses for p = 0.5, stars for p = 0.75 and squares for p = 1. On the other hand, in the right panels (b) and (d), circles correspond to p = 1, triangles to p = 1.25, crosses to p = 1.5, stars to p = 1.75 and squares to p = 2. Note that the discontinuity at α = 1 disappears for the averaged power loss Δ p /T .

Fig. 6 .

 6 Fig. 6. Total energy of the nonlinearly damped Duffing oscillator versus time for different values of p and α = √ 1 + 4E. Note that the behavior of E(t) follows that of Δ p when t is not too large, as is explained in the text: (a) α = 0.5, (b) α = 1.25, (c) α = 2 and (d) α = 5.

Fig. 7 .

 7 Fig. 7. Logarithmic decrement for successive cycles for different values of α = √ 1 + 4E and p. The dependence on this variable is discussed in the text. The symbols have the same meaning as in Fig. 2, i.e., circles stand for p = 0, triangles for p = 1, stars for p = 2 and inverted triangles for p = 3.

  Fig. 8. Value of the discontinuity of Δ p at α = 1 for several integer values of the damping exponents p. Symbols correspond to the numerical computations. The solid line joining the symbols is plotted as a simple visualization aid, and has no physical meaning.

Table 1

 1 This table shows the analytical results of the energy dissipation over a cycle Δ p for different integer values of the damping exponent p for the nonlinearly damped Duffing oscillator ẍ

Acknowledgements

This work has been supported by the Spanish DGES under project PB96-0123, and the Spanish Ministry of Science and Technology under project BFM2000-0967.

Appendix A. Analytical solutions of the Duffing oscillator

We want to explicitly obtain the exact solutions for the double-well Duffing oscillator ẍx + x 3 = 0.

(A.1)

Our objective is to look for analytical expressions for the solutions and the period of the motions of the system. The double-well potential V (x) = 1 4 x 4 -1 2 x 2 is shown in Fig. 3, and we have two different kinds of motion depending on the energy of the orbit. There is a motion inside a well, for values - 1 4 < E < 0, where the solution of the period is given by

where K(m) is the complete elliptic integral of first kind [START_REF] Lawden | Elliptic Functions and Applications[END_REF][START_REF] Milne-Thomson | Handbook of Mathematical Functions[END_REF], and m is the elliptic modulus. Inside one of the small well -1 4 < E < 0 and consequently 0 < α < 1, and the elliptic modulus m satisfies the condition 0 < m < 1, as expected. Fig. 4 shows the explicit dependence of the period of the Duffing oscillator on the energy.

The exact solution for an orbit in this region is given by

where α = √ 1 + 4E and -1 4 < E < 0. The period of an orbit for energies E > 0 is given by

where as E > 0, then α > 1, and consequently m satisfies the condition 1/ √ 2 < m < 1. The exact solution for an orbit in this region is given by

where α = √ 1 + 4E and E > 0.

Appendix B. List of integrals involving elliptic functions

We include a list of integrals that have been used in this work, based in [START_REF] Lawden | Elliptic Functions and Applications[END_REF][START_REF] Milne-Thomson | Handbook of Mathematical Functions[END_REF]. In the list, sn(x, m), cn(x, m) and dn(x, m) are the Jacobian elliptic functions, and K(m) and E(m) are the complete elliptic integrals of the first and second kind, respectively.


Sean Mcnamara

Jérôme Crassous

Axelle Amon

Eshelby inclusions in granular matter: theory and simulations

Keywords: numbers: 62.20.F-, 83.80.Fg, 62.20.D-

We present a numerical implementation of an active inclusion in a granular material submitted to a biaxial test. We discuss the dependence of the response to this perturbation on two parameters: the intra-granular friction coefficient on one hand, the degree of the loading on the other hand. We compare the numerical results to theoretical predictions taking into account the change of volume of the inclusion as well as the anisotropy of the elastic matrix.

I. INTRODUCTION

Despite decades of intensive research [START_REF]Critical State Soil Mechanics[END_REF][START_REF]Statics and Kinematics of Granular Materials[END_REF][START_REF] Davis | Plasticity and geomechanics[END_REF][START_REF]Bifurcation analysis in geomechanics[END_REF], the mechanics leading to strain localization in soils are not well identified. For geo-mechanicians, strain localization is the result of mechanical instabilities [START_REF]Statics and Kinematics of Granular Materials[END_REF][START_REF] Davis | Plasticity and geomechanics[END_REF][START_REF]Bifurcation analysis in geomechanics[END_REF]. When the stress exceeds a threshold (i.e. given by Mohr-Coulomb criterium for the stress component ratio), the material fails. Before such failure, soils flow plastically [START_REF]Critical State Soil Mechanics[END_REF][START_REF] Davis | Plasticity and geomechanics[END_REF][START_REF]Bifurcation analysis in geomechanics[END_REF], but the link between this plastic flow and failure is quite unclear [START_REF]Bifurcation analysis in geomechanics[END_REF]. Recent approaches to describe the rheology of granular materials propose non-local constitutive relations inspired by models coming from soft glassy materials [START_REF] Kamrin | [END_REF]6] and introducing a new variable to describe the local state of the material (the fluidity). The physical interpretation of this fluidity is still unclear and debated [6] Metallic glasses are another example of materials where plastic flow occurs before strain localization. It has been suggested that for these amorphous materials [7,8], the macroscopic plasticity is the result of numerous elementary plastic events. An elementary event may be viewed as a small and local reorganization at the molecular scale. Such events have been found in numerical simulations [9][10][11]. This description in term of elementary plastic events seems also relevant for soft glassy materials (such as dense emulsions or foams) or compressed granular materials [10]. Such elementary events have been indeed reported in some experimental studies [12][13][14][15]21].

A very important consequence of those elementary plastic events is that they are able to interact mechanically. Indeed, when they occur, they redistribute stress in the material [16]. For an elastic material, this redistribution is long-ranged and anisotropic, so that it may trigger other rearrangements further in the material along preferential directions. This coupling between plastic events is believed to be at the origin of formation of spatial heterogeneities [17][18][START_REF] Lin | Proc. Nat. Acad. Sci. Am[END_REF][START_REF] Lin | [END_REF] (such as shear band) at the macroscopic scale.

Consequently a good understanding of the stress redistribution around a single event is needed. It is a subject that has been already deeply studied, especially numerically [10,11,22], while experimental proofs of this redistribution are scarce [12,15,24]. Those studies show that the stress redistribution around a local event may be well described using a computation done by Eshelby [START_REF] Eshelby | Proc. R. Soc. Lond. A[END_REF] where an inclusion embedded in an elastic matrix spontaneously changes shape. The analytical solutions obtained from such a calculation provide a good description of the stress field measured numerically around plastic events, i.e. a quadrupolar stress redistribution. Using this theoretical representation of plastic events as Eshelby inclusions and considering the coupling between several events, predictions of the yield strain can be obtained [START_REF] Dasgupta | [END_REF]27], as well as of the inclination of shear bands [24,29]. Active rearrangements have also been studied numerically in order to investigate the dynamics of the response of the material to a local rearrangement [30] as well as the correlation between this elastic response and the subsequent plastic events [31].

It has to be noted that most of the numerical studies discussed above have been done in Lennard-Jones glasses for which a domain of truly elastic response exists. So, an elastic response of the matrix surrounding a small plastic event is not a surprise. For materials which do not behave as perfectly elastically, the pertinency of Eshelby's calculation can be questioned. Wu et al. show that the stress relaxation in liquid at times shorter than the Maxwell relaxation time may be described accordingly to Eshelby's stress tensors [32]. The fourfold pattern characteristic of Eshelby's stress tensors has observed in recent numerical studies of granular flows at constant volume [33]. Also, fluctuations of strains are observed experimentally [24] with inclinations related to the maximum of the deviatoric part of Eshelby's stress tensors. Similar oriented patterns were observed in numerical studies of compressed granular material [34], although the author does not link them to the Eshelby's tensors.

The aim of this study is to show that the stress redistribution due to a local reorganization in a granular material may be unambiguously described using the Eshelby formalism. For this, we focus on a realistic situation, where the granular material evolves at fixed pressure, allowing compaction or dilatation. We implement in a discrete element numerical model of a bidimensional frictional granular assembly an active inclusion with a new procedure very close to the spirit of the demonstration of Eshelby. We show that the observed response of the granular assembly corresponds to the one expected for an elastic material [START_REF] Eshelby | Proc. R. Soc. Lond. A[END_REF]. We study specifically the role of two parameters on the response of the system: the value of the friction coefficient, which governs the prevalency of sliding contacts in the system, and the loading state of the sample, i.e. its proximity to failure.

This study is organized in the following way:

In section II we explain the mode of loading of the material, which consists in a biaxial test, and we describe the preparation of the numerical experiments. We also describe the typical response obtained during the loading of the material. In section III, we give an overview of the principle of the calculation of Eshelby and we discuss our method of generating an Eshelby-type inclusion in our simulations. In section IV, we describe the stress field observed around an active inclusion at different stages of the loading and for different values of the friction coefficient. In section V we measure the angular distribution of the anisotropic response and we present the results of the evolution of the inclination of the maximum of the stress distribution as a function of the volumetric strain. We discuss the modification of the response of the system when the amount of plasticity in the granular matrix increases, either because of decreasing friction or because the system is closer to failure.

II. GEOMETRY AND PREPARATION OF THE NUMERICAL EXPERIMENT

In this section, we present the discrete element numerical simulations and explain the method of preparation of the samples. We consider a bidimensional sample submitted to a biaxial test. This classical configuration [START_REF] Lambe | Soil Mechanics[END_REF] consists in submitting a granular sample to well-controlled stress conditions without control of the volume. Experimentally, a uniaxial compression is exerted along one direction while a confining pressure is imposed on the lateral side of the granular sample. Those kind of tests are classically used in soil mechanics to study the strain-stress properties of granular materials [START_REF] Lambe | Soil Mechanics[END_REF]. For large enough initial volume fraction, localization of the deformation is observed above a threshold [START_REF] Desrues | [END_REF], with the formation of shear bands. Concerning the volumetric strain, still for initial volume fraction larger than a critical value, some compaction is observed at the beginning of the loading, followed by dilatancy at larger strain. Those kinds of test are unusual for physicists who prefer to work at constant volume and uncontrolled pressure, in order to study the behavior of a sample at a given volume fraction [37]. Yet, confining pressure controlled set-ups are the most natural configurations when trying to understand real life granular mechanical behavior. Because we are interested in this kind of configuration, we will give special care to local change of volume in the following. As a matter of fact dilatancy and plastic deformation are closely linked in granular materials.

In Sec. II A and II B, we discuss the numerical setup and test procedure. In Sec. II C we describe general prop- erties of our biaxial tests.

A. Numerical Setup

We perform numerical simulations of a twodimensional biaxial test.

The boundary conditions are sketched in Fig. 1. N grains of total mass M 0 are confined by four perpendicular walls. The left and bottom walls are fixed and motionless. The top wall moves downward with a constant velocity. These three walls can be considered to be of infinite mass, for their positions are not modified by the forces exerted by the grains. The right wall, however, is mobile, with a mass of M 0 /100. It is kept in place by a constant pressure p 0 .

The grains are polydisperse disks (r max /r min = 1.5) of uniform (two-dimensional) mass density ρ 0 . The quantities M 0 , ρ 0 , and p 0 establish the system of units. Thus length is measured in units of L 0 = M 0 /ρ 0 , time in units of T 0 = M 0 /p 0 , and velocity in units of V 0 = p 0 /ρ 0 .

The grains interact via traditional linear, damped springs in the normal and tangential directions (normal spring constant k N = 2000p 0 , tangential spring constant k T = k N /2). Weak linear damping is also included in these interactions. The grains also have a (microscopic) friction ratio µ m : the tangential force may not exceed µ m times the normal force. A weak rolling resistance is added so that all vibration modes will be damped. Finally, a weak gravitational force is added to gently push "rattlers" against the granular skeleton.

All four walls are frictionless, exerting only normal forces. The mobile wall has the same stiffness as the grains, but the three walls of infinite mass are soft (spring constant k wall = k N / √ N ) and dissipative. This novel boundary conditions efficiently removes vibrational energy from the system without resorting to global damping [START_REF] Mcnamara | III International Conference on Particle-based Methods: Fundamentals and Applications PARTICLES[END_REF].

B. Test procedure

The simulations presented below are done with N = 256 × 256 = 65536 grains. The initial state is formed by compressing a granular gas without friction (i.e. we set µ m = 0), yielding a relatively high density. Then, at the beginning of the test, friction is turned on (i.e. we set µ m to its final value), and the velocity of the top wall is set to 2 × 10 -5 H 0 /T 0 , where H 0 ≈ 1.1L 0 is the initial height of the system. We checked that decreasing further the velocity do not change the results of the simulation. As the simulation proceeds, the entire state of the system (particle positions and velocities, contact forces) are recorded at strain increments of ∆ε yy = 10 -5 , where the deformation ε yy is defined in II C

C. Global properties

In Fig. 2, we show the stress-strain curve for the simulations studied in this paper. The stress is deduced from the force exerted by the grains on the upper (constant velocity) wall. The strain is usually deduced in simulations from the position of this wall, but that does not work in this situation due to the softness of the wall. Indeed, more than half of the apparent deformation can occur at the wall, and not within the packing. Instead, we return to the definition of strain: ε yy = ∂u y /∂y, where u y is the y-component of the displacement of a material point relative to a reference state. We therefore calculate the displacement of each grain, and accumulate N pairs (y, u y ), and then do a linear regression and take the slope of as volumetric strain (%)

µ m =0.2 µ m =0.6 µ m =2
FIG. 3: Volumetric strain εxx + εyy for the tests studied in this paper.

∆ε yy , the deformation from the reference state. At the beginning of the simulation, the initial condition is the reference state, but each time the deformation increases by 10 -5 , the reference state is updated. The total deformation, averaged over the full sample, is ε yy in Fig. 2 and in the rest of the paper. Fig. 2 shows that we can control the peak stress by changing the friction ratio µ m . For µ m = 0.2 and µ m = 0.6, we obtain results similar to those obtained experimentally or in other numerical experiments. When µ m = 2, the peak stress exceeds 10 times the confining pressure. This quite high value do not correspond to experimental values of grain-grain friction coefficients. However, it is an interesting limit case for theoretical reasons.

In Fig. 3, we show the volumetric strain ε xx + ε yy during the tests. (The strain component ε xx is determined in the same way as ε yy strain component.) This result is in agreement with the experimental observations in biaxial and triaxial tests: the loading curve presents a maximum and the volumetric strain first decreases before continuously increasing.

III. ESHELBY'S INCLUSION: THEORY AND NUMERICAL IMPLEMENTATION

In this part we present our method for implementing an active inclusion in a granular material. Active inclusions in numerical simulation have been considered recently [30,31] in two-dimensional Lennard-Jones glasses by imposing displacements on the particles forming the inclusion. This approach amounts to imposing a strain on the inclusion. Here we propose a new procedure to implement an active inclusion where a stress is imposed. We do this for two reasons. First of all, we are dealing with granular materials, where changes of volume are very important. Thus we want to leave the inclusion free to choose the volume it wants. Second, the stress we impose appears at the beginning of the Eshelby's calculation [START_REF] Eshelby | Proc. R. Soc. Lond. A[END_REF]. We will explain this more precisely below.

Eshelby gave a general method based on the linear superposition of small deformation elasticity to solve a class of problems implicating an inclusion in an elastic matrix. He also gives analytical solutions of several typical problems. In the following, we will use the solution of this problem: consider a small part of an infinite elastic matrix which we will call the inclusion. Suppose that this inclusion has a spontaneous change of shape. The Eshelby solution of interest for us gives the stress field distribution in the matrix due to that local event when the unconstrained spontaneous change of shape of the inclusion is known.

The coordinate system used in the following is given in Fig. 1. For the vector positions, we use the polar coodinates r = r n = r cos θ sin θ . The origin is taken at the center of the inclusion. In all the considered cases, the imposed stress at the boundaries corresponds to biaxial stress-imposed configuration:

σ ∞ = σ ∞ xx 0 0 σ ∞ yy with σ ∞ xx and σ ∞ yy < 0, and σ ∞ xx -σ ∞ yy > 0.
In this theoretical section, we first recall the general principle of Eshelby's calculation. In the second part, we give the general form of the analytical solution far from the inclusion. In the last part our method of generating Eshelby-type inclusion in the numerical simulation is described.

A. Principle of the calculation

Consider an inclusion, within an infinite elastic matrix, which undergoes a spontaneous change of shape. If it had not been inside a matrix, this change of shape would correspond to the homogeneous strain tensor e * ij . The elastic stress field associated to this strain tensor is σ * ij = λe * kk δ ij + 2µe * ij , λ and µ being the Lamé coefficients [START_REF]Théorie de l'élasticité[END_REF]. As the inclusion is in a matrix, its change of shape will interact with the elastic surrounding medium, leading to a displacement field ũi ( r) in the matrix, generated by the rearrangement and given in 2D by: ũi

( r) = C G 2D ij ( r -r ′ )σ * jk dc ′ k (1)
where the integral covers the inclusion contour C, G 2D ij ( r -r ′ ) is the bidimensional Green function giving the displacement response awaited at r from a ponctual unitary force exerted in r ′ belonging to the contour on an element dc ′ k [START_REF] Eshelby | Proc. R. Soc. Lond. A[END_REF][START_REF]Théorie de l'élasticité[END_REF]. The following part discusses the far-field redistributed deviatoric stress for e * ij coaxial with the loading σ ∞ ij .

Complete solutions (with near-field terms) of the previous equations in the case of an elliptical inclusion can be found in [29].

B. Coaxial solutions

The displacement field in the matrix far from the inclusion can be approximated by (see Appendix A):

ũi ( r) = S 2πr λ + µ λ + 2µ e * jk g 2D ijk ( n) ( 2
)
where S is the area of the inclusion. The function g 2D ijk ( n) depends only on the direction n (see Eq. A3). The explicit determination of the solution depends on the expression of e * ij . In the present study, we restrict to the case of strain tensors characterizing the plastic rearrangement coaxial to the uniform imposed stress tensor σ ∞ , i.e.:

e * = e * xx 0 0 e * yy (3)
The geometry of the loading motivates this hypothesis. Furthermore, Dasgupta et al. [START_REF] Dasgupta | [END_REF]27] have shown that this orientation minimizes the energy of interaction between the external strain field and the inclusion at high enough strain. The components of the strain tensor ẽ induced by the rearrangement in the surrounding matrix can then be computed explicitly (see Appendix A).

In a biaxial experiment, the crucial quantity which characterizes the shear and which governs the response of the material is the deviatoric stress σ xx -σ yy . The significant component of the stress redistribution to determine is thus the redistributed deviatoric stress σxx -σ yy , which will add to the imposed stress field and will thus modify locally the total deviatoric stress. Here, we obtain from [START_REF]Statics and Kinematics of Granular Materials[END_REF] Figure 4(a) (resp. (b)) shows polar representations of the function f (θ) in a contracting (resp. dilating) case. We observe a dominant quadrupolar response with an increase of the deviatoric stress along directions close to 45 • from the principal stresses (red curve of Fig. 4). The change of volume induces a small modification of this inclination as well as an asymmetry in the negative part of deviatoric stress redistribution (blue curves of Fig. 4). This negative contribution counteract the macroscopic loading σ ∞ xx -σ ∞ yy .

(see Appendix A): σxx -σyy ∝ S πr 2 λ + µ λ + 2µ f (θ) (4)

C. Generating numerically an Eshelby inclusion

A careful reading of Sec. III A, especially of Eq. ( 1), shows that the stress redistribution due to the inclusion is obtained by integrating Green function for stress over the boundary of the inclusion. In the simulations, therefore, we will define a small inclusion, and then modify the stress at the boundary. But since the simulation is discrete, we must apply this stress in a discrete way.

We specify a perturbing stress σ * , and not a perturbing strain e * because in Eshelby's calculation, it is the stress that is propagated into the material [see Eq. ( 1)]. The perturbing stress is obtained from the perturbing strain, using the elastic coefficients. In our case, the elastic coefficients are not known, so that a more direct comparison can be obtained by directly imposing a stress.

Let us suppose that we want to estimate the stress in some region V. We have [28]

σ = 1 S α∈A F α ⊗ r α , (6)
where S = |V| is the area (in two-dimensions) of the region where we calculate the stress. The variable α labels the contacts and A is the set of contacts with at least one grain within V.

(A grain is in V if its center is in V).
The force between the grains is F α , r α is a vector along which the force is transmitted inside V, and ⊗ indicates a tensor (dyadic) product. Some contacts in A have both grains in V while others have only one -see Fig. 5. Let B be the set of all contacts with two grains in V while C = A\B is the set of all contacts with one grain in V. The vector r α depends on whether the contact α is in B or C. If α ∈ B (both grains in V), then r α = x jαx iα , where i α , j α label the two particles participating in contact α. On the other hand, if α ∈ C, then r α = x • αx iα , where x • α is the point where the line of centers crosses the boundary of V. We can therefore rewrite Eq. ( 6) as

σ = 1 S α∈B F α ⊗(x jα -x iα )+ 1 S α∈C F α ⊗(x • α -x iα ). ( 7
)
FIG. 5: Two contacts inside the region V where the stress is being calculated. Upper contact (α = 1): both grains lie within V, r1 is the line of centers. Lower contact (α = 2): only one grain lies within V, r2 extends from this grain, along the line of centers, to the boundary of V.

Instead of a sum over contacts, we can now write the stress as a sum over grains. We must define two new sets: let G be the set of all grains whose centers are in V, and let A i be the set of all contacts involving grain i.

Then we have

σ = 1 S i∈G α∈Ai χ iα F α ⊗ x i + 1 S α∈C F α ⊗ x • α . (8)
Here, χ iα = ±1, depending on whether we have i = i α or i = j α . Now the sum in brackets in Eq. ( 8) vanishes when the grains are in equilibrium, for it is simply the sum of forces exerted on grain i. The first term of Eq. ( 8) thus vanishes, leaving only those contacts that cross the boundary:

σ = 1 S α∈C F α ⊗ x • α . (9)
This equation tells us that the stress inside a region is determined just by the contacts that span the boundary.

Up to now, we have been considering these equations from a diagnostic point of view, that is, the contact forces and grain positions are extracted from numerical data, and then the stress tensor is calculated in view of passing over to a continuum description. Now we will take a different perspective: the stress tensor on the left hand side is a perturbing stress that we want to apply inside V, while the forces on the right hand side are modifications of the contact forces that we will introduce in order to apply that stress. Accordingly, we now write

σ * = 1 S α∈C F * α ⊗ x • α . ( 10
)
where we have written a star next to the stress tensor [as in Eq. ( 1)], and next to the forces, to indicate this new perspective.

Our method for generating an active inclusion consists of first choosing the stress σ * , and then choosing forces F * α so that Eq. ( 10) holds. The forces can then be added directly to the simuation.

We emphasize that our method begins with σ * , that appears in Eq. ( 1), at the beginning of Eshelby's calcution. The alternative of imposing displacements [30,31] corresponds to imposing ũi inside the inclusion. But Eshelby's calculation does not begin with ũi , but rather with σ * (or e * if the elastic moduli are known). Now let us return to the description of our method and consider the special case where V is a circle of radius R (S = πR 2 ), as shown in Fig. 5. The points x • α are all located on the circle. Placing the origin at the center of the circle, their coordinates can be written in terms of the R and the unit vector n α = cos θ α sin θ α . We have

σ * = 1 πR α∈C F * α ⊗ n α (11)
When we specify σ * , Eq. ( 11) contains four equations and 2|C| unknowns because each F * α is unknown and contains two components. The volume V is chosen to approximate the circular inclusion of the Eshelby calculation, and therefore must be larger than the grain size, leading to |C| ≫ 1. Thus the system ( 11) is always under-determined. To add extra conditions, we minimize α∈C |F * α | 2 . The meaning of this condition is that we look for the microscopic force distribution obeying Eq. ( 11) which disturbs the medium as little as possible. We therefore use Lagrange multipliers to minimize α∈C |F * α | 2 subject to the condition in Eq. (11). We thus minimize

L = α∈C |F * α | 2 -λ : πRσ * - α∈C F * α ⊗ n α . (12)
Here, λ is matrix of Lagrange multipliers, and the two points signify a total contraction between two matrices (A :

B = A ij B ij ).
Differentiating Eq. ( 12) by the forces, we get

F * α = 1 2 λn α . (13)
The Lagrange multipliers can be found by putting the forces in to the conditions Eq. ( 10):

2πrσ * = λ A C C B . ( 14
)
where A = α∈C cos 2 θ α , B = α∈C sin 2 θ α , C = α∈C cos θ α sin θ α . Before passing on to the results, we remark that a pair of grains can be considered as a contact in the set C even if the two grains are not touching. In that case the contact force vanishes, but that contact can be used to apply the perturbing stress.

IV. RESULTS OF THE NUMERICAL EXPERIMENT

In this section we study the response of the granular material to a provoked Eshelby inclusion.

A. Generation of a response to an Eshelby inclusion

We compute the response to an Eshelby inclusion in the following way. We first perform a compression of the granular material as described in II A up to a given deformation. We then identify a small circle region of radius R = L 0 /80 (approximately 6 grain radii) at the center of the sample. This small circle is the region V discussed above. For a region two times larger, the results presented below are unchanged. We also tested R = L 0 /160, and the main effect seems to be to rotate the far-field response by of order 10 • . If R is reduced another factor of two, there are not enough contacts to apply the perturbuation.

We then perform two simulations. In the first simulation, we simply continue the biaxial test without adding the perturbing force. The second simulation has exactly the same duration as the first, except that the perturbing forces are gradually "turned on" over a time long enough to avoid the generation of shock waves.

Then the stress σ i on each grain i is calculated as a sum over that grain's contacts. Using the notation of the previous section,

σ i = 1 πr 2 i α∈Ai χ iα (F α ⊗ r α ) r i |r α | . ( 15
)
where r i is the radius of grain i. This is an application of Eq. ( 6) to a single grain. The last factor arises because only part of the line joining the centers is inside the grain. Since we have two simulations, each grain has two stress values. Let σ (0) i be the stress on grain i in the unperturbed simulation, and σ

(1) i be the stress in the perturbed simulation. We then calculate σi = σ

(1) i -σ (0) i . (16)
Using this method, we calculate a perturbation stress for each grain that we can then use to visualize the effect of the inclusion.

B. Description of the stress response

We now focus on the spatial distribution of σyy,i -σxx,i inside the granular sample. In Fig. 6, we show the change in deviatoric stress σyy,i -σxx,i . The map is obtained by shading each grain accordingly to the value of the deviatoric stress. The perturbing stress applied in Fig. 6 is case of µ m = 2 simulation, the Eshelby pattern remains when the axial strain increases, but the angle of the bands increases. Meanwhile, the other two tests (µ m = 0.2 and 0.6) show a quite different behavior: the band angle remains near 45 • , but the texture of the image becomes much more grainy. We will discuss and interpret those different behaviors in section V.

C. Characterization of the angular response

We obtain quantitative information from Fig. 6 about the function f (θ). We divide the grains into classes according to their distance r from the center of V, and then into subclasses according to their angle θ. An average value of ∆σ i is calculated for each subclass. Then this average is multiplied by the r 2 . This allow us to test the 1/r 2 scaling of the redistributed stress. The resulting function for different classes of r is shown in Fig. 7. We observe a good collapse close to the perturbation (r ≤ 0.35L 0 ) showing that the far field approximation (leading to keep only the order 1/r 2 for the stress field) is sufficient to capture the radial response. Only the curves corresponding to the largest r differs significantly from the others, certainly because sampled region approaches the boundary of the simulation. (The simulation is approximately square with side length a bit larger than 1.1L 0 .)

The scaling in Fig. 7 is based on Eq. ( 4). This equation shows that the stress response should be proportional to S/r 2 , where S is the area of the inclusion (S = πR 2 in our case). Dividing by R 2 /r 2 removes this dependence. Futhermore, the response is proportional to the perturbing stress, so we divide by σ * yielding a function f (θ) of order unity.

We want to identify the position of the maxima of f (θ). To this end, we consider the grains with 0.1 < r/L 0 < 0.3, and sort them into classes according to their angle θ, and calculate the average change in stress in each class to obtain an estimate of f (θ). This function is then approximated by a Fourier series

f (θ) ≈ a 0 + n k=1 a k sin kθ + n k=1 b k cos kθ, (17)
with n 8. The maxima of this Fourier series are identified as the angle of the diagonal bands of the Eshelby.

We study the band angle as a function of loading. As the loading increases, a new problem appears: the packing becomes fragile, and the inclusion can trigger other events throughout the packing that swamp the change in stress caused directly by the inclusion. To get around this problem, we search for times where the packing is relatively stable (absence of kinetic energy peaks), and furthermore, we reduce the strain rate by a factor of 100, so that we will have time to "turn on" the forces in the inclusion.

In Fig. 8, we show the angle found in this way for the three biaxial tests of Figs. 2 and3. For each strain value, there are four angles, one for each of the diagonal bands giving its angle relative to the horizontal. These angles are averaged together to obtain a point on the graph. The error bars estimate the scatter of the points; they show the uncertainty of the mean. We observe that at the beginning of the loading all the angles are close to 45 • but with some fluctuations around this value. When the loading increase, the departure from 45 • becomes important. We will further discuss this figure in the next section.

To conclude, the mental construction of Eshelby underlying his computation can be implemented in a numerical simulation of dry frictional disks. In spite of the discrete nature of the material, the continuum mechanics solutions are clearly relevant as the observed response is very similar to the one predicted analytically. This picture is modified when the granular material if further loaded and the friction coefficient decreased.

V. DISCUSSION

In this part we discuss the departure from the quadrupolar response observed in Figure 6 either due to the value of the friction coefficient or to the loading of the material. As shown in the previous part, for an isotropically slightly pre-compressed granular material (initial preparation), the response of the material to a localized shear event is compatible with an elastic response (first column of Fig. 6). This picture is modified in two different ways when the material is submitted to an increasing deviatoric stress depending on the value of the friction coefficient. We discuss the different kind of departures from this ideal response in the following.

A. Beginning of the loading

As can be seen on the first column of figure 6, at the very beginning of the loading the response is very close to the theoretical prediction recalled in section III B and in Figure 4. Quantitatively, as discussed in the preceding section, the radial response is in 1/r 2 and the angular response is dominated by a quadrupolar dependence (Figure 7), with a maximum stress redistribution close to 45 • . The direct observation of the stress fields shown in figure 6 shows nevertheless some differences with the response of a continuous material: we can observe that a filamentary pattern underly the smooth response. Consequently, the observed stress redistribution is also compatible with a description of the granular material in term of force chains. The observed striations could indeed indicate that stress is concentrated on grains in force chains. Nevertheless, on a large scale, a smooth quadrupolar response is observed.

As the loading increases, two different behaviors are observed. For µ m = 2 (top row), the texture of the patterns remains the same, but the angles depart significantly from 45 • . For µ m = 0.6 and µ m = 0.2, the angles remain close to 45 • , but the texture changes significantly. We now discuss these changes.

B. Further loading: High value of µm

When the friction coefficient is artificially high to prevent sliding contacts (upper line of Figure 6), we observe a clear modification of the angles characterizing the stress resdistribution. Meanwhile, the increase in the fluctuations is much smaller than for the smaller friction coeffi-cients. A change of inclination can still be interpreted in an Eshelby picture and we will discuss two possible origins in the following: the effect of the volumetric strain and the effect of anisotropical elasticity.

Polar response and volumetric strain

As discussed in section III B, the quadrupolar response due to the shear part of the rearrangement is modified by the volumetric contribution. To make explicit this dependence, we study the polar function f (θ) of Eq. 5. The minima of the function are obtained for θ = 0

For the isovolumic case (e * xx + e * yy = 0), we obtain

θ E = 45 • [90 • ].
Our active inclusions correspond to an imposed traceless stress tensor Trace(σ * ) = 0. If the inclusion is made of an isotropic elastic material, the associated strain tensor will also be traceless. But the loading process introduces anisotropy into the contact network and thus into the elastic constants, as will be discussed in the next section. In such cases, Trace(e * ) may be non-zero.

If we suppose that the inclusion changes its volume, a reasonable hypothesis is that e * xx > 0 and e * yy < 0, i.e. that the local deformation follows the overall macroscopic deformation. We then obtain for the range of possible values θ E ∈ [ xx + e * yy > 0, i.e. to dilating plastic events, while the values lesser than 45 • correspond to contracting plastic events. Those two cases are shown on Fig. 4: the angle of the positive part of the deviatoric stress redistribution is indeed smaller than 45 • in the case of a contracting event (Fig. 4(a)) and larger in the case of a dilating one (Fig. 4(b)). The range of possible angles found here correspond to the direction of the maximum of the redistributed deviatoric stress from a single Eshelby inclusion. It differs from the range obtained in [29], where the range of directions of coupled Eshelby inclusion is determined by minimizing the elastic energy of interaction between the inclusions.

The result of this calculation is that only small departure from 45 • can be explained by non-isovolumic transformation. It could for example explain the fluctuations measured in Figure 8 for µ m = 0.2 and µ m = 0.6 but not the large increase observed in the µ m = 2 case.

Anisotropy of the elastic matrix

In this part, we will neglect possible change of volume of the inclusion (Trace(e * ) = 0) but we will consider the Eshelby problem in the case of anisotropic elasticity. Indeed, the loading process may induce some anisotropy of the elasticity of the material [START_REF] Mcnamara | [END_REF]41]. The computation of the Green's function in an anisotropic material is done in Appendix B. Equation B1 shows that in the case of an orthorhombic bidimensional material with a small anisotropy, terms in 2θ and 6θ appear. On the contrary, change of volume of the inclusion do not generate 6θ terms. So the presence of 6θ terms in f (θ) may be viewed as an indication of anisotropic elasticity. This can be checked by the Fourier analysis done in section IV C (fit of the measured response by equation 17). The different components of the fit are shown in Figure 9. We observe that indeed the 6θ component is the second most important harmonic in the case of µ m = 2. This supports the hypothesis that the observed modification of the angle is due to the effect of material anisotropy, and not to volume change of the inclusion. This effect is less clear for µ m = 0.6 and µ m = 0.2 where the second dominant harmonic is 2θ which can originate both from volumetric effects and anisotropical effects. In fact, for those values of µ m , the principal effect of the loading is not modification of the angles (that fluctuate close to 45 • -see Figure 8) but rather the increase of fluctuations (see Figure 6) as will be discussed in the next part.

C. Further loading: Low value of µm

For µ m = 0.6 and µ m = 0.2 at intermediate and large values of the strain (lower right four panels in Fig. 6) the filamentary striations visible in the other panels of Fig. 6) are replaced by a grainy strongly fluctuating texture. In spite of this, the response at large scales remains quadrupolar, with angles close to 45 • . This change in texture cannot be described in the Eshelby formalism, which supposes that the material is continuous. We do not think that this difference of texture is due to the disorder of the contacts between grains, but likely due to sliding contacts. Indeed, the µ m = 0.2 and µ m = 0.6, there are many sliding contacts (up to 35% and 15% of all contacts, respectively) but at µ m = 2, sliding contacts remain relatively rare (less than 2.5% of all contacts).

The permanence of the quadrupolar structure at large scale shows than an elastic description of a granular material can co-exist with different mesoscopic behaviors: either force chains (ε yy = 10 -5 or µ m = 2) or sliding contacts (large ε yy and µ m = 0.2, 0.6). When many contacts are sliding, the granular matrix does not behave as an elastic material. We then expect a stress redistribution which differs from the elastic Eshelby response. However, at large scales, we still see in Fig. 6 the Eshelby structure which is a characteristic of the material elasticity. It is unclear if it is possible to quantify objectively such a separation of scales.

VI. CONCLUSION

In this paper, we considered the effect of a local transformation on the mechanical stress distribution. For this, we applied the theoretical construction proposed by Eshelby to numerical simulations where extra forces are applied on some grains inside a granular material. In response to this local applied stress, the material relaxes stress at large distance. This stress relaxation may be, at least in some cases, successfully described by a theory where the granular material is treated as a continuous and elastic material. The presence of disorder of contact forces (such as the so called force-chains) does not disrupt this elastic response, at least at the macroscopic scale. If sliding contacts are supressed, we show that the material remains elastic at the macroscopic scale, even close to failure, but the elastic coefficients become anisotropic.

Our numerical simulations are performed with beads of finite stiffness. The ratio p 0 /k N = 5.10 -4 between the confining pressure and the contact stiffness, is the typical value of the relative deformations of the beads. It is interesting to compare it to a real 3D experiment. If we consider glass beads following Hertz law which are confined under a pressure P = 100 kP a, we may estimate a relative deformation of the beads of 3.10 -4 . So the value of the relative deformation is of the same order of magnitude as many biaxial compression experiments, such as experiments where stress relaxations along preferential directions are observed [23,24]. Some departures from the elastic behavior are clearly visible in our simulations. They may be linked to the occurrence of sliding contacts. The departure from an elastic response of the matrix is an interesting, but very challenging, problem. We may expect that the granular skeleton behaves elastically for deformations of the inclusion which are smaller than some plasticity threshold ε th . The variation of ε th with the microscopic friction coefficient µ m , with the reduced stiffness of the beads p 0 /k N , and with the proximity of the material failure are still unknown. This will be the subject of a future work.

Appendix A: Computation of the released stress

The bidimensional Green function is:

G (2D) ij ( r -r ′ ) = 1 4πµ λ + µ λ + 2µ - λ + 3µ λ + µ δ ij ln | r -r ′ | + (x i -x ′ i )(x j -x ′ j ) | r -r ′ | 2 (A1)
where r = x y = r cos θ r sin θ are the coordinates in the bidimensional plane.

The displacement field in the matrix can be obtained by integrating the Green function along the contour C of the inclusion: We consider the case of an orthorhombic bidimensional material which free energy is of the form [START_REF]Théorie de l'élasticité[END_REF]: To find the Green's function of this system, i.e. the solution of the system when b i = n i δ( r) where n is a vector of unit norm, we work in the Fourier space. Noting the Fourier transform of the components of the Green tensor G ij , we obtain the following relations:

ũi ( r) = 1 4πµ λ + µ λ + 2µ C - λ + 3µ λ + µ δ ij ln | r -r ′ | + (x i -x ′ i )(x j -x ′ j ) | r -r ′ | 2 σ * jk dc ′ k = σ * jk 4πµ λ + µ λ + 2µ S λ + 3µ λ + µ n k δ ij + (-δ ik n j -δ jk n i + 2n i n j n k ) dS ′ | r -r ′ | Where n = r-r ′ | r-r ′ | . Using σ * ij = λe * ll δ ij + 2µe
F = 1 2
λ xxxx q 2
x G xx + λ xyxy q 2 y G xx + (λ xxyy + λ xyxy )q x q y G yx = 1 4π 2 λ xxxx q 2

x G xy + λ xyxy q 2 y G xy + (λ xxyy + λ xyxy )q x q y G yy = 0 λ yyyy q 2 y G yy + λ xyxy q 2 x G yy + (λ xxyy + λ xyxy )q x q y G xy = 1 4π 2 λ yyyy q 2 y G yx + λ xyxy q 2 x G yx + (λ xxyy + λ xyxy )q x q y G xx = 0

The components G yx can be deduced from those relations. We compute them in the case of small anisotropy, i.e. when λ xxxx -λ xxyy -2λ xyxy and λ yyyy -λ xxyy -2λ xyxy are small [START_REF]Théorie de l'élasticité[END_REF]. We then write:

λ xxxx = λ + 2µ -ǫ λ yyyy = λ + 2µ + ǫ λ xxyy = λ λ xyxy = µ
with ǫ ≪ λ, µ. We obtain then

G xx = 1 4π 2 µ 1 q 2 - λ + µ λ + 2µ q 2
x q 4 1 + ǫ λ + 2µ q 2 x q 2 1 + (λ + µ)q 2 y (λ + µ)q 2 y + µq 2 G yy = 1 4π 2 µ 1 q 2 -λ + µ λ + 2µ q 2 y q 4 1 -ǫ λ + 2µ q 2 y q 2 1 + (λ + µ)q 2 x (λ + µ)q 2

x + µq 2 G xy = G yx = -1 4π 2 µ λ + µ λ + 2µ q x q y q 4 1 + ǫ λ + 2µ q 2 x -q 2 y q 2

We have ê * il = 4π 2 Sq k q l G ij σ * jk [16] with S the surface of the inclusion, from which we deduce:

êxx -êyy = -Sσ * µ 1 -λ + µ λ + 2µ (q 2 x -q 2 y ) 2 q 4 + ǫ λ + 2µ q 2 x -q 2 y q 2 -λ + µ λ + 2µ (q 2

x -q 2 y ) 3 q 6 (B1)

for σ * = σ * 0 0 -σ * .

The first isotropical term on the right hand side of eq. B1 Fourier transform leads to a Dirac term in r = 0. It is in fact due to the ponctual plastic event in the inclusion and should be removed when computing the elastic response [16]. The (q 2 x -q 2 y ) 2 q 4 leads to a cos 4θ r 2 term in the real space [16]. It differs from the usual form q 2 x q 2 y q 4 obtained in simple shear configuration because the pattern is rotate of 45 • in a biaxial loading compare to simple shear. We thus recover the quadrupolar response for an isovolumic inclusion in absence of anisotropy (ǫ = 0). The term of order ǫ gives the modification of the response due to the anisotropy of the matrix. From symetry arguments, we see that those terms will give rise in the real space to terms of the form cos 2θ r 2 and cos 6θ r 2 .

FIG. 2 :

 2 FIG.2: Stress strain curves for the simulations analyzed in this paper.

  with f (θ) = -(e * xx + e * yy ) cos(2θ) -(e * xx -e * yy ) cos(4θ) (5) This function has a quadrupolar part which originates from the shear part of the strain: (e * xx -e * yy ) and a bipolar part due to the volumetric strain (e * xx + e * yy ). The function f (θ) is independent of the elasticity coefficients because we are in two dimensions. In three dimension the corresponding function depends on the Poisson's ratio.

FIG. 4 :

 4 FIG. 4: Polar representation of the function f (θ) (Eq. 5). (a) Contractant event: e * xx = .3 and e * yy = -.7. (b) Dilatant event: e * xx = .7 and e * yy = -.3.

FIG. 6 :

 6 FIG. 6: Examples of generated Eshelbies, N = 256 × 256 grains. Top row: µm = 2, middle row: µm = 0.6, bottom row µm = 0.2. σ * is coaxial with the load; white grains correspond to a reduction in deviatoric stress. The value of the global strain is given under each panel.

2 FIG. 8 :

 28 FIG. 8: Angle of diagonal bands for simulations with different friction coefficients. For each value of the strain, four angles are averaged together to obtain a mean. The error bars give an idea of the scatter of the different angles measured. The horizontal line indicates an angle of 45 • .

FIG. 9 :

 9 FIG. 9: Harmonic components obtained from the fits obtained using equation 17. (a) µm = 2, (b) µm = 0.6 and (c) µm = 0.2.

  • ; 52.2 • ]. The values greater than 45 • correspond to e *

	[37.8	1 2 acos 1 4 ; 1 2 acos -1 4 ], i.e. θ E ∈

  (δ ij n k + δ ik n j -δ jk n i ) + 2n i n j n k dS ′ 2D ijk ( n) = µ λ+µ (δ ij n k + δ ik n j -δ jk n i ) + 2n i n j n k .For the particular case when e * = e * ẽyy = cos(2θ) (ẽ rr -ẽθθ ) -2 sin(2θ)ẽ rθ , we obtain for the deviatoric strain:

	using ẽxx -ẽxx -ẽyy =	S πr 2	λ + µ λ + 2µ	-(e * xx -e * yy ) cos 4θ -(e * xx + e * yy ) cos 2θ	(A4)
							* ij one obtains:
	ũi ( r) = Consequently, far from the inclusion, we can write: e * jk 2π λ + µ λ + 2µ S µ λ + µ	| r -r ′ |	(A2)
							ũi ( r) =	S 2πr	λ + µ λ + 2µ	e * jk g 2D ijk ( n)	(A3)
	with g xx 0 e * 0 yy	we have:
	ũr ( r) =	S 2πr	λ + µ λ + 2µ	λ + 2µ λ + µ	(e * xx -e * yy ) cos 2θ + (e * xx + e * yy )
	ũθ ( r) =	S 2πr	λ + µ λ + 2µ	-	µ λ + µ	(e * xx -e * yy ) sin 2θ
	The components of the strain tensor in polar coordinates are then:
	ẽrr =	S 2πr 2	λ + µ λ + 2µ	-	λ + 2µ λ + µ	(e * xx -e * yy ) cos 2θ -(e * xx + e * yy )
	ẽθθ =	S 2πr 2	λ + µ λ + 2µ	λ λ + µ	(e * xx -e * yy ) cos 2θ + (e * xx + e * yy )
	ẽrθ =	S 2πr 2	λ + µ λ + 2µ	-(e * xx -e * yy ) sin 2θ

Appendix B: Eshelby inclusion in an anisotropic material

  λ xxxx e 2 xx + 1 2 λ yyyy e 2 yy + λ xxyy e xx e yy + 2λ xyxy e 2 xy From the equilibrium equations ∂σ ik ∂x k + b i = 0, where b i is a body force, we can then deduce λ xxxx ∂ 2 u x ∂x 2 + λ xyxy ∂ 2 u x ∂y 2 + (λ xxyy + λ xyxy ) ∂ 2 u y ∂x∂y + b x = 0 λ yyyy ∂ 2 u y ∂y 2 + λ xyxy ∂ 2 u y ∂x 2 + (λ xxyy + λ xyxy )

	∂ 2 u x ∂x∂y	+ b y = 0