A robust monolithic solver for phase-field fracture integrated with fracture energy based arc-length method and under-relaxation

1]Ritukesh Bharali
2]Somdatta Goswami

3]Cosmin Anitescu
3]Timon Rabczuk
[1]Department of Industrial and Material Science, Chalmers University of Technology
[2]Department of Applied Mathematics, Brown University, U.S.A

[3]Institute of Structural Mechanics, Bauhaus-Universität Weimar

§ HIGHLIGHTS

* Robust monolithic solver with adaptive under-relaxation and arc-length method.

* Snap-back behaviour is captured with a phase-field fracture energy-based dissipation constraint.

* PHT-splines within IGA framework is utilized for adaptive mesh refinement.

The phase-field fracture free-energy functional is non-convex  with respect to the displacement and the phase field. This results in a poor performance of the conventional monolithic solvers like the Newton-Raphson method. In order to circumvent this issue, researchers opt for the alternate minimization (staggered) solvers. Staggered solvers are robust for the phase-field based fracture simulations as the displacement and the phase-field sub-problems are convex in nature. Nevertheless, the staggered solver requires very large number of iterations (of the order of thousands) to converge. In this work, a robust monolithic solver is presented for the phase-field fracture problem. The solver adopts a fracture energy-based arc-length method and an adaptive under-relaxation scheme. The arc-length method enables the simulation to overcome critical points (snap-back, snap-through instabilities) during the loading of a specimen. The use of an under-relaxation scheme stabilizes the solver by preventing the divergence due to an ill-behaving stiffness matrix. The efficiency of the proposed solver is further amplified with an adaptive mesh refinement scheme based on PHT-splines within the framework of isogeometric analysis. The numerical examples presented in the manuscript demonstrates the efficacy of the solver. All the codes and data-sets accompanying this work will be made available on GitHub (<https://github.com/rbharali/IGAFrac>).

§ INTRODUCTION

The seminal work of Francfort and Marigo [1] led to the emergence of the phase-field based fracture model, as an alternative fracture modelling technique. Therein, the Griffith fracture criterion was cast into a variational setting with certain limitations: no concept of internal length scale and maximum allowable stresses. Later, [2, 3] proposed a computationally convenient framework of the Francfort and Marigo model, adopting a scalar auxiliary variable that interpolates between fully fractured and intact material states. In this context, the Ambrosio-Tortorelli regularization of the Mumford-Shah functional [4] was utilized. Based on the minimization of the global energy function, the phase-field model eliminates the need for tedious tracking of the fracture path and remeshing techniques, frequently observed in the discrete fracture models like XFEM[5]. Furthermore, the phase-field model for fracture has proven its capabilities to handle topologically complex fracture patterns (branching, kinking and merging of cracks)[2].

Soon after the inception of phase field based fracture model, the concept was cast into a thermodynamically consistent framework in [6], adopting an energy-based fracture driving criterion. This work was later extended towards a generic stress-based fracture driving criterion in [7], ductile fracture with plasticity models in [8, 9], anisotropic fracture [10, 11], hydraulic fracture [12, 13], desiccation cracking [14, 15] in a non-exhaustive list of single-scale brittle fracture applications. In the context of multi-scale modelling, the overlapping domain decomposition techniques were adopted in [16, 17, 18, 19, 20], while [21, 22, 23] adopted the hierarchical modelling technique in the FE$^2$ sense [24].

Despite its popularity in several multi-physics domains, the phase-field model has its own set of computational challenges in its implementation in fracture analysis. They include,
(1) a non-convex free-energy functional with respect to the coupled field variables, (2) variational inequality due to fracture irreversibility constraint, and (3) the need to resolve the smeared fracture zone with an extremely fine mesh. The coupled fields can be solved using either a monolithic solver or a staggered solver. Provided that the non-linear solver converges, the monolithic solution scheme is more efficient and faster than the staggered one. However, the non-convex energy functional generally leads to poor convergence and loss of robustness of the monolithic solver. In order to circumvent this, [25] proposed a line search technique which included a negative search direction, [26] proposed convexification of the energy functional based on linear extrapolation of the phase-field for the momentum balance equation. Other methods developed in this context include the arc-length solvers [27, 28, 29], modified Newton-Raphson method [30], error-oriented Newton-Raphson method [31], and trust regions methods [32]. Nevertheless, the development of a robust monolithic solver still remains an active research area in the phase-field fracture community. As an alternative approach, [33] suggested the use of staggered solution scheme, since the energy functional is convex with respect to either of the coupled field, if the other one is held constant.

The second computational challenge associated with the monolithic solver pertains to the variational inequality formulation that arises from the fracture irreversibility constraint. In this context, [25, 34] opted for a simple penalization technique, [26] adopted the primal-dual active set method, whereas [31, 30] used an Augmented Lagrangian method based on the Moreau-Yoshida indicator function [31, 30]. More recent approaches include the micromorphic approach that transforms the phase-field into a local variable [35, 36], and the slack variable approach [37]. Alternatively, a heuristic approach was proposed by [33], replacing the fracture driving energy with its maximum value over the loading history. However, this method is not variationally consistent with the phase-field fracture energy functional [34, 38].

Finally, the phase-field model for fracture analysis requires an extremely fine mesh to resolve the smeared region. A simple and straight forward but computationally expensive way would be to use uniformly refined mesh. If the fracture path is known $\it{a\;priori}$, a certain part of the computational domain may also be pre-refined. The latter case is more applicable when it comes to benchmark models from the literature. However, in the scenarios where the fracture path is not known in advance, adaptive mesh refinement techniques is the preferred option. In this context, the elements are marked based on either a critical threshold value defined over the phase field parameter [26, 39], or by local increase in the tensile energy [40]. Other adaptive refinement schemes include the recovery-based error indicator [41], $\it{a\;posteriori}$ error estimation based on the dual-weighted residual method [42], the finite cell method [43], and the dual mesh concept for the two coupled fields with different mesh refinement indicators [44].

To overcome the issues discussed in the current state-of-the-art, we propose a novel monolithic solver, which is based on an adaptive under-relaxation scheme, and is integrated with fracture energy-based arc-length method. Although the under-relaxation strategies result in decreased rate of convergence of the solver, for the phase-field model, it ensures a guarantee to circumvent of divergence issues arising due to erratic behaviour of the jacobian [19]. The fracture energy-based arc-length method adopted in this work is displacement controlled, which provides the flexibility to take larger displacement steps while accurately capturing the brutal nature of the crack growth. In this work, we have studied displacement controlled fracture, where the load $\it{vs.}$ displacement curve in the post-peak behaviour encounters a discontinuity and the representative point drops onto the lower branch with negative slope, indicating that both load and displacement must decrease to obtain a controlled crack extension. Such observations are often neglected in staggered solvers, but these phenomena are captured accurately using our fast and efficient monolithic solver. Besides capturing the possible snap-back behaviour, the arc-length method also results in an adaptive time-stepping procedure, hence larger energy dissipation is permitted. The adaptive scheme enables a dynamically changing mesh which in turn allows the refinement to remain local at singularities and high gradients. The adaptive $\it{h}$-refinement technique is implemented using polynomial splines over hierarchical T-meshes (PHT-splines). The PHT-splines possess a very efficient local refinement algorithm and they also inherit the properties of adaptivity and locality of T-splines. Moreover, in all the examples, the crack is not initialized but it is allowed to nucleate naturally. The penalization approach is adopted to treat the variational inequality formulation.

This remainder of the article is structured as follows: Section <ref> introduces the reader to the phase-field model for fracture analysis, its underlying energy functional and the pertinent Euler-Lagrange equations. Subsequently, in Section <ref>, the isogeometric analysis framework and the discrete equations for the phase-field fracture model are introduced. Section <ref> presents the main contribution of this work, the robust monolithic solver. The numerical benchmark problems are addressed in Section <ref>, followed by concluding remarks in Section <ref>.

§ PHASE-FIELD MODEL FOR FRACTURE

§.§ The energy functional

Let $\dom \in \mathbb{R}^{\text{dim}}$ (dim $= 1,2,3$) be a domain occupied by a fracturing continuum. A discrete representation of fracture is shown in Figure <ref> where the fracture $\mathcal{C}$ may be represented by a cohesive zone fracture interface. Its diffused counterpart, obtained through the phase-field regularisation is presented in Figure <ref>. Here, the fracture is represented by an auxiliary variable, the phase-field parameter, $\pf \in [0,1]$ within a diffusive (smeared) zone of width $\l>0$, where $\l$ denotes a length scale parameter that controls the width of the diffused zone. The bounds over $\pf$, zero and one indicate the intact material state and total loss of integrity, respectively. Furthermore, the surface $\surf$ of both, the discrete and the diffused fracture continuum is decomposed into a Dirichlet boundary, $\surfarg{D}{u}$ and a Neumann boundary, $\surfarg{N}{u}$, such that $\surf = \surfarg{D}{u} \cup \surfarg{N}{u}$ and $\surfarg{D}{u} \cap \surfarg{N}{u} = \emptyset$.

(K) at (0,0);

[fill=black!10,line width=1pt] (K) plot [smooth cycle,tension=0.7]
coordinates (3,1) (7,1) (8,3) (7,4.475) (5,4.5) (2,4) (1.7,2.5);
[ ] at (3.25,2.15) $\dom$;

[line width=1.5pt,black] (5,0.75) to (5,2.5);
[ ] at (5.35,2.0) $\surfcrack$;

(K) [line width=2.5pt,black] plot [smooth, tension=0.8] coordinates (3,4.3) (1.7,3.7) (1.7,2.5);

[ ] at (1.15,4) $\surfarg{D}{u}$;

[ ] at (6.75,5.05) $\surfarg{N}{u}$;

discrete crack

(K) at (0,0);

[fill=black!10,line width=1pt] (K) plot [smooth cycle,tension=0.7]
coordinates (3,1) (7,1) (8,3) (7,4.475) (5,4.5) (2,4) (1.7,2.5);
[ ] at (3.25,2.15) $\dom$;

[black, path fading=fade out, draw=none] (5,2.5) circle (0.3);
[line width=0.1pt,black] (5,0.75) to (5,2.5);
[line width=0.1pt,black] (5,0.75) to (5,2.5);
[top color=black,bottom color=black!10,shading angle=90] (5,0.78) rectangle (5.3,2.5);
[top color=black!10,bottom color=black,shading angle=90] (4.7,0.78) rectangle (5,2.5);
[->,line width=1pt,black] (4.3,1.7) to (4.85,1.7);
[<-,line width=1pt,black] (5.15,1.7) to (5.7,1.7);
[ ] at (5.85,1.7) $\l$;

(K) [line width=2.5pt,black] plot [smooth, tension=0.8] coordinates (3,4.3) (1.7,3.7) (1.7,2.5);
[ ] at (1.15,4) $\surfarg{D}{u}$;

[ ] at (6.75,5.05) $\surfarg{N}{u}$;

diffused (smeared) crack

A solid, $\dom \in \mathbb{R}^2$, embedded with (a) discrete crack $\surfcrack$ and (b) diffused (smeared) crack, with Dirichlet and Neumann boundaries indicated as $\surfarg{D}{u}$ and $\surfarg{N}{u}$ respectively. Figure reproduced from [23].

A general form of the energy functional for the phase-field fracture model, shown in Figure <ref>, is given by,

\begin{equation}{\label{eqn:sec2:EFunc}}
\displaystyle E(\bdisp,\pf) = \int_{\dom}^{} g(\pf) \Psi^f(\strain(\bdisp)) \: \text{d}\dom + \int_{\dom}^{} \Psi^r(\strain(\bdisp)) \: \text{d}\dom + \int_{\dom}^{} \dfrac{\gc}{c_w \l} \left( w(\pf) + \l^2 |\bnabla \pf|^2 \right) \: \text{d}\dom,
\end{equation}

in the absence of any external loading (body and traction forces). Here, $g(\pf)$ is a monotonically decreasing stress-degradation function attached to the fracture driving strain energy $\Psi^f(\strain(\bdisp))$, and $\Psi^r(\strain(\bdisp))$ is the residual strain energy. Moreover, $\gc$ is the Griffith fracture energy, which is a material parameter, and $c_w$ is a normalisation constant associated with the choice of the locally dissipated fracture energy function, $w(\pf)$. Finally, $\strain(\bdisp)$ is the symmetric part of the deformation gradient, where $\bdisp$ is the displacement field. The choice of $w(\pf) = \pf^2$ ($c_w=1/2$) typically denotes a AT2 phase field model, while the choice $w(\pf) = \pf$ ($c_w=8/3$) is often referred to as the AT1 model.

The phase-field model for fracture allows great flexibility in the choice of degradation function $g(\pf)$ and locally dissipated energy function $w(\pf)$, albeit with some restrictions. The degradation function must satify the following criterion: $g(0) = 1$, $g(1) = 0$, and $g'(1) = 0$, to ensure that the fracture driving energy reaches zero for fully developed crack, $\textit{i.e.,}$ for $\pf = 1$.

Nevertheless, several researchers have proposed different combinations of degradation functions and locally dissipated fracture energy functions, some of which are presented in Table <ref>.

1cType     1c$g(\pf)$     1cContribution
Quadratic     $(1-\pf)^2$      Bourdin et. al. [2]

Cubic     $s((1-\pf)^3-(1-\pf)^2)+3(1-\pf)^2 - 2(1-\pf)^3$     Borden et. al. [45]

Rational     $\dfrac{(1-\pf)^p}{(1-\pf)^p + a_1\pf + a_1 a_2\pf^2 + a_1 a_2 a_3\pf^3}$     Wu [46]

Stress-degradation functions, popular in the phase-field fracture literature

1cAbbreviated name     1c$w(\pf)$     1c$c_w$
AT1     $\pf$     $8/3$

AT2     $\pf^2$     $2$

PFCZM     $2\pf - \pf^2$     $\pi$

Locally dissipated fracture energy functions, popular in the phase-field fracture literature

For AT1 and AT2 brittle fracture, the quadratic degradation function proposed in [2] is most commonly adopted. However, it is observed that the AT2 model lacks an initial elastic stage. In order to obtain a linear pre-peak response with the AT2 model, researchers opt for a cubic degradation function proposed in [45], with $0 < s \leq 1$, determining the slope of $g(\pf)$ in the undamaged state. For quasi-brittle fracture, [46] developed a rational degradation function with model parameters $a_1$, $a_2$, $a_3$, and $p$ is used. With these parameters, the different traction-separation laws are obtained. The reader is referred to [46] for more on this aspect. In this work, the AT2 model is adopted along with quadratic and cubic degradation functions.

Based on their choices, different traction-separation laws can be implemented, as shown in Table <ref>.

1cTraction-separation law     1c$a_1$     1c$a_2$     1c$a_3$     1c$p$
Linear     $\dfrac{4 E_0 \gc}{\pi \l f_t^2}$     $-0.5$     $0.0$     $2.0$

Exponential     "     $2^{5/3}-3$     $0.0$     $2.5$

Cornellisen     "     $1.3868$     $0.9106$     $2.0$

Traction-separation laws for quasi-brittle PFCZM

Furthermore, the choice of the fracture driving strain energy, $\Psi^f(\strain(\bdisp))$ and the residual, $\Psi^r(\strain(\bdisp))$ is also not unique. Table <ref> presents a few commonly adopted fracture driving and residual energy. The first model proposed in [2, 3] assumes the fracture to driven by the strain energy density $\Psi$, without accounting for tension-compressive asymmetry. Such a model predicts similar fracture in tension and compression. However, most researchers have adopted the notion that fracture is driven by the tensile strain energy density. In this context, [33] adopted a spectral decomposition of the strain energy density function. This yielded the tensile strain energy and the compressive strain energies as the fracture driving and residual energies respectively. In an alternative approach, [46] proposed the energy associated with the maximum principal stress $\sigma_1$ as the fracture driving energy. With reference to Table <ref>, $E$, $\lambda$ and $\mu$ are material constants corresponding to the Young's modulus, Lamé constant and shear modulus respectively. The trace operator is given by $tr$, while $\langle \cdot \rangle_{\pm}$ represents the positive/negative Macaulay brackets. Furthermore, $\strain^{\pm}$ indicates the tensile/compressive strain tensors, obtained through spectral decomposition of the strain tensor.

1cType    1c$\Psi^f$     1c$\Psi^r$     1cContribution
No Split     $\frac{1}{2} \lambda tr^2(\strain) + \mu \, \strain \colon \strain$     $0$     Bourdin et. al. [2, 3]

Spectral     $\frac{1}{2} \lambda \langle tr(\strain) \rangle_+^2 + \mu \, \strain^+ \colon \strain^+$     $\frac{1}{2} \lambda \langle tr(\strain) \rangle_-^2 + \mu \, \strain^- \colon \strain^-$     Miehe et. al. [33]

Rankine     $\frac{1}{2 E} \langle \sigma_1 \rangle^2_+ $     -     Wu [46]

Strain energy density decompositions in phase-field based fracture analysis.

The additive decomposition of the strain energy densitiy into fracture driving energy and a residual energy renders the displacement sub-problem non-linear. In order to preserve a linear displacement sub-problem, [47] proposed a `hybrid' approach. With this approach, the degradation function $g(\pf)$ is applied on the entire strain energy density $\Psi$ instead of $\Psi^f$ in the momentum balance equation. As a consequence, the variational consistency of the problem is lost. Nevertheless, the formulation is consistent w.r.t thermodynamic principles. This formulation, referred to as the `hybrid' phase-field fracture model, is adopted in this work.

§.§ Euler-Lagrange equations

The Euler-Lagrange equations for the phase-field model is obtained upon taking the first variation of the energy functional (<ref>) w.r.t. its solution fields, vector-valued displacements, $\bdisp$ and scalar-valued phase-field, $\pf$. Incorporating the hybrid formulation [47], and with appropriately defined test and trial spaces, the complete problem statement assumes the form:

Find ($\bdisp$, $\pf$) $\in \mathbb{U} \times \mathbb{P}$ with

\begin{align}
E'(\bdisp,\pf;\testbdisp) & = \int_{\dom}^{} \bigg( g(\pf) \dfrac{\partial \Psi(\strain(\bdisp))}{\partial \strain} \bigg) \colon \strain(\testbdisp)  \: \normalfont\text{d}\dom = 0 \hspace{9.15em}\forall \: \testbdisp \in \mathbb{U}^0, \label{eqn:sec:2:P1MomentumBalance}  \\
E'(\bdisp,\pf;\hat{\pf}) & = \int_{\dom}^{} \bigg( g'(\pf) \Psi^f(\strain(\bdisp)) + \dfrac{\gc}{c_w \, \l} w'(\pf) \bigg) (\hat{\pf} - \pf)  \: \normalfont\text{d}\dom \label{eqn:sec:2:P1PfEvolution} \\ \nonumber
& + \int_{\dom}^{} \dfrac{2 \gc \l}{c_w} \bnabla\pf \cdot \bnabla(\hat{\pf} - \pf)  \: \normalfont\text{d}\dom \geq 0 \hspace{11em} \forall \: \hat{\pf} \in \mathbb{P}.
\end{align}

considering the Dirichlet boundary conditions $\bdisp^p$ on $\surfarg{D}{u}$ and $\pf^p$ on $\surfarg{D}{\pf}$, and Neumann boundary condition $\trac{u}$ on $\surfarg{N}{u}$. Moreover, the trial and test spaces are given by

\begin{align}
\mathbb{U} & = \{ \bdisp \in [H^1(\dom)]^{\normalfont\text{dim}}| \bdisp = \bdisp^p \text{ on } \surfarg{D}{u} \},  {\label{eqn:sec:2:P1disp_trialspace}} \\
\mathbb{U}^0 & = \{ \bdisp \in [H^1(\dom)]^{\normalfont\text{dim}}| \bdisp = \mathbf{0} \text{ on } \surfarg{D}{u} \},  {\label{eqn:sec:2:P1disp_testspace}} \\
\mathbb{P} & = \{ \pf \in [H^1(\dom)] | \pf \geq {}^{n}\pf |\pf = \pf^p \text{ on } \surfarg{D}{\pf} \}.  {\label{eqn:sec2:P1pf_space}}
\end{align}

The left superscript $n$ in (<ref>) refers to the previous step in (pseudo) time. $\blacksquare$

Problem <ref> belongs to the variational inequality category (see Equation (<ref> and test/trial space (<ref>)). The treatment of variational inequality is not new in the phase-field fracture model literature. A review of the different approaches ensuring fracture irreversibility is presented in Section <ref>. Adopting the history-variable approach [33], in conjunction with appropriately defined test and trial spaces, the complete problem statement takes the form:

Find ($\bdisp$, $\pf$) $\in \mathbb{U} \times \mathbb{P}$ with

\begin{align}
E'(\bdisp,\pf;\testbdisp) & = \int_{\dom}^{} \bigg( g(\pf) \dfrac{\partial \Psi(\strain(\bdisp))}{\partial \strain}  \bigg) \colon \strain(\testbdisp)  \: \normalfont\text{d}\dom = 0 \hspace{9.15em}\forall \: \testbdisp \in \mathbb{U}^0, \label{eqn:sec:2-2:P2MomentumBalance} \\
E'(\bdisp,\pf;\testpf) & = \int_{\dom}^{} \bigg( g'(\pf) \mathcal{H} + \dfrac{\gc}{c_w \, \l} w'(\pf) \bigg) \testpf  \: \normalfont\text{d}\dom \label{eqn:sec:2-2:P2PfEvolution} \\ \nonumber
& + \int_{\dom}^{} \dfrac{2 \gc \l}{c_w} \bnabla\pf \cdot \bnabla\testpf  \: \normalfont\text{d}\dom = 0 \hspace{9.15em} \forall \: \testpf \in \mathbb{P}^0,
\end{align}

considering the Dirichlet boundary conditions $\bdisp^p$ on $\surfarg{D}{u}$ and $\pf^p$ on $\surfarg{D}{\pf}$. Moreover, the trial and test spaces are given by

\begin{align}
\mathbb{U} & = \{ \bdisp \in [H^1(\dom)]^{\normalfont\text{dim}}| \bdisp = \bdisp^p \text{ on } \surfarg{D}{u} \},  {\label{eqn:sec:2-2:P2disp_trialspace}} \\
\mathbb{U}^0 & = \{ \bdisp \in [H^1(\dom)]^{\normalfont\text{dim}}| \bdisp = \mathbf{0} \text{ on } \surfarg{D}{u} \},  {\label{eqn:sec:2-2:P2disp_testspace}} \\
\mathbb{P} & = \{ \pf \in [H^1(\dom)] | \pf = \pf^p \text{ on } \surfarg{D}{\pf} \},  {\label{eqn:sec2:pf_P2trialspace}} \\
\mathbb{P}^0 & = \{ \pf \in [H^1(\dom)] | \pf = 0 \text{ on } \surfarg{D}{\pf} \}.  {\label{eqn:sec2:P2pf_testspace}}
\end{align}

The history-variable $\mathcal{H}$ is defined as the maximum fracture driving energy $\Psi^f$ over the entire loading history. Mathematically,

\begin{equation}\label{eqn:sec2:P2hist}
\mathcal{H} = \operatorname{max} ({}^{n}\mathcal{H},\Psi^f).
\end{equation}

The left superscript $n$ in (<ref>) refers to the previous step in (pseudo) time. $\blacksquare$

The history-variable approach in Problem <ref> results in a variational equality problem, with relaxed test and trial spaces for the phase-field (cf. Problems <ref> and <ref>).

§ ISOGEOMETRIC ANALYSIS AND DISCRETE EQUATIONS

§.§ Isogeometric analysis and NURBS

Isgeometric analysis (IGA) [48] allows an exact representation of complex geometries, such as such as spheres, ellipsoids, paraboloids and hyperboloids. The representation is carried out using polynomial functions, with the Non-Uniform Rational B-splines (NURBS) most commonly adopted. The smoothness of the the NURBS basis functions is advantageous in problems with multi-faceted surface, that can trigger traction oscillations, when simulated using conventional geometry discretization. Furthermore, IGA also offers the ease in obtaining higher-order continuous basis functions with NURBS. As a consequence, it is appealing for higher-order Partial Differential Equations (PDEs). However, NURBS based modeling is often recognized to have significant flaws in constructing watertight geometries using tensor-product meshes. Also, the scale and scope of refining procedures causes the tensor product structure of NURBS to be inefficient leading to erroneous error estimation and improper implementation of adaptivity algorithms. In the context of phase-field fracture model, NURBS-based simulation was carried out in [49] for the fourth-order phase-field fracture model, albeit without mesh refinement.

The restrictions of the NURBS-based models were mitigated using T-Splines, while keeping the recognizable structure of NURBS algorithms. T-splines alleviate the deficiencies of NURBS by generating a single patch of watertight geometry that can be fine-tuned and coarsened locally. Implementation of T-splines within the framework of IGA has gained a lot of attention. The Bézier extraction [50] of the basis makes it suitable to be efficient integration into existing finite element programs. However, the linear independence of T-splines is not assured in generic T-meshes. The concept of analysis suitable T-splines was proposed in [51], which adopts the essential mathematical entities of NURBS, such as linear independence and partition of unity under certain restrictions on the T-mesh, while giving a highly localized and efficient refining algorithm. As an alternative to T-Splines, PHT-splines were proposed, which are generalization of B-Splines over hierarchical T-meshes. The local refinement algorithm for PHT-splines is extremely efficient, and easy to implement.

In this section, the basic concepts of PHT-splines based IGA is discussed, which is then used as a discretization technique to solve the phase-field fracture problem. In one-dimension, the PHT-spline representation takes the form,

\begin{equation}\label{eq:T_mesh}
\mathbb{T}(\xi) = \sum_{i=1}^{n_{cp}}\mathcal{N}_{i,p}(\xi)\mathbf{P}_i,
\end{equation}

where, $\mathcal{N}_{i,p}(\xi)$ indicates the cubic B-spline basis functions with $C^1$ continuity defined over the knot vector $\Xi$, and $n_{cp}$ is the total number of control points defined over the control mesh used to determine the scaffolding of the geometry. Furthermore, $p$ denotes the order of the polynomial, and $\mathbf{P}_i \in \mathcal{R}^d$ is the set of control points in $d$ dimensions for the B-spline curve with the knot vector $\Xi$. The initial set of knot vectors is denoted by the set of vertices $\Xi^{d}$ corresponding to each spatial direction in the parameter space, $\hat {\Omega} = \left[0,1\right]^d$, and is given by:

\begin{equation}\label{eq:knot_vector}
\Xi^{d} = \{{\xi_0^d, \xi_0^d},{\xi_1^d, \xi_1^d},{\xi_2^d},\ldots,{\xi_{n_i-1}^d, \xi_{n_i-1}^d},{\xi_{n_i}^d, \xi_{n_i}^d}\}.
\end{equation}

Here, ${\xi_i^d}\leq{\xi_{i+1}^d}$, and ${\xi_0^d}\leq{\xi_{1}^d} =0$, ${\xi_{n_i-1}^d}\leq{\xi_{n_i}^d}=1$. Moreover, $n_i =n_{cp} +p+1$ represents the number of elements in each parametric direction. The knot vector is uniform when the distance between the consecutive knots are constant. For each interior vertex $\xi_i$, there are two supporting basis functions, $\left[\xi_{i-1},\xi_{i+1}\right]$. For cubic polynomials, one of the distinguishing aspects of PHT-splines is that they maintain $C^1$ continuity, where the start and end knots are repeated $p+1$ times, while the interior knots are repeated only once. By repeatedly inserting vertices, all the knot spans can be obtained at the same refinement level, hence converting a PHT-spline to a B-spline.

B-splines are non-local, in the sense that, a B-spline typically encompasses more than more element. However, in a finite element framework, a local representation of the B-splines within each element is desired. This local representation of the B-spline is extracted using the Bézier decomposition technique. The Bézier extraction operator is generated using information from the knot vectors and does not rely on the control points or basis functions. Bernstein polynomials have an edge over NURBS basis functions in terms of implementation because the Bernstein basis functions are the same for all elements, as as observed from Figure <ref>. Following this idea, the Bézier extraction allows for the pre-computation of the Bernstein basis on the reference element. During the simulation, the Bernstein basis function can be mapped to each element, with minimal effort.

Bézier extraction operator implemented for a cubic B-spline curve with $\Xi$ = [0,0,0,0,1/3,2/3,1,1,1,1]. Here, a rational Bézier curve, $\mathbf{B}(\bar \xi) = \left\{\beta_{i,p}\right\}_{i=1}^{p+1}$ and its associated control points are $\mathbf{Q} = \left\{Q_i\right\}_{i=1}^{p+1}$ is defined in a reference space $\bar \Omega = [-1,1]$, where $\beta_{i,p}$ are rational Bernstein polynomial basis functions. The B-spline basis functions, $\mathbf{N}(\bar \xi)$ are obtained from the Bernstein basis functions, $\mathbf{B}(\bar \xi)$ using the linear Bézier extraction operator $\mathbf{C}$. In case of obtaining the Bézier control points, the mapping is reversed [52]..

The initial discretization, designated as level 0, is a tensor-product mesh. By splitting the components into $2d$ sub-elements using the cross-insertion approach, the coarse meshes at level $k$ are refined to obtain finer meshes at level $(k+1)$, following the principle of adaptivity. The basis functions for an element on the coarse mesh are replaced by a set of basis functions generated over the refined element in the hierarchical approach. For detailed implementation of the Bézier extraction operator and method to obtain the control points on the refined mesh, the reader is referred to [53].

§.§ Discrete phase-field fracture equations

The IGA proposed for the phase-field fracture model in this work, is similar to the classical Finite Element Analysis (FEA), the only difference being in the basis functions used. For the phase-field fracture model, the Euler-Lagrange equations in Problem <ref> is used as the starting point for discretization. Considering the displacement and the phase-field at control points ($\tilde{\bdisp}_i,\tilde{\pf}_i$) as the fundamental unknowns, the corresponding continuous fields ($\bdisp$,$\pf$) are approximated as,

\begin{equation}\label{eqn:discrete_disp_pf}
\bdisp = \sum\limits_{i = 1}^m {N_i^{\boldsymbol{u}}\tilde{\boldsymbol{u}}_i}\;\;\;\text{ , }
\pf = \sum\limits_{i = 1}^m N_i^{\varphi}\tilde{\varphi}_i.
\end{equation}

In the above equation, $N_i^{\boldsymbol{u}}$ and $N_i^{\pf}$ are the basis functions for the displacement and the phase-field, associated with the $i^{\text{th}}$ control point. The total number of control points is given by $m$. The spatial derivatives of the basis functions $N_i^{\boldsymbol{u}}$ and $N_i^{\pf}$ in a two-dimensional case are computed as,

\begin{equation}\label{eqn:Bmatrices}
\mathbf{B}_i^{\boldsymbol{u}} = \left[ {\begin{array}{*{20}{c}}
{{N_{i,x}}} \\
0 \\
\end{array}}&{\begin{array}{*{20}{c}}
0 \\
{{N_{i,y}}} \\
{{N_{i,x}}} \\
\end{array}}
\end{array}} \right] \text{ , }
\mathbf{B}_i^{\pf}  = \left[ {\begin{array}{*{20}{c}}
{{N_{i,x}}} \\
\end{array}} \right].
\end{equation}

Here, the subscripts $,x$ and $,y$ indicate spatial derivatives in $x$ and $y$ directions respectively. Using (<ref>), the strain $\strain$, and the gradient of the phase-field $\bnabla\pf$ are defined as,

\begin{equation}\label{eqn:strain_gradpf}
\boldsymbol{\epsilon} = \sum\limits_{i = 1}^m {\mathbf{B}_i^{\boldsymbol{u}} \tilde{\boldsymbol{u}}_i}\;\;\;\text{ , }
\bnabla\pf = \sum\limits_{i = 1}^m {\mathbf{B}_i^{\pf} \tilde{\pf}}.
\end{equation}

The discrete phase-field fracture problem is obtained upon inserting (<ref>-<ref>) in the Euler-Lagrange equations from Problem <ref>. Thereafter, (<ref>) and (<ref>) are assumed as the internal forces, and stiffness matrix derived from its derivative. This notation is consistent with [54]. This allows the presentation of the phase-field fracture problem within the incremental iterative framework as:

Compute the solution increment ($\tilde{\Delta\bdisp}$, $\tilde{\Delta\pf}$)$_{i+1}$ in the current iteration $i+1$ using

\begin{equation}\label{dproblem1:discSystem}
\underbrace{
\begin{bmatrix}
\mathbf{K}^{\boldsymbol{uu}} & \mathbf{K}^{\boldsymbol{u}\varphi} \\
\mathbf{K}^{\varphi\boldsymbol{u}} & \mathbf{K}^{\varphi\varphi}
\end{bmatrix}_{i}}_{\text{Stiffness matrix}} \begin{Bmatrix}
\Delta \tilde{\bdisp} \\
\Delta \tilde{\pf}
\end{Bmatrix}_{i+1}  =  \underbrace{\begin{Bmatrix}
\mathbf{f}^{ext,\boldsymbol{u}} \\
\mathbf{f}^{ext,\varphi}
\end{Bmatrix}_{i} - \begin{Bmatrix}
\mathbf{f}^{int,\boldsymbol{u}} \\
\mathbf{f}^{int,\varphi}
\end{Bmatrix}_{i}}_{\text{Residual}},
\end{equation}

and update the solution fields,

\begin{equation}\label{dproblem1:solUpdate}
\begin{Bmatrix}
\tilde{\bdisp} \\
\tilde{\pf}
\end{Bmatrix}_{i+1} = \begin{Bmatrix}
\tilde{\bdisp} \\
\tilde{\pf}
\end{Bmatrix}_{i} +
\begin{Bmatrix}
\Delta \tilde{\bdisp} \\
\Delta \tilde{\pf}
\end{Bmatrix}_{i+1},
\end{equation}

until convergence is achieved. The local element stiffness matrices are given by:

\begin{equation}\label{eq:stiff_comp}
\displaystyle
\begin{split}
\mathbf{K}^{\boldsymbol{uu}} & = \int_{\Omega} \left[\mathbf{B}^{\boldsymbol{u}}\right]^T \bigg( g(\pf) \dfrac{\partial^2 \Psi}{\partial\strain^2} \bigg) \left[\mathbf{B}^{\boldsymbol{u}}\right]\,d\Omega, \\
\mathbf{K}^{\boldsymbol{u}\varphi} & =  \int_{\Omega} \left[\mathbf{B}^{\boldsymbol{u}}\right]^T \bigg( g'(\pf) \dfrac{\partial\Psi}{\partial \strain} \bigg) \left[N^{\varphi} \right] \,d\Omega, \\
\mathbf{K}^{\varphi\boldsymbol{u}} &
= \int_{\Omega} \left[N^{\varphi}\right] \bigg( g'(\pf) \dfrac{\partial \mathcal{H}}{\partial \epsilon} \bigg) \left[\mathbf{B}^{\boldsymbol u}\right]\;d\Omega, \\
\mathbf{K}^{\varphi\varphi} & =
\int_{\Omega}\left\{\left[\mathbf{B}^{\varphi}\right]^T \left( \dfrac{2 \gc \l}{c_w} \right) \left[\mathbf{B}^{\varphi}\right] + \left[N^{\varphi}\right] \left(g''(\pf) \mathcal{H} +\dfrac{\gc}{c_w \l} w''(\varphi) \right) \left[N^{\varphi}\right] \right\} \;d\Omega, \\
\end{split}
\end{equation}

and the local internal force vectors are computed as

\begin{equation}\label{eq:fint_comp}
\begin{split}
\mathbf{f}^{int,\boldsymbol{u}} & = \int_{\Omega} \left[\mathbf{B}^{\boldsymbol{u}}\right]^T \bigg( g(\pf) \dfrac{\partial \Psi}{\partial\strain} \bigg) \,d\Omega,\\
\mathbf{f}^{int,\varphi} & =
\int_{\Omega}\left\{\left[\mathbf{B}^{\varphi}\right]^T \left( \dfrac{2 \gc \l}{c_w} \right) \left[\mathbf{B}^{\varphi}\right] \tilde{\pf} + \left[N^{\varphi}\right]^T \left(g'(\pf) \mathcal{H} +\dfrac{\gc}{c_w \l} w'(\varphi)  \right) \right\} \;d\Omega.
\end{split}
\end{equation}

The external force vectors $f^{ext,\boldsymbol{u}}$ and $f^{ext,\varphi}$ are considered equal to zero. $\blacksquare$

§ A MONOLITHIC SOLUTION TECHNIQUE

§.§ Incremental fracture energy-based arc-length method

Since the inception of the dissipation-based arc-length solver in [55], variants thereof have been utilized for the phase-field based fracture modeling [27, 28, 29]. Motivated by these studies, in particular, the fracture path controlled path following method proposed in [29], a fracture energy-based arc-length method is proposed in this work. From the phase-field fracture energy functional (<ref>), the energy associated with fracture is identified as:

\begin{equation}{\label{eqn:sec2:EFrac}}
\displaystyle G(\pf) = \int_{\dom}^{} \dfrac{\gc}{c_w \l} \left( w(\pf) + \l^2 |\bnabla \pf|^2 \right) \: \text{d}\dom,
\end{equation}

and its incremental form is given by

\begin{equation}{\label{eqn:sec2:EFracIncr}}
\displaystyle \Delta G(\pf,\Delta\pf) = \int_{\dom}^{} \dfrac{\gc}{c_w \l} \left( w'(\pf)\Delta\pf + 2 \l^2 \bnabla \pf \cdot \bnabla \Delta \pf \right) \: \text{d}\dom.
\end{equation}

Thereafter, an arc-length constraint equation $\Lambda(\pf,\Delta\pf)$ is devised, limiting the incremental phase-field fracture energy in (<ref>) to a certain value $\Delta \tau$. Mathematically,

\begin{equation}\label{eqn:arclength_cons}
\Lambda^{(\pf,\Delta\pf)} := \displaystyle \Delta G(\pf,\Delta\pf) - \Delta \tau = 0.
\end{equation}

Within a displacement-controlled loading scenario, an additive decomposition of the displacement is carried out as,

\begin{equation}
\tilde{\bdisp} = \mathbf{C} \tilde{\bdisp}_f + \tilde{\bdisp}_p + \tilde{\lambda} \hat{\bdisp}.
\end{equation}

Here, $\mathbf{C}$ is a constraint matrix [56], $\tilde{\bdisp}_f$, $\tilde{\bdisp}_p$ and $\hat{\bdisp}_p$ are free, prescribed and unit displacements. The load level $\tilde{\lambda}$ acts only on the prescribed displacements $\hat{\bdisp}_p$. With this setup, the discrete problem assumes the form:

Compute the solution increment ($\tilde{\Delta\bdisp}$, $\tilde{\Delta\pf}$, $\Delta\hat{\lambda}$)$_{i+1}$ in the current iteration $i+1$ using

\begin{equation}
\underbrace{
\begin{bmatrix}
\mathbf{K}^{\boldsymbol{uu}} & \mathbf{K}^{\boldsymbol{u}\varphi} & \mathbf{K}^{\boldsymbol{u}\lambda}  \\
\mathbf{K}^{\varphi\boldsymbol{u}} & \mathbf{K}^{\varphi\varphi} & \mathbf{K}^{\varphi\lambda}  \\
\mathbf{0} & \mathbf{K}^{\lambda\varphi} & \mathbf{0}
\end{bmatrix}_{i}}_{\text{Stiffness matrix}} \begin{Bmatrix}
\Delta \tilde{\bdisp} \\
\Delta \tilde{\pf} \\
\Delta \tilde{\lambda}
\end{Bmatrix}_{i+1}  = - \underbrace{\begin{Bmatrix}
\mathbf{f}^{int,\boldsymbol{u}} \\
\mathbf{f}^{int,\varphi} \\
\Lambda^{(\varphi,\Delta\varphi)}
\end{Bmatrix}_{i}}_{\text{Residual}}
\end{equation}

and update the solution fields,

\begin{equation}\label{eq:arcL_solUpdate}
\begin{Bmatrix}
\tilde{\bdisp} \\
\tilde{\pf} \\
\tilde{\lambda}
\end{Bmatrix}_{i+1} = \begin{Bmatrix}
\tilde{\bdisp} \\
\tilde{\pf}  \\
\tilde{\lambda}
\end{Bmatrix}_{i} +
\begin{Bmatrix}
\Delta \tilde{\bdisp} \\
\Delta \tilde{\pf}  \\
\Delta \tilde{\lambda}
\end{Bmatrix}_{i+1}.
\end{equation}

until convergence is achieved. The local element stiffness matrices $\mathbf{K}^{\boldsymbol{uu}}$, $\mathbf{K}^{\boldsymbol{u}\varphi}$, $\mathbf{K}^{\varphi\boldsymbol{u}}$ and $\mathbf{K}^{\varphi\varphi}$, and the local internal force vectors $\mathbf{f}^{int,\boldsymbol{u}}$ and $\mathbf{f}^{int,\varphi}$ remain same as that presented in Problem <ref>. The additional matrices are given by

\begin{equation}\label{eq:addl_stiff_comp}
\displaystyle
\begin{split}
\mathbf{K}^{\boldsymbol{u}\lambda} & = \mathbf{K}^{\boldsymbol{uu}} \: \hat{\bdisp}, \\
\mathbf{K}^{\varphi\lambda} & = \mathbf{K}^{\varphi\boldsymbol{u}} \: \hat{\bdisp}, \\
\mathbf{K}^{\lambda\varphi} & = \dfrac{\gc}{c_w \l}
\int_{\Omega}\left\{ \left[N^{\varphi}\right]^T \big( w''(\pf)\Delta\pf + w'(\pf) \big) + 2 \l^2 \left[\mathbf{B}^{\varphi}\right]^T \cdot \big( \bnabla\Delta\pf + \bnabla\pf  \big)  \right\} \;d\Omega,
\end{split}
\end{equation}

and the additional internal force vector is computed as

\begin{equation}\label{eq:addl_fint_comp}
\begin{split}
\Lambda^{(\varphi,\Delta\varphi)} & = \int_{\dom}^{} \dfrac{\gc}{c_w \l} \left( w'(\pf)\Delta\pf + 2 \l^2 \bnabla \pf \cdot \bnabla \Delta \pf \right) \: \text{d}\dom - \Delta\tau.
\end{split}
\end{equation}

In this work, the phase-field fracture problem is solved using both the conventional displacement-controlled solution scheme (Problem <ref>) and the Arc-length method (Problem <ref>). Algorithm <ref> presented in Appendix <ref> combines both within a single monolithic solution strategy. The time-stepping commences with a displacement-controlled solution scheme and the Newton-Raphson method. Upon convergence, the incremental dissipation $\Delta G$ is computed using (<ref>). If $\Delta G$ is greater than a certain user-defined switch energy, the method is switched to the Arc-length method. Additionally, $\Delta \tau$ is set to the switch energy, and $\Delta \lambda$ is set to zero. Upon achieving convergence in a certain step, an increment in $\Delta \lambda$ is carried out, subjected to a maximum value of $\Delta\tau_{max}$.

§.§ Adaptive under-relaxation scheme

An under-relaxation scheme introduces a scalar parameter $\beta \in (0,1)$ such that the solution field is updated using,

\begin{equation}\label{eq:underRelaxSolUpdate}
\begin{Bmatrix}
\tilde{\bdisp} \\
\tilde{\pf} \\
\tilde{\lambda}
\end{Bmatrix}_{i+1} = \begin{Bmatrix}
\tilde{\bdisp} \\
\tilde{\pf}  \\
\tilde{\lambda}
\end{Bmatrix}_{i} + \beta
\begin{Bmatrix}
\Delta \tilde{\bdisp} \\
\Delta \tilde{\pf}  \\
\Delta \tilde{\lambda}
\end{Bmatrix}_{i+1}.
\end{equation}

When $\beta$ is set to one, the Newton-Raphson method is recovered, otherwise, the method maybe referred to as a modified Newton approach. Under-relaxation schemes are usually robust, however, it may reduce the rate of convergence of the problem [57]. For the phase-field fracture problem, the under-relaxation scheme is adopted to prevent divergence due to an ill behaving stiffness matrix.

Algorithm <ref> in Appendix <ref> presents the under-relaxation adopted in this work. In this scheme, $\beta$ starts with a value one, corresponding to a full Newton-Raphson update. When the convergence is not achieved, the value to $\beta$ is reduced by a factor 1.25. This reduction is carried out twice before performing a reduction in the prescribed incremental dissipation $\Delta\tau$. The motivation behind this is to try the current dissipation step with smaller solution increments within the iterative process of the Newton-Raphson method. Such an approach prevents divergence due to an ill behaving stiffness matrix.

§.§ Adaptive mesh refinement and solution transfer

The phase-field fracture model requires an extremely fine mesh to resolve the crack zone in the computational domain, $\dom$. A sharp crack topology is recovered in the limit as $\l \to 0$ [2]. With fracture length-scales $\l$ very small compared to $|\dom|$, an uniformly refined mesh enormously increases the resources required in terms of computing power and storage. In this work, the novel monolithic solver is integrated with an efficient adaptive mesh refinement (AMR) scheme. The elements of the mesh are chosen for refinement based on a critical threshold value of $\pf$, $\pf_{\text{threshold}}$ [26, 44], which is typically referred as the the refinement indicator. To locally refine the crack path, polynomial splines over hierarchical T-meshes (PHT-splines) are used within the framework of IGA. The PHT-splines possess a very efficient and easy to implement local refinement algorithm. The hierarchical approach replaces the basis functions for an element on the coarse mesh with a set of basis functions constructed over the refined element. The refinement of an element originally defined on the coarse mesh is restricted by a pre-decided maximum number of allowable refinements, to avoid repeated refinements of the elements which are already in the cracked domain.

An adaptive h-refinement scheme is adopted in this work, in which the order of the basis functions remains constant throughout the refinement process. Within the simulation, a series of hierarchical meshes evolve. The mesh during the onset of the simulation is denoted as the base mesh or Level 0 mesh. At any hierarchical level, say `k', some (parent) elements are marked for refinement, following which they are sub-divided into $2^{\text{dim}}$ (children) elements. Once a parent element is refined, it becomes inactive and its children take the place in the computational domain as active elements. Finally, for computational efficiency, the basis functions are computed only for the children elements upon refinement, and not for the entire computational domain.

The field variables are projected from the coarser mesh to the finer mesh for each refined element, repeatedly during the mesh refinement. To avoid re-computation of the problem from the beginning, a variable transfer is required. The discretized variables include the field variables such as $\bdisp$ and $\pf$, computed at the control points that are required to be transferred to the new element. The projection of the field variables from a coarser mesh to a finer mesh is implemented using a similar technique to that described in [53]. For computing the new control points, instead of projecting the geometrical information at the basis vertex, we project $\pf$ on the finer mesh.

§ NUMERICAL EXPERIMENTS

In this section, numerical experiments are carried out on representative fracture problems. These include a tapered bar under tension, a specimen with an eccentric hole under tension, a single edge notched specimen under tensile and later shear loading. For all problems, the geometry, material properties as well as loading conditions are presented in the respective sub-sections. The load-displacement curves and the phase-field fracture topology at the final step of the analysis are also presented therein.

A residual-based convergence criterion is adopted in this work. More specifically, the iterations pertaining to the Newton-Raphson method is terminated, when the norm of the residual is less than $1e-3$.

§.§ Tapered Bar under Tension (TBT)

The first example in the numerical study section is a tapered bar under tensile loading, as shown in Figure <ref>. The bar has dimensions $5$ [mm] in length, $0.75$ [mm] and $2$ [mm] width of fixed end and the prescribed loading edges respectively. The loading in applied in the form of prescribed displacement increment $\tilde{\lambda}\hat{\bdisp}$, where $\hat{\bdisp}$ is a unit load vector and $\tilde{\lambda}$ is the load factor. When the analysis is started, the displacement-controlled approach is adopted and $\hat{\lambda}$ is incremented in steps of $1e-2$ [mm]. Following the switch to the arc-length method, $\hat{\lambda}$ becomes an unknown variable, and is solved with the arc-length constraint equation. The additional parameters required for the analysis are presented in Table <ref>.

(K) at (0,0);

[line width=0.75pt,black] (0,0) to (0,-0.75);
[line width=0.75pt,black] (0,-0.75) to (10,-2);
[line width=0.75pt,black] (10,-2) to (10,2);
[line width=0.75pt,black] (10,2) to (0,0.75);
[line width=0.75pt,black] (0,0.75) to (0,0);

[->,line width=1.5pt,black] (10.2,0.0) to (10.8,0.0);
[ ] at (11.5,0.0) $\tilde{\lambda}\hat{\bdisp}$;

[line width=0.75pt,black] (0,-1.25) to (0,1.25);

[line width=0.75pt,black] (0,1.25) to (-0.25,1.0);
[line width=0.75pt,black] (0,1.00) to (-0.25,0.75);
[line width=0.75pt,black] (0,0.75) to (-0.25,0.50);
[line width=0.75pt,black] (0,0.50) to (-0.25,0.25);
[line width=0.75pt,black] (0,0.25) to (-0.25,0.00);
[line width=0.75pt,black] (0,0.00) to (-0.25,-0.25);
[line width=0.75pt,black] (0,-0.25) to (-0.25,-0.50);
[line width=0.75pt,black] (0,-0.50) to (-0.25,-0.75);
[line width=0.75pt,black] (0,-0.75) to (-0.25,-1.00);
[line width=0.75pt,black] (0,-1.00) to (-0.25,-1.25);
[line width=0.75pt,black] (0,-1.25) to (-0.25,-1.50);

TBT experiment

Parameters     Value
Model     AT2

$\Psi^f$     No Split

$\lambda$     0.0 [MPa]

$\mu$     50.0 [MPa]

$\gc$     1.0 [N/mm]

$\l$     0.25 [mm]

$\pf_{\text{threshold}}$     0.2

$\Delta\tau_{max}$     0.0125 [N]

tableModel parameters

[thick,scale=0.95, every node/.style=scale=0.95]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.45, ymax=25, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/TaperBar/lodi_quad.txt;
[col sep = comma]./Data/TaperBar/lodi_cubics0_01.txt;
[col sep = comma]./Data/TaperBar/lodi_cubics0_1.txt;
[col sep = comma]./Data/TaperBar/lodi_cubics1.txt;

mark size=0.15pt,
x expr=1,
y expr=0
] ;

mark size=0.5pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 0.1$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 1.0$

[inner sep=0pt] () at (0,0)
[inner sep=0pt] () at (-1.3,-3.05)
hide axis,
scale only axis,

colorbar horizontal,
point meta min=0,
point meta max=1,
colorbar style=
xticklabel style = yshift=-0.075cm
[draw=none] coordinates (0,0);
[inner sep=0pt] () at (0,-3.75) $\pf$;

Figure (a) presents the load-displacement curves for the tapered bar under tension. The legend entries correspond to the choice of degradation functions. Figure (b) shows the distribution of the phase-field variable at the final step of the analysis.

Figure <ref> presents the load-displacement curves for the tapered bar under tension. For both, quadratic and cubic degradation function, the specimen exhibits a snap-back behaviour upon localisation. Furthermore, consistent with other studies in the literature, the cubic degradation function demonstrates a linear pre-peak behaviour for small values of $s$ ($<1$). It is also observed that the use of the cubic degradation function yields a higher peak load compared to the quadratic degradation function. The reason for this behaviour can be explained with analytical studies on 1D bar [58]. Finally, in Figure <ref>, the phase-field fracture topology at the final step of the analysis is presented, where the fracture is seen to appear on the fixed end.

§.§ Specimen with an eccentric hole under tension (EH)

A unit square (in mm) embedded with a eccentric hole [59], as shown in Figure <ref>, is considered for numerical study. The hole has a radius 0.2 [mm] and is centered at (0.6,0.0). A quasi-static loading is applied at the top boundary in the form of prescribed displacement increment $\tilde{\lambda}\hat{\bdisp}$, where $\hat{\bdisp}$ is a unit load vector and $\tilde{\lambda}$ is the load factor. The loading in applied in the form of prescribed displacement increment $\tilde{\lambda}\hat{\bdisp}$, where $\hat{\bdisp}$ is a unit load vector and $\tilde{\lambda}$ is the load factor. When the analysis is started, the displacement-control approach is adopted and $\hat{\lambda}$ is incremented in steps of $1e-4$ [mm]. Following the switch to the arc-length method, $\hat{\lambda}$ becomes an unknown variable, and is solved with the arc-length constraint equation. The bottom boundary remains fixed. The additional parameters required for the analysis are presented in Table <ref>.

(K) at (0,0);

[fill=black!20] (-2.5,-2.5) rectangle (2.5,2.5);
[fill=white] (0.5,0) circle (1.0);

[line width=0.75pt,black] (2.5,2.75) to (-2.5,2.75);
[->,line width=1.5pt,black] (0.0,3.0) to (0.0,3.5);
[ ] at (0,3.85) $\tilde{\lambda}\hat{\bdisp}$;

[line width=0.75pt,black] (-2.5,-2.5) to (-2.75,-2.75);
[line width=0.75pt,black] (-2.25,-2.5) to (-2.5,-2.75);
[line width=0.75pt,black] (-2.0,-2.5) to (-2.25,-2.75);
[line width=0.75pt,black] (-1.75,-2.5) to (-2.0,-2.75);
[line width=0.75pt,black] (-1.5,-2.5) to (-1.75,-2.75);
[line width=0.75pt,black] (-1.25,-2.5) to (-1.5,-2.75);
[line width=0.75pt,black] (-1.0,-2.5) to (-1.25,-2.75);
[line width=0.75pt,black] (-0.75,-2.5) to (-1.0,-2.75);
[line width=0.75pt,black] (-0.5,-2.5) to (-0.75,-2.75);
[line width=0.75pt,black] (-0.25,-2.5) to (-0.5,-2.75);
[line width=0.75pt,black] (0.0,-2.5) to (-0.25,-2.75);
[line width=0.75pt,black] (0.25,-2.5) to (0.0,-2.75);
[line width=0.75pt,black] (0.5,-2.5) to (0.25,-2.75);
[line width=0.75pt,black] (0.75,-2.5) to (0.5,-2.75);
[line width=0.75pt,black] (1.0,-2.5) to (0.75,-2.75);
[line width=0.75pt,black] (1.25,-2.5) to (1.0,-2.75);
[line width=0.75pt,black] (1.5,-2.5) to (1.25,-2.75);
[line width=0.75pt,black] (1.75,-2.5) to (1.5,-2.75);
[line width=0.75pt,black] (2.0,-2.5) to (1.75,-2.75);
[line width=0.75pt,black] (2.25,-2.5) to (2.0,-2.75);
[line width=0.75pt,black] (2.5,-2.5) to (2.25,-2.75);

EH experiment

Parameters     Value
Model     AT2

$\Psi^f$     Rankine

$\lambda$     121.154 [GPa]

$\mu$     80.769 [GPa]

$\gc$     2700 [N/m]

$\l$     2e-2 [mm]

$\pf_{\text{threshold}}$     0.2

$\Delta\tau_{max}$     0.05 [N/mm$^2$]

tableModel parameters

[thick,scale=0.95, every node/.style=scale=0.95]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.008, ymax=1600, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/EH/lodi_quad.txt;

[col sep = comma]./Data/EH/lodi_cubics0_1.txt;
[col sep = comma]./Data/EH/lodi_cubics1.txt;

mark size=0.25pt,
x expr=1,
y expr=0
] ;

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 0.1$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 1.0$

[inner sep=0pt] () at (0,0)
[inner sep=0pt] () at (-1.3,-3.05)
hide axis,
scale only axis,

colorbar horizontal,
point meta min=0,
point meta max=1,
colorbar style=
xticklabel style = yshift=-0.075cm
[draw=none] coordinates (0,0);
[inner sep=0pt] () at (0,-3.75) $\pf$;

Figure (a) presents the load-displacement curves for the single edge notched specimen under shear. The legend entries correspond to the choice of degradation functions. Figure (b) shows the distribution of the phase-field variable at the final step of the analysis.

Figure <ref> presents the load-displacement curves for the unit square specimen with an eccentric hole under tension. The different curves correspond to the choice of the degradation function, quadratic and cubic. It is observed that the specimen exhibits a linear pre-peak behaviour with the cubic degradation function for $s<1$. This linear stage is missing for the quadratic degradation function and the cubic degradation function with $s=1$. However, irrespective of the choice of the degradation function, two snap-back behaviours are observed. The first one occurring at the onset of the localization from the hole towards the right of the specimen. The second snap-back behaviour corresponds to the onset of the localization from the hole towards the left edge of the specimen. Next, in Figure <ref>, the phase-field fracture topology at the final step of the analysis is presented. The fracture topology similar to those observed in the literature [59]. The refined meshes corresponding to the different stages in the evolution of the phase-field is shown in Figure <ref>.

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

hide axis,
scale only axis,
point meta min=0.0,
point meta max=1.0,

colorbar right,
colorbar sampled,
colorbar style=
separate axis lines,
every colorbar/.append style=
height=/pgfplots/parent axis height,
ytick=0,0.5,1,yticklabels=0,0.5  $\pf$,1

[draw=none] coordinates (0,0);

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

Figures (a-d) shows the evolution of the phase-field, and the corresponding refined meshes are shown in Figures (e-h) for the eccentric hole. The colors in the latter figures represent the different patches in the IGA context.

§.§ Single Edge Notched specimen under Tension (SENT)

A unit square (in mm) embedded with a horizontal notch, midway along the height is considered, as shown in Figure <ref>. The length of the notch is equal to half of the edge length of the plate (shown in red). 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 $\tilde{\lambda}\hat{\bdisp}$, where $\hat{\bdisp}$ is a unit load vector and $\tilde{\lambda}$ is the load factor. When the analysis is started, the displacement-control approach is adopted and $\hat{\lambda}$ is incremented in steps of $1e-4$. Following the switch to the arc-length method, $\hat{\lambda}$ becomes an unknown variable, and is solved with the arc-length constraint equation. Furthermore, the bottom boundary remains fixed. The additional parameters required for the analysis are presented in Table <ref>.

(K) at (0,0);

[line width=0.75pt,black] (-2.5,-2.5) to (2.5,-2.5);
[line width=0.75pt,black] (2.5,-2.5) to (2.5,2.5);
[line width=0.75pt,black] (2.5,2.5) to (-2.5,2.5);
[line width=0.75pt,black] (-2.5,2.5) to (-2.5,-2.5);
[line width=1.5pt,red] (-2.5,0) to (0,0);

[line width=0.75pt,black] (2.5,2.75) to (-2.5,2.75);
[->,line width=1.5pt,black] (0.0,3.0) to (0.0,3.5);
[ ] at (0,3.85) $\tilde{\lambda}\hat{\bdisp}$;

[line width=0.75pt,black] (-2.5,-2.5) to (-2.75,-2.75);
[line width=0.75pt,black] (-2.25,-2.5) to (-2.5,-2.75);
[line width=0.75pt,black] (-2.0,-2.5) to (-2.25,-2.75);
[line width=0.75pt,black] (-1.75,-2.5) to (-2.0,-2.75);
[line width=0.75pt,black] (-1.5,-2.5) to (-1.75,-2.75);
[line width=0.75pt,black] (-1.25,-2.5) to (-1.5,-2.75);
[line width=0.75pt,black] (-1.0,-2.5) to (-1.25,-2.75);
[line width=0.75pt,black] (-0.75,-2.5) to (-1.0,-2.75);
[line width=0.75pt,black] (-0.5,-2.5) to (-0.75,-2.75);
[line width=0.75pt,black] (-0.25,-2.5) to (-0.5,-2.75);
[line width=0.75pt,black] (0.0,-2.5) to (-0.25,-2.75);
[line width=0.75pt,black] (0.25,-2.5) to (0.0,-2.75);
[line width=0.75pt,black] (0.5,-2.5) to (0.25,-2.75);
[line width=0.75pt,black] (0.75,-2.5) to (0.5,-2.75);
[line width=0.75pt,black] (1.0,-2.5) to (0.75,-2.75);
[line width=0.75pt,black] (1.25,-2.5) to (1.0,-2.75);
[line width=0.75pt,black] (1.5,-2.5) to (1.25,-2.75);
[line width=0.75pt,black] (1.75,-2.5) to (1.5,-2.75);
[line width=0.75pt,black] (2.0,-2.5) to (1.75,-2.75);
[line width=0.75pt,black] (2.25,-2.5) to (2.0,-2.75);
[line width=0.75pt,black] (2.5,-2.5) to (2.25,-2.75);

SENT experiment

Parameters     Value
Model     AT2

$\Psi^f$     No Split

$\lambda$     121.154 [GPa]

$\mu$     80.769 [GPa]

$\gc$     2700 [N/m]

$\l$     2e-2 [mm]

$\pf_{\text{threshold}}$     0.2

$\Delta\tau_{max}$     0.025 [N]

tableModel parameters

Figure <ref> presents the load-displacement curves for the single edge notched specimen under tension. The different curves correspond to the choice of the degradation function, quadratic and cubic. Similar to the previous section, it is observed that the specimen exhibits a linear pre-peak behaviour with the cubic degradation function for $s<1$. The quadratic degradation function does not exhibit a linear pre-peak behaviour. Moreover, beyond the first snap-back behaviour, the post-peak branches of all curves are similar. Next, in Figure <ref>, the phase-field fracture topology at the final step of the analysis is presented. The fracture topology is consistent with those observed in the literature, for instance, [33]. The refined meshes corresponding to the different stages in the evolution of the phase-field is shown in Figure <ref>.

[thick,scale=0.95, every node/.style=scale=0.95]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.0075, ymax=1300, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/SENT/lodi_quad.txt;
[col sep = comma]./Data/SENT/lodi_cubics0_01.txt;
[col sep = comma]./Data/SENT/lodi_cubics0_1.txt;
[col sep = comma]./Data/SENT/lodi_cubics1.txt;

mark size=0.25pt,
x expr=1,
y expr=0
] ;

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 0.1$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 1.0$

[inner sep=0pt] () at (0,0)
[inner sep=0pt] () at (-1.3,-3.05)
hide axis,
scale only axis,

colorbar horizontal,
point meta min=0,
point meta max=1,
colorbar style=
xticklabel style = yshift=-0.075cm
[draw=none] coordinates (0,0);
[inner sep=0pt] () at (0,-3.75) $\pf$;

Figure (a) presents the load-displacement curves for the single edge notched specimen under tension. The legend entries correspond to the choice of degradation functions. Figure (b) shows the distribution of the phase-field variable at the final step of the analysis.

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

hide axis,
scale only axis,
point meta min=0.0,
point meta max=1.0,

colorbar right,
colorbar sampled,
colorbar style=
separate axis lines,
every colorbar/.append style=
height=/pgfplots/parent axis height,
ytick=0,0.5,1,yticklabels=0,0.5  $\pf$,1

[draw=none] coordinates (0,0);

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

Figures (a-d) shows the evolution of the phase-field, and the corresponding refined meshes are shown in Figures (e-h) for the single edge notched specimen under tension. The colors in the latter figures represent the different patches in the IGA context.

Furthermore, sensitivity studies are carried out w.r.t to the choice of the prescribed maximum dissipation $\Delta\tau_{max}$ and the phase-field threshold for mesh refinement. The results in Appendix <ref> present the $\Delta\tau_{max}$ for which similar load-displacement curves are obtained. Next, it is also observed that the adaptive mesh refinement technique with different values of $\pf_{th.}$ yield similar load-displacement curve as that obtained on a fixed mesh. Finally, in Appendix <ref>, the proposed monolithic solver is compared with conventional alternate minimization solver and the quasi-Newton Raphson method [26]. For the SENT problem, it is observed that all methods/solvers yield similar peak load and pre-peak behaviour. However, the proposed monolithic solver augmented with the arc-length method captures snap-back behaviours, which is not possible with alternative minimization and quasi-Newton Raphson method.

§.§ Single Edge Notched specimen under Shear (SENS)

The shear test is carried out on the single edge notched specimen by loading in the horizontal direction, as shown in Figure <ref>. The model parameters remain same as presented in Table <ref>. Similar to the SENT model, a quasi-static loading is applied at the top boundary in the form of prescribed displacement increment $\tilde{\lambda}\hat{\bdisp}$, where $\hat{\bdisp}$ is a unit load vector and $\tilde{\lambda}$ is the load factor. The loading in applied in the form of prescribed displacement increment $\tilde{\lambda}\hat{\bdisp}$, where $\hat{\bdisp}$ is a unit load vector and $\tilde{\lambda}$ is the load factor. When the analysis is started, the displacement-control approach is adopted and $\hat{\lambda}$ is incremented in steps of $5e-4$ [mm]. Following the switch to the arc-length method, $\hat{\lambda}$ becomes an unknown variable, and is solved with the arc-length constraint equation. The bottom boundary of the specimen is fixed, while the left and the right edges are restricted in the vertical direction.

(K) at (0,0);

[line width=0.75pt,black] (-2.5,-2.5) to (2.5,-2.5);
[line width=0.75pt,black] (2.5,-2.5) to (2.5,2.5);
[line width=0.75pt,black] (2.5,2.5) to (-2.5,2.5);
[line width=0.75pt,black] (-2.5,2.5) to (-2.5,-2.5);
[line width=1.5pt,red] (-2.5,0) to (0,0);

[line width=0.75pt,black] (2.5,2.75) to (-2.5,2.75);
[->,line width=1.5pt,black] (0.1,3.25) to (0.8,3.25);
[ ] at (-0.5,3.25) $\tilde{\lambda}\hat{\bdisp}$;

[line width=0.75pt,black] (-2.5,-2.5) to (-2.75,-2.75);
[line width=0.75pt,black] (-2.25,-2.5) to (-2.5,-2.75);
[line width=0.75pt,black] (-2.0,-2.5) to (-2.25,-2.75);
[line width=0.75pt,black] (-1.75,-2.5) to (-2.0,-2.75);
[line width=0.75pt,black] (-1.5,-2.5) to (-1.75,-2.75);
[line width=0.75pt,black] (-1.25,-2.5) to (-1.5,-2.75);
[line width=0.75pt,black] (-1.0,-2.5) to (-1.25,-2.75);
[line width=0.75pt,black] (-0.75,-2.5) to (-1.0,-2.75);
[line width=0.75pt,black] (-0.5,-2.5) to (-0.75,-2.75);
[line width=0.75pt,black] (-0.25,-2.5) to (-0.5,-2.75);
[line width=0.75pt,black] (0.0,-2.5) to (-0.25,-2.75);
[line width=0.75pt,black] (0.25,-2.5) to (0.0,-2.75);
[line width=0.75pt,black] (0.5,-2.5) to (0.25,-2.75);
[line width=0.75pt,black] (0.75,-2.5) to (0.5,-2.75);
[line width=0.75pt,black] (1.0,-2.5) to (0.75,-2.75);
[line width=0.75pt,black] (1.25,-2.5) to (1.0,-2.75);
[line width=0.75pt,black] (1.5,-2.5) to (1.25,-2.75);
[line width=0.75pt,black] (1.75,-2.5) to (1.5,-2.75);
[line width=0.75pt,black] (2.0,-2.5) to (1.75,-2.75);
[line width=0.75pt,black] (2.25,-2.5) to (2.0,-2.75);
[line width=0.75pt,black] (2.5,-2.5) to (2.25,-2.75);

[fill=gray!50] (-2.75,1.25) – (-2.25,1.25) – (-2.5,1.75)– (-2.75,1.25);
[fill=black!75] (-2.5,1.05) circle (0.2);
[line width=1pt,black] (-2.25,0.8) to (-2.75,0.8);

[fill=gray!50] (-2.75,-1.5) – (-2.25,-1.5) – (-2.5,-1.)– (-2.75,-1.5);
[fill=black!75] (-2.5,-1.7) circle (0.2);
[line width=1pt,black] (-2.25,-1.95) to (-2.75,-1.95);

[fill=gray!50] (2.75,1.25) – (2.25,1.25) – (2.5,1.75)– (2.75,1.25);
[fill=black!75] (2.5,1.05) circle (0.2);
[line width=1pt,black] (2.25,0.8) to (2.75,0.8);

[fill=gray!50] (2.75,-1.5) – (2.25,-1.5) – (2.5,-1.)– (2.75,-1.5);
[fill=black!75] (2.5,-1.7) circle (0.2);
[line width=1pt,black] (2.25,-1.95) to (2.75,-1.95);

SENS experiment

Parameters     Value
Model     AT2

$\Psi^f$     Rankine

$\lambda$     121.154 [GPa]

$\mu$     80.769 [GPa]

$\gc$     2700 [N/m]

$\l$     2e-2 [mm]

$\pf_{\text{threshold}}$     0.1

$\Delta\tau_{max}$     0.025 [N]

tableModel parameters

Figure <ref> presents the load-displacement curves for the single edge notched specimen under shear.  The different curves correspond to the choice of the degradation function, quadratic and cubic. It is observed that the specimen exhibits a linear pre-peak behaviour with the cubic degradation function for $s<1$. This linear stage is not exhibited by the quadratic degradation function and the cubic degradation function with $s=1$. However, irrespective of the choice of the degradation function, two snap-back behaviours are observed. The first one occurring at the onset of the localization whereas the second snap-back behaviour appears when the crack has reached the bottom edge. Next, in Figure <ref>, the phase-field fracture topology at the final step of the analysis is presented. The fracture topology differs from that presented in [33], and the reason lies in the choice of fracture driving $\Psi^f$. In [33], the fracture is driven by the tensile strain energy obtained through spectral decomposition, whereas in this work, the Rankine criterion [46] is adopted. The refined meshes corresponding to the different stages in the evolution of the phase-field is shown in Figure <ref>.

[thick,scale=0.95, every node/.style=scale=0.95]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.014, ymax=800, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/SENS/lodi_quadn.txt;

[col sep = comma]./Data/SENS/lodi_cubics0_1n.txt;
[col sep = comma]./Data/SENS/lodi_cubics1n.txt;

mark size=0.25pt,
x expr=1,
y expr=0
] ;

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 0.1$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 1.0$

[inner sep=0pt] () at (0,0)
[inner sep=0pt] () at (-1.3,-3.05)
hide axis,
scale only axis,

colorbar horizontal,
point meta min=0,
point meta max=1,
colorbar style=
xticklabel style = yshift=-0.075cm
[draw=none] coordinates (0,0);
[inner sep=0pt] () at (0,-3.75) $\pf$;

Figure (a) presents the load-displacement curves for the single edge notched specimen under shear. The legend entries correspond to the choice of degradation functions. Figure (b) shows the distribution of the phase-field variable at the final step of the analysis.

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

hide axis,
scale only axis,
point meta min=0.0,
point meta max=1.0,

colorbar right,
colorbar sampled,
colorbar style=
separate axis lines,
every colorbar/.append style=
height=/pgfplots/parent axis height,
ytick=0,0.5,1,yticklabels=0,0.5  $\pf$,1

[draw=none] coordinates (0,0);

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

[inner sep=0pt] () at (0,0)

Figures (a-d) shows the evolution of the phase-field, and the corresponding refined meshes are shown in Figures (e-h) for the single edge notched specimen under shear. The colors in the latter figures represent the different patches in the IGA context.

§ CONCLUDING REMARKS

In this work, we have proposed a robust monolithic solver to accurately forecast the material behavior, which is essential for predicting damage progression and estimating possible failure paths. The literature on variational phase-field based fracture modeling still lacks a reliable, efficient, and simple monolithic solver capable of capturing the pre- and post-peak behaviours accurately. The proposed fully monolithic solver adopts a fracture energy-based arc-length method and an adaptive under-relaxation technique to bridge this gap. The proposed solver utilizes an adaptive mesh refinement scheme using polynomial splines over hierarchical T-meshes (PHT-splines) within the framework of IGA. The PHT-splines possess a very efficient and easy to implement local refinement algorithm, which makes it a right choice for capturing quantities of local interest. The combination of the proposed solver with an adaptive mesh refinement technique could facilitate the application of this approach to more complex structures and with sophisticated constitutive laws. Through the four test cases presented in this work, the crack is allowed to nucleate on its own, and also the solver captures the post-peak snap back effects which is not possible with the alternate minimization solver[3, 6] and the quasi-monolithic scheme [26].

Further extensions of this work may include complex multiphysics problems (e.g., porous media, corrosion), dynamic fracture, and the unified phase-field fracture model [46] for quasi brittle fracture. Also, plasticity models could be incorporated for ductile fracture.

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.

Acknowledgment: The first author (R.B) is thankful to Elias Börjesson for many helpful discussion on the arc-length method. The financial support from the Swedish Research Council for Sustainable Development (FORMAS) under Grant 2018-01249 and the Swedish Research Council (VR) under Grant 2017-05192 is gratefully acknowledged.

[1]
G.A. Francfort and J.-J. Marigo.
Revisiting brittle fracture as an energy minimization problem.
Journal of the Mechanics and Physics of Solids,
46(8):1319–1342, 8 1998.

[2]
B. Bourdin, G.A. Francfort, and J-J. Marigo.
Numerical experiments in revisited brittle fracture.
Journal of the Mechanics and Physics of Solids, 48(4):797 –
826, 2000.

[3]
B. Bourdin.
Numerical implementation of the variational formulation for
quasi-static brittle fracture.
Interfaces and Free Boundaries, 9:411–430, 2007.

[4]
David Bryant Mumford and Jayant Shah.
Optimal approximations by piecewise smooth functions and associated
variational problems.
Communications on pure and applied mathematics, 1989.

[5]
Stéphane PA Bordas, Sundararajan Natarajan, Pierre Kerfriden,
Charles Edward Augarde, D Roy Mahapatra, Timon Rabczuk, and Stefano Dal Pont.
On the performance of strain smoothing for quadratic and enriched
finite element approximations (xfem/gfem/pufem).
International Journal for Numerical Methods in Engineering,
86(4-5):637–666, 2011.

[6]
Christian Miehe, Fabian Welschinger, and Martina Hofacker.
Thermodynamically consistent phase-field models of fracture:
Variational principles and multi-field fe implementations.
International journal for numerical methods in engineering,
83(10):1273–1311, 2010.

[7]
Christian Miehe, Lisa-Marie Schänzel, and Heike Ulmer.
Phase field modeling of fracture in multi-physics problems. part i.
balance of crack surface and failure criteria for brittle crack propagation
in thermo-elastic solids.
Computer Methods in Applied Mechanics and Engineering,
294:449–485, 2015.

[8]
M. Ambati, T. Gerasimov, and L. De Lorenzis.
Phase-field modeling of ductile fracture.
Computational Mechanics, 55(5):1017–1040, May 2015.

[9]
C. Miehe, M. Hofacker, L.-M. Schänzel, and F. Aldakheel.
Phase field modeling of fracture in multi-physics problems. part ii.
coupled brittle-to-ductile failure criteria and crack propagation in
thermo-elastic–plastic solids.
Computer Methods in Applied Mechanics and Engineering,
294:486–522, 2015.

[10]
S. Teichtmeister, D. Kienle, F. Aldakheel, and M.-A. Keip.
Phase field modeling of fracture in anisotropic brittle solids.
International Journal of Non-Linear Mechanics, 97:1–21, 2017.

[11]
Jeremy Bleyer and Roberto Alessi.
Phase-field modeling of anisotropic brittle fracture including
several damage mechanisms.
Computer Methods in Applied Mechanics and Engineering,
336:213–236, 2018.

[12]
Zachary A. Wilson and Chad M. Landis.
Phase-field modeling of hydraulic fracture.
Journal of the Mechanics and Physics of Solids, 96:264–290,

[13]
Chukwudi Chukwudozie, Blaise Bourdin, and Keita Yoshioka.
A variational phase-field model for hydraulic fracturing in porous
Computer Methods in Applied Mechanics and Engineering,
347:957–982, 2019.

[14]
Tuanny Cajuhi, Lorenzo Sanavia, and Laura De Lorenzis.
Phase-field modeling of fracture in variably saturated porous media.
Computational Mechanics, 61(3):299–318, 2018.

[15]
Yousef Heider and WaiChing Sun.
A phase field framework for capillary-induced fracture in unsaturated
porous media: Drying-induced vs. hydraulic cracking.
Computer Methods in Applied Mechanics and Engineering,
359:112647, 2020.

[16]
RU Patil, BK Mishra, and IV Singh.
An adaptive multiscale phase field method for brittle fracture.
Computer Methods in Applied Mechanics and Engineering,
329:254–288, 2018.

[17]
Lam H Nguyen and Dominik Schillinger.
The multiscale finite element method for nonlinear continuum
localization problems at full fine-scale fidelity, illustrated through
phase-field fracture and plasticity.
Journal of Computational Physics, 396:129–160, 2019.

[18]
Savvas P Triantafyllou and Emmanouil G Kakouris.
A generalized phase field multiscale finite element method for
brittle fracture.
International Journal for Numerical Methods in Engineering,
121(9):1915–1945, 2020.

[19]
Tymofiy Gerasimov, Nima Noii, Olivier Allix, and Laura De Lorenzis.
A non-intrusive global/local approach applied to phase-field modeling
of brittle fracture.
Advanced modeling and simulation in engineering sciences,
5(1):1–30, 2018.

[20]
Nima Noii, Fadi Aldakheel, Thomas Wick, and Peter Wriggers.
An adaptive global–local approach for phase-field modeling of
anisotropic brittle fracture.
Computer Methods in Applied Mechanics and Engineering,
361:112744, 2020.

[21]
Francesca Fantoni, Andrea Bacigalupo, Marco Paggi, and Josè Reinoso.
A phase field approach for damage propagation in periodic
microstructured materials.
International Journal of Fracture, 223(1):53–76, 2020.

[22]
Bang He, Louis Schuler, and Pania Newell.
A numerical-homogenization based phase-field fracture modeling of
linear elastic heterogeneous porous media.
Computational Materials Science, 176:109519, 2020.

[23]
Ritukesh Bharali, Fredrik Larsson, and Ralf Jänicke.
Computational homogenisation of phase-field fracture.
European Journal of Mechanics - A/Solids, 88:104247, 2021.

[24]
Frédéric Feyel and Jean-Louis Chaboche.
Fe2 multiscale approach for modelling the elastoviscoplastic
behaviour of long fibre sic/ti composite materials.
Computer methods in applied mechanics and engineering,
183(3-4):309–330, 2000.

[25]
T. Gerasimov and L. De Lorenzis.
A line search assisted monolithic approach for phase-field computing
of brittle fracture.
Computer Methods in Applied Mechanics and Engineering,
312:276–303, 12 2016.

[26]
Timo Heister, Mary F. Wheeler, and Thomas Wick.
A primal-dual active set method and predictor-corrector mesh
adaptivity for computing fracture propagation using a phase-field approach.
Computer Methods in Applied Mechanics and Engineering, 2015.

[27]
Julien Vignollet, Stefan May, René De Borst, and Clemens V Verhoosel.
Phase-field models for brittle and cohesive fracture.
Meccanica, 49(11):2587–2601, 2014.

[28]
Stefan May, Julien Vignollet, and Rene De Borst.
A numerical assessment of phase-field models for brittle and cohesive
fracture: $\gamma$-convergence and stress oscillations.
European Journal of Mechanics-A/Solids, 52:72–84, 2015.

[29]
N Singh, CV Verhoosel, R De Borst, and EH Van Brummelen.
A fracture-controlled path-following technique for phase-field
modeling of brittle fracture.
Finite Elements in Analysis and Design, 113:14–29, 2016.

[30]
Thomas Wick.
Modified newton methods for solving fully monolithic phase-field
quasi-static brittle fracture propagation.
Computer Methods in Applied Mechanics and Engineering,
325:577–611, 2017.

[31]
Thomas Wick.
An Error-Oriented Newton/Inexact Augmented Lagrangian Approach for
Fully Monolithic Phase-Field Fracture Propagation.
SIAM Journal on Scientific Computing, 39(4):B589–B617, 2017.

[32]
Alena Kopaničáková and Rolf Krause.
A recursive multilevel trust region method with application to fully
monolithic phase-field models of brittle fracture.
Computer Methods in Applied Mechanics and Engineering,
360:112720, 2020.

[33]
Christian Miehe, Martina Hofacker, and Fabian Welschinger.
A phase field model for rate-independent crack propagation: Robust
algorithmic implementation based on operator splits.
Computer Methods in Applied Mechanics and Engineering,
199(45-48):2765–2778, 2010.

[34]
Tymofiy Gerasimov and Laura De Lorenzis.
On penalization in variational phase-field models of brittle
Computer Methods in Applied Mechanics and Engineering,
354:990–1026, 2019.

[35]
C Miehe, S Teichtmeister, and F Aldakheel.
Phase-field modelling of ductile fracture: a variational
gradient-extended plasticity-damage theory and its micromorphic
Philosophical Transactions of the Royal Society A: Mathematical,
Physical and Engineering Sciences, 374(2066):20150170, 2016.

[36]
Ritukesh Bharali, Fredrik Larsson, and Ralf Jänicke.
A micromorphic phase-field model for brittle fracture.
Preprint available at

[37]
Ritukesh Bharali, Fredrik Larsson, and Ralf Jänicke.
Phase-field fracture irreversibility using the slack variable
Preprint available at

[38]
Laura De Lorenzis and Corrado Maurini.
Nucleation under multi-axial loading in variational phase-field
models of brittle fracture.
International Journal of Fracture, 06 2021.

[39]
Somdatta Goswami, Cosmin Anitescu, and Timon Rabczuk.
Adaptive fourth-order phase field analysis for brittle fracture.
Computer Methods in Applied Mechanics and Engineering,
361:112808, 2020.

[40]
Markus Klinsmann, Daniele Rosato, Marc Kamlah, and Robert M McMeeking.
An assessment of the phase field formulation for crack growth.
Computer Methods in Applied Mechanics and Engineering,
294:313–330, 2015.

[41]
Chintan Jansari, K Kannan, RK Annabattula, S Natarajan, et al.
Adaptive phase field method for quasi-static brittle fracture using a
recovery based error indicator and quadtree decomposition.
Engineering Fracture Mechanics, 220:106599, 2019.

[42]
Thomas Wick.
Goal functional evaluations for phase-field fracture using pu-based
dwr mesh adaptivity.
Computational Mechanics, 57(6):1017–1035, 2016.

[43]
Sindhu Nagaraja, Mohamed Elhaddad, Marreddy Ambati, Stefan Kollmannsberger,
Laura De Lorenzis, and Ernst Rank.
Phase-field modeling of brittle fracture with multi-level hp-fem and
the finite cell method.
Computational mechanics, 63(6):1283–1300, 2019.

[44]
Somdatta Goswami, Cosmin Anitescu, and Timon Rabczuk.
Adaptive phase field analysis with dual hierarchical meshes for
brittle fracture.
Engineering Fracture Mechanics, 218:106608, 2019.

[45]
Michael J Borden, Thomas JR Hughes, Chad M Landis, Amin Anvari, and Isaac J
A phase-field formulation for fracture in ductile materials: Finite
deformation balance law derivation, plastic degradation, and stress
triaxiality effects.
Computer Methods in Applied Mechanics and Engineering,
312:130–166, 2016.

[46]
Jian-Ying Wu.
A unified phase-field theory for the mechanics of damage and
quasi-brittle failure.
Journal of the Mechanics and Physics of Solids, 103:72–99,

[47]
Marreddy Ambati, Tymofiy Gerasimov, and Laura De Lorenzis.
A review on phase-field models of brittle fracture and a new fast
hybrid formulation.
Computational Mechanics, 55(2):383–405, 02 2015.

[48]
Thomas JR Hughes, John A Cottrell, and Yuri Bazilevs.
Isogeometric analysis: Cad, finite elements, nurbs, exact geometry
and mesh refinement.
Computer methods in applied mechanics and engineering,
194(39-41):4135–4195, 2005.

[49]
Michael J. Borden, Thomas J.R. Hughes, Chad M. Landis, and Clemens V.
A higher-order phase-field model for brittle fracture: Formulation
and analysis within the isogeometric analysis framework.
Computer Methods in Applied Mechanics and Engineering,
273:100–118, 2014.

[50]
Michael A Scott, Michael J Borden, Clemens V Verhoosel, Thomas W Sederberg, and
Thomas JR Hughes.
Isogeometric finite element data structures based on bézier
extraction of t-splines.
International Journal for Numerical Methods in Engineering,
88(2):126–156, 2011.

[51]
Annalisa Buffa, Durkbin Cho, and Giancarlo Sangalli.
Linear independence of the t-spline blending functions associated
with some particular t-meshes.
Computer Methods in Applied Mechanics and Engineering,
199(23-24):1437–1445, 2010.

[52]
Somdatta Goswami.
Phase field modeling of fracture with isogeometric analysis and
machine learning methods.
PhD thesis, Bauhaus University-Weimar, Germany, 2021.

[53]
P Hennig, S Müller, and M Kästner.
Bézier extraction and adaptive refinement of truncated
hierarchical nurbs.
Computer Methods in Applied Mechanics and Engineering,
305:316–339, 2016.

[54]
René De Borst, Mike A Crisfield, Joris JC Remmers, and Clemens V Verhoosel.
Nonlinear finite element analysis of solids and structures.
John Wiley & Sons, 2012.

[55]
Clemens V. Verhoosel, Joris J. C. Remmers, and Miguel A. Gutiérrez.
A dissipation-based arc-length method for robust simulation of
brittle and ductile failure.
International Journal for Numerical Methods in Engineering,
77(9):1290–1321, 2009.

[56]
Mark S. Shephard.
Linear multipoint constraints applied via transformation as part of a
direct stiffness assembly process.
International Journal for Numerical Methods in Engineering,
20(11):2107–2112, 1984.

[57]
Erlend Storvik, Jakub Wiktor Both, Juan Michael Sargado, Jan Martin Nordbotten,
and Florin Adrian Radu.
An accelerated staggered scheme for variational phase-field models of
brittle fracture.
Computer Methods in Applied Mechanics and Engineering,
381:113822, 2021.

[58]
Kim Pham, Hanen Amor, Jean-Jacques Marigo, and Corrado Maurini.
Gradient damage models and their use to approximate brittle fracture.
International Journal of Damage Mechanics, 20(4):618–652,

[59]
Stefan May, Julien Vignollet, and René de Borst.
A new arc-length control method based on the rates of the internal
and the dissipated energy.
Engineering Computations, 33(1):100–115, 1 2016.

§ SOLUTION ALGORITHMS

$method \gets `displacement$'
$step \gets 0$
$\lambda \gets 0$
$\Delta\lambda \gets `stepsize$'

$step += 1$
$\lambda += \Delta\lambda$
not Converged
Solve Discrete Problem <ref>
Solve Discrete Problem <ref>

Compute energy dissipated ($\Delta G$) in the step using (<ref>)
$\Delta G \geq switchEnergy \:\&\&\: method == `displacement'$

$\Delta\lambda = 0$
$\Delta\tau = switchEnergy$

$iter < optiter$

$\Delta\tau *= \Delta\tau$
$\Delta\tau =  min(\Delta\tau,\Delta\tau_{max})$

$\Delta\tau = 0.5^{(0.25*(iter-optiter))}$
$\Delta\tau =  min(\Delta\tau,\Delta\tau_{max})$

Solution scheme without under-relaxation

$method \gets `displacement$'
$step \gets 0$
$\lambda \gets 0$
$\Delta\lambda \gets `stepsize$'

$step += 1$
$\lambda += \Delta\lambda$
not Converged

Solve Discrete Problem <ref>
Update solution fields with relaxation (<ref>)

Solve Discrete Problem <ref>
Update solution fields with relaxation (<ref>)

failed to Converge

Revert to old step solution and $\lambda$
method == `displacement'

$\Delta\lambda /= 10$

$fail += 1$
$fail > 2$

$\Delta\tau \,/= 2$
$\beta \,/=  1.25$

$\beta \,/=  1.25$

$fail \gets 0$
$\beta \,*= 1.05$
Compute energy dissipated ($\Delta G$) in the step using (<ref>)
$\Delta G \geq switchEnergy \:\&\&\: method == `displacement'$

$\Delta\lambda = 0$
$\Delta\tau = switchEnergy$

$iter < optiter$

$\Delta\tau *= \Delta\tau$
$\Delta\tau =  min(\Delta\tau,\Delta\tau_{max})$

$\Delta\tau = 0.5^{(0.25*(iter-optiter))}$
$\Delta\tau =  min(\Delta\tau,\Delta\tau_{max})$

Solution scheme with under-relaxation

§ SENSITIVITY W.R.T MAXIMUM DISSIPATION ALLOWED IN A STEP AND REFINEMENT INDICATOR

In this section, the sensitivity of the load-displacement curves w.r.t the maximum dissipated energy $\Delta\tau_{max}$ [N] and refinement indicator ($\pf_{th.}$) is investigated. To this end, the single edge notched specimen under tension (SENT) from Section <ref> is considered. All aspects of the model (geometry, material parameters and loading conditions) remain same, the only variation being in $\Delta \tau_{max}$ and the degradation function set to quadratic.

Figure <ref> presents the load-displacement curves obtained using the proposed solver for different values of $\Delta\tau_{max}$. For the values chosen in this study, the load-displacement curves are similar. This indicates that the relevant features (snap-back behaviour) could be sufficiently captured even with a larger dissipation steps, for instance $\Delta\tau_{max} = 0.05$ [N].

The sensitivity of the load-displacement curves obtained using the proposed solver with different refinement indicator is studied with three different threshold phase-field values. An element is marked for refinement once this threshold value is exceeded. Figure <ref> presents the load-displacement curves obtained for the different phase-field threshold values along with that obtained using a fixed mesh. It is observed that the curves are similar. This indicates that the proposed monolithic solver does not exhibit any bias w.r.t. adaptive mesh refinement techniques.

[thick,scale=0.95, every node/.style=scale=0.9]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.0065, ymax=1000, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/VaryDis/lodi_dTau0_00625.txt;
[col sep = comma]./Data/VaryDis/lodi_dTau0_0125.txt;
[col sep = comma]./Data/SENT/lodi_lby2pf0_2.txt;
[col sep = comma]./Data/VaryDis/lodi_dTau0_05.txt;

x expr=1,
y expr=0
] ;
$\Delta\tau_{max} = 0.00625$

x expr=1,
y expr=0
] ;
$\Delta\tau_{max} = 0.0125$

x expr=1,
y expr=0
] ;
$\Delta\tau_{max} = 0.025$

x expr=1,
y expr=0
] ;
$\Delta\tau_{max} = 0.05$

[thick,scale=0.95, every node/.style=scale=0.9]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.0065, ymax=1000, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/SENT/lodi_lby2pf0_1.txt;
[col sep = comma]./Data/SENT/lodi_lby2pf0_2.txt;
[col sep = comma]./Data/SENT/lodi_lby2pf0_4.txt;
[col sep = comma]./Data/SENT/lodi_FixMesh.txt;

mark size=0.25pt,
x expr=1,
y expr=0
] ;
$\pf_{th.} = 0.1$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
$\pf_{th.} = 0.2$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
$\pf_{th.} = 0.4$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Fixed Mesh

Figure (a) presents the load-displacement curves for the single edge notched specimen under tension. The legend entries correspond to the choice of prescribed maximum dissipation energy $\Delta\tau_{max}$. Figure (b) presents the load-displacement curves for the single edge notched specimen under tension. The legend entries correspond to the choice of solution techniques.

§ COMPARISON WITH ALTERNATIVE MINIMIZATION SOLVER AND QUASI-NEWTON RAPHSON METHOD [26]

In this section, a comparison of the proposed monolithic solver with alternative minimization solver [3, 33, 6] and the quasi Newton-Raphson method proposed by [26] is carried out. The comparison is based only on the load-displacement curves obtained from the respective solvers. The alternate minimization solver solves the displacement and the phase-field sub-problems alternatively until a certain tolerance criterion is met. The quasi Newton Raphson approach proposed in [26] adopts a linear extrapolation of the phase-field variable for the momentum balance equation. Although, the extrapolation strategy is questionable, it yields a convex energy functional, resulting in a better convergence behaviour of the Newton Raphson method.

For the comparative study, the SENT model from Section <ref> is chosen. For both the alternate minimization solver and the quasi Newton Raphson method, the fracture irreversibility is enforced using the penalization method [25, 34]. Figure <ref> presents the load-displacement curves obtained using the solver proposed in this work, and the alternate minimization solver and the quasi-Newton Raphson method. It is observed that all methods/solvers predict a similar peak load and pre-peak behaviour. However, the post-peak behaviours are different. The proposed monolithic solver equipped with the arc-length method is able to predict snap-back behaviour, which is not possible with conventional alternate minimization solver and the quasi-Newton Raphson method.

[thick,scale=1.15, every node/.style=scale=0.9]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.0065, ymax=1000, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/VaryDis/lodi_dTau0_00625.txt;
[col sep = comma]./Data/CompAMQNR/lodiAM.txt;
[col sep = comma]./Data/CompAMQNR/lodiQNR.txt;

x expr=1,
y expr=0
] ;
Present solver

x expr=1,
y expr=0
] ;
Alt. min.

x expr=1,
y expr=0
] ;
Quasi NR

Figure presents the load-displacement curves for the single edge notched specimen under tension. The legend entries correspond to the choice of solution techniques.

[thick,scale=0.95, every node/.style=scale=0.95]
[legend style=draw=none, legend columns = 2,
transpose legend, ylabel=Load [N],xlabel=Displacement [mm], xmin=0, ymin=0, xmax=0.014, ymax=800, yticklabel style=/pgf/number format/.cd,fixed,precision=2,
every axis plot/.append style=very thick]
[col sep = comma]./Data/SENS/lodi_quad.txt;

[col sep = comma]./Data/SENS/lodi_cubics0_1.txt;
[col sep = comma]./Data/SENS/lodi_cubics1.txt;

mark size=0.25pt,
x expr=1,
y expr=0
] ;

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 0.1$

mark size=0.25pt,
x expr=1,
y expr=0
] ;
Cubic: $s = 1.0$

[inner sep=0pt] () at (0,0)
[inner sep=0pt] () at (-1.3,-3.05)
hide axis,
scale only axis,

colorbar horizontal,
point meta min=0,
point meta max=1,
colorbar style=
xticklabel style = yshift=-0.075cm
[draw=none] coordinates (0,0);
[inner sep=0pt] () at (0,-3.75) $\pf$;

Figure (a) presents the load-displacement curves for the single edge notched specimen under shear. The legend entries correspond to the choice of degradation functions. Figure (b) shows the distribution of the phase-field variable at the final step of the analysis.