Yang Chen

Dmytro Vasiukov

Lionel Gélébart
email: lionel.gelebart@cea.fr

Chung Hae Park

A FFT solver for variational phase-field modeling of brittle fracture

Keywords: Phase-field model, FFT method, Brittle fracture, Fixed-point algorithm

Introduction

It is a challenging issue to simulate crack initiation and propagation within highly heterogeneous materials. Continuum damage models provide a powerful and easy-to-implement solution for this purpose, yet they face a well-known problem of mesh dependency. To address this issue, non-local damage models are usually employed with various regularization techniques [START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: Survey of progress[END_REF], such as spatially averaging state variables (integral-type, e.g. [START_REF] Pijaudier | Non local damage theory[END_REF][START_REF] Jirásek | Nonlocal models for damage and fracture: Comparison of approaches[END_REF]), or constructing the local constitutive relations with the addition of gradient-dependent terms (gradient-type, e.g. [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF][START_REF] Lorentz | Gradient constitutive relations: numerical aspects and application to gradient damage[END_REF]). Another non-local damage model that must be mentioned is the thick level set approach [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF], which deals with crack initiation and propagation with a level set function whose driving force is associated with the released fracture energy that is calculated from the integral over a transition zone between the undamaged and fully damaged zones.

The variational phase-field method is also included in the category of gradient-type models, and is drawing more and more attention, see e.g. [START_REF] Buliga | Energy minimization brittle crack propagation[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]. This approach is established based on the variational framework of [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] with a classical viewpoint of Griffith's theory. The variational approach uses a Mumford-Shah function [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] as the energy potential to be minimized, and the sharp crack surface is smeared using a Γ -convergent approximation, where a characteristic length is involved, see e.g. [START_REF] Ambrosio | Approximation of functionals depending on jumps by Elliptic functionals via gamma-convergence[END_REF][START_REF] Braides | Γ -Convergence for Beginners[END_REF]. In such crack regularization, a state variable is interpreted as a phase field that indicates the damage state of the solid body. A thermodynamically consistent framework has been later outlined by [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF], providing a compact formulation of phase-field model. The descriptive and comprehensive formulation of [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] promoted the application of variational phase-field models in the engineering community. Such damage models have proven capable of simulating crack initiation and propagation without any prescription of the crack geometry or any ad-hoc numerical treatment [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF]. They exhibit also a low mesh dependency and can handle very complex fracture processes, such as crack interaction, branching and coalescence. The efficiency of variational phase-field models has been demonstrated by various applications reported in the literature e.g. [START_REF] Nguyen | Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: direct comparison between in situ testing-microCT experiments and phase field simulations[END_REF][START_REF] Zhang | A modification of the phase-field model for mixed mode crack propagation in rock-like materials[END_REF][START_REF] Carollo | Modeling complex crack paths in ceramic laminates: a novel variational framework combining the phase field method of fracture and the cohesive zone model[END_REF][START_REF] Msekh | Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model[END_REF]. More recently, contributions have been made to extend this approach by considering anisotropic fracture properties [START_REF] Teichtmeister | Phase field modeling of fracture in anisotropic brittle solids[END_REF][START_REF] Bleyer | Phase-field modeling of anisotropic brittle fracture including several damage mechanisms[END_REF] and interface layers between constitutive materials [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF][START_REF] Paggi | Revisiting the problem of a crack impinging on an interface: a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model[END_REF]. All the aforementioned models have been solved using standard finite element method (FEM) solvers.

Although phase-field models are known as relatively mesh independent, a condition on element size must be ensured, i.e. the element size of the fractured zone needs to be smaller than the half of the characteristic length [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF]. This condition is problematic to some extent in the application of damage modeling of highly heterogeneous materials. In fact, the characteristic lengths of the constituents are usually much smaller than the size of the material heterogeneities embedded in the RVE (representative volume element) under investigation, which must be large enough to account for the statistical fluctuations of the microstructure. In addition, the microstructures of heterogeneous materials are usually quite complex, which makes it difficult or even impossible to anticipate the fractured zones even under simple loading conditions. Whilst adaptive meshing is an interesting way to reduce the problem size [START_REF] Zhou | Adaptive phase field simulation of quasi-static crack propagation in rocks[END_REF], a uniform refinement of the FE mesh is usually employed [START_REF] Nguyen | Initiation and propagation of complex 3D networks of cracks in heterogeneous quasi-brittle materials: direct comparison between in situ testing-microCT experiments and phase field simulations[END_REF][START_REF] Nguyen | Modeling of Complex Microcracking in Cement Based Materials by Combining Numerical Simulations Based on a Phase-Field Method and Experimental 3D Imaging[END_REF]. In this context, especially for 3D simulations, the use of standard FEM solvers for the damage phase-field modeling becomes very cumbersome and requires an efficient parallel implementation in order to overcome the computational limits. An example of parallel implementation based on MATLAB routines has been reported in [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro CT imaging[END_REF], where the largest simulation involving about 30 million elements was parallelized over 8 processors on a single workstation and completed within 357 h. This example demonstrates the need to use massively parallel solvers involving a distributed memory implementation in order to perform large-scale simulations in reasonable computation time.

As an alternative to standard FEM solvers, FFT based methods (e.g. [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid reenement[END_REF][START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF][START_REF] De Geus | Finite strain FFT-based non-linear solvers made simple[END_REF]) consist of matrix-free iterative algorithms for solving Partial Differential Equations with periodic boundary conditions: the problem is solved with no need to assemble the FE-type sparse matrix. The methods are well suited for parallelism and can be implemented quite easily in a distributed memory context. Moreover, a regular structured grid is used in FFT methods, so segmented digital images can be directly used as an input of the simulation without any meshing procedure. By virtue of this, FFT methods are quite convenient for image-based modeling, with virtual images or experimental images acquired from 3D Micro-Computed Tomography (µCT) for example or from any other 3D reconstruction technique. This aspect has been demonstrated by a recent work on ceramic matrix composites [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF], where the model was built from a synchrotron X-ray µCT image, involving several billions of elements. Furthermore, FFT methods have been successfully applied to various nonlinear constitutive material behaviors, such as hyper-elasticity [START_REF] Lahellec | Analysis of inhomogeneous materials at large strains using fast fourier transforms[END_REF], plasticity [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF], and visco-plasticity [START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF]. Continuum damage models have also been solved using the FFT technique, see [START_REF] Spahn | A multiscale approach for modeling progressive damage of composite materials using fast fourier transforms[END_REF] for an application under a multiscale framework combined with finite element method, and [START_REF] Wang | Progressive damage analysis of 3D braided composites using FFT-based method[END_REF] for an anisotropic damage model dedicated to textile composites. Although the FFT technique has been widely used in phase-field modeling for microstructural evolution problems (see [START_REF] Khachaturyan | Three-dimensional field model and computer modeling of martensitic transformations[END_REF][START_REF] Hu | A phase-field model for evolving microstructures with strong elastic inhomogeneity[END_REF][START_REF] Moelans | Phase field simulation of grain growth in three dimensional system containing finely dispersed secondphase particles[END_REF][START_REF] Chen | An integrated fast fourier transform-based phase-field and crystal plasticity approach to model recrystallization of three dimensional polycrystals[END_REF][START_REF] Kochmann | Two-scale FE-FFT-and phase-field-based computational modeling of bulk microstructural evolution and macroscopic material behavior[END_REF] among others), to the knowledge of the authors, its application to fracture phase-field model has never been reported.

In this paper, the phase-field model presented in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] is implemented to demonstrate how a FFT solver can be used to solve phase-field models for brittle fracture. The staggered scheme of [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] allows the solution of the mechanical and phase-field problems separately, both using the FFT technique. The paper is organized as follows. First, the phase-field model of Miehe et al. [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] is reviewed. This phase-field model will be then reformulated and solved within a FFT framework. This single edge notched specimen under tension will be selected as a benchmark to study the characteristics of the proposed method. Then, results of several representative numerical examples will be discussed. The capacity of the FFT solver in predicting different crack modes, as well as crack branching, coalescence and interaction, will be demonstrated. Finally, a simulation involving over 32 million voxels will be conducted to illustrate the efficiency of the method in solving large-scale problems.

Review of Miehe's phase-field model

The starting point of Miehe's phase-field model [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] is from smearing the sharp topology of a crack into a diffuse one with an exponential function i.e. d (x) = exp(-|x| /l c ). This diffuse regularization introduces a characteristic length l c and leads to a crack surface density γ expressed by the damage variable d and its spatial gradient ∇d:

γ (d, ∇d) = 1 2l c d 2 + l c 2 |∇d| 2 (1)
The diffuse crack topology can be determined by minimization of the crack surface Γ l :

d (x, t) = Arg{ inf d∈W Γ Γ l } with Γ l = ∫ Ω γ (d, ∇d)d V (2)
subject to the Dirichlet-type constraints W Γ = {d|d (x, t) = 1 at x ∈ Γ (t)}, where Γ (t) represents the crack surface with sharp topology at the instance t, which can be recovered from

Γ = lim l c →0 Γ l (3)
This description of crack topology is then incorporated into the variational framework of [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] to deduce the equations governing the evolutions of both damage and displacement fields. The variational approach introduces an energy functional for a cracked body:

Π (u, d) = E (ε (u) , d) + D(d) (4)
where D(d) is the energy dissipated by the created cracks, and it is defined using the diffuse crack topology based on Griffith's theory:

D (d) := ∫ Ω g c γ (d, ∇d)d V (5)
where g c represents the critical energy released by the creation of a unit crack surface. E(ε (u) , d) in Eq. ( 4) represents the elastic energy stored in the cracked body, and is considered to be degraded by the fracture evolution:

E (ε (u) , d) = ∫ Ω [ g (d) • ψ + 0 (ε(u)) + ψ - 0 (ε (u)) ] d V with g (d) = (1 -d) 2 + k ( 6
)
where k is a small positive parameter for stabilizing the numerical algorithm. Here, in order to take into account the unilateral condition, the strain energy density ψ 0 (ε(u)) of the unbroken body has been split into positive ψ + 0 (ε) and negative ψ - 0 (ε) parts corresponding to the contributions of tension and compression, respectively. The two parts are calculated for isotropic solids by

ψ ± 0 = λ 2 ⟨tr (ε)⟩ 2 ± + µ(ε ± : ε ± ) (7)
where λ and µ are Lamé coefficients, ⟨x⟩ ± = (x ± |x|)/2, and ε ± represents the spectral decomposition of the strain tensor

ε ± = ∑ 3 i=1 ⟨ ε i ⟩ ± n i ⨂ n i
, with ε i and n i being the eigenvalues and eigenvectors of the strain tensor. The strain decomposition is completed using two projection tensors. The demonstration of deducing an explicit form of the projection tensors is given in Appendix A.

The multi-field problem can be solved by minimizing the energy functional (Eq. ( 4)), leading to the local governing equations:

(a) :

g c l c [ d -l 2 c ∆d ] = 2(1 -d)H(ε) (b) : ⎧ ⎪ ⎨ ⎪ ⎩ div (σ (u, d)) = 0 σ = g (d) [ λ ⟨tr (ε)⟩ + I + 2µε + ] + [ λ ⟨tr (ε)⟩ -I + 2µε - ] ε = (∇u + ∇ T u)/2 (8)
where σ (u, d) is the local stress, which follows a material constitutive law depending on the displacement u (x) (or strain ε(x)) field and fracture phase d(x) field. I denotes identity matrix. H(ε) represents the so-called history field of maximum positive elastic energy:

H (ε (x, t)) : = max τ ∈[0,t] ψ + 0 (ε(x, t)) (9)
Table 1 Staggered scheme for solving multi-field problem using FFT method.

Initialization

The strain field ε 0 (x), history field H 0 (x) and phase field d 0 (x) are known Loop while t n+1 ≤ T , given ε tn (x), H tn (x) and d tn (x) 1.

Solve the phase-field problem (Eq. ( 8).a or [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF])

⇒ d t n+1 (x) 2.
Solve the mechanical problem (Eq. ( 8).b) ⇒ ε t n+1 (x) 3.

Update the history field (Eq. ( 9)) ⇒ H t n+1 (x)

Additionally, in order to improve the stability of the numerical solution, a viscous regularization can also be implemented as proposed in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] by adding a viscous term in Eq. ( 8).a:

g c l c [ d -l 2 c ∆d ] + η ḋ = 2(1 -d)H(ε) ( 10
)
where η is viscous parameter. It should be noted that this parameter has a physical unit of Pa s, which is the same as the dynamic viscosity in fluid mechanics. As a consequence, the order of magnitude of this viscous parameter must be chosen with respect to the physical units used in practical implementations.

As mentioned in the introduction, differently from FE implementations where matrix assembly is usually required, the FFT method iteratively solves the governing equations (Eqs. ( 8) and ( 10)) in a local fashion, in either real or Fourier space. The present work uses an existing FFT iterative procedure for the mechanical problem, while it proposes a new procedure to solve the phase-field equation.

Phase-field model under FFT framework

In this work the staggered scheme of [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] is employed to solve the phase-field (Eq. [START_REF] Bourdin | The variational approach to fracture[END_REF].a or Eq. ( 10)) and mechanical (Eq. [START_REF] Bourdin | The variational approach to fracture[END_REF].b) problems separately, both using the FFT technique. The staggered scheme consists of solving the phasefield problem by considering the mechanical field (strain, stress) to be known at a given time step, and vice-versa for the next time step. An overview of the staggered scheme is given in Table 1. FFT-based solvers are employed in both steps 1 and 2, respectively corresponding to the phase-field and mechanical problems.

Solution of the mechanical problem

The unknown in the mechanical problem (Eq. ( 8).b) is strain instead of displacement. Periodic boundary conditions (Eq. [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF].a) are prescribed to the unit cell Ω subjected to overall strain or stress load (Eq. [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF].b):

(a) : ε ( x + L p ) = ε (x) ; x ∈ Ω ∞ (b) : ⟨ε (x)⟩ = E; x ∈ Ω (11)
where L p stands for the three periodicity vectors, Ω ∞ denotes an infinite medium in which the unit cell Ω is periodically repeated. ⟨ * ⟩ represents the spatial mean over the unit cell Ω , and E is the prescribed overall strain tensor. Considering the damage field as "frozen" (constant) at a given time step, the mechanical problem (Eq. ( 8).b) becomes an elastic problem whose solution can be achieved by iteratively solving the Lippmann-Schwinger equation:

ε (x) = -Γ 0 * τ (x) + E with τ (x) = σ (x) -C 0 : ε(x) (12)
where τ (x) is called the polarization field, and Γ 0 is the fourth order Green operator for strain that is calculated in Fourier space as a function of the elastic constants of the chosen reference material C 0 . The stress field σ (x) is computed in real space using the constitutive equation in Eq. ( 8).b, hence the effects of fracture evolution and unilateral condition are taken into account.

Various FFT formulations are readily available for solving the mechanical problem, hence the implementation of this part is not presented in this paper. Readers are referred to the original papers, see e.g. [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Brisard | FFT-based Methods for the mechanics of composites: a general variational framework[END_REF][START_REF] Zeman | A finite element perspective on nonlinear FFT-based micromechanical simulations[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF]. In this work, the mechanical problem is chosen to be solved using a fixed-point algorithm as proposed in the basic scheme of [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] combined with a modified discrete Green operator [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF]. As demonstrated in [START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF], the use of this discrete Green operator is strictly equivalent to the use of hexahedral finite elements with reduced integration (i.e. 8 nodes and one Gauss point in the middle of the element) in a standard FEM code. In addition to this, a convergence acceleration procedure (see [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF]) is also employed to improve the performance of the fixed-point algorithm in mechanical solution. Even though these two improvements significantly reduce the dependence of convergence rate on the choice of the reference material C 0 , we still choose to follow the suggestion of [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF], i.e. λ 0 = (min (λ) + max (λ))/2 and µ 0 = (min (µ) + max (µ))/2, to achieve an optimal convergence rate of the mechanical solution.

In practice, instead of a strain-driven method as in Eq. ( 12), where the applied load is the average strain, it is usually more convenient to apply mixed loading conditions, i.e. average stress for some components and average strain for the others. This is made possible in the code AMITEX through a slightly modified algorithm, similar to that reported in [START_REF] Kabel | Mixed boundary conditions for FFT-based homogenization at finite strains[END_REF]. This consists of re-evaluating the complete average strain at each fixed-point iterate, so that the average stress conditions are fulfilled at convergence.

Fixed-point algorithm for the phase-field problem

Now we introduce a fixed-point iterative algorithm nested in the staggered scheme to solve the phase-field problem using the FFT technique (step 1 in Table 1). The non-viscous form can be recovered by setting η = 0, hence only the solution of the viscous form (Eq. ( 10)) is presented in the following.

The objective of the fixed-point iteration is to find the phase field d t n+1 (x) at the time step t n+1 based on the history field H t n (x) and the damage field d t n (x) that are known at the previous time step t n . For the sake of simplicity, the superscript * t n+1 will be hereafter omitted, i.e. d (x) := d t n+1 (x). In order to keep the formulation as the most generic possible, we consider that all the parameters (g c , l c , η) involved in Eq. ( 10) can vary from one point to another if the unit cell is heterogeneous. In other words, they are functions of the coordinates x. Replacing ḋ by d-d tn δt , with δt the time step of the simulation, Eq. ( 10) can be rewritten as

( 1 l 2 c (x) + 2H t n (x) g c (x) l c (x) + η(x) δtg c (x) l c (x) ) d (x) -∆d(x) = 2H t n (x) g c (x) l c (x) + η(x) δtg c (x) l c (x) d t n (x) (13)
and simplified into the following form:

A t n (x)d (x) -∆d(x) = B t n (x) (14)
with

A t n (x) = 1 l 2 c (x) + 2H t n (x) g c (x) l c (x) + η(x) δtg c (x) l c (x) ; B t n (x) = 2H t n (x) g c (x) l c (x) + η(x) δtg c (x) l c (x) d t n (x) (15)
These two field variables are calculated for each time step t n . By introducing a homogeneous parameter A n 0 , Eq. ( 14) becomes

A t n 0 d (x) -∆d(x) = χ(x) (16)
with χ (x) hereafter called polarization field for phase-field problem:

χ (x) = B t n (x) - ( A t n (x) -A t n 0 ) d(x) (17)
The fixed-point algorithm is now employed to solve Eq. ( 16) in order to find a new damage field d k+1 (x), by assuming the polarization fieldχ k (x) as known at the previous iteration step k. The Laplacian operator is computationally expensive in real space, yet in Fourier space it can be easily computed by local multiplication with wave (frequency) vector ξ

F (∆d(x)) = -(ξ • ξ ) • F (d(x)) (18)
where F( * ) represents the Fourier transform of a function * , and hereafter will be denoted by the symbol * . Thus, the solution of Eq. ( 16) can be expressed in Fourier space as:

dk+1 (ξ ) = χ k (ξ ) A n 0 + ξ • ξ (19)
Table 2 outlines the proposed algorithm that needs to be executed at each time step of the overall staggered scheme. We distinguish here the two superscripts * k and * t n : the former represents the index of fixed-point iteration step, while the latter is the index of time step. Indeed, this fixed-point loop is undertaken at each time step t n+1 to reach a converged phase field d k+1 (x). However, the fields of A t n (x), B t n (x) and d t n (x) remain known and unchanged during the fixedpoint loop. The loop is stopped once the convergence criterion is reached, i.e. the residual of the governing equation (Eq. ( 16)) is smaller than the predefined tolerance ϵ tol . In practice, we calculate the residual using the polarization

Initialization

Initialize the fracture phase field d k (x) = d tn (x) Compute the polarization field χ k (x) (Eq. ( 17)) Loop while ϵ > ϵ tol , A tn (x) , B tn (x) , A tn 0 , d tn (x) being known and unchanged 1.

FFT:

χ k (x) → χ(ξ ) 2.
Compute a new fracture phase field d(ξ ) in Fourier space (Eq. ( 19)) 3.

Inverse FFT:

dk+1 (ξ ) → d k+1 (x) 4.
Update the polarization field χ k+1 (x) (Eqs. ( 17)) 5.

Convergence test: Compute the residual ϵ of the governing equation (Eq. ( 20)).

fields of two sequential iterations:

ϵ =   χ k+1 (x) -χ k (x)   2   χ k+1 (x)   2 ( 20
)
where ∥ * ∥ 2 represents the L-2 norm over the unit cell. The tolerance is set to ϵ tol = 10 -6 in the present implementation.

Remarks. As mentioned in Section 3.1, a modified discrete Green operator is used for the mechanical problem.

The operator is built on the definition of discrete derivatives in real space. In practice, we use the finite-difference based definition that is reported in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] (the so-called 'rotated scheme') and in [START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] (the linear hexahedral FE with reduced integration in the Appendix of the reference). Moving from the classical discrete Green operator [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] to the modified Green operator [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF][START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] is simple: replacing the frequency vectors that define the Fourier space by the modified frequency vectors when applying the discrete Green operator. In the present work, the same finite-difference based definition and consequently the same modified frequency vectors are used for solving both the mechanical problem and the phase-field problem (Eq. ( 19)). It is useful to mention that for the mechanical problem, the modified discrete Green operator, compared to the classical one, drastically reduces spurious oscillations and allows for solving problems with infinite elastic contrasts.

A t n 0 in Eq. ( 16) has a great importance for the proposed fixed-point algorithm, because the choice of this parameter can significantly affect the convergence rate. By analogy to the basic scheme of [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] for the mechanical problem, the mean of the maximum and minimum values of A t n (x) (Eq. ( 21)) has been chosen in the present work, providing a relatively fast convergence rate.

A t n 0 = min ( A t n (x) ) + max ( A t n (x) ) 2 (21)
A detailed flowchart of the overall algorithm is given in Appendix B to provide a global view of the proposed FFT solver. Two fixed-point loops are nested within the staggered scheme, solving the mechanical and phase-field problems, respectively. Since fixed-point algorithms do not require computation of the tangent stiffness matrix, compared to tangent-based algorithms like Newton-Raphson's method, it will be easier to extend the FFT solver proposed herein to other models, such as damage with anisotropy [START_REF] Teichtmeister | Phase field modeling of fracture in anisotropic brittle solids[END_REF][START_REF] Bleyer | Phase-field modeling of anisotropic brittle fracture including several damage mechanisms[END_REF], or damage with plasticity [START_REF] Miehe | Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory[END_REF][START_REF] Alessi | Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples[END_REF].

Parallel implementation

In this FFT-based solver of the phase-field problem every step involves local operations only, except the FFT and inverse FFT. In addition, open sources of powerful FFT packages are available to complete the FFT and its inverse in a parallel fashion. Therefore, the proposed algorithm inherits the simplicity in parallel implementation of the classical FFT methods. The present implementation of the phase field problem has been introduced in a FFT-based code "AMITEX" solving mechanical problems at small and finite-strains in a massively parallel context. The code relies on the 2DECOMP&FFT library [START_REF] Li | 2DECOMP & FFT-A Highly Scalable 2D Decomposition Library and FFT Interface[END_REF] for the distributed memory implementation and readers interested in this aspect are invited to consult the documentation found on [START_REF]AMITEX_FFTP[END_REF] and some practical aspects such as the loading conditions and convergence criteria, can be found in [START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF].

Convergence acceleration

A convergence acceleration technique has been incorporated in both phase-field and mechanical problems to speed up the convergence of the fixed-point loops. This technique has readily been implemented in AMITEX [START_REF]AMITEX_FFTP[END_REF]. The main idea is to propose a better (or accelerated) solution of the strain or damage field, according to the four solutions and corresponding residuals obtained at previous steps. This "acceleration" procedure is applied every three steps of the fixed-point loop. A more detailed description is given in [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF].

Numerical examples

The proposed algorithm has been implemented for general 3D boundary value problems. The unit cells of the following "2D" examples (Sections 4.1-4.5) consist of a slab of one voxel thickness in the third dimension. Due to periodic boundary conditions, the out-of-plane strain components are equal in the whole slab, which is equivalent to a generalized plane strain condition. In Section 4.1, a single edge notched specimen is used as a benchmark to study the effects of different parameters. Similar specimens with inhomogeneous l c or g c are tested in Section 4.2. Then, mode II and mode III shear cracks are simulated in Section 4.3. Symmetric and asymmetric double edge notched specimens are studied in Section 4.4 to assess the capacity of predicting crack interaction. In Section 4.5, crack branching and coalescence are simulated within a notched bi-material specimen. Finally, a 3D example is given in Section 4.6 to demonstrate the potential of the proposed algorithm in large-scale simulations.

In order to conduct tests under tensile or shear loads, mixed loading conditions (see e.g. [START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF]) are employed. For example, a tensile test in X-axis is simulated by incrementally increasing the average strain component ε x x and keeping the average stress components σ yy , σ zz , σ x y , σ x z , σ yz equal to zero. The stability parameter k (see Eq. ( 6)) is set to 10 -6 for all the representative examples presented in the paper.

Benchmark study: single edge notched specimen under tension

The single edge notched specimen under tension is selected as a benchmark in order to examine the influences of various parameters. Fig. 1 depicts the geometry of the unit cell. To avoid the influence of periodic boundary conditions on the damage evolution, margins with null elastic properties are added to the two sides. Both the crack and margin voxels are defined as void, with elastic constants equal to zeros and l cr k c = 0.015 mm, g cr k c = 2.7 • 10 -3 kN/mm. Elastic constants of the solid material are λ = 121.15 kN/mm 2 and µ = 80.77 kN/mm 2 and fracture energy is g c = 2.7 • 10 -3 kN/mm. In the following macroscopic stress-strain curves, the macroscopic stress is evaluated as the average of the stresses in the solid voxels, i.e. notch voxels and margin voxels are excluded from this calculation.

Comparison with FEM solution

Using this benchmark specimen, the result from the FFT solver is compared with that published in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. The characteristic length is chosen as l c = 0.015 mm, which is the same as in the reference. In our simulations, a unit Fig. 2. Macroscopic stress-strain curve of the single edge notched tensile test, compared with that from [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and the result obtained using the ABAQUS implementation of [START_REF] Molnár | 2D and 3D Abaqus Implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF] . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) cell with voxel size of h = 0.005 mm is subjected to uniaxial tension in the Y-axis up to ε yy = 0.005 with initially a relatively large increment of ∆ε yy = 5 • 10 -5 and then the strain increment is refined to ∆ε yy = 5 • 10 -7 until the final failure. Note that Miehe's result was obtained from a 2D implementation. Therefore, in our simulation the unit cell is subjected to an equivalent loading condition, i.e. incrementally increasing ε yy , keeping σ x x , ε zz , σ x y , ε x z , ε yz = 0.

