Gilles Marckmann
email: gilles.marckmann@ec-nantes.fr

Grégory Chagnon

Matthieu Le Saux

Pierre Charrier

Experimental investigation and theoretical modelling of induced anisotropy during stress-softening of rubber

Keywords: Mullins effect, stress softening, constitutive equation, network alteration theory, anisotropy

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Industrial competition in automotive industry has led to a growing use of numerical simulation in the design process of elastomeric parts such as tires, seals, engine mounts, shock absorbers... A good understanding and a well-established modelling of the material behaviour are essential prerequisites for efficient numerical modelling as Finite Element analysis. Elastomers are known to exhibit a strongly non-linear elastic behaviour under static loading conditions, characterised by the well-known S-shaped uniaxial response. Inelastic features such as hysteresis, stress-softening and residual strain are observed under cyclic loading. The strain-induced stress-softening phenomenon, widely known as the Mullins effect [START_REF] Mullins | Effect of stretching on the properties of rubber[END_REF][START_REF] Mullins | Softening of rubber by deformation[END_REF], is characterised by a reduction of stiffness during the first loading cycles. Moreover, the material response is perceptibly more compliant at strain smaller than the maximum strain endured during this first cycle, as pointed out by many authors [START_REF] Mullins | Effect of stretching on the properties of rubber[END_REF][START_REF] Bonart | X-ray investigations concerning the physical structure of cross-linkingin segmented urethane elastomers[END_REF][START_REF] James | Strain energy functions of rubber. II-The characterization of filled vulcanizates[END_REF][START_REF] Muhr | Experimental determination of model for liquid silicone rubber: Hyperelasticity and Mullins effect[END_REF][START_REF] Robisson | Comportement mécanique d'un élastomère chargé en silice. Etude de l'influence des charges et modélisation par une loi viscohyperélastique endommageable[END_REF][START_REF] Pawelski | Softening behaviour of elastomeric media after loading in changing directions. Constitutive models for rubber[END_REF][START_REF] Chagnon | Modélisation de l'effet Mullins dans les élastomères[END_REF][START_REF] Dorfmann | A constitutive model for the Mullins effect with changes in material symmetry[END_REF][START_REF] Machado | Induced anisotropy by the mullins effect in filled silicone rubber[END_REF][START_REF] Merckel | A Mullins softening criterion for general loading conditions[END_REF][START_REF] Rebouah | Permanent set and stress softening constitutive equation applied to rubber like materials and soft tissues[END_REF], the level of softening is not identical in all directions and present induced anisotropy. Moreover, rubber-like materials subjected to cycling loading conditions do not return to their initial state when unloaded, and residual strain is generally observed after unloading. The Mullins effect does not find unanimous physical explanation at present [START_REF] Diani | A review on the Mullins effect[END_REF] and occurs both in unfilled and filled elastomers, though its influence is more pronounced in filled rubber-like materials. To describe the strain-induced anisotropy of the Mullins effect, different approaches were developed. [START_REF] Besdo | Directional sensitivity of Mullins effect constitutive models for rubber[END_REF] introduced a tensorial history function based on second-order tensors evaluated using ellipsoid representations of all right Cauchy-Green tensors previously encountered during the deformation process. [START_REF] Shariff | An anisotropic model of the Mullins effect[END_REF] proposed a phenomenological damage anisotropic model based on the assumption that stress-softening is related to principal stretches, and to maximum and minimum strain values of the principaldirection line elements encountered during the deformation history. [START_REF] Göktepe | A micro-macro approach to rubber-like materials. Part III: The micro-sphere model of anisotropic Mullins-type damage[END_REF] investigated the molecular damage concept developed by [START_REF] Govindjee | A micro-mechanically continuum damage model for carbon black filled rubbers incorporating Mullins' effect[END_REF], and decomposed the network structure into a particle-to-particle and a crosslink-to-crosslink part. They exclusively associated the anisotropic stress-softening phenomenon to the particle-to-particle contribution, which involves breakage of bonds between chains and filler particles, and associated a scalar damage parameter with each direction of the hyperelastic micro-sphere model [START_REF] Miehe | A micro-macro approach to rubberlike materials -Part I: the non-affine micro-sphere model of rubber elasticity[END_REF]. Thereby, the model inherently reproduces the strain-induced anisotropic distribution of the Mullins effect and naturally includes residual strain. Diani et al. (2006b) combined the [START_REF] Pawelski | Softening behaviour of elastomeric media after loading in changing directions. Constitutive models for rubber[END_REF] directional molecular-based hyperelastic model with the [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF] alteration theory and derived an anisotropic model for the Mullins effect and the inherent induced residual strain. [START_REF] Dargazany | A network evolution model for the anisotropic Mullins effect in carbon black filled rubbers[END_REF] proposed a similar approach by taking into account the existence of different chains with different lengths in each direction. In a similar approach, [START_REF] Itskov | A thermodynamically consistent phenomenological model of the anisotropic Mullins effect[END_REF] proposed three damage evolution functions for the three principal stretches. [START_REF] Dorfmann | A constitutive model for the Mullins effect with changes in material symmetry[END_REF] used a pseudo elastic approach to take into account the stress softening and changes in material symmetry. [START_REF] Merckel | Constitutive modeling of the anisotropic behavior of Mullins softened filled rubbers[END_REF] used the strain amplification factor concept introduced in the microsphere approach. [START_REF] Rebouah | Anisotropic Mullins softening of a deformed silicone holey plate[END_REF] and [START_REF] Machado | Theory and identification of a constitutive model of induced anisotropy by the Mullins effect[END_REF] proposed to develop a constitutive equation that depends both on the maximum and actual deformation in a full network approach.

A residual strain is generally observed after removal of the stress. The magnitude of this residual strain increases with the amount of fillers in rubber and with the maximum stretch endured before unloading [START_REF] Dorfmann | A constitutive model for the Mullins effect with permanent set in particule-reinforced rubber[END_REF]. To model both stress-softening and residual stretch, [START_REF] Ogden | A pseudo-elastic model for the Mullins effect in filled rubber[END_REF] used the previous pseudo-elastic approach to describe the Mullins effect, and incorporated a single additional damage variable function of the maximum strain previously achieved and for which no restriction is henceforth imposed. [START_REF] Holzapfel | Aspects of stress softening in filled rubbers incorporating residual strains[END_REF] extended this model to the inclusion of an anisotropic description of the damage mechanism governed by three damage variables related to the evolution of principal stretches. [START_REF] Dorfmann | A constitutive model for the Mullins effect with permanent set in particule-reinforced rubber[END_REF] also investigated the theory of pseudo-elasticity and developed a model combining stress-softening and residual strain by including two additional parameters respectively related to these inelastic phenomena. The authors pointed out that the anisotropic material response under uniaxial deformation (transversal isotropy) may be described by a strain energy function which is not symmetrical in terms of principal stretches (i.e. which does not satisfy the well-known [START_REF] Valanis | The strain-energy function of a hyperelastic material in terms of the extension ratios[END_REF] assumption). They as-sociated the induced anisotropy only with the residual strain variable but specified that a more general model, where the induced anisotropy is related to the softening variable, can be derived. Recently, [START_REF] Rebouah | Permanent set and stress softening constitutive equation applied to rubber like materials and soft tissues[END_REF] proposed a directional model were the stress becomes zero in some directions before the zero deformation is obtained to generate residual elongation.

The aim of the present study is to experimentally exhibit and model the evolution of stiffness of rubber in different directions depending on the history and the directions of first stretching. The residual strain and the loss of isotropy inside the material will be analysed through a large panel of loading conditions. In this way, we propose a phenomenological constitutive model established under micro-mechanical considerations and that principally accounts for the anisotropy induced by the Mullins effect. The first key contribution is concerned with the investigation of a set of original appropriate homogeneous experiments conducted under various deformation conditions.

Results obtained lead to a thorough understanding of the anisotropic nature of stress-softening and the residual strain and highlight the occurrence of a stiffening effect along contracted directions (stretch lower than 1) beyond the pre-softening stretch level. Pursuing the network alteration theory proposed by [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF], these inelastic phenomena are associated with micro-structural alteration mechanisms involving slippage and breakage of links between polymer chains. This rearrangement is reflected by the increase of chain contour length in tensile stretched regions and its slow decrease along contracted directions.

The paper is organised as follows. Section 2 reports experimental inves-tigations emphasing anisotropy of the material properties which takes place under homogeneous uniaxial simple tests: tensile test, pure shear, modified equibiaxial and "biaxial" extension. Section 3 is devoted to the derivation of a phenomenological model, based on micro-mechanical considerations, that incorporates the anisotropic evolution of stress-softening and stiffening effects. Section 4 is dedicated to the evaluation of predictive capabilities of the model. The model is fitted on experimental results presented in the paper.

Finally, some concluding remarks are given in Section 5.

Experimental observations

An important feature of the stress-softening effect is the occurrence of anisotropic material properties under stretching. This property has hitherto received little attention in the literature. A thorough understanding of this phenomenon and the evaluation of the anisotropic model proposed in the following require appropriate experiments. Within this framework, several quasi-static pure homogeneous experiments were performed according to different loading conditions. These experiments were conducted under controlled displacement conditions at a constant strain rate of 0.04167 s -1 by the French Research Department of the Trelleborg group. The material used in this work is filler-reinforced natural rubber with a mass percentage of carbon black of approximately 30%. The detailed composition of the material is given in Le [START_REF] Le Cam | Mechanism of fatigue crack growth in carbon black filled natural rubber[END_REF]. Note that stress values in the present paper are normalised. All experiments are repeated on four lots of samples.

All curves presented in this paper are the mean responses of the specimens.

Simple tension tests of previously softened specimens

The observation of the influence of anisotropic effects requires mechanical tests involving different stretching directions. One of them consists first in performing five pre-softening uniaxial tensile cycles up to a stretch equal to 3.25 on a virgin large flat dumbbell sample extracted from a calendered plate.

The later sample is 110 mm long and 40 mm large in the useful area. The thickness is equal to 2 mm. Smaller dumbbell samples are then cut from this large specimen with different orientations with respect to the direction of pre-softening stretch (reference direction)1 . The later are 10 mm long and 2.5 mm large in the useful area. Fig. 1(a) shows how the small dumbbell samples are extracted from the large one. Then, they are submitted to tensile stretch loading until breakage.

Small samples are also cut in a virgin large specimen and similar tensile tests were performed in order to evaluate the initial anisotropy of the plate. These experimental results exhibit an appreciable stiffer mechanical response when approaching the calendering direction.

As we only focus on the anisotropy induced by Mullins effect and not by the calendering process, the initial anisotropy observed above is taken into account by normalizing the mechanical response of specimens for each direction. The calendering direction is retained as the reference response (0 •direction) and a corrective factor is calculated for each angle. This factor is evaluated so that all virgin curves are superimposed with the reference one (cf. Table 1).

