An FFT‐based method for computing weighted minimal surfaces in microstructures with applications to the computational homogenization of brittle fracture

Cell formulae for the effective crack resistance of a heterogeneous medium obeying Francfort‐Marigo's formulation of linear elastic fracture mechanics have been proved recently, both in the context of periodic and stochastic homogenization. This work proposes a numerical strategy for computing the effective, possibly anisotropic, crack resistance of voxelized microstructures using the fast Fourier transform (FFT). Based on Strang's continuous minimum cut—maximum flow duality, we explore a primal‐dual hybrid gradient method for computing the effective crack resistance, which may be readily integrated into an existing FFT‐based code for homogenizing thermal conductivity. We close with demonstrative numerical experiments.


INTRODUCTION
Multiscale methods in mechanics use information about the mechanical behavior of the individual phases and their spatial arrangement on a particular length scale to draw conclusions about the mechanical behavior on a coarser length scale. For strain-hardening material behavior, multiscale methods have been established as a powerful and viable tool, see the recent survey of Matouš et al. 1 In the context of homogenization methods, continuum mechanical computations are performed either on a periodic cell 2 or a representative volume element. 3 Applying multiscale methods to fracture mechanics, however, is not so straightforward. Classically, linear elastic fracture mechanics, see Gross and Seelig 4 for an elementary exposition, concerns the propagation of preexisting cracks. However, in the context of homogenization, the stress singularity at the crack tip precludes scale separation. It is possible to conduct direct numerical simulations on volume elements, as done in the context of cohesive zone models, 5 local and nonlocal damage models [6][7][8][9][10] or phase-field fracture, [11][12][13] and to subsequently identify the parameters of a postulated macroscopic model. However, upon strain localization, the results may depend strongly on the chosen cell 14 and the imposed boundary conditions. [15][16][17] This contrasts with homogenization of strain-hardening materials, where a representative volume is characterized by the effective stress being independent of the imposed boundary conditions. 18 Even more to the point, Gitman et al 14 demonstrate by numerical simulations that the critical stress of their investigated continuum damage model approaches zero in the infinite volume limit.
To overcome these limitations, several approaches have been developed, including coupled volume methods, 19 second-order computational homogenization, 20 coarse-graining of failure modes, [21][22][23][24][25] and failure-zone averaging. [26][27][28] In contrast, for the Francfort-Marigo model, 29 which generalizes Griffith's energy criterion 30 and serves as continuum limit of some phase-field fracture models, [31][32][33] Braides et al 34 provide a homogenization result and corresponding cell formulae for heterogeneous elastic tensor and crack resistance in the context of periodic homogenization and antiplane shear. More precisely, they provide a formula for the effective stiffness tensor and a formula for the effective crack resistance (as a function of the unit normal), both of which are decoupled. We will recapitulate these formulae in Section 2.1.
This decoupling of elasticity and crack resistance may be startling at first. However, it is a result of the different scaling of the summands in the Francfort-Marigo model. The elastic energy is a volumetric quantity, whereas the crack resistance is multiplied by an area. For other models, such a decoupling through homogenization is well-known. For instance, in the context of fully coupled thermomechanics, the mechanical problem and the heat equation decouple almost entirely upon homogenization. 35 In a recent work, Cagnetti et al 36 extended the homogenization result and the cell formulae to stochastic homogenization. In particular, the existence of the representative volume element is ensured for homogenizing the crack resistance.
Apparently unaware of the homogenization results 34,36 for brittle fracture, Hossain et al 37 introduced a notion of effective fracture toughness and conducted computational studies by a phase-field fracture model. [31][32][33] We discuss the differences between the two approaches in Section 2.1.