According to the results shown in Fig. 2, the strain to failure of both simulations is quite consistent. However, the slope of the FFT solution (black solid curve) is different from the reference one [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] (red solid curve). This difference appears from the very beginning (even in the elastic regime). Therefore, it should not be due to an error in the solution of the phase-field model, but is more likely related to some other difference(s), for instance in the mesh configuration between the FFT and FEM solutions. In fact, in the FFT solution the notch has been considered as void voxels with a thickness of one voxel (0.005 mm), whereas it had no thickness in the FEM solution of [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. To check this point, a FEM simulation with a notch of thickness of 0.005 mm has been conducted using the FEM implementation of [START_REF] Molnár | 2D and 3D Abaqus Implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF] in ABAQUS/UEL. 4-node bilinear reduced integration with hourglass control elements of the size of h = 0.005 mm has been used in this simulation. The result is depicted by the blue dashed curve, whose slope agrees well with the FFT solution. To sum up, these two comparisons demonstrate the accuracy of the proposed method.

Convergence acceleration and viscous regularization

The efficiency of the convergence acceleration technique in solving the phase-field problem is examined for the same unit cell (h = 0.005 mm) as in Section 4.1.1. The macroscopic loading condition is applied by keeping σ x x , σ zz , σ x y , σ x z , σ yz equal to zero and incrementally increasing ε yy up to ε yy = 0.005 with a relatively large increment of ∆ε yy = 5 • 10 -5 and then the strain increment is refined to ∆ε yy = 5 • 10 -7 until the final failure. The viscous parameter η is set to zero. The characteristic length is l c = 0.015 mm. Two simulations are conducted, with and without the convergence acceleration technique. As shown in Fig. 3, the stress-strain curves are identical for both simulations, while the convergence acceleration technique reduces the maximum number of iterations for the phase-field solution by a magnitude of 10 (from ∼2300 to ∼250). This demonstrates that the convergence acceleration technique can significantly improve the convergence rate of the fixed-point algorithm that is proposed in the present paper for the phase-field problem.