We observe that this factor is linear with the angle. Then, mechanical responses of the pre-softened small samples are corrected by using the aforementioned corrective factor depending on the extraction angle.

Note furthermore that, due to residual strain, the dimensions of the small specimen are necessarily measured in a pre-softened configuration. approach is used in the numerical experiments detailed in Section 3. Thus, the unstrained virgin material configuration is adopted as the reference one.

Reference dimensions of the cut samples are determined with respect to this configuration. Thus, they are determined from the ratio of the small samples initial dimensions in the pre-softened configuration to the residual stretches in the corresponding principal directions of the secondary tensile deformation. Considering residual stretches, defined by the residual deformation gradient tensor F res , and the principal basis of deformation {e i } i=1...3 , residual stretches are defined in a new basis, e * i = Re i , which is the strain principal directions of the cut samples during the succeeding stretches (see ig. 3). and where R is the rotation tensor between e and e * . Then, the new reference deformation gradient tensor is given by:

F * res = R T F res R (1)
Fig. 4 shows stress-stretch responses of secondary specimens extracted at different angles (0 • , 22.5 • , 45 • , 67.5 • and 90 • ) from the reference direction.

The nominal stress and stretch are both expressed in the virgin unstrained configuration and the stretch-stress curves are corrected according to the initial anisotropy of the virgin material as described above. Mechanical responses of the pre-deformed specimens appear to be lower than the virgin material one. This reflects the strain-induced stress-softening phenomenon.

