Moving mesh finite element simulation for phase-field modeling of brittle fracture and convergence of newton's iteration

TL;DR

This paper introduces a moving mesh finite element method combined with regularization techniques to improve the numerical simulation and convergence of phase-field models for brittle fracture, effectively capturing complex crack propagation.

Contribution

The study develops a novel moving mesh finite element approach with regularization methods to enhance convergence and accuracy in phase-field modeling of brittle fracture.

Findings

Regularization improves Newton's iteration convergence for small parameters.

Moving mesh adapts to complex crack systems and propagating cracks.

Methods maintain solution accuracy while enhancing numerical stability.

Abstract

A moving mesh finite element method is studied for the numerical solution of a phase-field model for brittle fracture. The moving mesh partial differential equation approach is employed to dynamically track crack propagation. Meanwhile, the decomposition of the strain tensor into tensile and compressive components is essential for the success of the phase-field modeling of brittle fracture but results in a non-smooth elastic energy and stronger nonlinearity in the governing equation. This makes the governing equation much more difficult to solve and, in particular, Newton's iteration often to fail to converge. Three regularization methods are proposed to smooth out the decomposition of the strain tensor. Numerical examples of fracture propagation under quasi-static load demonstrate that all of the methods can effectively improve the convergence of Newton's iteration for relatively small values of the regularization parameter but without comprising the accuracy of the numerical solution. They also show that the moving mesh finite element method is able to adaptively concentrate the mesh elements around propagating cracks and handle multiple and complex crack systems.

Figures

Tables

Table 1. Table 1: Average CPU time (in seconds) for one time step for solving for d 𝑑 d , u 𝑢 u and mesh.

Mesh size	CPU for $d$		CPU for $u$		CPU for mesh		Total CPU time
	time	%	time	%	time	%
Moving mesh $41\times 41$	5.39	3.82	31.89	22.61	103.79	73.58	141.07
Uniform mesh $351\times 351$	47.9	3.58	1291.8	96.4			1339.7

Equations

$$ϵ=\frac{1}{2}\left(\nabla u+(\nabla u{)}^T\right),$$

$$W_e(ϵ)=\frac{\lambda }{2}{\left(\mathrm{tr}(ϵ)\right)}^2+\mu \mathrm{tr}({ϵ}^2),$$

$$\sigma :=\frac{\partial W_e}{\partial ϵ}=\lambda \mathrm{tr}(ϵ)I+2\mu ϵ,$$

$${\mathcal{W}}_e(ϵ)={\int }_{\Omega }{W}_e(ϵ)d\Omega .$$

$$\delta {\mathcal{W}}_e={\int }_{\Omega }\frac{\partial W_e}{\partial ϵ}:\delta ϵd\Omega ={\int }_{\Omega }\sigma :ϵ(\delta u)d\Omega ,$$

$$V_u=\left\{\phi |\phi \in H^1(\Omega ),\phi =\overline{u}\text{on}\partial {\Omega }_u\right\},$$

$$V_u^0=\left\{\phi |\phi \in H^1(\Omega ),\phi =0\text{on}\partial {\Omega }_u\right\},$$

$$H^1(\Omega )=\left\{\phi |{\int }_{\Omega }{\phi }^{2}d\Omega <+\infty ,{\int }_{\Omega }(\nabla \phi {)}^{2}d\Omega <+\infty \right\}.$$

$${\int }_{\Omega }\sigma :ϵ(\delta u)d\Omega ={\int }_{\partial {\Omega }_t}\overline{t}\cdot \delta udS+{\int }_{\Omega }f\cdot \delta ud\Omega ,\forall \delta u\in V_u^0.$$

$$\mathcal{W}(ϵ,\Gamma )={\mathcal{W}}_e(ϵ,\Gamma )+{\mathcal{W}}_c(\Gamma ):={\int }_{\Omega \setminus \Gamma }{W}_e(ϵ,\Gamma )d\Omega +{\int }_{\Gamma }{g}_{c}dS,$$

$${\mathcal{W}}_c^l={\int }_{\Omega }\frac{g_c}{4l}\left((d-1{)}^2+4l^2|\nabla d{|}^2\right)d\Omega .$$

$$f=f^{+}+f^{-},f^{+}=\frac{f+|f|}{2},f^{-}=\frac{f-|f|}{2}.$$

$$ϵ={ϵ}^{+}+{ϵ}^{-},{ϵ}^{+}=Q\text{diag}({\lambda }_1^{+},...,{\lambda }_n^{+})Q^T,{ϵ}^{-}=Q\text{diag}({\lambda }_1^{-},...,{\lambda }_n^{-})Q^T.$$

$$W_e=g(d){\Psi }^{+}(ϵ)+{\Psi }^{-}(ϵ),$$

$${\Psi }^{+}(ϵ)=\frac{\lambda }{2}{\left((\mathrm{tr}(ϵ){)}^{+}\right)}^2+\mu \mathrm{tr}(({ϵ}^{+}{)}^2),{\Psi }^{-}(ϵ)=\frac{\lambda }{2}{\left((\mathrm{tr}(ϵ){)}^{-}\right)}^2+\mu \mathrm{tr}(({ϵ}^{-}{)}^2).$$

