Computing the effective crack energy of heterogeneous and anisotropic microstructures via anisotropic minimal surfaces

A variety of materials, such as polycrystalline ceramics or carbon fiber reinforced polymers, show a pronounced anisotropy in their local crack resistance. We introduce an FFT-based method to compute the effective crack energy of heterogeneous, locally anisotropic materials. Recent theoretical works ensure the existence of representative volume elements for fracture mechanics described by the Francfort–Marigo model. Based on these formulae, FFT-based algorithms for computing the effective crack energy of random heterogeneous media were proposed, and subsequently improved in terms of discretization and solution methods. In this work, we propose a maximum-flow solver for computing the effective crack energy of heterogeneous materials with local anisotropy in the material parameters. We apply this method to polycrystalline ceramics with an intergranular weak plane and fiber structures with transversely isotropic crack resistance.


Introduction
Griffith [1] formulated an energetic crack criterion for preexisting cracks in the context of an isotropic, linear elastic body, laying the foundation of modern fracture mechanics [2]. He noted that cracks propagate under an incrementally increased load provided the newly formed crack surface is energetically more favorable than further elastic deformation. The identified material parameter indicating a critical energy is called critical energy release rate.
A further extension of linear elastic fracture mechanics was proposed by Irwin [3] who computed analytical solutions for a linear elastic pre-cracked body under certain loading scenarios, which he called modes. He identified stress intensity factors governing the behavior of the stress field at the crack tip and defined a crack propagation criterion based on the critical values of the stress intensity factors, referred to as fracture toughness. Note that in case of linear elastic, isotropic and homogeneous materials, Griffith's energy criterion may be expressed in terms of stress intensity factors [4], The mentioned approaches focus on the question when cracks propagate. To asses how cracks propagate, the concept of maximum energy-release [5] and the principle of local symmetry [6] are classical.
Traditionally, fracture mechanics focused on cracks in isotropic and linear elastic materials. Early extensions to elastoplastic fracture were proposed by Dugdale [7] and Barenblatt [8]. In case of anisotropic linear elasticity, Sih et al. [9] showed that the classical theory may result in complex-valued stress intensity factors. As no analytical solution is available for heterogeneous materials, their treatment requires different concepts.
The absence of analytical solutions may be remedied by computing the local stress fields numerically [10]. Since classical finite element discretizations have difficulties in resolving the stress at the crack tip accurately due to the stress singularity, extended finite elements [11] were introduced, which add ansatz functions with jumps in the displacement field. An alternative approach includes cohesive zone elements [12] where distinct traction-separation laws [13] may be utilized.
In 1998, Francfort and Marigo [14] reformulated Griffith's theory as a variational problem, introducing the Francfort-Marigo model of fracture. For a fixed domain and a fixed time discretization, they seek minimizers, i.e., the displacement field u and the crack surface S across which u may be discontinuous, of the energy functional The energy consists of a volume-energy part accounting for elastic deformations, and a surface-energy part quantifying the crack-surface energy. Physical assumptions of their model include S t ⊆ S t+1 at all time steps t, i.e., the crack surface may only grow. Notice that their formulation is based on global minimization for reasons of mathematical well-posedness, whereas Griffith seeks local critical points. The Francfort-Marigo model naturally includes heterogeneous material properties, as both the crack resistance γ and the stiffness tensor C may depend on the position. Furthermore, distinct material anisotropy may already be expressed in terms of an anisotropic stiffness tensor. Additionally, as noted by the authors themselves [14, eq. (17)] the surface energy may account for anisotropy by replacing the isotropic crack resistance in the surface energy by an anisotropic term γ (x, n(x)) depending on the unit normal n of the crack surface. The prevailing tool for treating the Francfort-Marigo model computationally is the phase-field model of fracture. The method, introduced by Bourdin, Francfort and Marigo [15], is motivated by the Ambriosio-Tortorelli approximation [16] of the Mumford-Shah functional [17], used in image segmentation. The phase-field model involves a length-scale parameter η and seeks minimizers u as well as d, namely the displacement field and the damage variable, of the functional Chambolle [18] showed that, for η → 0, the phase-field model −converges to the Francfort-Marigo model. Additionally, the phase-field model shows similarities to nonlocal damage models [19,20] and may be treated as such as long as η > 0 is regarded as a material parameter [21]. The phase-field fracture model is rather popular in the scientific community, including a variety of contributions over the last decade, see Wu et al. [22] for a recent overview. Of particular interest to the work at hand are contributions that account for material anisotropy. Investigations on modeling anisotropic fracture using an anisotropic stiffness but isotropic crack energy were carried out [23,24]. Van Dijk et al. [25] suggested an approach to incorpo-rate tension-compression anisotropy into the context of an anisotropic stiffness. Clayton-Knap [26] introduced a geometrically nonlinear phase-field model with an anisotropic fracture toughness, which takes, in case of small deformation elasticity, the form and a field of unit vectors p. They applied their model at small deformations to simulating cleavage in polycrystalline ceramics [27]. The positive cleavage parameter β is introduced to penalize crack propagation within the plane perpendicular to the unit vector p. Further extensions to incorporate crystal plasticity and ductile fracture were reported [28,29].
The tensor M p permits to model either one weak plane, or one tough direction. Introducing a general symmetric and positive definite tensor M, up to three different crack resistances, i.e, the eigenvalues of M, in the three eigendirections may be prescribed. More general approaches were proposed using a multi phase-field setting, see Nguyen et al. [30], or a higher order phase-field method, using fourth order tensors [31][32][33][34], to study polycrystalline materials. Pillai et al. [35] proposed an anisotropic cohesive phase-field model to simulate the behavior of anisotropic fiber structures. Incorporating weak interfaces via cohesive elements was proposed by Rezaei et al. [36].
To account for the influence of the microstructure of a material to its macroscopic behavior, multiscale methods may be used, see Matouš et al. [37] for an overview. For hardening material behavior, those multiscale approaches are well understood. Softening materials, in contrast, where distinct strain localisation may occur, are more challenging [38].
Classically, homogenization relies on a distinct scale separation. More precisely, the macroscopic displacement or stress fields should vary slowly compared to the local fields on the microscale. In the classical linear elastic fracture mechanics setting and for an evolving crack, the stress singularity at the crack tip typically prohibits such a scale separation. Indeed, in classical linear elastic homogenization, the displacement field is split into a macroscopic, smoothly varying part, and a microscopic, highly oscillating part. Upon homogenization, these two contributions decouple up to a rather simple, one-way coupling, which permits to transfer information from the microscale to the macroscale (via the effective stiffness). In the case of an evolving crack, the displacement field is no longer smooth at the crack tip, even for a homogeneous medium. In particular, the classical macromicro decomposition which was successful for linear elastic homogenization loses its promise. Moreover, the evolution of the crack needs to be accounted for at both scales.
In case of a fixed discretization in time (more precisely pseudotime, quasi-static problems are concerned), Braides et al. [39] established a periodic homogenization result for the Francfort-Marigo model under anti-plane shear. They showed that, for a distinct scale separation in the material parameters, the Francfort-Marigo model homogenizes to the effective functional with effective, possibly anisotropic stiffness tensor C eff and effective crack energy γ eff (n). Furthermore, they showed that the two effective quantities decouple, i.e., the local stiffness tensor has no influence on the effective crack energy and the local crack resistance does not influence the effective stiffness. A key ingredient for the homogenization result of Braides and coworkers was a formulation on a single unknown, namely the displacement field which is permitted to be discontinuous across specific crack surfaces. In this way, the headache concerning the distinction of a displacement field and the crack can be avoided. Indeed, the (possibly jumping) displacement field can be decomposed additively into a macroscopic and a microscopic part. The work of Braides et al. was further extended to stochastic homogenization by Cagnetti et al. [40], ensuring representative volume elements to exist for the Francfort-Marigo fracture model under anti-plane shear. Recently, the restriction to anti-plane shear was lifted by Friedrich et al. [41], so the homogenization result holds in general. To compute the effective quantities, Braides et al. [39] provide specific formulas. Computing the effective stiffness reduces to classical linear elastic homogenization [42], whereas the effective crack energy γ eff (n) associated to a crack with normal n may be computed as a γ -weighted minimal surface cutting the microstructure. Note that the homogenization results are based on the notion of −convergence. Thus, they only concern global minimizers of both the original and the effective functional.
To avoid confusion, we call the effective quantity γ eff (n) effective crack energy instead of effective crack resistance, as the latter term is used ambiguously in the literature. Hossain et al. [43] define the effective crack resistance as the maximum J-integral [44] evaluated along propagating cracks computed by phase-field fracture on heterogeneous structures. A different approach was pursued by Lebihain et al. [45], who investigated cracks bypassing inclusions, using a perturbative approach. They define the effective crack resistance as the maximum elastic energy release rate.
These two approaches differ from the ansatz proposed by Braides et al. [39] in the treatment of the scales. Both Hossain et al. and Lebihain et al. consider microstructures of a fixed size and compute crack propagation in continuous time, which is only discretized to enable a numerical treatment. Braides et al., on the other hand, fix a time discretization for the macroscopic model before evaluating the limit of the models as the microstructural length goes to zero. As a result, the microstructures considered in case of Braides are much smaller compared to Lebihain or Hossain. In case of Braides et al., the crack, which propagates incrementally on the macroscale, will cross the microstructure within one time step.
Minimizing the fracture energy has been considered even before the mathematical results by Jeulin [46][47][48], who studied crack propagation in two-dimensional micrographs. In fact, in two spatial dimensions, the problem of computing the effective crack energy simplifies drastically, as it reduces to the problem of computing minimum (weighted) geodesics, where efficient algorithms are available [49,50]. Schneider [51] proposed a method for computing the effective crack energy for heterogeneous, locally isotropic microstructures using a reformulation of the cell formula of Braides et al. [39] into a convex optimization problem. The transformation of the cell formula relies on Strang's [52] minimum cut-maximum flow duality. The established optimization problem was solved using FFT based algorithms. Recently, Schneider & Ernesti [53] proposed a discretization on a combinatorially consistent grid, which improves the solution quality, introducing an associated adaptive ADMM solver.

