Accounting for weak interfaces in computing the effective crack energy of heterogeneous materials using the composite voxel technique

We establish a computational methodology to incorporate interfaces with lower crack energy than the surrounding phases when computing the effective crack energy of brittle composite materials. Recent homogenization results for free discontinuity problems are directly applicable to the time-discretized Francfort-Marigo model of brittle fracture in the anti-plane shear case, and computational tools were introduced to evaluate the effective crack energy on complex microstructures using FFT-based solvers and a discretization scheme based on a combinatorially consistent grid. However, this approach only accounts for the crack resistance per volume and is insensitive to the crack resistance of the interface which is expected to play a significant role by considerations from materials science. In this work we introduce a remedy exploiting laminate composite voxels. The latter were originally introduced to enhance the accuracy of solutions for elasticity problems on regular voxel grids. We propose an accurate approximation of the effective crack energy of a laminate with weak interface where an explicit solution is available. We incorporate this insight into an efficient algorithmic framework. Finally, we demonstrate the capabilities of our approach on complex microstructures with weak interfaces between different constituents.


Introduction
Griffith [1] laid the foundations of modern fracture mechanics [2] by introducing an energetic criterion for the quasi-static propagation of a crack in a brittle medium.He considered a linear elastic, isotropic structure with a pre-existing crack subjected to an incrementally increasing load.In his model the crack will propagate if the energy release rate reaches a critical value, the so-called crack resistance.
A different perspective was taken by Irwin [3] who concentrated on the r −1/2 -singularity of the stress field in the r -vicinity of the crack tip.For three different crack modes he identified a fracture criterion dependent on the prefactors of the singular terms, the so-called stress intensity factors.In his model a crack will propagate if these stress intensity factors reach a critical value, the so-called fracture toughness.It can be shown that the criteria of Griffith and Irwin are actually equivalent for homogeneous isotropic materials [4], i.e., they can be converted into each other.A third approach to classic linear elastic fracture mechanics is given by the J-integral [5,6], which uses the divergence theorem to evaluate the energy release rate by a line integral in two dimensions and a surface integral in three dimensions, enclosing the crack tip in either cases.
Subsequently, extensions to homogeneous and anisotropic materials were proposed [7][8][9][10].Furthermore, nonlinear effects were incorporated to model small plastic effects near the crack tip [11].In particular, the work of Dugdale et al. [12], who accounted for elasto-plastic material behavior, laid the foundations for cohesive fracture models [13,14].Moreover, some early works on brittle fracture of heterogeneous materials investigate the crack behavior at the interface, i.e., whether an evolving crack penetrates the interface [15,16] or is deflected at the interface [17,18].
To deal with complex geometries, computational methods for fracture mechanics were introduced, see Sedmak [19] for a recent overview.For instance, extended finite element methods [20,21] were introduced to account for the jump in the displacement field within a finite element context.Furthermore, cohesive zone elements [22] were considered which allow to implement specific traction-separation laws [23].
Merely 25 years ago, Francfort and Marigo proposed a variational approach to fracture [24] based on Griffith's crack criterion [1].For a domain ⊂ R 3 and given stiffness tensor (field) C as well as crack resistance (field) γ , the Francfort-Marigo model is based on the functional where u : → R 3 denotes the displacement field which is assumed to be continuous except for across the crack surface S. For prescribed boundary conditions and in a time-incremental framework, the Francfort-Marigo model seeks minimizers of the functional (1.1) under the constraint that the crack surface is not allowed to shrink which prevents the material from healing.The Francfort-Marigo model naturally accounts for heterogeneous materials and an anisotropic stiffness tensor.Furthermore, the authors discuss how to incorporate a direction dependent crack resistance as well.Existence of solutions of the Francfort-Marigo were proven recently [25,26] based on global energy minimization.Although quite a few numerical methods for finding these minimizers were introduced [27][28][29][30], the most widely studied approach employs the phase-field model of brittle fracture [31] which was inspired by the Ambrosio-Tortorelli approximation [32] of the Mumford-Shah functional [33].The latter coincides with the Francfort-Marigo model in the case of anti-plane shear and a single load step except for the irreversibility constraint.The phase-field approach models the crack surface as a smeared interface of width in terms of a scalar damage variable D : → R. The energy functional reads In the quasi-static setting, minimizers u and D of this functional are sought for each time step.Passing to the limit of vanishing crack-width parameter → 0 in the phase field model recovers the Francfort-Marigo model in the sense of -convergence.To account for irreversibility in numerical computations, different numerical approaches were pursued [34][35][36].
Phase-field models for brittle fracture permit to treat both crack nucleation as well as crack propagation in a unified framework and do not require to pre-select the crack path a priori.However, the original model had a few shortcomings which where investigated by subsequent works.A first flaw in the original model was the inability to distinguish tensile and compressive stress/strain states.Tension-compression splittings of the elastic energy [37][38][39] offer an elegant way to overcome this issue.Extensions of the phase-field crack model to incorporate an anisotropic crack resistance, e.g., following the theoretical result of Focardi [40] for Mumford-Shah functionals [33], were investigated based on a second order tensor formulation [41][42][43][44].More complex types of anisotropy may be considered in the framework of higher-order phase-fields using fourth order tensor formulations [45][46][47][48].
Similar to the Francfort-Marigo model, phase-field fracture models naturally account for heterogeneous materials [49][50][51].Of particular interest to this work are contributions which specifically account for interface effects.Yoshioka et al. [52] considered interfaces of finite thickness, selecting the parameters for the crack energy in such a way that the resulting fracture energy of the phase field in the interface is equivalent to the energy of a sharp crack following the interface.Other approaches incorporated weak interfaces by combining cohesive zone elements at the interface with phase field fracture [53,54].Apart from using finite element and finite difference discretizations, phase-field fracture models may also be integrated into FFT-based micromechanics [55], see the pioneering work of Chen et al. [56], where cohesive zone elements are considered at the interface which they implement using the composite voxel approach [57][58][59].
Multi-scale methods [60] permit to account for micro-heterogeneous materials in simulations on the component scale.Homogenization results provide a mathematically rigorous way to provide such a scale transition between the microscopic and the component scale.In the context of fracture mechanics, Braides et al. [61] proved such a homogenization result for Mumford-Shah type functionals [33] and in the context of periodic homogenization, i.e., where both the stiffness tensor field and the crack resistance field are periodic with a period which is sent to zero in the process of homogenization.The setting of the Mumford-Shah functional arises as a special case of the Francfort-Marigo model under anti-plane shear loading.Moreover, the result of Braides et al. [61] considers a single time step only, i.e., does not account for irreversibility.To state the result, consider the functional (1.1) under anti-plane shear loading with periodic material parameters C and γ of periodicity η.Braides et al. [61] showed that this functional -converges for η → 0 to the homogeneous functional with effective stiffness C eff and effective crack energy γ eff , a scalar function of the normal n to the crack surface S. In addition to this homogenization result, Braides et al. [61] established explicit formulas for both effective quantities.The effective stiffness may be computed via well-established homogenization methods of linear elasticity [62,63] by solving the elastic cell problem [64].The formula for the effective crack energy is given by the periodic, γ -weighted minimum cut through one period of the structure [65,66].Particularly surprising about this result is the fact that the homogenization formulas for the stiffness and the crack resistance decouple, i.e., are completely independent from each other.In particular, the stiffness tensor field does not influence the effective crack energy and vice versa.Interestingly, this decoupling is no longer valid if certain degenerations of the stiffness and the crack energy are permitted, leading to what is called high-contrast energies in the mathematical literature [67,68].
Over the past 20 years, several extensions to the homogenization result of Braides et al. were published.Giacomini-Ponsiglione [69] proved that a similar homogenization result holds for the Francfort-Marigo model under anti-plane shear loading specifically accounting for the irreversibility constraint.Recently, Cagnetti et al. [70] demonstrated that essentially the same conclusions can be drawn for the case of random materials as in the periodic setting, i.e., they showed that the -limit has the form (1.3) and that appropriate corrector formulae holds for both the effective stiffness and the effective crack energy.In a different direction, Friedrich et al. [71] lifted the restriction to anti-plane shear for the two-dimensional setting.
These analytical results were complemented by approaches to compute the effective crack energy.Schneider [72] introduced an FFT-based solver for the minimum cut problem which relied on a duality result by Strang [73].To deal with artifacts in the solution fields and the accompanying convergence problems, the authors [74] established a novel solver and discretization method resulting in superior performance both in accuracy as well as in computation time.Furthermore, extensions to anisotropic materials [75] were proposed.For two-dimensional matrix-inclusion composite, the problem of finding a minimum cut through the microstructure simplifies (at least for some cases) to finding shortest paths through the microstructure.Inspired by the pioneering work of Jeulin [76][77][78] Ernesti et al. [79] used fast marching methods [80,81] to compute the effective crack energy.Similarities to limit load analysis [82] were pointed out by Michel-Suquet [83] who investigated this problem using FFT-based methods.Furthermore, they proposed a semi-analytical solution for the effective crack energy of a two-phase laminate.

