# Implementation of a phase field damage model with a nonlinear evolution equation in an FFT-based solver

Xiao Ma, Dmytro Vasiukov, Modesar Shakoor, Stepan Lomov, Chung Hae Park

<figure><img src="image_1.png" /></figure>

## To cite this version:

Xiao Ma, Dmytro Vasiukov, Modesar Shakoor, Stepan Lomov, Chung Hae Park. Implementation of a phase field damage model with a nonlinear evolution equation in an FFT-based solver. Engineering Fracture Mechanics, 2023, 290, pp.109518. 10.1016/j.engfracmech.2023.109518. hal-04172865

HAL Id: hal-04172865

https://hal.science/hal-04172865

Submitted on 28 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Implementation of a phase field damage model with a
nonlinear evolution equation in an FFT-based solver

Xiao Ma¹, Dmytro Vasiukov*¹, Modesar Shakoor¹, Stepan Lomov², and
Chung Hae Park¹

¹IMT Nord Europe, Institut Mines-Télécom, Univ. Lille, Centre for
Materials and Processes, F-59000 Lille, France

²Department of Materials Engineering, KU Leuven, Belgium

May 17, 2023

Abstract

This paper focuses on the numerical implementation of phase-field models of
fracture using the Fast Fourier Transform (FFT) based numerical method. Recent
studies on phase-field models focuses on the discussions of the choice of the value
of regularization length, which was proposed to smear the discontinuity of the
sharp crack. Some studies believe that should be considered as a material property
because it have a significant impact on the mechanical behavior of a material in
some phase-field models, for instance, the model proposed by Miehe. However,
our results in this study for heterogeneous materials have shown that the choice of
regularization length not only affects the macroscopic mechanical behavior but also
the local crack propagation patterns. As a result, it can be challenging to select
an appropriate value that produces both accurate macroscopic responses and local
crack patterns for certain phase-field models, such as Miehe's model. Thus, the

<footnote>*Corresponding author: dmytro.vasiukov@imt-nord-europe.fr</footnote>
phase-field model proposed by Wu, which has been proven to reduce the length sensitivity for homogeneous material, has been successfully implemented in an FFT-based solver with the application of Newton-Krylov algorithm.

The length sensitivity of Wu's phase-field model for heterogeneous materials is also addressed in this paper. Our tests show that Wu's phase-field model has partial length sensitivity for heterogeneous materials. However, the main reason for this sensitivity is that one phase enters the damage zone of another phase, which is different from Miehe's model. As a result, a set of criteria for safely choosing the regularization value has been established for Wu's model in this work. Meanwhile, it has also been found that Wu's model may be more suitable for brittle failure due to the introduction of an elastic stage compared to Miehe's model.

Keywords: Phase-field, Damage modeling, Fast Fourier Transform (FFT), Heterogeneous material

# 1 Introduction

A well-known problem with Continuum Damage Mechanics (CDM) models, which are frequently used to simulate the initiation and propagation of cracks, is their dependency on the mesh size. This is because the energy dissipated can vary depending on the size of the element [1]. To address this issue, many methods have been proposed. One technique that has gained popularity due to its mesh independence is the phase-field model [2-4], which is the main focus of this work.

It utilizes a damage variable ($d$) to describe the failure state, where sharp cracks are regularized as diffusive crack bands by a function of $d$ [5, 6] (as illustrated in Fig. 1(a)(b)). The damage variable $d$ is a continuous parameter that varies in the range [0, 1] and describes the degradation state of the material. A value of 0 represents an intact state and a value of 1 represents a fully degraded state. On the one hand, the evolution of the phase-field variable itself completely describes the initiation and propagation of cracks in the material. On the other hand, phase-field models are closely related to the variational approach to brittle fracture [7], which enables a simultaneous approximation
of the displacement field and the cracks by minimizing the total potential energy of the
solid.

<figure><img src="image_3.png" /><figcaption>Figure 1: A solid body Ω with the crack set Γ: (a) conventional sharp cracks and (b) diffuse crack bands in phase-field.</figcaption></figure>

While the phase-field approach is generally considered to be mesh-independent, there is a
condition on element size that must be met. Ref. [2] proved that the regularization length
of the phase-field model must be at least twice the element size of the fractured zone. In
this context, conventional finite element method (FEM) solvers can be cumbersome for
3D simulations and may require efficient parallel implementation [8] and/or automatic
mesh adaption [9] to extend computational limits. This will be even more important
when dealing with heterogeneous materials with complex microstructure cases. Thus, the
Fast Fourier Transform (FFT)-based method has gained popularity due to its intrinsic
parallelization capabilities. It was initially proposed by Moulinec and Suquet [10, 11]
as a voxel-based methodology that does not need stiffness matrix assembling, unlike
conventional FEM solvers.

In the FFT-based method, the local strain tensor is calculated by a convolution product
with a fourth-rank Continuous Green Operator (CGO) and a polarization term. This
convolution product is transformed into a simple multiplication in Fourier space. The
multiplication of the polarization term and the CGO are local operations that can be
easily parallelized [12]. Even if the FFT operation itself is not local, efficient parallel
implementations such as FFTW [13] are available. As a result, large-scale simulations
based on full-resolution images can be performed using the FFT-based method. Since
its introduction, FFT-based method has been improved and used to study a wide range
of physical phenomena in heterogeneous materials [14, 15]. Its high performance and comparisons with FEMs have been clearly demonstrated in Refs. [16–18].

As aforementioned, phase-field models require a minimum ratio between the minimum regularization length and mesh size to ensure mesh independence. Therefore, for complex microstructures, using phase-field methods necessitates a sufficiently refined mesh for better predictions. On the other hand, FFT is more conducive to massive parallelization. As a result, combining phase-field models with the FFT method for fracture problems can be advantageous. Some pioneering work in this area can be found in Refs. [19–22]. The present work also implements phase-field models in an FFT-based environment.

The phase-field model has two important parameters: the critical energy release rate ($G_c$) and the regularization length ($l_c$). $G_c$ is widely recognized as a material damage property, while $l_c$ describes the bandwidth of the diffusive crack, as shown in Fig. 1(b). The smaller $l_c$ is, the narrower the crack band becomes (Fig. 2). If $l_c = 0$, a sharp crack is fully recovered. Recently, there have been discussions in the literature about these parameters. In our latest study, we discovered that omitting certain terms in the phase-field evolution equation when dealing with heterogeneous materials can result in an abnormal damage diffusion between the different material components. Details can be found in Ref. [22], which will not be discussed in this work.

<figure><img src="image_4.png" /><figcaption>Figure 2: Illustration of diffusive damage profile of a crack at $x = 0$ for various length parameters $l_c$.</figcaption></figure>
Meanwhile, there have been discussions about the selection of the $l_c$ value. The primary focus of these discussions is the relationship between $l_c$ and the predicted mechanical behavior of the material. Initially, $l_c$ was intended to be used as a numerical parameter to describe the degree of crack diffusion (crack band-width). However, some studies [23, 24] have found that it can significantly affect the simulated mechanical behavior of materials in certain phase-field models, such as Miehe's phase-field model [2]. As a result, this parameter is now commonly treated as a material parameter [25, 26]. In the literature, there are two main approaches to determining $l_c$ values. One approach involves an analytical calculation of a 1D bar problem with a homogeneous crack phase-field [23], while the other involves an inverse analysis of experimental results [27].

Although these techniques can produce satisfactory results in some cases, they can also complicate the calculation. Additionally, Mandal et al. [28] have demonstrated that these $l_c$ estimates can sometimes be too large compared to the model size, resulting in some unusual phenomena. Our recent research [22] has shown that the impact of $l_c$ on mechanical response is more significant for heterogeneous materials. Furthermore, our results in this study for composite materials have shown that the choice of regularization length not only affects the macroscopic mechanical behavior but also the local crack propagation patterns. As a result, it can be challenging to select an appropriate $l_c$ that produces both accurate macroscopic responses and local crack patterns for certain phase-field models, such as Miehe's model. One of the solutions is to decouple the connection between $l_c$ and mechanical behavior and ensure that $l_c$ only fulfills its original role of describing the degree of crack diffusion. This is one of the main attractions of Wu's phase-field model [29], which introduces traction-separation laws from the cohesive zone model into the phase-field theory to make $l_c$ independent of mechanical response.

Mandal et al. [28] found that Wu's phase-field method demonstrates $l_c$ insensitivity when applied to homogeneous materials. As such, the main goal of this paper is to incorporate Wu's phase-field into the FFT solver and investigate its sensitivity when applied to heterogeneous materials such as composites. A short reminder of Miehe's and Wu's phase-field method will be presented in Sec. 2. The implementation of Wu's phase-field into the FFT solver is presented in Sec. 3. Its results and comparison with Miehe's
phase-field are shown in Sec. 4, and the final conclusions will be in Sec. 5.

# 2 Review of Miehe's and Wu's phase-field models

## 2.1 Regularization of sharp cracks

The starting point of the phase-field method is to smear the sharp crack using a diffusive crack band. The smear function can vary for different phase-field models. Those used in Miehe's (Eq. (1)(a)) and Wu's models (Eq. (1)(b)) are described below:

$$d^M = \exp\left(\frac{-|x|}{l_c}\right), \qquad (1a)$$

$$d^W = 1 - \sin\left(\frac{|x|}{l_c}\right) \qquad (1b)$$

The damage profiles of both models are clearly depicted in Fig. 3. From the figure, it is evident that Wu's phase-field has a narrower diffusion band-width that disappears at $x = \frac{\pi l_c}{2}$, while for Miehe's model this value is equal to $\infty$.

<figure><img src="image_5.png" /><figcaption>Figure 3: The damage profile of different diffusive crack functions [30].</figcaption></figure>

With Eq. (1), it is straightforward to calculate the diffusive crack surface (B) and its density ($\gamma$), where their generic forms can be expressed as follows,

$$B = \int_{-\infty}^{+\infty} \gamma dV, \qquad (2a)$$

$$\gamma = \frac{1}{c_0} [\alpha(d) + l_c (\nabla d)^2], \qquad (2b)$$

$$\alpha(d) = \xi d + (1 - \xi)d^2, \qquad (2c)$$
where $\xi$ and $c_0$ are parameters that vary depending on the specific phase-field model being used. For Miehe's and Wu's models, their values and expressions can be found in Table 1.

<table><caption>Table 1: The parameters of different phase-field models.</caption><thead><tr><th>Model name</th><th>ξ</th><th>α(d)</th><th>c<sub>0</sub></th><th>d(x)</th></tr></thead><tbody><tr><td>Miehe [2, 3]</td><td>0</td><td>d<sup>2</sup></td><td>2</td><td>exp(-|x|/l<sub>c</sub>)</td></tr><tr><td>Wu [4, 29]</td><td>2</td><td>2d - d<sup>2</sup></td><td>π</td><td>1 - sin(|x|/l<sub>c</sub>)</td></tr></tbody></table>

## 2.2 Variational approach

The diffusive crack surface (Eq. (2)) can then be incorporated into the variational approach to brittle fracture, which is based on Griffith's energy principle [31] for a cracked body. In this energetic approach, the damage evolution and crack propagation are a competition between the bulk energy stored in the body and the energy dissipation from the opening of the crack surface. From this point of view, in quasi-static loading, the total energy $\Pi$ can be expressed as:

$$ \Pi := \Phi_s + \Phi_d - P, \qquad (3) $$

where $\Phi_s$ represents the strain energy stored in the cracked body, $\Phi_d$ represents the energy dissipated by opening the crack surface, and $P$ represents the external loading.

In the phase-field method, the dissipated energy $\Phi_d$ can be approximated as the energy required to open the diffusive crack band, as described in Eq. (2). As a result, $\Phi_d$ can be converted from a surface integral to a volume integral. Meanwhile, an energetic degradation function, denoted as $g(d)$, is used to compute the value $\Phi_s$. We consider a domain $\Omega \subset \mathbb{R}^n$ ($n = 1, 2, 3$) that contains a crack set $\Gamma$ with $\Gamma \subset \mathbb{R}^{n-1}$ and is subject to volumetric loading by a body force $f^*$, a displacement $u^*$ on the boundary $\partial\Omega_u$, and a surface force $t^*$ on the complementary boundary $\partial\Omega_t$. The symbol $()^*$ indicates a
prescribed term. In summary, a detailed formulation of the components in Eq. 3 is provided below:

$$ \Phi_s = \int_{\Omega} \varphi(\varepsilon, \Gamma) dV = \int_{\Omega} g(d) \varphi_0(\varepsilon(u)) dV, \quad (4a) $$

$$ \Phi_d = \int_{\Gamma} G_c dS \approx \int_{B} G_c \gamma(d, \nabla d) dV, \quad (4b) $$

$$ P = \int_{\Omega} f^* \cdot u dV + \int_{\partial \Omega_t} t^* \cdot u dS + \int_{\partial \Omega_u} (\sigma \cdot n) \cdot u^* dS. \quad (4c) $$

In Eq. (4)(a), $\varphi$ denotes the elastic strain energy density stored in the cracked body, while $\varepsilon$ and $\varphi_0$ represent the strain tensor and initial strain energy, respectively. In Eq. (4)(b), $G_c$ stands for the critical energy release rate. Typically, the diffused crack domain $B$ is much smaller than the considered domain $\Omega$ with $B \subset \Omega$. In Eq. (4)(c), the symbols $\sigma$ and $n$ correspond to the Cauchy-stress tensor and the normal vector outside the boundary, respectively. Additionally, the symbols $\int dV$ and $\int dS$ represent volume and surface integration, respectively.

The strain and phase-fields $(\varepsilon, d)$ can then be determined by minimizing the total energy under the constraint of irreversibility, which is stated as $\dot{d} \ge 0$, and under the boundary conditions: $u(x) = u^*$ at $\partial \Omega_u$, and $\sigma \cdot n = t^*$ at $\partial \Omega_t$ with $n$ the outward unit normal vector to the external boundary $\partial \Omega$ ($\partial \Omega = \partial \Omega_u \cup \partial \Omega_t$). Eq. (3) can be written in a variational form:

$$ \delta\Pi(u, d) = \int_{\Omega} \sigma \delta\varepsilon dV + \int_{B} \frac{\partial\varphi}{\partial d} \delta dV + \int_{B} G_c \left( \frac{\partial\gamma}{\partial d} \delta d + \frac{\partial\gamma}{\partial \nabla d} \delta \nabla d \right) dV - \int_{\partial \Omega_t} t^* \cdot \delta u dS, \quad (5) $$

where $\sigma = \frac{\partial\varphi}{\partial\varepsilon}$ denotes the Cauchy stress, and $\delta u = 0$ for $\forall x \in \partial\Omega_u$. Eq. (5) can generally be split into two components: $\delta\Pi(u, d)^{P1}$ and $\delta\Pi(u, d)^{P2}$, which represent the mechanical and phase-field parts, respectively. This enables a separate solution for the mechanical and damage fields, and a coupling scheme would then be employed to link these two fields. Further details about the algorithms will be provided in Sec. 3. This section work only focuses the phase-field part $\delta\Pi(u, d)^{P2}$, which can be expressed as:

$$ \delta\Pi(u, d)^{P2} = \int_{B} g'(d) \varphi_0 \delta dV + \int_{B} \frac{G_c}{c_0 l_c} \alpha'(d) \delta dV + \int_{B} \frac{2G_c l_c}{c_0} \nabla d \nabla \delta dV, \quad (6) $$
where − ' means the first derivative. Eq. (6) allows us to easily derive the governing function for the damage field in phase-field models, and its expression is given in Eq. (7). Further details can be found in [32].

$$
g'(d) \mathcal{H} + \frac{G_c}{c_0 l_c} \alpha'(d) - \nabla \cdot \left( \frac{2G_c l_c}{c_0} \nabla d \right) = 0, \quad (7)
$$

Eq. (7) represents the generic form of the governing equation for the phase-field. The regularization parameters, α and c₀, are already known for each phase-field model and are shown in Table. 1. Therefore, the next step is to determine H and g(d) for each phase-field.

2.3 Energetic degradation functions

The energetic degradation function, g ( d ), varies for different phase-field models. For Miehe’s phase-field model, a quadratic energetic degradation function is used, as shown in Eq. (8)(a). In contrast, Wu’s model with a linear softening law, which is used in this work, employs a more complex degradation function as described in Eq. (8)(b).

$$
g(d)^M = (1 - d)^2, \tag{8a}
$$

$$
g(d)^W = \frac{(1-d)^2}{(1-d)^2 + a_1 d(1 - \frac{1}{2}d)}, \quad a_1 = \frac{EG_c}{\sigma_c^2} \cdot \frac{4}{\pi l_c}, \qquad (8b)
$$

where E stands for the Young’s modulus of the material in the case of isotropic material, while σ$_{c}$ represents the material’s strength. In contrast to the gradation function used in Miehe’s model, the formula for a$_{1}$ in Eq. (8)(b) includes σ$_{c}$ . As a result, Wu’s model has an elastic stage.

2.4 History variables

Damage models should reflect the physical properties of materials. A fully damaged material cannot be further damaged and the damage state cannot be reversed. As a result, the boundedness d ∈ [0 , 1] and irreversibility condition ˙ d ≥ 0 must be carefully considered. To ensure the irreversibility condition, the history field H is introduced into Eq. (7). This field represents the maximum initial strain energy from [0 , t$_{n}$ ], where t$_{n}$ is the current time. The two history fields for both Miehe’s (subscript M ) and Wu’s
(subscript W) models are defined as follows:

$$ \mathcal{H}^M (x, t_n) := \max_{0 \le t \le t_n} [\varphi_0^+(x, t)], \qquad (9a) $$

$$ \mathcal{H}^W (x, t_n) := \max_{0 \le t \le t_n} \left[ \varphi_0^+(x, t), \frac{1}{2} \frac{\sigma_c^2}{E} \right], \qquad (9b) $$

where $\varphi_0^+(x, t)$ represents the initial tensile strain energy at time $t$, as the compressive part does not contribute to damage in phase-field models. The $\varphi_0^+(x, t)$ of Wu's model is calculated as:

$$ \varphi_0^+ = \frac{1}{2} \frac{(\langle \bar{\sigma}_1 \rangle_+)^2}{E}, \qquad (10) $$

where $\bar{\sigma}_1$ represents the major principal value of the effective stress $\bar{\sigma}$, which is determined by the elastic relation $\bar{\sigma} = C : \varepsilon$. Here, $C$ refers to the material stiffness tensor without damage. The split operation, denoted by $\langle - \rangle_+$, is defined as $\langle x \rangle_+ = \{x, 0\}$. Since Miehe's model is not the main focus of this work, details about the initial strain energy and split formula can be found in Refs. [3, 19, 22], but will not be presented here.

As shown in Refs. [28, 30, 32], the damage value in Miehe's model inherently satisfies the boundedness $d \in [0, 1]$. However, this is not the case for Wu's model due to its elastic stage. As a result, an additional condition, the Rankine energy $(\frac{1}{2} \frac{\sigma_c^2}{E})$, must be included in the history field $\mathcal{H}$ for Wu's model to ensure that $d = 0$ during the elastic stage.

## 2.5 Phase-field governing equations

According to Eqs. (7) to (10), the governing equations for Miehe's and Wu's model are written as

$$ -2(1-d)\mathcal{H}^M + \frac{G_c}{l_c}d - \nabla \cdot (G_c l_c \nabla d) = 0, \qquad (11a) $$

$$ -\frac{4a_1(1-d)}{(2a_1d - 4d - a_1d^2 + 2d^2 + 2)^2} \mathcal{H}^W + \frac{2G_c}{\pi l_c}(1-d) - \nabla \cdot \left(\frac{2G_c l_c}{\pi} \nabla d\right) = 0, \qquad (11b) $$

where Eq. (11)(a) is specific to Miehe's model while Eq. (11)(b) applies to Wu's model.
### 3 Implementation of Wu's phase-field model into an FFT-based solver

In our previous work [18, 22], for Miehe's model, the derivative operations can be easily converted into simple multiplications. As a result, Eq. (11)(a) is a linear function that can be solved using a fixed-point algorithm. However, Wu's model introduces a new degradation function that introduces non-linearity. This creates a significant challenge when attempting to solve Wu's model using a fixed-point algorithm, as demonstrated in [32].

#### 3.1 An algorithm based on the Newton-Krylov method

Therefore, Newton-Raphson method, one of the most traditional method for solving nonlinear equations, is applied in the present work. This approach is based on the iteration:

$$d^{i+1} = d^i - \frac{R(d)}{R'(d)}, \qquad (12)$$

where the residual function $R(d)$ is

$$R(d) = g'(d)\mathcal{H} + \frac{2G_c}{\pi l_c}(1-d) - \nabla \cdot \left(\frac{2G_c l_c}{\pi}\nabla d\right), \qquad (13)$$

and its derivative

$$R'(d) = \frac{\partial R(d)}{\partial d}. \qquad (14)$$

Note that $\mathcal{H}$ in Eq. (13) is equal to $\mathcal{H}^W$ in Eq. (9)(b). Eq. 12 can be re-written as:

$$R'(d)\delta d = -R(d) \qquad (15)$$

with $\delta d = d_{i+1} - d_i$. Because most terms in $R(d)$ are linear, we only should compute $g''(d)$, where $-$ represents the second derivative. We get

$$\left\{ \begin{aligned} g''(d) &= \frac{a_1 \left(-3(1-d)^2 + a_1(1.5d^2 - 3d + 2)\right)}{\left[(1-d)^2 + a_1 d(1-0.5d)\right]^3} \\ R'(d)\delta d &= g''(d)\mathcal{H}^W \delta d - \frac{2G_c}{\pi l_c}\delta d - \nabla \cdot \left(\frac{2G_c l_c}{\pi}\nabla \delta d\right). \end{aligned} \right. \qquad (16)$$
This linear problem is solved using a Conjugate Gradient (CG) solver, which does not require to introduce a reference material or a polarization term. The non-local term is directly computed using

$$ \nabla \cdot \left( \frac{2G_c l_c}{\pi} \nabla d \right) = F^{-1} \left( J \cdot k \cdot F \left( Q'(x) \cdot F^{-1}(J \cdot k \cdot \hat{d}(x)) \right) \right), \quad (17) $$

where $J = \sqrt{-1}$, $\hat{d}$ is expression of $d$ in Fourier space, $k$ is the frequency vector, $F()$ and $F^{-1}()$ denotes the Fast Fourier transform (FFT) and its inverse.

The FFT and its inverse only need to be computed for the calculation of this non-local term. Ma et al. [18] have shown that Willot's rotated scheme, proposed in Ref. [33], produces fewer numerical artifacts than Moulinec and Suquet's original collocation scheme proposed in [34]. However, we were unable to implement Willot's rotated scheme for the phase-field solver in this work. To avoid using different schemes for the mechanical and phase-field parts, we used Moulinec and Suquet's original collocation scheme for both solvers. This could be an interesting area for future research.

Regarding the CG solver, we need to introduce $LP(d)$ and $RP(\delta d)$ as follows:

$$ \left\{ \begin{aligned} LP(\delta d) &= -\frac{2 G_c}{\pi l_c} \delta d + g''(d) \mathcal{H} \delta d - \nabla \cdot \left( \frac{2 G_c l_c}{\pi} \nabla \delta d \right), \\ RP(d) &= -g'(d) \mathcal{H} - \frac{2 G_c}{\pi l_c} (1-d) + \nabla \cdot \left( \frac{2 G_c l_c}{\pi} \nabla d \right), \end{aligned} \right. \quad (18) $$

where $LP(d)$ in Eq. (18) represents the left-hand side (LHS) of Eq. (15), while $RP(\delta d)$ is the right-hand side (RHS) of that equation. The pseudo-code of the CG solver for Wu's phase-field is presented in Appendix 1. It is a standard CG solver except for operation $f_1$ which is a special treatment to remove the accumulation of floating point error. The convergence criterion of the CG solver (tol) is set as $1.00 \cdot 10^{-6}$.

The output of the CG solver is $\delta d$, so the new damage field is calculated as $d_{new} = d_{old}+\delta d$, and convergence is checked for this new field, and the tolerance of this Newton algorithm is also set to $1.00 \cdot 10^{-6}$. Unlike the conventional convergence criteria used in Refs. [19, 22], our implementation uses the residual $RP(d)$ of the phase-field evolution equation as the convergence criterion for Newton algorithm. This is because $RP(d)$ is derived from the phase-field governing equation, so every value of $d$ is checked against this equation. The algorithm is summarized in pseudo-code in Eq. (19).
Initialization, if $t_n = 0$ :

$$ (a_0) \quad d^{t_n=0}(x) = 0 $$

If $t_n \neq 0$ :

$$ (a_1) \quad d^{t_n} \text{ is known} $$

Time $t_n + 1$ : first input $d_{i=0}^{t_n+1} = d^{t_n}$

While convergence criterion is not met:

$$ (a_2) \quad \text{Call CG solver, with input } d_{i=0}^{t_n+1} $$

$$ (b_2) \quad d_{i+1}^{t_n+1} = d_i^{t_n+1} + \delta d $$

$$ (c_2) \quad \text{Re-calculate } RP(d_{i+1}^{t_n+1}, x) $$

$$ (d_2) \quad \text{if } ||RP||_2 < \text{tolerance } (1.00 \cdot 10^{-6}), \text{ convergence is met,} \tag{19} $$

where $|| - ||_2$ represents the $L_2$ norm. The next step is to input this new $d_{i+1}^{t_n+1}$ into the mechanical part $(\delta\Pi(u, d)^{P1})$. Then the new mechanical fields $(t_{n+1})$ can be calculated. In order to simplify the calculation, this implementation adopts a hybrid formulation where the history variable used for damage evolution Eq. (9) is split, but there is no tension/compression split in the mechanical part. The possibility of applying a hybrid formulation has been confirmed in Ref. [30], and more details can be found in Ref. [32].

Thus, the equations for the mechanical part can be written as:

$$ \left\{ \begin{array}{l} \text{div}(\sigma(u, d)) = 0, \\ \sigma = g(d) \frac{\partial \varphi_0}{\partial \varepsilon}, \\ \varepsilon = \frac{\nabla u + \nabla^T u}{2}. \end{array} \right. \tag{20} $$

## 3.2 The choice of coupling scheme

This workflow can be solved using the same weak coupling approach used for Miehe's model (Refs. [19, 22]), where the damage and strain fields are solved only once per loading step. However, to obtain converged results, this approach requires the use of a sufficiently fine loading step. In this work, to use larger loading steps, we opted for a
stronger coupling and switched to an alternating minimization algorithm.

<figure><img src="image_6.png" /><figcaption>Figure 4: The flow chart of Wu's phase-field method based on the Newton-Krylov algorithm in FFT solver. (remark: CV is damage field convergence check. AltCV is alternating minimization convergence check.)</figcaption></figure>

The flow chart of the alternating minimization is shown in Fig. 4. As depicted in this figure, the mechanical solve is performed first, followed by the phase-field solve, which includes the details of the Newton-Krylov algorithm. Outside these two loops is another loop called the alternating minimization loop, which ensures that all damage and mechanical fields are updated in the current time increment. This loop is controlled by the alternating minimization iteration number << alt >> and the convergence check << AltCV >>, which is based on the change of damage variable (d):

$$e = ||d_{\text{alt}} - d_{\text{alt}-1}||_2, \qquad (21)$$

where || - ||₂ represents the L₂ norm over the model. Due to this convergence check, the damage field computed by the phase-field solver can be used to update the degradation function and return to the mechanical solver. The new strain field would then update the history variable and enter the phase-field solver again. This process can be repeated several times until convergence of the alternating minimization scheme is achieved, thus ensuring a strong coupling between the strain field and the phase-field.
The advantage of this algorithm is that we do not need a very refined loading step to get accurate results. Based on our tests, $\delta\langle\epsilon\rangle = 1.00 \cdot 10^{-4}$ is acceptable, which is 200 times larger than the value used with the weak coupling. Details will be discussed in Sec. 4.

## 3.3 Summary

Generally speaking, this new algorithm involves three nested loops to compute the strain field and phase-field at each increment. The top loop is the alternating minimization, in which we alternately solve the phase-field solver where the degradation function is updated at each iteration and the mechanical solver where the history variable is updated at each iteration. Due to its nonlinearity, the phase-field evolution equation is solved using a Newton-Krylov loop, of which each iteration involves a linear system of equations that is solved using a CG solver. The last loop is the CG solver itself, where at each iteration, the FFT and its inverse are used to compute the non-local term with the correct formulation for heterogeneous materials. Same as the fixed-point algorithm, most operations of this algorithm are still in local that does not require a global assembly hence can be easily parallelized. Considering that the efficiency of FFT has been largely proven in the literature [16], the comparison of computational time between FFT and FEM will not be presented in this work.

# 4 Numerical experiments and analysis

## 4.1 Introduction

The algorithm we propose is designed to solve general 3D problems. However, to reduce computation time in our tests, we used models with only one voxel thickness in the third dimension. Due to the periodic boundary conditions, the out-of-plane strain components are constant in the third direction, which is equivalent to a generalized plane strain condition.

In these analyses, the notch models are subjected to a macroscopic strain along the y-direction with a loading steps of $\delta\langle\epsilon\rangle_{yy} = 1.00 \cdot 10^{-4}$ until final failure, where $\langle - \rangle$
means the average operation over the whole domain. Stress-free conditions are imposed
in the other directions ($\langle\sigma\rangle_{xx} = \langle\sigma\rangle_{zz} = \langle\sigma\rangle_{xy} = \langle\sigma\rangle_{xz} = \langle\sigma\rangle_{yz} = 0$). These mixed-type
loading conditions are applied using the method presented in Ref. [35]. Due to stress-free
boundaries, Poisson's effect is active during loading. For the two-fibers composite models,
the loading direction is changed to $x$ with $\delta\langle\epsilon\rangle_{xx} = 1.00 \cdot 10^{-4}$ and $\langle\sigma\rangle_{yy} = 0$, while the
other conditions remain unchanged.

## 4.2 A single-notch plate

### 4.2.1 Influence of the loading step

As shown in Ref. [19], the fixed-point algorithm used in Miehe's model is highly sensitive to the loading step size and requires a very small increment size. Based on our previous work [22], $\delta\langle\epsilon\rangle = 5.00 \cdot 10^{-7}$ is a good value for Miehe's model and is also used in this work. For Wu's model, thanks to alternating minimization, it is not as sensitive as the fixed-point algorithm, but using too large an increment size can still lead to inaccurate results. Therefore, the first numerical experiment is to find an acceptable increment using a single notch plate model.

This example is a widely used benchmark test to verify phase-field models for brittle fracture [19, 29, 36]. As shown in Fig. 5, it consists of a square plate with a length of 1mm. A straight horizontal notch of 0.5mm is located at the mid-height of the left edge and has a width of one voxel size. Due to the intrinsic periodicity of the FFT solver, lateral bands denoted as M2 (purple color) are added to prevent spurious effects from the boundary conditions. Material M2 is also used for the notch. In order to minimize the effect of the lateral bands, its dimension should be as small as possible. Based on our tests, a dimension of 5 voxels is the minimum value to get convergence in the FFT solver and break the periodicity.

In the single notch model tests, the following material parameters are used for M1: Young's modulus of $E = 2.1 \cdot 10^5$ MPa, Poisson's ratio of $\nu = 0.3$, failure strength of $\sigma_c = 2.45 \cdot 10^3$ MPa, and critical energy release rate of $G_c = 2.7$ N/mm, which are from Ref. [29]. These properties result in an internal length of $l_{ch} = \frac{EG_c}{\sigma_c^2} = 0.095$ mm. Accord-
<figure><img src="image_7.png" /><figcaption>Figure 5: An illustration of the geometry of single notch plate. Material 1 (M1) is in yellow and Material 2 (M2) in purple.</figcaption></figure>

ing to Ref. [30], $a_1 \ge \frac{3}{2}$ is an exigence to ensure the convexity of degradation function of Wu's phase-field model. That means $l_c \le 0.08\text{mm}$ must be satisfied for this material. On the other side, to ensure the non-locality property, Miehe's criterion regarding the ratio between $l_c$ and the voxel size $h$ still needs to be verified for Wu's model. The $l_c/h$ ratio is set to be larger than 5 in some references [37] for Wu's model rather than 2 for Miehe's model. But, our next example will show that 2 is still acceptable for Wu's model. Thus, we have an interval for the choice of $l_c$: $2h < l_c \le 0.08\text{mm}$.

Additionally, the choice of void-like properties (M2) is arbitrary. In this work, they are set as follows: $E_0^e = 21\text{MPa}$, $\nu^e = 0.3$, $\sigma_c^e = 24.45\text{MPa}$, and $G_c^e = 2.7\text{N/mm}$, which result in the same $l_{ch}$ value as for M1.

In this work, three loading steps are used: $\delta\langle\epsilon\rangle = 1.00 \cdot 10^{-4}$, $5.00 \cdot 10^{-5}$ and $2.00 \cdot 10^{-5}$. The model has a resolution of $N = 251$ with a voxel size of $h = 4.00 \cdot 10^{-3}$ mm and the value of $l_c$ is set to $l_c = 0.02\text{mm}$. The averaged stress-strain relations are shown in Figure 6. The figure shows that the averaged stress-strain curves are coincident even when the loading step is decreased by a factor of five. Thus, a loading step of $1.00 \cdot 10^{-4}$ can be considered acceptable.
<figure><img src="image_8.png" /><figcaption>Figure 6: Averaged stress-strain relation curves (remark: the curve thickness of different loading steps are different, which is for distinguishing the curves between each other).</figcaption></figure>

### 4.2.2 Influence of the loading step for high contrast models

In our previous study [22], we discovered that the number of iterations when using an FFT-based solver for Miehe's phase-field evolution equation is sensitive to the contrast in damage properties. As a result, it is crucial to examine the impact of the loading step for Wu's model under varying contrast ratios. The critical energy release rate ($G_c$) of material M1 remain fixed while those of void-like material M2 vary from 2.7N/mm to 2.7.10³N/mm to represent different material contrast ratios. It is important to note that these properties do not represent the actual physics and are only used for numerical studies. Three contrast values are presented in Table 2.

<table><caption>Table 2: The critical energy release rates for different contrast ratios.</caption><tr><td>$G_c^e$ (N/mm)</td><td>2.70</td><td>27</td><td>2700</td></tr><tr><td>Contrast</td><td>1</td><td>10</td><td>1000</td></tr><tr><td>$G_c^e/G_c$</td><td></td><td></td><td></td></tr></table>

The results of contrast 1 are shown in Figure 6. Figure 7 only shows the results of contrasts 10 and 1000. This figure clearly shows that the loading step convergence is not affected by different contrast ratios, which confirms the loading step used in this work:
$\delta\langle\epsilon\rangle = 1.00 \cdot 10^{-4}$. During this study, the authors also observed that the computation time for a high contrast model is significantly longer with same computational hardware. Therefore, further quantitative examination of the number of interactions of Newton-Krylov, CG, and alternating minimization schemes under varying contrasts in damage properties remains an interesting point of study. However, this will be left for future work.

<figure><img src="image_9.png" /><figcaption>Figure 7: Averaged stress-strain relation curves of (a) contrast 10, (b) contrast 1000. (remark: the curve thickness of different loading steps are different, which is for distinguishing the curves between each other).</figcaption></figure>

### 4.2.3 Analysis of the sensitivity to the regularization length

In this part, the material properties stays unchanged as in Sec. 4.1.1, and the contrast ratio $G_c^e/G_c = 1$. Fig. 8 displays stress-strain curves for various $l_c$ values. The FEM results, which lack lateral bands and notch material, are also shown as reference. The resolution of FEM model is set to $N = 1000$. It's important to note that the model and solver differ between FFT and FEM, so their results cannot be compared directly. This figure indicates that there is always sensitivity to $l_c$ for the single notch model. However, there is a clear trend: as $l_c$ decreases, the results approach those of FEM, where the lateral bands are absent. Meanwhile, it should be noted that the residual stress values shown in Fig. 8 for the FFT-based solver may be primarily due to the presence of lateral bands.
<figure><img src="image_10.png" /><figcaption>Figure 8: Results of Wu's phase-field model for single notch model with different $l_c$ values. ($G_c^e/G_c = 1$.)</figcaption></figure>

In this scenario, cracks always initiate at the notch location, which can also be viewed as an interface between two phases. Due to the non-locality criteria, the $l_c$ value cannot be smaller than twice the mesh size. This means that two materials always exist in the damage region of M1, which may contribute to the $l_c$ sensitivity. From a physical perspective, adding M2 material to the single notch model is similar to adding an adhesive or glue layer that strengthens the model's damage behavior. Reducing $l_c$ decreases the area of M2 material considered and weakens this effect, bringing the results closer to those of FEM.

Based on this, it can be inferred that for a given $l_c$ value, the smaller the $G_c^e$ of M2 material, the closer the results are to FEM. When $G_c^e = 0$, a full void notch is obtained. However, due to differences in models and solvers used by FFT and FEM, a direct comparison of their results is not possible. Nonetheless, Fig. 7 can be used to illustrate this trend. By selecting a loading step value of $1.00 \cdot 10^{-4}$ and plotting three contrast results in Fig. 9, a clear trend can be observed that validates our earlier hypothesis.

## 4.3 Two-fibers unit cell composite model

In the next step, a composite unit-cell model with two fibers is analyzed, as shown in Figure 10. Because, cracks will initiate in the matrix if there is no debonding model. The model size is 0.02mm (20µm) with a fiber volume fraction of $V_f = 0.15$ with one voxel size in thickness. The resolution of the model is $N = 225$. The blue part represents the fiber,
<figure><img src="image_11.png" /><figcaption>Figure 9: Stress-strain relation of single notch model with different contrasts.</figcaption></figure>

which is made of E-glass with properties: $E_f = 7.4 \cdot 10^4$ MPa, $\nu_f = 0.2$, $G_f = 60$ N/mm and $\sigma_f^c = 3.40 \cdot 10^3$ MPa. The red part represents the matrix, which is made of Epoxy with properties: $E_m = 4.65 \cdot 10^3$ MPa, $\nu_m = 0.35$, $G_m = 3.00 \cdot 10^{-3}$ N/mm and $\sigma_m^c = 80$ MPa. Additionally, the inter-fiber distance is much smaller than twice the distance from the fiber-matrix interface to the border. Considering the periodic boundary conditions, the strain concentration point should be at the center of the model.

<figure><img src="image_12.png" /><figcaption>Figure 10: Two-fibers unit-cell composite model. (blue part: fiber and red part: matrix.)</figcaption></figure>
4.3.1 Study with a fixed inter-fiber distance

In this study, the value of $l_c$ is varied to study its impact, and an inter-fiber distance ($l_f$) is set to $2.75 \cdot 10^{-3}$mm. The variation of $l_c$ values ranges from $1.6 \cdot 10^{-3}$mm to $2.3 \cdot 10^{-4}$mm.

<figure><img src="image_13.png" /><figcaption>Figure 11: Results of Miehe's phase-field model for two-fibers unit-cell composite model ($l_f = 2.75 \cdot 10^{-3}$mm): (a) Averaged stress-strain curves; (b) Local crack patterns for $l_c = 1.60 \cdot 10^{-3}$mm; (c) Local crack patterns for $l_c = 8.00 \cdot 10^{-4}$mm; (d) Local crack patterns for $l_c = 5.00 \cdot 10^{-4}$mm. (Remark: The orange dash circles are fiber contours.)</figcaption></figure>

Fig. 11 shows the results of Miehe's phase-field model with different $l_c$ values. As shown in the Fig. 11(a), the results of the stress-strain relation confirm the conclusion in Refs. [23, 24]. Increasing $l_c$ reduces the maximum stress and accelerates damage evolution. This is why $l_c$ is linked to the mechanical response. Furthermore, the crack patterns in Fig. 11(d) are in the center zone, while those in Figs. 11(b) and (c) are on the model borders. Sakata et al. [38] have shown that stress/strain concentration can depend on the inter-fiber distance ($l_f$). The smaller the inter-fiber distance, the higher the
concentration. In this case, the inter-fiber distance of the periodic unit-cell implies that strain concentration is not expected to be on the boundaries for crack initiation position. Thus, the crack patterns in Figs. 11(c) and (d) may not be physical. However, as mentioned earlier, $l_c$ is linked to the mechanical response and cannot be freely chosen. Therefore, for Miehe's model, it is difficult to choose an appropriate $l_c$ that produces both proper macro-responses and local crack patterns.

<figure><img src="image_14.png" /><figcaption>Figure 12: Results of Wu's phase-field model for two-fibers unit-cell composite model: (a) Averaged stress-strain curves ($l_f = 2.75 \cdot 10^{-3}$mm); (b) Local crack patterns for $l_c = 8.00 \cdot 10^{-4}$mm; (c) Local crack patterns for $l_c = 5.00 \cdot 10^{-4}$mm; (d) Local crack patterns for $l_c = 2.30 \cdot 10^{-4}$mm. (remark: the curves for $8.00 \cdot 10^{-4}$mm to $l_c = 2.30 \cdot 10^{-4}$mm are superimposed. The orange dash circles are fiber contours.)</figcaption></figure>

In contrast, the results of Wu's model in Fig. 12 show much less sensitivity to the $l_c$ values. In some cases ($l_c$ ranges $8.00 \cdot 10^{-4}$mm to $2.30 \cdot 10^{-4}$mm), Wu's model shows full $l_c$ insensitivity. In Wu's phase-field model, the $l_c$ values and material mechanical response
are uncoupled. However, based on our studies, especially for heterogeneous models, there
is still a criterion for choosing $l_c$. Because heterogeneous models are generally multi-phase
systems, it is necessary to ensure that other phases do not enter the fracture zone of the
damaged phase. If it is the case where $l_c$ value is too large, the macroscopic mechanical
response can be affected (e.g., $l_c = 1.60 \cdot 10^{-3}\text{mm}$). As shown in Fig. 3, the fracture zone
of Wu's phase-field model is limited at a value of $\pi l_c$. Thus, for a composite model, a
criterion $l_c \le \frac{l_f}{\pi}$ should be satisfied.

Another difference between Wu's and Miehe's models was found in this study. Wu's
phase-field model introduces an elastic stage where the damage at a point is always zero
until it exceeds the limit, while in Miehe's model the damage state becomes positive from
the start of the load. This suggests that Wu's model may be more suitable for brittle
failure. An example illustrating this difference is as follows.

Fig. 13(a) shows the damage field of Miehe's model, while Fig. 13(b) shows that of Wu's model. Fig. 13(c) represents the damage profiles of center line for both models. In this figure, both Miehe's and Wu's models are set as $l_c = 5.00 \cdot 10^{-4}$. As can be seen in Fig. 13(c), the strain concentration is in the center zone for both models. But because Miehe's model does not have an elastic part, the damage on the borders is much more significant than in Wu's model, where that region stays intact because the stresses are not beyond the limit.

### 4.3.2 Study with a fixed regularization length

Moreover, if the $\frac{l_f}{l_c}$ ratio is too small, the local crack patterns may not be physical even for Wu's model. To demonstrate this, a new series of simulations have been carried out. In these tests, $l_c$ is kept constant at $l_c = 5.00 \cdot 10^{-4}\text{mm}$ while the inter-fiber distance varies from $2.67 \cdot 10^{-4}\text{mm}$ to $2.75 \cdot 10^{-3}\text{mm}$, which means that the $l_f/l_c$ ratio varies from 0.53 to 5.5.

Fig. 14 shows three representative results. When the ratio is too small, such as 0.53 when $l_f = 2.67 \cdot 10^{-4}\text{mm}$, both mechanical and damage behaviors become unrealistic as shown in Figs. 14 (a) and (b). In this case, the crack should be at the center zone
<figure><img src="image_15.png" /><figcaption>Figure 13: (a) Damage fields of Miehe's model; (b) Damage fields of Wu's models; (c) Damage profiles along the center line. (Remark: $l_f = 2.75 \cdot 10^{-3}$ mm and $l_f = 5 \cdot 10^{-4}$ mm)</figcaption></figure>

and damage should initiate earlier due to the much more significant strain concentration effect compared to a larger inter-fiber distance. However, if the ratio is larger than $\pi$, both mechanical and damage behaviors remain good and show more realistic behavior. For example, in Figs. 14 (a), (c) and (d), the results for $l_f = 1.80 \cdot 10^{-3}$ mm show more concentration effect and earlier damage initiation than $l_f = 2.75 \cdot 10^{-3}$ mm. Additionally, both cracks remain at the center zone which is more realistic according to theory. Therefore, a too small $\frac{l_f}{l_c}$ ratio is unacceptable. However, $l_f$ is considered a geometric property of a model. Therefore, for Wu's model, $l_c$ should be selected as small as possible.

As mentioned earlier, $a_1$ must be greater than or equal to 1.5, which means that $l_c \le \frac{4EG_c}{\pi\sigma_c^2}$ to ensure that Wu's phase-field equation can be solved. All of the above tests provide upper limits for the value of $l_c$. Meanwhile, $l_c$ also has a lower limit to ensure the non-locality property of the phase-field. In Ref. [2], the $l_c/h$ ratio should be larger than 2,
<figure><img src="image_16.png" /><figcaption>Figure 14: Results of Wu's phase-field model for two-fibers unit-cell composite model ($l_c = 5.00 \cdot 10^{-4}$mm): (a) Averaged stress-strain curves; (b) Local crack patterns for $l_f = 2.67 \cdot 10^{-4}$mm; (c) Local crack patterns for $l_f = 1.80 \cdot 10^{-3}$mm; (d) Local crack patterns for $l_f = 2.75 \cdot 10^{-3}$mm. (Remark: The orange dash circles are fiber contours.)</figcaption></figure>

while in Ref. [4] it should be larger than 5. Based on our studies shown in Figure 12, a proper result can also be obtained when $l_c/h = 2.5$. Therefore, this work prefers the proposition in Ref. [2] in terms of the $l_c/h$ ratio criterion. In summary, the proposed criteria for Wu's phase-field are:

$$ \left\{ \begin{array}{l} l_c > 2h, \\ l_c \le \frac{l_f}{\pi}, \\ l_c \le \frac{4EG_c}{\pi\sigma_c^2}. \end{array} \right. \qquad (22) $$

# 5 Conclusion

In the present work, the Wu's phase-field model has been successfully implemented into a FFT-based solver, which is our first main contributions. In Wu's phase-field model,
degradation function $g(d)$ is nonlinear, which makes it quite difficult to solve Wu's phase-field evolution equation using the FFT solver with a fixed-point algorithm. Thus, we developed a new phase-field solver based on the FFT method. We took advantage of this new implementation to consider a strong coupling of the phase-field with the strain field to alleviate restrictions on the loading step.

This new algorithm involves three nested loops to compute the strain field and phase-field at each increment. The top loop is the alternating minimization, in which we alternately solve the phase-field solver where the degradation function is updated at each iteration and the mechanical solver where the history variable is updated at each iteration. Due to its nonlinearity, the phase-field evolution equation is solved using a Newton-Krylov loop, of which each iteration involves a linear system of equations that is solved using a CG solver. The last loop is the CG solver itself, where at each iteration, the FFT and its inverse are used to compute the non-local term with the correct formulation for heterogeneous materials.

Our second contribution is a detailed investigation of the $l_c$ sensitivity of Wu's phase-field model for heterogeneous materials. This is an important topic for phase-field modeling that has not been studied much. Our tests show that Wu's phase-field model has partial $l_c$ sensitivity for heterogeneous materials. But the main reason for this sensitivity is different from Miehe's model. In this respect, Wu's phase-field model has some advantages over Miehe's model.

* In Miehe's model, the main reason for $l_c$ sensitivity is that $l_c$ is coupled with the mechanical behavior. The $l_c$ value also plays a crucial role in forming crack patterns. This makes choosing an appropriate $l_c$ value for Miehe's model very difficult.

* In Wu's model, $l_c$ and mechanical behavior are completely decoupled. The main reason for $l_c$ sensitivity is that one phase enters the damage zone of another phase. As a result, a set of criteria for safely choosing $l_c$ value has been established in this work.

However, if the crack starts at the interface, Wu's model always shows $l_c$ sensitivity. This is mainly because the damage zones always contain multiple phases. But this effect
can be reduced by decreasing the $l_c$ value while satisfying the criteria set in this work. The authors also note that there has been increased attention to interface modeling in phase-field, with proposals for multi-phase-field models or CZMs. These techniques could also be used to reduce $l_c$ sensitivity at the interface, but that is beyond the scope of this paper. The comparison of computational time between Wu's and Miehe's models, which is not included in this work, would still be valuable for the future work.

Another difference between Wu's and Miehe's models was found in this study. Wu's phase-field model introduces an elastic stage where the damage at a point is always zero until it exceeds the limit, while in Miehe's model the damage state becomes positive from the start of the load. This suggests that Wu's model may be more suitable for brittle failure.
# 6 Appendix 1

Initialization: $d(x)$ is already known for previous iteration

(a₀) Define the CG solver tolerance (tol)

(b₀) $i = 0$ (Iteration number of CG solver)

(c₀) $\delta d^{i=0} = 0$

(d₀) Calculate $RP(d)$ and $LP(\delta d^{i=0})$

(e₀) Calculate $r = RP - LP$

(f₀) Define $a = r$

(g₀) Define $ERR_{new} = r^T \cdot r$

While ($ERR_{new} \ge \text{tol}^2$):

(a₁) Calculate $i = i + 1$

(b₁) Calculate $LP(a)$

(c₁) Define $q = LP(a)$

(d₁) Calculate $\alpha = \frac{ERR_{new}}{a^T \cdot q}$

(e₁) Calculate $\delta_d^i = \delta_d^{i-1} + \alpha \cdot a$

(f₁) if $i$ is divisible by 10

Calculate $LP(\delta_d^i)$

Calculate $r = RP - LP$

else

Calculate $r = r - \alpha \cdot q$

end

(g₁) Define $ERR_{old} = ERR_{new}$

(h₁) Calculate $ERR_{new} = r^T \cdot r$

(i₁) Define $\beta = \frac{ERR_{new}}{ERR_{old}}$

(j₁) Calculate $a = r + \beta \cdot a$

End

(23)
References

[1] G. Ernst et al. "Multiscale progressive failure analysis of textile composites". In: *Composites Science and Technology* 70.1 (2010), pp. 61–72. ISSN: 0266-3538. DOI: https://doi.org/10.1016/j.compscitech.2009.09.006.

[2] C. Miehe, F. Welschinger, and M. Hofacker. "Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations". In: *International Journal for Numerical Methods in Engineering* 83.10 (2010), pp. 1273–1311. DOI: https://doi.org/10.1002/nme.2861.

[3] C. Miehe, M. Hofacker, and F. Welschinger. "A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits". In: *Computer Methods in Applied Mechanics and Engineering* 199.45 (2010), pp. 2765–2778. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2010.04.011.

[4] J.-Y. Wu. "A unified phase-field theory for the mechanics of damage and quasi-brittle failure". In: *Journal of the Mechanics and Physics of Solids* 103 (2017), pp. 72–99. ISSN: 0022-5096. DOI: https://doi.org/10.1016/j.jmps.2017.03.015.

[5] B. Bourdin, G. Francfort, and J.-J. Marigo. "Numerical experiments in revisited brittle fracture". In: *Journal of the Mechanics and Physics of Solids* 48.4 (2000), pp. 797–826. ISSN: 0022-5096. DOI: https://doi.org/10.1016/S0022-5096(99)00028-9.

[6] B. Bourdin, G. A. Francfort, and J.-J. Marigo. "The Variational Approach to Fracture". In: *Journal of Elasticity* 2008 91:1 91 (1 Mar. 2008), pp. 5–148. ISSN: 1573-2681. DOI: 10.1007/S10659-007-9107-3.

[7] G. Francfort and J.-J. Marigo. "Revisiting brittle fracture as an energy minimization problem". In: *Journal of the Mechanics and Physics of Solids* 46.8 (1998), pp. 1319–1342. ISSN: 0022-5096. DOI: https://doi.org/10.1016/S0022-5096(98)00034-9.

[8] T. Nguyen et al. "Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro CT imaging". In:
*Extreme Mechanics Letters* 17 (2017), pp. 50–55. ISSN: 2352-4316. DOI: https://doi.org/10.1016/j.eml.2017.09.013.

[9] S. Zhou and X. Zhuang. "Adaptive phase field simulation of quasi-static crack propagation in rocks". In: *Underground Space* 3.3 (2018). Computational Modeling of Fracture in Geotechnical Engineering Part I, pp. 190–205. ISSN: 2467-9674. DOI: https://doi.org/10.1016/j.undsp.2018.04.006.

[10] H. Moulinec and P. Suquet. "A fast numerical method for computing the linear and nonlinear mechanical properties of composites". In: *Comptes rendus de l'Académie des sciences. Série II, Mécanique, physique, chimie, astronomie* 318.11 (1994), pp. 1417–1423. ISSN: 1251-8069.

[11] H. Moulinec and P. Suquet. "A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures". In: *IUTAM Symposium on Microstructure-Property Interactions in Composite Materials*. Springer Netherlands, 1995, pp. 235–246. DOI: 10.1007/978-94-011-0059-5\_20.

[12] Y. Chen et al. "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". In: *International Journal of Solids and Structures* 161 (2019), pp. 111–126. ISSN: 00207683. DOI: 10.1016/j.ijsolstr.2018.11.009.

[13] *FFTW*. http://www.fftw.org.

[14] W. Leclerc et al. "A numerical investigation of effective thermoelastic properties of interconnected alumina/Al composites using FFT and FE approaches". In: *Mechanics of Materials* 92 (2016), pp. 42–57. ISSN: 0167-6636. DOI: https://doi.org/10.1016/j.mechmat.2015.09.002.

[15] B. Wang et al. "Progressive damage analysis of 3D braided composites using FFT-based method". In: *Composite Structures* 192 (2018), pp. 255–263. ISSN: 0263-8223. DOI: https://doi.org/10.1016/j.compstruct.2018.02.040.
[16] P. Eisenlohr et al. "A spectral method solution to crystal elasto-viscoplasticity at finite strains". In: *International Journal of Plasticity* 46 (2013). Microstructure-based Models of Plastic Deformation, pp. 37-53. ISSN: 0749-6419. DOI: https://doi.org/10.1016/j.ijplas.2012.09.012.

[17] S. Lucarini et al. "Adaptation and validation of FFT methods for homogenization of lattice based materials". In: *Computer Methods in Applied Mechanics and Engineering* 388 (2022), p. 114223. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2021.114223.

[18] X. Ma et al. "Numerical artifacts of Fast Fourier Transform solvers for elastic problems of multi-phase materials: their causes and reduction methods". In: *Computational Mechanics 2021* 67:6 67 (6 Apr. 2021), pp. 1661-1683. ISSN: 1432-0924. DOI: 10.1007/S00466-021-02013-5.

[19] Y. Chen et al. "A FFT solver for variational phase-field modeling of brittle fracture". In: *Computer Methods in Applied Mechanics and Engineering* 349 (2019), pp. 167-190. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2019.02.017.

[20] F. Ernesti, M. Schneider, and T. Böhlke. "Fast implicit solvers for phase-field fracture problems on heterogeneous microstructures". In: *Computer Methods in Applied Mechanics and Engineering* 363 (2020), p. 112793. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2019.112793.

[21] Y. Cao et al. "A novel FFT-based phase field model for damage and cracking behavior of heterogeneous materials". In: *International Journal of Plasticity* 133 (2020), p. 102786. ISSN: 0749-6419. DOI: https://doi.org/10.1016/j.ijplas.2020.102786.

[22] X. Ma et al. "Simplified and complete phase-field fracture formulations for heterogeneous materials and their solution using a Fast Fourier Transform based numerical method". In: *Engineering Fracture Mechanics* 279 (2023), p. 109049. ISSN: 0013-7944. DOI: https://doi.org/10.1016/j.engfracmech.2023.109049.
[23] T. T. Nguyen et al. "On the choice of parameters in the phase field method for simulating crack initiation with experimental validation". In: *International Journal of Fracture* 197 (2 Feb. 2016), pp. 213-226. ISSN: 15732673. DOI: 10.1007/S10704-016-0082-1/METRICS.

[24] X. Zhang et al. "Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale". In: *Computational Mechanics* 2017 59:5 59 (5 Jan. 2017), pp. 737-752. ISSN: 1432-0924. DOI: 10.1007/S00466-017-1373-8.

[25] J. Espadas-Escalante, N. van Dijk, and P. Isaksson. "A phase-field model for strength and fracture analyses of fiber-reinforced composites". In: *Composites Science and Technology* 174 (2019), pp. 58-67. ISSN: 0266-3538. DOI: https://doi.org/10.1016/j.compscitech.2018.10.031.

[26] T. Guillén-Hernández et al. "A micromechanical analysis of inter-fiber failure in long reinforced composites based on the phase field approach of fracture combined with the cohesive zone model". In: *International Journal of Fracture* 220 (2 Dec. 2019), pp. 181-203. ISSN: 15732673. DOI: 10.1007/S10704-019-00384-8/METRICS.

[27] T. Nguyen et al. "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". In: *Journal of the Mechanics and Physics of Solids* 95 (2016), pp. 320-350. ISSN: 0022-5096. DOI: https://doi.org/10.1016/j.jmps.2016.06.004.

[28] T. K. Mandal, V. P. Nguyen, and J.-Y. Wu. "Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture". In: *Engineering Fracture Mechanics* 217 (2019), p. 106532. ISSN: 0013-7944. DOI: https://doi.org/10.1016/j.engfracmech.2019.106532.

[29] J.-Y. Wu and V. P. Nguyen. "A length scale insensitive phase-field damage model for brittle fracture". In: *Journal of the Mechanics and Physics of Solids* 119 (2018), pp. 20-42. ISSN: 0022-5096. DOI: https://doi.org/10.1016/j.jmps.2018.06.006.
[30] J.-Y. Wu et al. "Phase-field modeling of fracture". In: *Advances in Applied Mechanics*. Ed. by S. P. Bordas and D. S. Balint. Vol. 53. Elsevier, 2020, pp. 1-183.
DOI: https://doi.org/10.1016/bs.aams.2019.08.001.

[31] A. Griffith. "The phenomena of rupture and flow in solids". In: *Philosophical Transactions of the Royal Society of London* 221 (1920), pp. 163-198.

[32] X. Ma. "The elastic and damage modeling of heterogeneous materials based on the Fast Fourier Transform". PhD thesis. IMT Nord Europe - University of Lille, 2022.

[33] F. Willot. "Fourier-based schemes for computing the mechanical response of composites with accurate local fields". In: *Comptes Rendus - Mecanique* 343.3 (2015), pp. 232-245. ISSN: 16310721. DOI: 10.1016/j.crme.2014.12.005.

[34] H. Moulinec and P. Suquet. "A numerical method for computing the overall response of nonlinear composites with complex microstructure". In: *Computer Methods in Applied Mechanics and Engineering* 157.1-2 (Apr. 1998), pp. 69-94. ISSN: 00457825. DOI: 10.1016/S0045-7825(97)00218-1.

[35] M. Kabel, S. Fliegener, and M. Schneider. "Mixed boundary conditions for FFT-based homogenization at finite strains". In: *Computational Mechanics* 57.2 (2016), pp. 193-210. ISSN: 1432-0924. DOI: 10.1007/s00466-015-1227-1.

[36] M. Ambati, T. Gerasimov, and L. De Lorenzis. "A review on phase-field models of brittle fracture and a new fast hybrid formulation". In: *Computational Mechanics* 55 (2015), pp. 383-405. DOI: 10.1007/s00466-014-1109-y.

[37] P. Zhang et al. "Modelling progressive failure in multi-phase materials using a phase field method". In: *Engineering Fracture Mechanics* 209 (2019), pp. 105-124. ISSN: 0013-7944. DOI: https://doi.org/10.1016/j.engfracmech.2019.01.021.

[38] S. I. Sakata and T. Sakamoto. "A Local Sensitivity-Based Multiscale Stochastic Stress Analysis of a Unidirectional Fiber-Reinforced Composite Material Considering Random Location Variation of Multifibers". In: *ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part B: Mechanical Engineering* 5 (3 Sept. 2019). ISSN: 23329025. DOI: 10.1115/1.4043400.