Contributions
This work extends the approach of the authors [51,53] to account for a locally anisotropic crack resistance, enabling to treat matrix-inclusion composites with anisotropic phases or polycrystalline materials. In Sect. 2.1, we introduce the cell formula for the effective anisotropic crack energy and anisotropic phases. We describe the anisotropy in terms of a tensor, similar to approaches applied in phase-field fracture [26]. We transform the anisotropic minimum cut problem into an anisotropic maximum flow problem in Sect. 2.2, which gives rise to a convex optimization problem.
Powerful solution methods for maximum flow problems arose for maximum-flow problems on graphs [54]. Due to metrication artifacts, these are not directly applicable to the continuous maximum flow problem at hand [55]. This may be overcome by a combinatorial continuous maximum flow (CCMF) discretization [56], recently applied to computing the effective crack energy of solids [53]. We propose a way to account for anisotropic crack resistance within the CCMF discretization in Sect. 3.1. As the anisotropy of the crack resistance is described in terms of a (symmetric positive definite) second-order tensor, the maximum-flow formulation involves the projection onto an ellipsoid. We express the governing projection operator as the solution of a constrained optimization problem in Sect. 3.3. Finally, we apply our anisotropic minimum cut-maximum flow approach to brittle materials with a distinct anisotropy. In Sect. 4.2 we investigate polycrystalline brittle materials, such as ceramics with distinct cleavage in each grain. Last but not least, we investigate fracture of carbon fiber reinforced composites, where each fiber itself shows strong anisotropy of transversely isotropic kind.
2 Minimum cut-maximum flow with anisotropic crack resistance

