Anthony Hun

Johann Guilleminot
email: johann.guilleminot@duke.edu

Julien Yvonnet

D A Hun

| J Guilleminot

| J Yvonnet

| M Bornert

Stochastic Multi-Scale Modeling of Crack Propagation in Random Heterogeneous Media

Keywords: Crack propagation, phase-field method, mesoscale, stochastic analysis, random field, uncertainty quantification

come

Stochastic multiscale modeling of crack propagation in random heterogeneous media

INTRODUCTION

Modeling the damage and strength properties in heterogeneous brittle or quasi-brittle materials is of major interest in civil and mechanical engineering, biomechanics and many other areas of engineering. In such complex materials like concrete, cementitious materials, cortical bone, or fiber/particle-reinforced composites, the main damage mechanisms are due to the initiation and propagation of crack networks at the microscale. However, in most applications, modeling the crack propagation at the scale of constituents up to the macroscale is, in general, not tractable. In the recent years, a few attempts have been devoted to multiscale modeling of damage from microscale crack propagation, including, among many others: FE 2 -approaches considering separated scales and microscale nonlinear damage [START_REF] Kulkarni | Coupled multi-scale cohesive modeling of failure in heterogeneous adhesives[END_REF][START_REF] Nguyen | Homogenization-based multiscale crack modelling: from micro-diffusive damage to macro-cracks[END_REF], extended FE 2 techniques with discontinuities at the macroscale [START_REF] Loehnert | A multiscale projection method for macro/microcrack simulations[END_REF][START_REF] Liu | Multiscale simulation of major crack/minor cracks interplay with the corrected XFEM[END_REF][START_REF] Coenen | A multi-scale approach to bridge microscale damage and macroscale failure: a nested computational homogenization-localization framework[END_REF][START_REF] Oliver | Continuum approach to computational multiscale modeling of propagating fracture[END_REF], domain decompositions methods embedding discontinuities [START_REF] Rudoy | Domain decomposition method for crack problems with nonpenetration condition[END_REF][START_REF] Teichtmeister | Phase field modeling of fracture in anisotropic brittle solids[END_REF][START_REF] Waisman | An adaptive domain decomposition preconditioner for crack propagation problems modeled by XFEM[END_REF], or more recently methodologies enabling the identification of regularized damage models at the macroscale from microscale calculations [START_REF] Hossain | Effective toughness of heterogeneous media[END_REF][START_REF] Nguyen | Identification of fracture models based on phase field for crack propagation in heterogeneous lattices in a context of non-separated scales[END_REF].

In the mentioned works, only deterministic geometries and crack paths were considered. Since microstructural randomness strongly impacts the macroscopic response of (quasi-)brittle materials in various ways, ranging from size effects [START_REF] Bažant | Fracture and Size Effect in Concrete and Other Quasibrittle Materials[END_REF][START_REF] Bažant | Probabilistic Nonlocal Theory for Quasibrittle Fracture Initiation and Size Effect. I: Theory[END_REF][START_REF] Vliet | Experimental investigation of size effect in concrete and sandstone under uniaxial tension[END_REF][START_REF] Genet | Scaling strength distributions in quasi-brittle materials from micro-to macro-scales: A computational approach to modeling Nature-inspired structural ceramics[END_REF] to high stochasticity in failure patterns and ultimate properties [START_REF] Tregger | Modeling the effect of mesoscale randomness on concrete fracture[END_REF][START_REF] Daphalapurkar | Predicting variability in the dynamic failure strength of brittle materials considering pre-existing flaws[END_REF][START_REF] Rahman | Stochastic multiscale fracture analysis of three-dimensional functionally graded composites[END_REF], the development of approaches incorporating multiscale and probabilistic ingredients all together is a natural path to extend the predictive capabilities in fracture simulations. The main objective of this work is to propose a stochastic, multiscale-informed phase-field approach to model crack propagation in heterogeneous media. In the proposed framework, the parameters involved in the elasticity-phase-field formulation are specifically defined through multiscale analysis with non-separated scales. This particular setting ensures consistency with critical subscale information, and allows for the propagation of stochasticity at the macroscopic level. Similar ideas were pursed in the very recent work [START_REF] Acton | Voronoi tessellation based statistical volume element characterization for use in fracture modeling[END_REF], with a few noticeable differences though. First, the approach developed in the above reference is concerned with dynamical fracture, solved using an asynchronous spacetime discontinuous Galerkin method, and is focused on fracture strength random fields. A phase-field approach to brittle fracture modeling is alternatively considered and extended in this paper, in which validation is further assessed on a macroscopic quantity of interest. Second, and while both contributions invoke information theory as a rationale to define probability measures, stochastic modeling aspects and related methodological issues are addressed more extensively hereinafter. Note also that crack paths are simulated in the sequel by propagating a pre-existing crack, whereas crack nucleation sites are identified, for each sample of the microstructure, as the weakest material points in [START_REF] Acton | Voronoi tessellation based statistical volume element characterization for use in fracture modeling[END_REF].

This paper is organized as follows. The computational approach enabling the description of crack propagation at the microscopic scale is first detailed in Section 2. The phase-field formulation is introduced, and reference stochastic computations are performed on a prototypical random microstructure. The description of crack propagation using mesoscopic descriptors is next addressed in Section 3. Methodologies to define and subsequently identify these descriptors are introduced, and applied to the aforementioned microstructure. The relevance of the framework is finally assessed by comparing macroscopic predictions based on either the reference, microscopic model or the proposed mesoscopic approach.

STOCHASTIC MODELING OF CRACK PROPAGATION AT MICROSCALE

Formulation

