Phase-field modeling of brittle fracture in heterogeneous bars

We investigate phase-field modeling of brittle fracture in a one-dimensional bar featuring a continuous variation of elastic and/or fracture properties along its axis. Our main goal is to quantitatively assess how the heterogeneity in elastic and fracture material properties influences the observed behavior of the bar, as obtained from the phase-field modeling approach. The results clarify how the elastic limit stress, the peak stress and the fracture toughness of the heterogeneous bar relate to those of the reference homogeneous bar, and what are the parameters affecting these relationships. Overall, the effect of heterogeneity is shown to be strictly tied to the non-local nature of the phase-field regularization. Finally, we show that this non-locality may amend the ill-posedness of the sharp-crack problem in heterogeneous bars where multiple points compete as fracture locations.


Introduction
Although at some scale all materials are heterogeneous, i.e. their properties vary in space, the adoption of homogeneous effective properties is often sufficient for mechanical modeling in engineering. However, the effect of heterogeneity cannot be ignored for a large class of mechanical problems such as those involving composite materials, biological tissues and metamaterials. The evolution of cracks in these materials follows complex patterns that challenge many modeling and computational approaches.
Phase-field modeling of brittle fracture was proposed by Bourdin et al. [1] as the regularization of the variational fracture formulation by Francfort and Marigo [2] and was later re-interpreted as a special family of gradient damage models [3]. It provides a remarkably flexible variational framework to describe the nucleation and propagation of cracks with arbitrarily complex geometries and topologies in two and three dimensions [4].
The original phase-field modeling approach is based on the assumption of homogeneous elastic and fracture properties of the material throughout the domain. Previous studies addressing phase-field modeling in heterogeneous materials adopt a pragmatic approach, by simply substituting the constant fracture toughness of the original model with a fracture toughness depending on the material point [5][6][7][8]. Natarajan et al. [5] propose a phase-field formulation for fracture in functionally graded materials. The approach is further developed by Kumar et al. [6], where it is shown that the peak stress of a functionally graded material remains bounded between the values pertaining to the single constituents in homogeneous conditions. Hossain et al. [7] propose a technique based on phase-field modeling to evaluate the effective fracture toughness of heterogeneous media, while Shen et al. [8] show that the introduction of a spatially variable fracture toughness in phase-field models is a promising tool to model fracture in bones. However, to the best of our knowledge the implications of heterogeneous material properties on the key predictions of the phase-field model have never been the object of a fundamental investigation. Thus, the relationship between local material properties and observed behavior as predicted by the phase-field model remains unclear, which in turn may prevent the proper calibration of the model and the proper interpretation of its results.
In this work, we perform such investigation for the one-dimensional case. We revisit the fundamental mathematical analysis in [3] by assuming that the elastic and/or fracture material properties are heterogeneous with different possible profiles of spatial variations. We aim at quantitatively assessing how the heterogeneity in the material properties influences the observed behavior, and especially the peak stress and the fracture toughness, in the context of phase-field modeling.
The paper is structured as follows. Section 2 formulates the one-dimensional phase-field model of brittle fracture and the related evolution problem. Section 3 defines the classes and profile shapes of heterogeneity adopted in the subsequent sections. Section 4 briefly reviews the solution of the evolution problem for the homogeneous bar. The core of the study is Section 5, where the evolution problem is solved for the heterogeneous bar. The analysis is first carried out in the one-dimensional space and then extended to a bar in the three-dimensional space. In Section 6, we discuss the consequences of heterogeneity on the fracture behavior as predicted by the phase-field approach in contrast to the predictions of the sharp-crack model. The main conclusions are drawn in Section 7.
In the following, the dependence on the pseudo-time t of the quasi-static setting is denoted with the subscript t, e.g. α t is the damage variable at pseudo-time t; the prime symbol denotes the derivative with respect to either the spatial coordinate x, e.g. u ′ t = ∂u t /∂x, or the damage variable α, e.g. w ′ = dw/dα; the symbols ∇(•) and ∆(•) represent respectively the gradient and the Laplacian of the vectorial quantity (•) with respect to the spatial coordinates, while the divergence is indicated as div(•); the dot symbol denotes the derivative with respect to the pseudo-time, e.g.α t = ∂α t /∂t; vectors and second-order tensors are denoted with bold symbols, e.g. u is the displacement vector in three-dimensions; the 2nd-order identity tensor is indicated with I and the symmetric part of the 4th-order identity tensor with I s ; the superscript T denotes the transpose of a matrix.