Cell formula for an anisotropic minimum cut
and a given field of heterogeneous and direction-dependent crack resistances γ : Y × R 3 → R. There are several expressions for determining the effective material properties governing brittle fracture [43,45]. We follow the result established Braides et al. [39] and recent extensions [40,41]. In particular, the heterogeneous crack resistance field homogenizes to an effective crack energy γ eff (ξ), which depends on the macroscopic crack normalξ and is computed via a γ -weighted minimal surface with mean surface normalξ . Suppose that a field γ : Y × R 3 → R of anisotropic crack resistances is given. We assume that, for any microscopic point x ∈ Y , the association defines a norm on the vector space R 3 , i.e., it is onehomogeneous, satisfies the triangle inequality and is nondegenerate. Moreover, we assume that there are positive constants γ − , γ + , s.t. the inequalities Under these assumptions, we define the effective crack energy as a function on the unit sphere S 2 via where the infimum is evaluated over all smooth periodic fields p. Upon one-homogeneous extension, the effective crack energy γ eff gives rise to a norm on R 3 , as well.
In this article, we specialize the form of the microscopic crack resistances considered to those which describe an anisotropic Euclidean norm. More precisely, for a field G : Y → Sym(3) of symmetric, positive definite crack resistance tensors, we consider the local crack resistance field Then, Eq. (2.2) becomes 3) The isotropic case [51,53], with isotropic crack resistance field γ : Y → R, is recovered via G(x) = γ (x) Id.