Contributions
The work at hand takes a closer look at interface properties in the context of the effective crack energy of a microstructure [72,75,83].Originally, the effective crack energy arises naturally from homogenization results in the periodic [61] and the stochastic setting [70].We revisit primal and dual formulations in section 2.1.These results account for the crack resistance in the bulk phases only.However, for a number of composite materials, like fiber-reinforced materials it is well known that the adhesion between the individual phases is crucial for the effective strength of the resulting material.
To take a closer look at interface properties, we follow the traditional route and introduce a small positive parameter δ and thicken the interface between the different phases to an interphase whose thickness is proportional to δ. Assigning a crack resistance to this interphase, we are interested in the limit of the effective crack energy as δ goes to zero.For a two-phase laminate, see Sect.2.2, we obtain an analytic expression for the effective flow set which differs from a pure bulk laminate provided the crack resistance of the introduced interphase is strictly lower than the bulk crack resistances.
Comparing this observation to weak elastic interfaces [84], we notice that a different scaling of the interphase properties is natural to obtain an appropriate weak interphase.Indeed, if we repeated the previous exercise with an elastic three-phase laminate, and set C 3 ≡ C int for a fixed stiffness C int assigned to the interphase, its influence would vanish in the limit δ → 0. Rather, a scaling C 3 ≡ δ C int is necessary to introduce weak interface effects [84], more precisely a spring-type weak interface.
Physically speaking, the interpretation of weak elastic interfaces and weak interfaces in the context of the effective crack resistance is thus completely different.For elastic interfaces to be weak, the elastic stiffness of the interphase needs to become infinitesimally soft, as present, for instance, in an adhesive material.For the crack resistance, in contrast, no such degeneration is necessary to result in interface effects.
However, both weak interfaces share the same interpretation: they arise not only from homogenization but in the limit as two different parameters, the correlation length in homogenization and the interface thickness related to the interphase property, tend to zero.
After identifying the need to consider a more general and robust definition of the effective crack energy, we provide such an extension in Sect.2.3 following the work of Baldi [85] on weighted functions of bounded variation (BV).We demonstrate that the extended definition encompasses the previously identified limit laminate naturally and provide a simplified approximation of the flow set which is convenient for numerics, see Sect.2.4.
These theoretical developments are complemented by an associated numerical technique.More precisely, we extend the combinatorial continuous maximum flow discretization (CCMF) [74], which we revisit in Sect.3.1, to account for weak interfaces in terms of the composite voxel technique [57,86], more precisely laminate composite voxels [58,59,87].The composite-voxel approach, originally formulated for voxel-based computational micromechanics, introduces specific constitutive laws in voxels which are intersected by the interface between bulk materials.Within these composite voxels, instead of the material response of a single material, the material behavior of a laminate with an appropriate normal and suitable volume fractions of the present materials is used to account for the heterogeneity of the material within the voxel.The original purpose of the composite voxel method was to enhance the accuracy of FFT-based methods which enables the use of a coarser grid and thus reduces the overall computational cost.
We mimic the composite-voxel strategy in the context of the effective crack energy by using appropriate laminates in the composite voxels.Moreover, we account for the weak interface by using the effective crack energy of a laminate with weak interface based on the previously developed (semi-)analytical description.The integration into the CCMF-framework is discussed in Sect.3.2.The CCMF-discretization is complemented by an efficient solver framework [75] based on the alternating direction method of multipliers (ADMM), see Sect.3.3.To handle composite voxels, it is necessary to project onto the effective flow set of laminates efficiently.Such a procedure is described in Sect.3.4.
We investigate the proposed methods for complex computational examples in Sect. 4. Notation: We denote the absolute value of a real number q with |q| = q 2 and write for the norm of a vector in R n .
2 The effective crack energy of a microstructure

