An iterative staggered scheme for phase field brittle fracture propagation with stabilizing parameters

TL;DR

This paper introduces a new iterative staggered scheme with stabilization for simulating brittle fracture propagation using phase field models, improving robustness and efficiency, and validated through numerical benchmarks.

Contribution

A novel iterative scheme with stabilization terms for phase field brittle fracture models, including convergence analysis and numerical validation.

Findings

Scheme ensures bounded elastic energy under natural conditions.

Monotonic convergence achieved with sufficiently thick diffusive zones.

Validated effectiveness through benchmark problems.

Abstract

This paper concerns the analysis and implementation of a novel iterative staggered scheme for quasi-static brittle fracture propagation models, where the fracture evolution is tracked by a phase field variable. The model we consider is a two-field variational inequality system, with the phase field function and the elastic displacements of the solid material as independent variables. Using a penalization strategy, this variational inequality system is transformed into a variational equality system, which is the formulation we take as the starting point for our algorithmic developments. The proposed scheme involves a partitioning of this model into two subproblems; phase field and mechanics, with added stabilization terms to both subproblems for improved efficiency and robustness. We analyze the convergence of the proposed scheme using a fixed point argument, and find that under a natural condition, the elastic mechanical energy remains bounded, and, if the diffusive zone around crack surfaces is sufficiently thick, monotonic convergence is achieved. Finally, the proposed scheme is validated numerically with several bench-mark problems.

Figures

Equations

$$\left(g(\phi )\text{amsmathbb}C\mathbf{e}(u),\mathbf{e}(v)\right)$$

$$G_c\epsilon (\nabla \phi ,\nabla \psi )-\frac{G_c}{\epsilon }(1-\phi ,\psi )+(1-\kappa )(\phi |\text{amsmathbb}C\mathbf{e}(u){|}^2,\psi )$$

$${\lambda }_m\Vert \mathbf{u}\Vert \le (\text{amsmathbb}C\mathbf{u},\mathbf{u}{)}^{\frac{1}{2}}\le {\lambda }_M\Vert \mathbf{u}\Vert ,\text{ for }\mathbf{u},\mathbf{v}\in (L^2(B){)}^{d\times d},\mathbf{u},\mathbf{v}\ne 0.$$

$$u^n$$

$${\phi }^n$$

$$\left(g({\phi }^n)\text{amsmathbb}C\mathbf{e}(u^n),\mathbf{e}(v)\right)=(b^n,v),$$

$$G_c\epsilon (\nabla {\phi }^n,\nabla \psi )-\frac{G_c}{\epsilon }(1-{\phi }^n,\psi )+(1-\kappa )({\phi }^n|\text{amsmathbb}C\mathbf{e}(u^n){|}^2,\psi )$$

$$+([\Xi +\gamma ({\phi }^n-{\phi }^{n-1}){]}^{+},\psi )=0,$$

$$\Vert g(\psi )-g(\eta )\Vert \le 2(1-\kappa )\Vert \psi -\eta \Vert ,\forall \psi ,\eta \in W.$$

$$(u^{n,i-1},{\phi }^{n,i-1},b^n)$$

$$a_u(u^{n,i},v):=L_u(u^{n,i}-u^{n,i-1},v)+\left(g({\phi }^{n,i-1})\text{amsmathbb}C\mathbf{e}(u^{n,i}),\mathbf{e}(v)\right)=(b^n,v),$$

$$({\phi }^{n,i-1},u^{n,i},{\phi }^{n-1})$$

$$a_{\phi }({\phi }^{n,i},\psi ):=L_{\phi }({\phi }^{n,i}-{\phi }^{n,i-1},\psi )+G_c\epsilon (\nabla {\phi }^{n,i},\nabla \psi )-\frac{G_c}{\epsilon }(1-{\phi }^{n,i},\psi )$$

$$+(1-\kappa )({\phi }^{n,i}|\text{amsmathbb}C\mathbf{e}(u^{n,i}){|}^2,\psi )+({\eta }^i(\Xi +\gamma ({\phi }^{n,i}-{\phi }^{n-1})),\psi )=0,$$