Anisotropic maximum flow formulation
We define the energy functional and, for fixedξ ∈ R 3 , the set of kinematic constraints With the energy functional and the kinematic constraints at hand, we call the anisotropic minimum cut problem. For fixed normalξ , the minimum effective crack energy (2.3) computes as the minimum value of this minimization problem. Treating this problem numerically is challenging, since the energy functional is homogeneous of degree one and thus non-differentiable. Extending the isotropic case [51], we introduce a dual formulation, i.e., a corresponding anisotropic maximum flow problem [52]. The formal dual problem to the minimization problem (2.5) is given by where f * denotes the Legendre-Fenchel dual of f , given by and the minimum (2.6) is evaluated over all periodic solenoidal fields v. Since the tensor G is symmetric and Thus, the Legendre-Fenchel dual f * equals the indicator function +∞ otherwise, of the closed set Since G is symmetric and positive definite, the set C is a convex domain. More precisely, the constraint in the definition of the set C describes an ellipsoid centered at the origin, see Fig. 1. Combining (2.7) with this expression for C, we arrive at the anisotropic maximum flow problem

Numerical treatment
Ernesti and Schneider [53] proposed to solve the maximum flow problem in the combinatorial continuous maximum flow (CCMF) discretization [56] by FFT-based methods [57]. The formulation is based on doubling the degrees of freedom.
We present an extension of this strategy to account for an anisotropic crack resistance expressed via a heterogeneous field of crack resistance tensors G : Y → Sym(3). Thus we refer to Ernesti and Schneider [53] as a general reference throughout this section

The CCMF discretization for anisotropic maximum flow
We discretize the unit cell Y = [0, L 1 ]×[0, L 2 ]×[0, L 3 ] on a regular cubic voxel grid Y N with N i , (i = 1, 2, 3) voxels in each direction under the assumption that each voxel is cubic with edge length h = L i /N i , (i = 1, 2, 3). We evaluate the crack resistance tensor at the voxel center, i.e., This placement of the flow field enables to encode the conservation of mass via the discrete backwards divergence To enforce the constraint G −1 v ≤ 1 in the context of the CCMF discretization, we introduce the nonlocal shift operator S The constraint is then enforced via the inequality To integrate the CCMF discretization of the anisotropic maximum flow problem into an FFT-based homogenization solver for heat conductivity, we introduce the extension operator A, given by .
With this notation at hand, accompanied by the inner product v, w = 1 we may rewrite the discrete maximum flow problem as For the CCMF discretization, the set C from Eq. (2.8) reads The derivation of the associated discrete primal problem follows the same steps as the isotropic case, described in Ernesti and Schneider [53, Sec. 3.1]. The discrete minimum cut problem with a tensorial crack resistance is given by with set of discretely compatible fields The operator A * denotes the left inverse of A, i.e., A * A = Id holds, and ∇ + refers to the discrete forward gradient operator.