$$\{\begin{array}{ll}g(0)=0:& \text{Damage occurred for d=0 and this part should vanish};\\ g(1)=1:& \text{No damage occurs for d=1};\\ g^{\prime }(0)=0:& \text{No more changes after the fully broken state};\\ g^{\prime }(1)\ne 0:& \text{The damage has to be initiated at the onset.}\end{array}$$

$${\mathcal{W}}^l={\mathcal{W}}_e+{\mathcal{W}}_c^l={\int }_{\Omega }\left((d^2+k_l){\Psi }^{+}(ϵ)+{\Psi }^{-}(ϵ)+\frac{g_c}{4l}\left((d-1{)}^2+4l^2|\nabla d{|}^2\right)\right)d\Omega ,$$

$$\delta {\mathcal{W}}^l={\int }_{\Omega }\frac{\partial W}{\partial d}\delta dd\Omega +{\int }_{\Omega }\frac{\partial W}{\partial \nabla d}\cdot \nabla \delta dd\Omega +{\int }_{\Omega }\frac{\partial W}{\partial ϵ}:ϵ(\delta u)d\Omega$$

$$W=(d^2+k_l){\Psi }^{+}(ϵ)+{\Psi }^{-}(ϵ)+\frac{g_c}{4l}\left((d-1{)}^2+4l^2|\nabla d{|}^2\right).$$

$${\int }_{\Omega }\left(\left(2d\mathcal{H}+\frac{g_c(d-1)}{2l}\right)\delta d+2g_{c}l\nabla d\cdot \nabla \delta d\right)d\Omega =0,\forall \delta d\in V_d$$

$${\int }_{\Omega }\sigma :ϵ(\delta u)d\Omega ={\int }_{{\Omega }_t}\overline{t}\cdot \delta udS+{\int }_{\Omega }f\cdot \delta ud\Omega ,\forall \delta u\in V_u^0$$

$$\sigma :=\frac{\partial W}{\partial ϵ}=(d^2+k_l)\left(\lambda (\mathrm{tr}(ϵ){)}^{+}I+2\mu {ϵ}^{+}\right)+\left(\lambda (\mathrm{tr}(ϵ){)}^{-}I+2\mu {ϵ}^{-}\right).$$

$$\mathcal{H}=\underset{s\le t}{\max }{\Psi }^{+}(ϵ)(s),$$

$$V_d^h=\left\{{\phi }_h|{\phi }_h\in C^0(\overline{\Omega })\cap V_d;{\phi }_h{|}_K\in P_1(K),\forall K\in {\mathcal{T}}_h(\Omega )\right\},$$

$$V_u^h=\left\{{\phi }_h|{\phi }_h\in C^0(\overline{\Omega })\cap V_u;{\phi }_h{|}_K\in P_1(K),\forall K\in {\mathcal{T}}_h(\Omega )\right\},$$

$$V_u^{0,h}=\left\{{\phi }_h|{\phi }_h\in C^0(\overline{\Omega })\cap V_u^0;{\phi }_h{|}_K\in P_1(K),\forall K\in {\mathcal{T}}_h(\Omega )\right\},$$

$${\int }_{\Omega }\left(\left(2d_h\mathcal{H}+\frac{g_c(d_h-1)}{2l}\right){\phi }_h+2g_{c}l\nabla d_h\cdot \nabla {\phi }_h\right)d\Omega =0,\forall {\phi }_h\in V_d^h$$

$${\int }_{\Omega }\sigma (u_h):ϵ({\phi }_h)d\Omega ={\int }_{{\Omega }_t}\overline{t}\cdot {\phi }_{h}dS+{\int }_{\Omega }f\cdot {\phi }_{h}d\Omega ,\forall {\phi }_h\in V_u^{0,h}.$$

$${\lambda }_{\alpha }^{+}=\frac{\lambda +\sqrt{{\lambda }^2+{\alpha }^2}}{2},{\lambda }_{\alpha }^{-}=\frac{\lambda -\sqrt{{\lambda }^2+{\alpha }^2}}{2},$$

$${\delta }_{\alpha }(\lambda )=\frac{1}{\sqrt{2\pi }\alpha }{e}^{-\frac{{\lambda }^2}{2{\alpha }^2}}.$$

$${\lambda }_{\alpha }^{+}={\int }_{-\infty }^{\infty }{\lambda }^{+}(\eta ){\delta }_{\alpha }(\lambda -\eta )d\eta =\frac{\lambda }{2}\left(1+\text{erf}\left(\frac{\lambda }{\sqrt{2}\alpha }\right)\right)+\frac{\alpha }{\sqrt{2\pi }}{e}^{-\frac{{\lambda }^2}{2{\alpha }^2}},$$

$${\lambda }_{\alpha }^{-}={\int }_{-\infty }^{\infty }{\lambda }^{-}(\eta ){\delta }_{\alpha }(\lambda -\eta )d\eta =\frac{\lambda }{2}\left(1+\text{erf}\left(-\frac{\lambda }{\sqrt{2}\alpha }\right)\right)-\frac{\alpha }{\sqrt{2\pi }}{e}^{-\frac{{\lambda }^2}{2{\alpha }^2}}.$$

$${\delta }_{\alpha }(\lambda )=\{\begin{array}{ll}\frac{1}{\alpha }\left(\frac{3}{4}-\frac{{\lambda }^2}{{\alpha }^2}\right),& \text{ for }|\lambda |⩽0.5\alpha \\ \frac{1}{\alpha }\left(\frac{9}{8}-\frac{3}{2}\frac{|\lambda |}{\alpha }+\frac{1}{2}\frac{{\lambda }^2}{{\alpha }^2}\right),& \text{ for }0.5\alpha ⩽|\lambda |⩽1.5\alpha \\ 0,& \text{ for }|\lambda |⩾1.5\alpha .\end{array}$$

Full text

AMS 2010 Mathematics Subject Classification. 65M50, 65M60, 74B99

Key Words. brittle fracture, phase-field model, Newton’s iteration, moving mesh, mesh adaptation, finite element method

Abbreviated title. Moving Mesh Finite Element Simulation of Brittle Fracture

1 Introduction

Brittle fracture is the fracture of a metallic object or other elastic material where plastic deformation is strongly limited. It usually occurs very rapidly and can be catastrophic in engineering practice; e.g., see Pokluda and Šandera [48]. Understanding the initiation and propagation of brittle fracture and preventing fracture failure are vital to the engineering design, where numerical simulation of fracture processes has become a powerful tool. Computational approaches for studying brittle fracture can be roughly categorized into two groups, discrete crack models and smeared crack models. In the former group, discontinuous fields are introduced into the numerical model and cracks are described as moving boundaries. One major challenge for those models is to track moving boundaries. A commonly used strategy is to change the mesh geometry by introducing new boundaries at each time step together with adaptive remeshing; e.g., see [6, 12, 30, 47]. The mesh regenerating and boundary updating not only increase computational cost but also further complicate the implementation of boundary conditions. In order to avoid complex remeshing, Moës et al. [43, 44] propose the extended finite element method, which enriches the finite element spaces with discontinuous fields based on the partition-of-unity concept and allows the propagation of cracks along element interfaces. A drawback for the method is that it requires explicit description of crack patterns and thus has difficulty in dealing with complex cracks and unforeseen patterns of crack propagation.

In the second group of computational approaches, smeared crack models approximate cracks with continuous fields and do not rely on explicit description of cracks. The phase-field model based on the variational approach proposed by Francfort and Marigo [17] is a commonly used type of smeared crack model. In the phase-field modeling, a phase-field variable $d$, which depends on a parameter $l$ describing the actual width of the smeared cracks, is introduced to indicate where the material is damaged. One of the major advantages of this model is that the initiation and propagation of cracks are completely determined by a coupled system of partial differential equations based on the energy functional. Another advantage is that the generation and propagation of fracture networks do not require explicitly tracking fracture interfaces. The phase-field modeling is used in this work.

The phase-field modeling has been successfully applied in many other fields including image segmentation [2], dendritic crystal growth [31, 56], and multiple-fluid hydrodynamics [34, 50, 51, 59]. Since its first application in brittle fracture simulation by Bourdin et al. [13], significant progress has been made in this area; e.g., see [3, 9, 11, 32, 35, 38, 40, 42, 45, 53]. However, there still exist challenges. In particular, the strain tensor has to be decomposed along eigen-directions into tensile and compressive components in the presence of cracks, with only the former component contributing to generation and propagation of cracks. This decomposition of the strain tensor is introduced by Miehe et al. [40] to account for the reduction of the stiffness of the elastic solid by cracks and to rule out unrealistic branching. Unfortunately, it also makes the elastic energy non-smooth and increases the nonlinearity of the governing equation. As a consequence, the Jacobian matrix of the governing equation may not exist at places and Newton’s iteration can often fail to converge [35, 45].

The phase-field modeling is governed by a coupled system for the phase-field variable $d$ and the displacement $u$ which can be solved in the monolithic or staggered approach. In the monolithic approach (e.g., see [19, 57, 58]), the system is solved simultaneously for both $d$ and $u$, and it is often challenging to obtain a convergent solution due to the non-convexity nature of the energy functional. A damped Newton method with linear search [19] and an error-oriented Newton method [57] have been developed to overcome convergence issues while a modified Newton scheme with Jacobian modification has been proposed by Wick [58]. The staggered approach solves the coupled system sequentially for $d$ and $u$ and has been used by a number of researchers; e.g., see [1, 3, 18, 40, 42]. By fixing one variable such as the phase-field $d$, the underlying problem becomes convex in the other unknown variable $u$. Since $d$ and $u$ are not coupled, the procedure becomes simpler and more robust. The main disadvantage of this approach is that at a loading step many staggered iterations are required to reach convergence, which can be costly. The effects of the number of staggered iterations in relation to the size of the loading increment has been studied by Ambati et al. [1]. Their results show that an insufficient number of $d$-$u$ iterations can lead to inaccurate results when large loading increments are used. However, the number of staggered iterations usually has relatively little performance impact on the shape of the load-displacement curves when loading increments are sufficiently small. This implies that the staggered approach without iteration can be used as long as small loading increments are taken. Since our focus in this work is on mesh adaptation and convergence of Newton’s iteration, we use the staggered approach without iteration for quasi-static brittle fracture problems with small loading increments.

It is worth mentioning that for the staggered approach with/without iteration, the equation for $u$ remains highly nonlinear due to the decomposition of the strain tensor, which can make Newton’s iteration fail or be slow to converge; see [1] or the numerical results in Section 4. On the other hand, like implicit schemes versus explicit schemes for ordinary differential equations, the staggered approach with iteration or the monolithic approach can be more robust than the staggered approach without iteration which is used in the current work. Comparison of their performances with mesh adaptation and regularization of the decomposition of the strain tensor (see discussion below) can be an interesting topic for future investigations.

The model parameter that describes the width of smeared cracks is needed to be small for the phase-field model to be a reasonably accurate approximation of the original problem. This in turn requires that the mesh elements be small at least in the crack regions, meaning that mesh adaptation is necessary to improve the computational accuracy and efficiency. The mesh adaptation should be dynamical too since cracks can propagate under continuous load.

Moving mesh methods are well suited for the numerical simulation of the phase-field modeling of brittle fracture. Although they have been successfully applied to phase-field models for other applications, e.g., see [16, 36, 49, 54, 60, 61], they have not been employed for brittle fracture simulation, which has distinguished challenges associated with the above mentioned decomposition of the strain tensor. As a matter of fact, mesh adaptation has rarely been employed in brittle fracture simulation so far and there are only a few published studies on the topic. Noticeably, Heister et al. [20] develop a predictor-corrector local mesh adaptivity scheme that allows the mesh to refine around cracks. Artina et al. [5] present an a posteriori error estimator for anisotropic mesh adaptation that generates thin, anisotropic elements around cracks and isotropic elements away from the cracks, but they use an early model that does not decompose the strain tensor and therefore does not distinguish between fracture caused by tension and compression.

The objective of this paper is twofold. The first is to study the MMPDE (moving mesh partial differential equation) moving mesh method [14, 26, 27, 28] for the phase-field modeling of brittle fracture. The MMPDE method is a type of dynamic mesh adaptation method specially designed for time dependent problems. It employs a mesh PDE to move the mesh continuously in time to follow and adapt to evolving structures in the solution. A new formulation of the MMPDE method was developed recently in [24], which provides a simple, compact analytical formula for the nodal mesh velocities (cf. (23) below), and this makes its implementation relatively easy. It is also shown in [25] that the mesh governed by the underlying mesh equation stays nonsingular if it is nonsingular initially. The MMPDE moving mesh method is combined with a piecewise linear finite element discretization to solve the governing equation in a quasi-static condition, with the quasi-time being introduced to represent the load increments.

It is noted that a number of other moving mesh methods have also been developed in the past and there is a large literature in the area. The interested reader is referred to the books or review articles [7, 8, 14, 28, 52] and references therein.

The second goal of the paper is to investigate the convergence of Newton’s iteration used for solving the nonlinear algebraic system resulting from the finite element discretization of the displacement equation.

As mentioned above, Newton’s iteration often fails to converge due to the non-existence of the Jacobian matrix and nonlinearity of the governing equation caused by the decomposition of the strain tensor. To avoid this difficulty, we propose to smooth out the decomposition with regularization and discuss three regularization methods. Numerical examples show that all of them can effectively improve the convergence of Newton’s iteration for relatively small values of the regularization parameter but without comprising the accuracy of the numerical solution. Note that no special treatment is needed for Newton’s iteration to guarantee its convergence. This is in contrast to previous works such as [19, 57, 58] where modified Newton’s schemes are used along with some modification in the computation of the Jacobian matrix. Our numerical results also show that the proposed moving mesh finite element method improves the computational accuracy and efficiency and is able to handle multiple and complex crack systems.

The rest of this paper is organized as follows. Section 2 provides a brief introduction of the phase-field model for brittle fracture. The moving mesh finite element method and the three regularization methods for the decomposition of the strain tensor are described in Section 3. Numerical results obtained for three two-dimensional examples are presented in Section 4, and conclusions are drawn in Section 5.

2 Phase-field models for brittle fracture

2.1 Variational approach to elasticity models without cracks

We first describe small strain isotropic elasticity models without cracks. We consider an elastic body occupying a bounded domain $\Omega\subset\mathbb{R}^{2}$ with the boundary $\partial\Omega=\partial\Omega_{t}\cup\partial\Omega_{u}$, where the surface traction $\overline{t}$ is specified on $\partial\Omega_{t}$ and the displacement $\overline{u}$ is given on $\partial\Omega_{u}$. The strain tensor is given by

[TABLE]

where $\nabla{u}$ is the displacement gradient tensor. For isotropic material without damage, the elastic energy per unit volume, or strain energy density, is given by Hooke’s law as

[TABLE]

where $\lambda$ and $\mu$ are the Lamé constants. The stress tensor is given by the derivative of the strain energy with respect to the strain tensor, i.e.,

[TABLE]

where the symbol $:=$ stands for “definition”. Then the total strain energy stored in the elastic solid is given by

[TABLE]

The variation of $\mathcal{W}_{e}$ is

[TABLE]

where $A:B$ is the inner product of tensors $A$ and $B$, i.e., $A:B=\sum\limits_{i,j}A_{i,j}B_{i,j}$. Define the function spaces

[TABLE]

where $H^{1}(\Omega)$ is a Sobolev space, viz.,

[TABLE]

If we add the boundary traction $\overline{t}$ and body force $f$, then the variational formulation of the elasticity model is to find $u\in V_{u}$ such that

[TABLE]

2.2 Phase-field approach to elasticity models with cracks

We now consider the situation with cracks in the elastic body. We use the approach proposed by Francfort and Marigo [17] where the total energy of the body with a given crack $\Gamma$ is given by

[TABLE]

where $\mathcal{W}_{e}({\epsilon},\Gamma)$ represents the energy stored in the bulk of the elastic body, $\mathcal{W}_{c}(\Gamma)$ is the energy required to create the crack according to the Griffith criterion, and $g_{c}$ is the fracture energy density (also referred to as the fracture toughness) which is the amount of energy needed to create a unit area of fracture surface.

In the phase-field modeling, the fracture surface is approximated by a phase-field variable $d(x,t)$, which depends on a parameter $l$ describing the width of the smooth approximation of the crack. This function is smooth with the value 0 or close to 0 near the crack and 1 away from the crack (see Fig. 1). The fracture energy $\mathcal{W}_{c}(\Gamma)$ is approximated by the smeared total fracture energy [13] as

[TABLE]

The elastic energy $\mathcal{W}_{e}({\epsilon},\Gamma)$ needs to be modified to reflect the loss of material stiffness in the damage zone. We follow the approach by Miehe et al. [40]. Define the decomposition of a scalar function $f$ as

[TABLE]

Assuming that the strain tensor has the eigen-decomposition $\epsilon=Q\text{diag}(\lambda_{1},...,\lambda_{n})Q^{T}$, we define

[TABLE]

The tensile strain component, $\epsilon^{+}$, contributes to the damage process resulting in crack initiation and propagation whereas the compression strain component, $\epsilon^{-}$, does not contribute to the damage process. A commonly used damage model is given as

[TABLE]

where $g(d)$ is a degradation function that describes the reduction of the stiffness of the bulk of the solid and

[TABLE]

The degradation function $g(d)$ is required to satisfy the property (e.g., see [33])

[TABLE]

We take a commonly used quadratic degradation function $g(d)=d^{2}$ in our computation. Combining (6) and (8), we obtain the total energy as

[TABLE]

where $k_{l}\ll l$ is the (small) regularization parameter for avoiding degeneracy. The variation of the energy is

[TABLE]

where

[TABLE]

The weak formulation is to find $d\in V_{d}=H^{1}(\Omega)$ (with $d$ satisfying a homogeneous Neumann boundary condition) and $u\in V_{u}$ such that

[TABLE]

where $\mathcal{H}=\Psi^{+}(\epsilon)$ and

[TABLE]

To ensure crack irreversibility in the sense that the cracks can only grow, we replace $\mathcal{H}=\Psi^{+}(\epsilon)$ in (10) by

[TABLE]

where $t$ is the quasi-time corresponding to the load increments.

In addition to the above described decomposition by Miehe et al. [40], there is another often-used decomposition by Amor et al. [3] that is based on the split of the strain tensor into volumetric and deviatoric parts. Discussion of the advantages and disadvantages of both decompositions can be found in [1, 9].

3 The moving mesh finite element approximation

3.1 Finite element discretization and solution procedure

We consider a simplicial mesh $\mathcal{T}_{h}$ for the domain $\Omega$ and denote the number of its elements and vertices by $N$ and $N_{v}$, respectively. The piecewise linear approximations of the function spaces $V_{d}$, $V_{u}$ and $V_{u}^{0}$ are given by

[TABLE]

where $P_{1}(K)$ is the set of polynomials of degree less than or equal to 1 defined on $K$. The linear finite element approximation of the phase-field problem (10) and (11) is to find $d_{h}\in V_{d}^{h}$ and $u_{h}\in V_{u}^{h}$ such that

[TABLE]

Note that $u_{h}$ and $d_{h}$ are strongly coupled through (14) and (15). Common procedures for solving these equations can be roughly categorized into two groups, simultaneous (or called monolithic) solution and alternating (or called staggered) solution. In the former group (such as see [19, 20, 57, 58]) the phase-field variable and displacement are solved simultaneously often by Newton’s method. This procedure has to deal with a large system and the highly nonlinear coupling between $d_{h}$ and $u_{h}$. Moreover, the decomposition of the strain tensor further increases the nonlinearity (see Section 3.2 for detailed discussion). It is still challenging to obtain a convergent solution using the simultaneous solution procedure. The second solution procedure (such as see [1, 3, 40, 42]) is to solve (14) for $d_{h}$ and (15) for $u_{h}$ alternatingly. One of the advantages of this approach is that $d_{h}$ and $u_{h}$ are not coupled. Moreover, (14) is linear about $d_{h}$ and its solution does not need Newton’s iteration.

In our computation, we use the alternating solution procedure and consider the problem in a quasi-static condition. The quasi-time $t$ is introduced to represent the load increments. The solution procedure from $t^{n}$ to $t^{n+1}$ is described as follows.

1. (i)

We assume that the mesh $\mathcal{T}_{h}^{n}$ at time $t^{n}$ and the history field $\mathcal{H}_{h}^{n}$ in (13) (defined on $\mathcal{T}_{h}^{n}$) are given. 2. (ii)

Compute the phase-field variable $d_{h}^{n+1}$ and new mesh $\mathcal{T}_{h}^{n+1}$ as follows.

• •

Set $\mathcal{T}_{h}^{n+1,1}=\mathcal{T}_{h}^{n}$;

• •

For $k=1:kk$

Compute $\mathcal{H}$ on $\mathcal{T}_{h}^{n+1,k}$ using linear interpolation of $\mathcal{H}_{h}^{n}$;

• -

Compute $d_{h}^{n+1,k}$ using (14) and $\mathcal{H}$ on $\mathcal{T}_{h}^{n+1,k}$;

• -

If $k<kk$, compute the new mesh $\mathcal{T}_{h}^{n+1,k+1}$ by the MMPDE moving mesh method based on $\mathcal{T}_{h}^{n+1,k}$ and $d_{h}^{n+1,k}$. (See Section 3.3)

• •

Set $\mathcal{T}_{h}^{n+1}=\mathcal{T}_{h}^{n+1,kk}$ and $d_{h}^{n+1}=d_{h}^{n+1,kk}$. 3. (iii)

Compute the displacement field $u_{h}^{n+1}$ by solving the nonlinear system (15) based on $d_{h}^{n+1}$ and $\mathcal{T}_{h}^{n+1}$. (See Section 3.2) 4. (iv)

Compute $\Psi^{+,n+1}_{h}(\epsilon(u_{h}^{n+1}))$ and interpolate $\mathcal{H}_{h}^{n}$ from the old mesh $\mathcal{T}_{h}^{n}$ to the new mesh $\mathcal{T}_{h}^{n+1}$ using linear interpolation and denote it by $\tilde{\mathcal{H}}_{h}^{n+1}$. Let $\mathcal{H}_{h}^{n+1}=\max\{\Psi^{+,n+1}_{h}(\epsilon(u_{h}^{n+1})),\tilde{\mathcal{H}}_{h}^{n+1}\}$.

In the above procedure, the parameter $kk$ can affect the adaptivity of the mesh $\mathcal{T}_{h}^{n+1}$ to the phase-field variable $d_{h}^{n+1}$. From our computational experience, we choose $kk=5$, for which we have found that $\mathcal{T}_{h}^{n+1}$ is sufficiently adaptive to $d_{h}^{n+1}$ but without compromising the computational efficiency too much.

3.2 Convergence of Newton’s iteration

We recall that the strain tensor is decomposed into tensile and compressive components (cf. (7)), with only the former component contributing to generation of cracks. The tensile and compressive components are nonsmooth functions of the displacement, which has two major effects on Newton’s iteration applied to the solution of the equation (15), the existence and computation of the Jacobian matrix and the convergence of Newton’s iteration. As can be seen in Fig. 2(a) and Fig. 2(b), the Jacobian matrix of equation (15) does not exist in the places where any of the eigenvalues vanishes. In other places where the eigenvalues do not vanish, the computation of the Jacobian matrix can also be tricky. In principle, the analytical expressions of the derivatives of $\epsilon^{+}$ and $\epsilon^{-}$ with respect to $\epsilon$ can be obtained through the derivatives of the eigenvalues and eigenvectors with respect to $\epsilon$ (e.g, see [37]). However, these expressions are too complicated to be used in practical computation. Special algorithms are needed to compute these derivatives using analytical formulas; see Miehe [39] and Miehe and Lambrecht [41].

Another effect of the decomposition of the strain tensor is that Newton’s iteration often fails to converge; see the numerical examples in Section 4.1.1. The degeneracy of the equation (15) can further complicate the situation. In principle, the regularization parameter $k_{l}$ can be chosen sufficiently large to make Newton’s iteration to converge. However, a large value of $k_{l}$ will smear the crack and lead to an unphysically overestimated bulk energy.

To overcome the above mentioned difficulties, we propose to smooth out the decomposition (7) and compute the Jacobian matrix of the equation (15) using finite differences. The latter is straightforward so we will not elaborate it here. For the former, there are a variety of regularization methods. We consider here three of these methods where the positive and negative eigenvalue functions are smoothed using a switching technique or convolution with a smoothed delta function (or a mollifier).

• •

The sonic-point regularization method. We first consider the so-called eigenvalue-switching technique which is used to obtain a smooth transition of the solution through the sonic point [4] in computational fluid dynamics. The regularized positive and negative eigenvalue functions are defined as

[TABLE]

where $\alpha>0$ is the regularization parameter.

• •

The exponential convolution method. The second regularization method is to take the convolution of $\lambda^{+}$ and $\lambda^{-}$ with the exponential delta function

[TABLE]

Then the regularized positive and negative eigenvalue functions are given by

[TABLE]

• •

The smoothed 2-point convolution method. The third method is similar to the exponential convolution method but with a localized, smoothed 2-point delta function

[TABLE]

We then have

[TABLE]

It is remarked that all of the above regularization methods satisfy

[TABLE]

and thus

[TABLE]

Moreover, as can be seen in Fig. 2(c), all of the methods provide a smooth transition around zero for both positive and negative eigenvalue functions. Furthermore, the sonic-point and exponential convolution methods have global effects whereas the smoothed 2-point convolution method only changes the values near $\lambda=0$. Finally, for a fixed value of the regularization parameter, $\alpha=3\times 10^{-3}$, in terms of the closeness of the curve of $\lambda_{\alpha}^{+}$ ($\lambda_{\alpha}^{-}$) to that of $\lambda^{+}$ (or $\lambda^{-}$) the best is the smoothed 2-point convolution and then the exponential convolution and sonic-point methods.

The effects of these regularization methods on the the convergence of Newton’s iteration and on the numerical solution will be discussed in Section 4.1.1. It is worth mentioning that we have tried a few other regularization methods which only modify the local behavior of the eigenvalue function near the origin. Since they do not satisfy (16) in general and are much less effective in making Newton’s iteration convergent than the above methods, we choose not to report them here to save space.

3.3 The MMPDE moving mesh method

In this section we describe the MMPDE moving mesh method [26, 28] for generating the new mesh $\mathcal{T}_{h}^{n+1}$ (as mentioned in Section 3.1). The method takes the $\mathbb{M}$-uniform mesh approach where a nonuniform mesh is viewed as a uniform one in the metric specified by a tensor $\mathbb{M}$. In our computation, we choose the metric tensor as

[TABLE]

where $H(d_{h})$ is a recovered Hessian of $d_{h}$ and $|H(d_{h})|=Q\text{diag}(|\lambda_{1}|,...,|\lambda_{2}|)Q^{T}$, assuming that the eigen-decomposition of $H(d_{h})$ is $Q\text{diag}(\lambda_{1},...,\lambda_{2})Q^{T}$. The recovered Hessian of $d_{h}$ at a vertex is obtained by twice differentiating a local quadratic polynomial fitting in the least-squares sense to the nodal values of $d_{h}$ at the neighboring vertices. The form of (18) is known optimal in terms of the $L^{2}$ norm of linear interpolation error on triangular meshes (e.g., see [29]). Notice that $\mathbb{M}$ is symmetric and uniformly positive define on $\Omega$. It will be used to control the shape, size, and orientation of mesh elements through the so-called equidistribution and alignment conditions (to be discussed below). Since (18) is based on the Hessian of the phase-field variable $d$, we may expect that the mesh elements are concentrated in the crack regions where the curvature of $d$ is large.

Let $x_{1},....,x_{N_{v}}$ be the coordinates of the vertices of $\mathcal{T}_{h}$. We choose the very initial physical (simplicial) mesh as the reference computational mesh which is denoted by $\hat{\mathcal{T}}_{c,h}=\{\hat{\xi}_{1},...,\hat{\xi}_{N_{v}}\}$. For the purpose of mesh generation, we need an intermediate simplicial mesh which we refer to as a computational mesh, $\mathcal{T}_{c,h}=\{\xi_{1},...,\xi_{N_{v}}\}$. We assume that $\mathcal{T}_{h}$, $\hat{\mathcal{T}}_{c,h}$, and $\mathcal{T}_{c,h}$ have the same number of elements and vertices and the same connectivity. As a consequence, there exists a unique element $K_{c}\in\mathcal{T}_{c,h}$ corresponding to any element $K\in\mathcal{T}_{h}$. The affine mapping between $K_{c}$ and $K$ is denoted by $F_{K}$ and its Jacobian matrix by $F_{K}^{\prime}$. Let the coordinates of the vertices of $K$ and $K_{c}$ be $x_{0}^{K},\,x_{1}^{K},\,x_{2}^{K}$ and $\xi_{0}^{K},\,\xi_{1}^{K},\,\xi_{2}^{K}$, respectively. Then we have

[TABLE]

From this we get

[TABLE]

where $E_{K}$ and $\hat{E}_{K}$ are the edge matrices of $K$ and $K_{c}$ defined as

[TABLE]

The objective of the MMPDE moving mesh method is to generate an adaptive mesh as a uniform one in the metric $\mathbb{M}$. Such an $\mathbb{M}$-uniform mesh requires that (i) the ratio of the area of $K$ in the metric $\mathbb{M}$ to the area of $K_{c}$ in the Euclidean metric stay constant for all elements and (ii) $K$ measured in the metric $\mathbb{M}$ be similar to $K_{c}$ measured in the Euclidean metric for all elements. These two requirements can be expressed mathematically as the equidistribution and alignment conditions (e.g., see [23, 28]),

[TABLE]

where $|K|$ and $|K_{c}|$ denote the area of $K$ and $K_{c}$, respectively, $\mathbb{M}_{K}$ is the average of $\mathbb{M}$ over $K$, $\det(\cdot)$ and $\operatorname{tr}(\cdot)$ denote the determinant and trace of a matrix, and

[TABLE]

An energy functional based on these conditions has been proposed by Huang [21] as

[TABLE]

where $\mathbb{J}_{K}=(F^{\prime}_{K})^{-1}$ and

[TABLE]

Here, $0<\theta\leq\frac{1}{2}$ and $p>1$ are two dimensionless parameters. We use $\theta=\frac{1}{3}$ and $p=\frac{3}{2}$, which are known experimentally to work well for most problems.

Our goal is to find a new physical mesh $\mathcal{T}_{h}^{n+1}$ by minimizing $I_{h}$. We use an indirect approach with which we take $\mathcal{T}_{h}=\mathcal{T}_{h}^{n}$ and then minimize $I_{h}$ with respect to $\mathcal{T}_{c,h}$. The minimization is carried out by integrating the MMPDE (see (22) below) from $t^{n}$ to $t^{n+1}$ with the reference computational mesh $\hat{\mathcal{T}}_{c,h}$ as the initial mesh. The obtained computational mesh is denoted as $\mathcal{T}_{c,h}^{n+1}$. Notice that $\mathcal{T}_{h}^{n}$ is kept unchanged during the integration and $\mathcal{T}_{c,h}^{n+1}$ and $\mathcal{T}_{h}^{n}$ form a correspondence, i.e., $\mathcal{T}_{h}^{n}=\Phi_{h}(\mathcal{T}_{c,h}^{n+1})$. The new physical mesh $\mathcal{T}_{h}^{n+1}$ is then defined as

[TABLE]

which can be computed readily using linear interpolation. This procedure is explained in Fig. 3.

We now describe the MMPDE moving mesh method [26, 27] for minimizing $I_{h}$. The MMPDE is defined as a gradient system of $I_{h}$, i.e.,

[TABLE]

where $\tau>0$ is a parameter used to adjust the time scale of mesh movement and $P_{j}=\det(\mathbb{M}(x_{j}))^{\frac{p-1}{2}}$ is chosen such that (22) is invariant under the scaling transformation of $\mathbb{M}$. The derivative of $I_{h}$ with respect to $\xi_{j}$ is considered as a row vector and can be found analytically by the notion of scale-by-matrix differentiation; see [24]. Using the analytical formula, we can rewrite (22) as

[TABLE]

where $\omega_{j}$ is the element patch associated with the $j$-th vertex, $j_{K}$ is its local index of the vertex on $K$, and the local velocities ${v}_{j_{K}}^{K}$ are given by

[TABLE]

where the derivatives of the function $G$ are given by

[TABLE]

MMPDE (23), with proper modifications for the boundary vertices to allow them only to slide on the boundary, is integrated from $t^{n}$ to $t^{n+1}$ with the initial mesh $\hat{\mathcal{T}}_{c,h}$. Notice that the integration of the MMPDE is equivalent to performing steepest descent for minimizing $I_{h}$. In our computation, we use the Matlab®  function ode15s, a Numerical Differentiation Formula based integrator, for the integration.

4 Numerical results

In this section we present numerical results obtained with the moving mesh finite element method described in the previous section for three examples. The first two examples are benchmark problems commonly used in the existing literature to examine mathematical models for brittle fracture and related numerical algorithms. The last example is chosen to test the ability of our method to handle multiple and complex cracks. Special attention is paid to the demonstration of the effectiveness of the MMPDE method to track crack propagation and the effects of the regularization methods on the convergence of Newton’s iteration. In the results presented in this section, an adaptive mesh of size $N=6,400$ and the sonic-point regularization method with $\alpha=1\times 10^{-3}$ are used, unless stated otherwise.

4.1 Example 1. Single edge notched tension test

We first consider a single edge notched tension test from Miehe et al. [40], with the domain and boundary conditions shown in Fig. 4(a). For the boundary conditions, the bottom edge of the domain is fixed and the top edge is fixed along $x$-direction while a uniform $y$-displacement $U$ is increased with time to drive the crack propagation.

The solid is assumed to be homogeneous isotropic with elastic bulk modulus $\lambda=121.15$ kN/mm2 and shear modulus $\mu=80.77$ kN/mm2. The fracture toughness is $g_{c}=2.7\times 10^{-3}$ kN/mm. Two displacement increments have been prescribed for the computation, $\Delta U=1\times 10^{-5}$ mm is chosen for the first 500 time steps and $\Delta U=1\times 10^{-6}$ mm afterwards. We consider two values of $l$, $0.00375$ mm and $0.0075$ mm. An initial triangular mesh is constructed from a rectangular mesh by subdividing each rectangle into four triangles along the diagonal directions.

Typical adaptive meshes and contours of the phase-field variable during crack evolution for $l=0.00375$ mm and $0.0075$ mm are shown in Fig. 5 and Fig. 6, respectively. As can be seen, the mesh points concentrate in the region around the crack, which demonstrates the effectiveness of the mesh adaptation strategy for this tension test. Closer views around the crack and crack tip are shown in Fig. 7(b) and Fig. 7(c). It can be observed that the mesh stays symmetric near the crack for all time. For comparison purpose, we compute the surface load vector on the top edge as

[TABLE]

where $n$ is the unit outward normal to the top edge. Particularly, we are interested in $F_{y}$ for tension test and $F_{x}$ for shear test. The load-deflection curves are shown in Fig. 8.

4.1.1 Effects of the regularization methods

We now investigate the effects of the regularization methods (cf. Section 3.2) on the convergence of Newton’s iteration and load-deflection curves. For convergence, we consider the first displacement increment (where $\Delta U=1\times 10^{-5}$ mm is applied on the top edge) with $k_{l}=0$. The convergence history of Newton’s iteration is shown in Fig. 9. As can be seen, Newton’s iteration fails to converge without regularization on the decomposition of the strain tensor. On the other hand, convergence is reached using the sonic-point regularization with $\alpha\geq 1\times 10^{-4}$ and exponential convolution and smoothed 2-point convolution methods with $\alpha\geq 4\times 10^{-4}$. Moreover, Newton’s iteration converges faster for larger $\alpha$ for all of the three methods.

Next, we examine the effects of the regularization methods on load-deflection curves. We first compare the three regularization methods ($\alpha=1\times 10^{-3}$, $l=0.0075$ mm and $k_{l}=1\times 10^{-3}$). For the computation without regularization ($\alpha=0$), $k_{l}=1\times 10^{-3}$ is used to obtain convergent results. The results are shown in Fig. 10(a). As we can see, there is no significant difference between the regularization methods with $\alpha=1\times 10^{-3}$ and without regularization, except at the beginning of the load time.

We then investigate the effects of the regularization parameter. The load-deflection curves for the three regularization methods with various values of $\alpha$ ($4\times 10^{-4},1\times 10^{-3},5\times 10^{-3}$, and $1\times 10^{-2}$) are shown in Fig. 10(b), 10(c), 10(d), respectively. We can see that, for all three regularization methods, when $\alpha$ is small ($\alpha\leq 1\times 10^{-3}$), influence on the numerical solution is almost invisible. For large $\alpha$ values ($\alpha\geq 5\times 10^{-3}$), all three regularization methods underestimate the load before the crack starts propagating and overestimates afterwards. Moreover, the effects using the smoothed 2-point convolution method are less significant than those with the sonic-point method. The effects of the exponential convolution method lie between those of the smoothed 2-point convolution and sonic-point methods. Since the sonic-point method has the simplest form and is more economic to compute than the other two methods, and since all three methods give almost identical results when $\alpha$ is small, we use the sonic-point method with $\alpha=1\times 10^{-3}$ for later computations.

4.1.2 Effects of mesh adaptivity

We now study the effects of mesh adaptivity on the numerical solution. We take $l=0.0075$ mm in the computation and compare the results with uniform meshes. As mentioned before, the mesh should be sufficiently fine ($h\ll l$) to resolve the cracks. A uniform triangular mesh is obtained by first partitioning the domain into a uniform rectangular mesh and then dividing each sub-rectangle into four triangles using the diagonal lines and is represented by the number of sub-rectangles in $x-$ and $y$-directions. We consider different mesh sizes as $131\times 131$ ($N=67,600$), $151\times 151$ ($N=90,000$), $171\times 171$ ($N=115,600$), $211\times 211$ ($N=176,400$), $251\times 251$ ($N=250,000$), and $351\times 351$ ($N=490,000$). For adaptive meshes, we start with uniform meshes of size $31\times 31$ ($N=3,600$), $41\times 41$ ($N=6,400$), and $51\times 51$ ($N=10,000$). Fig. 11 shows the load-deflection curves for different meshes. It clearly demonstrates the spatial convergence in terms of the number of mesh elements for both uniform and adaptive meshes. Moreover, moving meshes use significantly fewer elements to achieve the same accuracy, as can be seen from Fig. 11(b) and Fig. 12.

To verify the accuracy of the numerical solutions, the load-deflection curves with a uniform mesh of $351\times 351$ ($N=490,000$) and an adaptive mesh of $41\times 41$ ($N=6,400$) are plotted in Fig. 12. They are comparable with each other and agree very well with the result obtained by Miehe et al. [40] where a mesh pre-refined around the regions of the crack and its expected propagation path is used. The computation (with a Matlab®  implementation of the method) takes about 1,339.7 seconds of CPU time (for one time step) on the machine with a single AMD Opteron 6386 SE $2.80$ GHz processor for the uniform mesh and about 141.07 seconds of CPU time for the adaptive mesh. It is clear that the moving mesh method improves computational efficiency significantly. The relative cost between moving mesh and solving $u$ and $d$ can be clearly seen in Table 1. The CPU time for mesh movement takes a large proportion ($73.58\%$) when the moving mesh method described in the previous section was used. This can further be improved by employing a two-level mesh movement strategy; e.g., see Huang [22].

4.2 Example 2. Single edge notched shear test

In this example, we consider a single edge notched shear test, with the domain and boundary conditions shown in Fig. 4(b). The bottom edge of the domain is fixed and the top edge is fixed along $y$-direction while a uniform $x$-displacement $U$ is increased with time to drive the crack propagation. The elastic bulk modulus is $\lambda=121.15$ kN/mm2, the shear modulus is $\mu=80.77$ kN/mm2, and the fracture toughness is $g_{c}=2.7\times 10^{-3}$ kN/mm. The displacement increment is chosen as $\Delta U=1\times 10^{-5}$ mm for the computation.

Typical adaptive meshes and contours of the phase-field variable during crack evolution for $l=0.00375$ mm and $0.0075$ mm are shown in Fig. 13 and Fig. 14, respectively. The mesh concentrates dynamically around the crack as it evolves under continuous load. Close views of the meshes around the crack and crack tip are plotted in Fig. 15. The mesh concentration is adequate especially during the turning process of the shear crack. This property is important for handling large fracture networks with complex crack propagation such as joining, branching and nonplanar propagation. Fig. 16 shows the load-deflection curves for the shear test with $l=0.00375$ mm and $0.0075$ mm using a moving mesh of $N=6,400$.

Next, we compare the regularization methods for the shear test problem. Fig. 17 shows the convergence of Newton’s iteration using the three regularization methods. As we can see, all of the methods make Newton’s iteration convergent for $\alpha\geq 4\times 10^{-4}$. The convergence improves as $\alpha$ increases. Moreover, Newton’s iteration converges even for $\alpha=1\times 10^{-4}$ when the sonic-point and exponential convolution methods are used.

Lastly, the load-deflection curves for the shear test with $l=0.0075$ mm and $k_{l}=1\times 10^{-3}$ are shown in Fig. 18(a) for the regularization methods and without regularization. We can see that the curves for the regularization methods with $\alpha=1\times 10^{-3}$ are almost the same as that without regularization. The load-deflection curves obtained using the three regularization methods with various values of $\alpha$ are presented in Fig. 18(b), 18(c) and 18(d), respectively. For the sonic-point regularization method, the load-deflection curve becomes unphysical when $\alpha\geq 1\times 10^{-2}$. It is also interesting to observe that the load for the exponential convolution and smoothed 2-point convolution regularization is overestimated after crack starts propagating for large values of $\alpha$ ($\alpha\geq 5\times 10^{-3}$). The effects of $\alpha$ are less significant with smoothed 2-point convolution than the other two methods.

4.3 Example 3. Test with multiple cracks

In this example, we first consider two cracks in a square plate of width $2$ mm, with the domain and boundary conditions shown in Fig. 19(a), to test the modeling of the junction between two cracks. Crack 1 is centered at ($-0.2,0$) with length $0.6$ mm and polar angle (relative to the horizontal direction) $9^{\circ}$, while Crack 2 is centered at ($0.46,0$) with length $0.8$ mm and polar angle $65^{\circ}$. The bottom edge of the domain is fixed and the top edge is fixed along $x$-direction while a uniform $y$-displacement $U$ is increased with time to drive the crack propagation. The material parameters are the same as previous examples, that is, elastic constants $\lambda=121.15$ kN/mm2, $\mu=80.77$ kN/mm2, and the fracture toughness is $g_{c}=2.7\times 10^{-3}$ kN/mm. The displacement increment is chosen as $\Delta U=1\times 10^{-4}$ mm for the computation. The mesh consists of $10,000$ triangular elements and the length scale parameter is chosen as $l=0.00375$ mm. A similar configuration with different boundary conditions has been used by Budyn et al. [15].

Typical adaptive meshes and contours of the phase-field variable during crack evolution are shown in Fig. 20. As can be seen, the mesh adapts dynamically to capture the junction process of the two cracks. As the load $U$ increases, the evolution can be described as follows.

• (a)

$U=8\times 10^{-3}$ mm: the tip of Crack 1 activates;

• (b)

$U=8\times 10^{-3}$ to $1.0\times 10^{-2}$ mm: Crack 1 propagates along a curved path heading towards Crack 2;

• (c)

$U=1.1\times 10^{-2}$ mm: Crack 1 connects to Crack 2, and one of the tip for Crack 2 activates;

• (d)

$U=1.6\times 10^{-2}$ mm: both Crack 1 and Crack 2 have propagated to the edge of the plate.

Since there is no analytical solution for this problem, we cannot compare the computed solutions with the exact solution. Nevertheless, they can be justified in physics. In fact, it is known (e.g., see [46]) that the critical stress increases as the polar angle increases and cracks with large polar angles are harder to initiate. Therefore, Crack 1, which is oriented with a smaller polar angle, is the first to initiate and propagate. Moreover, the frequency of crack coalescence is strongly related to the crack density (the crack size and their relative location) within the solid; e.g., see [55]. The higher the crack density is, the more likely two cracks will merge. Since being close to Crack 2, Crack 1 begins to merge with Crack 2 as it grows. At the same time, the other side of Crack 1 has propagated to the edge of the plate, which causes the left of the plate to lose strength. At later stages, the right side of the tip for Crack 2 activates and propagates to the edge of the plate.

Next, we consider five initial cracks in the same domain and with the same boundary conditions as for the two-crack problem, see Fig. 19(b). The material parameters are the same as previous examples except $g_{c}=2.7\times 10^{-4}$ kN/mm. The lengths of Crack 1, 2, 3, 4, and 5 are $0.3$ mm, $0.35$ mm, $0.35$ mm, $0.5$ mm, and $0.5$ mm, centered at ($-0.6,0.3$), ($0,0.5$), ($0.6,0.5$), ($-0.5,-0.4$), and ($0.5,-0.2$), with polar angle $30^{\circ}$, $45^{\circ}$, $17^{\circ}$, $28.6^{\circ}$, and $9^{\circ}$, respectively.

Typical adaptive meshes and contours of the phase-field variable during crack evolution are shown in Fig. 21. The evolution can be described as follows.

• (a)

$U=4\times 10^{-3}$ mm: the tips of Crack 3, 4 and 5 activate;

• (b)

$U=6\times 10^{-3}$ mm: Crack 4 connects to Crack 5, and Crack 3 has propagated to the right edge of the plate;

• (c)

$U=1\times 10^{-2}$ mm: Crack 3 joins Crack 2, and Crack 4 has propagated to the left edge of the plate;

• (d)

$U=1.7\times 10^{-2}$ mm: Crack 2 joins Crack 1.

Notice that the tips of Crack 3, 4, and 5 activate earlier due to their smaller polar angles. Compared with other cracks, Crack 4 and 5 have longer length and a closer distance, which results in the early crack merging. As the load increases, Crack 4 connects to Crack 5 and Crack 3 from both ends. At later stages, Crack 2 begins to propagate after Crack 3 joining Crack 2. Throughout the entire process, Crack 1 is limited to propagation due to its higher polar angle and the interaction of other cracks.

Finally, we consider a more complex situation with ten initial cracks and with the same domain and boundary conditions as for the two-crack problem; see Fig. 19(c). The material parameters are the same as in the previous examples except $g_{c}=2.7\times 10^{-4}$ kN/mm. All of Crack 1 through Crack 10 have the same length $0.1$ mm. They are centered at ($-0.5,0.8$), ($0.2,0.8$), ($-0.3,0.3$), ($0.5,0.5$), ($0,0$), ($-0.7,-0.2$), ($-0.5,-0.5$), ($-0.1,-0.8$), ($0.5,-0.75$), and ($0.7,-0.2$), with polar angle $40^{\circ}$, $45^{\circ}$, $109^{\circ}$, $132^{\circ}$, $143^{\circ}$, $40^{\circ}$, $45^{\circ}$, $120^{\circ}$, $40^{\circ}$, and $115^{\circ}$, respectively. Typical adaptive meshes and contours of the phase-field variable during crack evolution are shown in Fig. 22. The evolution can be described as follows.

• (a)

$U=6\times 10^{-3}$ mm: the tips of Crack 4, 5, 6, and 9 activate;

• (b)

$U=8\times 10^{-3}$ mm: Crack 5 connects to Crack 6, and Crack 5 has propagated to the left edge of the plate;

• (c)

$U=2.4\times 10^{-2}$ mm: Crack 4 joins Crack 3, and Crack 4 has propagated to the right edge of the plate.

The results for both the two-, five- and ten-crack problems also demonstrate that the MMPDE moving mesh method captures the crack propagation successfully for multiple and complex crack systems.

5 Conclusions

In the previous sections we have studied the moving mesh finite element solution of phase-field models for brittle fracture and investigated regularization methods to improve the convergence of Newton’s iteration for solving the nonlinear system resulting from finite element discretization. The MMPDE moving mesh method has been used to track the crack propagation and improve the computational efficiency. Numerical examples have been presented to demonstrate the effectiveness of the moving mesh method to dynamically concentrate mesh elements around propagating cracks and its ability to handle multiple and complex crack systems.

A distinguished feature in the numerical simulation of brittle fracture using phase-field modeling is the decomposition of the strain tensor in the elastic energy which is necessary to account for the reduction in the stiffness of an elastic solid due to the presence of cracks. The decomposition makes the energy functional non-smooth and increases the nonlinearity of the governing equation. Computationally, this causes the non-existence of the Jacobian matrix of the nonlinear discrete equation (15) at some places and the failure of Newton’s iteration to converge (cf. Figs. 9(a) and 17(a)). Three methods, the sonic-point, exponential convolution, and smoothed 2-point convolution, have been proposed in Section 3.2 to regularize the decomposition of the strain tensor. Numerical examples have demonstrated that all of these methods can effectively improve the convergence of Newton’s iteration for relatively small values of the regularization parameter $\alpha$ but without compromising the accuracy of numerical simulation. Generally speaking, the larger $\alpha$ is, the faster Newton’s iteration converges. The sonic-point method works with smaller $\alpha$ than the other two methods while the smoothed 2-point convolution method has slightly smaller effects on the load-deflection curves for the same $\alpha$ than the other two methods.

Another benefit of the regularization is that the finite-difference approximation of the Jacobian matrix, which is known to work for smooth problems, also works well in the current situation with regularization.

Acknowledgment. The authors would like to thank Dr. Junbo Cheng for bringing their attention to the sonic-point regularization method in computational fluid dynamics. F.Z. was supported by China Scholarship Council (CSC) and China University of Petroleum - Beijing (CUPB) for his research visit to the University of Kansas from September of 2015 to September of 2017. F.Z. is thankful to Department of Mathematics of the University of Kansas for the hospitality during his visit. The authors are grateful to the anonymous referees for their valuable comments in improving the quality of the paper.