The anisotropic distribution of stress-stretch properties is also clearly ob- servable. Although the prior tensile stretch produced stress-softening in all selected directions, the softening level significantly decreases as the direction deviates from the one of the initial stretch. It also appears to be nearly null in the orthogonal direction (90 • in the figure). Furthermore, all secondary specimens are conditioned to the same stretch ratio. The material gets stiffer when the stretch level approaches the maximum stretch endured during the pre-softening stretch level (i.e. λ = 3.25). The material does not completely recover its virgin behaviour beyond the pre-softening stretch level (cf. table 2). Indeed, due to residual strain, the secondary loading curves should rejoin the virgin curve for stretch levels slightly higher than the maximal previously angle of extraction ( endured one [START_REF] Johnson | The Mullins effect in uniaxial extension and its influence on transverse vibration of rubber string[END_REF]. As [START_REF] Muhr | Experimental determination of model for liquid silicone rubber: Hyperelasticity and Mullins effect[END_REF] pointed out, the softening generally remains above the prior maximum endured stretch level [START_REF] Diani | A review on the Mullins effect[END_REF]. The authors specified that application of the residual strain correction, usually adopted for the idealisation of the Mullins effect, would reduce this drop below the virgin curve.

0 • 22.5 • 45 • 67.5 •
Mechanical responses of the 0 • , 22.5 • and 45 • -small samples exhibit a slope close to the virgin curve one beyond the maximum stretch previously reached. On the other hand, relatively to the virgin material response, significant stiffening appears for angles of 67.5 • and 90 • when the pre-softening stretch level is exceeded. The latest samples have submitted stretches lower than 1 during the pre-softening cyclic test. Consequently, it seems that stresssoftening phenomenon is coupled with a stiffening effect in these directions.

Even if the observed effect is rather small in comparison with the experimental dispersal, the same phenomenon is observed for each lot of samples.

We may presume that stiffening is linked to the compression loading conditions. This assumption is consistent with the common perception of the Mullins effect which is generally ascribed to the quasi-irreversible rearrange-ment of the molecular network under stretching. In that case, the material is assumed to have a selective memory of the maximum stretch ratio encountered during its loading history and to adapt a new equilibrium configuration softer than the virgin one as long as this maximum stretch is not exceeded.

Finally note that, in conformity with the incompressibility of the material, samples cut at 0 • and 22.5 • reveal a residual stretch greater than 1, those oriented at 67.5 • and 90 • manifest a contraction residual stretch (< 1) whilst the specimen extracted at 45 • does not show any.

Pure shear tests from previously softened specimens

Similar experiments are conducted with a pure shear pre-softening test.

Here, the specimen used for the pure shear test is moulded and is assumed to be initially homogeneous and isotropic unlike the calendered plate. The test is clearly as lower as the direction deviates from the pure shear loading direction. Observations similar to those relative to the previous uniaxial tensile tests can be made about the recovering of the virgin material behaviour beyond the pre-softening stretch level (λ = 4). However, a stiffening effect which occurs in directions close to the direction orthogonal to the presoftening uniaxial tensile stretch loading is not clearly observable here. Since the material is not stretched nor contracted in the direction orthogonal to the pure shear loading one, this observation coincides with the previous analysis where the stiffening phenomenon appears when stretches are lower than 1.

The smaller the angle is (in comparison with the reference one), the greater the amount of residual stretch is. The sample oriented to 90 • does not exhibit any residual strain.

Modified biaxial tension tests from previously softened specimens

Tests inspired by the so-called modified biaxial tension test proposed by [START_REF] Pawelski | Softening behaviour of elastomeric media after loading in changing directions. Constitutive models for rubber[END_REF] ones used for the uniaxial tensile test of section 2.1. A correction is performed to take into account the initial anisotropy of the plates. The responses of the material are reported in Fig. 8. As expected, the 45 • -secondary sample exhibits the smaller amount of softening since the maximum pre-stretch endured in this direction is the lowest one. However, a remaining significant softening occurs beyond the pre-softening stretch level. The 22.5 • and 0 •samples are slightly stiffer than the virgin sample beyond the prior maximum stretch ratio, while 67.5 • and 90 • -samples exhibit a significant drop below the virgin curve. Consequently, it seems that, under the maximum stretch ratio endured during the pre-softening test, the succession of compressive and tensile loading affects the material response. Indeed, due to the testing procedure, this experiment differs from a pure equi-biaxial test. Load histories in the two pre-stretching directions are different. The 22.5 • and 0 • -directions successively experience contraction followed by tensile stretch, whilst the 67.5 • and 90 • -directions undergo tensile stretch followed by compressive one.

These observations confirm the previously assumed existence of a directional chain network rearrangement which takes place during stretching. Mechanisms involved in this process may generate a new stress-free equilibrium configuration stiffer than the virgin one after compression.

An analogous test is conducted by applying a pre-stretching of 3.25 in the first direction of the cross-shaped specimen, perpendicular to the calendering one, whereas the reference direction is stretched up to a factor 2.25.

Responses are depicted in Fig. 9. The respective contributions of the two pre-softening contributions, as well as the induced anisotropy, are clearly distinguishable. Remarks similar to those enunciated above may be stated concerning the influence of the loading history undergone in each material direction.

Conclusions

All the experimental observations made in the present section clearly highlight the anisotropic nature of the Mullins effect. The angular distribution of stress-softening occurring under various homogeneous deformation modes apparently depends on the maximum stretches undergone by the material. However, an induced stiffening phenomenon along contracted directions was exhibited and the influence of the loading history is pointed out. This large variety of experimental results will be used to validate the model derived in the next section.

Constitutive model

The aim of this section is to build a constitutive model for rubber which is able to take into account anisotropy induced by pre-stretches. Mathematical symbols are defined at the end of the paper.

The previous experimental observations are consistent with a network alteration as proposed by [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF] where alteration depends on the stretch endured by chains and thus on the direction of stretching.

Description of the model

Rubber elasticity can be described through micro-mechanically based models. As they are explicitly expressed in terms of these two parameters, the chain density n and the average number of monomers per chains N, micromechanical network models based on non-Gaussian chain statistics are well-adapted to the use of stress-softening.

In the paper, the affine "Micro-sphere" model of [START_REF] Miehe | A micro-macro approach to rubberlike materials -Part I: the non-affine micro-sphere model of rubber elasticity[END_REF] is chosen as the reference elastic model with a set of 42 directions. For convenience, denote

λ c = λ (d c ) = F • d c = √ d c T • C • d c the stretch of a single
chain oriented in the direction of the unit vector d c . The representation of the macroscopic Cauchy stress tensor under the incompressibility assumption (det F = 1) is given by : Besides the alteration processes previously stated by Marckmann et al.

σ = C R √ N r c=1 ω c λ c L -1 λ c √ N (F • d c ) ⊗ (F • d c ) -pI (
(2002), we can presume that reversible mechanisms should be involved. In particular, we can mention the slippage of chains over filler particles, drawing of chains from coils or from their layers coating the surface of filler particles, slippage of macromolecular chains at topological entanglement junctions [START_REF] Svistkov | Modeling of relaxation processes in filled elastomeric materials[END_REF], etc. This latest mechanism may explain the stiffening phenomenon observed in Section 2 along contracted directions. With a view to simplification, material directions of the idealized network are assumed to be independent one on another. Stress-softening and stiffening phenomena in a direction are supposed to only depend on the alteration process in this direction, involving both breakage and slippage mechanisms.

The alteration process is modelled by considering an alteration function of both a stress-softening and a stiffening variable, respectively denoted α c and β c , which affect the chain oriented in the direction d c of the discrete network.

The alteration function is identical for all directions (isotropy of the mechanism) and describes the evolution of the average number of monomers per chain upon loading. The material is assumed to have a directional selective memory of the alteration process such that for c = 1, ..., r:

N c (t) = N 0 f αc , βc , α c , β c , (3)
where N 0 = N c (t = 0) denotes the average number of monomers per chain within the virgin material and where alteration variables α c (t) and β c (t) are defined along the considered chain direction c by:

α c (t) = max 0≤τ ≤t [λ c (τ )] (4)
and

β c (t) = min 0≤τ ≤t [λ c (τ )] (5)
ẋ denotes the rate of x.

The present model differs from those proposed by Diani et al. (2006a) or [START_REF] Göktepe | A micro-macro approach to rubber-like materials. Part III: The micro-sphere model of anisotropic Mullins-type damage[END_REF] because the alteration function does not only depend on the maximum stretch. To simplify, we assume that the evolution of the number of monomers by chain only depends on the history of the each direction. We also assume that the material presents a continuous and discontinuous evolution law as proposed by [START_REF] Miehe | Discontinuous and continuous damage evolution in Ogden type large strain elastic materials[END_REF]. A term of continuous damage D c is introduced. D c is supposed to increase whatever the chain increases or decreases, and only depends on the stretch rate | λc |. A term of discontinuous evolution is added to take into account the dependence upon the maximum and the minimum stretches (respectively α c and β c ) and the associated rates (respectively αc and βc ).

(N c C c R ) = N 0 C 0 R (1 -D c ) g αc , βc , α c , β c (6)
with D c defined by Ḋc = D(| λc |).

Thermodynamical aspects

Within this framework, the amended macroscopic strain energy density is defined by

W F, α, β, α, β, D ≈ r c=1 ω c (1 -D c )n c w (λ c , N c , C c r ) -pI (7)
The alteration variables α = (α c ) c=1,...,r , β = (β c ) c=1,...,r and D = (D c ) c=1,...,r may be considered as internal variables. Then the evolution of N c must be consistent with the Clausius-Planck inequality that guarantees a positive macroscopic dissipation. The dissipation D int is defined for isothermal deformation by:

D int = Π : Ḟ -Ẇ 0 ∀ F with det F = 1, (8)
According to the model of [START_REF] Miehe | A micro-macro approach to rubberlike materials -Part I: the non-affine micro-sphere model of rubber elasticity[END_REF] and Eqs. ( 6), the strain energy density W can be writen as:

W = C 0 R N 0 r c=1 ω c (1 -D c )g αc , βc , α c , β c λ c √ N c ξ c + ln ξ c sinh ξ c +W 0 (9)
with

ξ c = L -1 (λ c √ N c ). We note χ(ξ) = λ c √ N c ξ c + ln ξ c
sinh ξ c . N c depends on α c and β c while (N c C c R ) depends on α, β, α and β, and finally Ḋ depends on | λc |. Therefore the expression of the dissipation D int is obtained by developping partial derivatives of each term.

D c int = + ∂g ∂ αc ∂ αc ∂t -ω c C 0 R N 0 (1 -D c )χ(ξ) ∂g ∂α c αc + ∂g ∂β c βc + ∂g ∂ βc ∂ βc ∂t -ω c C 0 R N 0 (1 -D c )gN 0 ∂χ(ξ) ∂N c ∂f ∂α c αc + ∂f ∂ αc ∂ αc ∂t + ∂f ∂β c βc + ∂f ∂ βc ∂ βc ∂t + ω c C 0 R N 0 gχ(ξ) Ḋc = D c(g) int + D c(f ) int + D c(D) int (10)
D int = r c=1 (D c(g) int + D c(f ) int + D c(D) int ) = D (g) int + D (f ) int + D (D) int (11)
It is to note that the damage always increases, as the bounds of the stretches, thus Ḋc 0, D c 1, αc 0 and βc 0. Nevertheless the stiffening phenomenon imposes ∂N c ∂β c 0. So, the only way to satisfy the positivity of D int is that positive terms in Eq. ( 10) balance negative ones.

For a given load (for instance a general biaxial extension), stretches λ i∈1,2,3 are linked by the kinematic constraint of incompressibility (λ 1 λ 2 λ 3 = 1) unlike (α c ) c=1,...,r and (β c ) c=1,...,r . Reported on the Micro-sphere representation, α and β form yield surfaces whose evolution are directional and history dependant (see Fig. 11). Unfortunately, these surfaces are not regular and not convex in the undeformed configuration. We can only notice that the lower value of β is generally reached in a direction orthogonal to a direction where α is the greater. The model induces a complex expression of D int which does not allow a simple relation between terms to enforce its positiveness. We propose to verify it a posteriori in section 4.2.

Comparison of the model with experimental data

Proposed model and identification

The aim of this section is to determine the alteration function involved in the model and to estimate the material parameters in order to reproduce the trend of experimental data. The ability of the formulation to predict stress-softening, stiffening and residual strain phenomena is assessed.

Since evolution of the average chain contour length cannot be experimentally observed, the evolution of alteration during stretching is described by the following empirical alteration functions:

f αc , βc , α c , β c = exp(a 1 (α c -1)) + b 1 (β c -1), (12)
g αc , βc , α c , β c = a 3 exp(a 3 (α c -1)) αc + b 3 exp(b 4 (α c -1)) βc ( 13
)
where a 1 , b 1 , a 3 , b 3 , and b 4 are material parameters. These equations must satisfy conditions Eq. ( 10). Note that, insofar as it clearly depends on the nature of the material and the forming process, this form is specific to the material used in the present study. The current average number of monomers per chain N c (t) then follows Eq. ( 3). The larger is the value of N c , the greater is the chain locking stretch value.

The damage function is chosen as:

Ḋ = ν| λ| ( 14
)
where ν is a material parameter. That means that less and less macromolecular chains are concerned by the mechanical response of the material.

Nevertheless, by assuming ν = 0, the density of active macromolecules stay constant. The material parameters were fitted on all the experimental tests, the values are summarized in Table 3.

An illustration of the abilities of the model is presented in First, the uniaxial tension test with secondary specimens is numerically reproduced by firstly simulating kinematics of the preconditioning tensile test in order to evaluate the new stress-free configuration (average chain length and density, prior maximum and minimum chain stretches) in each discretized direction. While considering these characteristics as initial ones, analytical simulation of the subsequent tensile tests on cut samples is then achieved for the retained set of loading directions (oriented at 0 • , 22.5 • , 45 • , 67.5 • and 90 • from the pre-softening loading direction). As shown in Fig. 13, the proposed model provides a good agreement with experimental data depicted in Fig. 4, more precisely, it successfully predicts the anisotropic distribution of stress-softening, stiffening and residual strain phenomena under uniaxial tensile stretching. to t 1 , the micro-sphere is stretched in direction 90 • to λ = 3.25. Then from t 1 to t 2 the microsphere is unstained, then stretched in direction 0 • to λ = 4.5 (reaching value λ = 3.25 at time t 3 ). Finally, the microsphere is unstrained.

Thermodynamics

As seen in section 3.2, the only way to verify the positiveness of the energy dissipation is to verify it a posteriori. For the proposed equation, the derivatives to introduce in eq. 10 are:

∂f ∂α c = a 1 exp(a 1 (α c -1)), ∂g ∂α c = a 3 exp(a 3 (α c -1)) αc + b 3 b 4 exp(b 4 (α c -1)) βc ∂f ∂ αc = 0, ∂g ∂ αc = a 3 exp(a 3 (α c -1)) ∂f ∂β c = b 1 , ∂g ∂β c = 0 ∂f ∂ βc = 0, ∂g ∂ βc = b 3 exp(b 4 (α c -1)) (15)
We can notice that, for the present model, the positiveness of the energy dissipation can not be ensured mathematically. This invalidates it for further predictive simulations. But, in this paper, the model is only used to validate the assumption of a stiffening phenomenon involved in the stretching of rubber. In the following, we only try to validate the positiveness of the dissipation of the present simulation, with the material parameters given in table 3.

We first verified the positiveness of the dissipation for the above experiments. This process is repeated for many different load paths and the second principle of thermodynamics is systematically verified. Besides the generally observed stress-softening and residual strain phenomena, a stiffening effect has been stated in contracted regions beyond the maximum stretch level previously undergone by the material. Within this framework, we proposed a new network mechanism for the anisotropic description of these inelastic effects, which are considered to result from the breakdown and slippage of junctions bridging polymer chains and/or filler particles. This alteration scenario leads to an increase of the average chain contour length in tensile loaded regions and a decrease in contracted zones. This approach is implemented in a hyperelastic network-based formulation by describing the evolution of the material parameters through an alteration function involving both a softening and a stiffening variable. Since the structural rearrangements are assumed to generate a new equilibrium stress-free configuration, this evolution is described in an incremental way. Discretized integration formula, well-adapted to the representation of initial or induced anisotropy, is retained to numerically perform the homogenisation procedure.

It was shown that, our model reproduces the residual stretch as a natural consequence of the anisotropic character of the alteration process. Suitability of the resulting formulation is demonstrated on the basis of available experimental results highlighting the anisotropic material properties generated by prior loading. Finally, the proposed formulation permits to fit globally a large set of experimental data that are very different in term of strain states with few material parameters.

Fig. 1

 1 Fig. 1(b) shows the evolution of the nominal stress (first Piola-Kirchhoff stress) π as a function of the principal stretch λ along the loading direction.

Figure 1 :

 1 Figure 1: Anisotropic distribution of stress-softening developed under tensile stretching. (a): geometry of the large and the secondary small dumbbells specimens. (b): initial anisotropy of the large virgin dumbbell samples extracted from a calendered plate. Stressstretch responses of small test pieces cut at various angles from the reference direction are presented.

Fig

  Fig. 2).

Figure 2 :

 2 Figure 2: Cinematics of the deformation

Figure 3 :

 3 Figure 3: Reference basis of the small specimen extracted from the large pre-softened specimen.

Figure 4 :

 4 Figure 4: Anisotropy induced by the Mullins effect under cyclic tensile conditions up to a stretch of 3.25. Comparison of tensile stress-stretch uniaxial responses of a virgin sample and pre-softened samples cut at different angles from the reference direction. Error bars shows experimental dispersal on four samples.

  specimen is 300 mm long, 40 mm width and 2 mm thick. The width inside fixation rolls is equal to 35 mm. It is submitted to five cycles of pre-softening pure shear loading up to a stretch of 4 (Fig.5: the width of the specimen is stretched of 400%).

Figure 5 :

 5 Figure 5: The pure shear sample and the secondary small dumbbells samples used for the characterisation of the anisotropic distribution of stress-softening developed under pure shear stretching

Figure 7 :Figure 8 :

 78 Figure 6: (a) Anisotropy induced by the Mullins effect under cyclic pure shear conditions up to λ = 4. Comparison of tensile stress-stretch uniaxial responses of a virgin sample and samples cut at different angles from the reference direction, (b) close-up to residual stretches

Figure 9 :

 9 Figure 9: Anisotropy induced by the Mullins effect under cyclic modified biaxial conditions up to a stretch of 2.25 in the first direction and 3.25 in the second one Comparison of tensile stress-stretch uniaxial responses of a virgin sample and samples cut at different angles from the reference direction.

  2) where d c represents the direction vector of the chain c, λ(d c ) is the stretch of the chain c, ω c is the weight factor associated to the chain c, C R = nkT represents the material stiffness expressed in terms of the chains density n, I is the unit tensor and p a pressure classically introduced to impose incompressibility. This hyperelastic constitutive equation contains two physicallymotivated material parameters N and C R .

From a statistical point

  of view, chains embedded in the macromolecular network are improbably oriented along a straight line. Considering a chain passing through an entanglement loop formed by another molecule, as shown in Fig. 10(a), ends of this chain are certainly oriented in different directions. Application of an uniaxial tensile external force induces the slippage of the chain through this entanglement junction. As a chain is a sequence of monomer segments between two subsequent junctions, the length and consequently the number of monomers of the chain oriented along the tensile direction increases (i.e. the part of the chain situated above the junction point in Fig. 10(b)), whilst the one of the other chain (i.e. the part of the chain situated below the junction point) decreases. We can presume that interactions with the neighbouring molecular network or eventually the quasi-irreversible tightening of the entanglement loop make this change of conformation not fully recoverable during unloading (Fig. 10(c)). Another consequence is that the number of monomers in each direction changes and the mass conservation cannot be considered in each direction, but globally (if any).

Figure 10 :

 10 Figure 10: Schematic representation of a chain, depicted in bold line, slipping through an entanglement loop formed by another macromolecule (in grey line): (a) initial undeformed configuration, (b) deformed configuration under uniaxial tensile loading and (c) stress-free configuration after application of the stretching.