The effective crack energy with bulk crack resistance
Suppose a cuboid domain ) is given, together with a periodic scalar field γ : Y → R of crack resistances, which we assume to be Lebesgue measurable and to obey the bounds for positive numbers γ ± independent of x ∈ Y .For a prescribed mean crack normal ξ ∈ S d−1 , the effective crack energy [72,74] in direction ξ is defined as follows where the infimum is taken over all functions p which are integrable together with their gradient.Here and in the following, we restrict to periodic functions and denote the corresponding function spaces with a #-symbol.
Intuitively, the fact that the integrand in the expression (2.2) grows linearly in the gradient of the field p fosters the minimizers to concentrate on a surface, i.e., a minimal surface weighted by the crack resistance γ .In particular, it appears evident that the minimum of the variational problem (2.2) is not attained in the class W 1,1  # (Y ).Rather, it is imperative to work with functions of bounded variation p ∈ BV # (Y ), i.e., Lebesgueintegrable functions p ∈ L 1 (Y ), s.t. the total variation is also a BV function, and we have the identity The functions of bounded variation form a Banach space under the natural norm Moreover, the total variation (2.3) is lower semicontinuous w.r.t.L 1 -convergence.A key property of BVfunctions is that their gradient defines a vector-valued Radon measure, i.e., they represent a rather general class of functions which admit a reasonable calculus.We refer to Evans-Gariepy [88, §5] for further details.The problem at hand (2.2) may be interpreted as a weighted version of the total variation in view of the equivalence (2.4).The term inside the infimum of the minimum cut problem (2.2) may be extended to all BV # (Y ) functions based on the identity ( valid for p ∈ W 1,1 # (Y ) and established via integration by parts.Leveraging the compactness and lower semicontinuity properties of the BV space, one obtains the representation Let us return to the BV-formulation of the minimum-cut problem (2.7), and discuss a dualization technique introduced by Strang [73] which is of vital importance to compute the effective crack energy (2.2) numerically.
We change the order of minimization and maximization, assuming that this operation is permitted, (2.9) The identity implies the representation of the effective crack energy.The problem (2.11) is called maximum-flow problem due to the divergence constraint characterizing incompressible flow fields.It maximizes the average flow under a point-wise constraint on the length of the flow vector.
In fact, Michel-Suquet [83] pointed out that the maximum-flow problem (2.11) may be conveniently recast in the form with the set of attained average flows To reveal the potential advantage of using the maximum-flow formulation (2.11), notice that the original minimum-cut problem (2.2) involves an integrand that is homogeneous of degree one.In particular, the functional is not differentiable at the origin.In contrast, the maximum-flow problem involves a linear objective function, a linear constraint encoding incompressibility and the point-wise norm constraint which may be recast in the form of quadratic constraints In particular, the formulation (2.11) is amenable to classical techniques of differentiable optimization [89,90].