Using the same unit cell, material parameters and loading conditions, another simulation has also been conducted with the convergence acceleration disabled for the mechanical problem. The number of iterations for the mechanical solution becomes very large (larger than several thousands) even from the early stage of the load, making the computation time too long, hence the computation was interrupted before its completion. This confirms the efficiency of the technique for accelerating the fixed-point algorithm that is commonly used for the mechanical problem.

Next, the effect of the viscous regularization is investigated. Two simulations are performed using different viscous parameters, i.e. η = 0 and 10 9 Pa s. Other parameters (∆ε yy , l c , h) are the same as in the previous study. The convergence acceleration technique is disabled in the phase-field problem in this study. Fig. 4 shows that the viscous regularization with η = 10 9 Pa s slightly retards the aggressive crack growth in the post critical regime, as already presented in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. On the other hand, the viscous regularization helps to reduce the number of iterations for the phase-field solution from the early stage of loading (compared to the convergence acceleration technique, see Fig. 3), and the maximum number of iterations is reduced from ∼2300 to ∼600 by virtue of the viscous regularization. More simulations using different values of the viscous parameter η have also been tested and the stress-strain curves are gathered in Fig. 5. The retarding effect of the viscous regularization on the aggressive crack growth becomes more and more obvious as η increases. However, a too high value of η (e.g. 10 12 Pa s) can lead to unrealistic prediction of the critical stress, hence the viscous parameter must be chosen carefully.