The alternating direction method of multipliers
To solve Eq. (3.6) with the alternating direction method of multipliers (ADMM), we rewrite the discrete minimum cut problem equivalently as with the two convex functions f (ξ ) = ι Kξ (ξ ) and where ι Kξ is the indicator function of the set Kξ described in (3.7). The operator-splitting approach starts by rewriting the problem as a constrained optimization problem The alternating direction method of multipliers (ADMM) [58,59] was introduced into the context of FFT-based methods by Michel et al. [60,61], see also the application to nonsmooth optimization by Willot [62]. Applied to our problem, we investigate the augmented Lagrangian function involving a penalty factor ρ > 0 and the Lagrange multiplier, i.e., our flow field, v : Y N → R 6 . The damped ADMM recursion [53,Sec. 3.2] with damping factor δ is given by (3.11) More explicitly, the ADMM algorithm with adaptive penalty factor ρ k is given by (3.12) where P Kξ and P C G denote the orthogonal projectors onto the sets Kξ and C G , respectively.

The anisotropic projector problem
Within the ADMM iterations (3.12), evaluating two projection operators is required. The projection onto the compatible fields may be efficiently performed with the help of the FFT [53,Eq. (3.16)]. Additionally, the orthogonal projection onto the set C G is required, expressed via the projection operator P C G and illustrated in Fig. 2. In the isotropic case, the set of constraints C comprises a sphere per voxel. Thus, the projection onto C involves orthogonal projections onto spheres with radii γ (x) and is straightforward. In the anisotropic case, however, the the set C G comprises an ellipsoid for each voxel, with the eigenvalues γ i of G as semi-axes. Following Kiseliov [63], we write this projection as an optimization problem. Consider a vector v ∈ R n and a symmetric, positive definite tensor Q. We seek the projection w ∈ R n , such that (3.13) Introducing the Lagrange parameter λ, the governing KKTconditions, see for instance [64,Thm. 12.1], read From conditions (3.14) we find If the vector v satisfies v T Q v ≤ 1, the problem (3.13) is trivially solved for w = v (and λ = 0). We therefore focus on the case v T Q v > 0. Since the tensor Q is symmetric and positive definite, the inequality (3.15) describes a convex domain. Thus, the projection w lies on the boundary of the admissible set. Thus, the inequality (3.15) is satisfied as an equality. We insert the representation (3.17) into the conditions w T Q w − 1 = 0 and solve the resulting equation for the scalar λ, using Newton's method and initial value λ 0 = 0, following Kiseliov [63]. In the mentioned reference, global convergence of this algorithm is proved.
For the application at hand, where n = 6, we set for each x ∈ Y . Therefore, the projection operator summarizes as with λ solving (3.18).

Setup
The presented computational approach was integrated into an existing FFT-based code for computing effective crack energies on microstructures [53], which is embedded into an FFT-based computational homogenization code for thermal conductivity [65]. The software is written in Python with Cython extensions (and OpenMP). All computations were performed using the ADMM solver presented in Ernesti-Schneider [53] with either the Barzilai-Borwein adaptivity or the averaging adaptivity, and a damping factor of 0.25. If not mentioned otherwise, the governing equations were solved for a prescribed tolerance tol = 10 −4 and the convergence criterion [66] (4.1) The computational experiments were run on a desktop computer with 32GB RAM and six 3.7GHz cores and on a workstation with 512 GB RAM and two Intel Xeon(R) Gold 6146 processors (12 × 3.20 GHz), respectively.