The effective crack energy of a two-phase laminate with a weak interface
The purpose of this section is to investigate weak interfaces in microstructures for the example of a laminate, i.e., a material which is layered in a specific direction n and which contains a number K of different phases.
Michel-Suquet [83, Appendix A] provided an analytical expression for the effective flow set (2.13) of laminates.More precisely, they considered the case of two phases.However, their arguments carry over directly to the case of K phases with volume fractions φ i and crack resistances γ i .Notice that the unit normal vector n permits to decompose any vector v ∈ R d uniquely into n-parallel and n-transverse Then, the effective flow set of a K -phase laminate is given by Let us use this result to study a two-phase laminate with a weak interface, see Fig. 1b.To do so, we will actually study an interphase problem and take an appropriate limit.More precisely, suppose three crack resistances γ i (i = 1, 2, 3) are given, together with a volume fraction φ ∈ (0, 1) of the first phase.For the parameter δ ∈ (0, min(φ, 1 − φ)), we consider a three-phase laminate with volume fractions and corresponding crack resistances γ i (i = 1, 2, 3).The setup is illustrated in Fig. 1a.For any fixed value δ > 0, the effective flow set takes the form (2.18) Fig. 1 Illustration of two-and three-phase laminates with unit normal n In the limit δ → 0, we obtain the effective flow set In case the crack resistance γ 3 of the interphase exceeds the bulk crack resistances, i.e., the inequality In particular, the two sets do not coincide, and a distinct interface effect is noticable.Put differently, the limiting effective crack resistance takes the reduced interface crack resistance γ int , as shown in Fig. 1b, into account.

The effective crack energy with a weak interface
The previous section revealed that the notion of the effective crack energy as discussed in Sect.2.1 is insensitive to interface effects.The purpose of this section is to show that this defect is not intrinsic to the minimum cut/maximum flow problems but easily repaired by choosing the proper framework.Recall from Sect.2.1 that the minimum cut problem (2.2) could be interpreted as a weighted total variation problem.However, due to the appearance of the weight γ in the integral, interface crack resistances do not matter.Indeed, the Lebesgue integral of a function is impervious to values on a set of Lebesgue measure zero, e.g., on a (sufficiently regular) surface.
Baldi [85] studied functions of bounded variation w.r.t. a general positive weight function ω.More precisely, she considered a locally integrable weight function ω ∈ L 1 #, loc (Y ), belonging to the Muckenhoupt A 1 (Y )-class of weight functions [91,92], i.e., functions which are characterized by the inequality valid for almost every x ∈ Y and a universal constant c > 0. Here, B r (x) denotes the ball of radius r around x ∈ Y in the Y -periodic distance and the symbol − refers to the mean integral.For such a weight function ω, Baldi [85] defines the weighted total variation in complete analogy to the unweighted case (2.3).Please notice that the constraint of the flow field v in the definition (2.24) is pointwise everywhere, not pointwise almost everywhere.
The corresponding set of weighted BV functions BV # (Y ; ω) consists of (classes of) functions p : Y → R, which are integrable w.r.t. the ω-weighted Lebesgue measure, i.e., which satisfy the condition and whose weighted total variation (2.24) is finite.Similar to their unweighted counterparts, the weighted BV functions form a Banach space w.r.t. the norm s.t. the weighted total variation (2.24) is lower semicontinuous w.r.t.ω-weighted L 1 -convergence.Moreover, existence of ω-weighted minimal surfaces can be established.We refer to Baldi [85] for details.
Let us revisit the weighted minimum-cut problem (2.2).We would like to use the crack resistance γ as a suitable weight ω.Notice, for a start, that the upper and lower bounds (2.1) imply that the crack resistance γ defines an element of the Muckenhoupt A 1 (Y )-class of weight functions due to the estimate γ (y) dy. (2.27) Thus, we may consider the space BV # (Y ; γ ) to be a reasonable candidate to contain minimizers of the minimum cut problem.Then, we may consider the effective crack resistance where we used that the assumptions (2.1) imply that the weighted and unweighted BV -spaces actually coincide, for instance evident from the norm equivalence Let us compare the novel formulation (2.28) with the previous formulation (2.7) where we attached a subindex prev to be able to distinguish the two formulations.The difference between the defintions (2.28) and (2.30) is that the constraint v(x) ≤ γ (x) is enforced everywhere for the former and almost everywhere for the latter.In particular, the estimate follows by a mollification argument.
However, some care has to be taken.Baldi [85,Lemma 3.1] points out that one may assume without loss of generality that the weight function ω is lower semicontinuous, i.e., it satisfies the condition for any sequence (x k ) in Y converging to a point x ∈ Y as k → ∞.Indeed, for any weight function ω, consider the lower semicontinuous envelope  Let us put this observation into a physical perspective.Suppose we have a two-phase medium comprising two disjoint domains Y 1 and Y 2 whose union covers the cell, with associated constant values of crack resistance γ 1 = γ 2 .A priori, the crack resistance of the interface I = Y 1 ∩ Y 2 is not uniquely defined.However, the previous discussion entails that the only consistent extension on the interface satisfies the condition (2.38) Indeed, suppose the crack resistance γ would be higher in a neighborhood of x ∈ I .Then, one may displace a crack through the interface at x, represented by a surface, a little bit to either side of the interface and use the lower value min(γ 1 , γ 2 ).As this displacing shift can be as small as desired, such an elevated value γ (x) cannot be perceived upon relaxation.
If only the bulk values γ 1 and γ 2 are prescribed, relaxing the minimum-cut problem (2.28) via its lower semicontinuous envelope (2.33) automatically determines the crack resistance of the interface γ (x) = min(γ 1 , γ 2 ), x ∈ I. (2.39) Interestingly, the theory (2.28) developed by Baldi [85] permits the interface to have a lower, yet positive crack resistance.Such weak interfaces are of primary interest for the article at hand.With similar arguments as in Sect.2.1 for bulk crack resistances, we obtain the maximum flow representation of the effective crack energy.Also, the convenient form holds with the set of attained average flows where the overline refers to the set-theoretic closure.