We choose to conduct the further numerical tests only with the convergence acceleration technique, but with the viscous regularization disabled, i.e. η = 0 Pa s.

Load increment, characteristic length and voxel size

Due to the specificities of FFT method compared to standard FE method (e.g. loading and boundary conditions) and in order to convince the readers who are more familiar with standard FE implementations, we check how various parameters of the phase field model influence the prediction of material fracture in the FFT implementation. In this section, the effects of strain increment ∆ε, characteristic length l c and voxel size of the unit cell h will be discussed. Simulation results similar to those from FE solutions will be obtained.

In order to assess the influence of strain increment on the solution, the unit cell is loaded in the tensile Y -axis (ε yy ) with a rate of ∆ε yy = 5 • 10 -5 up to ε yy = 0.005, and then the strain increment is refined with different values of ∆ε yy . The other stress components (σ x x , σ zz , σ x y , σ x z , σ yz ) are kept equal to zero. The unit cell has a voxel size of h = 0.005 mm and the characteristic length is chosen as l c = 0.015 mm. Fig. 6 shows the stress-strain graphs of these simulations. A large strain increment retards the aggressive softening behavior, and even affects the prediction of critical stress (e.g. ∆ε yy > 5 • 10 -6 in the present case). This effect is similar to the effect of the viscous regularization (Fig. 5) and has also been reported in [START_REF] Molnár | 2D and 3D Abaqus Implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF]. Therefore, the strain increment must be chosen with care, keeping in mind the fact that a too small increment makes the computation slow, while a too large one results in unrealistic prediction.