Contributions
In this work, we propose a numerical strategy for computing the effective crack resistance of Braides et al 34 and Cagnetti et al 36 on large-scale digital volume images based on the fast Fourier transform (FFT). For prescribed average cracksurface normal, we introduce a convenient cell formula for computing the effective crack resistance which permits arbitrarily periodic crack surfaces. The proposed formula is related to the original cell formula of Braides et al 34 in the same way that imposing periodic boundary conditions is related to imposing displacement boundary conditions for mechanical problems. 18 In particular, in the infinite volume limit, both approaches are expected to coincide.
Based upon Strang's work, 38 we recast the weighted minimal surface problem central to the cell formula as a convex optimization problem (the minimum cut formulation). The latter program, being convex, has the advantage that any critical point is a global minimizer. Again building upon Strang's insights, 38 we exploit dual and primal-dual formulations of the minimumcut problem, leading to the maximumflow problem. The details comprise Section 2.2.
Based upon the primal-dual formulation, we propose a numerical algorithm for computing the effective crack resistance. More precisely, we investigate the primal-dual hybrid gradient (PDHG) method, 39,40 which has been successfully applied to image segmentation problems 41 , and point out similarities to polarization-based methods 42-44 popular in FFT-based computational micromechanics, see Section 3 for details.
Due to formal similarities to the homogenization of thermal conductivity problems, we investigate two popular discretization schemes used in FFT-based methods: the under-integrated Fourier-Galerkin discretization, pioneered by Moulinec and Suquet 45,46 and the rotated staggered grid discretization, 47 which may also be interpreted as an under-integrated finiteelement discretization. 48 The approach in this article differs from other FFT-based computational approaches to phase-field crack propagation, 11-13 as we do not compute stress-strain curves for applied mechanical loading. Instead, the effective, possibly anisotropic crack resistance is our quantity of interest. The latter may be used in a phase-field crack model 49,50 on the macroscale.
We demonstrate the power of the proposed numerical schemes for problems of industrial complexity, cf. Section 4.