One-dimensional phase-field model for brittle fracture
In this section, we formulate the phase-field model of brittle fracture for a one-dimensional domain (Figure 1), along the lines of [3,[9][10][11] but with some generalizations to prepare for the later developments. The primary unknowns, both functions of the spatial coordinate x, are the displacement u and the phase-field or damage variable α. The latter is an internal variable which describes the material damage level. Its magnitude is bounded between α = 0, corresponding to a sound material, and α = 1, denoting a fully damaged material.

Energetic quantities
As follows, we introduce some definitions which will be used in the remainder of this paper, especially concerning important energetic quantities. The total energy density W is defined as where ϕ el is the elastic energy density and ϕ d is the dissipation density. The elastic energy density is given by where E 0 > 0 (considered here as a continuous function of x to account for possible heterogeneity in the elastic properties of the material) is the undamaged elastic modulus and a(α) is the degradation function. The latter describes the degradation of the elastic modulus due to damage, thus it is a monotonically decreasing function such that a(0) = 1 and a(1) = a ′ (1) = 0. We also introduce the compliance modulation function s(α) as the reciprocal of a(α), The dissipation density reads hence it consists of a local term, depending on the damage variable, and a non-local term, depending on its spatial derivative. In the local term, w(α) is the local dissipation function, a monotonically increasing function of α such that w(0) = 0 and w(1) = 1. There are two common options for the definition of the functions a(α) and w(α) [3,12]: AT2: a(α) = (1 − α) 2 and w(α) = α 2 .
Throughout this paper, we will focus on the AT1 model. In the non-local term, the dependency on the spatial derivative of the damage variable calls for the introduction of the internal length parameter ℓ. The local magnitude of the dissipation density is modulated by a specific fracture energy w 1 , which is considered a continuous function of x to account for possible heterogeneity in the fracture properties of the material. The total energy functional reads and it is the sum of the elastic energy functional E el (u, α) and of the dissipation functional D(α), with (9) and

Evolution problem
Let us now consider a bar clamped at the left end and loaded with a prescribed displacement at the opposite end (Figure 1), where U t is a positive smooth function of the pseudo-time t. At pseudo-time t, a displacement field v and a damage field β are admissible if they respectively belong to U t and A, with A more precise definition of these functional spaces is out of the scope of this work. We limit ourselves to require for them a sufficient regularity so that the total energy functional remains finite. The evolution of the system can be studied as a quasi-static process parameterized through the pseudo-time t ≥ 0 and described with the function t → (u t , α t ) and can be characterized variationally by means of total energy minimization. Following the variational approach in [9][10][11], the evolution problem is governed by the principles of irreversibility, local stability and energy balance and can be formulated as follows: Problem 1 (Evolution problem). Given the initial state (u 0 , α 0 ) = (0, 0) at the pseudo-time t = 0, find t → (u t , α t ) ∈ U t × A fulfilling the following conditions: 1. irreversibility: t → α t is a non-decreasing function, 2. local stability: 3. energy balance: where is the work made by external actions in the pseudo-time interval [0, t].
In Eq. 16, σ s denotes the Cauchy stress at pseudo-time s:

First-order evolution problem
Upon a first-order expansion of Eq. 14 (under the assumption of sufficient smoothness), an evolution t → (u t , α t ) is a solution of Problem 1 only if it is solution of the first-order evolution problem [9][10][11] Problem 2 (First-order evolution problem). Given the initial state (u 0 , α 0 ) = (0, 0) at the pseudo-time t = 0, find t → (u t , α t ) ∈ U t × A sufficiently smooth fulfilling the following conditions: Remark 1. Eqs. 21, 23, 24 involve the spatial derivative of E 0 (x) and w 1 (x) which do not appear in the classical homogeneous formulation.
Remark 2. The difference between the KKT conditions Eqs. 22-24 for the general case of the heterogeneous bar problem and the analogous conditions for the special case of homogeneous bar is the presence of the terms containing the spatial derivative w ′ 1 . This contribution to the strong form of the governing equations is not included in previous literature dealing with heterogeneous materials, see e.g. [5,6]. However, these studies perform numerical finite element analyses based on the weak form associated to Eqs. 22-24; in the weak form the term containing w ′ 1 is compensated for by a similar term with opposite sign appearing after integration by parts, hence it does not influence results. As will be shown later, in the present study the same term is essential to understand the role played by heterogeneity on qualitative and quantitative aspects of the solution of the evolution problem.

Homogeneous and heterogeneous bars
Thus far, the local elastic and fracture material properties have been characterized through E 0 (x) and w 1 (x), respectively. In the following, we distinguish between two cases: • homogeneous bar : the special case in which E 0 and w 1 are constant along the bar; • heterogeneous bar : the more general case in which E 0 and/or w 1 vary along the bar, as described by the functions E 0 (x) and/or w 1 (x), assumed to be sufficiently regular.
The spatial distribution of the material properties for the heterogeneous bar problem is further defined as where the constantsĒ 0 andw 1 are reference values of the undamaged elastic modulus and the specific fracture energy, respectively, and the functions f E (x) and f w (x) define the corresponding spatial variation profiles. The material with E 0 (x) =Ē 0 and w 1 (x) =w 1 is denoted as the homogeneous material associated to a given heterogeneous material. In the following, all the quantities referred to the associated homogeneous material are denoted with a bar. For the future developments, we now define three shapes of the heterogeneity profile, denoted by h i (x) with i = lin, par, exp, each of which depends on a length ℓ f termed the characteristic length of the heterogeneity ( Table 1). The magnitude of ℓ f characterizes how rapidly the material properties vary along the axis of the bar. This variation is associated to the slope of the profile (i = lin), to its curvature (i = par), or to both (i = exp).
We also define three classes of heterogeneous materials as summarized in Table 2, each class assigning a profile shape to f w (x) and/or to f E (x). Accordingly, we distinguish between heterogeneity in the specific fracture energy (hw), heterogeneity in the undamaged elastic modulus (hE) and full heterogeneity (hwE).

Remark 3.
We assume the material properties to be minimum at the midpoint cross-section of the bar and we choose symmetric increasing profiles of three different shapes ( Table 1). As a result, the midpoint cross-section is the weak location where we expect damage to start and develop first. We will compare the behavior of these heterogeneous bars with that of homogeneous bars where the properties are everywhere equal to the minimum values, i.e. E 0 = min x E 0 (x) andw 1 = min x w 1 (x).
Finally, we introduce the dimensionless coordinatex := x/ℓ. With a slight abuse of notation, we denote , respectively. In particular, the profile shapes are written in terms ofx as follows: where r is the characteristic ratio: The limit case r → 0 is obtained for: • the sharp-crack model: ℓ → 0 • the associated homogeneous material: ℓ f → ∞ Since we are interested in the diffusive approximation of cracks in heterogeneous materials, the relevant range of values for the applications we have in mind is 0 < r < 1, i.e. the variation of the material properties occurs over lengths that are sufficiently larger than the intrinsic length scale of the phase-field model. Values of r (much) larger than 1 would represent a variation of material properties occurring at a length scale below the size of the characteristic length. Due to the nature of the regularization, such variations would not be "seen" by the model and are therefore not relevant for the present study.

Shape
Expression Plot Tab. 1: Shapes of heterogeneity: linear shape, parabolic shape and exponential shape. Type with i = lin, par, exp Tab. 2: Classes of heterogeneity: heterogeneity in specific fracture energy (hw), heterogeneity in undamaged elastic modulus (hE), full heterogeneity (hwE).

Solution of the evolution problem for the homogeneous bar
In this section, we briefly summarize the solution of the evolution problem formulated in Section 2 for the special case of homogeneous bar with undamaged elastic modulusĒ 0 and specific fracture energyw 1 . This solution has been thoroughly analyzed in the literature, see e.g. [9][10][11] . As mentioned earlier, we limit ourselves to the case of the AT1 model (Eq. 5). This model satisfies the strain hardening condition, i.e. −w ′ (α) /a ′ (α) is increasing with respect to α (this implies that the elastic domain in the strain space expands for increasing damage), and the stress softening condition, i.e. w ′ (α) /s ′ (α) is decreasing with respect to α for all values of α in [0, 1] (this implies that the elastic domain in the stress space shrinks for increasing damage) [3].
Let us consider an initially unstrained and undamaged bar, i.e. (u 0 , α 0 ) = (0, 0), loaded with an imposed end displacement U t as introduced in Section 2.2. The starting point of the analysis is the construction of a homogeneous solution, i.e. a solution characterized by a constant value of the damage variable along the bar under a monotonically increasing prescribed displacement. Since the stress is constant due to equilibrium and the elastic properties are homogeneous, the strain field is also constant and given by u ′ t = U t /2L. The bar being initially undamaged, the solution of the evolution problem is characterized by an initial elastic phase, where the damage criterion in Eq. 23, which in the case at hand simplifies to is a strict inequality with α t = 0. This phase continues until the applied displacement U t reaches its value at the elastic limit corresponding to the elastic limit stress or yield stress We denote the corresponding pseudo-time as t e . For t > t e (hence U t > U e ), the damage criterion (30) becomes an equality and damage can grow. For the ensuing homogeneous solution in the damaging phase the homogeneous value of the damage variable can be computed from the prescribed displacement U t through and the corresponding stress is Note that the validity of the strain hardening condition guarantees that the functional relationship α t → U t in Eq. 32 is monotonically increasing, hence there is a unique α t solution for a given U t , whereas the stress softening property implies that the stress in Eq. 33 decreases with the damage level and hence with the applied displacement. Thus, the peak stress of the homogeneous response is reached for α t = 0, hence Denoting the pseudo-time at peak stress as t p , for the homogeneous bar it is thus t p = t e . A stability analysis [9][10][11] demonstrates that for a sufficiently long bar (2L ≫ l) the homogeneous state is unstable for any U t ≥ U e and a damage localization necessarily arises at the end of the elastic phase. In this localized solution, damage is only non-zero within an open interval S t ∈ (−L, L) where the damage criterion holds as an equality, while it vanishes in the remainder of the domain. Infinite localized solutions are possible based on the position of the damage localization region within the domain. Without loss of generality, we assume here that this region is centered at x * = 0, i.e. at the midpoint of the bar. A thorough analysis of the localization phase can be found in [9][10][11] and references therein and is not repeated here. During localization, the maximum value of the damage variable, i.e. α t (x * ), increases monotonically whereas the stress decreases, hence σ p is the peak stress not only of the homogeneous response but of the overall stress-displacement response. To follow this phase, control is switched from increasing prescribed displacement to decreasing stress or increasing α t (x * ), as the corresponding U t may no longer be monotonically increasing depending on the length of the bar.
At the end of the localization phase, α t (x * ) reaches the value 1 leading to failure of the bar. We denote the corresponding pseudo-time as t u , at which the stress σ u = 0 and the fully localized damage profile α u , symmetric about x * = 0, reads [12]: where δ u = 2ℓ is the half-support width.
The fracture toughness is defined as the dissipated energy at failure: and in the present case is given by Remark 4. According to Eq. 37, for a given ℓ, knowledge of the local quantitȳ w 1 is sufficient to determine the global quantity G c .

Solution of the evolution problem for the heterogeneous bar
In this section, we derive the solution of the evolution problem for the heterogeneous bar, including the homogeneous and the localized solutions, and especially focusing on the effect of the heterogeneity on peak stress and fracture toughness. As in the case of the homogeneous bar, we carry out the analysis for the AT1 model.

Governing equations on half domain
We can take advantage of symmetry and study the problem on half (e.g. on the left half) of the domain. For later reference, we rewrite here the equilibrium equation and the KKT conditions 1. irreversibility:α 2. damage criterion: 3. loading-unloading conditions: The governing equations above do not require the existence of f ′ w (0) but only of the left-hand derivative. The natural boundary conditions read: We now need additional boundary conditions atx = 0. These can be easily retrieved from the variational approach as for the case of the homogeneous bar [9], and read as follows depending on t:

Homogeneous solution
Using Eqs. 5 and 17, the damage criterion Eq. 40 takes the form With the dimensionless coordinate defined in Section 3, the dimensionless damage criterion reads whereσ t is the dimensionless stresš andσ e =σ p = Ē 0w1 is the yield stress, equal to the peak stress, for the associated homogeneous material.
As in the analysis for the homogeneous bar, we first look for a homogeneous solution in the elastic phase, where the damage criterion is satisfied as a strict inequality with α t = 0. It is straightforward to determine the dimensionless elastic limit stressσ e and the positionx * of the first point of the bar reaching the elastic limit as followš x * = arg miň Remark 5. According to Eq. 49, the elastic limit stress for the heterogeneous bar is the same as for the bar made of the associated homogeneous material and is equal to Next, we look for a homogeneous solution in the damaging phase, where the damage criterion is satisfied as an equality with α t = 0 and uniform along the bar. Assuming uniform damage delivers Eq. 52 only admits a solution for the special case where f E (x) · f w (x) is constant along the bar, which is excluded a priori by our choice of the heterogeneity profiles in Section 3. Hence, the evolution problem for the general case of a heterogeneous bar does not admit a homogeneous solution in the damaging phase.

Localized solution
After reaching the elastic limit, i.e. for t > t e , the heterogeneous bar problem admits only a localized solution with α t (x) = 0. Hence, within the left half-domain there exists an interval (−δ t , 0) where the damage criterion holds as an equality while the remainder of the bar is undamaged, i.e.
withδ t := δ t /ℓ, where δ t is the half-support width (a priori unknown) at pseudo-time t. The regularity of the functions a(α) and w(α) implies that α t and α ′ t are continuous within (−L/ℓ, 0) [9], hence while the boundary conditions Eqs. 44, 45 continue to hold. The damage profile α t (x) is assumed to be monotonically increasing over (−δ t , 0) (an assumption that can be easily verified a posteriori), therefore the maximum damage value is the value of the damage variable atx * = 0, i.e. α t (0).

Peak stress and stress-displacement curve during localization
Next, we show that the boundary value problem in α t constituted by Eq. 53 along with the boundary conditions in Eqs. 54, 44 for t ∈ [t e , t u ) admits solutions for increasing values of the stress, up to a peak value that, as in the case of the homogeneous bar, we denote as the peak stress. After the peak, the stress decreases. As follows, we devise a semi-analytical scheme to solve the problem and determine the peak stress.
The boundary value problem can be reformulated as an initial value problem through the spatial coordinate transformation : whereδ in t is a guess for the a priori unknown half-support widthδ t . The damage criterion is rewritten in terms ofx as The initial value problem defined by Eqs. 56, 57 is solved via a Runge-Kutta scheme using the algorithm ODE 45 of Matlab [13] starting from the initial pointx = 0 ( Figure 2). The integration is stopped at the target pointx =δ out t = 0 such that For a givenσ t with t > t e , we assign toδ t the values corresponding to the conditionδ in t =δ out t along the piecewise linear interpolation of the pairs (δ in t ,δ out t ) obtained from Eqs. 56-58. When t = t e , the bar is undamaged anď δ t = 0. For 1 ≤σ t ≤σ p , two values are assigned toδ t . The peak stressσ p is found as the stress for which these two values coincide (Figure 3). Details about the numerical implementation are reported in Appendix A.
Since the computation is based on Eqs. 56-58, the result, i.e. the dimensionless peak stressσ p , depends on the functions f w (x) and f E (x), i.e. it depends on the heterogeneity class and profile shape and, for a given class and profile shape, it depends on the characteristic ratio r only. Results for all the considered heterogeneity classes and profile shapes are illustrated in Figure 4.
Remark 6. For the heterogeneous bar, the peak stress σ p is larger than the elastic limit stress σ e (Figure 4). This is due to the presence of a short hardening phase during the initial damage localization process.
Computing the support extension for different values ofσ t with t ∈ [t e , t u ) enables also the definition of the full stress-displacement curve during the damage localization phase. In this case, the numerical integration is performed via the stiff equation solver ODE 23s of Matlab [13], see Appendix A for the detailed algorithm. To eachδ t we can associate a damage profile α t (x), its maximum α * t =α t (δ t ), and the stressσ t . Also, recalling Eq. 17, by knowing the current dimensionless stressσ t and the damage profile α t , the dimensionless applied displacement U t /2L can be computed through the integral Sorting these quantities based on an ascending order of α * t yields the stressdisplacement curve during the localization phase ( Figure 5).

Remark 7.
In the spirit of previous studies on phase-field modeling of brittle fracture [3,14] and in compliance with the Γ-convergence arguments at the root of the approach, we performed here path-independent energy minimization, i.e. we do not account for the irreversibility condition. As noted in [9], the solution stemming from an incremental procedure and compatible with the irreversibility condition corresponds to the upper envelope of the set of localization profiles obtained for t ∈ (t e , t u ]. for anyδ in t . The blue and red lines represent the profile shape h i and the damage variable α t , respectively. Numerical integration starts fromx = 0 and is stopped atx =δ out t . Fig. 3: Representation of the piecewise linear interpolation of pairs (δ in t ,δ out t ) for increasing values ofσ t starting fromσ e = 1 (dotted lines). The plot is obtained for the hE class with profile shape f (x) = 1 − r · (x −δ in t ). Points e and p correspond to the solutions for pseudo-times t = t e and t = t p , respectively.

Fracture toughness
The fracture toughness is defined as the dissipated energy at failure, and its determination requires the computation of the fully localized damage profile. At t = t u it is σ t = 0 and the damage criterion within the halfsupport width reads It is immediate to notice that, since function f E is only contained in the stress term which is now zero, the damage profile at t = t u for the heterogeneous bar only depends on function f w . This implies that the fracture toughness for classes hw and hwE depends on the profile shape of the heterogeneity and on the characteristic ratio r, whereas a heterogeneity of class hE leaves the fracture toughness unchanged. As follows, we exemplify the computations for the linear heterogeneity profile shape, whereas the analogous computations for the parabolic and exponential profile shapes follow similar lines and are reported in Appendix C and Appendix D, respectively. For the linear profile shape, Eq. 60 becomes The analytical solution for α u depends on the unknown coefficients c 1 and c 2 : The two unknown coefficients can be obtained as functions of the unknowň δ u combining Eq. 62 with the boundary conditions Eq. 54 leading to By the following substitutions z = 8 r 2 − 1 exp(1) and τ = 8 and using Eqs. 62, 63, 64, we can rewrite the remaining boundary condition Eq. 45 as (Appendix E) z = τ · exp(τ ).
Eq. 66 has solution where W k is the Lambert function whose definition is provided in Appendix B. Substituting backwards and prescribingδ u ≥ 0 we finď which gives the half-support widthδ u as a function of the characteristic ratio r (Figure 6a). The dimensionless fracture toughnessǦ c , defined aš whereḠ c = 8/3 ℓw 1 is the fracture toughness for the associated homogeneous material, can be obtained recalling Eqs. 36 and 10 aš Combining Eq. 62 with Eqs. 63 and 64 and using f w (x) = h lin (x) in Eq. 28 1 , G c can be expressed in terms ofδ u aš A further substitution of Eq. 68 in Eq. 71 leads to an analytical expression forǦ c that is plotted as solid line in Figure 6b. The expression forǦ c can be simplified using the polynomial approximation of the Lambert function proposed by Veberič [15] along with a Taylor expansion about r = 0. Although the approximation order can be freely selected (Appendix F), here the Veberič approximation is truncated at order 6 and the Taylor expansion at order 5 givinǧ The exact dimensionless fracture toughness (Eq. 71) and its polynomial approximation (Eq. 72) are compared in Figure 6b. Finally, combining Eqs. 62, 63, 64, 68 we can plot the damage profile at failure, α u , for different values of the characteristic ratio r (Figure 7). The results on the dimensionless fracture toughness for the three different profile shapes are summarized in Figure 8.

Remark 8.
For the heterogeneous bar, the damage profile at failure is narrower and the fracture toughness G c is larger than for the bar made of the associated homogeneous material (Figures 7, 8).

Extension to heterogeneous bars in three dimensions
In this subsection, we extend the previous one-dimensional analysis to the case of bars with properties varying along the longitudinal axis embedded in the three-dimensional space. The main aim is to investigate the effects of this more realistic setting on the peak stress and fracture toughness and highlight the differences with the one-dimensional case (Sections 5.3.1 and 5.3.2). Figure 9 illustrates the geometry and the boundary conditions of the problem. The bar has length 2L and square cross-section with dimensions 2H × 2H, and a monotonically increasing displacement U t is prescribed at x = L. Homogeneous Neumann boundary conditions are enforced to the damage variable on the whole boundary.

Problem setting
The three-dimensional version of the total energy functional in Eq. 7 reads where u is the displacement vector, C 0 is the undamaged 4th-order elasticity tensor for linear elastic isotropic materials, ε(u) = 1 2 (∇u + ∇ T u) is the 2ndorder infinitesimal strain tensor, Ω is the spatial domain and x = {x, y, z} is the spatial coordinate vector with x corresponding to the longitudinal axis of the bar.
The distribution of the elastic and fracture properties along the major axis is defined again through the functions f E and f w as whereC 0 andw 1 are independent of x. In Appendix G it is shown that, for isotropic material properties, Eq. 74 implies a heterogeneous distribution of the undamaged elastic modulus, namely E 0 (x) =Ē 0 · f E (x), and a homogeneous Poisson's ratio ν(x) =ν. The main hypotheses adopted in the one-dimensional case are preserved. In particular, we perform a path-independent energy minimization (Remark 7), we introduce the dimensionless spatial coordinate vectorx = x/ℓ and the dimensionless stress tensorσ t = σ t / w 1Ē0 . We limit the analysis to a linear type of heterogeneity with slope r (Table 1).
Equilibrium prescribes div(σ t ) = 0, while the damage criterion is written as whereŠ t is the support of the damage variable.Š 0 is the dimensionless 4th-order undamaged compliance tensor, which readš and depends on the Poisson's ratio ν only. Note that Eq. 75 depends exclusively on the parameters r and ν and on the class of heterogeneity (hw, hE, hwE). Thus, in a three-dimensional setting the damage criterion depends on the Poisson's ratio and the criterion defining the elastic limit t e becomes maxx(Š 0σe (x) ·σ e (x)) = 1.

Peak average stress
In three dimensions the stress can vary within the cross-section. However, equilibrium imposes a constant axial force, hence a constant cross-sectional average normal stress σ t , which is defined as where σ xx t is the xx-component of the Cauchy stress tensor σ t = a(α t ) C 0 ε(u t ). Therefore, in the three-dimensional context the concept of peak stress used in Section 5.3.1 is replaced by the average peak stress σ p = sup t σ t . Thus, let us study the peak average dimensionless stress σ p = σ p / w 1Ē0 . We perform a series of numerical experiments for different value of r ∈ [0, 1] and four values of the Poisson's ratio ν = {0, 0.15, 0.3, 0.45}. We assume to have a long bar with L/ℓ = 15, aspect ratio L/H = 6.25, ℓ = 0.1,Ē 0 = 1 and w 1 = 1. The finite element analyses are carried out with our code GRIPHFiTH [16] using a structured mesh with 8-node brick elements. Along y and z, the mesh size is 0.02. The mesh is refined in the central zone along the x-axis for a length of 4 ℓ, i.e. where the damage is expected to localize. In this central region, the element size is 0.01, while outside it is 1.3/7. The load U t is applied in steps with increments of ∆U t = 5 · 10 −4 . The non-negativity of the damage variable is enforced through the penalty method proposed by Gerasimov et al. [12] and the problem is solved by means of an alternate minimization scheme [17] which is stopped when the norm of the residuals becomes lower than a predefined tolerance of 10 −4 .
The results in terms of σ p for the three classes of heterogeneity are summarized in Figure 10 and are qualitatively very similar to those obtained in the one-dimensional setting, showing an increase of the peak stress with the characteristic ratio r. Quantitatively we observe that the values obtained for ν = 0 overlap with those of the one-dimensional bar (as expected), while increasing the Poisson's ratio leads to an increase of the peak average stress σ p . For the classes hE and hwE, this is due to a varying lateral contraction introduced by the Poisson's effect whose magnitude is modulated by the elastic modulus along the bar. As a consequence, a necking deformation is introduced already in the elastic regime whose magnitude varies along the x-axis and has its maximum in correspondence of the minimum value ofĒ 0 (x). It is worth noting that this necking deformation is not present in homogeneous bars, at least before the onset of the failure after reaching the peak stress σ p = σ e . The change in lateral contraction introduces parasitic shear stresses that break up the uniaxial stress state, leading to a variation of the termŠ 0σt (x) ·σ t (x) within the cross-section and, thus, to a point-wise different local elastic limit. This effect is more pronounced for larger values of the Poisson's ratio and is further exacerbated by the progressive evolution of the damage during the hardening regime (see Section 5.3.1), which makes the heterogeneity in the elastic parameters more pronounced. As a result of this mechanism the damage does not evolve homogeneously within the cross-section, and this ultimately leads to an increase of the peak average stress with respect to the one-dimensional case.
Although to a limited extent, this effect is also observed in absence of heterogeneity in the elastic properties, as demonstrated by the results for class hw (Figures 10 and 11a). In this case, the elastic limit is the same for all the points within the cross-section where localization is expected, however the uniaxial stress state is altered by the evolution of the damage during the hardening regime, which promotes a spatial variation of the elastic parameters and an inhomogeneous distribution of the damage variable within the cross-section (Figure 11b).

Fracture toughness
At failure, the left-hand side of Eq. 75 vanishes since σ u = 0, thus, the dependence of the problem on the Poisson's ratio disappears and the threedimensional problem simplifies in a scalar one-dimensional-like problem. As a consequence the dimensionless fracture toughness depends only on r and on the function f w and the same observations reported in Section 5.3.2 apply also here.
This result is also verified numerically. After the peak load, the system enters a softening branch characterized by a snap-back where the damage evolves abruptly until reaching the complete failure of the bar, i.e. α u = 1 at x = 0. At this point the damage is homogeneous within the cross-section and its profile obtained numerically coincides with the analytical result derived in the one-dimensional case (Figure 12).

Discussion
In this subsection we discuss some implications of the obtained results.

Non-locality
In the present section we have investigated how heterogeneity in the elastic and fracture material properties affects the observed behavior of a bar. We have assumed the material properties to be minimum at the midpoint cross-section of the bar (we have taken these as reference material properties). As a result, the midpoint cross-section is the location where the elastic limit is reached first and around which damage localization starts and develops, finally leading to failure of the bar.
We have concluded that, for a bar with a given heterogeneity class and profile shape, the peak stress σ p and the fracture toughness G c can be expressed as whereσ p andḠ c denote the peak stress and the fracture toughness of the bar made of the homogeneous reference material, and p(r) : r →σ p and g(r) : r →Ǧ c have been determined to be always larger than 1 and increasing with r. Thus the peak stress and fracture toughness of the heterogeneous bar are both larger than those of the homogeneous reference bar. This increase is a non-local effect resulting from the elastic modulus and/or the specific fracture energy being larger than those of the reference homogeneous material in the neighborhood of the midpoint cross-section, thus it is a consequence of the choice made for the reference homogeneous material. The non-locality is naturally induced by the phase-field model through its intrinsic length scale, so that the macroscopic behavior of the bar does not simply result from the local properties at the midpoint cross-section but involves its neighborhood. Accordingly, the increase in peak stress and fracture toughness is a function of the ratio r between the internal length of the phase-field model and the length characterizing the speed of variation of the material properties. The non-local effect vanishes, i.e. p(r) and g(r) approach 1, when the limit case r → 0 is approached, i.e. for the sharp-crack model (or, trivially, for homogeneous material properties).

Model calibration
For calibration of the phase-field model, different options are possible depending on which properties can be realistically assumed to be known. Let us assume that the shape of the heterogeneity in the material properties is known upfront, e.g. through computed tomography by correlation with the density profile. This is a common practice in many fields, e.g. bone biomechanics [18,19] (however correlation is typically assumed between the density and the value of the elastic modulus, whereas the analogous correlation with the specific fracture energy is less investigated). Under this assumption, ℓ f is also known. To fix ideas, let us further assume that the heterogeneity shape corresponds to one of the profile shapes we considered in this study. Quantities which can be realistically measured on the bar geometry are the initial stiffness k 0 , the elastic limit stress σ e and the peak stress σ p . The initial stiffness can be expressed as from which the value of the elastic modulus at the midpoint cross-section,Ē 0 , can be deduced. From Eq. 51,w 1 can be computed from the measurement of σ e . Finally, from the measurement of σ p , the intrinsic length ℓ of the phasefield model can be calibrated using Eq. 78 1 and recalling thatσ p =σ e = σ e .

Sharp-crack model vs. phase-field model
Let us now explore further the consequences of the non-local nature of the phase-field model, as opposed to the locality of the sharp-crack model, in bars made of heterogeneous materials.

Sharp-crack model
The sharp-crack model, put forth by Francfort and Marigo [2] as a variational reformulation of Griffith's brittle fracture criterion, relies on the global minimization of a total energy functional. In the one-dimensional case, this functional can be expressed as [12] E sc (u, Γ c ) := where Ω is the problem domain, Γ c is the crack set, i.e. the set including the cracked points of the bar, and H 0 (Γ c ) is its Hausdorff measure which, in the one-dimensional case, returns the number of points belonging to Γ c .
The productḠ c f w (x) can be regarded as the specific fracture energy for the sharp-crack model. Accordingly, also in this case f w (x) plays the role of spatial variation profile of the fracture property andḠ c represents the minimum value of the specific fracture energy. While the phase-field model is based on local energy minimization (see the local stability condition in Problem 1), the sharp-crack model needs the minimization of the total energy functional to be global, otherwise crack nucleation would not be predicted since the undamaged elastic solution is always locally stable. Within this global minimization framework, an initially uncracked bar under tension remains sound until the creation of a crack becomes energetically more convenient than the storage of additional elastic energy. At critical conditions, the stored elastic energy is released completely and the fracture energy takes the value corresponding to a single crack formed in correspondence of the minimum value of f w (x) [20]. As claimed in [2], this result marks a contrast with classical fracture mechanics based on Griffith's criterion, as it amends its inability to predict crack initiation.

6.2.
Multiple minima for f w An interesting benchmark showcasing the differences between sharp-crack and phase-field models is the case of bars where multiple points compete as initiation sites, namely when multiple minima for f w are present. In the following, we present two different examples, one with two equal minima and one with two different minima.

Two equal minima for f w
We consider heterogeneity in specific fracture energy with a profile f w (x) possessing two equal minima in x 1 and x 2 . We further assume that the profile can be split into two parts, one symmetric about x 1 and the other symmetric about x 2 , and we consider linear and parabolic profiles as in Figure 13. We also introduce the function f w,1 (x − x 1 ) describing the right-half of the first part of the profile and f w,2 (x − x 2 ) describing the right-half of the second part of the profile. For the examples in Figure 13, we have: Let us first study this problem using the sharp-crack model. Gerasimov et al. [21] showed that, for this profile, the problem with the sharp-crack model is ill-posed due to the competition between the two possible crack locations x 1 and x 2 . They proposed a stochastic relaxation to transform the ill-posed deterministic problem into a well-posed stochastic problem formulated in terms of fracture probability. The stochastic solution was found by introducing a random perturbation to the specific fracture energy profile, in form of a white noise with magnitude controlled by the small parameter η > 0, and then letting η approach 0. Denoting with P i the probability that the crack forms at x i (i = 1, 2), it was concluded that, for the linear and parabolic examples in Figure 13, P 1 = 1/3 and P 2 = 2/3 [21].
In [21], the mentioned probabilities are calculated numerically but it is possible to obtain the same results in closed form. For a given parameter η controlling the magnitude of the noise, the number of favorable cases for fracture at x i is proportional to f −1 w,i (1 + η), whereas the number of possible cases is proportional to 2 j=1 f −1 w,j (1 + η). Therefore, Eq. 82 can be easily extended to the case where the profile is not symmetric and to an arbitrarily large number of minima. Both for linear and parabolic examples, Eq. 82 together with Eq. 81 yields P 1 = 1/3 and P 2 = 2/3 which coincide with the results in [21]. Let us now solve the same problem using the phase-field modeling approach. Exploiting the results in the present study, it is straightforward to demonstrate that the problem becomes well-posed. For the profiles in Figure 13, for a given value of ℓ, the characteristic ratio r 1 about x 1 is larger than the characteristic ratio r 2 about x 2 . Therefore, according to the result in Figure 4, fracture at x 1 requires a larger stress than fracture at x 2 , hence, fracture can only occur at x 2 . This is confirmed by a one-dimensional finite element analysis performed assuming ℓ = 0.1,Ē 0 = 1,w 1 = 1 and no heterogeneity in the elastic modulus. We use a uniform mesh with linear elements of size ℓ/30. At every time step, the prescribed end displacement is increased by a constant increment ∆U t = 0.1. Irreversibility is enforced through the penalty method [12] and the same alternate minimization scheme mentioned in the threedimensional analysis is adopted. In Figure 14, the results at failure are shown. As expected, in both cases the crack forms at point x 2 .
Thus, in this example the non-local regularization inherent to the phasefield model amends the ill-posedness of the problem. Since the non-local model "sees" not only the minimum in the specific fracture energy but also its neighborhood, equal minima surrounded by different neighborhoods (which reflects in unequal values of r) lead to different peak stresses and the phasefield regularized problem becomes well-posed.

Two different minima for f w
Let us now study an example with two different minima at x 1 and x 2 ( Figure 15). This time, according to the sharp-crack model the problem is well-posed. Indeed, the crack is predicted to initiate at x 1 , at which f w (x) possesses its global minimum [20].
In order to solve this problem with the phase-field approach, we assume that the material is heterogeneous in w 1 only and that ℓ = 0.1. Accordingly, r 1 = 0.5 and r 2 = 0.1125. We can then estimate the peak stress for a fracture at x i with the relationship σ p i =σ p i Ē 0w1 f w (x i ), whereσ p i is the dimensionless peak stress for linear shape and class hw extracted from Figure 4 at r = r i . We find that σ p 1 ≃ 1.08 Ē 0w1 is lower than σ p 2 ≃ 1.11 Ē 0w1 and hence fracture is predicted at x 2 . To further demonstrate the validity of this result, we perform the finite element analyses illustrated in Section 6.2.1 with the same numerical and material parameters as in the analytical computations. The finite element results, given in Figure 16, are in agreement with the theoretical prediction of a crack at x 2 .
Further, Figure 16 illustrates the normalized peak stresses corresponding to the two minima as functions of the regularization length of the phase-field model ℓ. As ℓ decreases, the two normalized peak stresses decrease at a different rate, so that below a threshold value ℓ * ≃ 0.056 the failure location is shifted from x 2 to x 1 , i.e. to the location predicted by the sharp-crack approach.
Thus, in this example the non-local regularization inherent to the phasefield model may lead to different predictions on the failure location than the sharp-crack model, for which only the minimum value of the specific fracture energy counts. The predicted failure location of the phase-field model depends not only on the minimum value but also on the characteristics of its neighborhood (embodied in the parameter r). The amount of sampled neighborhood depends on the value of ℓ, and the predicted failure location shifts to the one of the sharp-crack model as ℓ becomes sufficiently small. Interestingly, in this case it is the phase-field approach to be ill-posed for ℓ = ℓ * .    17: Peak stresses for fracture at x 1 (blue curve) and at x 2 (red curve) vs. internal length ℓ for the example with two different minima. The actual peak stress corresponds to the minimum among these two options and fracture occurs at the respective point. There is a value of length ℓ * which marks a transition in the predicted failure location, from the one of the sharp-crack model x 1 (ℓ < ℓ * ) to the other one x 2 (ℓ > ℓ * ).

Conclusions
We investigated phase-field modeling of brittle fracture in a heterogeneous one-dimensional bar. We assumed the material properties to be minimum at the midpoint cross-section (taking these minimum values as reference material properties), and chose continuously and symmetrically increasing profiles of different shapes along the axis of the bar. Our main goal was to quantitatively assess how the heterogeneity in elastic and fracture material properties influences the observed tensile strength and fracture toughness of the bar, as obtained from the phase-field modeling approach.
The main findings can be summarized as follows: • The elastic limit stress for the heterogeneous bar is the same as for the bar made of the reference homogeneous material; • The evolution problem for the heterogeneous bar does not admit a homogeneous solution in the damaging phase; • Heterogeneous bars show a hardening branch after the elastic limit that leads to a peak stress larger than the elastic limit stress. The value of the peak stress is influenced by both elastic and fracture properties; for a given class and profile shape of heterogeneity, it only depends on the ratio between the internal length of the phase-field model and the length characterizing the speed of variation of the material properties (characteristic ratio); • The fracture toughness for the heterogeneous bar is larger than for the bar made of the reference homogeneous material and is influenced by the fracture properties only; for a given profile shape of heterogeneity, it only depends on the characteristic ratio.
• The results obtained in the one-dimensional space are also valid for the bar embedded in the three-dimensional space. The only difference is that in the three-dimensional case the peak stress of the bar is also influenced by the Poisson's ratio, with a larger Poisson's ratio making the effect of heterogeneity more pronounced.
The observed effects of heterogeneity are direct consequences of the nonlocal nature of the phase-field model. This becomes evident through the comparison between sharp-crack and phase-field model predictions. Within the sharp-crack modeling framework, fracture in bars with heterogeneous specific fracture energy featuring multiple equal minima is an ill-posed problem and can only be addressed via stochastic relaxation. However, the same problem becomes well-posed with phase-field modeling if the heterogeneity profile features different values of the characteristic ratio at the equal minima. It is also shown that more complex cases of heterogeneity can be easily addressed by directly exploiting the findings in this study. In particular, the critical location among competing minima corresponding to "defects" of different sizes can be easily identified, which seems very relevant to the study of fracture in heterogeneous materials.
localization phase. In both, we collect the input lengthsδ in t in a vectorδ in and the corresponding output lengthsδ out Eq. 70 yields the dimensionless fracture toughnesš G c = 1 720 r 3 (−δ u r (60 + 25δ 2 u r 2 + + 3δ 4 u r 4 ) + 15 1 +δ 2 u r 2 2 4 +δ 2 u r 2 arctan(δ u r)). (C.6) Combining the numerical solution forδ u and Eq. C.6 we obtain the curve of the dimensionless fracture toughness vs. the characteristic ratio in Fig Appendix D. Fracture toughness for f w with exponential profile For f w with the exponential profile shape, the dimensionless damage criterion at t = t u within the half-support width reads We simplify the expression for the dimensionless fracture toughnessǦ c using the approximation of the Lambert function proposed by Veberič [15] and the Taylor expansion about r = 0. Veberič's approximation is truncated at order 6 (Appendix F) and the Taylor expansion at order 5 (Figure D that can be rearranged as follows: Proceeding with the further substitutions z = C exp(1) and τ = C/y 2 we retrieve Eq. 66.