The effect of the characteristic length l c is analyzed in a unit cell with voxel size of h = 0.005 mm. For this set of simulations, the increment of ∆ε yy = 5 • 10 -7 is fixed for the second loading stage. According to Fig. 7.a, a smaller characteristic length leads to a higher strength but a smaller strain to failure. This confirms that the characteristic length should be considered as a material parameter instead of an algorithmic one, as also mentioned in [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF]. Note that a continuous softening of the stress-strain curve is observed for each simulation and the phenomenon becomes more obvious as l c increases. This is attributed to the evolution of the damage parameter over a large area, even quite far from the crack tip (see the light blue color on the right side of each damage field, Fig. 7.b). Close to the final rupture, a small jump is observed on the curves. This could be attributed to the spurious oscillations of the damage field observed around the crack tip (see Fig. 7.b). However, the effect of this oscillation becomes negligible when the crack starts to propagate. This issue would possibly be eliminated in the future by considering a regularization of the predefined crack notch.

Fig. 8 illustrates the effect of the voxel size h. Unit cells with different resolutions are tested using the same characteristic length l c = 0.015 mm and the same strain increment ∆ε yy = 5 • 10 -7 for the second loading stage. The stress-strain curves are almost identical for the different simulations, except that with h = 0.01mm. This demonstrates that the phase-field model is independent on the grid resolution if the condition of h ≤ l c /2 is satisfied.