A simplified approximation of a two-phase laminate with weak interface
The purpose of this section is twofold.On the one hand, we would like to show that the effective crack energy of a two-phase laminate with weak interface according to the new definition (2.28) coincides with the limit (2.19) of a three-phase laminate with vanishing interphase, as derived in Sect.2.2.On the other hand, we will derive a simplified approximation of the flow set (2.19) of a two-phase laminate with weak interface that is amenable to efficient numerical treatment.Suppose a two-phase laminate with a weak interface is given, see Fig. 1b.More precisely, we consider the first phase Y 1 to occupy a volume fraction of φ ∈ [0, 1], whereas the second phase Y 2 comprises the remaining volume 1 − φ.The direction of lamination n is given in terms of a unit vector in R d .
We assume that the two bulk phases are endowed with crack resistances γ i (i = 1, 2), whereas the interface I = Y 1 ∩ Y 2 is characterized by a crack resistance γ int .Thus, the local crack resistance field is given explicitly in the form The function γ is lower semicontinuous provided the interface crack resistance γ int does not exceed both bulk resistances, i.e., the inequality γ int ≤ min(γ 1 , γ 2 ) (2.44) holds.To determine the effective crack energy of such a laminate we follow Michel-Suquet [83] and determine the set of realized average flows (2.42) Then, the following representation may be deduced, which coincides with the limit (2.19) of the flow set of a three-phase laminate with vanishing interphase thickness.To maintain a fluid reading experience, the derivation is outsourced into Appendix A. With a more effective computational resolution in mind, we will also consider the following simplification.Instead of the admissible flow set (2.46), we consider the following simplification which is motivated by considering the case of a homogeneous laminate with weak interface, i.e., the case It is readily seen that the set inclusion volume averaging leads to the second inequality in the definition of the simplified admissible flow set (2.47).We investigated the difference between the exact laminate (2.46) and the approximation (2.47) numerically using the ECOS solver [93] for quadratic cone problems.An exemplary solution for the effective crack energy of a laminate with direction n = e x is shown as a polar plot in Fig. 2a.We chose the parameters φ = 0.5, γ 1 = γ , γ 2 = 10 γ and γ int = 0.5 γ .We observe that the exact solution is almost indistinguishable from the approximation.This is also quantified by the relative error, see Fig. 2b which is below 1% for the investigated material properties.We would also like to highlight that the error vanishes completely both in lamination direction n and perpendicular to this direction.

The combinatorial consistent maximum flow discretization
To solve the maximum flow problem (2.40) it is necessary to select an appropriate discretization scheme together with a suitable solution method.Due to the complexity of industrially used microstructures, discretizations on regular voxel grids [62,63] have been demonstrated to lead to powerful and computationally efficient schemes for a variety of physical problems [55, §5].However, it is not uncommon that such regular grid based discretizations lead to oscillations in the vicinity of material interfaces [55, §2].Thus, it does not come as a surprise that these traditional techniques, e.g., based on trigonometric polynomials [62,63], do not show optimal performance for the maximum flow problem (3.1) due to the appearance of the pointwise constraint on the flow field.The combinatorial consistent maximum flow (CCMF) discretization [74,94] was introduced as a remedy, building upon the advantages of finite-volume discretizations.More precisely, let us consider a voxel discretization Y N with N i (i = 1, 2, 3) voxels in each coordinate direction and cubic voxels with the edge length h, i.e., the identity h = L i /N i holds for i = 1, 2, 3. Notice that we restrict to the physically relevant dimension d = 3, but typically use d = 2 for visualization.
In the classical setting of scalar bulk crack resistances, we suppose that a crack resistance γ is given for each voxel, see Fig. 3a.Then, the CCMF discretization associates flow degrees of freedom to the faces of the voxels, i.e., each voxel has six (d = 3) resp.four (d = 2) adjacent flow variables to account for.Notice that the flow variables and the crack resistances live on different grid points, see Fig. 3a for two and Ernesti-Schneider [74,Fig. 4] for three dimensions.These flow variables are not independent, as each face is shared by two adjacent voxels.For this reason, we attach the flow variables to a specific voxel.This procedure is not unique, and we use the association for each voxel with integer coordinates [i, j, k].Then, we may express the incompressibility constraint by the equation with the divergence operator where we consider the indices to by (N 1 , N 2 , N 3 )-periodic.The discrete incompressibility condition (3.3) ensures the balance of in-and outflow for each voxel, as usual for finite-volume methods.
To discretize the maximum flow problem (3.1), it remains to discretize the pointwise constraint Recall that the crack resistance γ is associated to each voxel, whereas the flow field lives on the voxel faces, see Fig. 3a.Therefore, we use the operator A, defined for each voxel [i, j, k] by