We consider an open bounded domain Ω ⊂ ℝ , = 2, 3, with boundary Ω , subjected to Dirichlet and Neumann boundary conditions and on the corresponding boundaries Ω and Ω respectively, with Ω ∪ Ω = Ω and Ω ∩ Ω = ∅ (see Fig. 1 (a)). The total energy of the cracked solid is given by

( , Γ) = ∫ Ω∖Γ ( , Γ) Ω +  -1 (Γ) -∫ Ω ⋅ Γ , (1)
where denotes the displacement vector, is the elastic strain density function, is the fracture toughness, or critical energy release rate in the sense of Griffith and  the Hausdorff measure of dimension -1 representing the discontinuous crack shape Γ. In this work, we adopt the variational approach to fracture as proposed by Bourdin, Francfort and Marigo [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] and developed in an efficient computational framework in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], called in the literature the phase field approach to fracture. The above energy form can be replaced by a regularized one [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Tortorelli | Approximation of functional depending on jumps by elliptic functional via t-convergence[END_REF], given by:

( , ) = ∫ Ω ( , ) Ω + ∫ Ω ( ) Ω -∫ Ω ⋅ Γ , ( 2
)
where is a crack density function, whose model can be chosen among several possible forms, leading to a class of shapes for the regularized damage field near the crack (see, e.g., [START_REF] Borden | A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework[END_REF]). In this regularized framework, the cracks are no more described by surfaces but by a smooth field  → ∈ [0, 1] (see Fig. 1 (b)). For example, the so-called first-order model is associated with

( ) = 1 2 2 + 2 ∇( ) ⋅ ∇( ) , (3)
where ∈ ℝ + is a regularization length associated to the regularized description of the cracks (see Fig. 1 (c)) such that ≪ Ω. Discussions about the choice of can be found in, e.g., [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF]. In Eq. ( 2), the strain density function can be decomposed as follows to involve damage in tension only and then avoiding self interpenetration while taking into account automatically self-contact within the crack [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]:

Ψ = (( ) + )Ψ + ([ + ]) + Ψ -([ -]) , (4)
where  → ( ) denotes the (stiffness) degradation function. Here, the quadratic function ( ) = (1-) 2 is used and 0 < ≪ 1 is a small regularizing parameter ensuring the well-posedness of the boundary value problem (see e.g. ( 23) and a discussion for the choice of degradation functions in [START_REF] Kuhn | On degradation functions in phase field fracture models[END_REF]. Furthermore, Ψ + and Ψ -denote parts of the strain density related to tensile and compressive parts of the strain tensor, respectively. For an isotropic medium, it leads to (see [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]):

Ψ ± ([ ]) = 2 ⟨ ([ ])⟩ ± 2 + [ ± ] , (5)
where (.) is the trace operator and  → ( ) and  → ( ) are the field of Lamé's constants. In Eq. ( 5), [ + ] and [ -] are such that [ ] = [ + ] + [ -] and are defined by the spectral decomposition

[ ± ] = ∑ =1 < > ± ( ) ⊗ ( ) , ( 6
)
in which {( , ( ) )} =1 are the pairs of associated eigenvalues and eigenvectors of the strain tensor [ ], and < ⋅ > ± is the operator given by Variation of (2) with respect to and leads to the coupled equations:

< > ± = 1 2 ( ± | |) , ∀ ∈ ℝ . ( 7
⎧ ⎪ ⎨ ⎪ ⎩ ⋅ [ ( , )] = , ( -2 Δ ) -2(1 -)([ ]) = 0, (8)
where ∇ ⋅ (.) is the divergence operator and  is a strain density history function, used to prescribe damage irreversibility [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]:

( , ) = max ∈ [0, ] Ψ + ( , ) . ( 9
)
The stress tensor for the damaged material reads, for an isotropic material, as:

[ ( , )] = Ψ [ ] = (( ) + ) ⟨ ([ ])⟩ + [1] + 2 [ ] + + ⟨ ([ ])⟩ -[1] + 2 [ ] -.
The system of equations ( 8) is complemented by the following boundary conditions

⎧ ⎪ ⎨ ⎪ ⎩ = on Ω , [ ] = on Ω , ⋅ n = 0 on Ω , (10)
where ∇(.) is the gradient operator, and are prescribed vector fields of displacements and tractions, and is the outwardpointing normal vector on Ω (see Fig. 1 (a)).

It should be noticed at this point that Eq. ( 5) encapsulates the stochastic aspect of the propagation at microscale, since the field  → ( ( ), ( )) corresponds to a realization of the elasticity tensor random field { (( ( ), ( ))) , ∈ Ω} given by

( ) = ∑ =1 1 Ω ( ) ( ( ), ( )) , ( 11
)
where {1 Ω ( ), ∈ Ω} and ( ( ), ( )) are the indicator function and elasticity tensor of phase (assumed to be isotropic here), occupying domain Ω , and the domain Ω is assumed to contain constitutive phases. For random microstructures, the indicator functions are, indeed, non-Gaussian random fields: in practice, it is thus required to proceed to, e.g., Monte Carlo simulations of these fields, and to solve the coupled elasticity-phase-field problem for each realization of Ω. This strategy allows the variability in crack paths (and consequently, in the nonlinear part of the macroscopic response) to be simulated at microscale. This point is specifically illustrated in the next section.

Computing Realizations of Crack Paths at the Microscale

In this section, we describe the Monte Carlo approach used to generate realizations of crack paths at the micro scale. These realizations will be used, in Section 3.3, to identify a mesoscale stochastic model constructed in Section 3.2.1. In order to illustrate the methodology throughout this paper, a prototypical stationary, isotropic random microstructure made up of a homogeneous matrix and monodisperse spheres is selected hereinafter ( = 2). A two-dimensional square domain Ω = (]0, [) 2 is considered (that is, = 2), with = 1 mm, and the radius of the inclusions is set to = 0.04 × . Plane strain conditions are assumed. A set of = 1000 realizations was generated using the molecular-dynamics-type algorithm proposed in [START_REF] Zamponi | Comment on "packing hyperspheres in high-dimensional Euclidean spaces[END_REF]. Each sample is periodized and contains = 50 non-overlapping heterogeneities. Independent realizations of this microstructure are shown in Fig. 2 .

FIGURE 2

Four independent realizations of the periodized random microstructure.

The constitutive materials are assumed isotropic, and the bulk and shear moduli are denoted as ( , ) and ( , ) for the matrix and inclusions, respectively. Accordingly, , and , denote the toughnesses of the matrix and inclusions. These properties are chosen such that the mechanical contrast satisfies = ∕ = ∕ = , ∕ , , where the properties of the matrix are taken as = 175 [GPa],

= 81 [GPa] and , = 2.7×10 -3 [kN.mm -1 ]. Dirichlet boundary conditions are applied in the form ( 1 , 0) = and ( 1 , ) = (1) for 0 ⩽ 1 ⩽ , in which (1) = (1, 0) is the first vector of the canonical basis in ℝ 2 (note that Ω = ∅) and ∈ [0, ]; see Fig. 3 . An initial crack is positioned as described in Fig. 3 . These boundary conditions correspond to a pure shear loading (see Fig. 3 ), and the evolution of the damage field and displacementforce curve (associated with the microstructural sample shown in Fig. 3 ) can be seen in Figs. 4 and5 , for = 0.0075 within the matrix phase. Since the crack path Γ is, by definition, identified as the collection of points Γ for which ≈ 1, and upon restricting the analysis to configurations containing a single crack, the variability in the crack propagation generated by the underlying microstructural randomness can be observed by considering the stochastic process

{ Γ 1 ( 2 ), 2 ∈ 2 }, in which 2 ⊆ [0, 0.5].
Likewise, the stochasticity induced on the macroscopic response can be characterized by computing the horizontal force on the top edge ( 2 = ) of the samples, denoted by . The mean and standard deviations for these quantities of interest, together with sample-based envelopes are shown in Fig. 6 . In these figures,

 → ( ) = { ( )} and  → ( ) = √ { ( ) 2 } -{ ( )} 2
represent the mean and variance functions for the macroscopic force, and the statistical estimators and envelopes are obtained using 100 independent realizations of the microstructure. Similar notations are used for studying the second-order properties of the process { Γ 1 ( 2 ), 2 ∈ 2 }.

STOCHASTIC MODELING OF CRACK PROPAGATION AT MESOSCALE

In this section, we propose a simplified model of crack propagation at a mesoscale, i.e. at a level of scale where the fully detailed microstructure is replaced by a smoother approximation to reduce the computational costs in a stochastic framework. Note that (i) this mesoscale model differs from a fully homogenized one, as it maintains some statistical fluctuations raised by microstructural randomness; (ii) the formulation remains predictive to study crack propagation, in contrast with an approach that would describe cracks through first-order, averaged characteristics (such as crack density).

Formulation for Stochastic Crack Propagation at Mesoscale

In order to characterize crack propagation using coarse descriptors, we introduce an upscaled version of the elasticity-phase-field problem as follows:

⎧ ⎪ ⎨ ⎪ ⎩ ⋅ [ ( , )] = , ̃ ̃ ( -̃ 2 Δ ) -2(1 -) Ψ+ ([ ]) = 0 , ( 12
)
where the stress tensor

[ ( , )] = (( ) + ) Ψ+ ([ ( )]) [ ( )] + Ψ-([ ( )]) [ ( )] , ( 13
)
is here expressed as a function of a new, mesoscopic stored energy function

Ψ± ([ ]) = [ ± ] ∶ ̃ ∶ [ ± ] , ( 14
)
in which  → ̃ ( ) is a mesoscopic elasticity tensor field, ̃ represents an equivalent toughness for the mesoscale medium and ̃ is the characteristic length associated with the regularized description at the mesoscale. In what follows, the definition of these quantities is investigated through a two-step methodology:

• First, the definition of the elasticity field  → ̃ ( ) is achieved using a moving-window upscaling approach under different types of boundary conditions. This point is discussed in Section 3.2.1.

• Second, the definition of the toughness ̃ is addressed in Section 3.2.2, where a statistical inverse problem involving the peak force at the macroscopic scale is introduced.

The results from the first step will be used, in Section 3.3, to construct a stochastic surrogate for the elasticity field. This model will enable us to draw additional samples of the mesoscopic elasticity field without having recourse to the homogenization solver.

Definition of the Mesoscopic Coefficients 3.2.1 Homogenization-based Definition of Mesoscopic Elasticity Tensor Fields

In this section, we define the technique used to construct a mesoscopic (smooth) definition of the heterogeneous elastic medium from fully detailed realizations of microstructures. Such an approach has been extensively discussed over the past two decades. Within a multiscale setup, this can be achieved by using a local homogenization (see [START_REF] Ostoja-Starzewski | Microstructural Randomness and Scaling in Mechanics of Materials[END_REF] for a survey) or a filtering approach [START_REF] Yvonnet | A consistent nonlocal scheme based on filters for the homogenization of heterogeneous linear materials with non-separated scales[END_REF][START_REF] Bignonnet | Macroscopically consistent non-local modeling of heterogeneous media[END_REF][START_REF] Tran | Stochastic modeling of mesoscopic elasticity random field[END_REF]. For averaging-type upscaling, kinematic and static uniform boundary conditions (which are denoted by KUBC and SUBC hereinafter) can be considered [START_REF] Ostoja-Starzewski | Random field models of heterogeneous materials[END_REF] and provide bounds for the field of apparent tensors [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF][START_REF] Hazanov | Order relationships for boundary condition effect in heterogeneous bodies smaller than the representative volume[END_REF]. Alternatively, periodic boundary conditions (PBC) can be invoked, especially when a fast convergence toward the effective properties is sought; PBC were employed in the so-called moving-window approach (39), for instance. In the sequel, kinematic uniform boundary conditions (KUBC) and static uniform boundary conditions (SUBC) are selected. While square-shaped domains are typically used in the literature of homogenization, a circular moving window is considered to prevent the generation of spurious anisotropic features at the mesoscale created by the corners of a square bow (see, e.g., [START_REF] Bignonnet | Macroscopically consistent non-local modeling of heterogeneous media[END_REF]). In this work, an original method using a moving window for obtaining a smoothed homogenized field with fluctuations is described in the following.

Let Ω be a circular domain of radius ̃ , centered at ∈ Ω, with boundary Ω . For one realization of the microstructure, { ̃ ( ) , ∈ Ω} is obtained by performing a homogenization locally in Ω , for both KUBC and SUBC (see Fig. 7 ). As moves within Ω, we obtain a continuous, smooth homogenized medium characterized by a wavelength associated with the radius ̃ (see illustrations of this process in Fig. 10 for different radii ̃ . We recall that KUBC correspond to the following boundary conditions:

( ) ( ) = [ ( ) ] , ∀ ∈ Ω , ( 15
)
where in the present 2D plane strain context the indices and run over {1, 2}, leading to:

[ (11) ] = 1 0 0 0 , [ (22) ] = 0 0 0 1 , [ (12) ] = 1 2 0 1 1 0 . ( 16
)
The strategy is schematically depicted in Fig. 7 .

For SUBC, traction vectors are applied in the form

( ) = [Σ ( ) ] ( ) , ∀ ∈ Ω , ( 17
)
where ( ) is the outward-pointing normal vector at point ∈ Ω . Combinations of indices similar to those introduced for KUBC are considered, with [Σ (11) ] = 1 0 0 0 , [Σ (22) ] = 0 0 0 1 , [Σ (12) ] = 1 2

0 1 1 0 . ( 18
)
Let { ̃ ( ) , ∈ Ω} and { ̃ ( ) , ∈ Ω} be the random fields of mesoscopic elasticity tensors defined under the aforementioned boundary conditions. Note that when approaches the boundary Ω, the realization of the microstructure is virtually replicated, by periodicity, and the mesoscopic tensor is still well defined.

Denoting by the radius of the inclusions, the ratio ̃ ∕ plays an important role in defining a continuous transition from the microscale ( ̃ ∕ → 0 + ) to the macroscale ( ̃ ∕ → +∞). Additionally, it specifies the level of anisotropy exhibited by the local apparent elasticity tensor, ranging from microscopic isotropy to mesoscopic anisotropy, and then to macroscopic isotropy (in the present case). The definition of the isotropic approximation for the mesoscale elasticity tensor field is addressed in the FIGURE [START_REF] Rudoy | Domain decomposition method for crack problems with nonpenetration condition[END_REF] Defining mesoscale fields of stiffness through homogenization (case of KUBC).

following.

Let { ̃ ( ) , ∈ Ω} denote the isotropic approximation of the field { ̃ ( ) , ∈ Ω}, where the subscript BC refers to the type of boundary conditions under consideration (KUBC or SUBC). In two-dimensional elasticity, and assuming plane strain conditions, the Voigt-type matrix representation {[ ̃ ( )], ∈ Ω} of the aforementioned field is given by

[ ̃ ( )] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ̃ ( ) + 4 3 ̃ ( ) ̃ ( ) -2 3 ̃ ( ) 0 ̃ ( ) -2 3 ̃ ( ) ̃ ( ) + 4 3 ̃ ( ) 0 0 0 ̃ ( ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ∀ ∈ Ω , ( 19
)
where { ̃ ( ), ∈ Ω} and {̃ ( ), ∈ Ω} are the random fields of three-dimensional bulk and shear moduli defining the isotropic approximation, in plane strain elasticity, of the actual (anisotropic) elastic tensor ̃ ( ). To obtain these coefficients, we minimize the distance (in the sense of the metric defined as follows) between ̃ ( ) and [ ̃ ( )] (see [START_REF] Guilleminot | Stochastic modeling of anisotropy in multiscale analysis of heterogeneous materials: A comprehensive overview on random matrix approaches[END_REF][START_REF] Tran | Stochastic modeling of mesoscopic elasticity random field[END_REF] and the references therein):

( ̃ ( ), ̃ ( )) = argmin >0, >0 ‖[ ̃ ( )] -[ ̃ ( )]‖ 2 , ( 20
)
in which

[ ̃ ] = ⎛ ⎜ ⎜ ⎝ ̃ 11 ̃ 12 ̃ 13 ̃ 12 ̃ 22 ̃ 23 ̃ 13 ̃ 23 ̃ 33 ⎞ ⎟ ⎟ ⎠ , ( 21
)
where the Voigt's notation has been used for the different components of the tensor and ‖ ⋅ ‖ is the Frobenius norm.

The optimization problem defined by Eq. ( 20) can readily be solved by a direct differentiation of the cost function, and the mesoscopic moduli of the isotropic approximation are found as

̃ ( ) = 1 60 (11 ̃ 11 ( ) + 11 ̃ 22 ( ) -4 ̃ 33 ( ) + 38 ̃ 12 ( )) , ̃ ( ) = 1 5 ( ̃ 11 ( ) + ̃ 22 ( ) + ̃ 33 ( ) -2 ̃ 12 ( )) , ( 22
)
where the subscript BC has been dropped for notational convenience (this convention will be used in the sequel when no confusion is possible). In order to assess the relevance of the isotropic approximation, the following random field { ( ), ∈ Ω} is introduced (40, 35):

 ( ) = ‖[ ̃ ( )] -[ ̃ ( )]‖ F ‖[ ̃ ( )]‖ F , ∀ ∈ Ω . ( 23
)
The graphs of the fields of mean and standard deviation evaluated for a coarse mesoscopic grid (with 20 points along each direction) are shown in Figs. 8 and9 , for both KUBC and SUBC.

It is seen that the error between the homogenization-based random fields and their isotropic approximations remain small in mean and variance, and that the error is larger in the case of KUBC. The effect of the mesoscopic resolution on the isotropic approximation (for KUBC) is qualitatively shown on a single realization in Fig. 10 . As expected, the field becomes more homogeneous as the ratio ̃ ∕ increases. The approximation in the set of isotropic tensors is satisfactory and allows the phase-field approach to be readily applied in the isotropic case. In this rest of this paper, we will consider the characterization The graphs of the first-order marginal probability density functions for the bulk and shear moduli are shown in Fig. 11 . The well-known ordering with respect to boundary conditions is observed almost surely (that is, ̃ ⩽ ̃ and ̃ ⩽ ̃ for each microstructural sample), and it is seen that the level of statistical fluctuations associated with KUBC is larger than for SUBC.

FIGURE 11

Graphs of the first-order marginal probability density functions for the bulk (left) and shear (right) moduli.

The estimated normalized correlation functions along (1) and (2) for the fields of bulk and shear moduli (for SUBC) are shown in Fig. 12 . In these figures, the notation  → ̃ ( ) indicates that the correlation function of the bulk modulus random field is evaluated along the unit vector ( ) (a similar notation is used for the shear modulus). It is seen that the differences between the

FIGURE 12

Graph of the correlation function along (1) and (2) , estimated from the simulated data, for the random fields of bulk (left) and shear (right) moduli.

correlation functions for the two random fields are almost indistinguishable, due to the very strong cross-correlation between the two properties. Moreover, it can be observed that the correlation first decreases over the range [0, ∕2] (with = 1) and then starts increasing on [ ∕2, ], in accordance with the periodicity of the underlying background medium. This information will be used to select an appropriate form of the correlation functions for the random field models, constructed in Section 3.3.1.

Identification of the Mesoscopic Toughness

Let us now turn to the identification of the fracture and phase-field parameters ̃ and ̃ at the mesoscale (see Section 3.1). In this work, the characteristic length ̃ of the mesoscopic medium is set to be equal to the characteristic length at the microscale, previously denoted by . This choice is supported by unreported numerical convergence results, as well the results provided in [START_REF] Nguyen | Identification of fracture models based on phase field for crack propagation in heterogeneous lattices in a context of non-separated scales[END_REF]. The mesoscopic toughness parameter ̃ is next identified by solving an inverse problem involving the peak force = max ( ). More precisely, ̃ is calibrated by imposing a match between the mean value of , estimated with 500 independent microstructural samples and the fine-scale elasticity-phase-field simulations (detailed in Section 2.2), and the mean value determined with the mesoscopic description introduced in Section 3.1, denoted by ̃ . In the latter description, samples of the elasticity random field are estimated through the mowing-window homogenization procedure, and the mesoscopic toughness then appears as the unique unknown parameter. Let ̃ ( ̃ * ) be the mean peak force associated with the candidate value ̃ * for the mesoscopic toughness (that is, by substituting ̃ * for ̃ in Eq. ( 12)). An optimal value can thus be defined by minimizing the relative error function

 ( ̃ * ) = | -̃ ( ̃ * )| (24)
over the admissible set [ , , , ]:

̃ = argmin ̃ * ∈ [ , , , ]  ( ̃ * ) . ( 25
)
Since the mesoscopic elasticity field depends on the boundary conditions applied, the above optimization problem must be solved independently for KUBC and SUBC. The graph of the cost function obtained for a resolution parameter = 3 is shown in the left panel in Fig. 13 . The optimal values are obtained as

̃ = 1.022 × , ≈ 2.75 × 10 -3 [kN.mm -1 ] ( 26
)
for KUBC and

̃ = 1.255 × , ≈ 3.375 × 10 -3 [kN.mm -1 ] (27)
for SUBC. A comparison of the macroscopic responses obtained with the mesoscopic formulations (for KUBC and SUBC), parametrized with the identified values, and the reference computations is shown in the right panel in Fig. 13 . This figure shows that the mesoscale formulation identified under SUBC provides a fairly accurate estimate of the mean macroscopic response. This conclusion similarly holds for KUBC, although this type of boundary conditions leads to a stiffening of the response, in accordance with the fact that apparent tensors obtained under KUBC constitute upper bounds for the mesoscopic elasticity. In both cases, the variability is underestimated due to filtered elasticity fluctuations. These effects can clearly observed in Fig. 13 , where the mean and standard deviation on crack paths are reported for the two types of boundary conditions. It is seen that the cracks paths obtained with the mesoscopic formulations are localized near the mean crack path at microscale, for both KUBC and SUBC. Not surprisingly, the mesoscopic-based crack paths exhibit a variability that is much smaller than the one obtained at microscale (where the crack trajectory is constrained by the radius of the inclusions). While the results support the relevance micro KUBC of the formulation, the latter necessitates solving a very large number of homogenization problems to represent the fluctuations of the elasticity field at the mesoscopic scale. In the next section, we address the construction of a stochastic model that enables the elasticity tensor random field to be sampled in a robust manner.

Γ micro Γ d SU BC Γ g SU BC

Stochastic Modeling and Identification of Mesoscale Elasticity

Construction of Random Field Models

The construction of stochastic models for random fields of elasticity tensors with arbitrary material symmetries has been investigated in [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF][START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF][START_REF] Staber | Stochastic modeling and generation of random fields of elasticity tensors: A unified information-theoretic approach[END_REF] using an information-theoretic formulation [START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial derivative operators[END_REF] (see [START_REF] Guilleminot | Non-Gaussian Random Fields in Multiscale Mechanics of Heterogeneous Materials:1-9[END_REF] for a survey), and in (46) using a spectral expansion. For a given type of boundary conditions, let { ̃ ( ), ∈ Ω} and {̃ ( ), ∈ Ω} be the random fields of three-dimensional bulk and shear moduli as defined in Eq. [START_REF] Acton | Voronoi tessellation based statistical volume element characterization for use in fracture modeling[END_REF].

Let { ( ) = (Ξ 1 ( ), Ξ 2 ( )), ∈ Ω} be a bivariate Gaussian field with statistically independent, normalized components. These components are defined by the correlation functions ( , )  → Ξ 1 ( , ; (1) ) and ( , )  → Ξ 2 ( , ; (2) ). The non-Gaussian random fields of elastic moduli are then defined through the nonlinear transformations

̃ ( ) = -1 ( ̃ , ̃ ) (0,1) (Ξ 1 ( )) (28)
and

̃ ( ) = -1 ( ̃ , ̃ ) (0,1) ( Ξ 1 ( ) + √ 1 -2 Ξ 2 ( )) , (29)
where -1 ( , ) is the inverse cumulative distribution function of the Gamma law with (shape and scale) parameters and (note that the dependence of these parameters on the boundary conditions is not reported for notational convenience), (0,1) is the cumulative distribution function of the standard Gaussian law and denotes the correlation coefficient between ̃ ( ) and ̃ ( ).

Given the stationarity and the form of the correlation functions estimated for the random fields of elastic moduli (see Fig. 12 ), the following separable form is retained:

Ξ ( , ; ( ) ) = ( 1 ; ( ) 1 ) × ( 2 ; ( ) 2 ) , ∀ ∈ ([0, ]) 2 , ∈ {1, 2} , (30)
where = | -| is the lag distance along ( ) , ∈ {1, 2}, and the one-dimensional normalized correlation function  → ( ; ) is defined as

( ; ) = exp - 2 2 sin 2 . ( 31
)
In Eq. ( 31), is a model parameter related to the internal length

 = ∕2 ∫ 0 | ( ; )| , (32)
which is interpreted, in the periodic setting under consideration, as the spatial correlation length of the Gaussian random field along the associated basis vector (e.g., along (1) if the function

1  → ( 1 ; ( ) 1 ) is considered). It can be shown that the correlation length  < ∕2 reads as  = 2 exp{--2 } 0 ( -2 ) , ( 33
)
where  0 denotes the zero-order modified Bessel function. It should be noted that by construction, one has < ∕2. The graph of  → ( ; ) is shown in Fig. 16 for different values of . The following properties can easily be deduced.

• The first-order marginal probability measure is a bivariate Gamma law [START_REF] Patrick | Statistical inference with bivariate gamma distributions[END_REF][START_REF] Arnst | An overview of nonintrusive characterization, propagation, and sensitivity analysis of uncertainties in computational mechanics[END_REF], which is consistent with previous results derived within the framework of information theory (see [START_REF] Guilleminot | Non-Gaussian Random Fields in Multiscale Mechanics of Heterogeneous Materials:1-9[END_REF] and the references therein).

• The mean values ̃ = { ̃ ( )} and ̃ = {̃ ( )} read as

̃ = ̃ × ̃ , ̃ = ̃ × ̃ , (34)
and the coefficients of variation are given as

̃ = 1 √ ̃ , ̃ = 1 √ ̃ . ( 35
)
Recall that these properties are independent of location , owing to the stationarity of the random fields.

• The fields of stiffness and compliance tensors are of second-order:

{‖[ ̃ ( )]‖ 2 } < +∞ , {‖[ ̃ ( )] -1 ‖ 2 } < +∞ , ∀ ∈ Ω , ( 36
)
hence ensuring that the stochastic linear elastic boundary value problem is well posed [START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial derivative operators[END_REF].

• The random fields { ̃ ( ), ∈ Ω} and {̃ ( ), ∈ Ω} are mean-square continuous and mean-square differentiable.

From a computational standpoint, the underlying Gaussian fields are sampled using a truncated Karhunen-Loève expansion. In order to reduce the associated computation time, the random fields are sampled on a grid that is coarser than the one used to solve the elasticity-phase-field problem at the mesoscale. Realizations of fields are then obtained by interpolating on the fine mesoscopic grid. To that end, the coarse mesoscopic grid is specifically defined so that the correlation structure is properly discretized. In the results presented hereinafter, the coarse mesh includes six Gauss points per correlation length, in each direction.

Identification of the Elasticity Random Field

The probabilistic model involves two sets of parameters controlling the joint probability density function of the elastic moduli at a given location and the correlation structure of the underlying Gaussian fields. The first set of parameters gathers ( ̃ , ̃ ) and ( ̃ , ̃ ) (or equivalently ( ̃ , ̃ ) and (̃ , ̃ ), in view of Eq. ( 34)), as well as the the coefficient of cross-correlation . These hyperparameters can be estimated from the database using standard statistical estimators, here with 500 sample realizations:

̃ = 258.
) 38
The kernel density estimations of the first-order marginal and joint distribution for the bulk and shear moduli obtained with the data and the model-based samples are shown in Figs. 17 and 18 (recall that the resolution is fixed by ̃ ∕ =3 here). A very good agreement is observed between the probability density functions corresponding to the homogenization-based data and those estimated with model-based samples. Since the transformations given by Eqs. ( 28) and ( 29) are nonlinear, the correlation functions associated with the random fields of elastic moduli cannot be inferred explicitly. In this case, the identification of the vector-valued hyperparameters (1) and (2) is performed through a two-step procedure, described below.

Let the correlation functions of { ̃ ( ), ∈ Ω} and {̃ ( ), ∈ Ω} be written as  → ̃ ( ; (1) ) and  → ̃ ( ; ( (1) , (2) )), respectively: this notation emphasizes the underlying dependence on the parameters of the Gaussian fields (see Eqs. ( 28) and ( 29)). By a slight abuse of notation, these correlation functions will also be denoted as  → ̃ ( ; (1) )

and  → ̃ ( ; ( (1) , (2) )) when evaluated along ( ) , ∈ {1, 2}. In a first stage, the components of (1) are identified, for a given type of boundary conditions, as

= argmin > 0  (1) 1 ( ) , (1) 1
= argmin > 0  (1) 2 ( ) , ( (1) 2
) 39
where the cost functions are given by and

 (1) 1 ( ) = ⎛ ⎜ ⎜ ⎝ ∕2 ∫ 0 ( ̃ ( 1 ) -̃ ( 1 ; )) 2 1 ⎞ ⎟ ⎟ ⎠ 1∕2 (40)
 (1) 2 ( ) = ⎛ ⎜ ⎜ ⎝ ∕2 ∫ 0 ( ̃ ( 2 ) -̃ ( 2 ; )) 2 2 ⎞ ⎟ ⎟ ⎠ 1∕2 . ( 41
)
In a second stage, the hyperparameters controlling the correlation structure of the Gaussian random field {Ξ 2 ( ), ∈ Ω} are identified as

(2) 1 = argmin > 0  (2) 1 ( ) , (2) 2
= argmin > 0  (2) 2 ( ) , (42)
in which

 (2) 1 ( ) = ⎛ ⎜ ⎜ ⎝ ∕2 ∫ 0 ( ̃ ( 1 ) -̃ ( 1 ; ( (1) 1 , )) 2 1 ⎞ ⎟ ⎟ ⎠ 1∕2 (43)
and

 (2) 2 ( ) = ⎛ ⎜ ⎜ ⎝ ∕2 ∫ 0 ( ̃ ( 2 ) -̃ ( 2 ; ( (1) 2 , )) 2 2 ⎞ ⎟ ⎟ ⎠ 1∕2 . ( 44
)
In Eqs. ( 43) and ( 44), the values of (1) 1 and (1) 2 are those obtained within the first step of the methodology (see Eq. ( 39)). Below, the estimations of the correlation functions associated with the stochastic model are performed by combining statistical and ergodic averaging on a set of 500 independent realizations of the random fields (for fixed values of the hyperparameters). The optimal values are found as (1) = (0.4624, 0.4574) , (2) = (0.4043, 0.4014)

for the elasticity random fields identified under KUBC, and

(1) = (0.4654, 0.4604) , (2) = (0.4694, 0.4654)

for the case of SUBC.

(τ 2 ) R model k (τ 2 ) R gen k (τ 2 )
FIGURE [START_REF] Acton | Voronoi tessellation based statistical volume element characterization for use in fracture modeling[END_REF] Normalized correlation function of the bulk modulus random field along (1) (left) and (2) (right), estimated from the multiscale data (red line) and the calibrated stochastic model (black line).

that the calibrated model allows the decays of the correlation functions to be accurately reproduced, which is key to mimicking the mesoscopic elasticity (and in particular, the frequency of sample path oscillations that has a substantial impact on the crack paths in the phase-field formulation at mesoscale).

Comparison Between the Microscopic and Mesoscopic Formulations

The variability in crack paths can be observed in Fig. 21 for the reference microscale-based computations and the mesoscopic formulation defined with either locally-homogenized microstructural samples or the elasticity field stochastic model. It can be observed that the crack paths corresponding to a description at microscale present larger statistical fluctuations. As previously indicated, this result is indeed expected, since the propagation only occurs in the matrix phase then (hence forcing the crack to get around inclusions). In contrast, the mesoscopic formulations exhibit much smaller variability but capture quite accurately the mean crack path. Interestingly, it is seen that the stochastic model for the elasticity tensor random field performs well in delivering crack paths that are consistent those computed with homogenization-based fields. The predictions of the macroscopic force-displacement curve are shown in Fig. 22 . In accordance with the results presented in Section 3.2.2 (see the right panel in Fig. 13 ), where microscopic and homogenization-based mesoscopic formulations were compared, it is observed that the mesoscopic elasticity-phase-field formulation defined under SUBC delivers predictions in better agreement with the reference solution, in terms of both the mean elastic response and mean peak force. In contrast, the The error measure is given by ≈ 0.084% for KUBC and ≈ 0.037% for SUBC, showing that an accurate prediction of the mean peak force can be obtained with the two types of boundary conditions.

CONCLUSION

A stochastic approach to model crack propagation in random media has been proposed in this paper. The formulation relies on a phase-field formulation where material coefficients are defined and identified through multiscale computations. Monte-Carlo simulations were first performed using a description at the microscopic scale. These computations enables the characterization of subscale-induced randomness on the macroscopic response of the domain and were subsequently used as reference results to assess the relevance of the framework. The definition of the mesoscopic parameters was then addressed. The elasticity field at mesoscale was specifically defined as the isotropic approximation of spatially dependent homogenized tensors, obtained by means of a moving-window upscaling approach (under kinematically and statically uniform boundary conditions). The (deterministic) mesoscopic toughness was identified by solving an inverse problem related to the mean peak force. It is shown that the formulation under statically uniform boundary conditions allows for an accurate prediction of the mean elastic response and mean peak force. In contrast, kinematically uniform boundary conditions generate a stiffening of mesoscale elasticity, in accordance with theoretical results derived elsewhere. An information-theoretical probabilistic model for the elasticity random field was then constructed and allows for a fast, robust sampling of mesoscopic elasticity. The results obtained by feeding this stochastic surrogate model into the phase-field formulation were finally compared with those corresponding to the full-scale, microscopic model. It is shown, in particular, that the model-based, mesoscopic elasticity-phase-field formulation associated with statically uniform boundary conditions allows for an accurate prediction of both the mean elastic response and mean peak force.

FIGURE 1

 1 FIGURE 1 Panel (a): sharp description of a cracked solid. Panels (b), (c) and (d): smeared description within the phase field framework, showing damage profile  → ( ) for different values of the regularization length (in this figure, denotes the characteristic length of the domain, here of the minor axis of the ellipsoidal shape).

FIGURE 3 FIGURE 4

 34 FIGURE 3 Boundary conditions applied to the domain and initial crack.

FIGURE 5 Force 2 FFIGURE 6

 526 FIGURE 5 Simulated displacement-force response for the microstructure shown in the left panel in Fig. 3 .

FIGURE 8 FIGURE 9

 89 FIGURE 8 Graphs of the mean functions for the random fields { ( ), ∈ Ω} (left) and { ( ), ∈ Ω} (right).

FIGURE 10

 10 FIGURE 10 One realization of the random field {‖[ ( )]‖ , ∈ Ω} for ̃ ∕ ∈ {0.2, 0.6, 0.8, 1, 2, 4} (from left to right).

FIGURE 13

 13 FIGURE[START_REF] Bažant | Probabilistic Nonlocal Theory for Quasibrittle Fracture Initiation and Size Effect. I: Theory[END_REF] Graphs of the cost function for KUBC and SUBC (left), and envelopes of the macroscopic response for the reference microscale model and the mesoscopic formulations (right).

FIGURE 14

 14 FIGURE 14 Comparison of crack paths obtained with the microscopic description (black dashed/solid lines) and the mesoscopic formulation. Left panel: case of KUBC. Right panel: case of SUBC.

FIGURE 15

 15 FIGURE 15 Crack path for the same microstructure: microscale description (left), KUBC-based description (middle), SUBCbased description (right).

2 FIGURE 16

 216 FIGURE 16 Plot of the correlation function  → ( ; ) for different values of .

7 [

 7

FIGURE 18

 18 FIGURE[START_REF] Rahman | Stochastic multiscale fracture analysis of three-dimensional functionally graded composites[END_REF] Kernel density estimates for the joint probability density function of the bulk and shear moduli, for homogenization-based (left) and model-based (right) samples (for SUBC).

FIGURE 22

 22 FIGURE 22 Envelopes of the macroscopic response for the reference microscale model ( ) and the mesoscopic formulations parametrized either by upscaled coefficients ( ̃ , ̃ ) or model-based coefficients ( ̃ , ̃ ).

  Kernel density estimates for the probability density function of the bulk (left) and shear (right) moduli, for the two types of boundary conditions KUBC (blue) and SUBC (red).

	1.8	10 -11					3.5	10 -11
	1.6						3
	1.4
							2.5
	1.2
	1						2
	0.8						1.5
	0.6
							1
	0.4
	0.2						0.5
	0						0
	1.5	2	2.5	3	3.5	4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	2.2
						10 11									10 11
	FIGURE 17

  The graphs of the normalized correlation functions estimated with the data and with the stochastic model thus identified (under SUBC) are shown in Figs.19 and 20for the bulk and shear moduli random fields, respectively. It is seen

	1										1
	0.8					R data k k R gen k (τ 1 ) (τ 1 ) (τ 1 ) R model				0.8						R data k
	0.6										0.6
	0.4										0.4
	0.2										0.2
	0										0
	-0.2										-0.2
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
					Lag distance									Lag distance

FIGURE [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] Normalized correlation function of the shear modulus random field along (1) (left) and (2) (right), estimated from the multiscale data (red line) and the calibrated stochastic model (black line).

FIGURE 21

Comparison of crack paths obtained with the microscopic description (black dashed/solid lines), the mesoscopic formulation where the elasticity is obtained microstructural samples (blue lines), and the mesoscopic formulation involving the stochastic model for elasticity tensors (red lines). Left panel: case of KUBC. Right panel: case of SUBC.

use of KUBC at the mesoscale generates a stiffer response in the elastic regime, while still allowing for a good prediction of the mean peak force. In addition, the use of the stochastic model in lieu of homogenization-based samples does not introduce any significant bias in the predictions, regardless of the type of boundary conditions. Upon interpreting the macroscopic force as a stochastic process indexed by the prescribed displacement, the error generated by the proposed model-based, mesoscopic formulation can be characterized as

where the stochastic process { ̃ ( ), ∈ [0, ]} is implicitly depending on the boundary conditions applied at mesoscale. The right-hand side term in Eq. 47 can be estimated through Monte Carlo simulations, and the relative error remains small for the two types of boundary conditions, with ≈ 1.6% for KUBC and ≈ 2.2% for SUBC. The error for the prediction of the mean peak force can be characterized by )