Furthermore, the numbers of iterations and computation time for the simulations with unit cells of different resolutions are collected in Table 3. The time per iteration for either phase-field or mechanical problem increases with the number of voxels in a quasi-linear trend, as shown in Fig. 9. This is a promising property of the FFT solver compared to standard FEM solvers. Unfortunately, the numbers of iterations also increase as the grid is refined, making the total computation time increase with the problem size by an order slightly higher than one. In addition, it is useful to mention that the time per iteration is shorter (about 10 times) for the phase field problem, which is not surprising since the number of unknowns is 6 (strain components) for the mechanical problem and only one for the phase-field. On the other hand, the number of iterations is larger (2∼3 times) for phase-field problem than for mechanical problem. Overall, more time is spent for solving the mechanical problem than for solving the phase-field problem. In the phase-field model, the fracture behavior is primarily determined by two parameters: the characteristic length l c and the fracture energy g c . Now, we check the capability of the proposed method for solving different cases with inhomogeneous fracture properties. Single edge notched specimens are used for this purpose. Fig. 10.a depicts the geometry, which is similar to the one used in previous simulations (Fig. 1). The two parts of the solid body have identical elastic constants (λ = 121.15 kN/mm 2 , µ = 80.77 kN/mm 2 ), but different characteristic lengths (l ci ) or fracture energies (g ci ). The predefined crack and margin voxels have zero elastic constants with l cr k c = 0.02 mm, g cr k c = 2.7•10 -3 kN/mm. The voxel size of the unit cell is h = 0.01mm. The loading condition is the same as that in Section 4.1.2. Fig. 10.b shows the macroscopic stress-strain curves of the simulations with different configurations on l c and g c , together with the corresponding crack patterns shown in Fig. 10.c-f.

The characteristic length l c clearly controls the diffuse crack topology as shown in Fig. 10.c and d. A greater l c produces a higher regularization of crack surface, hence a more diffuse crack topology. As for the tests with different g c , when the crack grows from high g c to low g c material phase, the smeared zone becomes larger, which is unexpected and the reason of this is unclear yet. When the crack enters into the phase with higher g c from a lower g c , its growth rate decreases, which explains the second peak of the red dash curve in Fig. 10.b.

Cracked body under shear loads

In order to check the prediction of mode II and mode III cracks, the unit cell of a body with periodically distributed cracks is used as depicted in Fig. 11.a. No void margins are added, so the periodicity is kept in all directions. As a result, the crack can be located at any arbitrary position within the unit-cell as shown in Fig. 11.b. This configuration illustrates the specificity of FFT method, i.e. periodic boundary conditions are intrinsically applied to the unit cell. The voxel size is h = 0.005 mm. The crack voxels are given null elastic properties with l cr k c = 0.015 mm, g cr k c = 2.7 • 10 -3 kN/mm. The solid material parameters are λ = 121.15 kN/mm 2 , µ = 80.77 kN/mm 2 , l c = 0.015 m and g c = 2.7 • 10 -3 kN/mm. The macroscopic mode II (mode III resp.) shear load is applied by incrementally increasing the strain component ε x y (ε yz resp.) and keeping the stress components σ x x , σ yy , σ zz , σ x z , σ yz (σ x y resp.) equal to zero. Once again, using a one voxel thick slab, these 3D simulations are equivalent to 2D simulations with generalized plane strains (with constant out of plane strain components). At the first loading stage, the shear strain is increased up to ε x y = 0.008 (ε yz = 0.015 resp.) with a load increment of ∆ε x y = ∆ε yz = 5 • 10 -5 . Then at the second loading stage, the load increment is reduced to ∆ε x y = 5 • 10 -6 (∆ε yz = 5 • 10 -6 resp.).

As shown in Fig. 12.a, two cracks appear due to periodic boundary conditions. The mode II crack path first follows a straight line with an angle of 45 • with respect to the plane of the predefined notch. Then their orientations change as they approach the boundaries. This is explained by the increasing interaction between the periodic cracks. On the other hand, two mode III cracks also periodically nucleate at the ends of the predefined notch, and they grow in the plane of the notch (Fig. 12.b).

Symmetric and asymmetric double edge notched tensile specimens

In this section, crack interaction is tested in double edge notched tensile specimens. Both specimens with different relative positions of edge notches, as depicted in Fig. 13, are studied. The macroscopic stress-strain curves are shown in Fig. 14 for both specimens. They are almost identical before the final stage (point C). At the final stage, the symmetric specimen reaches its total failure, which corresponds to the connection of the two coplanar cracks (see Fig. 15.a). In contrast, the asymmetric specimen continues to sustain a small value of tension because the two non-coplanar cracks avoid their connection (Fig. 15.b). As shown in Fig. 15, when the two cracks are relatively far from each other, there is no obvious difference in crack interaction between the two specimens, so the macroscopic responses are identical. However, when the two cracks approach each other, their path depends on their relative position, leading to the difference in macroscopic behavior of the two specimens. The crack patterns of both specimens are in excellent agreement with the literature [START_REF] Molnár | 2D and 3D Abaqus Implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF][START_REF] Melin | Why do cracks avoid each other?[END_REF].

Crack branching and coalescence: notched bi-material specimen under tension

We examine here the capability of the proposed algorithm for predicting crack branching and coalescence. The geometry of the unit cell and the material parameters for the two solid parts used in this test are depicted in Fig. 16.a. The elastic constants of margin and crack voxels are set to zero with l cr k c = 0.3 mm, g cr k c = 1 • 10 -3 kN/mm. The upper part is 10 times stiffer and tougher than the lower part. Macroscopic tension in the X -axis is applied up to ε x x = 0.002 with an increment of ∆ε x x = 2 • 10 -5 and then the increment is refined to ∆ε x x = 5 • 10 -6 until failure. The other stress components (σ yy , σ zz , σ x y , σ x z , σ yz ) are kept equal to zero. Fig. [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF].c shows that the simulation captures both crack branching and coalescence. A crack initiates from the notch (Fig. 16.c-1), corresponding to the slope change at point 1 on the stress-strain curve. Then at point 2, it bifurcates into two branches (Fig. 16.c-2) when it is close to the interface of the two solid parts. This bifurcation does not induce an obvious change of slope on the stress-strain curve. Finally, as the load continues to increase, the two branches penetrate into the upper part and then join each other (Fig. 16.c-3-4, crack coalescence), leading to the final failure.