$${\eta }^i(x)=\{\begin{array}{l}1,\text{ if }\Xi (x)+\gamma ({\phi }^{n,i}(x)-{\phi }^{n-1}(x))\ge 0,\\ 0,\text{ if }\Xi (x)+\gamma ({\phi }^{n,i}(x)-{\phi }^{n-1}(x))<0.\end{array}$$

$$e_u^i$$

$$e_{\phi }^i$$

$$L_u(e_u^i-e_u^{i-1},v)+(g({\phi }^n)\text{amsmathbb}C\mathbf{e}(e_u^i),\mathbf{e}(v))+((g({\phi }^{n,i-1})-g({\phi }^n))\text{amsmathbb}C\mathbf{e}(u^{n,i}),\mathbf{e}(v))=0,$$

$$L_{\phi }(e_{\phi }^i-e_{\phi }^{i-1},\psi )+G_c\epsilon (\nabla e_{\phi }^i,\nabla \psi )+\frac{G_c}{\epsilon }(e_{\phi }^i,\psi )+\gamma ({\eta }^i{e}_{\phi }^i,\psi )+(1-\kappa )(e_{\phi }^i|\text{amsmathbb}C\mathbf{e}(u^{n,i}){|}^2,\psi )$$

$$+(1-\kappa )\left({\phi }^n\text{amsmathbb}C\mathbf{e}(e_u^i):\mathbf{e}(u^{n,i}+u^n),\psi \right)=0,$$

$$\mathrm{ess}\underset{x\in B}{\sup }|\mathbf{e}(u^n(x))|\le M,\forall n.$$

$$\mathrm{ess}\underset{x\in B}{\sup }|\mathbf{e}(u^{n,i}(x))|\le M,\forall (n,i).$$

$$\left(\frac{L_{\phi }}{2}+\frac{G_c}{\epsilon }+\frac{G_c\epsilon }{c_P}-8\xi \frac{(1-\kappa {)}^2}{\kappa }\right)\Vert e_{\phi }^i{\Vert }^2+\left(\frac{L_u}{2}+\frac{\kappa {\lambda }_{\text{min}}^2}{2c_P{c}_K}\right)\Vert e_u^i{\Vert }^2$$

$$\le \left(\frac{L_{\phi }}{2}+8\xi \frac{(1-\kappa {)}^2}{\kappa }\right)\Vert e_{\phi }^{i-1}{\Vert }^2+\frac{L_u}{2}\Vert e_u^{i-1}{\Vert }^2,$$

$${\epsilon }^2-16\xi \frac{(1-\kappa {)}^2}{\kappa }\frac{c_P}{G_c}\epsilon +c_P>0,$$

$$\left(\frac{L_{\phi }}{2}+\frac{G_c}{\epsilon }\right)\Vert e_{\phi }^i{\Vert }^2+\frac{L_{\phi }}{2}\Vert e_{\phi }^i-e_{\phi }^{i-1}{\Vert }^2+G_c\epsilon \Vert \nabla e_{\phi }^i{\Vert }^2+\gamma ({\eta }^i{e}_{\phi }^i,e_{\phi }^i)$$

$$+(1-\kappa )(e_{\phi }^i|\text{amsmathbb}C\mathbf{e}(u^{n,i}){|}^2,e_{\phi }^i)+\frac{L_u}{2}\Vert e_u^i{\Vert }^2+\frac{L_u}{2}\Vert e_u^i-e_u^{i-1}{\Vert }^2+(g({\phi }^n)\text{amsmathbb}C\mathbf{e}(e_u^i),\mathbf{e}(e_u^i))$$

$$=\frac{L_{\phi }}{2}\Vert e_{\phi }^{i-1}{\Vert }^2+\frac{L_u}{2}\Vert e_u^{i-1}{\Vert }^2-(1-\kappa )({\phi }^n\text{amsmathbb}C\mathbf{e}(e_u^i):\mathbf{e}(u^{n,i}+u^n),e_{\phi }^i)$$

$$-((g({\phi }^{n,i-1})-g({\phi }^n))\text{amsmathbb}C\mathbf{e}(u^{n,i}),\mathbf{e}(e_u^i)),$$

$$2(x-y,x)=\Vert x{\Vert }^2+\Vert x-y{\Vert }^2-\Vert y{\Vert }^2.$$

$$\left(\frac{L_{\phi }}{2}+\frac{G_c}{\epsilon }\right)\Vert e_{\phi }^i{\Vert }^2+G_c\epsilon \Vert \nabla e_{\phi }^i{\Vert }^2+\frac{L_u}{2}\Vert e_u^i{\Vert }^2+\kappa (\text{amsmathbb}C\mathbf{e}(e_u^i),\mathbf{e}(e_u^i))$$

$$\le \frac{L_{\phi }}{2}\Vert e_{\phi }^{i-1}{\Vert }^2+\frac{L_u}{2}\Vert e_u^{i-1}{\Vert }^2+(1-\kappa ){\int }_B|\text{amsmathbb}C\mathbf{e}(e_u^i):\mathbf{e}(u^{n,i}+u^n)e_{\phi }^i|\text{d}x$$

$$+{\int }_B|(g({\phi }^{n,i-1})-g({\phi }^n))\text{amsmathbb}C\mathbf{e}(u^{n,i}):\mathbf{e}(e_u^i)|\text{d}x$$

$$\displaystyle\qquad\quad\leq\frac{L_{\varphi}}{2}\lVert e_{\varphi}^{i-1}\rVert^{2}+\frac{L_{u}}{2}\lVert e_{u}^{i-1}\rVert^{2}+2(1-\kappa)\lambda_{\textnormal{max}}M\Big{(}\lVert e_{\varphi}^{i}\rVert+\lVert e_{\varphi}^{i-1}\rVert\Big{)}\lVert\mathbf{e}(e_{u}^{i})\rVert,$$

$$\left(\frac{L_{\phi }}{2}+\frac{G_c}{\epsilon }-2(1-\kappa ){\lambda }_{\text{max}}M\frac{1}{2{\delta }_1}\right)\Vert e_{\phi }^i{\Vert }^2+G_c\epsilon \Vert \nabla e_{\phi }^i{\Vert }^2$$

$$\displaystyle\quad+\frac{L_{u}}{2}\lVert e_{u}^{i}\rVert^{2}+\Big{(}\kappa\lambda_{\textnormal{min}}^{2}-2(1-\kappa)\lambda_{\textnormal{max}}M(\delta_{1}+\delta_{2})\Big{)}\lVert\mathbf{e}(e_{u}^{i})\rVert^{2}$$

$$\le \left(\frac{L_{\phi }}{2}+2(1-\kappa ){\lambda }_{\text{max}}M\frac{1}{2{\delta }_2}\right)\Vert e_{\phi }^{i-1}{\Vert }^2+\frac{L_u}{2}\Vert e_u^{i-1}{\Vert }^2,$$

$$\left(\frac{L_{\phi }}{2}+\frac{G_c}{\epsilon }-8\xi \frac{(1-\kappa {)}^2}{\kappa }\right)\Vert e_{\phi }^i{\Vert }^2+G_c\epsilon \Vert \nabla e_{\phi }^i{\Vert }^2+\frac{L_u}{2}\Vert e_u^i{\Vert }^2+\frac{\kappa {\lambda }_{\text{min}}^2}{2}\Vert \mathbf{e}(e_u^i){\Vert }^2$$

$$\le \left(\frac{L_{\phi }}{2}+8\xi \frac{(1-\kappa {)}^2}{\kappa }\right)\Vert e_{\phi }^{i-1}{\Vert }^2+\frac{L_u}{2}\Vert e_u^{i-1}{\Vert }^2.$$

$$\Vert e_{\phi }^i{\Vert }^2\le c_P\Vert \nabla e_{\phi }^i{\Vert }^2\text{ and }\Vert e_u^i{\Vert }^2\le c_P{c}_K\Vert \mathbf{e}(e_u^i){\Vert }^2,$$

$$\left(\frac{L_{\phi }}{2}+\frac{G_c}{\epsilon }+\frac{G_c\epsilon }{c_P}-8\xi \frac{(1-\kappa {)}^2}{\kappa }\right)\Vert e_{\phi }^i{\Vert }^2+\left(\frac{L_u}{2}+\frac{\kappa {\lambda }_{\text{min}}^2}{2c_P{c}_K}\right)\Vert e_u^i{\Vert }^2$$

Full text

Key words: phase field; quasi-static; brittle fracture; fracture propagation; $L$-scheme; fixed stress; iterative algorithm; linearization; convergence analysis; fixed point; finite element;

1 Introduction

Fracture propagation is currently an important topic with many applications in various engineering fields. Specifically, phase-field descriptions are intensively investigated. The theory of brittle fracture mechanics goes back to the works of A. Griffith [23], wherein a criterion for crack propagation is formulated. Despite a foundational treatment on the subject of brittle fracture, Griffith’s theory fails to predict crack initiation. This deficiency can however be overcome by a variational approach, which was first proposed in [10, 20]. Using such a variational approach, discontinuities in the displacement field $u$ across the lower-dimensional crack surface are approximated by an auxiliary phase-field function $\varphi$. The latter can be viewed as an indicator function, which introduces a diffusive transition zone between the broken and the unbroken material. The enforcement of irreversibility of crack growth finally yields a variational inequality system, of which we seek the solution $\{u,\varphi\}$.

In this work, we concentrate on improvements of the nonlinear solution algorithm, which is still a large bottleneck of phase-field fracture evolution problems. Specifically, high iteration numbers when the crack initiates or is further growing are reported in many works [21, 29, 44, 45]. However, in most studies iteration numbers are omitted. Both staggered (splitting) schemes and monolithic schemes are frequently employed. Important developments include alternating minimization/staggered schemes [9, 11, 12, 29, 30], quasi-monolithic scheme with a partial linearization [25], and fully monolithic schemes [21, 44, 45].

The goal of this work is to propose a linearized staggered scheme with stabilizing parameters. In particular, the proposed scheme is based on recent developments on iterative splitting schemes coming from poroelasticity [13, 26, 32, 33]. Iterative splitting schemes are widely applied to problems of coupled flow and mechanics, where at each iteration step either of the subproblems (i.e., flow or mechanics) is solved first, keeping some physical quantity constant (e.g., fixed stress or fixed strain), followed by solving the next subproblem with updated solution information. This procedure is then repeated until an accepted tolerance is reached. Further extensions of this technique involves tuning some artificial stabilization terms according to a derived contraction estimate in energy norms. Here, the quantity held constant during solving of the subproblems need not represent any physical quantity present in the model. This is the central idea in the so-called ‘$L$-scheme’, which has proven to perform robustly for Richards equation [28, 38], for linear and nonlinear poroelasticity [7, 8], and for nonlinear thermo-poroelasticity [27].

We propose here a variant of the $L$-scheme, adapted to phase field brittle fracture propagation models. This scheme is based on a partitioning of the model into two subproblems; phase field and mechanics. Here, the $L$-scheme acts both as a stabilization and as a linearization (as a linearization scheme, the stabilization parameters mimics the Jacobian from Newton iteration). Assuming that the mechanical elastic energy remains bounded during the iterations, and that the diffusive zone around crack surfaces is sufficiently thick, we give a proof of monotonic convergence of the proposed scheme by employing a fixed point argument.

The efficiency and robustness of the proposed scheme is demonstrated numerically with several bench-mark problems. Moreover, we compare the number of iterations needed for convergence with ‘standard’ staggered schemes (i.e., without stabilizing terms), and monolithic schemes in which the fully-coupled system is solved all-at-once. Furthermore, it is well known that when reaching the critical loading steps during the computation of brittle fracture phase field problems (i.e., when the crack is propagating), spikes in iteration numbers appear. For this reason, and thanks to the monotonic convergence property of the proposed scheme, we show that a (low) upper bound on the number of iterations may be enforced, while the computed results are still in very good agreement with the non-truncated solutions. Thus, using this ‘truncated $L$-scheme’, we effectively avoid the iteration spikes at the critical loading steps at the cost of negligible loss of accuracy. We mention that this strategy is not available with e.g. Newton iteration, as the iterate solutions may behave erratically for any number of iterations before finally converging. Moreover, the assumption that the mechanical elastic energy remains bounded during the iterations is verified numerically for all tests cases.

The main aims of this work are three-fold: Under a natural assumption, we prove the convergence of a novel iterative staggered scheme, optimized for phase field brittle fracture propagation problems. Based on these theoretical findings, we design a robust solution algorithm with monotonic convergence properties. Finally, several numerical tests are presented in which our variants of the $L$-scheme are tested in detail.

The outline of this paper is as follows: In Section 2 we present the model equations and coefficients, in Section 3 we introduce the partitioned scheme and derive a convergence proof, in Section 4 we describe in detail our numerical algorithm in pseudo-code, and in Section 5 we provide several numerical experiments, in particular the single edge notched tension test, the single edge notched shear test, and the L-shaped panel test. Finally, in Section 6 we provide some conclusions and summary of the work.

1.1 Preliminaries

In this section we explain the notation used throughout this article, see e.g. [18, 47] for more details. Given an open and bounded set $B\subset{}^{d}$, $d\in\{2,3\}$, and $1\leq p<\infty$, let $L^{p}(B)=\{f:B\rightarrow\real:\int_{B}|f(x)|^{p}\textnormal{d}x<\infty\}$. For $p=\infty$, let $L^{\infty}(B)=\{f:B\rightarrow\real:\operatorname{ess}\sup_{x\in B}|f(x)|<\infty\}$. In particular, $L^{2}(B)$ is the Hilbert space of square integrable functions with inner product $(\cdot,\cdot)$ and norm $\lVert f\rVert:=(f,f)^{\frac{1}{2}}$. For $k\in{\amsmathbb N}$, $k\geq 0$, we denote by $W^{k,p}(B)$ the space of functions in $L^{p}(B)$ admitting weak derivatives up to $k$’th order. In particular, $H^{1}(B):=W^{1,2}(B)$ and we denote by $H_{0}^{1}(B)$ its zero trace subspace.

Note that we reserve the use of bold fonts for second order tensors. Hence, if $u,v\in L^{2}(B)$, their inner product is $(u,v):=\int_{B}u(x)v(x)\textnormal{d}x$, and similarly, if $u,v\in(L^{2}(B))^{d}$ then we take their inner product to be $(u,v):=\int_{B}u(x)\cdot v(x)\textnormal{d}x$. Finally, if $\mathbf{u},\mathbf{v}\in(L^{2}(B))^{d\times d}$ then their inner product is $(\mathbf{u},\mathbf{v}):=\int_{B}\mathbf{u}(x):\mathbf{v}(x)\textnormal{d}x$.

We will also frequently apply several classical inequalities, in particular: Cauchy-Schwarz, Young, Poincaré, and Korn. See e.g. [15, 24] for a detailed description of these.

2 Governing equations

What follows is a brief description of the phase field approach for quasi-static brittle fracture propagation, see e.g. [20, 30] for more details. Consider a (bounded open) polygonal domain $B\subset{}^{d}$, wherein $\mathcal{C}\subset{}^{d-1}$ denotes the fracture, and $\Omega\subset{}^{d}$ is the intact domain, and a time interval $(0,T)$ is given with final time $T>0$. By introducing the phase field variable $\varphi:B\times(0,T)\rightarrow[0,1]$, which takes the value [math] in the fracture, $1$ in the intact domain, and varies smoothly from [math] to $1$ in a transition zone of (half-)thickness $\varepsilon>0$ around $\mathcal{C}$, the evolution of the fracture can be tracked in space and time. Using the phase field approach, the fracture $\mathcal{C}$ is approximated by $\Omega_{F}\subset{}^{d}$, where $\Omega_{F}:=\{x\in{}^{d}:\varphi(x)<1\}$.

Introducing the displacement vector $u:B\times(0,T)\rightarrow{}^{d}$, the model problem we consider arises as a minimization problem: An energy functional $E(u,\varphi)$ is defined according to Griffith’s criterion for brittle fracture [23], which is then sought to be minimized over all admissible $\{u,\varphi\}$. From this minimization problem, the Euler-Lagrange equations are obtained by differentiation with respect to the arguments, yielding a variational equality system. Finally, a crack irreversibility condition must be enforced (the crack is not allowed to heal), which takes the form $\partial_{t}\varphi\leq 0$. Thus, the variational equality system, which is the previously mentioned Euler-Lagrange equations, is transformed into a variational inequality system, which reads as follows: Find $(u(t),\varphi(t))\in V\times W:=(H_{0}^{1}(B))^{d}\times W^{1,\infty}(B)$ such that for $t\in(0,T]$ there holds

[TABLE]

where $G_{c}>0$ is the critical elastic energy restitution rate, $0<\kappa<<1$ is a regularization parameter, the purpose of which is to avoid degeneracy of the elastic energy (equivalent with replacing the fracture with a softer material), and $g(\varphi):=(1-\kappa)\varphi^{2}+\kappa$ is a standard choice for the degradation function (see e.g. [39, 45]. Note that $g(\varphi)\rightarrow\kappa$ when approaching the fracture zone). The body force acting on the domain $B$ is $b:B\times(0,T)\rightarrow{}^{d}$, and $|\amsmathbb{C}\mathbf{e}(u)|^{2}:=\amsmathbb{C}\mathbf{e}(u):\mathbf{e}(u)$ is the elastic mechanical energy, where $\mathbf{e}(\cdot):=(\nabla(\cdot)+\nabla(\cdot)^{\top})/2$ is the symmetric gradient, and $\amsmathbb{C}=[C_{ijkl}]_{ijkl}$ is the fourth order tensor containing the elastic material coefficients, where each $C_{ijkl}\in L^{\infty}(B)$. We assume that $\amsmathbb{C}$ satisfies the usual symmetry and positive definiteness properties, i.e., $(\amsmathbb{C}\mathbf{u},\mathbf{v})=(\mathbf{u},\amsmathbb{C}\mathbf{v})$, and $(\amsmathbb{C}\mathbf{u},\mathbf{u})^{\frac{1}{2}}$ defines an $L^{2}$-equivalent norm, i.e., there exists constants $\lambda_{m},\lambda_{M}>0$ such that

[TABLE]

In order to facilitate the following developments we assume continuity in time for $\{u,\varphi,b\}$. Let now $0=t^{0}<t^{1}<\cdots<t^{N}=T$ be a partition of the time interval $(0,T)$, with time step $\delta t:=t^{n}-t^{n-1}$, and denote the time discrete solutions by

[TABLE]

The irreversibility condition now becomes $\varphi^{n}\leq\varphi^{n-1}$ (using a backward Euler method), and the time-discrete version of the problem (2.1a)–(2.1b) reads as follows: Find $(u^{n},\varphi^{n})\in V\times W$ such that

[TABLE]

where $b^{n}:=b(\cdot,t^{n})$. The last term in the phase field equation (2.5b) is a penalization to enforce the irreversibility condition, thus transforming the variational inequality (2.1b) into a variational equality, with penalization parameter $\gamma>0$, and where $\Xi\in L^{2}(B)$ is given (in practice $\Xi$ will be obtained by iteration, cf. Section 4). Note that we also used the notation $[x]^{+}:=\max(x,0)$. From here on, we shall refer to (2.5a) as the mechanics subproblem, and to (2.5b) as the phase field subproblem. Regarding the degradation function $g$, it is easily seen to satisfy the following Lipschitz condition:

[TABLE]

The time-discrete system (2.5a)-(2.5b) was analyzed in [36], and there it was shown that at least one global minimizer $(u^{n},\varphi^{n})\in V\times W$ exists, provided $b^{n}\in(L^{2}(B))^{d}$, for each $n$. We mention also that the analysis of a pressurized phase field brittle fracture model can be found in [34, 35].

3 Iterative scheme

In this section we introduce the iterative staggered solution procedure for the fully discrete formulation of (2.5a)-(2.5b). To this end, let $\mathcal{T}_{h}$ be a simplicial mesh of $B$, such that for any two distinct elements of $\mathcal{T}_{h}$ their intersection is either an empty set or their common vertex or edge. We denote by $h$ the largest diameter of all the elements in $\mathcal{T}_{h}$, i.e., $h:=\max_{K\in\mathcal{T}_{h}}\operatorname{diam}(K)$, and let $V_{h}\times W_{h}\subset V\times W$ be appropriate (conforming) discrete spaces. We continue now with the same notation for the variables and test functions as before (omitting the usual $h$-subscript), since we will from here on mostly deal with the discrete solutions.

For each $n$, the iterative algorithm we propose defines a sequence $\{u^{n,i},\varphi^{n,i}\}$, for $i\geq 0$, initialized by $\{u^{n-1},\varphi^{n-1}\}$. The iteration is then done in two steps: First, the mechanics subproblem is solved, with the degradation function held constant. Then, the phase field subproblem is solved, with the elastic energy held constant. Note that there are also artificial stabilizing terms which are held constant during solving of the subproblems. Introducing the stabilization parameters $L_{u},L_{\varphi}>0$ (to be determined later), the iterative algorithm reads as follows:

[TABLE]

where, in order to avoid the $[\cdot]^{+}$-bracket, we also introduced the function $\eta^{i}\in L^{\infty}(B)$ defined for a.e. $x\in B$ by

[TABLE]

3.1 Convergence analysis

We now proceed to analyze the convergence of the scheme (3.1a)-(3.1b). Our aim is to show a contraction of successive difference functions in energy norms, which implies convergence by the Banach Fixed Point Theorem (see e.g. [14]). To this end we define the following difference functions

[TABLE]

where $\{u^{n},\varphi^{n}\}$ denotes the (exact) solutions to (2.1a)-(2.1b) at time $t^{n}$. Using the symmetry properties of $\amsmathbb{C}$, the following set of difference equations are then obtained by subtracting (3.1a)-(3.1b) solved by $\{u^{n},\varphi^{n}\}$ from the same equations solved by the iterate solutions:

[TABLE]

Furthermore, we introduce the following assumption related to the elastic mechanical strain.

Assumption 1** (Boundedness of elastic strain).**

We assume there exists a constant $M>0$ such that

[TABLE]

Moreover, we assume that $M$ is large enough such that the above bound holds also for the iterate elastic strain, i.e.,

[TABLE]

Note that $M$ is nothing else than an upper bound for the elastic strain in the system for the converged solution, which is arguably finite for any reasonable problem. Note also that with sufficient regularity of the domain, coefficients, source terms, and initial data, the above assumption is satisfied, i.e., the problem (2.5a)-(2.5b) admits a solution $u^{n}\in(W^{1,\infty}(B))^{d}$, thus implying the existence of $M$. Alternatively to introducing the constant $M$, we could introduce instead a so-called ‘cut-off operator’ in the iterate equations (3.1a)-(3.1b), as seen in e.g. [40, 41]. Note that in all numerical tests to be done in the next sections, we provide figures validating the second part of this assumption (cf. Section 5.4). With the above definitions, we state our main theoretical result.

Theorem 3.1** (Convergence of the scheme).**

The scheme (3.1b)–(3.1a) defines a contraction satisfying

[TABLE]

if $L_{u},L_{\varphi}>0$, and if the model parameter $\varepsilon>0$ is sufficiently large such that

[TABLE]

where $\xi:=(M\lambda_{\textnormal{max}}/\lambda_{\textnormal{min}})^{2}>0$, and where $c_{P},c_{K}>0$ are the Poincaré and Korn constants, respectively, depending only on the domain $B$ and spatial dimension $d$.

Proof.

We begin by taking $v=e_{u}^{i}$ and $\psi=e_{\varphi}^{i}$ in (3.5a) and (3.5b), respectively, add the resulting equations together and obtain

[TABLE]

where we used the following inner product identity

[TABLE]

Discarding some non-negative terms from the left hand side of (3.10), using the fact that $\operatorname{ess}\sup_{x\in B}\varphi^{n}(x)\leq 1$, in addition to the Lipschitz property of the degradation function $g$ (2.6), yields

[TABLE]

where we also invoked the Assumption 1 in the last line, and applied the Cauchy-Schwarz inequality. Using the Young inequality, the properties of elastic tensor (2.2), and rearranging, leads to

[TABLE]

for some constants $\delta_{1},\delta_{2}>0$. Choosing $\delta_{1}=\delta_{2}=\kappa\lambda_{\textnormal{min}}^{2}/8(1-\kappa)\lambda_{\textnormal{max}}M$ yields (3.13) as

[TABLE]

Next, by applying the Poincaré inequality on $\lVert e_{\varphi}^{i}\rVert$, and by applying successively the Poincaré and Korn inequalities on $\lVert e_{u}^{i}\rVert$, we obtain

[TABLE]

where $c_{P},c_{K}$ are the (squares of the) Poincaré and Korn constants, respectively (depending only on the domain $B$ and spatial dimension $d$). Finally, employing these bounds on the left hand side of (3.14) yields

[TABLE]

Thus, for (3.16) to be a contraction estimate, $\varepsilon$ must satisfy the following second order inequality

[TABLE]

Setting the left hand side of (3.17) equal to zero yields a second order polynomial, the discriminant of which must satisfy one of the following three statements:

If

[TABLE]

then $P(\varepsilon)=0$ has two distinct positive real roots $\varepsilon_{1},\varepsilon_{2}>0$, in which case (3.16) is a contraction for $\varepsilon\in(0,\varepsilon_{1})\cup(\varepsilon_{2},\infty)$. 2. 2.

If

[TABLE]

then $P(\varepsilon)=0$ has one positive real root, $\varepsilon_{0}>0$, of multiplicity two, in which case (3.16) is a contraction for all $\varepsilon\neq\varepsilon_{0},\varepsilon>0$. 3. 3.

If

[TABLE]

then $P(\varepsilon)=0$ has two complex roots, in which case (3.16) is a contraction for all $\varepsilon>0$.

∎

Remark 3.1** (Convergence rate).**

According to the above proof, if the scheme is not converging for a given value of $\varepsilon$, then a larger or a smaller value may be chosen to rectify the situation. However, since crack surfaces become singular as $\varepsilon\rightarrow 0$ (thus necessitating finer meshing, i.e., $h\rightarrow 0$), we choose to state Theorem 3.1 with the condition that $\varepsilon$ be large enough. We note also that due to some unknown constants in the convergence rate (3.1), it is not known whether this rate is optimal. Furthermore, working with a large $\varepsilon$ is substantiated by the theory of phase field fracture being based on $\Gamma$ convergence [2, 3]. Applying this to phase field fracture was first done in [10]. Specifically, the setting is suitable when $h=o(\varepsilon)$; namely when $\varepsilon$ is sufficiently large.

4 Algorithm

In practice, we apply the stabilizations and penalizations proposed in the previous sections as outlined below. It is well-known (e.g., [37]) that the choice of $\gamma$ is critical. If $\gamma$ is too low, crack irreversibility will not be enforced. On the other hand, if $\gamma$ is too large, the linear equation system is ill-conditioned and influences the performance of the nonlinear solver. For this reason, $\gamma$ is updated in at each iteration step. Better, in terms of robustness, is the augmention in such an iteration by an additional $L^{2}$ function $\Xi$, yielding a so-called augmented Lagrangian iteration going back to [19, 22]. For phase-field fracture this idea was first applied in [42]. Thus, combining the staggered iteration for the solid and phase-field systems with the update of the penalization parameter $\Xi$ yields the following algorithm:

Algorithm 1**.**

*At the loading step $t^{n}$. *

*Choose initial $\Xi^{0}$. Set $\gamma>0$. *

repeat

*     Iterate on $i$ (augmented Lagrangian loop) *

*     Solve two-field problem, namely *

*          Solve elasticity in Problem (3.1a) *

*          Solve the nonlinear phase-field in Problem (3.1b) *

*     Update *

[TABLE]

until* *

[TABLE]

[TABLE]

*Set: $(u^{n},\varphi^{n}):=(u^{n,i},\varphi^{n,i})$. *

*Increment $t^{n}\rightarrow t^{n+1}$. *

For the stabilization parameters $L_{u},L_{\varphi}$, we have the following requirements (somewhat similar to $\gamma$): If the stabilization is too small, the stabilization effects vanish. If the stabilization is too large, we revert to an unacceptably slow convergence, and potentially, may converge to a solution corresponding to an undesirable local minimum of the original problem. In order to deal with these issues, we employ here a simple, yet effective strategy: We draw $L:=L_{u}=L_{\varphi}$ from a range of suitable values and compare the results, i.e., $L\in\{1.0e-6,1.0e-3,1.0e-2,1.0e-1\}$. Moreover, we include also the configurations $L_{u}=0$, $L_{\varphi}>0$ and $L_{u}=L_{\varphi}=0$ in all the numerical tests to be done in the following.

Remark 4.1**.**

In this paper we use $\operatorname{TOL}=10^{-6}$.

4.1 Nonlinear solution, linear subsolvers and programming code

Both subproblems (phase field and mechanics) may be nonlinear. In our theory presented above, we assumed a standard elasticity tensor. However, the model (3.1a)–(3.1b) is too simple for most mechanical applications. More realistic phase-field fracture applications require a splitting of the stress tensor (based on an energy split) in order to account for fracture development only under tension, but not under compressive forces. Consequently, we follow here [31] and split $\boldsymbol{\sigma}$ into tensile $\boldsymbol{\sigma}^{+}$ and compressive parts $\boldsymbol{\sigma}^{-}$ :

[TABLE]

and

[TABLE]

where the elasticity tensor $\amsmathbb{C}$ has been replaced by the Lamé parameters, $\mu_{s}$ and $\lambda_{s}$. Moreover, $\mathbf{I}$ is the $d\times d$ identity matrix, and $<\cdot>$ is the positive part of a function. In particular, for $d=2$, we have

[TABLE]

where $\lambda_{1}(u)$ and $\lambda_{2}(u)$ are the eigenvalues of the strain tensor $\mathbf{e}:=\mathbf{e}(u)$, and $v_{1}(u)$ and $v_{2}(u)$ the corresponding (normalized) eigenvectors. Finally, the matrix $\mathbf{P}$ is defined as $\mathbf{P}:=\mathbf{P}(u):=[v_{1}|v_{2}]$; namely, it consists of the column vectors $v_{i},\ i=1,2$. We notice that another frequently employed stress-splitting law was proposed in [4].

The modified scheme reads:

[TABLE]

These modifications render the displacement system (4.2a) nonlinear, for which we use a Newton-type solver. The phase field equation is also nonlinear due to the penalization term and the stress splitting. Our version of Newton’s method is based on a residual-based monotonicity criterion (e.g., [17]) outlined in [45][Section 3.2]. Inside Newton’s method, the linear subsystems are solved with a direct solver; namely UMFPACK [16]. All numerical tests presented in Section 5 are implemented in the open-source finite element library deal.II [5, 6]. Specifically, the code is based on a simple adaptation of the multiphysics template [43] in which specifically the previously mentioned Newton solver is implemented.

5 Numerical experiments

In this section, we present several numerical tests to substantiate our algorithmic developments. The goals of all three numerical examples are comparisons between an unlimited number of staggered iterations (although bounded by $500$) denoted by ‘$L$’, and a low, fixed number, denoted by ‘$LFI$’, where we use $30$ (Ex. 1 and Ex. 2), and $20$ (Ex. 3) staggered iterations, respectively. These comparisons are performed in terms of the number of iterations and the correctness of the solutions in terms of the so-called load-displacement curve, measuring the stresses of the top boundary versus the number of loading steps.

5.1 Single edge notched tension test

This test was applied for instance in [31]. The configuration is displayed in Figure 1. We use the system (3.1a)-(3.1b). Specifically, we study our proposed iterative schemes on different mesh levels, denoted as refinement (Ref.) levels 4, 5, 6 (uniformly refined), with $1024$ elements ($2210$ Dofs for the displacements, $1\,105$ Dofs for the phase-field, $h=0.044$), $4\,096$ elements ($8\,514$ Dofs for the displacements, $4257$ Dofs for the phase-field, $h=0.022$), and $16\,384$ elements ($33\,410$ Dofs for the displacements, $16\,705$ Dofs for the phase-field, $h=0.011$).

Specifically, we use $\mu_{s}$ = $80.77\text{\,}\mathrm{k}\mathrm{N}\mathrm{/}\mathrm{m}\mathrm{m}^{2}$, $\lambda_{s}$ = $121.15\text{\,}\mathrm{k}\mathrm{N}\mathrm{/}\mathrm{m}\mathrm{m}^{2}$, and $G_{c}$ = $2.7\text{\,}\mathrm{N}\mathrm{/}\mathrm{m}\mathrm{m}$. The crack growth is driven by a non-homogeneous Dirichlet condition for the displacement field on $\Gamma_{\textnormal{top}}$, the top boundary of $B$ . We increase the displacement on $\Gamma_{\textnormal{top}}$ over time, namely we apply non-homogeneous Dirichlet conditions:

[TABLE]

where $t$ denotes the current loading time. Furthermore, we set $\kappa=10^{-10}$ [mm] and $\varepsilon=2h$ [mm]. We evaluate the surface load vector on the $\Gamma_{\textnormal{top}}$ as

[TABLE]

with normal vector $\nu$, and we are particularly interested in $F_{y}$ for Example 1 and $F_{x}$ for Example 2 (Section 5.2). Graphical solutions are displayed in the Figures 2 and 3 showing the phase-field variable and the discontinuous displacement field. Our findings of using different stabilization parameters $L$ are compared in the Figures 4, 5, 6, 7, and 8. Different mesh refinement studies are shown in the Figures 7 and 8. Here, the number of staggered iterations does not increase with finer mesh levels, which shows the robustness of our proposed methodology.

5.2 Single edge notched shear test

The configuration of this second setting is very similar to Example 1 and was first proposed in a phase-field context in [31]. We now use the model with strain-energy split (4.2a)–(4.2b). The parameters and the geometry (see Figure 1) are the same as in the previous test case. The boundary condition is changed from tensile forces to a shear condition (see also again Figure 1):

[TABLE]

As quanitity of interest we evaluate $F_{x}$ in (5.2). Our findings are shown in the Figures 9, 10, 11, 12, 13, 14, and 15. The major difference to Example 1 is that the scheme is converging even with $L_{u}=0$, as computationally justified in Figure 11. As in Example 1, the load-displacement curves are close to the published literature and, again, the proposed $L$ scheme is robust under mesh refinement (see Figures 12 - 15).

The results are very comparable to the published literature. In particular, it is nowadays known that the proposed Miehe et al. stress splitting does not release all stresses once the specimen is broken (see [1]) and it is also known that we do not see convergence of the curves when both $h$ and $\varepsilon$ are refined (see [25]).

5.3 L-shaped panel

For the configuration of this third example we refer to [1, 29, 44], which are based on an experimental setup [46]. We use again the model with strain-energy split; namely (4.2a)-(4.2b). Moreover, in this test a carefully imposed irreversibility constraint is important since the specimen is pushed, pulled, and again pushed (see Figure 16 for the loading history on the small boundary part $\Gamma_{u}$). In the pulling phase the fracture vanishes if the penalization is not strong enough.

The geometry and boundary conditions are displayed in Figure 1. In contrast to the previous examples, no initial crack prescribed. The initial mesh is $1,2$ and $3$ times uniformly refined, leading to $300,1200,4800$ mesh elements, with $h=$ $29.1548\text{\,}\mathrm{m}\mathrm{m}$, $14.577\text{\,}\mathrm{m}\mathrm{m}$, $7.289\text{\,}\mathrm{m}\mathrm{m}$, respectively.

We increase the displacement $u_{D}:=u_{y}=u_{y}(t)$ on $\Gamma_{u}:=\{(x,y)\in B|\,$470\text{,}\mathrm{m}\mathrm{m}$\leq x\leq$500\text{,}\mathrm{m}\mathrm{m}$,\,y=$250\text{,}\mathrm{m}\mathrm{m}$\}$ over time, where $\Gamma_{u}$ is a section of $30\text{\,}\mathrm{m}\mathrm{m}$ length on the right corner of the specimen. We apply a loading-dependent, non-homogeneous Dirichlet condition (see also Figure 16):

[TABLE]

where $t$ denotes the total loading time. Due to this cyclic loading the total displacement at the end time $T$ = $2\text{\,}\mathrm{s}$ is $1\text{\,}\mathrm{m}\mathrm{m}$.

We use $\mu_{s}$ = $10.95\text{\,}\mathrm{k}\mathrm{N}\mathrm{/}\mathrm{m}\mathrm{m}^{2}$, $\lambda_{s}$ = $6.16\text{\,}\mathrm{k}\mathrm{N}\mathrm{/}\mathrm{m}\mathrm{m}^{2}$, and $G_{c}$ = $8.9\text{\times}{10}^{-5}\text{\,}\mathrm{k}\mathrm{N}\mathrm{/}\mathrm{m}\mathrm{m}$. The time (loading) step size is $\delta t$ =${10}^{-3}\text{\,}\mathrm{s}$. Furthermore, we set $k=10^{-10}h$[mm] and $\varepsilon=2h$. As before, we observe the number of Newton iterations and we evaluate the surface load vector

[TABLE]

with normal vector $\nu$, and now we are particularly interested in $F_{y}$. The crack path at the chosen time step snapshots in Figure 17 corresponds to the published literature [44, 29, 1]. The load-displacement curves and the number of iterations for different $L$ and corresponding mesh refinement studies are displayed in the Figures 18, 19, 20, 21 and 22.

5.4 Verification of Assumption 1

In this last set of computations, we verify whether Assumption 1 holds true in our computations. We choose some prototype settings, namely on the coarest mesh level Ref. 4 and $L_{u}=L_{\varphi}=1e-6$. In Figure 23, we observe that $\operatorname{ess}\sup_{x\in B}|\mathbf{e}(u^{n}(x))|$ varies, but always can be bounded from above with $M>0$. The value of $\operatorname{ess}\sup_{x\in B}|\mathbf{e}(u^{n}(x))|$ is the final strain when the $L$-scheme terminates. The minimum and maximum values shows that there are no significant variations in $\operatorname{ess}\sup_{x\in B}|\mathbf{e}(u^{n}(x))|$ during the $L$-scheme iterations with respect to the finally obtained value.

6 Conclusions

We have proposed a novel staggered iterative algorithm for brittle fracture phase field models. This algorithm is employing stabilization and linearization techniques known in the literature as the ‘$L$-scheme’, which is a generalization of the Fixed Stress Splitting algorithm coming from poroelasticity. Through theory and numerical examples we have investigated the performance of our proposed variants of the $L$-scheme for brittle fracture phase field problems.

Under natural constraints that the elastic mechanical energy remains bounded, and that the model parameter $\varepsilon$ is sufficiently large (i.e., that the diffusive zone around crack surfaces must be sufficiently thick), we have shown that a contraction of successive difference functions in energy norms can be obtained from the proposed scheme. This result implies the algorithm is converging monotonically with a linear convergence rate. However, in the convergence analysis there appears some unknown constants which makes the precise convergence rate, as well as the precise lower bound on $\varepsilon$ unknown.

We provide detailed numerical tests where our proposed scheme is employed on several phase field brittle fracture bench-mark problems. For each numerical example we provide findings for different values of stabilization parameters. For most cases we let $L_{u}=L_{\varphi}>0$, but for comparison we include also for the stabilization configurations $L_{u}=0$ with $L_{\varphi}>0$, and $L_{u}=L_{\varphi}=0$. For the test cases presented here, there is only Example 1 where $L_{u}=0$ does not work. This might be due to the very rapid crack growth, which sets Example 1 apart from Examples 2 and 3. In this regard, we conclude that further work is needed to find an optimal configuration of $L_{u}$ and $L_{\varphi}$. For all numerical test we also provide computational justification for the assumption of bounded elastic mechanical energy. Furthermore, a slight dependency on $h$ in the iteration counts is observed in the numerical tests, but this is expected since we use $\varepsilon=2h$, and as our analysis demonstrates, the convergence rate is dependent on $\varepsilon$. The variation in iteration numbers with mesh refinement is in any case sufficiently small enough that we conclude our algorithm is robust with respect to mesh refinement.

Moreover, due to the iteration spikes at the critical loading steps, we have included, for comparison, several results in which the iteration has been truncated (labeled $LFI$ in Examples 1-3). Due to the monotonic convergence of the scheme, this strategy still produces acceptable results, while effectively avoiding the iteration spikes. We therefore conclude, at least for the particular examples presented here, that a truncation of the $L$-scheme can be employed for greatly improved efficiency with only negligible (depending on the situation at hand, of course) loss of accuracy.

Acknowledgements

This work forms part of Research Council of Norway project 250223. The authors also acknowledges the support from the University of Bergen. The first author, MKB, thanks the group ‘Wissenschaftliches Rechnen’ of the Institute of Applied Mathematics of the Leibniz University Hannover for the hospitality during his research stay from Oct - Dec 2018. The second author, TW, has been supported by the German Research Foundation, Priority Program 1748 (DFG SPP 1748) named Reliable Simulation Techniques in Solid Mechanics. Development of Non-standard Discretization Methods, Mechanical and Mathematical Analysis. The subproject within the SPP1748 reads Structure Preserving Adaptive Enriched Galerkin Methods for Pressure-Driven 3D Fracture Phase-Field Models (WI 4367/2-1).