Figure 11 :

 11 Figure 11: Representation of the internal state variables in the Micro-sphere. Vectors d c are scaled respectively with 1 (white), with α c (red) and with β c (green) for a biaxial extension (a) and an uniaxial extension (b). The three surfaces are concentric; offsets on axis 3 have been applied for illustration.

Figure 12 :

 12 Figure 12: Example of a cyclic loading with the "extension" component and the "path" component of the stess softening model.

Figure 13 :

 13 Figure 13: Model prediction (right) of the anisotropy induced by the Mullins effect under cyclic tensile conditions up to a stretch of 3.25, compared with experiments (left).

Fig. 14

 14 Fig. 14 presents the evolution of variables involved in the proposed formulation and highlights the influence of history of each chain in the global response of the specimen. Loading is decomposed in time ranges. From t = 0

Figure 14 :

 14 Figure 14: Evolution of internal variables during unixial stretch for chains oriented at 0 • and 90 • respectively.

  When the both chains reach the same value of α c and β c (i.e. α 90 = α 0 and β 90 = β 0 ) the parameters N c are identical (N 90 = N 0 = 12.21) but the damage parameters are different (D 0 = 0.029 and D 90 = 0.018) as they depends on the loading paths. That explains the difference of the responses observed in Fig. 13. Pure shear tests with secondary specimens are then simulated in a similar way. As illustrated in Fig. 15, model predictions correctly describe the dependence of the material response on the direction, this dependence being induced by the prior pure shear loading. Nevertheless, the remaining softening effect, significant on experimental results above the prior maximum stretch level endured, is not reproduced by the constitutive equation. predictions (right) of the anisotropy induced by the Mullins effect under cyclic pure shear conditions up to a stretch of 4, compared with experiments (left).Results of calculations for the modified equibiaxial test with secondary samples are shown in Fig.16. The response of the model is corrected by the initial anisotropy in order to be compared with experimental data of Fig.8. Despite a small discrepancy about the amount of softening, the general trend of data is fairly well described below and beyond the pre-softening stretch level. While comparing experimental and theoretical results for the modified biaxial respectively shown in Fig.9and 17, a similar conclusion can be drawn. Model predictions (right) of the anisotropy induced by the Mullins effect under cyclic modified "equibiaxial" conditions up to a stretch of 3.25 in both directions, compared with experiments (left).

  Model predictions (right) of the anisotropy induced by the Mullins effect under cyclic modified "biaxial" conditions up to a stretch of 2.25 in the first direction and 3.25 in the second one, compared with experiments (left).

Figure 18 :

 18 Figure 18: Dissipation inside the material during the uniaxial tensile test. A first extension to λ = 3.25 is applied in direction 0 • then a second extension of λ = 6 in a different direction.

Fig. 19

 19 Fig. 19(b) shows that of negative part D g of the dissipation due to function g of equation (9) is counterbalanced by the other terms of the dissipation. It is also to note the part D D of the damage function is negligible with respect to the other terms.

  Figure 19: (a) Arbitrary biaxial load path given as function of λ 1 and λ 2 , (b) coresponding dissipation given as function of λ 2

  dot product of A by b A : B double-dot product of A by B A ⊗ B dyadic product of A by B ≈ is approximated by...

Table 1 :

 1 Value of corrective factors depending on the direction of extraction of the small

	Since

specimen relatively to the calendered direction of the large specimen

Table 2 :

 2 Residual strain of the small specimens in their principal axes

	λ 2res	λ 3res

• ) λ 1res

All the samples cut in a specimen which was already stretched will be called "secondary specimens".

Acknowledgment

We would like to thank Odile Duroc of the Advanced Engineering Department of Trelleborg Automotive for her enthusiasm and help to carry out the unconventional experiments reported in of this paper.

We also thanks Pr. Erwan Verron and Pr. Laurent Gornet for their help during the Master-thesis intership of M. Le Saux.


Jonathan Barés

Marina Barlet

Cindy L Rountree

Luc Barbier

Daniel Bonamy
email: daniel.bonamy@cea.fr

Mikko Alava

Eytan Katzav

Nominally brittle cracks in inhomogeneous solids: from microstructural disorder to continuum-level scale

Keywords: fracture, solid mechanics, crackling, fractals, scaling laws www.frontiersin.org

We analyze the intermittent dynamics of cracks in heterogeneous brittle materials and the roughness of the resulting fracture surfaces by investigating theoretically and numerically crack propagation in an elastic solid of spatially-distributed toughness. The crack motion splits up into discrete jumps, avalanches, displaying scale-free statistical features characterized by universal exponents. Conversely, the ranges of scales are non-universal and the mean avalanche size and duration depend on the loading microstructure and specimen parameters according to scaling laws which are uncovered. The crack surfaces are found to be logarithmically rough. Their selection by the fracture parameters is formulated in term of scaling laws on the structure functions measured on one-dimensional roughness profiles taken parallel and perpendicular to the direction of crack growth.

INTRODUCTION

Understanding how solids break continues to pose significant fundamental challenges. For brittle solids broken under tension, Linear Elastic Fracture Mechanics (LEFM) tackles the difficulty by reducing the problem to the destabilization and subsequent growth of a dominant pre-existing crack (see e.g., [START_REF] Lawn | Fracture of brittle solids (Cambridge Solide State Science)[END_REF] for an introduction to LEFM). The theory is based on the idea that, in an elastic material, all dissipative and damaging processes are localized in a small zone around the crack tip, so-called fracture process zone, FPZ. Crack destabilization and further motion are then governed by the balance between the flux of mechanical energy released in the FPZ from the surrounding material and the dissipation rate into this zone. The former is computable within linear elasticity theory and connects to the stress intensity factor, which characterizes the near-tip stress field and depends on the external loading and specimen geometry only. The dissipation rate is quantified by the fracture energy required to expose a new unit area of cracked surfaces. Equivalently, it can be quantified by the fracture toughness, i.e., the onset value for the stress intensity factor above which crack starts to grow. Within LEFM theory, both the fracture energy and the fracture toughness are assumed to be material constants, to be measured experimentally.

LEFM framework provides a coherent framework to describe fracture in homogeneous solids. In contrast, heterogeneous solids remain unclear. Stress concentration at the crack tip makes the behavior observed at the continuum-level scale extremely sensitive to material's small scale inhomogeneities. Consequences include crackling dynamics for fracture with discrete pulses (avalanches) of a variety of sizes and erratic crack paths (see [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF] for review). A generic observation in this field is the existence of scale free statistics: (i) Direct imaging of the crack motion in paper sheets [START_REF] Santucci | Subcritical statistics in rupture of fibrous materials: experiments and model[END_REF][START_REF] Stojanova | High frequency monitoring reveals aftershocks in subcritical crack growth[END_REF] or along heterogeneous interfaces [START_REF] Måløy | Local waiting time fluctuations along a randomly pinned crack front[END_REF][START_REF] Tallakstad | Local dynamics of a randomly pinned crack front during creep and forced Frontiers in Physics[END_REF] has revealed power-law distribution for the size of the crack jumps; (ii) the acoustic [START_REF] Petri | Experimental evidence for critical dynamics in microfracturing processes[END_REF][START_REF] Garcimartin | Statistical properties of fracture precursors[END_REF][START_REF] Baro | Statistical similarity between the compression of a porous material and earthquakes[END_REF] or seismic [START_REF] Bak | Unified scaling law for earthquakes[END_REF][START_REF] Corral | Long-term clustering, scaling, and universality in the temporal occurrence of earthquakes[END_REF] events going along with fracture are characterized by power-law distribution for the energy; and iii) fracture surfaces are found to exhibit self-affine morphological features [START_REF] Mandelbrot | Fractal character of fracture surfaces of metals[END_REF][START_REF] Bouchaud | Fractal dimension of fractured surfaces: a universal value?[END_REF][START_REF] Måløy | Experimental measurements of the roughness of brittle cracks[END_REF][START_REF] Ponson | 2d scaling properties of experimental fracture surfaces[END_REF][START_REF] Ponson | Failure mechanisms and surface roughness statistics of fractured fontainebleau sandstone[END_REF] or logarithmic [START_REF] Dalmas | Crack propagation through phase separated glasses: effect of the characteristic size of disorder[END_REF] roughness.

These observations, by essence, cannot be addressed by conventional LEFM continuum approaches. In this context, over the past 25 years, promising alternative approaches have emerged from statistical and non-linear physics: Statistical lattice models like fiber bundle models (see [START_REF] Pradhan | Failure processes in elastic fiber bundles[END_REF] for review) or random fuse models (see [START_REF] Alava | Statistical models of fracture[END_REF] for review), for example, were found to reproduce, with a minimal set of ingredients, crackling dynamics in failure [START_REF] Hansen | Burst avalanches in bundles of fibers -local versus global load-sharing[END_REF][START_REF] Zapperi | Crack roughness and avalanche precursors in the random fuse model[END_REF][START_REF] Danku | Temporal and spacial evolution of bursts in creep rupture[END_REF] and self-affine fracture roughness [START_REF] Hansen | Roughness of crack interfaces[END_REF][START_REF] Nukala | Crack surface roughness in threedimensional random fuse networks[END_REF][START_REF] Nukala | Fracture roughness in three-dimensional beam lattice systems[END_REF], in qualitative agreement with some of the experimental observations. However, these approaches remain too minimal to provide quantitative predictions for situations of engineering interests.

A different approach consists in considering the crack propagation in a solid with spatially-distributed toughness. The availability of asymptotic formulas [START_REF] Rice | 1st-order variation in elastic fields due to variation in location of a planar crack front[END_REF][START_REF] Movchan | On perturbations of plane cracks[END_REF][START_REF] Leblond | Crack paths in three-dimensional elastic solids. i: two-term expansion of the stress intensity factors -application to crack path stability in hydraulic fracturing[END_REF][START_REF]Perturbation approaches of a planar crack in linear elastic fracture mechanics: a review[END_REF] for the variations of stress intensity factors along a slightly distorted crack can explicitly take into account the microstructural disorder in a continuum-level scale LEFM-like description [START_REF] Gao | A first order perturbation analysis on crack trapping by arrays of obstacles[END_REF][START_REF] Schmittbuhl | Interfacial crack pinning: effect of nonlocal interactions[END_REF][START_REF] Larralde | The shape of slowly growing cracks[END_REF][START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF]. Within this framework, referred thereafter to as the Random-Toughness Continuum-Mechanics (RT-CM) approach, the fracture onset can be mapped to a critical depinning transition [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF][START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF][START_REF] Ponson | Depinning transition in the failure of inhomogeneous brittle materials[END_REF][START_REF] Bonamy | Failure of heterogeneous materials: a dynamic phase transition?[END_REF]. Depending on the situation, this approach yields logarithmically rough [START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF] or self-affine [START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF] fracture surfaces. It can also explain the crackling dynamics sometimes observed as a consequence of a selfadjustment of the driving force around the depinning value [START_REF] Bonamy | Crackling dynamics in material failure as the signature of a self-organized dynamic phase transition[END_REF]. The main advantage here is that the external parameters involved in the depinning model can be mapped to conventional LEFM parameters. This has been used to relate effective toughness and microstructural disorder [START_REF] Patinet | Quantitative prediction of effective toughness at random heterogeneous interfaces[END_REF][START_REF] Démery | From microstructural features to effective toughness in disordered brittle solids[END_REF], or to unravel the specific conditions required to observed crackling in brittle fracture [START_REF] Barés | Crackling versus continuumlike dynamics in brittle failure[END_REF].