3D example: Continuous fiber reinforced composite with void

To check the capability of the FFT solver for large-scale simulations, we conduct a numerical test on a unit cell of a continuous fiber composite. Being periodic in three axes, the idealistic geometrical model consists of a unidirectional composite with periodically distributed fibers oriented at 45 • as shown in Fig. 17 = 60 • 10 -3 kN/mm, g matri x c = 1 • 10 -4 kN/mm. The unit cell is loaded in tension along the X -axis up to ε x x = 0.035 with rate of ∆ε x x = 10 -4 and then the increment is refined to ∆ε x x = 5 • 10 -6 until the final failure, with other stress components (σ yy , σ zz , σ x y , σ x z , σ yz ) being zero. The macroscopic stress-strain curve is shown in Fig. 18, together with the numbers of iterations for solving either the phase-field or the mechanical problem at each load increment. Both the phase-field and mechanical solutions require more iterations when the material failure occurs. This trend is similar to that observed previously in Section 4.1.2 (Figs. 3 and4). After the material failure (ε x x > 0.0412), the unit cell is still able to withstand load. More particularly, the slope of the last part of the stress-strain curve remains constant, and its extension passes through the axis origin. This suggests that no damage evolution occurs at this last loading stage, which can be further confirmed by the crack patterns shown in Fig. 19,i.e. no crack growth is observed between the loading levels 5 and 6. As a consequence, the numbers of iterations required at this loading stage decrease drastically. It should be noted that this last stiffening stage is due to the fact that the fibers have a much higher fracture energy making them undamaged within the unit cell. The growth of crack topology is shown in Fig. 19. The simulation captured the crack initiation and the influence of fibers on the crack propagation.

This simulation with over 32 million voxels has been repeated three times on a cluster using one, four or eight nodes (28 cores per node). In other words, the simulation has been parallelized over 28, 112 and 224 cores. Each node consists of two Intel Xeon 6132 processors (14 cores per processor, 2.6 GHz, Slylake R ⃝ technology). An intel OmniPath 100 network is used for inter-node communications. The code has been compiled with the intel 15 compiler together with the intel MPI library, version 5.

The computation time (wall time) and efficiency are summarized in Table 4, the efficiency being defined as E(N ) = T 1 /(N T N ), with T 1 and T N the wall time on 1 node and N nodes respectively. According to Table 4, even if parallelized on 28 cores (1 node), the computation time is still high (36 h), making more nodes necessary. By virtue of a distributed memory parallel implementation, using 4 nodes and 8 nodes has reduced the computation time to 11.5 h (3 times less) and 7 h (5 times less), respectively. Indeed, the efficiency is quite satisfactory when using 4 nodes (0.78) and still acceptable on 8 nodes (0.65). It must be mentioned that the decrease of the efficiency is related not only to the increasing time spent during the communication between processors, but also to the increasing number of iterations. For instance, the number of iterations for the mechanical problem has increased from 40008 on one node to 47667 on eight nodes. This increase has been essentially observed during the stress jump due to the unstable propagation of the crack. The reason for this increase is still under investigation. One idea is that during the instable crack propagation, the algorithm is more sensitive to numerical errors, and furthermore, the amount of numerical errors is more important in a parallel implementation.

To conclude, even if it cannot be really called as a "large-scale" simulation, this "medium-scale" simulation demonstrates the ability of AMITEX to use distributed memory architectures, and to take advantage of High Performance Computing to run massively parallel simulations as reported in [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF] with a problem of more than 7 billion voxels.

Conclusions

A FFT solver of a variational phase-field model for brittle fracture has been proposed in this paper. Relying on the staggered update scheme of [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], the phase-field and elastic problems have been separately solved both using the FFT method with fixed-point algorithm together with a convergence acceleration procedure. Benefitting from the FFT technique, the proposed method is easy to parallelize and convenient for image-based modeling. A convergence acceleration technique has been added to the fixed-point algorithms to improve the computational performance. The influences of this technique, as well as other parameters involved in the phase-field model, have been analyzed using a single edge notched specimen.

Numerical examples have demonstrated that the proposed algorithm is stable and convergent for predicting not only simple mode I, II and III cracks, but also complex cases such as crack interaction, branching and coalescence. Introduced in the AMITEX code and taking advantage of its distributed memory parallel implementation, a simulation with over 32 million voxels has demonstrated the attractive potential of the solver for large-scale simulations. As a conclusion, the FFT solver proposed in this paper is an efficient alternative to standard FEM solvers of phase-field models for brittle fracture, especially in the context of image-based modeling where a large number of voxels are usually involved.

ε + = ( Q t • s • Q) • ε (A.5)
Thus, the projection tensors P + , P -are constructed by

P + = Q t • s • Q; P -= I -P + (A.6)
where I is unitary 3 × 3 matrix. The explicit expression of the projection tensor P + is:

P + = ⎛ ⎜ ⎜ ⎝
s 1 q 11 q 11 + s 2 q 21 q 21 + s 3 q 31 q 31 s 1 q 11 q 12 + s 2 q 21 q 22 + s 3 q 31 q 32 s 1 q 11 q 13 + s 2 q 21 q 23 + s 3 q 31 q 33 s 1 q 11 q 12 + s 2 q 21 q 22 + s 3 q 31 q 32 s 1 q 12 q 12 + s 2 q 22 q 22 + s 3 q 32 q 32 s 1 q 12 q 13 + s 2 q 22 q 23 + s 3 q 32 q 33 s 1 q 11 q 13 + s 2 q 21 q 23 + s 3 q 31 q 33 s 1 q 12 q 13 + s 2 q 22 q 23 + s 3 q 32 q 33 s 1 q 13 q 13 + s 2 q 23 q 23 + s 3 q 33 q 33 ⎞ ⎟ ⎟ ⎠

(A.7)
where q i j is the component of Q. In practice, strains are not stored in the form of tensor but of vector: Therefore, the projection tensor must also be adjusted:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ε 1 = ε 11
P + = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (A.9)
where p i j is the component of P + . The negative projection tensor is calculated by P -= I -P + , where I is unitary 6 × 6 matrix. Thus, the strain vector can be decomposed into positive and negative parts by multiplying with the projection tensors:

⎧ ⎨ ⎩ ε ± = P ± • ε in the form of tensor ε ± = P ± • ε in the form of vector

Fig. 1 .

 1 Fig. 1. Geometry of the unit cell used in the single edge notched tensile test.

Fig. 3 .

 3 Fig. 3. Efficiency analysis of the convergence acceleration technique: stress-strain curves together with the numbers of iterations used for phasefield solution for the simulations with and without convergence acceleration technique.