.6)
The operator A collects the flows on all six faces of the voxel and multiplies them with a normalization factor [74].With this operator at hand, the CCMF discretization replaces the constraint (3.5) by the voxel-wise condition In view of the definition (3.6), the latter condition may also be written in explicit form The original constraint (3.5) features three terms in the norm on the left hand side, the CCMF discretization involves six terms on the left hand side of the constraint (3.8).To ensure consistency for homogeneous flows, a factor two needs to be included into the right hand side.To retain a formal similarity of the CCMF-constraint (3.7) with the non-discretized equivalent (3.5), the prefactor 1/ √ 2 is included in the definition (3.6) of the A-operator.
With these operators at hand, we may write down the discrete maximum flow problem in the form with the constraint set (3.10)

Combinatorial consistent maximum flow and composite voxels
To account for weak interfaces (2.40) in the CCMF discretization, we take a second look at the maximum flow problem (3.9) by writing the constraint set (3.10) in the form with the flow set associated to the voxel with index [i, j, k].In this formulation (3.11)-(3.12),more general, in particular nonisotropic, flow sets V[i, j, k] may be considered naturally for computing the effective crack energy.Such an idea was already exploited for locally anisotropic crack resistance [75].To account for weak interfaces, we make use of the composite voxel method [57,86], which was originally introduced to improve the accuracy of FFT-based computational micromechanics on regular grids.Let us consider a two-phase microstructure comprising a matrix material with embedded inclusions, see Fig. 3b.For FFT-based computational micromechanics, the microstructure is typically encoded by a color field given on this voxel image.Each color uniquely identifies the phase, and the color of each voxel determines to Fig. 4 Subvoxeling approach and substitute laminate which phase the corresponding voxel belongs to.For the two-phase microstructure shown in Fig. 3b, color 1 corresponds to the matrix and color 2 encodes the inclusions.
For the composite voxel method, one assumes that a voxel image with a much finer resolution is given.Based on this representation, a third color is introduced for voxels which are not homogeneous in this fine-scale representation, see Fig. 4a.
The classical composite voxel method in mechanics [57][58][59] assumes the interface between the materials to be accurately described by a linear surface within the voxel and furnishes the composite voxel with the effective mechanical properties of an equivalent laminate, see Fig. 4b.More precisely, based on the finely resolved voxel image, estimates for the volume fractions of the phases and for the normal of the interface between the phases is computed.Following Kabel et al. [57,59], we evaluate the volume fractions of the individual phases via quadrature on the fine grid.The unit normals are computed by normalizing the vector pointing from the centroid of one phase to the center of the composite voxel.
Thus, we end up with three phases: the sets Y N 1 and Y N 2 of matrix and inclusion voxels and the composite voxel phase Y N  3 , where for each voxel with indices (i, j, k) the volume fraction φ[i, j, k] of phase 1 and the unit normal n[i, j, k] are available.
For use in the maximum flow problem (3.11), we furnish any composite voxel (i, j, k) ∈ Y N 3 with the flow set (3.14) based on the approximation (2.47) of the equivalent laminate (2.46) and where we defined n CCMF ∈ R 6 by Together with the definitions for the isotropic phases, we may formulate the maximum flow problem (3.11) & (3.12) with weak interface by defining the local flow sets (3.17) Notice that sets V 1 and V 2 both correspond to spheres of radius γ 1 and γ 2 .The set V 3 [i, j, k] describes an anisotropic shape of a sphere with radius φ[i, j, k] γ 1 + (1 − φ[i, j, k]) γ 2 and spherical caps removed, see Fig. 5a.This contrasts with the investigation [75] where ellipsoidal flow sets were used to model an anisotropic crack resistance.

The ADMM solver for the maximum flow problem
In this section, we recall the alternating direction method of multipliers (ADMM) solver tailored to solving the maximum flow problem (3.11) & (3.12) (3.18) The corresponding discrete minimum cut problem is derived similarly to both the isotropic case [74,Sec. 3.1] and the ellipsoidal case [75, Sec.3.1], i.e., we obtain the problem with the voxel-wise (1-homogeneously extended) crack resistance and the set K ξ of compatible fields with the discrete forward gradient operator ∇ + and the left inverse A * of A, i.e., characterized by the condition A * A = Id.We solve the discrete maximum flow problem using the alternating direction method of multipliers (ADMM) [95,96], first introduced into FFT-based methods by Michel et al. [97,98].To derive the algorithm, we first define the two convex functions (3.26) Following the same section, line 1 and 3 may expressed explicitly using the orthogonal projectors P K ξ and P C onto the sets K ξ and C N (where we suppress the dependence on the discretization parameter N ).
The projection P K ξ onto the compatibility set is efficiently computed using FFT [74, eq. ( 33)], whereas the projection onto the constraint set P C may be evaluated for each voxel independently by projecting onto the corresponding flow set.