In the present work, depinning models unravel how the loading rate, specimen geometry, and microstructural disorder quantitatively select the statistics of continuum-level scale avalanches (as measured in conventional experimental fracture mechanics) and of the crack roughness (as recorded in conventional fractography). Section 2 presents the derivation of the depinning model from LEFM. This section details the successive assumptions and their implications. Section 3 assesses the statistics of the continuum-level scale avalanches, i.e., bursts evidenced in the time evolution of the spatially-averaged crack velocity. Statistics on the avalanche's size and duration demonstrate a power-law characterized by universal exponents. Conversely, the associated cutoffs do depend on the loading, microstructure, and specimen parameters according to scaling laws which are uncovered. Two regimes can be distinguished: A regime of pseudo-isolated avalanches not so different from what is predicted in the quasi-static limit (vanishing loading rate); and a regime where avalanches coalesce with each other. Section 4 analyzes the statistics of fracture roughness. The structure functions measured along profiles parallel and perpendicular to the direction of crack growth exhibit logarithmic scaling with prefactors and characteristic length-scales depending on the Poisson's ratio and microstructure parameters according to scaling which are uncovered.

MATERIALS AND METHODS

In LEFM theory, the crack velocity v is governed by the balance between the mechanical energy G and the (material constant) fracture energy (Griffith criterion). For slow cracks, v/μ = Gwhere μ = c R / relates the mobility μ to the Rayleigh wave speed c R . The principle of local symmetry (PLS) governs the growth direction [START_REF] Gol'dstein | Brittle fracture of solids with arbitrary cracks[END_REF]. This imposes a zero stress intensity factor condition for mode II (K II = 0) all along the fracture. As consequences, an initially straight crack in an ideal homogeneous material gently loaded in mode I would continuously grow within its plane, without any jerky motion or roughness. But inhomogeneities at the microstructure scale yield distortions of the crack front, which, in turn, induce perturbations δK I , δK II , and δK III in the local loading of the front.

In the RT-CM approach, the effect of material inhomogeneities is taken into account by introducing a random spatiallydistributed component for the fracture energy and for the shear part of the loading. Then, Griffith criterion and PLS combined with an asymptotic estimation of the perturbations in the loading induced by the front distortion describes the crack growth. This approach was pioneered by Gao and Rice [START_REF] Gao | A first order perturbation analysis on crack trapping by arrays of obstacles[END_REF] and subsequently developed [30-32, 37, 40, 42-44]. Here, the model is reviewed with an emphasis on the relation between the model parameters and experimentally measurable quantities. We also pay special attention to list the various assumptions and discuss the implied limits.

FROM FRONT DISTORTIONS TO THE PERTURBATIONS IN THE CRACK-TIP LOADING

Let one consider a crack embedded in an isotropic elastic solid of size L × H × W under tension. In the following, we adopt the usual convention of fracture mechanics and the axes e x , e y and e z align with the direction of crack propagation (L direction), tensile loading (H direction), and mean crack front (W direction), respectively. For a straight crack, the mode I stress intensity factor K 0 fully characterizes the near-tip stress field loading. Let one now consider the situation depicted in Figure 1: left where the presence of inhomogeneities yields small in-plane and out-of-plane crack distortions, f (z, t) and h(x = f (z, t), z) respectively. The first step in the RT-CM approach is to estimate how these distortions perturb locally the loading. To make the problem tractable, simplifying hypothesis are needed:

Hyp. 1 The Young modulus E and Poisson ratio ν of the solid are considered as homogeneous. Hyp. 2 The spatial distribution of (x, y, z) and the PLS noise are narrow enough such that a first-order perturbation analysis can characterize the coupling of the local perturbations δK I (z), δK II (z) and δK III (z) to the front distortions {f , h}.

Then, Rice [START_REF] Rice | 1st-order variation in elastic fields due to variation in location of a planar crack front[END_REF] and Movchan et al's [START_REF] Movchan | On perturbations of plane cracks[END_REF] formula relate δK i (z) to the front distortions:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ δK I (z, {f }, {h}) = K 0 2π front f (ζ,t) -f (z,t) (ζ -z) 2 dζ δK II (z, {f }, {h}) = K 0 2 ∂h ∂x x = f (z, t), z -K 0 2 -3ν 2π(2 -ν) front h(x = f (z,t),ζ )-h(x = f (z,t),z) (ζ -z) 2 dζ δK III (z, {f }, {h}) = K 0 2(1 -ν 2 ) 2 -ν ∂h ∂x x = f (z, t), z (1)
where denotes a principal value integral. Note that δK I depends only on the in-plane distortions {f }, and δK II , δK III depend on the out-of-plane distortions {h}. Here, several additional assumptions were required:

Hyp. 3 In the expressions of δK II and δK III , we have omitted the terms brought by the non-singular stresses near the reference straight crack front (the so-called T-stress and the third term A of the asymptotic expansion of neartip stress field). This assumption remains valid for stable crack paths (T ≤ 0 [START_REF] Cotterell | Slighly curved or kinked cracks[END_REF]) when the roughness wavelength is small with respect to the system size and the characteristic distances defined by the external loading [START_REF] Larralde | The shape of slowly growing cracks[END_REF][START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF].

A complete form of the omitted terms can be found in Movchan et al. [START_REF] Movchan | On perturbations of plane cracks[END_REF]; Hyp. 4 Rice's formula for δK I assumes a sample with infinite width W and thickness H. It remains valid as long as one considers front corrugations the wavelength of which are small with respect to H and W. A more accurate formula proposed in Legrand et al. [START_REF] Legrand | Coplanar perturbation of a crack lying on the mid-plane of a plate[END_REF] addresses the case of finite H. To the very best of our knowledge, there does not exist any formula taking into account the effect of finite width, W.

GRIFFITH CRITERION AND EQUATION OF MOTION

The above expression for the local perturbations δK i leads to a first order expression of the energy release rate G(z) as a function of the front distortions:

G(z) = G 0 1 + δK I (z) K 0 , ( 2
)
where G 0 (t, f ) = K 0 (t, f ) 2 /E is the reference energy release rate which, as K 0 (t, f ), would result from the same loading with a straight crack front within the mean plane y = 0 at the mean position x = f (t). Henceforth, the operator a indicates averaging of the variable a(z) over the spatial coordinate z. The two quantities G 0 and K 0 coincide with the ones that would have been defined in the conventional LEFM approach when ignoring front distortions. They only depend on the specimen geometry and imposed load (of which both evolve with t) and can be determined using continuum mechanics (e.g., finite element analysis). Slow cracking (as considered herein) implies that the solid is loaded by imposing external displacements rather than external stresses (the opposite would yield dynamic fracture with crack at speed on the order of the sound speed). This induces additional constraints: G 0 should decrease with f (t) (specimen compliance always increases with crack length) and increase with t (to drive a crack, external displacement can only increase with time). Without loss of generality, we set t = 0 at a time when the crack has just stopped and x = 0 is the crack tip's position at this time. G 0 (t = 0, f = 0) = where = (x, y, z) is the mean value of fracture energy. Considering the subsequent variation f (z, t) to be small with respect to the crack length at t = 0, we can write:

G 0 (t, f ) = + Ġt -G f (t), (3)
where Ġ = ∂G 0 /∂t (driving rate) and G = -∂G 0 /∂f (unloading factor) are positive constants characterizing the loading rate and the specimen geometry, respectively. Inserting Equations 1 and 3 in Equation 2, and the resulting expression for G(z) into Griffith criterion, the equation of motion for the in-plane front displacement is:

1 μ ∂f ∂t (z, t) = Ġt -G f + π front f (ζ, t) -f (z, t) (ζ -z) 2 dζ +γ x = f (z, t), y = h(x = f (z, t), z), z , ( 4
)
where γ (x, y, z) = (x, y, z)is the fluctuating part of the fracture energy. The solution of this equation provides the spacetime dynamics of the in-plane projection of the crack front. Subsequently, it gives the time evolution of the fracture velocity at the continuum-level scale, v(t) = df /dt. This is the relevant observable characteristics of the crack dynamics in standard experiments of fracture mechanics.

PRINCIPLE OF LOCAL SYMMETRY AND EQUATION OF PATH

We now derive an equation of path by making use of the principle of local symmetry [START_REF] Gol'dstein | Brittle fracture of solids with arbitrary cracks[END_REF]. To do so, we take inspiration from the work of Katzav et al. [START_REF] Katzav | Stability and roughness of tensile cracks in disordered materials[END_REF][START_REF] Katzav | Fracture surfaces of heterogeneous materials: a 2d solvable model[END_REF] to model crack path in a model 2D situation and extend it for three-dimensional solids. The idea is to decompose the front propagation along x-axis into infinitesimal straight kinks of length , identified with the microstructural length-scale characterizing the spatial distribution of (x, y, z) (see Figure 1:right). Consider now the kink occurring at the front location z between time t and t + dt. Leblond et al's formulas [START_REF] Leblond | Crack paths in three-dimensional elastic solids. ii: three-term expansion of the stress intensity factors -applications and perspectives[END_REF] relate the stress intensity factors after kinking K i (z) (i = I, II, III) to the stress intensity factors before kinking K i (z) and to the kink

angle δθ = θ (z, t + dt) -θ(z, t) = ∂ 2 h/∂x 2 (x = f (z, t), z) via: K i (z) = j = I,II,III F i,j (δθ)K j (z), (5)
where F i,j (δθ) are universal functions that were computed in Leblond [START_REF] Leblond | Crack paths in three-dimensional elastic solids. i: two-term expansion of the stress intensity factors -application to crack path stability in hydraulic fracturing[END_REF] for three-dimensional solids. The principle of local symmetry implies that K II must be uniformly zero, which leads to:

F II,I (δθ)K I (z) + F II,II (δθ)K II (z) = 0. (6)
Now, to first order in δθ, F II,I (δθ) and F II,II (δθ) are given by Leblond [START_REF] Leblond | Crack paths in three-dimensional elastic solids. i: two-term expansion of the stress intensity factors -application to crack path stability in hydraulic fracturing[END_REF]:

F II,I (δθ) = δθ + O(δθ 2 ) and F II,II (δθ) = 1 + O(δθ 2 )
. The kink angle is then related to K 0 , δK I (z) and δK II (z) via: δθ ≈ -δK II (z)/(K 0 + δK I (z)). Finally, by keeping only the first order loading perturbations , it writes: To introduce the effect of microstructural disorder, a spatially distributed noise term K 0 ϑ(x, y, z) is added to the right-handed side of the equation. In the situation of inter-granular crack growth in a material made of sintered grains (e.g. sandstone), for example, this additional term would translate the difference between the kink angle predicted by the principle of local symmetry and that truly selected so that the crack propagates between the cemented grains. Finally, combining the resulting equation with the asymptotic formula for stress intensity factors (Equation 1), one gets:

∂ 2 h ∂x 2 (x = f (z, t), z) = - δK II (z) K 0 . ( 7
2 ∂ 2 h ∂x 2 (x, z) = - ∂h ∂x (x, z) + 1 π 2 -3ν 2 -ν front h(x, ζ ) -h(x, z) (ζ -z) 2 dζ + ϑ(x, y = h(x, z), z) ( 8 )
The solution of this equation provides the topography of the postmortem fracture surfaces, h(x, z). This is the quantity of interest in fractography science. Note that this equation differs from that given in Larralde [START_REF] Larralde | The shape of slowly growing cracks[END_REF] and Bonamy et al. [START_REF] Bonamy | Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics[END_REF] as it includes an additional curvature term ∂ 2 h/∂x 2 .

RELEVANT PARAMETERS AND NUMERICAL ASPECTS

Equations 4 and 8 predict deterministically the fracture dynamics and the morphology of fracture surfaces using the following inputs:

Input 1 The loading rate Ġ; Input 2 The specimen geometry (unloading factor G and specimen width W) Input 3 The material constants (fracture energy , Poisson ratio ν, and mobility μ), Input 4 The microstructure disorder (microstructure lengthscale and the spatial distribution of the two noise terms γ (x, y, z) and ϑ(x, y, z)).

Statistically, the two noise terms are characterized by the probability density functions p(γ ) and p(ϑ) and the spatial correlation functions C γ ( r) = γ ( r 0 + r)γ ( r 0 ) r 0 and C ϑ ( r) = γ ( r 0 + r)γ ( r 0 ) r 0 . Additional assumptions simplify the problem:

Hyp. [START_REF] Måløy | Local waiting time fluctuations along a randomly pinned crack front[END_REF] The two noise terms are not correlated; At first glance, such an assumption may appear odd since both terms originate from the material heterogeneities. But due to the tensorial nature of elasticity, these heterogeneities will affect the equation of motion and the equation of path independently. As an illustrative example, let us consider again the situation of intergranular crack growth in a solid made of sintered grains. Two distinct space-dependent noise terms are required to describe the microstructural texture: The first one quantifies the local variations of adhesion between grains and mainly affects the equation of motion; the second one describes the local variations of joint orientation and affects the equation of path. Hyp. 6 Both probability density functions are Gaussian of zero mean and standard deviation γ and θ, respectively;

Hyp. 7 In both cases, the correlation function C decreases linearly with |r| over the microstructure distance and is zero above.

Note that, in elastic interfaces equations as that proposed in Equations 4 and 8, the scaling properties remain unaffected upon changes in the probability functions (p(γ ) and p(ϑ)) and in the correlation functions (C γ (|r|) and C ϑ (|r|)) [START_REF] Vandembroucq | Universal depinning force fluctuations of an elastic line: application to finite temperature behavior[END_REF]. Thanks to this series of assumptions, the two spatially distributed noise terms are fully characterized with 3 parameters: γ , θ and . At this point, the Equations 4 and 8 are coupled via the two 3D noise terms γ (x, y, z) and ϑ(x, y, z), since they both depend on the in-plane and out-of-plane positions of the crack front. A last assumption permits the decoupling of the two equations:

Hyp. 8 The 3D spatially distributed terms γ (x, y, z) and ϑ(x, y, z) reduce to their in-plane projection γ (x, z) and ϑ(x, z).

As shown in Ramanathan et al. [START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF], Equation 8with a 2D projection for the noise term yields logarithmically smooth fracture surfaces. Hence, zooming out on the fracture surfaces makes them appear flatter and flatter. When the zooming out is sufficient, fluctuations in h become negligible when compared to in-plane length scales. The two noise terms reduce to γ (x, y = h(x, z), z) ≈ γ (x, y = 0, z) and ϑ(x, y = h(x, z), z) ≈ ϑ(x, y = 0, z). This validates Hyp. 8. At this point, the selection of the fracture behavior brings into play 8 parameters: μ, , Ġ, G , γ , θ, and W. The introduction of dimensionless time t → t/ /(μ ) and length {x, y, z, f , h} → {x/ , y/ , z/ , f / , h/ }) reduces this number. Equations 4 and 8 become:

∂f ∂t (z, t) = ct -kf (t) + 1 π × front f (ζ, t) -f (z, t) (ζ -z) 2 dζ + η γ (x = f (z, t), z), (9a) 2 ∂ 2 h ∂x 2 (x, z) = - ∂h ∂x (x, z) + A π front h(x, ζ ) -h(x, z) (ζ -z) 2 dζ + ϑ(x, z). ( 9b
)
where the dimensionless noise term η is characterized by a standard deviation η = γ / and a spatial correlation length equal to unity.

The final two decoupled dimensionless equations require only 6 parameters:

• The dimensionless driving rate c = Ġ /μ 2 ,
• the dimensionless unloading factor k = G / ,

• the Poisson-ratio dependent parameter A = (2 -3ν)/(2 -ν),

• the parameters η and θ characterizing the disorder strength,

• the continuum-level vs. microstructure scale ratio N = W/ .

The following sections subsequently address the continuum-level scale statistics of crack dynamics v(t) and the morphology of the post-mortem fracture surfaces h(x, z). They also reveal their dependencies on the external parameters {c, k, A, η, θ, N}. The numerical scheme is the following. For each set of parameters, we built two uncorrelated random N × pN maps η(x, z) and ϑ(x, z) with Gaussian distribution of variance η and θ, respectively. At time t = 0, the crack front is a straight line at x = 0: f (z, t = 0) = 0 and h(x = 0, z) = 0 ∀z ∈ [0, N -1]. Solving Equation 9a gives the space-time evolution of f (z, t). Solving Equation 9b provides h(x, z). Both cases invoke a fourth order Runge-Kutta scheme and periodic boundary conditions along z: f (z = 0, t) = f (z = N, t) and h(x, z = 0) = h(x, z = N) for all t and x. This speeds up the computation time as the two integral terms in the equations can be solved in the Fourier domain.

RESULTS

CRACKLING DYNAMICS AND AVALANCHE STATISTICS

We first look at the crack dynamics. Equation 9a describes the motion of a long-range elastic chain driven in a frozen random pinning potential with a driving force F = ctkf (t) selfadjusting around the depinning value [START_REF] Bonamy | Crackling dynamics in material failure as the signature of a self-organized dynamic phase transition[END_REF]. For {c, k} → 0, the chain propagates while remaining at the critical depinning point and the crack moves through irregular jumps, or avalanches. The size S and duration D of the avalanche 1 follow a universal powerlaw distribution [START_REF] Bonamy | Crackling dynamics in material failure as the signature of a self-organized dynamic phase transition[END_REF][START_REF] Ertas | Critical dynamics of contact line depinning[END_REF]: P(S) ∝ S -τ and P(D) ∝ D -α where the exponents τ and α can be estimated using renormalization [START_REF] Ertas | Critical dynamics of contact line depinning[END_REF][START_REF] Chauve | Renormalization of pinned elastic systems: how does it work beyond one loop?[END_REF][START_REF] Dobrinevski | Statistics of avalanches with relaxation, and barkhausen noise[END_REF] or numerical [START_REF] Rosso | Roughness at the depinning threshold for a longrange elastic string[END_REF][START_REF] Duemmer | Depinning exponents of the driven long-range elastic string[END_REF] methods. The most precise estimations yields τ = 1.280 ± 0.01 and α = 1.500 ± 0.010 [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF]. These scale free statistics extend to upper cutoffs S 0 and D 0 set by the system 1 In problems of depinning of interfaces, the definition of the avalanche size differ depending on the publication. Sometimes, the size is defined as the area A swept by the front between two successive depinning configurations. Here, the size S is defined as the integral of the continuum-level scale velocity v(t) between the two successive depinning configurations. The two definitions are related by S = A/N when N is the system size.

size N: S 0 ∝ N ζ and D 0 ∝ N κ where ζ and κ are the roughness exponent and dynamic exponent, respectively. The most precise estimations of their value yield [START_REF] Rosso | Roughness at the depinning threshold for a longrange elastic string[END_REF][START_REF] Duemmer | Depinning exponents of the driven long-range elastic string[END_REF]: ζ = 0.385 ± 005 and κ = 0.770 ± 0.005.

Here, both c and k are finite. The effects of a finite unloading factor k keeping c → 0 is now fairly well documented [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF][START_REF] Laurson | Avalanches and clusters in planar crack front propagation[END_REF]: It modifies the upper cutoffs S 0 and D 0 to the preceding scale-free distributions:

P(S|N, k, c → 0) = S -τ f (S/S 0 ) with S 0 ∝ N 1/σ -1 /k 1/σ (10a) P(D|N, k, c → 0) = D -α g(D/D 0 ) with D 0 ∝ (N/k) 1/ (10b)
Here, the form of the functions f (u) and g(u) is expected to be universal. The exponents 1/σ and 1/ are also predicted to be universal: 1/σ = 0.69 ± 0.010 and 1/ = 0.385 ± 0.010 [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF][START_REF] Laurson | Avalanches and clusters in planar crack front propagation[END_REF].

Conversely, the effect of a finite c is not uncovered yet. By yielding some overlap between the avalanches [START_REF] White | Driving rate effects on crackling noise[END_REF], it can significantly alter the dynamics. In particular, a recent work [START_REF] Barés | Crackling versus continuumlike dynamics in brittle failure[END_REF] has evidenced a transition line between a crackling-like dynamics made of irregular power-law distributed jumps and a continuum-like dynamics ruled by the conventional LEFM theory. between the pulse, but becomes equal to a small value proportional to c (prefactor dependent of the Runge Kunta scheme).

The avalanches are then identified with the bursts where v(t) is above a prescribed reference level v c = 10 -3 . Their duration D is defined as the interval between the two successive intersections of v(t) with v c , and their size S is defined as the integral of v(t)v c between the same points. to depend on the reference level v c and, thus, are to be associated with the procedure to extract the avalanches. We will consider in the following only the part of the distributions above these lower cutoffs.

The power-law distributions observed in Figure 3 also exhibit upper cutoffs S 0 and D 0 . These cutoffs decrease with k in all cases. The effects of c is of two types:

• For large k/small c the cutoff does not depend on c (only on k).

• At small k/large c, the cutoff increases with c and the distribution also displays a bump at large sizes and durations. 4B). At large enough k, both S and D are independent of c. This large k regime is attributed to a regime of pseudo-isolated (pi) avalanches. The distributions are then expected to take forms similar to that of Equation 10. As a result, S pi is expected [START_REF] Ledoussal | Size distributions of shocks and static avalanches from the functional renormalization group[END_REF] to take the form

