A micromorphic phase-field model for brittle and quasi-brittle fracture

The phase-field model for fracture, despite its popularity and ease of implementation comes with its set of computational challenges. They are the non-convex energy functional, variational inequality due to fracture irreversibility, the need for extremely fine meshes to resolve the fracture. In this manuscript, the focus is on the numerical treatment of variational inequality. In this context, the popular history-variable approach suffers from variationally inconsistency and non-quantifiable nature of the error introduced. A better alternative, the penalisation approach, has the potential to render the stiffness matrix ill-conditioned. In order to circumvent both aforementioned issues, a micromorphic approach towards phase-field fracture modelling is proposed in this manuscript. Within this approach, a micromorphic extension of the energy functional is carried out. This transforms the phase-field into a local variable, while introducing a micromorphic variable that regularises the fracture problem. This reduction the regularity requirements for the phase-field enables an easier implementation of the fracture irreversibility constraint through simple 'max' operation, with system level precision. Numerical experiments carried out on benchmark brittle and quasi-brittle problems demonstrate the applicability and efficacy of the proposed model for a wide range of fracture problems.

The efficacy of the novel micromorphic phase-field fracture model is demonstrated with benchmark problems in both brittle and quasi-brittle fracture. These include a single edge notched specimen under tension, the Winkler L-panel experiment [47], and a three-point bending experiment carried out in [48]. This manuscript is structured as follows: Section 2 introduces the reader to the phase-field model for fracture, its underlying energy functional and pertinent Euler-Lagrange equations. Subsequently, in Section 3, the micromorphic approach for phase-field fracture is introduced. The numerical benchmark problems are addressed in Section 4, followed by concluding remarks in Section 5.
2 Phase-field fracture model 2.1 The energy functional Figure 1a and 1b illustrates a discrete fracture and a phase-field regularised fracture in a 2D continuum respectively. The salient feature of the a discrete fracture model is the explicit representation of the fracture through displacement discontinuities. The phase-field fracture model belongs to the smeared representation of fracture, where the phase-field, ϕ ∈ [0, 1] interpolates between intact and fully broken (fractured) material states. The fractured continuum in both figures is assumed to occupy a domain Ω ∈ R dim (dim = 2 in this case). The external boundary Γ comprises of the Dirichlet and Neumann boundaries, represented by Γ u D and Γ u N respectively. Note that Γ = Γ u D ∪ Γ u N and Γ u D ∩ Γ u N = ∅.  . diffused (smeared) crack, with Dirichlet and Neumann boundaries indicated as Γ u D and Γ u N respectively. Figure adopted from [32] .
The energy functional pertaining to the phase-field fracture model [33] is given by, where, a degradation function g(ϕ) is attached to the fracture driving strain energy density Ψ + , with Ψ − and t u p being the residual strain energy density and the prescribed traction respectively. The last integral in the above equation represents the fracture energy, where G c and l are the Griffith fracture energy and the fracture length-scale respectively, while c w is a normalisation constant associated with the choice of the locally dissipated fracture energy function w(ϕ). Some commonly adopted phase-field fracture models are presented in the Table 1. The quasi-brittle phase-field fracture model requires additional parameters p, a 1 , a 2 and a 3 , as observed from the Table 1. These parameters are chosen in way, such that different traction-separation laws may be obtained. For a detailed derivation of these parameters, the reader is referred to [22]. Following [22], the constant a 1 is given by, Quasi-brittle [22] 2ϕ − ϕ 2 π (1 − ϕ) p (1 − ϕ) p + a 1 ϕ + a 1 a 2 ϕ 2 + a 1 a 2 a 3 ϕ 3 Table 1: Brittle and quasi-brittle phase-field fracture models.
where, the newly introduced parameters E 0 and f t represent the Young's Modulus and the material tensile strength respectively. The other parameters vary for different traction-separation laws, and given in Table 2.  Table 2: Quasi-brittle phase-field fracture model parameters for different traction-separation laws [22].
As demonstrated in Equation (1), the strain energy density is split into fracture driving and residual components, Ψ + and Ψ − . In this manuscript, the common consensus that a fracture occurs under tensile loading is adopted. Therefore, the fracture driving and residual strain energy densities are computed using the Amor split [49]. With this approach, the fracture is driven by the deviatoric part of the strain energy density together with a positive volumetric part of the strain energy density. The negative volumetric part corresponds to the residual energy. Mathematically, and Here, K and µ represents the bulk and the shear modulus of the material respectively, tr is a trace operator, • ± represents the positive/negative Macaulay brackets, and dev is the deviatoric strain. The corresponding Cauchy stresses are given by, and adopting the heaviside function H.
In order to simulate fracture initiation and propagation in continuum Ω, minimisation of the energy functional (1) is required, w.r.t. its solution variables u and ϕ. This results in the Euler-Lagrange equations. Along with appropriately defined test and trial spaces, they result in the following problem: using pertinent time-dependent Dirichlet boundary conditions u p on Γ u D and ϕ p on Γ ϕ D , and and Neumann boundary condition t u p on Γ u N . The trial and test spaces are defined as Note that the requirement ϕ ≥ n ϕ in (8c) ensures fracture irreversibility, with n referring to the previous time-step.
From Problem 1, it is clear that the fracture irreversibility constraint manifests in the form of a variationally inequality Euler-Lagrange equation (7b) with restrictive trial and test set (8c) for the phase-field. As mentioned in Section 1, several researchers have proposed different methods to treat the variational inequality problem. In the next section, a micromorphic phase-field fracture model is proposed, where the phase-field variable becomes a local quantity, while a micromorphic variable ensures regularisation of the problem. This enables a simpler 'point-wise' treatment of the fracture irreversibility constraint.

Energy functional and variants
In this sub-section, a micromorphic extension of the phase-field fracture energy functional (1) is carried out in the spirit of [40]. This results in, Here, d is the 'new ' micromorphic variable that regularises the fracture problem. From the above equation, it is observed that the micromorphic variable d is solely associated with the gradient term in the fracture energy integral, and the phase-field ϕ becomes a local quantity. Consequently, the regularity requirements on the phase-field w.r.t. the existence of derivatives is circumvented. Moreover, the micromorphic approach also introduces an additional energy term associated with the difference between the phase-field and the micromorphic variable. Theoretically, for the continuous problem, in the limit, the interaction parameter, α → ∞, the original energy functional (1) is recovered.

Euler-Lagrange equations
The set of Euler-Lagrange equations for the micromorphic phase-field fracture model is obtained upon minimising the energy functional (9) w.r.t. its solution variables u, ϕ and d. Along with appropriately defined test and trial spaces, it results in the following problem: using pertinent time-dependent Dirichlet and Neumann boundary conditions, u p on Γ u D , and t u p on Γ u N respectively. The trial and test spaces are given by, In Problem 2, the phase-field evolution equation (10b) is local. As such, ϕ can be computed 'point-wise' in the computational domain. In particular, the root(s) of the possibly nonlinear scalar equation, yield(s) the phase-field ϕ, if n ϕ < ϕ < 1. Using the locally computed phase-field, the complete problem for micromorphic phase-field fracture is assumes the form: using pertinent time-dependent Dirichlet boundary conditions u p on Γ u D , and Neumann boundary condition t u p on Γ u N . The trial and test spaces are defined as with the local phase-field ϕ computed using (12).
It is worth mentioning that the local phase-field evolution (12) is linear for brittle AT1 and AT2 fracture models. This is due to the quadratic nature of the degradation function g(ϕ) = (1 − ϕ) 2 coupled with linear/quadratic locally dissipated fracture energy function w(ϕ) = ϕ and ϕ 2 for the AT1 and AT2 models respectively. This yield explicit expressions for the local phase-field variable, for AT1, and ϕ = min max 2Ψ + + αd for AT2 model respectively. In the case of quasi-brittle phase-field fracture model, the rational nature of the degradation function (see Table 1) results in nonlinear scalar equation for the local phase-field variable, warranting the need for the Newton-Raphson method.

Discrete equations
The work in this manuscript is executed within the framework of the finite element method [51,52], using triangular (T3) elements for both, displacement and micromorphic fields. This allows assuming the displacement and the micromorphic fields at the nodes (ũ i ,d i ) as the primary unknowns, with the corresponding continuous fields (u,d) approximated as, In the above equation, N u i and N d i are the interpolation functions for the displacement and the micromorphic phase-field, associated with the i th node. The spatial derivatives of the interpolation functions N u i and N d i in a two-dimensional case are given by, Here, the subscripts , x and , y indicate spatial derivatives in x and y directions respectively. Using (18), the strain , and the gradient of the micromorphic variable ∇d are defined as, The discrete phase-field fracture problem is obtained upon inserting (17)(18)(19) in the Euler-Lagrange equations from Problem 3. Thereafter, (13a) and (13b) are assumed as the internal forces, and stiffness matrix derived from its derivative. This notation is consistent with [53], and allows the presentation of the phase-field fracture problem in the incremental iterative framework as: Discrete Problem 1. Compute the solution increment (∆ũ, ∆d) i+1 in the current iteration i + 1 using and update the solution fields, until a certain convergence measure is fulfilled. The local element stiffness matrices are computed as: and the local internal force vectors are computed as The external force vectors f int,u and f int,d are considered equal to zero. Computation of the material stiffness D in Voigt notation is presented in Appendix A.
The discrete problem 1 is non-convex. As such, conventional incremental solution techniques, like the Newton-Raphson method fails to achieve convergence in the softening regime, possibly due to an indefinite stiffness matrix. Efforts to circumvent this issue include the novel line search technique proposed in [34], the use of arc-length solvers [54][55][56], trust-region methods [57] and convexification via extrapolation of the phase-field for the momentum balance equation [36]. These aforementioned techniques are within the framework of monolithic solution techniques. In [39], the alternate minimisation solution techniques was proposed, since the phase-field energy functional is convex with respect to either of the coupled field, if the other one is held constant. In this manuscript, a monolithic solution technique is desired. As such, the convexification strategy proposed by [36] is adopted. With this choice, the phase-field for the momentum balance equationφ, is computed using, withd being the linearly extrapolated micromorphic variable from the last two converged (time) steps. It is important to note that this extrapolation-based computation of the phase-field applies only to the momemtum balance equation and not the micromorphic variable evolution equation. For the micromorphic variable evolution equation, the phase-field is computed using (12), restated below for clarity, Consequently, the discrete problem assumes the form: Discrete Problem 2. Compute the solution increment (∆ũ, ∆d) i+1 in the current iteration i + 1 using and update the solution fields, until a certain convergence measure is fulfilled. The local element stiffness matrices are computed as: and the local internal force vectors are computed as The local phase-fieldsφ and ϕ are computed using (21) and (22) respectively. The external force vectors f int,u and f int,d are considered equal to zero. Computation of the material stiffness D in Voigt notation is presented in Appendix A.

Numerical Study
In this section, numerical experiments are carried out on brittle and quasi-brittle phase-field fracture problems. The problems in brittle fracture include the single edge notched specimen under tension and shear, and the three-point bending test [39], whereas for quasi-brittle fracture, the Winkler concrete L-panel experiment [47] and the concrete three point bending experiment [48] are analysed. For each problem, the geometry, loading conditions as well as the additional model parameters are presented in the respective sub-sections. The phase-field fracture topology in the final step of the analysis, and the load-displacement curves are also shown therein. All problems are solved in a fully coupled (monolithic) sense, adopting the Newton-Raphson method. The iterative procedure is terminated when an error measure defined as ratio of the norm of the residual in the current iteration to that of the first iteration is less than a certain tolerance, tol. In this manuscript, tol = 1e − 3. The linear problem within each iteration is solved using the GMRES solver, preconditioned with an incomplete LU preconditioner. Moreover, for all numerical experiments, the interaction parameter α is parametrised as with β being a user-defined non-dimensional scalar.

Single Edge Notched specimen under Tension (SENT)
The single edge notched specimen [39] has been studied extensively under tensile and shear loading in the phase-field fracture literature. The geometry consists of a unit square (in mm) embedded with a horizontal notch, midway along height and equal to half of the edge length as shown in Figure 2. The notch is modelled explicitly in the finite element mesh. A quasi-static loading is applied at the top boundary in the form of prescribed displacement increment ∆u = 1e − 4[mm] for the first 55 steps, following which it is changed to 1e − 6[mm]. The bottom boundary remains fixed. The additional model parameters are presented in Table 3.

Single Edge Notched specimen under Shear (SENS)
In order to perform a shear test, the single edge notched specimen is loaded horizontally along the top edge as shown in Figure 3. The material properties remain same as presented in Table 3. A quasi-static loading is applied to the top boundary in the form of prescribed displacement increment ∆u = 1e − 4[mm] for the first 85 steps, following which it is changed to 5e − 6[mm]. Furthermore, the bottom boundary remains fixed, and roller support is implemented in left and right edges thereby restricting the vertical displacement. Figure 6a shows the load displacement curves obtained using β = 10, 100, 250. The curves obtained using β = 100, 250 is similar to that obtained in [59]. Moreoever, the phase-field topology at the final step of the analysis (see Figure 6b) is similar to that presented in Figure 5 of [59]. For β = 10, the non-local behaviour of the micromorphic model is not sufficient. This exaplains the loss of integrity of specimen at an earlier stage in the loading process (cf. Figures 6a and 4a).

Winkler L-panel
The concrete L-shaped panel studied by [47,60] is considered in this sub-section. Figure 7 shows the geometry as well as the loading conditions. The longer edges of the panel are 500 [mm] and the smaller edges are 250 [mm]. The loading is applied on the edge marked in blue, 30 [mm] in length, and is in the form of displacement increments of ∆u = 1e − 3 [mm]. The additional model parameters are presented in Table 4. Miehe [39] Ambati [58] 100 Kristensen [20] 10 250        Figure 9. Displacement-based load increments of ∆u = 1e − 3 [mm] is enforced throughout the simulation. The model parameters are presented in Table 5. Figure 10a presents the load-displacement curves using the micromorphic phase-field fracture model, with β = 100, 250. The range of the experimentally obtained load-displacement curves is exhibited by the shaded region. For both values of β, the load-displacement curves exhibit a good agreement with the experimental region. Moreover, the phase-field topology in the final step of the simulation in Figure 10b is similar to that from the literature [48].

Conclusions and Outlook
A novel phase-field fracture model is proposed in this manuscript, based on the micromorhic extension of the phase-field fracture energy functional. In this model, the phase-field variable is local, and a micromorphic variable is introduced for regularization. Combined with an irreversibility criterion on the phase-field variable, this results in a local variational inequality problem. The problem is then trivially solved pointwise (at integration points) throughout the computational domain. For brittle AT1 and AT2 models, an explicit expression for the local phase-field exists, however, for quasi-brittle fracture models, a nonlinear scalar equation needs to be solved. Furthermore, the local nature of the phase-field variable also provides the ease in implementing bounds, i.e., ϕ ∈ [0, 1] using the trivial 'min' and 'max ' operations. Unlike the penalization technique proposed in [35], fracture irreversibility is achieved with system precision without the need for a high penalty parameter. As regards to computational effort, since the phase-field variable is solved locally, the size of the finite element problem is identical to that of a standard phase-field implementation.  that ϕ ≈ d. This has been demonstrated in this manuscript with numerical experiments on the single edge notched specimen loaded in tension. When the interaction parameter is sufficiently high (≈ 250G c /l), the phase-field and the micromorphic variable profiles are similar for arbitrary sections in a computational domain. This essentially recovers the original phase-field fracture energy functional, while enabling the ease of local phase-field fracture irreversibility. Numerical experiments were also carried out in the context of quasi-brittle fracture, the Winkler L-shaped panel and the concrete three-point bending tests. It is observed that the micromorphic phase-field fracture model results in fracture topologies and load-displacement curves within the experimentally observed range. This demonstrates capability of the micromorphic approach to demonstrate both brittle and quasi-brittle phenomenon. Finally, the novel micromorphic phase-field fracture model opens a plethora of future research extensions, particularly, in multi-physics applications, composite laminates [61]. Other studies may include the implementation of a dissipation-based arc-length method [55,56] or quasi-Newton methods [62,63] for addressing the non-convexity of energy functional.

Software Implementation and Data Availability
The numerical study in Section 5 is carried out using an in-house finite element code, inspired by ofeFrac [64], and based on the Jem and Jive C++ libraries from the Dynaflow Research Group, the Netherlands. The simulation data would be made available in Github repository of the corresponding author (https: //github.com/ritukeshbharali).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.