A polycrystalline microstructure
As our first example, we consider a polycrystalline microstructure generated synthetically based on Laguerre tessellations and the algorithm described in Kuhn et al. [67]. Following a similar approach as Nguyen et al. [30] within the context of phase-field fracture, we distinguish between 2D and 3D structures. In a 2D microstructure with distinct anisotropy, we identify one weak and one tough direction, whereas in 3D, we identify one weak and two tough directions, resulting in one weak plane. Since elastic deformation and thus elastic material constants do not play a role in our model, we normalize the crack resistance tensor to a value γ and consider a cleavage anisotropy factor β ∈ [1, 100] as proposed in Clayton-Knap [26,27] for polycrystalline silicon carbide. The crack resistance tensor for 2D and 3D structures is given by respectively, where R grain is a rotation matrix encoding the grain orientation.

Two-dimensional structures
We start with a 2D Laguerre tessellation and planar grain orientations. In our first study, we vary the number of grains as well as the cleavage anisotropy factor β. For each structure, we vary the loading angle in seven equidistant steps from 0 to 90 degrees. For each loading angle, we additionally consider the case where each grain is rotated by 90 degrees. This results in 14 computations per structure. We used the Barzilai-Borwein adaptivity within the ADMM solver for most computations. For some cases, the averaging adaptivity converged faster and we switched to the latter solver choice in those cases. We extracted the mean value as well as the standard deviation of the 14 computations per structure. The crack path for various loading angles is shown in Fig. 3. We see that the crack changes its direction for each grain in order to minimize its surface energy. To close the crack path for a non-axis aligned normal, the crack has to pass the cell several times. Furthermore, we observe local similarities in the crack path for different loading angles. Figure 4a shows the influence of the number of grains on the effective crack energy. For a low number of grains, we observe a rather large standard deviation and no clear trend in the mean value. From 16 2 to 32 2 grains, the standard deviation decreases significantly. This trend continues for 64 2 grains, were the standard deviation decreases even further while the mean value remains within the same range, indicating representativity. We thus find a structure of 32 2 grains to be sufficiently large.
Secondly, we investigate the influence of the anisotropy factor on the effective crack energy for the structure containing 1024 grains, see Fig. 4b. Note that the case β = 1 gives rise to a homogeneous ξ -field, since no local differences arise in the crack resistance. In particular, the effective crack energy becomes γ eff ≡ γ . The effective crack energy increases with an increasing anisotropy factor until a certain threshold is reached at β = 50, beyond this internal contrast, no significant increase is visible. At this threshold, the crack favors the weakest direction in each grain. This is also reflected in Fig. 5. Clear differences in the crack path may be observed between β = 5 and β = 10, see Fig. 5a, b, respectively. For increasing anisotropy, the differences become more subtle, such that the crack paths for β = 20 and β = 50 are almost indistinguishable, see Fig. 5c and see Fig. 5d, respectively.

Three-dimensional structures
We consider generated 3D Laguerre tessellations [67] with random grain orientation in each grain. Similar to the 2D case, we vary both the number of grains and the anisotropy factor. Since 3D simulations are much more costly compared to 2D simulations, we perform only three simulations per microstructure for the size study by investigating a normals in e x , e y and e z -direction, respectively, and only one simulation per cleavage anisotropy factor β for the internal contrast study. Figure 6 shows two cracked microstructures with 216 and 1 728 grains, respectively. Similar to the 2D case, the crack surface changes its orientation in each grain. Figure 7a shows the results of the size study for an anisotropy factor β = 10. The effective crack energy increases for increasing number of grains. For a very small structure containing 27 grains, the deviation between the three loading directions is rather large. For an increasing number of grains, this standard deviation is reduced to a negligible amount. Two main distinctions between the 2D and the 3D case emerge in case of the size study. Firstly, the range of the effective crack energies differ significantly: In the 2D case the effective crack energy ranged around γ eff 1.4 γ , whereas for the 3D case we find γ eff 5.2 γ for β = 10. Secondly, the effective crack energy increases with increasing number of grains in the 3D, at least within the range of our observation, whereas in the 2D case we observed saturation. This indicates that even more than the considered 1 3824 grains could be necessary to reach representative results.
Shifting our focus to the contrast study, see Fig. 7b, we observe a nearly linear correlation between the cleavage factor β and the effective crack energy. This clearly differs from the 2D case, which displayed a saturation. This effect has a geometric origin. In the 2D case, a continuous crack path can be found which passes each grain in the energetically most favorable direction. In 3D, this is not the case, as planar cracks within each grain cannot be combined into a continuous, global crack surface, in general. Hence, in the two-and three-dimensional context, the effective crack resistance of polycrystalline materials differs fundamentally in the context of anisotropic intergranular fracture, modeled with one cleavage parameter. In two spatial dimensions, the parameter β plays a subordinate role (at least if it is sufficiently large), in the three-dimensional case we see clearly that β is an important material parameter, which needs to be identified.