The laminate projection problem
To use ADMM for solving the minimum cut/maximum flow problem, evaluating orthogonal projectors (3.28) onto the sets V c (c = 1, 2, 3) in equation (3.17) is required.In the isotropic case, the set V c describes a sphere with radius γ c (c=1,2), and the projection proceeds via rescaling the norm if necessary.
For a fixed composite voxel (i, j, k), i.e., for the set the following projection problem needs to be solved.Given a vector w ∈ R 6 , we seek a vector v ∈ R 6 solving the problem For sake of readability, we set γ = φ γ 1 + (1 − φ) γ 2 , suppress the dependence on the voxel (i, j, k) for both the volume fraction and the normal, and omit the superscript CCMF for the normal, as well.This projection problem may be interpreted as follows.The first constraint describes a sphere with radius γ .The second constraint introduces two half-planes, described by v • n ≤ γ int and v • n ≥ −γ int .Hence, the set we project onto is a sphere where two spherical caps were removed, see Fig. 5a.To project onto this set, see Fig. 5a, which has a piecewise smooth surface, we distinguish three cases and follow a geometric argument.First, we decompose the vector w in a normal part w = (n • w) n and a tangential part w ⊥ = w − w .Based on the three vectors w, w and w ⊥ we describe three cases which characterize the domain outside of V 3 , see Fig. 5b.
• Domain D I is described via the conditions w > γ int and w ⊥ < γ 2 − γ 2 int .(3.30)In this case, only the normal part of w has to be reduced.• Domain D II is given by the constraint w > γ int and and is dependent on the relation between w ⊥ and w .• Domain D III is defined via the constraints Based on the decomposition the projection onto the set V 3 is readily achieved by the projectors P I , P II and P III valid for their respective domains ) In domain I, only the normal part of w is scaled, whereas the tangential part remains the same.In domain II and the displayed quarter of the planar domain in Fig. 5b, all points are projected onto the same point.Thus, the projector P V 3 reads We implemented the governing equations (3.27) into an in-house FFT-based micromechanics solver written in Python with Cython extensions and we parallelized the code using OpenMP.We used the CCMF discretization and the damped alternating direction method of multipliers with an adaptive parameter strategy [74].We used a damping factor of δ = 0.25 and the averaging adaptivity approach proposed by Lorenz-Tran-Dinh [99].We solved the governing equations up to a tolerance of 10 −4 , measured via the convergence criterion [74, eq. ( 36)].The simulations were performed on a desktop computer with 32GB RAM and 6 cores of each 3.7GHz.As a first step we investigate multi-grid convergence of our approach on a two-dimensional structure containing a rotated square shown in Fig. 6a.The square inclusion is located in the center of the microstructure, has a edge length of L/2, where L denotes the length of the microstructure, and is rotated by 45 degrees.The discretized composite voxel structures of N × N pixels were generated via downsizing a (16N ) × (16N ) image, see Fig. 6b for the discretized structure with N = 128.We consider a crack resistance γ for the matrix material, 3 γ for the inclusion and γ int = 0.5 γ in the interface.The resulting minimum cut for ξ = e y for N = 128 is shown in Fig. 6c.Due to the material contrast between the interface and both the inclusion and the matrix, it is energetically more favorable for the cut to follow the interface.Due to the symmetry of the structure, two cuts are equally likely.The effective crack energy may be computed exactly from the geometry e l s e (4.1) In the first case the cut follows the interface of the inclusion whereas for γ int > γ / √ 2 the cut is a straight line passing the inclusion.Using the exact solution, we may compute the relative error for each resolution and study multi-grid convergence.We compare our composite voxel approach with a naive method of inserting voxels at the interface with isotropic crack resistance γ int , i.e., assuming V 3 = {v ∈ R 6 | v ≤ γ int } within the composite voxels.A comparison of the error is shown in Fig. 6d.We see that even for a low resolution of 16 2 pixels the composite voxel approach yields an error below 1%.Furthermore, multi-grid convergence is observed upon refinement and the accuracy of the composite voxel approach is one magnitude higher than the naive approach of inserting single voxels with a lower crack resistance.Even for the finest resolution of 128 2 pixels the naive approach does not reach the accuracy of the composite voxel approach of the coarsest grid of 16 2 pixels.