S pi ≈ S τ -1 min S 2 -τ 0 ∝ k -(2 -τ )/σ , D pi ≈ D α-1 min D 2 -α 0 ∝ k -(2 -α)/
. These two scaling are compatible with the observations at large k (dash line in Figures 4A,B). As a synthesis, the mean avalanche size and duration are found to take the following form at large k:

S (c, k c) = S pi (k) ≈ 10 -2 k -(2 -τ )/σ (11a) D (c, k c) = D pi (k) ≈ 11k -(2 -α)/ (11b)
Let us try now to characterize the effects of the avalanche overlap when k becomes small or c becomes large. Previous work [START_REF] Barés | Crackling versus continuumlike dynamics in brittle failure[END_REF] evidenced a transition between the crackling dynamics studied here and a continuum-like dynamics when c becomes large enough or k small enough. This transition is believed [START_REF] Barés | Crackling versus continuumlike dynamics in brittle failure[END_REF] to coincide with the point where the avalanche overlap percolates throughout the entire system. At constant η and N, this transition was shown [START_REF] Barés | Crackling versus continuumlike dynamics in brittle failure[END_REF] to be fully driven by the ratio c/k. We hence plotted, in Figure 5, S / S pi and D / D pi (first order estimation of the number of individual avalanches having merged together to form the bursts detected from the signal v(t)), as a function of the control parameter c/k. A coarse collapse is observed. As expected, the master curves diverge at the transition value between crackling and continuum-like dynamics (materialized by the vertical dash lines in the main panels of Figures 5A,B).

Note that significant deviations to the collapse are observed. They are believed to stem from a qualitative change in the distribution shape as k decreases/c increases and the avalanche overlap increases. Actually, for small k/large c a bump develops at large scales, for both the distribution in size and in duration. Due to this change in shape, several length-scales (resp. time-scales) intervenes in the distributions of S (resp. D). Thus, the recording of the mean values alone is not sufficient to capture their evolutions with c and k. Ongoing work aims at characterizing these effects more accurately.

To make the analysis complete, we looked at the effects of the system size N and disorder strength η. Figure 6 plots the probability density function of avalanche size P(S|N, η) and avalanche duration P(D|N, η) as measured for different values of N and η at fixed values of c and k. For both size and duration, the decrease of N yields an increase of the lower-cutoffs S min and D min (main panel of Figures 6A,A ). The two dependencies are well fitted by power-laws: S min ∝ N -a SN and D min ∝ N -a DN with a SN ≈ 1.7 and a DN ≈ 0.6 (inset of Figures 6A,A ). 11and the collapse obtained in Figure 5, we anticipate the following form for mean size and duration: where the two functions F(u) and G(u) exhibit a plateau at small u, and both diverge at the same value u c . The value of the exponents τ , α, 1/σ and 1/ are well known [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF]: They can e.g., be related to the so-called roughness exponent ζ and dynamic exponent κ classically defined in the realm of critical depinning transition: In this context, it is of interest to rewrites the equation with the original variables, before the non-dimensionalization procedure:

S ≈ ηb S η N b SN k (2 -τ )/σ F cN c N k ηc η , D ≈ 1 ηb D η N b DN k (2 -α)/ G cN c N k ηc η , ( 12
)
τ = 2 -1/(1 + ζ ), α = 1 + ζ /κ, 1/σ = (1 + ζ )/
S ≈ γ b S η × b SN -(2 -τ )/σ G (2 -τ )/σ × b S η -(2 -τ )/σ × W b SN F Ġ × c η -1 × W c N μG × γ c η × c N , ( 13a
)
D ≈ b S η +(2 -α)/ × b DN -(2 -α)/ G (2 -α)/ × W b DN G Ġ × c η -1 × W c N μG × γ c η × c N , ( 13b
)
where the mean avalanche size S and duration D are expressed with real length and time units, and Ġ, G , W, , μ, and γ are recalled to be the loading rate, the unloading factor, the specimen width, the fracture energy, the mobility, the microstructure length-scale, and the contrast in local fracture energy. The value of the exponents are recalled to be τ = 1.280 ± 0.01 (predicted), α = 1.500 ± 0.01 (predicted), 1/σ = 0.69 ± 0.01 (predicted),

ROUGHNESS OF FRACTURE SURFACES

We turn now to the topography h(x, z) of the post-mortem fracture surfaces as predicted by equation 9b. Figure 8 reports typical topographies for different values of the two external parameters A and θ. When A is close to 1, the surface seems to be statistically isotropic while as A gets smaller, the surface appears more elongated in the direction of z. Conversely, the parameter θ only affects the range swept by the roughness. Note that in almost all elastic solids, Poisson ratio ν lies between 0 and 0.5, which impose a finite interval for A = (2 -3ν)/(2 -ν): 1/3 ≤ A ≤ 1. Herein, only A within this interval are considered.

To characterize quantitatively the spatial distribution of fracture roughness, we adopted the classical procedure [START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF] and computed the structure function S( r) = (h( r + r)h( r)) 2 . Here, the operator denotes averaging over all positions r = (x, z). First, we computed the structure function S z ( z) along 1D profiles taken parallel to z (mean direction of the crack front). The procedure is the following: (i) an initially straight front was first propagated over a distance equal to 10N to obtain a statistically stationary regime; (ii) 10000 subsequent profiles h(x i , z) separated by a distance x i+1x i = 1 were recorded; (iii) the structure function S i z ( z) was computed for each of these profiles; and (iv) finally, these 10000 individual structure functions were averaged to get S z ( z).

Figure 9 depicts typical examples of S z ( z) curves for different values of N, θ and A. S z goes as S z = p z log ( z/λ z ) up to an upper cutoff set by the system size N. This logarithmic scaling is anticipated to extend over the whole range of length-scales as N → ∞. Note that logarithmically rough crack surfaces were also predicted in earlier theoretical works [START_REF] Larralde | The shape of slowly growing cracks[END_REF][START_REF] Ramanathan | Quasistatic crack propagation in heterogeneous media[END_REF] analyzing crack propagation through a three-dimensional heterogeneous solid. As a plus, the present model allows relating the prefactor p z and characteristic length-scale λ z with the fracture parameters: p z is found to scale as θ2 /A (Figure 10A), while λ z is independent of both A and θ (Figure 10B). Finally, the structure function along z is:

S z ( z) = C θ2 A log ( z/λ z )
, with C = 0.32 ± 0.01 and λ z = 0.24 ± 0.03 [START_REF] Måløy | Experimental measurements of the roughness of brittle cracks[END_REF] We now look at the structure function S x ( x) along x. A direct computation of S x following the standard procedure proposed for S z was found to give a large scattering even for a constant set of parameters {N, A, θ}. Hence, the computation procedure was modified as follows: (i) an initially straight front was first propagated over a distance equal to 10N to obtain a statistically stationary regime; (ii) the x evolution of h(x, z i ) = h(x, z i )h(x) was subsequently recorded at 100 locations z i uniformly distributed along the specimen width N (h(x) denotes averaging over the specimen width N); (iii) the structure function 2 was computed for each of these profiles; (iv) these 100 individual structure functions were averaged to get the S x ( x) for a single specimen; and (v) the so-obtained structure functions were further averaged over 100 specimens. This procedure produces accurate and reproducible curves S x vs. x.

S i z ( x) = ( h(x + x, z i ) -h(x, z i ))
Figure 11 depicts typical examples of S x ( x) curves for different values of N, θ and A. The behavior resembles that observed for S z , with a logarithmic scaling S x = p x log ( x/λ x ) up to an upper cutoff set by the system size N. Note that S x saturates above the cutoff, and does not decrease down to zero as was observed for S z due to periodic boundary conditions. As for S z , the prefactor (now referred to as p x ) goes as θ (Figure 12A). The characteristic length-scale (λ x ) is independent of θ (Figure 12B In this context, it is of interest to rewrites them with the original variables, before the non-dimensionalization procedure: S

z ≈ 0.32 2 -ν 2 -3ν θ2 2 log (4.2 z ), S x ≈ 0.32 2 -ν 2 -3ν θ2 2 log (4.8 2 -3ν 2 -ν x ), (16)
where S x , S z , x and z are expressed with real length units, and ν, and θ are recalled to be the Poisson ratio, the microstructure scale, and disorder contrast.

CONCLUDING DISCUSSION

Stress enhancement at crack tips makes the macroscale failure behavior observed extremely sensitive to the presence of disorder at the microstructure scale. This translates into crackling dynamics and rough fracture surfaces, which, by essence, cannot be addressed within the conventional LEFM framework. In this paper, we have used the RT-CM approach to obtain quantitative relations between some statistical observables characteristic of these two aspects and the fracture parameters: Loading rate (time derivative of the energy release rate), specimen geometry (specimen thickness and unloading factor), conventional mechanical constants (fracture energy, Poisson ratio), and microstructural disorder (microstructure scale and disorder strength).

Over a certain range of the fracture parameters, this RT-CM approach predicts crackling dynamics [START_REF] Barés | Crackling versus continuumlike dynamics in brittle failure[END_REF]: The crack growth splits up into discrete jumps, which are power-law distributed in size and duration. The characteristic exponents associated to these power-laws are universal. Conversely, the scales covered by these scale-free features are non-universal and, in particular, the mean size and duration of the crack jumps are found to depend on the fracture parameters according to scaling laws that are uncovered. These scaling laws can be understood over a certain range of the fracture parameters, in the regime of pseudo-isolated avalanches addressable via standard functional renormalization theory [START_REF] Ertas | Critical dynamics of contact line depinning[END_REF][START_REF] Dobrinevski | Statistics of avalanches with relaxation, and barkhausen noise[END_REF][START_REF] Rosso | Roughness at the depinning threshold for a longrange elastic string[END_REF][START_REF] Duemmer | Depinning exponents of the driven long-range elastic string[END_REF]. Conversely, the effect of the avalanche overlapping is not understood. On-going work aims at analyzing the distribution of the local avalanches as detected in the space-time diagrams of the front dynamics, in order to understand the coalescence process.

Also, this RT-CM approach predicts rough fracture surfaces. The fracture roughness can be characterized by computing the structure function, which exhibits logarithmic scaling. The associated prefactor and characteristic length-scale are found to depend on the Poisson ratio, microstructure length-scale, and disorder strength according to laws that were uncovered. This may have interesting applications: It allows one to infer the microstructure parameters (the access of which could be made difficult otherwise, due to the smallness of the length-scales involved) from the analysis of postmortem fracture surfaces at larger scale. Work in progress aims at testing the scaling predicted here for the structure functions against fractography experiments achieved in oxide glasses.

Note finally that the RT-CM model studied here call upon a variety of assumptions (see Section 2). An interesting perspective would be to measure to which extend these assumptions can be released. Work in this direction is currently under progress. The model is also limited to nominally brittle fracture, with a single macroscopic crack propagating in an otherwise intact material. Promising alternative approaches have emerged from statistical and non-linear physics [START_REF] Alava | Statistical models of fracture[END_REF][START_REF] Pradhan | Failure processes in elastic fiber bundles[END_REF] and may permit to address quasi-brittle fracture, with many microcracking events interacting with each-others. A central challenge in the field would be to bridge the gap between these approaches and engineering damage mechanics.

FIGURE 1 |

 1 FIGURE 1 | Left: Sketch and notation of a crack front propagating in a heterogeneous solid. Right: Decomposition of the crack path into infinitesimal straight kinks.

  ) www.frontiersin.org November 2014 | Volume 2 | Article 70 | 3