Fig. 4 .

 4 Fig. 4. Viscous regularization effect for phase-field solution with different values of η: stress-strain curves together with the number of iterations.

Fig. 5 .

 5 Fig. 5. Effect of the viscous parameter η on macroscopic stress-strain curves.

Fig. 6 .

 6 Fig. 6. Effect of stain increment ∆ε yy on macroscopic stress-strain curves.

Fig. 7 .

 7 Fig. 7. Effect of the characteristic length l c : (a) stress-strain curves of the simulations with different values of l c ; (b) damage field at the three load levels marked on the stress-strain curve for the simulation with l c = 0.060 mm.

Fig. 8 .

 8 Fig. 8. Voxel size h effect on macroscopic stress-strain curves.

Fig. 9 .

 9 Fig. 9. Time per iteration for either the phase-field or mechanical problem as a function of number of voxels involved in the simulation.

Fig. 10 .

 10 Fig. 10. Crack patterns in the single edge notched tensile specimen with inhomogeneous fracture energy or characteristic length: (a) geometry of the unit cell; (b) macroscopic stress-strain curves; (c-f) crack patterns for unit cells with different configurations of l c and g c .

Fig. 11 .

 11 Fig. 11. (a) Geometry of the unit cell used for shear tests. (b) A solid medium with periodic cracks: the red square indicates the unit cell chosen in the present study and it is equivalent to any choice with the same dimensions, e.g. the black square . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

  . The volume fraction of fibers has been chosen as 20%, with fiber diameter of 15 µm. An initial defect, represented by a pennyshaped void, is inserted between the fibers. The unit cell has dimensions of 401 × 401 × 201 voxels (32 321 001 voxels in total) with a voxel size of h = 0.25 µm. The elastic constants of the defect are set to zero with l cr k c = 2 µm, g cr k c = 1 • 10 -4 kN/mm. The characteristic lengths and fracture energies of fiber and of matrix are chosen as l f iber c = l matri x c = 2 µm, g f iber c

Fig. 12 .

 12 Fig. 12. Crack patterns under (a) mode II (in-plane shear) and (b) mode III (anti-plane shear) loads, together with the corresponding macroscopic stress-strain curves.

Fig. 13 .

 13 Fig. 13. Geometries of double edge notched tensile specimens: (a) with symmetric crack-like notches, (b) with asymmetric crack-like notches.

Fig. 14 .

 14 Fig. 14. Macroscopic stress-strain curves of the symmetric and asymmetric double edge notched specimens under tension. The letters (A-D) indicate the loading levels at which the crack patterns are given in Fig. 15.

Fig. 15 .

 15 Fig. 15. Crack patterns of symmetric (a) and asymmetric (b) double edge notched specimens at different loading levels. The letters (A-D) refer the loading levels shown in Fig. 14.

Fig. 16 .

 16 Fig. 16. Crack branching test: (a) geometry of the unit cell used in the simulation; (b) stress-strain curve; (c) crack patterns at different loading levels showing crack propagation in soft material, crack branching, crack coalescence and final failure.

Fig. 17 .

 17 Fig. 17. 3D visualization of the unit-cell microstructure of the continuous fiber composite with void. Matrix is hidden for the sake of clarity.

Fig. 18 .

 18 Fig. 18. (a) Macroscopic stress-strain curves of the tensile test on the unit cell of continuous fiber composite, with the red straight line being the extension of the last part (ε > 0.0412) of the stress-strain curve. (b) Numbers of iterations for phase-field and mechanical solutions; the macroscopic stress-strain curve is also plotted to provide a reference of loading level.

Fig. 19 .

 19 Fig.[START_REF] Zhang | A modification of the phase-field model for mixed mode crack propagation in rock-like materials[END_REF]. Crack patterns of the continuous fiber composite under tension in the X -axis at different loading levels (indicated by the numbers on the macroscopic stress-strain curve in Fig.18.a). For the sake of clarity the matrix is hidden and only the voxels with damage variable greater than 0.9 are shown.

ε 2 = ε 22 ε 3 = ε 33 ε 4 = 2ε 12 ε 5 = 2ε 13 ε 6 = 2ε 23 (

 22333412513623

(A. 10 )

 10 Appendix B. Flowchart of the overall algorithm See Fig. B.1.

Fig. B. 1 .

 1 Fig. B.1. Flowchart of the overall algorithm: two fixed-point loops are nested in the staggered scheme, solving the mechanical and the phase-field problems, respectively. CV. check: convergence check for either the mechanical problem or the phase-field problem; ACV: convergence acceleration procedure.

Table 2

 2 Fixed-point algorithm for solving the phase-field problem at a time step t n+1 .

Table 3

 3 Computation time used for the single edge notch simulations with unit cells of different resolutions. All the simulations were undertaken using 8 processors.

	h (mm)	Number of	Total time	Number of	Number of	Approximate time	Approximate time
		voxels	(minutes)	iterations for	iterations for	per iteration for	per iteration for
				PF problem	mechanical	PF problem	mechanical
					problem	(second)	problem (second)
	0.0050	44 823	6	55 656	18 897	0.0015	0.015
	0.0033	100 233	13	64 169	26 947	0.0019	0.024
	0.0025	177 643	41	68 688	32 320	0.0059	0.064
	0.0020	277 053	71	73 025	36 213	0.0090	0.099
	4.2. Single edge notched tensile specimen with inhomogeneous l c or g c

Table 4

 4 Computation time and efficiency for simulations conducted on 1, 4 and 8 nodes (i.e. 28, 112 and 224 cores).

	Number of nodes (Number of cores)	1 (28)	4 (112)	8 (224)
	Wall time	36 h	11.5 h	7 h
	Efficiency	1	0.78	0.65

Acknowledgment

The authors would like to acknowledge the European Union (European Regional Development Fund, ERDF/FEDER), the French state and the Hauts-de-France Region Council for partly funding the ELSAT2020 by CISIT project (POPCOM action). Our gratitude is to Dr. Julien Derouillat from Maison de la Simulation for his contribution to the improvement of AMITEX code. We also would like to thank Prof. James Marrow from Oxford university for the language proofreading of the manuscript.

Appendix A. Projection tensor for strain decomposition

In our implementation, the projection tensors for strain decomposition are constructed as follows. The definition of ε + can be written into matrix form:

where ε is the strain tensor in the principal coordinates. The rotation matrix Q is computed from the eigenvectors (v 1 , v 2 , v 3 ) of the strain tensor, and it is orthogonal, so

Now we replace ⟨x⟩ + = (x + |x|)/2 by a sign matrix s:

The sign matrix s is defined by the eigenvalues (λ 1 , λ 2 , λ 3 ) of the strain tensor: