Ye Lu
email: ye.lu@cea.fr

Thomas Helfer
email: thomas.helfer@cea.fr

Benoît Bary
email: benoit.bary@cea.fr

Olivier Fandeur
email: olivier.fandeur@cea.fr

An efficient and robust staggered algorithm applied to the quasi-static description of brittle fracture by a phase-field approach

Keywords: Phase field, Brittle fracture, Discontinuous crack propagation, Semi-implicit algorithm, Irreversibility implementation

The phase field method has been widely adopted in brittle fracture analysis for its ability to handle complex crack topology. This paper presents a novel efficient and robust phase field algorithm for quasi-static brittle fracture analysis. This algorithm overcomes two major issues that affect significantly the numerical cost of the method: the treatment of unstablediscontinuous crack propagation and the inequality constraint associated with the irreversibility of the damage evolution. To handle unstablediscontinuous crack propagation, a semi-implicit scheme, which combines the usual explicit and implicit schemes, is proposed. Different from explicit schemes that require small time steps and purely implicit schemes that looselose immediately efficiency when encountering unstable crackdiscontinuous propagation, the proposed method can releasealleviate the steps constraint while keeping a good robustness with discontinuous cracking. Concerning the irreversibility constraint, this work proposes a practical and easy-to-implement method. It is shown that this method is extremely efficient and robust without any supplementary numerical coefficient. The efficiency of the method is demonstrated by means of representative numerical examples.

Introduction

Numerical simulations of fracture processes play an important role in engineering designs. The prediction of crack propagation in brittle materials is usually based on Griffith's theory [START_REF] Griffith | Vi. the phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london[END_REF], which compares the energy released by the crack propagation to a critical energy release rate. A general concept in these models is that a necessary condition for crack propagation is that the energy release rate reaches a critical value. Discontinuous crack propagation happens when the elastic energy stored in the body is greater than this value.

One of the major difficulties in fracture simulations is the intrinsic spatial singularity, i.e. the displacement field is discontinuous across the crack point the crack tip. The numerical treatment of this singularity in finite element (FE) models can be considered, either by embedding discontinuity lines by means of remeshing strategies [START_REF] Yang | A 3d h-adaptive local remeshing technique for simulating the initiation and propagation of cracks in ductile materials[END_REF][START_REF] Gibert | A 3d automatic mesh refinement x-fem approach for fatigue crack propagation[END_REF] and cohesive elements [START_REF] Elices | The cohesive zone model: advantages, limitations and challenges[END_REF][START_REF] Blal | Artificial compliance inherent to the intrinsic cohesive zone models: criteria and application to planar meshes[END_REF],

or by enriching the displacement field with discontinuities using the partition of unity method [START_REF] Babuška | The partition of unity method[END_REF], as introduced in extended finite element methods (XFEM) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. However, tracing the evolution of complex crack paths, including crack initiation, propagation, merging and branching, in such models has proven to be a tedious task, especially in three dimensional cases [START_REF] Moës | Non-planar 3d crack growth by the extended finite element and level sets-part i: Mechanical model[END_REF][START_REF] Gravouil | Non-planar 3d crack growth by the extended finite element and level sets-part ii: Level set update[END_REF].

An alternative approach to deal with the discontinuous crack topology consists in incorporating a smooth auxiliary (phase) field. This auxiliary variable describes continuously the transition between fully broken and intact material phases within a small band. As a consequence, discontinuities are not directly introduced into the model, the evolution of fracture surfaces is provided by the problem solution on a fixed mesh. This is particularly advantageous for handling complex crack paths. Following the terminology of [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], such regularized models are referred to as phase field models.

Phase field models of quasi-static brittle fracture are based on the variational formulation of Griffith's-type methods, proposed by Francfort and Marigo [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and implemented for the first time by Bourdin et al. [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]. The energy functional in the phase field model resembles closely the potential used in image segmentation [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF], which relies on the regularization concept through Γ-convergence [START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via t-convergence[END_REF]. More recently, Miehe et al. [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field fe implementations[END_REF] presented a quasi-static phase field formulation aligned with thermodynamic arguments. This model possesses several practical features for its numerical implementation and has been widely adopted for various applications (see e.g. [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] Xia | Topology optimization for maximizing the fracture resistance of quasi-brittle composites[END_REF]).

A higher-order extension within the isogeometric analysis framework can be found in [START_REF] Borden | A higherorder phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework[END_REF].

The major limit of the phase field method for engineering-size applications is its expensive computational cost. Phase field models usually require very fine meshes around crack paths for convergence reasons. Although parallelization [START_REF] Ziaei-Rad | Massive parallelization of the phase field formulation for crack propagation with time adaptivity[END_REF] or adaptive meshing strategies [START_REF] Yulong | Adaptive consistent element-free galerkin method for phase-field model of brittle fracture[END_REF][START_REF] Burke | An adaptive finite element approximation of a variational model of brittle fracture[END_REF][START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF][START_REF] Li | Variational h-adaption method for the phase field approach to fracture[END_REF][START_REF] Patil | An adaptive multiscale phase field method for brittle fracture[END_REF][START_REF] Zhang | Moving mesh finite element simulation for phase-field modeling of brittle fracture and convergence of newton's iteration[END_REF] can be applied to accelerate the simulations, the development of efficient solution schemes remains a key point to make the method more attractive for problems having real engineering interests.

Depending on whether the displacement and the crack field are computed simultaneously or alternativelyalternately, two families of phase field solution schemes can be distinguished: monolithic and staggered schemes. Monolithic schemes are expected to have higher convergence rates, since the two field solutions are solved simultaneously in a unique Newton-Raphson loop.

However, this type of approaches suffers from the well-known discontinuous crack propagation in quasi-static simulations. The loss of convexity of the problem makes the solutions difficult to converge. Additional numerical treatments are usually required to improve its robustness (see e.g. [START_REF] Wick | Modified newton methods for solving fully monolithic phasefield quasi-static brittle fracture propagation[END_REF]). Furthermore, this approach cannot be directly implemented into a commercial code without supplementary developments.

Staggered schemes are more robust and easier to implement in an existing commercial code [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]. They are based on an operator split algorithm in which the mechanical problem and the damage problem are solved alternativelyalternately. Explicit or implicit strategies can then both be considered. Explicit schemes [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][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] consider that the damage is constant over an incremental time step when computing the displacement field. The damage field is updated for the next time step after the convergence of the mechanical problem. They were found to be extremely robust, even with discontinuous crack propagation. Nevertheless, this kind of approaches usually requires very small time steps, which are not accessible for many engineering applications. Implicit schemes can overcomealleviate this time step dependency (see e.g. [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF][START_REF] Helfer | Modélisation par champ de phase de la fissuration des matériaux fragiles: Aspects numériques et applications au combustible nucléaire oxyde, 13ème colloque national en calcul des structures[END_REF]). The implicit method proposed by [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF] is usually known as an alternate minimization scheme, in which two minimization problems for the displacement and damage fields are solved independently by fixing one of these two fields. Similarly to explicit schemes, this method computes the displacement field at constant damage, but instead updates the latter at current time step and then computes again a new displacement. Therefore, different from the explicit scheme, the convergence of damage field can be ensured for each time step in the alternate minimization procedure. This method is very stable but usually converges extremely slowly, although relaxation methods can be employed to accelerate it [START_REF] Farrell | Linear and nonlinear solvers for variational phasefield models of brittle fracture[END_REF]. Experience [START_REF] Helfer | Modélisation par champ de phase de la fissuration des matériaux fragiles: Aspects numériques et applications au combustible nucléaire oxyde, 13ème colloque national en calcul des structures[END_REF] shows that the implicit methods work independently ofwith relatively large time steps but usually require numerous iterations to find the converged solutions, particularly when encountering unstablediscontinuous crack propagation.

Nevertheless, this kind of approaches is appealing for many engineering problems subjected to long-term loading (e.g. nuclear fuel simulations [START_REF] Helfer | Modélisation par champ de phase de la fissuration des matériaux fragiles: Aspects numériques et applications au combustible nucléaire oxyde, 13ème colloque national en calcul des structures[END_REF]), in which small time steps are not acceptable. This work falls within this perspective. In order to have a solution scheme independent ofallowing large time steps while keeping a good efficiency and robustness for unstable propagationdiscontinuous cracking, we present here a novel staggered scheme, referred to as semi-implicit method, for quasi-static brittle fracture simulations. The first important ingredient of the semi-implicit method relies on an alternativea modified one-loop implicit scheme [START_REF] Helfer | Modélisation par champ de phase de la fissuration des matériaux fragiles: Aspects numériques et applications au combustible nucléaire oxyde, 13ème colloque national en calcul des structures[END_REF], in which we suggest to integrate the damage update into the mechanical Newton loop so as to improve the damage convergence rate.

Since the damage field is updated after computing each displacement increment estimate, the convergence of the mechanical problem may be strongly perturbed but we can expect a faster damage convergence, compared to the alternate minimization scheme [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF]. Indeed, as shown in the examples, damage field can converge relatively faster than displacement field, and even small variations of the damage field can affect significantly the mechanical equilibrium. Particularly in unstablediscontinuous propagation steps, globally stable damage field cannot ensure immediately that equilibrium.

We observed that this problem causes a significant number of iterations only needed for achieving that mechanical convergence, after the convergence of damage field. A simple way for overcoming this problem is to switch to a less strict convergence criterion for the mechanical problem. However, this is not considered in this work. Instead, on the level of solution scheme, this work proposes to apply an explicit-type resolution for the mechanical convergence by fixing the damage variables, once the damage convergence is found by the suggested implicit procedure. Therefore, this semi-implicit scheme, which combines the implicit and explicit schemes, presents a twofold advantage: large time steps independence and robustness for discontinuous cracking.

Another important issue in phase field implementation to deal with is the irreversibility condition. The thermodynamics arguments demand that the crack field should not be reversible in any case. This results in an inequalityconstrained phase field problem. Different methods to impose the irreversibility can significantly affect the computational cost of phase field solutions.

The simplest way, proposed by Miehe et al. [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], is to introduce a monotonically increasing history field energy function in the phase field equation to replace the original loading-induced reference energy. This method seems very efficient, since no explicit constraints are introduced in the phase field equation. However, the introduction of the history field makes the phase field solution differ from the original variational framework and its equivalence to the original minimization problem can not be proven [START_REF] Gerasimov | On penalization in variational phase-field models of brittle fracture[END_REF]. From this viewpoint, the penalty methods [START_REF] Gerasimov | On penalization in variational phase-field models of brittle fracture[END_REF], which rely on equality-based formulations and can keep the variational nature of the original problem, seem an appealing option. The main drawback of these methods lies on the introduction of the penalty coefficient, which may cause ill-conditioning of the problem.

Although analytic derivations of the 'optimal' coefficients have been proposed, the choice of these penalty coefficients remains a delicate task. The most standard method to solve the constrained phase field problem is the Lagrangian method or its augmented version [START_REF] Wheeler | An augmented-lagrangian method for the phase-field approach for pressurized fractures[END_REF]. This kind of methods may result in a large size system and therefore an extremely high numerical cost, with the introduction of a large amount of extra-variables. Hence, this kind of methods needs an efficient implementation. This work proposes an efficient way to impose the irreversibility condition, which can be viewed as an efficient variant of Lagrangian method for the inequality-constrained phase field problem. Hence, this method can keep the original variational nature of the phase field solution. The implementation relies on an iterative procedure with only equality constraints on a reduced and irreversible active subset, which limits the size of the augmented system.

Particularly, a vanishing energy driving force is applied to the constrained subset, which ensures the positivity of Lagrange multipliers. It is shown that the proposed method is very efficient and can lead to a similar result to that of Lagrangian method. Together with the proposed semi-implicit solution scheme, the novel phase field algorithmic implementation seems very robust for quasi-static brittle fracture simulations. Several representative numerical examples will be presented for demonstrating the efficiency and robustness of the method. This paper is organized as follows. Section 2 introduces the phase field formulation for brittle fracture, in which the novel method for imposing the irreversibility condition is described. Then, the semi-implicit solution scheme is exposed in Section 3. Section 4 will present some numerical examples.

Finally, the paper will be closed by some concluding remarks.

Problem formulation

Phase field modeling of crack topology

Let us consider a n dimensional domain Ω ⊂ R n with a fully open crack Γ ⊂ R n-1 . The crack field can be then described by a damage function d(x) which is equal to 0 everywhere, except for Γ. For a fully damaged point, we consider that d equals to 1. Figure 1(a) shows a sharp crack topology. Diffuse approximation consists in introducing a continuous transition between the crack and undamaged zone. This leads to a partially damaged domain around the crack, in which d ∈ (0, 1) (see e.g. Figure 1(b)). In phase field methods [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], the width of the damaged zone is controlled by an internal characteristic length l c (called also regularization length) and the crack surface density can be defined as As a consequence, the total crack surface is approximated by its integral over the domain

γ(d, ∇d) = 1 2l c (d 2 + l 2 c |∇d| 2 ) (1)
Ω Γ d = 0 d = 1, ∀x ∈ Γ (a) Discrete crack Ω d > 0 d = 0 d = 1, ∀x ∈ Γ (b) Diffuse crack
Γ ≈ Γ(d) = Ω γ(d, ∇d)dV = Ω 1 2l c (d 2 + l 2 c |∇d| 2 )dV (2)

Variational formulation and governing equations

The variational formulation of brittle fracture is given by the work of Francfort and Marigo [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. Let Ω denote again the elastic body with a crack Γ, a prescribed surface loading t is applied on its boundary ∂Ω, the energy functional in quasi-static setting reads

Π(u, Γ) = Ω ψ e (ε e (u))dV Πe + Γ G c dΓ Π d - ∂Ω tudA Πext (3)
where u denotes the displacement field, G c is the critical release energy per unit crack surface, and Π e , Π d , Π ext denote respectively the elastic bulk energy, dissipation potential and external work. In this work, no external applied forces are taken into account. Hence, we can simply consider t = 0 without changing the formulation. According to Griffith's theory, cracks should propagate along a path of the least energy and satisfy the irreversibility condition: Γ(s) ⊆ Γ(t), ∀s < t, which leads to a minimization problem for crack propagation

(u, Γ) = Arg inf u * ,Γ * Π(u * , Γ * ) (4)
We remark that the volume body forces are neglected here. Under infinitesimal strain assumption, the elastic strain tensor is defined as

ε e = ε = ∇ s u (5)
where ∇ s denotes the symmetric gradient operator.

The dissipation functional is defined as the work needed to create the corresponding crack surface Γ which is approximated by Γ(d) in diffuse approximation, i.e.

Π d = Γ G c dΓ ≈ Γ G c γ(d, ∇d)dΓ (6)
In order to be consistent with the second thermodynamic principle, we demand a positive crack dissipation rate: Πd ≥ 0. This results in

δ d γ ≥ 0 and ḋ ≥ 0 (7)
where δ d γ = 1 lc (d -l 2 c ∆d). Note that this positive increment constraint [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] ensures the global irreversibility condition of Γ. In addition, the first condition in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] is automatically satisfied in phase field approximation.

Taking into account the rigidity degradation of material due to the fracture and a spectral decomposition of strain tensor [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] for preventing the crack opening in compression, we can define the elastic energy functional as

Π e = Ω g(d)ψ + e (ε) + ψ - e (ε) dV (8)
where the degradation function is chosen as: g(d) = (1 -d) 2 + k, with a small value k appearing to avoid singular problems. In our work, this value is set to 10 -10 . The spectral decomposition of strain is written as

ε = ε + + ε -with ε ± = i=1 ε i ± n i ⊗ n i (9)
where (•) + and (•) -denote respectively the tensile and compressive modes.

ε i and n i are the eigenvalues and eigenvectors in different dimensions. The Macaulay brackets are defined as:

• ± = 1 2 (•±|•|).
With the Lamé constants λ and µ, the elastic energy relating to each part of strain is defined as

ψ ± e (ε) = λ 2 ( tr[ε] ± ) 2 + µε ± : ε ± (10)
Therefore, considering that Π e = Ω Ψ e dV = Ω (σ : ε -f d)dV , the definition of stress is given by

σ = ∂Ψ e ∂ε = g(d) λ tr[ε] + I + 2µε + + λ tr[ε] -I + 2µε - (11)
and the energetic driving force of damage reads

f = - ∂Ψ e ∂d = 2(1 -d)ψ + e ( 12
)
Taking the variation of the total energy functional, the optimization problem (4) with inequality constraint becomes [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]:

div σ = 0 in Ω σ • n = t on ∂ σ Ω (13)
with a Kuhn-Tucker-type (KT) condition

ḋ ≥ 0, f -G c δ d γ ≤ 0, ḋ(f -G c δ d γ) = 0 (14)
in the whole domain Ω and ∇d • n = 0 on the boundary ∂Ω. Note that the equation ( 14) 3 allows to compute the damage d when ḋ > 0, i.e.

f -G c δ d γ = 0 → 2(1 -d)ψ + e = G c l c d -l 2 c ∆d for ḋ > 0 (15)

Finite element discretization

The implementation of phase field method can be performed in different manners (monolithic or staggered), as mentioned in introduction. For the reason of stability and robustness, this work is based on a operator split scheme, i.e. staggered implementation [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF].

The finite element method is used for the spatial discretization.

u(x) = N u (x)U, d(x) = N d (x)d (16)
where U and d are respectively the usual nodal displacement and damage vectors. N u and N d are corresponding FE shape functions. Noting by B the gradient of shape functions, the FE discretization leads to the two following equations ( 17) et ( 18):

R u (u, d, t) = ∂Π ∂u = ∂Π e ∂u - ∂Π ext ∂u = F int (u, d) -F ext (t) (17)
where

F int = Ω B T u σdV , F ext = ∂Ω N T u t dA. R d (u, d) = ∂Π ∂d = ∂Π e ∂d + ∂Π d ∂d = K d (u)d -F d (u) (18)
where

K d = Ω (2ψ + e + Gc lc )N T d N d dV + Ω G c l c B T d B d dV , F d = Ω 2N T d ψ + e dV .
The solution (u, d) to the original variational problem (4) should be computed through minimizing the residuals R u and R d under the irreversibility condition ḋ ≥ 0. Due to the splitting of the energy [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], the mechanical problem (17) remains nonlinear even when fixing the damage variable at a constant stage. The standard Newton-Raphson method is usually applied for the mechanical problem.

Irreversibility condition

As the reference elastic energy ψ + e can decrease with external loading, the desired condition: ḋ ≥ 0 in Ω is not automatically satisfied in the method.

Hence, this point needs some additional treatment. To this end, different methods exist in the literature.

Brief review of previous methods

The most thorough way from the mathematical point of view is to apply an inequality constraint when solving the phase field equation

2(1 -d)ψ + e = G c l c d -l 2 c ∆d subject to ḋ ≥ 0 (19)
As mentioned in the introduction, this inequality constraint can be enforced using the Lagrange multiplier (see e.g. [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Wheeler | An augmented-lagrangian method for the phase-field approach for pressurized fractures[END_REF]). However, the standard Lagrangian method results in a large size augmented system, which is extremely expensive to solve. Its application to an engineering size problem is inappropriate. An alternative option is to use the penalty method [START_REF] Gerasimov | On penalization in variational phase-field models of brittle fracture[END_REF] which does not introduce extra-variables in the problem. However, the introduction of penalty coefficients may result in ill-conditioned systems.

A more efficient way [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] is to solve the phase field equation ( 18) with a constraint on the damage variation only when a point reaches a critical value, i.e.

ḋ = 0 for d ≈ 1 (20)
Therefore, this method only prevents the cure of fully broken points. The irreversibility condition is not ensured in the whole domain Ω.

An alternative efficient method [START_REF] Lancioni | The variational approach to fracture mechanics. a practical application to the french panthéon in paris[END_REF] consists in solving, at first, only the unconstrained phase field equation [START_REF] Borden | A higherorder phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework[END_REF] for the whole domain Ω, and then doing an a posteriori projection, which reads

2(1 -d * )ψ + e = Gc lc (d * -l 2 c ∆d * ) d(t) = max(d * , d(s)), ∀s < t (21)
This method enforces the exact irreversibility condition, but the projection makes the solution deviate from the original variational framework.

More recently, Miehe et al. [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] proposed to introduce a compact history field function H to replace the reference elastic energy ψ + e in phase field equation, which results in

H = max (ψ + e (t), ψ + e (s)) , ∀s < t 2(1 -d)H = Gc lc (d -l 2 c ∆d) (22)
In this method, the increasing of damage is driven by the history field function

H. The following KT condition is expected to be satisfied in Ω ḋ ≥ 0, ψ + e -H ≤ 0, ḋ(ψ + e -H) = 0 (23)
As commented by [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], the damage will not decrease with the monotonically increasing function H, therefore the irreversibility ḋ ≥ 0 can be always ensured. However, nothing can guarantee the consistency condition [START_REF] Li | Variational h-adaption method for the phase field approach to fracture[END_REF] 3 : the equation ( 22) does not guarantee that the damage will not increase with a decreasing reference energy ψ + e . Indeed, the introduction of the history field function, which can be viewed as an a priori projection, makes the solution lose its original variational nature. Regardless of this drawback, this method has been widely adopted for its implementation simplicity and efficiency.

An efficient irreversibility implementation

In order to keep at maximum the variational nature of the phase field solutions without losing much efficiency, we propose here a novel implementation of the inequality condition. Let us consider a time discrete formulation of the constrained damage problem [START_REF] Ziaei-Rad | Massive parallelization of the phase field formulation for crack propagation with time adaptivity[END_REF], the damage at instant t is given by

2(1 -d t )ψ + e = G c l c d t -l 2 c ∆d t subject to d t ≥ d t-1 (24)
The use of standard Lagrangian method leads to the following equations with equality constraints on a set of points

     2(1 -d t )ψ + e + γ = Gc lc (d t -l 2 c ∆d t ) γ(d t -d t-1 ) = 0 γ ≥ 0 (25)
where γ denotes the Lagrange multiplier, which should not be negative: γ ≥ 0. The solution of the above equations can be found by an iterative procedure, in which the values of Lagrange multipliers are updated according to the current estimate and successively until all the constrained points have positive damage increments. This procedure is extremely expensive for two reasons: introduction of a large number of additional variables (which is equal to the degree number of original system) and numerous iterations required for finding the a priori unknown values of Lagrange multipliers. Although active-set algorithms (e.g. [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]) can be used to reduce the additional variables, the iterations number may remain important.

In this work, we propose an efficient solution procedure with only equality constrains on a subset of the global system, which is based on the following modified formulation

         2(1 -d t )(1 -p)ψ + e + pγ = Gc lc (d t -l 2 c ∆d t ) γ(d t -d t-1 ) = 0 γ ≥ 0 γ(1 -p) = 0 (26)
where p is an indicator function which equals to 1 for the constrained subset and equals to 0 for the unconstrained subset. The first equation implies that the energetic driving force f = ∂Ψe ∂d [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] should be vanishing for the constrained points. This supplementary condition: f = 0 (ψ + e = 0) has a minor effect on those unconstrained points. Hence, the final phase field solution is expected to keep the original variational nature. This formulation enables the following solution procedure:

• Given a precomputed energy ψ + e and the previous damage d t-1

• Initiation: p(x) = 0, ∀x ∈ Ω

• For iteration i

* = ∅ or d * i -d * i-1 ∞ ≤ c , then d t = d * i , END. Otherwise, repeat the iteration.
This procedure is computationally efficient, since only equality constraints are imposed on a subset of the global system. In addition, it should be noticed that the constrained set is irreversibly updated in the sense that one already constrained point cannot become unconstrained in a solution loop. The feasibility of this concept is ensured by the vanishing energetic driving force, which guarantees the positivity of the Lagrange multiplier.

This way, the iterative procedure usually converges very fast within several iterations. The irreversibility condition: ḋ ≥ 0 is automatically satisfied with the converged solution. Compared to those projection-based approaches (e.g. [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Lancioni | The variational approach to fracture mechanics. a practical application to the french panthéon in paris[END_REF]), this method is more expensive but is expected to lead to more robust variationally-consistent phase field solutions.

The proposed implementation differs from active-set algorithms in the following aspects. The constrained set is increasingly updated for a given energy state:

D i ⊂ D i+1
, where i stands for the iteration step. However, in active-set algorithms, active constrained points may be removed from the constrained set (e.g. according to the possitivity of the Lagrange multiplier), hence D i ⊂ D i+1 . In addition, the energy is removed for the constrained set in proposed algorithm, which is not the case for active-set algorithms. The impact of this irreversibility implementation is numerically analyzed by a 1D example in Appendix.

3. A robust solution scheme for unstable crack propagation

UnstableDiscontinuous crack growth

In quasi-static crack modeling, unstable crack propagation is present as a discontinuous crack evolution in timediscontinuous crack evolution happens as the jump of the crack field. It remains a computationally challenging problem with standard iterative Newton algorithms, since the jump of the solution (e.g. displacement field) is difficult to capture with the tangent operator due to the loss of convexity of the problem. One possible way to overcome this issue is to switch to truly dynamic simulations [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Larsen | Existence of solutions to a regularized model of dynamic fracture[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], in order to capture the loss of kinetic energy which is not taken into account in quasi-static simulations. However, the numerical integration scheme requires very small time steps for capturing the stress waves. Although this issue can be overcome using for example mass scaling techniques [START_REF] Ye | A selective mass scaling method for shear wave propagation analyses in nearly incompressible materials[END_REF], this kind of dynamic approaches is not considered in this work. Hence ,the development of efficient solution schemes is essential to deal with the unstable discontinuous crack propagation in quasi-static simulations.

Explicit scheme

An efficient phase field solution scheme has been proposed by [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and improved by [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]. This kind of approaches, referred to as explicit scheme, consists in solving the mechanical equation at a constant damage for each time step. Once the converged mechanical solution obtained, the damage is then updated for the next step. Since the mechanical and damage problems are decoupled, the problem to be solved becomes convex. This method shows an excellent convergence property even when unstable discontinuous cracking occurs. However, as the convergence on damage is not checked, small time steps are usually demanded. The solution procedure is summarized as follows:

• At time step t, d * = d t-1
• Loop for mechanical equilibrium: find u by minimizing R u (u, d * , t)

• If convergence, update u t = u and internal variables

• Loop for the damage problem: find d by minimizing R d (u t , d) with ḋ ≥ 0

• Update the damage field:

d t = d • Pass to next time step
We remark here the irreversibility condition is assumed to be considered by one of the previously mentioned methods in section 2.4.

This remark holds for the following sections.

3.3. Implicit scheme

Alternate minimization scheme

In order to release the time step constraint, the convergence on damage field must be additionally checked. The alternate minimization scheme [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF] is similar to the explicit scheme, but has an additional loop on the damage field before passing to the next step:

• At time step t, d 0 = d t-1
• Loop on j = 0, 1, 2, , . . . Some variants with different convergence criterion for damage field can also be found in the literature (e.g. [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF][START_REF] Gerasimov | On penalization in variational phase-field models of brittle fracture[END_REF]). This scheme is stable in the sense that the mechanical loop always converges within limited iterations, thanks to the decoupling of the two problems. The main drawback is the low convergence rate of the global damage field.

A modified one-loop implicit scheme

Alternatively, we can consider an implicit method by integrating the damage update in the mechanical loop:

• At time step t, u 0 = u t-1 , d 0 = d t-1
• Mechanical Newton loop on i = 0, 1, 2, , . . . Generally speaking, the iterative Newton procedure should be repeated until the convergence of both mechanical and damage problems. In this one-loop scheme, a necessary condition for the global mechanical equilibrium is the convergence of the damage field. Therefore, the residual R u can be used as a global convergence criterion. This algorithm is particularly attractive, since the integration of damage update into the mechanical loop makes the two problems coupled in a stronger manner and should be helpful for the global convergence. In addition, the acceleration techniques [START_REF] Ramière | Iterative residual-based vector methods to accelerate fixed point iterations[END_REF] available for the mechanical convergence can be easily applied to accelerate the global convergence.

1. Compute a displacement increment δu i+1 = -K -1 u R u (u t-1 , d i , t) 2. Update u i+1 = u t-1 + δu i+1
This approach has been implemented for a phase field modeling of nuclear fuel [START_REF] Helfer | Modélisation par champ de phase de la fissuration des matériaux fragiles: Aspects numériques et applications au combustible nucléaire oxyde, 13ème colloque national en calcul des structures[END_REF]. For overcoming the discontinuous propagation steps, a fictive path loading method [START_REF] Michel | A new numerical methodology for simulation of unstable crack growth in time independent brittle materials[END_REF] is employed in that implementation. This method consists in accepting unstable states during the propagation and re-starting the iteration by considering the displacement increment δu 0 equals to zero.

It is shown that this implicit scheme can find a solution of the non-convex problem with an important number of iterations for the discontinuous crack propagation. Furthermore, we observed that many iterations are performed after the convergence of damage field, if a large jump of crack occurs. The reason is multi-fold. First, due to the jump of crack, the displacement solution is far from the current estimate, this challenges the standard Newton procedure. Second, as mention previously, the splitting of the strain energy increases the nonlinearity of the problem. Even if the damage field remains stable, numerous iterations are still needed. In addition, a large crack may lead to ill-conditioned mechanical systems, the integration of the damage update can strongly perturb the mechanical convergence. Numerically, small damage variations can lead to a significant change of mechanical state.

From a physical point of view, the total dissipated energy in a discontinuous cracking step should be larger than the amount due to the crack growth.

Therefore, the crack stability may be reached before the system falls into its equilibrium state. In order to solve the convergence problem, especially for the mechanical equilibrium, a semi-implicit scheme is proposed.

A semi-implicit scheme

The basic idea is to start with a purely implicit solution for the damage prediction, and then switch to an explicit solution once the damage converged. As mentioned earlier, many iterations are required for the mechanical equilibrium after the damage convergence, and the decoupling of these two problems should help to accelerate this procedure. Therefore, appropriate convergence criteria are needed for detecting the instabilities as well as the convergence of each problem (mechanical and damage). We can summarize the solution scheme as follows:

• At time step t, u 0 = u t-1 , d 0 = d t-1
• Mechanical Newton loop on i = 0, 1, 2, , . . .

1. Compute a displacement increment δu i+1 = -K -1 u R u (u t-1 , d i , t) 2. Update u i+1 = u t-1 + δu i+1
If d i+1 -d i > d and R u (u i+1 , d i , t) > u , then (•) i ← (•) i+1 , return to step 1 for next iteration. If d i+1 -d i ≤ d and R u (u i+1 , d i , t) ≤ u , then (•) t = (•) i+1 , pass to next time step. If d i+1 -d i ≤ d and R u (u i+1 , d i , t) > u , then d t = d i+1 , end
of iterations and pass to the explicit solution step.

• Explicit solution at the constant damage d t : find u by minimizing

R u (u, d t , t)
• If convergence, then u t = u, pass to next time step and start from the implicit loop ContrarilyContrary to purely explicit schemes, the semi-implicit solution is obtained at a converged damage state. Hence, some characteristics (e.g. large time steps) of the implicit scheme are conserved in this semiimplicit method. In the meantime, the explicit solution at a constant damage allows to efficiently overcome the difficulties due to the discontinuous propagation. This will be illustrated in numerical experiments with comparison to the modified implicit method.

Remark that the solution schemes presented in this section are independent of the irreversibility implementation. They can be used for any staggered phase field model.

Discussion on the accuracy of different solution schemes

Assuming the spatial discretization error is small enough, a necessary condition to accurate solutions is full-filling the convergence criteria for both mechanical and damage problems: R u = 0 and

R d = 0 for ḋ > 0 and R d > 0 for ḋ = 0.
As shown previously, the explicit schemes only update displacement and damage fields at staggered steps. The convergence criteria on R u (u, d, t) and R d (u, d) are never strictly verified . In general, explicit schemes only approximate the accurate solutions with sufficiently small time steps.

Fully implicit schemes, both alternate minimization and the modified one-loop scheme, should be able to provide accurate results as they can strictly full-fill the necessary condition at each time step.

The semi-implicit scheme can be regarded as a trade-off between them. At each time step, the scheme starts from a purely implicit solution, but activates an explicit solution as long as the damage field is converged. In discontinuous cracking steps, the converged solution always verifies the R u (u, d, t) = 0, but not necessarily the one on R d (u, d). Indeed, the explicitly computed displacement may introduce a new energy state that can affect the previously computed damage field. However, the influence is expected to be small and limited, since the semi-implicit solution is obtained with a con- Proof. Let u ∈ V ⊂ H 1 (Ω) denote the exact solution,

u h ∈ V h ⊂ H 1 (Ω) the discretized one. u Expl h , u Impl h , u Semi-Impl
h are defined to be the solutions computed respectively by the explicit, implicit and semi-implicit schemes for the same spatial and time discretization h. Assuming u = 0, the following relation can be obtained by Hölder's inequality

n t=1 u h (t) -u(t) L 1 (Ω) = n t=1 ( u h (t) u(t) -1)u(t) L 1 (Ω) ≤ ( u h (t) u(t) -1) L 2 (Ω×Ω T ) u(t) L 2 (Ω×Ω T ) ≤ C h u(t) L 2 (Ω×Ω T ) (27)
where we introduce an error discrepancy factor C h which only depends on the discretization factor h, Ω T denotes the time domain. Assuming the spatial discretization error is neglectable, the convergence of u h with respect to time refinement is defined as: lim h→0 C h = 0. Hence, the convergence of explicit and implicit solutions can be defined as

n t=1 u Expl h (t) -u(t) L 1 (Ω) ≤ C 1 h u(t) L 2 (Ω×Ω T ) n t=1 u Impl h (t) -u(t) L 1 (Ω) ≤ C 2 h u(t) L 2 (Ω×Ω T ) (28)
By definition, we have

n t=1 u Impl h (t) -u(t) L 1 (Ω) ≤ n t=1 u Semi-Impl h (t) -u(t) L 1 (Ω) ≤ n t=1 u Expl h (t) -u(t) L 1 (Ω) (29)
Therefore, if lim h→0 C 1 h = 0 and lim h→0 C 2 h = 0, then

lim h→0 n t=1 u Semi-Impl h (t) -u(t) L 1 (Ω) = 0 (30)
Analogically, the above equality holds for u = 0. The proof is closed.

Implementation in Cast3M

The proposed method can be easily implemented in the code Cast3M without supplementary developments. Particularly, the material behavior generator MFront [START_REF] Helfer | Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the pleiades fuel element modelling platform[END_REF] is used here for implementing the softening mechanical response due to the cracking. For summarizing, the overall phase field algorithmic implementation is illustrated in Algorithm 1. In addition, this work makes use of the convergence acceleration tools available in Cast3M, i.e. the fix point acceleration techniques [START_REF] Ramière | Iterative residual-based vector methods to accelerate fixed point iterations[END_REF], when solving the mechanical problem.

Discussion on the choice of convergence criteria

The proposed semi-implicit solution scheme requires several convergence checks. Generally speaking, different criteria can be used. For example, for the mechanical loop, we can check also the convergence on displacement Algorithm 1: Semi-implicit staggered phase field implementation Input: Solution at previous instant: u t-1 , d t-1

Output: Solution at current instant: u t , d t 1 Initiation: δu 0 = 0, δd 0 = 0 2 for i = 0, . . . , i max do

if d i+1 -d i ≤ d & R u (u i+1 , d i ) ≤ u then 8 (•) t ← (•) i+1 9 End Loop 10 else if d i+1 -d i ≤ d & R u (u i+1 , d i ) > u then 11 d t ← d i+1 12
Solve u t at constant damage d t // Algorithm 3

Comparison of different irreversibility implementations

In order to demonstrate the capability of the proposed irreversibility implementation, we compare three different methods: the method introducing H [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] (referred to as H-model), standard Lagrangian method [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] (referred to as variational model), and the proposed method of this work. Particularly, the explicit solution scheme is employed here for its robustness to discontinuous crack propagation. In the following experiments, relatively fine loading steps have to be used, which are considered as follows: δu = 6.71 × 10 -5 mm for the first 80 steps, then δu = 6.71 × 10 -6 mm. The Lagrangian method is implemented in Cast3M using an active-set method. expected, the proposed method can produce a very similar result to the variational model, which confirms the equivalence between them. However, the proposed method is much more efficient than the variational model using the existing implementation for the variational inequality, as shown in Table 1.

We remark here that the validity of these solutions is out of the scope of this discussion. For obtaining a realistic crack pattern with an explicit scheme, a much finer time discretization is needed. However, this is not a problem for the comparison of different numerical behaviors conducted by different phase field models. The use of explicit scheme is for the purpose of giving similar constant input energies to different models at each time step.

At the end of the first comparison, we illustrate additionally the numerical dissipation energy of the proposed model during the crack propagation (see Figure 4). The dissipation is computed as follows Ĝc = δ(Π ext -Π e ) 2δ crack length [START_REF] Farrell | Linear and nonlinear solvers for variational phasefield models of brittle fracture[END_REF] we recall that Π ext and Π e stand for the external work and the elastic energy stored in the cracked body respectively. In a continuous cracking step, the numerical dissipation should be equal to the theoretical one that is prescribed in the phase field formulation. As shown in the figure, the initiation of cracking needs a much higher energy. Then, due to the discontinuous crack propagation, the computed dissipation differs significantly from the theoretical one at the beginning of propagation. As the cracking becomes continuous, the numerical dissipation converges to the theoretical reference. This confirms the energetic aspect of the proposed irreversibility implementation.

Comparison of different solution schemes

In the second part of this test-case, we compare different solution schemes: explicit, implicit (the modified one) and semi-implicit schemes. The proposed irreversibility implementation is adopted, as an alternative implementation for the variational inequality. The internal characteristic length l c = 0.015 mm is used here. The loading increment is considered as follows: δu = 6.1 × 10 -4 mm for the first 8 steps, δu = 6.1 × 10 -5 mm for the remaining 20 steps until u = 6.1 × 10 -3 mm. It should be noticed that these loading steps are usually considered too large for an explicit solution in such experiments.

Suitable loading steps for accurate explicit solutions should be at the most of the order of 10 -6 mm [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF].

As shown in Figure 5, these loading steps are too large for an explicit scheme. The final crack patterns are strongly dependent of the time discretization. For obtaining a full crack, the loading steps have indeed to be 100-times smaller. On the contrary, the implicit scheme does not impose such requirement for the loading steps. The full crack is obtained without refining time discretization. The proposed semi-implicit scheme has similar properties. Figure 6 illustrates the crack evolution in the specimen. It is shown that the crack initiates around u = 5.612 × 10 -3 mm and goes through the entire specimen within one time step. This discontinuous propagation appears in explicit solutions only when the time step decreases sufficiently. We can see that the final crack pattern is completely independent of loading steps with the semi-implicit method in this example. Larger loading steps can be used without perturbing the final crack pattern.

However, if we take a look at Figure 7, the reaction force is still sensible to the time discretization, even with implicit or semiimplicit schemes. In order to have a good representation of loading history, the time steps should not be too large. In general, semiimplicit solutions should converge with time refinement at a better rate than explicit solutions. Taking the maximal reaction force of the implicit solution as a reference, we can illustrate the convergence trends (see Figure 8). It is shown that the semi-implicit solutions converge much faster than the explicit solutions by refining the time steps.

In this example, the discontinuous propagation causes a sudden drop of external forces and a big jump of state on displacement field. Therefore, numerous iterations are needed to find mechanical equilibrium, and sometimes, no convergence can be found within a limited time with the purely implicit scheme.

The semi-implicit scheme is notably more efficient than the implicit one.

As shown in Table 2, the semi-implicit solutions take only several hours in this experiment, whereas the implicit solutions encounter severe difficulties of convergence and need more than 20 hours to find global mechanical equilibrium. This implicit scheme is able to bypass the discontinuous cracking step with numerous iterations, but still hardly converges in the following steps. This happens as well in other experiments, this point will be more discussed in the second example. The semi-implicit scheme is able to overcome these difficulties by fixing the damage variable. For the same reason, explicit solutions are also very robust for the discontinuous crack propagation. However, as shown in Table 2, the computational cost of explicit solutions increases exponentially as the time step decreases and can rapidly become un-affordable if a high accuracy of results is demanded.

Asymmetric double notched tensile specimen

The second experiment consists in the well studied asymmetric double notched specimen [START_REF] Melin | Why do cracks avoid each other?[END_REF][START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF][START_REF] Schwaab | Interacting cracks obey a multiscale attractive to repulsive transition[END_REF]. Figure 10 illustrates the crack evolution obtained with the proposed phase field implementation and semi-implicit scheme. The final crack pattern shows an excellent agreement with the experimental observation in many brittle materials [START_REF] Schwaab | Interacting cracks obey a multiscale attractive to repulsive transition[END_REF]. Similar numerical results can be obtained with explicit schemes using small time steps, as reported in [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]. Physically, this repulsive behavior of two parallel cracks is completely possible within linear elastic fracture mechanics theory, as explained by the work [START_REF] Schwaab | Interacting cracks obey a multiscale attractive to repulsive transition[END_REF], although the attraction or repulsion depends strongly on the geometry condition of the two approaching cracks. One explanation is that the propagation direction is altered by the interaction between the stress fields around the crack tips, as they get closer.

The crack angle θ with respect to initial direction intends to promote the pure opening mode, i.e. stress intensity factor K II (θ) = 0.

Interested readers can refer to [START_REF] Schwaab | Interacting cracks obey a multiscale attractive to repulsive transition[END_REF] for more details.

In the numerical experiment, initial cracks start propagating around u = 4.11 × 10 -2 mm and the en passant cracks appear within one time step. This discontinuous propagation requires many iterations for the convergence of solutions, as shown in Figure 11. Table 3 summarizes the computational cost for different solution cases.

The implicit scheme did not converge after a long time for the discontinuous cracking steps. It is shown that the semi-implicit scheme is much more robust and efficient than the purely implicit scheme with a significant speedup. Figure 13 shows the evolution of crack for a maximum displacement loading u = 6 × 10 -2 mm. As expected, the crack grows vertically and stops near the elastic region. Figure 13 shows the crack evolution for two internal lengths.

Similar results have been reported in the literature (see e.g. [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]).

In our numerical experiments, a smaller internal length leads to stronger numerical instabilities, which causes therefore a higher computational cost (see Table 4). One of the reasons may relate to the mesh size effect, as the mesh becomes closer to the internal length when decreasing the latter. Regardless of the strong numerical instabilities, the semi-implicit scheme is able to converge at a limited time cost. crease. Therefore, we can adopt large time steps without modifying final crack patterns, using the implicit or semi-implicit schemes. However, this is not the case for this heterogeneous material, loading steps have to be carefully chosen for accurately representing the loading history. Numerically, it can be noticed that the irreversibility condition is implemented in a time-discrete manner. which clearly shows the history-dependence nature of the underlying problem. We remark that this dependency is true for any kind of solution schemes. The semi-implicit scheme can accept relatively larger time steps in any cases, compared to explicit ones. This example emphasizes the adequate choice of time steps, even with the implicit or the proposed semi-implicit schemes.

However, we can expect better convergence rates against the explicit schemes. Indeed, explicit schemes have been tested with the finest time steps: δu/10, the solution is far from convergence.

The computational cost of each solution is summarized in Table 5. This

shows again the efficiency of the proposed phase field implementation. conditions for numerical analysis is illustrated in Figure 17. The plate is initially subjected to a uniform temperature T 0 . From time t > 0, a colder temperature T 1 is prescribed on the upper side. All the exposed surfaces are considered adiabatic. Assuming the length of plate is sufficiently long, the temperature field at t > 0 can be analytically given by

T (x, y) = T 0 -(T 1 -T 0 )f c ( y 2 √ k c t ), ∀t > 0 (32)
where k c is the thermal conductivity, f c is the complementary error function:

f c (x) = 2 √ π ∞ x e -s 2 ds.
Therefore, the temperature is uniform in the direction x, while a high temperature gradient appears in the direction y. The elastic strain induced by the thermal expansion reads then

ε e = ε -ε th = ε -α(T -T 0 )I ( 33
)
where α is the thermal expansion coefficient, I is the second order identity matrix.

In this work, the loading parameters and material properties are considered as the same as in [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF]. Some important loading parameters and material properties are given in Table 6. In order to have results comparable with those reported in [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF][START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF], we do not differentiate the compression and traction effect in the phase field model. The FE model is generated using linear triangular elements with a plane stress assumption. In order to well represent the thermal shock history, the time steps are considered as follows:

δt = 5 × 10 -4 for the first two steps, then δt = 1 × 10 -3 for the remaining steps.

As shown in Figure 18, the crack starts being homogeneous in the direction parallel to the surface of the thermal shock. At some critical time, the homogeneous solution bifurcates towards a periodical solution with equal distance cracks penetrating inside the specimen. The wave length, which stands for the distance between two neighboring cracks, is initially equal to 7.4 times of the internal characteristic length. After some time, some cracks stop to propagate whereas the others continue with a wave length approximately two times higher than the original one. These results show an excellent agreement with both analytical analyses [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF] and experimental results [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF]. In terms of computational cost, this experiment takes only several hours with the proposed phase field implementation. However, full Lagrangian methods will take several days for this kind of computations.

Conclusion

A novel phase field method for quasi-static brittle fracture analysis has been developed. This method is based on two novel algorithmic implementations: a novel efficient algorithm for imposing the irreversibility condition and a robust staggered semi-implicit solution scheme for overcoming the discontinuous propagation and time step constraints.

The irreversibility implementation is based on an efficient implementation of the inequality constrained optimization procedure with vanishing energetic driving force. Unlike the H field based model, the proposed method can keep the original variational nature of the phase field solution. Moreover, this method does not introduce any supplementary numerical coefficient which may result in ill-conditioned systems. This method can be considered as a variant to conventional variational phase field models.

The proposed semi-implicit staggered scheme relies on two ingredients: the integration of phase field solution into the mechanical loop, and the combination of the purely implicit and explicit solution schemes. This method allows to alleviate the time step constraints, while being very robust with the numerical instabilities associated with the discontinuous propagation. In the above procedure, the mapping introduces naturally an interpolation error, but we can expect this error is controllable by refining the mesh and choosing appropriate interpolation methods.

When applying the vanishing energy condition to a constrained point, a part of the reference energy ψ + e will be removed from its surrounding elements. Hence, the computed damage can be less important than that of full Lagrangian method.

In order to illustrate the impact of the method, let us consider a This energy intends to create a thinner crack at x = 3 mm with a higher peak value and initiate a new crack at x = 7 mm. Hence, the damage field around the initial crack is decreased if the irreversibility condition is not imposed. We use the proposed method to compute the corresponding damage fields.

Figure 1 :

 1 Figure 1: Approximation of crack topology

1 .

 1 Compute the damage d * i with d * i = d t-1 on the constrained set D = {x ∈ Ω| p(x) = 1} 2. Determine the decreasing points D * = {x ∈ Ω| d * i < d t-1 } 3. Update the constrained set D = D D * 4. Update the indicator function p(x) = 1, ∀x ∈ D 5. Check the convergence: If D

1 . 2 . 3 .≥ 0 4 . 5 .

 12345 Loop for mechanical equilibrium: find u by minimizing R u (u, d j , t) If convergence, update u t = u and internal variables Loop for the damage problem: find d by minimizing R d (u t , d) with ḋ Update the damage field: d j+1 = d If d j+1 -d j ≤ d , then d t = d j+1 , pass to next time step 6. Otherwise, repeat the j-loop

and internal variables 3 .≥ 0 4 . 5 .

 345 Loop for the damage problem: find d by minimizing R d (u i+1 , d) with ḋ Update the damage field: d i+1 = d If convergence, then u t = u i+1 , d t = d i+1 , pass to next time step 6. Otherwise, repeat the Newton loop where K u stands for the tangent stiffness matrix or the initial elastic stiffness matrix for the modified Newton method.

and internal variables 3 .≥ 0 4 . 5 .

 345 Loop for the damage problem: find d by minimizing R d (u i+1 , d) with ḋ Update the damage field: d i+1 = d Convergence check:

Proposition 1 .

 1 verged damage field at each time step. In other words, the criterion on R d (u, d) is expected to be verified if the damage field has well converged. The semi-implicit solution should globally approximate the accurate fully implicit solution. In particular, we have the following result. If the explicit and implicit solutions can both converge to the exact solution with time refinement, then the semi-implicit solution converges to the exact solution with time refinement.

15 (4. Numerical experiments 4 . 1 .

 1541 •) i ← (•) i+116 Return u t , d t variation: δu i+1 in addition to the global equilibrium R u = 0, or just use one of these two criteria instead. In principle, the global equilibrium ensures automatically the convergence on displacement variation. Conversely, it is not true. The influence of different criteria on computational cost can be studied, but is out of the scope of this work. In any cases, the relative tolerance for mechanical and damage convergence check, respectively related to u and d (Algorithm 1 and 3), should be at least of order 10 -4 for a good accuracy. The irreversibility tolerance c used in Algorithm 2 should be at least of order 10 -2 . In the following numerical examples, both displacement stability and equilibrium residual are used for the mechanical convergence. The damage convergence is considered by the crack stability, as shown in the algorithms. Single edge notched tensile testThe first test concerns the well known single edge notched tensile test[START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. Figure2illustrates the geometry and prescribed boundary conditions.The displacement loading is imposed on the top side of the specimen while keeping the bottom side fixed. The material properties are set to the same as[START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]: Young's modulus E = 210 kN/mm 2 , Poisson's ratio ν = 0.3, critical energy release rate G c = 2.7×10 -3 kN/mm. The FE mesh is generated using linear triangular elements and refined around the expected crack path for a size about twice smaller than the concerned regularization length.

Figure 2 :

 2 Figure 2: Geometry and boundary conditions for single edge notched test [mm]

Figure 3

 3 Figure 3 illustrates the crack patterns for different models. Due to the loss of equivalence to the original variational formulation, the H-model leads to a crack profile different from the variational model. Although the crack length at the current loading stage is overall the same in this example, this difference in crack profiles shows that the different models can result in completely different cracking behaviors different numerical behaviors and potentially lead to different local minima and convergence rates. As

Figure 3 :

 3 Figure 3: Crack patterns at u = 6.71 × 10 -3 mm provided by the explicit solution scheme for different models

Figure 4 :

 4 Figure 4: Numerical dissipation during the crack propagation with the proposed irreversibility implementation

  (a) Explicit, δu (b) Explicit, δu/10 (c) Explicit, δu/100 (d) Implicit, δu (e) Semi-implicit, δu (f) Semi-implicit, δu × 2

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Final crack patterns in the single edge notched tensile test with different solution schemes

Figure 8 :

 8 Figure 8: Convergence studies with respect to time step refinement

Figure 9

 9 illustrates the geometry and prescribed boundary conditions. The used materials properties are the same as[START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]:E = 210 kN/mm 2 , ν = 0.3, G c = 2.7 × 10 -3 kN/mm, l c = 0.2 mm. TheFE mesh is generated using linear triangular elements and refined around the expected crack paths for a mesh size up to 0.1 mm, which is twice inferior to the internal characteristic length l c . This model contains 26297 nodes and 52520 elements. The displacement increment is considered as follows: δu = 5.01 × 10 -3 mm for the first 8 steps, δu = 5.01 × 10 -4 mm for the remaining 20 steps until u = 5.1 × 10 -2 mm.

Figure 11 (

 11 a) depicts the iteration numbers of each time step with the one-loop implicit and the semi-implicit schemes. The most difficult time step for the implicit scheme turns out to be the one (step 11) right after the discontinuous cracking. This may be explained by the fact that the system becomes highly nonlinear with respect to the displacement due to the splitting of strain energy and ill-conditioned with a large crack. The strong coupling of displacement and damage fields amplifies the perturbation of damage variation on the mechanical convergence. Experience shows that even small variations in the damaged region can severely perturb the iterative Newton procedure. As shown in Figure11(b), the mechanical problem converges very slowly with the purely implicit scheme, despite of the globally stable damage field. However, by decoupling the mechanical and damage problems and considering that the damage solution has been found, the solution can quickly converge to a local minimum. This is shown by the excellent convergence behavior of the semi-implicit solution.

ure 12 )Figure 9 :Figure 10 :

 12910 Figure 9: Geometry and boundary conditions for double notched tensile specimen [mm]

Figure 11 :

 11 Figure 11: Iteration number and iterative residual

Figure 12 :

 12 Figure 12: Geometry and boundary conditions for three points bending test [mm]

Figure 13 :

 13 Figure 13: Crack evolution in the three points bending test

Figure 15

 15 Figure 15 illustrates the final crack patterns obtained with different loading steps. The evolution of the reaction forces is shown in Figure 16. In this example, a loading step dependency is observed. This is associated with the intrinsic history-dependence nature of the underlying fracture problem. In previous examples, where the materials are considered homogeneous, stresses induced by the applied loading are uniform and monotonically in-

Figure 14 :

 14 Figure 14: Geometry and boundary conditions for the two-phase concrete material [mm]

10 Figure 15 :

 1015 Figure 15: Final crack patterns of the concrete material obtained with different step sizes.

Figure 16 :

 16 Figure 16: Evolution of the reaction force in the concrete material

Figure 17 :

 17 Figure 17: Geometry and boundary conditions for thermal chock test [mm]

Figure 18 :)

 18 Figure 18: Crack evolution in the thermal shock test. (a) (b) (c): Simulated cracks with l c = 5 × 10 -2 mm, (d): Experimental results [45].

2 (A. 3 ) 2 (

 232 1D bar with G c = 2.7 × 10 -3 kN/mm and l c = 0.015 mm. We assume a piece-wise constant element-based reference energy ψ + e is given in different time steps. Initially, a predefined crack is computed with the following energy distributionψ + e (x, t = 0) = 0.5e -|x-3| 0.At the first time step t = 1, the strain energy is changing to ψ + e (x, t = 1) = 5e -|x-3| 0.05 + e -|x-7| 0.

  Figure A.19 depicts the crack profiles of different time steps. We can see that the damage field did not decrease with the proposed irreversibility implementation. Compared to the full Lagrangian method, the damage field is indeed degraded with the proposed method, as shown in

Figure A. 19

 19 Figure A.19(b). Fortunately, the difference is limited by refining the mesh (Figure A.19(c)). The final mesh size is reasonable, as it is only slightly inferior to the internal crack length.

Table 1 :

 1 Computational cost for different models in the tensile test

	Model	Solver l c (mm) CPU Time
	Variational model Expl.	0.0075 0.00375	225 min 703 min
	This work	Expl.	0.0075 0.00375	55 min 160 min

Table 2 :

 2 Computational cost for the single edge notched tensile test

	Model	Solver	Step size Full crack Iterations CPU Time
			δu	No	229	2 min
		Expl.	δu/2 δu/10	No No	756 6462	4 min 43 min
	This work		δu/100	Yes	97691	506 min
		Impl.	δu	Yes	-	>1000 min
		Semi-impl.	δu δu × 2	Yes Yes	17467 14329	348 min 306 min

Table 3 :

 3 Computational cost for the double notched tensile test

	Model	Solver	Iterations CPU Time
	This work	Impl. Semi-impl.	>17685 >1000 min 7673 269 min
	4.3. Symmetric three points bending test
	Next, we investigate the performance of proposed methods in a different
	loading case. The symmetric three points bending test is used (see Fig-

Table 4 :

 4 Computational cost for the bending test

	Model	Solver	L c (mm) Iterations CPU Time
	This work Semi-impl.	0.06 0.03	6653 70959	137 min 780 min
	4.4. Crack nucleation and propagation in a two-phase concrete material
	Then, we consider a heterogeneous specimen. Figure 14 illustrates a two-

phase simplified concrete material without any initial crack. Particularly,

Table 5 :

 5 Computational cost for the two-phase concrete specimen

	Model	Solver	Step size Iterations CPU Time
			δu	1002	132 min
	This work Semi-impl.	δu/5	4344	391 min
			δu/10	8649	387 min
	4.5. Thermal shock test
	The final test concerns a thermal shock problem [44, 45], which is well
	studied both experimentally and analytically. The geometry and boundary

Table 6 :

 6 Loading parameters and material properties for the thermal chock test

T 0 (K) T 1 (K) E (GPa) ν α (K -1 ) G c (kN/

mm) 673 293 370 0.22 8.4 × 10 -6 12.16×10 -6

Acknowledgement

This research was conducted in the framework of the 'PLEIADES' project, which is supported financially by the CEA (Commissariat à l'Energie Atomique et aux Energies Alternatives), EDF (Electricité de France) and Framatome. In addition, the authors would like to acknowledge the financial support of the Cross-Disciplinary Program on Numerical Simulation of CEA (PTC FUEL).

Detect new constraint set:

End Loop

11 Return u t Appendix A. Numerical study of the impact of the irreversibility implementation

The important features of the proposed implementation are:

irreversible updating of the constrained set and vanishing energetic force. The impact of the first feature is easier to see. Due to the updating strategy, the resulting constraint set may not be optimal for a given energy state. Consequently, the computed damage field may be somehow degraded.

The impact of vanishing energy is less obvious. Unlike the analytic formulation, the vanishing energy ψ + e = 0 can not be directly imposed with discretized formulation for the constrained nodes, since the strain energy is usually computed inside an element for integration points. To do this, a mapping of the energy between its