Frontiers

  in Physics | Interdisciplinary Physics November 2014 | Volume 2 | Article 70 | 4

Figure 2

 2 FIGURE 2 | (A) Evolutions of the continuum-level scale velocity v(t) for three typical parameter sets within the crackling phase: {c = 10 -5 , k = 10 -3 , η = 1, N = 1024} (top), {c = 10 -5 , k = 10 -2 , η = 1, N = 1024} (middle), {c = 10 -4 , k = 10 -2 , η = 1, N = 1024} (bottom). (B) Zoom in on an avalanche, identified as a pulse where v(t) is above a prescribed value v c = 10 -3 (horizontal red dash line). The avalanche duration D is the interval between the two successive intersections of v(t) with v c . The avalanche size S is the integral of v(t)v c between the same points.

Figure 3

 3 reports the probability density function of avalanche size P(S|c, k) and duration P(D|c, k) for different values of c and k. Power-laws are observed. The exponents τ and α associated with the power-law decrease are independent of c and k. Moreover, they compare well with the values τ = 1.280 ± 0.01 and α = 1.500 ± 0.010 expected for {c, k} → 0. On the other hand, the valid region of the power-law is observed depends on c and k. Both the lower cutoffs S min and D min are roughly independent of c and k. However, the precise values S min and D min were observed

  Direct computation of the cutoff is quite imprecise. Hence, the selection of the typical length and time scales is studied via variations of the mean values S (c, k) = ∞ S min S × P(S|c, k)dS

FIGURE 3 | 2 . 6 (

 326 FIGURE 3 | Top: Distribution of the avalanche size measured (A) for various values of k at constant c = 10 -4 and (B) for various values of c at constant k = 10 -2 .In both cases, N = 1024 and η = 1 and the axes are logarithmic. The power-law exponent τ is found to be independent of the parameters and compatible to the universal value τ = 1.28 (black dashed line) predicted for c → 0 and k → 0. The lower cutoff is found to be independent of the parameters: S min ≈ 10 -3 (vertical dash line). The upper-cutoff is found to decrease with k and to increase with c (resp. to be independent of c) when c is large enough (resp. when c is small enough).Bottom: Distribution of the avalanche duration measured (A ) for various values of k at constant c = 10 -4 and (B ) for various values of c at constant k = 10 -2 . In both cases, N = 1024 and η = 1 and the axes are logarithmic. The power-law exponent α is found to be independent of the parameters and compatible to the universal value α = 1.50 (black dashed line) predicted for c → 0 and k → 0. The lower cutoff is found to be independent of the parameters: D min ≈ 3 (vertical dash line). The upper-cutoff is found to decrease with k and to increase with c (resp. to be independent of c) when c is large enough (resp. when c is small enough).

  η does not seems to affect S min . Conversely, D min decreases with η as D min ∝ η-a D η with a D η ≈ 1.2 (inset of Figure 6B ). The effects of N and η of S and D are analyzed in Figure 7, by plotting the curves S vs c and D vs c for different N at constant k and η (insets in Figures 7A,A ), and for different η at constant k and N (insets in Figures 7B,B ). Increasing N yields a decrease of the low-c plateau and a leftward shift of the divergence location for both S and D . Increasing η yields an increase of the low-c plateau for S , a decrease of the low-c plateau for D , and a rightward shift of the divergence location for both S and D . All curves can then be superimposed by making { S → S /N -b SN , c → c/N -c N } in the main panel of Figure 7A, { S → S / ηb S η , c → c/ ηc η } in the main panel of Figure 7B, { D → D /N -b DN , c → c/N -c N } in the main panel of Figure 7A , { D → D / η-b D η , c → c/ ηc η } in the main panel of Figure 7B . By combining this with Equation

FIGURE 4 | 7 FIGURE 5 |

 475 FIGURE 4 | (A) Evolution of the mean avalanche size S as a function of the unloading factor k for different loading rates c. The curves collapse for large k onto a c-independent power-law curve. Black thick dash line depicts a fit over the collapsed region using Equation 11a: S pi ≈ 10 -2 k -(2 -τ )/σ with τ = 1.28 and 1/σ = 0.69. (B) Evolution of the mean avalanche duration D as a function of the unloading factor k for different value of the loading rate c. At large k, the curves collapse onto a c-independent power-law curve. Black thick dash line is a fit over the collapsed region using Equation 11b : D pi ≈ 11k -(2 -α)/ with α = 1.50 and 1/ = 0.385. In both (A,B) N = 1024 and η = 1 and the axes are logarithmic. The different symbols correspond to different values of c, given in the legends.

2 and 1 /

 1 = κ/(1 + ζ ). Conversely, the precise origin of the exponents b SN , b S η, b DN , b D η, c N , and c η and their link with ζ and κ remain to be uncovered. Equation 12 quantitatively relate the material parameters to quantities that are accessible in conventional experimental mechanics, namely the mean size and duration of the avalanches.

1 /

 1 = 0.385 ± 0.01 (predicted), b SN = 1.3 ± 0.1 (fitted), b S η = 0.7 ± 0.1 (fitted), b DN = 0.45 ± 0.1 (fitted), b D η = 0.65 (fitted), c N = 0.65 ± 0.1 (fitted), and c η = 1.05 ± 0.1 (fitted).

Frontiers 8 FIGURE 6 |

 86 FIGURE 6 | Effect of the system size N and disorder strength η on the avalanche statistics. Main panel, top: Distribution of the avalanche size measured (A) for different N at constant c = 2 × 10 -5 , k = 10 -2 and η = 1, and (B) for different η at constant c = 2 × 10 -5 , k = 10 -2 and N = 1024. The axes are logarithmic. Main panel, bottom: Distribution of the avalanche duration measured (A ) for different N at constant c = 2 × 10 -5 , k = 10 -2 and η = 1, and (B ) for different η at constant c = 2 × 10 -5 , k = 10 -2 and N = 1024. All axes are logarithmic. In all panels, the power-law exponents τ and α are compatible with the universal values τ = 1.28 and α = 1.50 (inclined dashed line) predicted for c → 0 and k → 0. Regarding the size, the lower cutoff S min is observed to be independent of η and to decrease with N. Regarding the duration, the lower cutoff D min decreases with both N and η. S min (resp. D min ) is quantitatively defined as the intersection of the power-law regime with exponent τ (resp. α) and the saturation value for P observed at small S (resp. small D). The red and blue dash lines in (A,A ,B ) present illustrations for N = 128 and N = 1024, respectively. The variation of S min with N is shown in the inset of a. The variations of D min with N and η are shown in the insets of (A ,B ). The different symbols correspond to different values for c (from 2 × 10 -6 to 10 -4 ). The red line are fits: S min ∝ N -a SN with a SN = 1.7 ± 0.1, D min ∝ N -a DN with a DN = 0.6 ± 0.1, and D min ∝ η-a D η with a D η = 1.2 ± 0.1.

9 FIGURE 7 | 1 FIGURE 8 |

 9718 FIGURE 7 | Effect of the system size N and disorder strength η on the mean size and duration of avalanches. Top, inset: Mean avalanche size S as a function of c for constant k = 10 -2 and (A) different N at constant η = 1 and (B) different η at constant N = 1024. Bottom, inset: Mean avalanche duration D as a function of c and constant k = 10 -2 and (A ) different N and constant η = 1 and (B ) different η and constant N = 1024. In all cases, the axes are logarithmic. Main panels: Curve collapse obtained by plotting S /N -b SN vs c/N -c N (A), S / η-b S η vs c/N -c η (B), D /N -b DN vs c/N -c N (A ), and D / η-b D η vs c/N -c η . The fitted exponents are found to be b SN = 1.3 ± 0.1, b S η = 0.7 ± 0.1, b DN = 0.45 ± 0.1, b D η = 0.65 ± 0.1, c N = 0.65 ± 0.1, c η = 1.05 ± 0.1

  :inset). But contrary to what is observed for S z , this characteristic length λ x depends on A: λ x ∝ 1/A. This dependency is responsible for the apparent stretching along s of the images in Figure 10 observed as A decreases. Finally, the structure function along x is: S x ( x) = C θ2 A log (A z/D), with C = 0.32 ± 0.01 and D = 0.21 ± 0.02 (15) Equations 14 and 15 quantitatively relate the material parameters (microstructure and Poisson ratio) to quantities accessible in conventional fractography analysis.

FIGURE 9 |FIGURE 10 |

 910 FIGURE 9 | Structure function S z ( z) measured along z: (A) for various system sizes N at A = 1 and θ = 1, (B) for various θ at A = 1 and N = 1024, and (C) for various A at θ = 1 and N = 1024. The ordinate axis is logarithmic. For all figures, the straight lines correspond to the solutions of Equation 14.

FIGURE 11 |

 11 FIGURE 11 | Structure function S x ( x) measured along x (A) for various system sizes N at A = 1 and θ = 1, (B) for various θ at A = 1 and N = 1024, and (C) for various A at constant θ = 1 and N = 1024. The ordinate axis is logarithmic. For all figures, the straight lines correspond to the solutions of Equation 15.

FIGURE 12 |

 12 FIGURE 12 | (A)Slope p x associated with the curve S x vs. x as a function of A at constant θ = 1 (main) and θ at constant A = 1 (inset). In the inset, the axes are logarithmic. In both graphs, the red lines are fits p 1 = C x /A (main) and p 1 = C θ 2 where the fitted parameter is found to be C = 0.31 ± 0.02 (95% confident interval). (B) Characteristic length-scale λ x associated with the curve S x vs. x as a function of A at constant θ = 1 (main) and θ at constant A = 1 (inset). In both graphs, the red lines are fits λ x = 0.21 ± 0.02/A (main) and λ x = 0.21 ± 0.02 (inset). Here, ± indicates a 95% confident interval.

Frontiers inPhysics | Interdisciplinary Physics November 2014 | Volume 2 | Article 70 | 10

ACKNOWLEDGMENTS

The authors would like to thank Alberto Rosso for enlightening discussions and the "Agence Nationale de la Recherche" (ANR) for financial support of the project MEPHYSTAR (ANR-09-SYSC-006-01).

Conflict of Interest Statement:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.