Homogenization of brittle fracture
In the context of quasi-static small-strain brittle fracture mechanics, let Ω be a domain in R d with sufficiently smooth boundary, and let a linear elastic tensor C ∈ L(Sym(d)) and a positive crack resistance * be given. Herein, Sym(d) denotes *In linear elastic fracture mechanics, the terms fracture toughness, crack resistance, and critical energy release rate are used synonymously. In this work, we consistently use "crack resistance". the vector space of symmetric d × d matrices. For prescribed boundary displacements u 0 ∶ Ω × [t 0 , t 1 ] → R d , depending on (pseudo-)time, the Francfort-Marigo model 29 seeks a displacement field u ∶ Ω × [t 0 , t 1 ] → R d , such that, for any t ∈ [t 0 , t 1 ], u(⋅, t) minimizes the energy where S u denotes the set of jump points of u (away from the boundary Ω) and ∇ s stands for the symmetrized gradient, subjected to prescribed boundary displacements u(⋅, t) = u 0 (⋅, t) on Ω and the constraint that cracks may only grow, and an initial crack set S u(⋅,t 0 ) . To render the model Equation (1) well-defined for each instant of time, the displacement field u needs to be chosen in a suitable function space of discontinuous functions, see Chambolle et al. 51,52 To treat the Francfort-Marigo model numerically, it is customary to discretize the quasi-static evolution 53 in time by a backward Euler scheme, that is, to subdivide the time interval [t 0 , t 1 ] into K+1 increasing time instants t 0 = t 0 , t 1 , t 2 , … , t K−1 , t K = t 1 and to minimize, for each k = 1, 2, … , K, the functional (1) subjected to the boundary condition u k = u 0 (⋅, t k ) on Ω and the constraint S u k ⊇ S u k−1 , where u k and u k−1 denote minimizers of the functional (1) corresponding to the time steps t k and t k−1 , respectively.
For a (periodically) microstructured material, let a rectangular unit cell Y ⊆ R d be given, together with a heterogeneous elastic tensor field C ∶ Y → L(Sym(d)) and a lower semicontinuous crack resistance ∶ Y → R + , both of which we consider as Y -periodically extended functions on R d . For > 0 and a given time discretization, we consider the heterogeneous problem subjected to similar boundary conditions and constraints as above, see the left-hand side of Figure 1 for a schematic. For ≪ 1, the functional W in Equation (2) is very difficult to solve numerically (in addition to the "usual" computational challenges associated to the Francfort-Marigo model, see Bourdin et al 31 ). Thus, for ≪ 1, it is convenient to approximate Equation (2) by a suitable homogenized model, as is common for other types of multiscale problems. 1 In the context of antiplane shear, Braides et al 34 have shown a periodic homogenization result for problems of the type (2). To be more precise, a special case of their result concerns a given rectangular unit cell Y ⊆ R d , together with periodic, possibly discontinuous, coercive elastic tensor field C ∶ Y → L(R d×d ), and lower semicontinuous crack resistance ∶ Y → R + , bounded from above and from below. Then, under suitable conditions on C as well as on and as → 0, the functionalsW Ω Ω Y converge, in the sense of Γ-convergence, to the homogenized functional for homogeneous, effective "stiffness" C eff and crack resistance eff ∶ S d−1 → R, a (symmetric) function of the unit vector n normal to the crack surface. Furthermore, cell formulae for both C eff (the formula appearing in classical works like in Babuška 2 ) and eff are provided. More precisely, the effective crack resistance eff ∶ S d−1 → R, is given by and, for L > 0, eff where Q L (n) is the cube ] d rotated in such a way that the e 1 -axis is mapped onto n, and the infimum is taken over all surfaces S lying inside Q L , whose boundary is constrained to Q L ∩ {x ∈ R d |x ⋅ n = 0}. Geometrically, Equation (6) seeks the -weighted minimal surface with fixed boundary conditions on the halved cube Q L . The latter definition is visualized in Figure 2A, where the case d = 2 is considered. The figure shows nine copies of the unit cell, with attaining two distinct values 1 and 2 . Q L is centered at the origin, such that one of the face normals corresponds to n. Furthermore, we see two competitors for the "shortest line," shown in green and blue, both with identical, fixed boundary conditions. The blue line is shortest if 2 < 1 , whereas the green line is (almost) optimal if 2 ≫ 1 . As the basis for this work, we assume that the homogenization result of Braides et al 34 continues to hold beyond the restrictive case of antiplane shear, that is, we assume that the limiting model, as → 0 and in the sense of Γ-convergence, of the multiscale problem (2), is given by where the effective stiffness C eff is determined by the usual formula 54 the infimum is taken over all smooth, periodic fields u ∶ Y → R d (|Y| denotes the volume of the cell Y ), and the effective crack resistance eff ∶ S d−1 → R may be computed as detailed in Equation (5).
To conclude this section, several comments are in order.
F I G U R E 2 Schematics for eff L (n), cf. Equation (6), and for determining the effective fracture toughness following Hossain et al 37 1. From a mechanical point of view, the result of Braides et al 34 may be surprising, as the "elastic" and the fracture part decouple during homogenization, that is, C eff solely depends on C and eff can be determined from alone. On second thought, the result becomes clearer taking into account the different scaling of the two summands entering the Francfort-Marigo energy (1). Indeed, for more general energies of Mumford-Shah type, these two summands may interact. For instance, Braides et al 55 consider the relaxation of energies of the form where the surface part is 1-homogeneous in the displacement jump [u] = u + − u − across the interface. Due to this 1-homogeneity, the bulk and surface parts have the same scaling, and may interact upon relaxation. Also, for functionals of the type (3) with high-contrast coefficients, interaction upon homogenization can be proved. More precisely, Barchiesi et al 56  Equation (5), that was established by Braides et al 34 and Cagnetti et al. 36 Hossain et al consider mode-I-crack propagation in a plane-stress setting, cf. Figure 2B. More precisely, they consider a square cell Ω centered at the origin, with a preexisting crack located at the negative x-axis, and filled by the periodic material. They investigate the Francfort-Marigo model (1) subjected to "surfing" displacement boundary condition, that is, those which solve mode-I-crack propagation in a homogeneous isotropic setting, and determine the time-dependent energy release rate by evaluating the J-integral 59,60 on the boundary Ω. As a result, they obtain a function J eff , depending on (pseudo-)time, and define eff HHBB ∶= max This choice is motivated by the interpretation of the crack resistance as the critical energy releaserate, that is, as the minimum value of the energy release-rate which permits unhindered crack growth. With this definition at hand, Hossain et al 37 show that heterogeneity in the elastic moduli as well as in the crack resistance may lead to an increase in the crack resistance. These conclusions appear to contradict the homogenization theory of Braides et al 34 and Cagnetti et al. 36 These differences in predictions originate from the different scales of the boundary-condition increments Δu 0 for a discretization in time. Hossain et al 37 consider increments which are only a fraction of the microstructural length. Otherwise, crack propagation would not be able to stagnate between the microstructural phases. In the setup of Equation (2), the increments Δu 0 are macroscopic, that is, the crack growth between successive time steps may traverse a large number of unit cells.
Notice that the setup of Hossain et al, 37 cf. Figure 2B, constitutes a special case of the heterogeneous Francfort-Marigo model (2), cf. Figure 1. In particular, upon → 0, the homogenized functional (7) is recovered, together with the corresponding cell formulae. Put differently, eff arises from a suitable averaging procedure of J eff as → 0.
From a physical point of view, the toughening effect of heterogeneous media predicted by Hossain et al 37 may be observed experimentally if is small, but not exceedingly so, that is, for microtensile experiments 61 if the length scale of the heterogeneities is small. In contrast, for purely macroscopic loading conditions, the homogenization theory of Braides et al 34

A fully periodic cell formula and the maximum flow formulation
To work with the homogenized model (7) the effective stiffness C eff and the effective crack resistance eff ∶ S d−1 → R need to be determined, either analytically or by computational means. The cell formula (8) for C eff is the usual one, and classical techniques may be applied, see Matouš et al 1 for a recent overview.
In contrast, the cell formulas (5) and (6), and n ∈ S d−1 , are less standard, essentially for three reasons. First, the minimal surface problem with prescribed boundary condition is nonconvex in its original formulation. We will see below that the minimal surface area may be computed using convex functionals, however. Second, in contrast to the cell formula (8), the definition (5) of eff (n) involves a limiting procedure on cells of growing size. Such formulae are well-known in mathematical, asymptotic, and computational homogenization of random materials or in the nonconvex setting. As eff may be computed by minimizing a (sequence of) convex functional(s), the level of difficulty is comparable to the computational homogenization of random strain-hardening materials. The third difficulty in treating the cell formulas (5) and (6) is the intrinsic degeneracy of the functional. Indeed, there may be several (or even a continuum of) minimizers of the functional in Equation (6), and the computational method needs to account for that.
Last but not least, by engineering experience it appears to be smarter to use periodic boundary conditions for the "cutting surface," although we do not provide a mathematical proof that our formulation is equivalent to the original formulation. With these considerations in mind, we fix the cell Q L ≡ Q and assume that ∶ Q → R is defined on Q per se, that is, we do not care whether arises from periodization or not, and seek, for prescribed , the infimum of the variational problem Herein, plays the role of the unit normal (typically normal to one of the faces of Q) and |Q| denotes the volume of Q. It can be shown, cf. Strang, 38 that Equation (9) gives rise to a well-defined problem on the space of (periodic) functions of bounded variation, and the interfaces of sublevel sets {x ∈ Y| (x) < t} of the minimizer , for almost any t ∈ (0, 1), are periodic -weighted minimal surfaces. Actually, Strang only considers the case of continuous . However, Hintermüller et al 62 provide a corresponding result for crack resistances with jumps using a Fenchel predualization technique.
From a computational point of view, the minimization problem (9) appears somewhat innocent. If the norm term in the integrand of Equation (10) was raised to a power of two, we would face a homogenization problem that typically arises, for instance, in thermal conductivity. However, the exponent one gives rise to some difficulties. Although, the existence of minimizers for Equation (10) is easy to establish † , the solution may not be unique. For instance, if is homogeneous, any cut across Q normal to gives rise to a minimizer. Also, the functional in Equation (10) is not classically differentiable. Still, Equation (10) involves a convex function, that is, on the one hand we may use subdifferential calculus to characterize critical points, and, on the other hand, know that any critical point is global minimzer. Thus, all we need to find is some critical point. Naive subgradient algorithms are comparatively slow, cf. Section 4.3 in Chambolle and Pock. 41 Thus, following similar approaches for image processing problems, we investigate dual and primal-dual formulations of Equation (10).
For fixed ∈ R d , introduce the (convex) set  of kinematic constraints and denote by ⟨⋅⟩ Q the averaging operator on Q. Then, for a convex function f , the formal dual of the problem is given by see Hintermüller et al, 62 where f * denotes the Legendre-Fenchel dual of f , that is, The primal minimizationproblem (10) may be put in the abstract primal form (11) by setting f( ) = ⟨ || ||⟩ Q . Then, the Legendre-Fenchel dual function according to Equation (13) computes as the indicator function  of the convex set , Thus, the dual formulation of the problem (10) seeks a vector field v ∶ Q → R d , such that, holds. Put differently, the formal dual of the minimum cut problem (10) is the maximumflow problem (15), which seeks to maximize the flow ⟨v⟩ Q in direction subject to an incompressibility constraint and a bound constraint on the flow v.
The strong duality of these problems, that is, was first shown by Strang 38 for continuous and extended to piece-wise continuous by Hintermüller et al. 62 This minimum cut-maximum flow duality immediately gives rise to simple bounds on the effective crack resistance. Indeed, inserting = into the primal problem, and v = (min ) into the dual problem provides the two-sided bound Last but not least, we will discuss briefly a primal-dual formulation of the minimal cut problem (16). As f is convex, its bidual f * * equals f , and we may write f as the Legendre-Fenchel dual of its Legendre-Fenchel dual Thus, we may write

A PDHG METHOD FOR COMPUTING THE EFFECTIVE CRACK RESISTANCE
To find a saddle point of the primal-dual formulation (18) Arrow et al 63 proposed to alternate a (proximal) descent in the variable and an ascent in the dual variable v with respective step sizes s and t where P  and P  denote the L 2 -orthogonal projections onto the closed convex sets  and , respectively. The iterative scheme (20) is not convergent; in general. Pock et al 40 and, independently, Esser et al 39 noticed that, for st ≤ 1, a simple overrelaxation in one of the variables, that is, in our context, improves the convergence behavior of this primal-dual hybrid gradient method (PDHG). More precisely, Chambolle and Pock 64 showed an ergodic 1 k -convergence rate for the PDHG method (21) if the step sizes satisfy st ≤ 1. In the context of FFT-based methods, the Hilbert space under consideration, both for and v, is the space L 2 (Q; R d ), equipped with the L 2 inner product Then, for any ∈ L 2 (Q; R d ), the projector P  becomes in terms of the Green's operator, a linear operator Γ ∶ L 2 (Q; R d ) → L 2 (Q; R d ) given by Γ = ∇Δ −1 div, see for instance, Vondřejc et al. 65 Also, P  may be evaluated explicitly on any ∈ L 2 (Q; R d ) via || (x)|| (x), otherwise.
Applied to the PDHG iteration (21), we rewrite the first line in the form Introducing the "reference material" 0 = 1 s , the associated Green's operator Γ 0 = 1 s Γ, and rewriting st ≤ 1 in the form t = 0 for ∈ (0, 1] the PDHG iteration becomes where evaluating P  is specified in Equation (23). Thus, the PDHG method alternates a single step of the basic scheme in thermal conductivity, cf. Eyre and Milton, 42 and a local projection operation. Several remarks are in order.

←5
: v ←6 : old ← 7: while k < maxit and residual> tol do 8: k ← k + 1 9: ← v − 0 ‖v(x)‖ , otherwise. 16: residual old ← 18: end while 19: return , v, residual 1. 0 has the same dimensions as the crack resistance , whereas the factor has dimension 1. In the absence of theoretical hints for choosing 0 , we investigate several rational choices of 0 in terms of the function ∶ Q → R in Section 4. 2. We initialize the algorithm by 0 = and v 0 = . 3. As the convergence criterion, we choose that is, we measure div v in an appropriate norm, cf. Schneider et al 44 for elaboration in the context of small-strain problems. 4. For = 1, the PDHG iterative scheme (24) may be equivalently rewritten as a Douglas-Rachford method (with relaxation factor 1 2 ), see Section 5.4 in Chambolle and Pock. 41 In the context of FFT-based computational homogenization methods, the Douglas-Rachford scheme (with relaxation factor 1 2 ) was introduced by Michel et al 43 where ( * , v * ) denotes a saddle point, that is, a fixed point of the iterative scheme (24). Some caution has to be taken in interpreting estimate (26), because will only be a bounded Radon measure, in general. Thus, its L 2 norm will not be finite. However, we are mostly interested in this estimate for the discretized problems. In that case, the L 2 norm will be finite, but not uniformly so. Rather, the norm will blow up as the mesh size gets refined.

Implementation
The PDHG method detailed in Algorithm 1 was implemented into an in-house FFT-based computational homogenization code, based upon previous work on the computational homogenization of thermal conductivity. 66 Python with Cython extensions was used as the programming language of choice, relying upon OpenMP for parallelization. The examples were run on a laptop with two 2.70 GHz cores and 16 GB RAM as well as a desktop computer with six 3.70 GHz cores and 32 GB RAM, respectively. We consider two classical discretization schemes used in FFT-based computational homogenization: 1. The discretization by trigonometric polynomials of Moulinec and Suquet, 45,46 where, in addition, the "energy" is approximated by the trapezoidal rule, cf. Vondřjec et al. 65 2. The rotated staggered grid of Willot et al, 47 which may be interpreted as a discretization by Q1-(trilinear hexahedral)-finite elements with reduced integration. 48 Formulae for the Γ-operator of these schemes are readily available in the mentioned references. For both discretizations, we set the Nyquist frequencies of the Γ-operator to zero, that is, we force the Nyquist frequencies of to be zero (We could have alternatively chosen to force the Nyquist frequencies of v to zero).
In general, we set the parameter entering Algorithm 1 equal to 0.99. We use tol = 10 −4 in the convergence criterion Equation (25), unless mentioned otherwise.

Validation and calibration
As a warm-up, we consider simple examples with analytically known effective crack resistances. We use these examples to study the accuracy of the discretizations, and to determine the sensitivity to the algorithmic parameters.

A two-phase laminate
First, we investigate a simple two-phase laminate with volume fractions 1 = 1 3 and 2 = 2 3 , with lamination direction e x , and phase-wise constant crack resistances 1 and 2 , respectively. Then, we have eff (e x ) = min( 1 , 2 ) and eff (e y ) = ⟨ ⟩ Q = 1 The microstructure was discretized by 3 3 = 27 voxels. We test the PDHG solver for different contrasts 2 ∕ 1 and different choices of the reference material, that is, 3. Arithmetic mean : 0 = − + + 2 , 4. Maximum: 0 = + , where − = min( ) and + = max( ). Due to the initialization, the case = e y is solved in a single iteration. Thus, we only investigate = e x . Both the Moulinec-Suquet and the rotated staggered grid discretization give rise to identical iteration counts (for both solvers). Thus, the number of iterations is collected in a single Table 1. We see that choosing the minimum for the reference material is best. Furthermore, increasing the contrast also increases the iteration count. However, the iteration count is not symmetric with respect to multiplicative inversion.

A single-ball inclusion
Next, we investigate a single spherical inclusion in a cubic box with crack resistance 1 . The inclusion has a volume fraction of 12.89%, and is equipped with a crack resistance of 2 . For studying the convergence behavior of the algorithms for different discretizations and reference-material choices, we work on a 32 3 microstructure, and solve up to a tolerance of 10 −4 . Iteration counts for the Moulinec-Suquet discretization are collected in Table 2. Consistently, a contrast of 100 required more iterations than a contrast of 10, which in turn, was more expensive than a 2 ∕ 1 -ratio of 0.1. For the minimum and the geometric mean reference-material choice strategy and a contrast of 0.01, less iterations were required than for a contrast of 0.1. For the other two investigated choices, the reverse was true. Still, the "minimum" choice outperformed all other reference material choices by a large margin. For Willot's discretization, cf. Table 2, similar conclusions as for Moulinec-Suquet's discretization hold. The number of iterations differs slightly, but is roughly on the same level.
Last but not least, let us compare the effective crack resistances computed for different discretizations, cf. Table 2. All obtained values coincide for the first three significant digits.
A local comparison of the solution fields is shown in Figure 3. The crack resistance 2 of the fibers is chosen as 2 = 10 1 , where 1 denotes the crack resistance of the matrix. For = e y , it is optimal to crack linearly through the matrix without touching the inclusion. For the Moulinec-Suquet discretization and the discretization on a rotated staggered grid, the fields are not homogeneous, as might be expected for the analytical solution. Instead, they exhibit ringing and checkerboard artifacts, as usual for these types of discretizations, cf. Dorn and Schneider. 66 This may also be observed for the flow fields v, cf. Figure 3A. However, these visual artifacts have no influence on the predictive quality of the computed effective quantities.  Next, we consider a resolution study for the single-ball microstructure. The same geometry was remeshed, using 32 3 , 64 3 , 128 3 , and 256 3 voxels. Using the PDHG iterative scheme, and 0 = min , the iteration counts are listed in Table 3. As a general trend, we see that the iteration count increases with increasing resolution. This contrasts with other iterative algorithms used in FFT-based computational micromechanics, such as the basic scheme 45,46 or the conjugate gradient method. 67 However, the mentioned algorithms are applied to less degenerate problems, and the logarithmic convergence rate Equation (26) of the PDHG scheme implies that the iteration count is very sensitive to the initial condition. Still, the total iteration counts are on a reasonable level for practical applications.
Turning our attention back to Table 3, we see that the iteration counts of the Moulinec-Suquet discretization and the discretization on the rotated staggered grid are on a similar level. For 2 ∕ 1 = 10 and 2 ∕ 1 = 100, the iteration counts are comparable for fixed voxel count. For 2 ∕ 1 < 1, the iteration counts are lower, on average. However, the iteration counts for 2 ∕ 1 = 0.01 remain below 400, whereas for 2 ∕ 1 = 0.1 they are much larger, with a peak for 64 3 voxels.
In general, if 2 ∕ 1 ≫ 1, the inclusions act as obstacles for the minimal surfaces, and the algorithm needs to account for these hindrances. In contrast, for 2 ∕ 1 ≪ 1, traversing the inclusions is "free," simplifying the problem.
Last but not least, we investigate the computed effective crack resistances, listed in Table 4. The discretizations induce only differences in the fourth significant digit for fixed resolution. Also, the resolution dependence is negligible.

Continuously fiber-reinforced composites
In this section, we investigate the effective transverse crack resistance of continuously fiber-reinforced composites. We consider unidirectional microstructures with up to 50% fiber-volume content, which were generated by the adaptive shrinking-cell algorithm of Torquato and Jiao 68 modified to not shear the cell. We fix the crack resistance matrix of the matrix and set fiber = 10 matrix for the fibers.
First, we investigate the size of a representative volume element, 3,14 and subsequently discuss the effective transverse crack resistance as a function of the fiber-volume fraction.

Representativity study
Due to the transverse isotropy of continuously fiber-reinforced composite microstructures, the effective transverse crack resistance has to be isotropic. The rectangular shape of the computational cell Y and the finiteness of the number of fibers breaks this symmetry and induces artifacts we wish to quantify in this section. We fix the x-axis to coincide with the fibers' principal axis and consider an average crack normal in the y-z-plane, that is, perpendicular to the fibers. We consider different angles for the crack normal, that is, = cos e y + sin e z for angles ranging from 0 o to 90 o in 15 o -steps. The procedure is shown in Figure 4 for a reference microstructure with 50% fiber-volume fraction, 100 fibers and a resolution of 256 2 . We see that, for varying angle , the crack paths are changing. Consistently for all angles, several crack paths with roughly identical lengths emerge. We wish to determine the unit-cell size for which the computed effective transverse crack resistance is representative. Keeping the investigation of the fiber-volume fraction influence of Section 4.3.2 in mind, we consider a fiber-volume fraction = 50%, as the statistical fluctuation should be smaller for smaller fiber-volume fractions, as well.
Based on the reference cell with 100 = 10 2 fibers and a resolution of 256 2 pixels, we successively increase both the fiber count and the resolution, maintaining the accuracy level of the results. More precisely, we consider (10k) 2 F I G U R E 5 Differently sized unit cells for the representative volume element study for 50% filler content and = e y , including (10k) 2 fibers discretized by (256k) 2 pixels. We show || ||, clipped to [0, 10], with axes and scale as in Fig. 4  fibers resolved by (256k) 2 pixels and k = 1, 2, 3, 4. To get an impression on the cell sizes, we refer to Figure 5, where the results for = e y are shown. For increasing volume-element size, the number of pure "matrix islands" increases. Otherwise, the structures appear statistically similar. For the four considered microstructures, the angle-dependent effective crack resistances are listed in Table 5. To arrive at an effective scalar quantity, taking into account the expected isotropy, we computed the mean (and the standard deviation) of the results, as well. We see that the standard deviation decreases roughly by a factor of two if the edge length of the microstructure is doubled. If only the absolute value eff is of interest, the standard deviations are exceedingly small. Indeed, the standard deviation for the smallest cell is already below 1%. However, if the relative increase eff − matrix serves as the reference quantity, the standard deviation is at about 26% for the smallest microstructure. For the largest microstructure, the standard deviation is about 6%. This dependence is visualized in Figure 6A, where we also see that the effective transverse crack resistance is not

Dependence on the volume fraction
In the previous section, a representative volume-element size with 20 2 fibers was determined, and we wish to investigate the dependence of the effective transverse crack resistance on the fiber-volume fraction . For this purpose, volume elements with different fiber-volume fractions were generated, increasing in 5% steps up to 50% volume fraction. We fixed the fiber count to 400 and the resolution to 512 2 pixels. The resulting unit cells are shown in Figure 7, together with the crack paths for = e y . As a general trend, the number of possible cracks decreases for increasing fiber-volume fraction. Indeed, for low filler content, a variety of (almost) straight cracks is possible, whereas, for high filler content, due to the close packing, only few shortest paths exist. With the same procedure as in the previous section, the mean and standard deviation for the considered crack-normal angles were computed. The results comprise Figure 6B. The effective transverse crack resistances increase with the fiber-volume fraction. Also, the standard deviation increases. To get an increase by 1% in transverse crack resistance, about 30% fibers (in volume) are necessary, whereas a 2% increase necessitates slightly more than 40% fibers. Overall, the increase in transverse crack resistance is not excessive.

Industrial-scale examples
In this section, we consider more realistic microstructures.

A porous microstructure
A sand-binder aggregate, a typical microstructure of a blown sand core used for foundry applications, is shown in Figure 8A, discretized by 256 3 voxels. The microstructure was generated by the mechanical-contraction method described in Schneider et al, 69 with a sand-volume fraction of 58.6%, held together by 1.28% anorganic binder. These microstructures are modeled based on typical sieve lines of casting sands and microcomputed tomography images. For our investigation, we fix the crack resistance of the sand and the binder phase to a value , and set the crack resistance of the pore space to 0.01 . For the described discretization methods, the PDHG method was used to solve the saddle point problem (18) with = e y up to a precision of 10 −4 . The iteration counts and effective crack resistances are listed in Table 6. For both the Moulinec-Suquet and the rotated staggered grid discretization, the iteration counts are similar, approximately 2350. Also, the computed effective crack resistances are similar for the Moulinec-Suquet and the rotated staggered grid discretization.
To gain a visual impression, we take a look at Figure 8, where the resulting crack paths are shown. More precisely, we thresholded the norm of at a level of 10 to produce the images. As expected, on average, the crack paths are normal to F I G U R E 7 Influence of increasing volume fraction for 400 fibers, 512 2 pixels and = e y . We show || ||, clipped to [0, 10], with axes and scale as in Figure 4A [Colour figure can be viewed at wileyonlinelibrary.com] the prescribed average normal = e y . For the considered discretizations, cracking takes place almost exclusively in the binder phase. This might be a result of the geometric necks induced by the binder. From a mechanical point of view, this appears plausible, as stress concentrations are expected to appear in these necks. The predicted transverse crack paths are identical for both considered discretizations.

A short-fiber reinforced composite
As our next example, we consider a short-fiber reinforced plastic featuring a brittle matrix. The structure we consider is composed of 18% short fibers with a length of 275 m and a diameter of 10 m, dispersed in a cubic box with an edge

TA B L E 7
Results for the PDHG method with 0 = min and the fiber-filled microstructure shown in Figure 9A Moulinec-Suquet Rotated staggered grid length of 512 m. The fiber orientation tensor is chosen as perfectly isotropic, and the structure is discretized by 256 3 and 512 3 voxels, corresponding to voxel resolutions of 2 m and 1 m, respectively. The microstructure was generated by the Sequential Addition and Migration method, 70 and contains 1119 fibers. We set the crack resistance of the fibers to 10 times the crack resistance of the matrix, which we denote by matrix . The PDHG method was used to solve the saddle point problem (18) to a precision of 3 × 10 −4 .
The effective crack resistances computed using the two discretization methods are listed in Table 7. For a resolution of 256 3 , the predicted effective crack resistances differ in the third significant digit. Thus, we have rerun the computations on 512 3 voxels. For this resolution, the computed effective crack resistances agree to three significant digits for both discretization schemes, and predict an increase by about 62% in the effective crack resistance compared with the pure matrix. In contrast to the transverse crack resistance of a continuously fiber-reinforced composite, the influence of the reinforcing fibers clearly manifests for this isotropic short-fiber reinforced example.
The iteration counts are listed in Table 7, as well. For 256 3 , Moulinec-Suquet's discretization requires about 2000 iterations, which is about 50% of the iterations required by the rotated staggered grid. In contrast, for 512 3 , the iteration count for the Moulinec-Suquet discretization is roughly tripled, whereas, for the discretization on a rotated staggered grid, the iteration count is almost the same as for the low resolution.
In Figure 9, both the original microstructure and the computed crack field is shown for the Moulinec-Suquet discretization and a resolution of 256 3 . For the rotated staggered grid, the predicted crack field is similar.

CONCLUSION
This work was devoted to studying FFT-based solvers for computing the effective crack resistance of a heterogeneous microstructure, according to a variant of the cell formula of Braides et al. 34 Based on ideas of Strang, 38 we relied upon a convex reformulation of the cell problem for determining the effective crack resistance. For the computational resolution, we proposed a PDHG method in an FFT-based context. The PDHG is striking in its simplicity, and may be readily integrated into an existing FFT-based thermal homogenization code. For Moulinec-Suquet's discretization and the discretization on a rotated staggered grid, we demonstrated the applicability of the introduced solver on problems of academic interest and also of industrial scale.
A simple finite-volume type discretization 66 was also investigated, but led to strong mesh artifacts, whence we decided to not report on this discretization. These artifacts may be overcome by employing more sophisticated cut metrics. 71,72 We leave this as future work.
For the industrial-size problems we investigated, the PDHG combined with the rotated staggered-grid discretization performed best, and may-in the presented form-prove extremely valuable for applications to questions of materials science.
The biggest limitation we see is the low convergence rate of the PDHG algorithm. However, extensive numerical experiments with other solvers did not improve upon the simple-to-implement PDHG. This might be intrinsic, that is, as a result of the nonsmoothness and degeneracy of the functional under consideration, or significantly improved by dedicated solvers or formulations. For instance, connections to the limit-load problem 73 could be exploited.