A polycrystalline microstructure with weak interfaces
Next, we consider a polycrystalline brittle material with weak interfaces between the grains.The individual grains are considered to be isotropic with crack resistance γ .The interface is considered weaker with crack resistance γ int , modeled via composite voxels.Notice that in this case, where the crack resistance on each side of the interface is equal, the simplified laminate problem is equal to the original laminate problem and thus the affiliated projector computes the exact laminate projection.The microstructure under consideration contains 15 grains, see Fig. 7 and was generated as Laguerre tesselations [100] with a uniform grain size distribution, i.e., each grain takes up 6.67% of the total volume.The composite voxel image of resolution 256 3  was generated via downsizing a 1024 3 microstructure.We consider two different interface crack resistances, Fig. 8 Polycrystalline microstructure and minimum cut for different interface crack resistances see Tab. 1 for the material parameters and the resulting effective crack energies for ξ = e x .Figure 8a shows the minimum cut for γ int = 0.5 γ .We notice two sections of this minimum cut.In some section it follows the interface between the grains, whereas in others it cuts a grain in a straight line.This contrasts with Fig. 8b which depicts the minimum cut for γ int = 0.25 γ .In this case the cut is entirely restricted to the interface.In the first case the effective crack energy compared to the crack resistance of the grains is reduced by 12%, whereas the second case shows a reduction of 45%.

A fiber reinforced composite
As our final example, we consider the microstructure of a fiber reinforced composite shown in Fig. 9a which was generated using sequential addition and migration [101,102].All fibers are of equal length and diameter with an aspect ratio of 10 and a total volume fraction of 30%.The fiber orientation tensor of second order [103] was chosen as diag(0.6,0.2, 0.1).In order to accurately resolve the surface of the fibers, we chose a resolution of 20 voxels per fiber diameter in the composite voxel image, see Fig. 9b, resulting in a 256 3 structure.The composite voxel image was generated by downsizing a 1024 3 microstructure.We set the crack resistance of the matrix material to γ and the crack resistance of the fibers to 10 γ .Furthermore, we vary the interface crack resistance ranging from 10 −3 γ to γ .The effective crack energy in the favored fiber direction for varying interface crack resistance is shown in Fig. 9c.We notice a reduction of the effective crack energy as the interface crack resistance decreases.If the interface is 90% weaker than the matrix material, the effective crack energy will be reduced by 19%.However, even for an almost vanishing interface crack resistance, the observed effective crack energy is higher than the matrix crack resistance.In other words, a toughening of the material due to cylindrical inclusions is observed even in the presence of a weak interface.
Figure 9d-f show the cut surfaces with mean normal e x for γ int = γ , γ int = 0.5 γ and γ int = 0.1 γ .We notice that for a lower interface crack resistance the overall surface is less smooth and adapts closer to the interface.In particular, for a lower interface crack resistance, the effect of fiber pullout becomes more prominent.

Conclusion
In this work we studied the effective crack energy of microstructured materials with a weak interface between the constituents.As a first step, we established the solution for minimum cut/maximum flow for a laminate with a weak interface.Due to a challenging numerical treatment of the resulting model, we proposed a simplified Fig. 9 Fiber microstructure, influence of the interface crack resistance and minimum cut model that returns similar results.Taking this simplified model, we followed Kabel et al. [57,59] who introduced laminate composite voxels to FFT-based micromechanics.We used their approach, combined with our laminate formula for the effective crack energy, to equip the interfaces in complex microstructures with a lower crack resistance and computed the effective crack energy of the composite.We investigated multi-grid convergence in comparison to a naive approach, i.e., setting the crack resistance in interface voxels to a lower value.In a consecutive study we investigated polycrystalline materials with isotropic grains of a higher crack resistance than the interface.We noticed that above a certain threshold the resulting minimum cut is entirely restricted to the interface.Last but not least, we computed the effective crack energy of a complex fiber reinforced composite material with weak interfaces.We investigated the influence of the crack resistance of the interface and found that even for an almost vanishing crack resistance, a toughening of the effective crack energy in contrast to the matrix material may be observed.
Our approach establishes the use for composite voxels in FFT-base micromechanics for computing the effective crack energy.Originally introduced to ensure high accuracy in computing effective properties while reducing the overall computational cost, we proposed a way to use them to implement a different material behavior at the interface, similar to Chen et al. [56] for phase-field fracture.In a consecutive step, we may combine the proposed approach with a minimum cut/maximum flow formula for locally anisotropic crack resistances [75].This may require a laminate formula with an anisotropic crack resistance for the bulk phases of the laminate.

Fig. 3
Fig. 3 Illustration of the CCMF discretization and a two-phase microstructure

Fig. 5
Fig. 5 Domain V 3 and orthogonal projectors w), w ∈ D I , P II (w), w ∈ D I I , P III (w), w ∈ D I I I , w,

Fig. 6
Fig. 6 Rotated square microstructure, minimum cut and error versus resolution

Table 1 Fig. 7
Fig. 7 A polycrystalline microstructure With this remark at hand, one observes that the two formulations (2.28) and (2.33) may both be represented in Baldi's formulation.Indeed, starting from a given crack resistance field γ , the interface-aware formulation (2.28) arises by considering the weight ω = γ 34) holds for any p ∈ BV # (Y ).* , (2.35) involving the lower semicontinuous envelope (2.33).In contrast, the original approach (2.30) emerges by setting ω = Mγ (2.36) in terms of the Hardy-Littlewood maximal function Mγ (x) = sup γ (y) dy, x ∈ Y. (2.37) Indeed, for any weight γ in the Muckenhoupt A 1 (Y )-class, the function Mγ is lower semicontinuous, an element of A 1 (Y ) and coincides with the function γ almost everywhere, see Baldi [85, §1] for a discussion.