A carbon fiber reinforced polymer
In this section, we investigate a carbon fiber reinforced composite. Carbon fibers are used due to their favorable strength-density ratio. The individual fibers have a strong anisotropy in both elastic and strength properties. As a resulting of to the manufacturing process, they show higher Young's modulus and higher strength in fiber direction compared to the plane perpendicular to it. We use a crack resistance tensor for a fiber oriented in (unit) direction p as assuming an internal contrast of 50 for the crack resistance, as suggested by Pillai et al. [35]. Furthermore, we model the matrix material in our composite with the isotropic crack resistance γ , assuming that the fibers transverse crack resistance is lower than that of the matrix material. The microstructure under consideration contains 290 fibers of 450μm length and 7μm diameter, oriented in an almost unidirectional manner in x-direction with second order orientation tensor diag(0.9, 0.5, 0.5) and a total of filler content of 15%, see Fig. 8a. The structure was generated using sequential addition and migration [68]. Figure 8b shows the crack surface oriented in x-direction. We notice fiber pullout and matrix failure, as well as fiber damage that appears to be non-perpendicular to the fiber direction. The effective crack energy is increased by 50% compared to the matrix material. The crack surfaces in y-and z− direction are shown in Fig. 8c, d, respectively. Both crack surfaces look qualitatively similar. We notice both matrix failure and inter-fiber debonding, since we assumed the fibers  Fig. 7 Influence of the number of grains and of the material contrast on the effective crack energy in the three-dimensional case perpendicular to their orientation to be weaker than the matrix material. For both cases, the effective properties are lower compared to the matrix material and almost equal.

Conclusion
In this work, we presented a numerical method for computing the effective crack energy of a heterogeneous medium with distinct anisotropy via weighted minimal surfaces. We derived the anisotropic maximum flow problem and discussed the implementation into an FFT-based homogenization tool for isotropic fracture. We saw that both the solver framework and the discretization established for the isotropic case serves as a firm foundation for the anisotropic case. Indeed, the extension requires evaluating the projection onto ellipsoids in an efficient manner. This anisotropic extension of the homogenization of the fracture energy enables to treat of additional material classes compared to previous works. Applications were presented for polycrystalline ceramics and carbon-fiber reinforced composites.
In the literature on anisotropic phase-field fracture, 2D polycrystalline materials are often investigated in addition to the 3D case. The anisotropy may be encoded by the cleavage parameter β, which penalizes crack propagation in certain directions and forces the crack path to follow a weak plane. In our homogenization framework, we observed the behavior for the 2D case to fundamentally differ from the 3D case. In two dimensions, a crack is not geometrically restricted to follow the weakest plane through each grain. Therefore, the cleavage parameter β has no further meaning beyond a certain threshold. In the 3D case, on the other hand, neighborhood relations between different grains prohibit the crack to form freely in order to follow the weakest plane for each grain. Therefore, we observe a strong dependence of the effective crack energy on the cleavage parameter β, emphasizing its importance from the viewpoint of materials science.
Additionally, we saw that our framework allows to investigate carbon fiber structures, which show distinct anisotropy within each fiber. This enables modeling additional effects compared to isotropic fibers, studied in previous contributions [53]. With isotropic inclusions either weakening, similar to pores, or toughening with respect to the matrix material can be modeled. The inclusion of anisotropic fibers allows for toughening in one direction, for instance the fiber direction, and weakening in the transverse direction, broadening the spectrum for material design.