Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies

TL;DR

This paper presents a straightforward and efficient implementation of the phase-field model for crack propagation in COMSOL Multiphysics, validated through various 2D and 3D case studies demonstrating its effectiveness.

Contribution

The authors develop a simple, coupled implementation of the phase-field fracture model in COMSOL, enabling effective simulation of quasi-static and dynamic cracks without complex coding.

Findings

Good agreement with previous models confirms validity.

Implementation effectively handles complex crack patterns.

Simulation results demonstrate computational efficiency.

Abstract

The phase-field model (PFM) represents the crack geometry in a diffusive way without introducing sharp discontinuities. This feature enables PFM to effectively model crack propagation compared with numerical methods based on discrete crack model, especially for complex crack patterns. Due to the involvement of \phased field", phase-field method can be essentially treated a multifield problem even for pure mechanical problem. Therefore, it is supposed that the implementation of PFM based on a software developer that especially supports the solution of multifield problems should be more effective, simpler and more efficient than PFM implemented on a general finite element software. In this work, the authors aim to devise a simple and efficient implementation of phase-field model for the modelling of quasi-static and dynamic fracture in the general purpose commercial software developer, COMSOL Multiphysics. Notably only the tensile stress induced crack is accounted for crack evolution by using the decomposition of elastic strain energy. The width of the diffusive crack is controlled by a length-scale parameter. Equations that govern body motion and phase-field evolution are written into different modules in COMSOL, which are then coupled to a whole system to be solved. A staggered scheme is adopted to solve the coupled system and each module is solved sequentially during one time step. A number of 2D and 3D examples are tested to investigate the performance of the present implementation. Our simulations show good agreement with previous works, indicating the feasibility and validity of the COMSOL implementation of PFM.

Figures

Equations

$${\Psi }_{opt}(\mathbf{u},\Gamma )={\int }_{\Omega }{\psi }_{\epsilon }(\mathbf{\epsilon })\mathrm{d}\Omega +{\int }_{\Gamma }{G}_c\mathrm{d}S-{\int }_{\Omega }\mathbf{b}\cdot \mathbf{u}\mathrm{d}\Omega -{\int }_{\partial {\Omega }_{h_i}}\mathbf{f}\cdot \mathbf{u}\mathrm{d}S$$

$${\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$

$${\psi }_{\epsilon }(\mathbf{\epsilon })=\frac{1}{2}\lambda {\epsilon }_{ii}{\epsilon }_{jj}+\mu {\epsilon }_{ij}{\epsilon }_{ij}$$

$$\gamma (\varphi ,▽\varphi )=\frac{{\varphi }^2}{2l_0}+\frac{l_0}{2}\frac{\partial \varphi }{\partial x_i}\frac{\partial \varphi }{\partial x_i}$$

$${\int }_{\Gamma }{G}_c\mathrm{d}S={\int }_{\Omega }{G}_c\left[\frac{{\varphi }^2}{2l_0}+\frac{l_0}{2}\frac{\partial \varphi }{\partial x_i}\frac{\partial \varphi }{\partial x_i}\right]\mathrm{d}\Omega$$

$$\{\begin{array}{r}{\mathbf{\epsilon }}_{+}=\sum _{a=1}^d⟨{\epsilon }_a{⟩}_{+}{\mathbf{n}}_a\otimes {\mathbf{n}}_a\\ {\mathbf{\epsilon }}_{-}=\sum _{a=1}^d⟨{\epsilon }_a{⟩}_{-}{\mathbf{n}}_a\otimes {\mathbf{n}}_a\end{array}$$

$$\{\begin{array}{r}{\psi }_{\epsilon }^{+}(\mathbf{\epsilon })=\frac{\lambda }{2}⟨tr(\mathbf{\epsilon }){⟩}_{+}^2+\mu tr\left({\mathbf{\epsilon }}_{+}^2\right)\\ {\psi }_{\epsilon }^{-}(\mathbf{\epsilon })=\frac{\lambda }{2}⟨tr(\mathbf{\epsilon }){⟩}_{-}^2+\mu tr\left({\mathbf{\epsilon }}_{-}^2\right)\end{array}$$

$${\psi }_{\epsilon }(\mathbf{\epsilon })=\left[(1-k)(1-\varphi {)}^2+k\right]{\psi }_{\epsilon }^{+}(\mathbf{\epsilon })+{\psi }_{\epsilon }^{-}(\mathbf{\epsilon })$$

$${\psi }_{kin}(\dot{\mathbf{u}})=\frac{1}{2}{\int }_{\Omega }\rho {\dot{u}}_i{\dot{u}}_i\mathrm{d}\Omega$$

$$\begin{gathered} L=\frac{1}{2}{\int }_{\Omega }\rho {\dot{u}}_i{\dot{u}}_i\mathrm{d}\Omega -{\int }_{\Omega }\left\{\left[(1-k)(1-\varphi {)}^2+k\right]{\psi }_{\epsilon }^{+}(\mathbf{\epsilon })+{\psi }_{\epsilon }^{-}(\mathbf{\epsilon })\right\}\mathrm{d}\Omega -{\int }_{\Omega }{G}_c\left[\frac{{\varphi }^2}{2l_0}+\frac{l_0}{2}\frac{\partial \varphi }{\partial x_i}\frac{\partial \varphi }{\partial x_i}\right]\mathrm{d}\Omega +{\int }_{\Omega }{b}_i{u}_i\mathrm{d}\Omega + \\ {\int }_{\partial {\Omega }_{h_i}}{f}_i{u}_i\mathrm{d}S \end{gathered}$$

$$\{\begin{array}{r}\frac{\partial {\sigma }_{ij}}{\partial x_j}+b_i=\rho {\ddot{u}}_i\\ \left[\frac{2l_0(1-k){\psi }_{\epsilon }^{+}}{G_c}+1\right]\varphi -l_0^2\frac{{\partial }^2\varphi }{\partial x_i^2}=\frac{2l_0(1-k){\psi }_{\epsilon }^{+}}{G_c}\end{array}$$

$${\sigma }_{ij}=\left[(1-k)(1-\varphi {)}^2+k\right]\frac{\partial {\psi }_{\epsilon }^{+}}{\partial {\epsilon }_{ij}}+\frac{\partial {\psi }_{\epsilon }^{-}}{\partial {\epsilon }_{ij}},$$

$$\mathbf{\sigma }=\left[(1-k)(1-\varphi {)}^2+k\right]\left[\lambda ⟨tr(\mathbf{\epsilon }){⟩}_{+}\mathbf{I}+2\mu {\mathbf{\epsilon }}_{+}\right]+\lambda ⟨tr(\mathbf{\epsilon }){⟩}_{-}\mathbf{I}+2\mu {\mathbf{\epsilon }}_{-}$$

$$H(\mathbf{x},t)=\underset{s\in [0,t]}{\max }{\psi }_{\epsilon }^{+}\left(\mathbf{\epsilon }(\mathbf{x},s)\right)$$

$$\{\begin{array}{r}\frac{\partial {\sigma }_{ij}}{\partial x_j}+b_i=\rho {\ddot{u}}_i\\ \left[\frac{2l_0(1-k)H}{G_c}+1\right]\varphi -l_0^2\frac{{\partial }^2\varphi }{\partial x_i^2}=\frac{2l_0(1-k)H}{G_c}\end{array}$$

$$\{\begin{array}{rl}& {\sigma }_{ij}{m}_j=f_i\mathrm{on}\partial {\Omega }_{h_i}\\ & \frac{\partial \varphi }{\partial x_i}{m}_i=0\mathrm{on}\partial \Omega \end{array}$$

$$\{\begin{array}{rlr}& \mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0(\mathbf{x})& \mathbf{x}\in \Omega \\ & \mathbf{v}(\mathbf{x},0)={\mathbf{v}}_0(\mathbf{x})& \mathbf{x}\in \Omega \\ & \varphi (\mathbf{x},0)={\varphi }_0(\mathbf{x})& \mathbf{x}\in \Omega \end{array}.$$

$${\sigma }_{cr}=\frac{9}{16}\sqrt{\frac{E{G}_c}{3l_0}}$$

$$l_0=\frac{27E{G}_c}{256{\sigma }_{cr}^2}$$

$$\mathbf{D}=\frac{\partial \mathbf{\sigma }}{\partial \mathbf{\epsilon }}=\lambda \left\{\left[(1-k)(1-\varphi {)}^2+k\right]H_{\epsilon }(tr(\mathbf{\epsilon }))+H_{\epsilon }(-tr(\mathbf{\epsilon }))\right\}\mathbf{J}+2\mu \left\{\left[(1-k)(1-\varphi {)}^2+k\right]\frac{\partial {\mathbf{\epsilon }}_{+}}{\partial \mathbf{\epsilon }}+\frac{\partial {\mathbf{\epsilon }}_{-}}{\partial \mathbf{\epsilon }}\right\}$$

$$\bm{D}_{e}=\left[\begin{array}[]{cccccc}D_{1111}&D_{1122}&D_{1133}&D_{1112}&D_{1123}&D_{1113}\\ D_{2211}&D_{2222}&D_{2233}&D_{2212}&D_{2223}&D_{2213}\\ D_{3311}&D_{3322}&D_{3333}&D_{3312}&D_{3323}&D_{3313}\\ D_{1211}&D_{1222}&D_{1233}&D_{1212}&D_{1223}&D_{1213}\\ D_{2311}&D_{2322}&D_{2333}&D_{2312}&D_{2323}&D_{2313}\\ D_{1311}&D_{1322}&D_{1333}&D_{1312}&D_{1323}&D_{1313}\end{array}\right]$$

$${\bar{D}}_{ijkl}=\lambda \left\{\left[(1-k)(1-\varphi {)}^2+k\right]H_{\epsilon }(tr(\mathbf{\epsilon }))+H_{\epsilon }(-tr(\mathbf{\epsilon }))\right\}{\delta }_{ij}{\delta }_{kl}$$

$${\tilde{D}}_{ijkl}=2\mu \left\{\left[(1-k)(1-\varphi {)}^2+k\right]P_{ijkl}^{+}+P_{ijkl}^{-}\right\}$$

$$P_{ijkl}^{+}=\sum _{a=1}^3\sum _{b=1}^3{H}_{\epsilon }({\epsilon }_a){\delta }_{ab}{n}_{ai}{n}_{aj}{n}_{bk}{n}_{bl}+\sum _{a=1}^3\sum _{b\ne a}^3\frac{1}{2}\frac{⟨{\epsilon }_a{⟩}_{+}-⟨{\epsilon }_b{⟩}_{+}}{{\epsilon }_a-{\epsilon }_b}{n}_{ai}{n}_{bj}(n_{ak}{n}_{bl}+n_{bk}{n}_{al})$$

$$P_{ijkl}^{-}=\sum _{a=1}^3\sum _{b=1}^3{H}_{\epsilon }(-{\epsilon }_a){\delta }_{ab}{n}_{ai}{n}_{aj}{n}_{bk}{n}_{bl}+\sum _{a=1}^3\sum _{b\ne a}^3\frac{1}{2}\frac{⟨{\epsilon }_a{⟩}_{-}-⟨{\epsilon }_b{⟩}_{-}}{{\epsilon }_a-{\epsilon }_b}{n}_{ai}{n}_{bj}(n_{ak}{n}_{bl}+n_{bk}{n}_{al})$$

$$\{\begin{array}{rlr}& {\epsilon }_1={\epsilon }_1(1+\delta )& \mathrm{if}{\epsilon }_1={\epsilon }_2\\ & {\epsilon }_3={\epsilon }_3(1-\delta )& \mathrm{if}{\epsilon }_2={\epsilon }_3\end{array}$$

$$\dot{H}=\{\begin{array}{rlr}& {\dot{\psi }}_{\epsilon }^{+},& {\psi }_{\epsilon }^{+}>0\mathrm{and}H={\psi }_{\epsilon }^{+}\\ & 0,& \mathrm{else}\end{array}$$

$$H_0(\mathbf{x})=\{\begin{array}{rlr}& \frac{B{G}_c}{2l_0}\left(1-\frac{2d(\mathbf{x},l)}{l_0}\right),& d(\mathbf{x},l)\le \frac{l_0}{2}\\ & 0,& d(\mathbf{x},l)>\frac{l_0}{2}\end{array}$$

$$B=\frac{\varphi }{(1-k)(1-\varphi )}$$

$${\int }_{\Omega }\left(-\rho \ddot{\mathbf{u}}\cdot \delta \mathbf{u}-\mathbf{\sigma }:\delta \mathbf{\epsilon }\right)\mathrm{d}\Omega +{\int }_{\Omega }\mathbf{b}\cdot \delta \mathbf{u}\mathrm{d}\Omega +{\int }_{{\Omega }_{h_i}}\mathbf{f}\cdot \delta \mathbf{u}\mathrm{d}S=0$$

Full text

1 Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Vietnam

2 Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Vietnam

3 Institute of Structural Mechanics, Bauhaus-University Weimar, Weimar 99423, Germany

4 Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, P.R. China

5 Institute of Continuum Mechanics, Leibniz University Hannover, Hannover 30167, Germany.

• Corresponding author: timon.rabczuk@tdt.edu.vn

Keywords: phase-field, multi-field, fracture mechanics, COMSOL, crack propagation, Quasi-static, dynamic fracture

1 Introduction

Fracture induced failure has obtained extensive concern in engineering designs because of the potential serious risks for structures and machines being used [1]. The research on crack initiation and propagation in solids has therefore become very important [2]. Particularly, when experiments are difficult, or even impossible to perform for studying certain type of crack propagation, researchers have to employ numerical approaches to predict complicated crack paths [3] such as those in multiple scales [4, 5, 6, 7, 8]. Consequently, a great number of numerical methods have been proposed to deal with crack problems in recent years.

Most of these methods have to describe complex crack geometry in the discrete setting, such as the discrete crack models [9], the extended finite element method (XFEM) [10, 11], generalized finite-elements method (GFEM) [12], and the phantom-node method [13, 14]. These methods all enrich the displacement field with discontinuities. Particularly, the discrete crack model [9] introduces new boundaries for the freshly created crack surfaces by an adaptive reconstruction of the mesh. XFEM [10] enriches the cracked elements by adding a set of discontinuous shape functions to the standard parts of FEM. Another common option to model cracks is the so-called cohesive elements [15, 16, 17] that allow displacement jumps on element boundaries and cracks are therefore restricted to penetrate along the corresponding element edges. In addition, the element-erosion methods [18, 19, 20] also succeeds in dealing with the fracture surfaces by setting the stresses of the elements, which meet the fracture criterion, as zero. However, the element-erosion methods have the disadvantage that they cannot simulate crack branching correctly [21]. Therefore, the complicated and special treatments for complex crack topologies have made these numerical approaches not so easy to implement and apply in practical engineering.

A recently emerged and developed approach, the phase-field method (PFM) [22, 23, 24, 25, 26], has attracted a lot of attention because of its relatively easier numerical implementation for fracture. The phase-field models utilize a scalar field (so-called phase-field) to represent the discrete cracks. The phase-field does not describe the crack as a physical discontinuity and just smoothly transits the intact material to the thoroughly broken one. The shape and propagation of the crack depend on the evolution equations of the phase-field. Thus, implementation of the phase-field does not require additional work to track the fracture surfaces algorithmically [24]. This results in that the phase-field methods have a large advantage over the discrete fracture models for modeling multiple and crack branching and merging in materials with arbitrary 2D and 3D geometries.

The phase-field models for quasi-static brittle crack started from Bourdin et al. [27] and improved by Miehe et al. [22, 23]. All these models are regarded as extension of the classical Griffith fracture theory and then extended to dynamic problems by Borden et al. [24]. In addition, Landau-Ginzburg type evolution equations [28] instead of the Griffith type have also been proposed and developed for the phase-field description of dynamic fracture. The progress in the phase-field models for quasi-static and dynamic crack problems has made PFM successfully applied in different problems, such as cohesive fractures [29], ductile fractures [30, 31], large strain problems [25], hydraulic fracturing [32], thermo-elastic problems [33, 34], electrochemical problems [35], thin shell [36], and stressed grain growth in polycrystalline metals [37, 38, 39]. These attempts imply that the application of the phase-field methods is quite beyond purely mechanical problems. This naturally requires a much easier implementation approach for the phase-field models. Otherwise, extensive application of the phase-field models will be restricted, especially in multi-physical problems.

Due to the smooth characteristics of the phase-field, the phase-field method can be implemented in any existing standard finite element to model complex crack patterns as shown in [22, 23]. Therefore, to reduce the efforts in implementation, it is desirable to implement phase-field method to an extensively used FEM code or commercial software. In fact, Msekh et al. [40] and [2] have implemented the phase-field method for brittle cracks in Abaqus. However, the phase-field modeling itself is essentially a multi-field problem even in the case of pure mechanical problem [22, 23]. From the authors’ experience, it is laborious and time consuming to implement a multifield problem in Abaqus. Therefore, a general purpose programme developer that especially supports the programming of multifield problem such as COMSOL has the potential to become a better solution than Abaqus.

In this paper, the possibility of simple and fast implementation of phase-field method is exploited for fracture modelling in a multifield programme developer, namely COMSOL Multiphysics. The phase-field modeling in COMSOL can be easily extended to problems that have more coupled fields by just adding suitable modules and coupling terms. It will be quite easy for readers to use this first-step implementation and augment it by other physical phenomena to solve multiphysics problems involving crack propagation. For example, the phase field implementation in COMSOL can be extended and applied to hydraulic fracturing, or compressed air energy storage [41, 42], which involves fluid pressure field, temperature, and cyclic effects [43, 44, 45, 46]. In this work, one phase-field model presented by Miehe et al. [22, 23] for a quasi-static crack problem and another one presented by Borden et al. [24] for dynamic problems are implemented in COMSOL in a staggered scheme. The elastic strain energy density is decomposed into two individual parts resulting from compression and tension, respectively. Thus, the fractures only due to tension can be obtained. In COMSOL, we use an implicit time integration scheme to enable the simulation. We also calculate some 2D and 3D benchmarks for quasi-static and dynamic crack propagation to show the feasibility of our approach for modeling fracture.

The paper is organized as follows. We begin with a short introduction of the phase-field model for brittle fractures based on the variational approach in Section 2. Subsequently, we present the numerical implementation of the phase-field model in COMSOL in Section 3. In Section 4, we examine some 2D and 3D numerical examples for cracks under quasi-static and dynamic loading. Finally, we end with conclusions regarding our findings in Section 5.

2 Phase-field model for fracture

2.1 Theory of brittle fracture

Let us consider an arbitrary body $\Omega\subset\mathbb{R}^{d}$ ($d\in\{1,2,3\}$) as shown in Fig. 1. The body $\Omega$ has an external boundary $\partial\Omega$ and internal discontinuity boundary $\Gamma$. The displacement of body $\Omega$ at time $t$ is denoted by $\bm{u}(\bm{x},t)\subset\mathbb{R}^{d}$ where $\bm{x}$ is the position vector. The displacement field satisfies the time-dependent Dirichlet boundary conditions, $u_{i}(\bm{x},t)=g_{i}(\bm{x},t)$, on $\partial\Omega_{g_{i}}\in\Omega$, and also the time-dependent Neumann conditions on $\partial\Omega_{h_{i}}\in\Omega$. We also consider a body force $\bm{b}(\bm{x},t)\subset\mathbb{R}^{d}$ acted on the body $\Omega$ and a traction $\bm{f}(\bm{x},t)$ on the boundary $\partial\Omega_{h_{i}}$.

A variational approach for fracture problems according to Griffith’s theory has been proposed in [47]. It states that the required energy to create a fracture surface per unit area is equal to the critical fracture energy density $G_{c}$, which is also commonly referred to as the critical energy release rate. The total potential energy $\Psi_{opt}(\bm{u},\Gamma)$ can be expressed in terms of the elastic energy $\psi_{\varepsilon}(\bm{\varepsilon})$, fracture energy and energy due to external forces:

[TABLE]

with the linear strain tensor $\bm{\varepsilon}=\bm{\varepsilon}(\bm{u})$ given by

[TABLE]

Isotropic linear elasticity is assumed and the elastic energy density $\psi_{\varepsilon}(\bm{\varepsilon})$ is given by [22]

[TABLE]

where $\lambda$ and $\mu$ are Lamé constants.

In addition, the variational approach [47] states that initiation, propagation and branching of the crack $\Gamma(\bm{x},t)$ at the time $t\in(0,T)$ for a point $x\in\Omega$ occur when the potential reaches the minimum value and the irreversible condition $\Gamma(\bm{x},s)\in\Gamma(\bm{x},t)(s<t)$ is satisfied. The irreversible condition means that a crack cannot be recovered to the uncracked state after its formation.

2.2 Phase filed approximation for fracture energy

The variational approach for brittle fracture was successfully implemented by Miehe et al. [22, 23], Borden et al. [24] and Bourdin et al. [27] with the introduction of a scalar phase filed. In this paper, we define a phase-field $\phi(\bm{x},t)\in[0,1]$ to approximate the fracture surface, $\Gamma$ (see also Fig. 1). $\phi(\bm{x},t)\in[0,1]$ represents naturally a diffusive crack topology. $\phi=1$ represents the crack and $\phi=0$ means that the body is uncracked. Thus, the crack surface density per unit volume of the solid is given by [22],

[TABLE]

where $l_{0}$ is a parameter that controls the transition region of the phase-field from 0 to 1. We call $l_{0}$ the length scale parameter that reflects the shape of a crack. The crack region will have a larger width with a larger $l_{0}$ and vice versa, see Fig. 1.

Applying Eq. (4), the fracture energy is approximated by

[TABLE]

The crack surface energy is transformed from the elastic energy as shown in [47], indicating that the elastic energy drives the evolution of the phase-field. In order to ensure that the crack is only driven by tensile load, it is important to decompose the elastic energy into tensile and compressive parts [27]. Here, the decomposition approach in Miehe et al. [22] is adopted to ensure the evolution of the phase-field will only be induced by the tensile part of the elastic energy density while compressive stress will not contribute to the propagation of crack. Therefore, the strain tensor $\bm{\varepsilon}$ is decomposed as follows

[TABLE]

where $\bm{\varepsilon}_{+}$ and $\bm{\varepsilon}_{-}$ are the tensile and compressive strain tensors, respectively. $\varepsilon_{a}$ is the principal strain and $\bm{n}_{a}$ is its direction vector. The operators $\langle\centerdot\rangle_{+}$ and $\langle\centerdot\rangle_{-}$ are defined as [22]: $\langle\centerdot\rangle_{+}=(\centerdot+|\centerdot|)/2$, $\langle\centerdot\rangle_{-}=(\centerdot-|\centerdot|)/2$.

By applying the decomposed strain tensor, the elastic energy density is represented as follows:

[TABLE]

It is assumed that only the tensile part of the elastic energy density is affected by the phase-field and then the following equation is used to model the stiffness reduction [24]:

[TABLE]

where $k$ is a model parameter that prevents the positive part of the elastic energy density from disappearing and the numerical singularity when phase-field $\phi$ tends to 1. In addition, it is required that $k>0$ and $k\ll 1$.

2.3 Governing equations

We consider also the kinetic energy of body $\Omega$:

[TABLE]

with $\bm{\dot{u}}=\frac{\partial{\bm{u}}}{\partial t}$ and $\rho$ being the density of body $\Omega$.

The total Lagrange energy functional can be expressed by the sum of the phase-field approximation for the fracture energy (5), the elastic energy (8), the kinetic energy (9) and the external potential energy by the external loads:

[TABLE]

The variation of the functional $L$ can be derived and its first variation should be zero, which leads to the following governing equations:

[TABLE]

where $\ddot{u}_{i}=\frac{\partial^{2}u}{\partial t^{2}}$ and $\sigma_{ij}$ is the component of Cauchy stress tensor given by

[TABLE]

and it can be rewritten as

[TABLE]

where $\bm{I}$ is a unit tensor $\in\mathbb{R}^{d\times d}$.

In order to ensure a monotonically increasing phase-field, the irreversibility condition is needed during compression or unloading. Thus, we introduce a strain-history field $H(\bm{x},t)$ [22, 23] defined by

[TABLE]

Replacing $\psi_{\varepsilon}^{+}$ by $H(\bm{x},t)$ in Eq. (11), the strong form is obtained by

[TABLE]

In addition, the zero first variation of the functional $L$ also achieves the Neumann conditions,

[TABLE]

with $m_{j}$ the component of the outward-pointing normal vector of the boundary.

For dynamic problems, the following initial conditions must be fulfilled,

[TABLE]

Here, the initial phase-field $\phi_{0}(\bm{x})$ can be used to model pre-existing cracks or geometrical features by setting it locally equal to 1 [24].

2.4 Estimation of $l_{0}$

Selecting the length scale parameter $l_{0}$ remains an open topic in phase field models for fracture. An analytical solution for the critical tensile stress $\sigma_{cr}$ that a 1D bar can sustain under tension has been derived by Borden et al. [24]:

[TABLE]

where $E$ is the Young’s modulus and $G_{c}$ the critical energy release rate. There is an apparent singularity when $l_{0}$ tends to zero, i.e. in case of a sharp crack which leads to a phyiscally meaningless infinite tensile strength. However, when all other parameters except $l_{0}$ are known, Eq. (18) can be solved for $l_{0}$ yielding

[TABLE]

In Eq. (19), the critical energy release rate $G_{c}$, Young’s modulus $E$, and critical stress $\sigma_{cr}$ can be estimated through some regular tests. Though the extension of this approach into higher-order dimensions and complex mixed-mode fracture is difficult, it gives at least some estimate how to choose $l_{0}$. It should be noted here again that no external fracture criterion is needed in the phase field method. The crack path is obtained by the evolution equation of phase field.

3 Implementation method in COMSOL

3.1 Overall framework

The phase-field modeling in this work is naturally a two-field problem ($\bm{u}$ and $\phi$). The phase-field model is implemented into the software COMSOL Multiphysics. We establish three main modules namely, Solid Mechanics Module, History-strain Module and phase-field Module. These modules are used to solve the three fields, $\bm{u}$, $H$ and $\phi$, respectively. These modules are all written in strong forms and solved based on standard finite element discretization in space domain and finite difference discretization in time domain. We also establish a preset Storage Module to calculate and store the internal field variables during each time step, such as the principal strains and their corresponding direction vectors.

3.2 Module setup

The Solid Mechanics Module is set up based on a linear elastic material library. The boundary and initial conditions in the Solid Mechanics Module are implemented as shown in Section 2. However, the elasticity matrix during each time step requires modification. The elasticity matrix is calculated on basis of the elasticity tensor of four order $\bm{D}$:

[TABLE]

where $H_{\varepsilon}\langle x\rangle$ is a Heaviside function with $H_{\varepsilon}\langle x\rangle=1$ for $x>0$ and $H_{\varepsilon}\langle x\rangle=0$ for $x\leq 0$ and $J_{ijkl}=\delta_{ij}\delta_{kl}$ with $\delta_{ij}$ and $\delta_{kl}$ are Kronecker deltas. Finally, the elasticity matrix $\bm{D}_{e}$ is rewritten as following:

[TABLE]

with $D_{ijkl}=\bar{D}_{ijkl}+\tilde{D}_{ijkl}$.

The component $\bar{D}_{ijkl}$ is expressed as

[TABLE]

According to the algorithm for fourth-order isotropic tensor [48], the component $\tilde{D}_{ijkl}$ is calculated based on the following:

[TABLE]

with

[TABLE]

and

[TABLE]

in which $n_{ai}$ denotes the $i$-th component of vector $\bm{n}_{a}$.

It can be seen from Eqs. (24) and (25) that Eq. (23) cannot be evaluated if $\varepsilon_{a}=\varepsilon_{b}$. We therefore adopt a “perturbation” technology for the principal strains [49] and make a change for better application in COMSOL:

[TABLE]

with the perturbation $\delta=1\times 10^{-9}$ for this paper. The second principal strain remains unchanged.

The Solid Mechanics Module has an inertial term in the governing equation as shown in Eq. (15). Thus, the governing equation is automatically suitable for a dynamic crack problem. For a quasi-static problem, the inertia term vanishes.

The phase-field Module is established by revising a module governed by the Helmholtz equation. The governing equation in Eq. (15), boundary condition in Eq. (16) and initial condition (17) are implemented in this module. The History-strain Module is set up based on the Distributed ODEs and DAEs Interfaces in COMSOL, which provide the possibility to solve distributed ODEs and DAEs in domains. The history strain field $H(\bm{x},t)$ is obtained by solving the following equation:

[TABLE]

Initial conditions are also required for the History-strain Module. Commonly $H_{0}(\bm{x})=0$ unless pre-existing cracks are modelled as the induced ones by the following expression [24]:

[TABLE]

where B is a scalar that controls the magnitude of the induced history field.

Letting $d=0$ and substituting $H_{0}$ into the second equation of (15), one will get:

[TABLE]

$B$ will become quite large if $\phi$ is close to 1, the value of the phase-field for the initial crack. Here, we chose $B=1\times 10^{6}$ for the simulation in this paper.

3.3 Finite element method and discretization

In COMSOL, the finite element method is used with the weak form of the governing equations given by

[TABLE]

and

[TABLE]

The standard vector-matrix notation is used with $\bm{u}_{i}$ and $\phi_{i}$ being the nodal values of the displacement and phase field. Then, we let the discretization as

[TABLE]

where $\bm{d}$ and $\hat{\bm{\phi}}$ are the vectors consisting of node values $\bm{u}_{i}$ and $\phi_{i}$. $\bm{N}_{u}$ and $\bm{N}_{\phi}$ are shape function matrices:

[TABLE]

where $n$ is the node number in one element and $N_{i}$ is the shape function at node $i$. Assuming that the test functions have the same discretization, we obtain

[TABLE]

where $\delta\bm{d}$ and $\delta\hat{\bm{\phi}}$ are the vectors consisting of node values of the test functions.

The gradients are thereby as follows

[TABLE]

where $\bm{B}_{u}$ and $\bm{B}_{\phi}$ are the derivatives of the shape functions:

[TABLE]

The equations of weak form (30) and (31) are then written as

[TABLE]

[TABLE]

For admissible arbitrary test functions, Eqs. (37) and (38) produces the discretized weak form as

[TABLE]

[TABLE]

where $\bm{F}_{u}^{ine}$, $\bm{F}_{u}^{int}$, and $\bm{F}_{u}^{ext}$ are the inertial, internal, and external forces for the displacement field and $\bm{F}_{\phi}^{int}$ and $\bm{F}_{\phi}^{ext}$ are the internal and external force terms of the phase field. Additionally, the mass and stiffness matrices follow

[TABLE]

3.4 Staggered method

The relationship between all the modules established is shown in Fig. 2. The “Storage Module” stores the results obtained from the “Solid Mechanics Module”, such as the magnitude of principal strains, direction of principal strain and elastic energy. The positive part of the elastic energy is then calculated and imported into the “History-strain Module” to solve and update the local history strain field. Then the updated history strain is used to solve the phase-field. The updated phase-field solution and the previously stored principal strains and their corresponding directions are used to modify the stiffness in the “Solid Mechanics Module” and then update the solution for the mechanical field. Fig. 2 shows the coupling for the solution of each module. Based on this, we employ a staggered scheme to solve the coupled system of equations as indicated in Fig. 3. Thus, the Newton-Raphson approach is adopted to obtain the residual of the discrete equations $\bm{R}_{u}=\bm{F}_{u}^{ext}-\bm{F}_{u}^{ine}-\bm{F}_{u}^{int}=0$ and $\bm{R}_{\phi}=\bm{F}_{\phi}^{ext}-\bm{F}_{\phi}^{int}=0$, respectively.

For the staggered time integration scheme, the equations of displacement, history strain and phase-field are solved independently. To obtain unconditional stability for the calculation, we use the implicit Generalized-$\alpha$ method [24]. When the time comes to a new value $t_{i}$, a new guess for the three field variables ($\bm{u}_{i}^{j=0}$, $H_{i}^{j=0}$ and $\phi_{i}^{j=0}$) is made first based on the results that have been solved in previous time steps. That is, linear extrapolation of the previous solution is used to construct the initial guess for the nonlinear system of equations to be solved at the present time step. For the given time step $i$ and iteration step $j+1$, the displacement $\bm{u}_{i}^{j+1}$ is first solved based on one Newton-Raphson iteration by using the results ($\bm{u}_{i}^{j}$, $H_{i}^{j}$ and $\phi_{i}^{j}$) of previous iteration step $j$. Using the updated displacements $\bm{u}_{i}^{j+1}$, the equation concerning history strain is then solved based on another Newton-Raphson iteration. Subsequently, the phase-field $\phi_{i}^{j+1}$ is obtained by the updated $\bm{u}_{i}^{j+1}$, $H_{i}^{j+1}$ and also a Newton-Raphson iteration. We finally compare the total relative error between the solution in previous and present iteration steps. If the error is less than the tolerance $\varepsilon_{t}$, the calculation is finished for current time step and will switch to the next step. Otherwise, the calculation will go through another iteration process until the tolerance requirement is satisfied. We choose the tolerance $\varepsilon_{t}=1\times 10^{-6}$ for our simulation. Thus, we succeed in obtaining all the solutions in the whole time domain by the implicit staggered time integration scheme.

It should be noted here that the iteration is slow to converge and more iteration steps are required when the material starts to fracture. Therefore, Anderson acceleration, a nonlinear convergence acceleration method that uses information from previous Newton iterations, is used to accelerate convergence [50]. The dimension of iteration space field is chosen as more than 50 to control the number of iteration increments in our work. In addition, we take standard Lagrangian elements (see the examples in the following section) to discretize the space domain for the three physical fields. Finally, Fig. 4 gives the flow chart of our implementation of phase-field method for crack problems in COMSOL. Our original codes can be downloaded from ”https://sourceforge.net/projects/phasefieldmodelingcomsol/”.

4 Numerical examples

In this section, we present several quasi-static and dynamic benchmark problems testing the influence of the length-scale parameter $l_{0}$, the element size, load and time step sizes as well as the critical energy release rate $G_{c}$ on the numerical results.

4.1 2D notched square plate subjected to tension

Consider a square plate with an initial notch subjected to static tension loading. This benchmark test has been calculated and analyzed by Miehe et al. [22, 23] and Hesch and Weinberg [25]. The geometry and loading condition are shown in Fig. 5. A vertical displacement $u_{y}$ is applied on the upper boundary of the plate with $u_{x}=0$. The material parameters are: $E$ = 210 GPa, $\nu$ = 0.3, and $G_{c}$ = 2700 J/m2. We choose $k=1\times 10^{-9}$ to avoid a singular stiffness matrix. The length parameter $l_{0}=7.5\times 10^{-3}$ mm and $1.5\times 10^{-2}$ mm, respectively. Plane strain conditions are assumed. The domain is discretized with 64516 Q4 (4 node quadrilateral) elements (with bi-linear shape functions). The element size $h$ is around $3.96\times 10^{-3}$ mm yielding $l_{0}=2h$ and $l_{0}=4h$ for $l_{0}=7.5\times 10^{-3}$ mm and $1.5\times 10^{-2}$ mm, respectively.

We apply an displacement increment of $\Delta u=1\times 10^{-5}$ mm for the first 450 time steps. Then, a displacement increment of $\Delta u=1\times 10^{-6}$ mm is chosen for the remaining time steps. We obtain the crack patterns at different displacements for the two fixed length scale parameters $l_{0}$, as shown in Fig. 6. As expected, the crack is less diffused for a smaller length scale parameter. The presented crack patterns are the same as those reported by Miehe et al. [22, 23], Hesch and Weinberg [25], Liu et al. [2]. The cracks propagate in horizontal direction. The load-displacement curves on the top boundary of the plate are shown in Fig. 7 in comparison with the results by Hesch and Weinberg [25]. The loads obtained by this work are in good agreement with those by Hesch and Weinberg [25] with the increase in the vertical displacement. A minor difference exists due to the different algorithm used in both methods. For a total of 453,520 degrees of freedom, COMSOL required 4 h 20 min ($l_{0}=1.5\times 10^{-2}$ mm) and 4 h 11 min ($l_{0}=7.5\times 10^{-3}$ mm) on two I5-6200U CPUs.

We also perform the simulation by changing the mesh size $h$ to $7.92\times 10^{-3}$ and $1.98\times 10^{-3}$ mm and the displacement increment $\Delta u$ to $2\times 10^{-6}$ and $5\times 10^{-7}$ mm in the remaining time steps. The results show that mesh size and displacement increment have no effect on the crack pattern. Figure 8 represents the load-displacement curves for different mesh sizes and displacement increments. A larger mesh size leads to a larger peak load. As Miehe et al. [22] have suggested a mesh size $h\leq 0.5l_{0}$ to obtain a precise crack topology, the mesh size $h=7.92\times 10^{-3}$ mm achieves a much larger peak load than the other mesh sizes for $l_{0}=7.5\times 10^{-3}$ mm. In addition, a much steeper post-peak stage can be seen for a smaller displacement increment.

4.2 3D notched square plate subject to tension

We now extend the benchmark problem in section 4.1 to 3D. The geometry of the plate in $x-y$ plane is the same as that in Section 4.1 with a thickness of 0.05 mm in $z$ direction. The material parameters are identical to the 2D example and we present results for a length scale parameter of $l_{0}=1.5\times 10^{-2}$, $1.25\times 10^{-2}$, $1.0\times 10^{-2}$ and $7.5\times 10^{-3}$ mm. All boundaries of the plate are fixed in the normal direction except the top boundary, which is subjected to a displacement of $u$ in the $y$ direction and is fixed in the $x$ and $z$ direction ($u_{x}=0$ and $u_{z}=0$).

The plate is discretized with 8-node Lagrangian elements of the same size $h=7.5\times 10^{-3}$ mm without special refinement in the expected path for crack propagation. For the staggered scheme, we apply the displacement increment $\Delta u=1\times 10^{-5}$ mm for the first 400 time steps and then adopt the displacement increment as $\Delta u=1\times 10^{-6}$ mm for the remaining time steps. Our simulation show that different length scale parameters $l_{0}$ achieve the same crack pattern. Crack patterns for $\phi>0.95$ in the 3D simulation with $l_{0}=1.5\times 10^{-2}$ and $7.5\times 10^{-3}$ mm are shown in Fig. 9. The load-displacement curve for the top boundary of the plate is shown in Fig. 10. As observed, the crack patterns and load-displacement curve are quite similar to those of the 2D case.The peak load of the plate increases as the length scale parameter decreases.

4.3 2D notched square plate subjected to shear loading

We now test the 2D notched square plate from Section 4.1 under shear loading. The geometry and boundary conditions of the plate subjected to shear is depicted in Fig. 11. For the first 80 time steps, we take the displacement increment $\Delta u=1\times 10^{-4}$ mm, afterwards, $\Delta u=1\times 10^{-5}$ mm. We calculate the crack patterns under shear loading for two length scale parameters $l_{0}=1.5\times 10^{-2}$ mm and $l_{0}=7.5\times 10^{-3}$ mm as shown in Fig. 12. As expected, the crack propagates under the shear loading and the crack has wider crack width when $l_{0}=1.5\times 10^{-2}$ mm. Figure 12 shows a curved crack path under shear, which is the same as those reported by Miehe et al. [22, 23], Hesch and Weinberg [25], Liu et al. [2].

The load-displacement curves for the top edge of the plate are depicted in Fig. 13 in comparison with the results by Hesch and Weinberg [25]. A close observation is shown in Fig. 13. Thus, the loads obtained by this work and Hesch and Weinberg [25] are exactly matching as the displacement increases, particularly for a smaller length scale $l_{0}=7.5\times 10^{-3}$ mm. The excellent agreement in the crack pattern and load-displacement curves indicates the feasibility and practicability of the presented phase field modeling approach in COMSOL. We then show the influence of mesh size and step size on the load-displacement curves at a fixed $l_{0}$ in Fig. 14. Larger mesh size and displacement increment achieve larger peak load.

4.4 2D and 3D notched square plate subjected to tension and shear

Nooru-Mohammed carried out experiments on double-edge-notched plates subjected to both tension and shear loading [51]. Figure 15 shows the geometry and boundary conditions of their experiments. The length, height and thickness of the plate are 200, 200 and 50 mm, respectively. Two horizontal notches of 25 mm $\times$ 5 mm exist in the middle of the left and right edges of the plate. The shear force $P_{s}$ is applied as Fig. 16. The tensile load is zero with the increase in $P_{s}$ and the vertical displacement increment $\delta$ is then prescribed on the upper and lower boundaries when $P_{s}=5$ kN.

We conduct both 2D and 3D simulations of the Nooru-Mohammed experiment [51]. In the simulations, these elastic constants are used: Young’s modulus $E$ = 32.8 GPa and Poisson’s ratio $\nu=0.2$. We fix the length scale of $l_{0}=2.5$ mm and the plate is discretized into uniform elements with size $h=1.25$ mm ($l_{0}=2h$) for 2D and $h=2.5$ mm for 3D. We choose constant displacement increment $\Delta u=5\times 10^{-6}$ mm for each time step. In this work, we test different critical energy release rates, i.e. $G_{c}=$ 25, 50, 75 and 100 J/m2, respectively.

Figure 17 shows the crack patterns of the 2D notched plate when the displacement is $\delta=0.026$ mm. The crack for lager $G_{c}$ has a smaller inclination angle to the horizontal direction and is more curved. This results in a smaller distance between the two parallel cracks. Figure 18 presents the load-displacement curves of the 2D notched plated under shear and tension for different critical energy release rates. A quick drop of the load after a nearly linear increase is observed. The peak and residual values of the load increases with the increase in the critical energy release rate. Figure 19 shows the crack propagation in the 3D notched square plate for $G_{c}=75$ and 100 J/m2. The crack patterns of the 2D and 3D simulations are in good agreement. The load-displacement curves of the 3D plate are depicted in Fig. 20, which are less steep after the peak compared with the 2D plate. The reason is that the 2D plate section suffers from tension perpendicular to the section under the plane strain assumption, which accelerates the drop in the load bearing capacity of the plate.

4.5 3D three-point bending test

Let us consider a 3D three-point bending test shown in Fig. 21. 2D results were presented by Miehe et al. [22, 23]. The thickness of the beam is 0.4 mm and the following material parameters are used: $\lambda=12$ kN/mm2, $\mu=8$ kN/mm2, and $G_{c}=0.5$ N/mm. Prism elements are used to discretize the beam with a maximum element size of $h=6\times 10^{-2}$ mm except $h=1.5\times 10^{-2}$ mm in the region where the crack is expected to propagate. The specimen is loaded displacement-driven with a constant displacement increment $\Delta u=5\times 10^{-5}$ mm is applied.

We choose the length scale $l_{0}=0.06$ mm and 0.03 mm and show the crack patterns of the simply supported notched beam in Fig. 22. The simulation shows a larger crack width when $l_{0}=0.06$ mm. Figure 23 presents the reaction force at the top of the beam versus the applied displacement. The presented results in 3D are then compared with those 2D results proposed by Miehe et al. [22]. As observed, the 3D results by the present results are in good agreement with the 2D simulations by Miehe et al. [22].

4.6 2D dynamic shear loading of Kalthoff experiment

We next test our method for dynamic fracture by taking advantage of the Kalthoff-Winkler experiments [52] which has been studied by several other researchers [14, 53, 54, 55, 56, 57, 58, 59]. We adopt the symmetry condition to reduce the computational cost and the dimensions and loading conditions are shown in Fig. 24.

The impactor is modelled by applying the following velocity $v$:

[TABLE]

with $v_{0}=16.5$ m/s and $t_{0}=1$ $\mu$s. Moreover, the initial crack is assumed to be traction free.

The material parameters are from [24]: $\rho=8000$ kg/m3, $E=190$ GPa, $\nu=0.3$, $G_{c}=2.213\times 10^{4}$ J/m2, and $k=1\times 10^{-9}$. The Rayleigh wave speed of the plate is $v_{R}=2803$ m/s. The length scale parameter $l_{0}$ is fixed as $3.9\times 10^{-4}$ m. The initial crack is modeled as a notch with a width of $l_{0}$. Here, we use Q4 elements to discretize the plate and we simulate the crack patterns for two mesh levels: Mesh 1 with element size $h=3.9\times 10^{-4}$ m ( $l_{0}=h$) and Mesh 2 with element size $h=1.95\times 10^{-4}$ m ($l_{0}=2h$). The time steps are set as $\Delta t=0.04$ $\mu$s (for Mesh 1, $\Delta t=0.288h/v_{R}$ and for Mesh 2, $\Delta t=0.576h/v_{R}$) and $\Delta t=0.01$ $\mu$s (for Mesh 1, $\Delta t=0.072h/v_{R}$ and for Mesh 2, $\Delta t=0.144h/v_{R}$), respectively.

The phase-field of our simulation at 90 $\mu$s is shown in Fig. 25 for different mesh levels and time steps. The crack starts to propagate at 26 $\mu$s. The crack angle versus the horizontal axis varies from $63^{\circ}$ to $67^{\circ}$ which matches well the experimental results and other numerical results [24]. As shown in Fig. 25, the crack tip by the coarser mesh (Mesh 1) has larger distance from the upper boundary than that by the finer mesh (Mesh 2). The crack has a larger angle and goes more close to the upper boundary for the smaller time step $\Delta t=0.01$ $\mu$s.

Figure 26 shows the contour plot of the maximum principal tensile stress for Mesh 2 at $t=75$ $\mu$s under different time steps. Note that the deformation is scaled by a factor of 5. The region with $\phi>0.95$ is also removed from the figure to see the broken geometry of the plate.

We also obtain the velocity of crack tip as shown in Fig. 27. The velocity $\bm{v}_{n}$ at the crack tip is calculated as follows:

[TABLE]

with $\bm{x}_{n}$ the position of current crack tip at time $t_{n}$. The position of crack tip is determined from the iso-curve of the phase-field $\phi=0.75$ according to Borden et al. [24]. The crack speed increases to a velocity of $0.6v_{R}$ and remains nearly constant until the end of the simulation. Figure 27 also shows that the results are independent for two different time steps.

Figures 28 and 29 present the elastic strain energy and dissipated energy curves, respectively. The elastic strain energy is calculated as

[TABLE]

while the dissipated energy is obtained by

[TABLE]

As shown in Fig. 28, the elastic strain energy curves for Mesh 1 and Mesh 2 are in good agreement. The elastic strain energy increases quickly before reaching a maximum. After that, the elastic strain energy starts to decrease slowly. In Fig. 29, the dissipated energy increases as the time increases. .

In the end of this example, we test the influence of the critical energy release rate $G_{c}$. We present the results of $G_{c}=5\times 10^{3}$, $1\times 10^{4}$, $2.213\times 10^{4}$, and $3\times 10^{4}$ J/m2 for Mesh 1 and $\Delta t=0.04$ $\mu$s. Figure 30 gives the phase field at 90 $\mu$s for different $G_{c}$. More complex crack patterns are observed for smaller $G_{c}$. Crack branching occurs when $G_{c}=5\times 10^{3}$ and $1\times 10^{4}$ J/m2. Secondary cracks occur in the bottom right corner of the model because of wave reflections [21]. However, for a larger $G_{c}$, only a single crack is seen. The crack is more hard to reach the upper boundary of the plate when $G_{c}$ becomes larger. The variation in the critical energy release rate $G_{c}$ does not change the main crack pattern. The main cracks initiate from the tip of the pre-existing crack.

Figure 31 presents the crack-tip velocity for different $G_{c}$. The velocity is calculated on the longest cracks in Fig. 30. Figure 31 shows that the maximum crack-tip velocity decreases and the time for crack initiation increases with the increase in $G_{c}$. Particularly, for $G_{c}=5\times 10^{3}$ J/m2, two crack branching successively occur when the crack-tip velocity reaches the maximum.

4.7 2D dynamic crack branching under tension

The last example is another classical benchmark problem for dynamic fracture: a pre-notched rectangular plate subjected to uni-axial traction. Figure 32 gives the geometry of the plate along with the boundary conditions. This benchmark test has been calculated by Song et al. [21], Liu et al. [2] and Rabczuk and Belytschko [60, 61].

The material parameters are adopted from [24] as $\rho=2450$ kg/m3, $E=32$ GPa, $\nu=0.2$ and $G_{c}=3$ J/m2. Plane strain conditions are assumed. The length scale parameter $l_{0}$ is fixed as $5.0\times 10^{-4}$ m. The pre-existing notch is modeled by introducing an initial history strain field as explained in Section 3.2. Q4 elements are used to discretize the plate with a uniform mesh; $k=1\times 10^{-9}$ is also picked to avoid the singularity during calculation. We conduct the simulation using two different mesh levels: Mesh 1 with size $h=2.5\times 10^{-4}$ m ($l_{0}=2h$) and Mesh 2 with $h=1.25\times 10^{-4}$ m ($l_{0}=4h$), respectively. For each mesh, the time step is chosen as 0.1 $\mu$s, 0.05 $\mu$s and 0.025 $\mu$s.

Figure 33 shows the results of the phase field at 80 $\mu$s. Mesh 1 and Mesh 2 have similar crack patterns for different time steps and the crack fail to reach the boundary because the staggered scheme is adopted. As shown in [2] and [24], the monolithic scheme is easy to reach the boundary. Thus, the staggered scheme in this work seems to delay the crack branching. In addition, Fig. 33 shows that a coarser mesh and a larger time step have wider cracks and the cracks are much more difficult to reach the boundary than a finer mesh and a smaller time step.

Figure 34 gives the maximum tensile stress for Mesh 1 and Mesh 2 at $t=70$ $\mu$s. In the figure, we scale the displacement field by a factor of 100 and the regions where the phase field is more than 0.95 are also removed from the figure to display the broken geometry of the plate. Figure 34 shows a tensile stress concentration at the crack tip and the results of both meshes are similar and in good agreement with the results of Liu et al. [2].

Figure 35 presents the elastic strain energy curves for all the meshes and time steps. All the curves are in good agreement before 60 $\mu$s. However, some discrepancies exist after 60 $\mu$s. For a smaller time step, the elastic strain energy approximately decreases with the increase in time. But for a larger time step, the elastic strain energy increases as the time goes. Additionally, the elastic strain energy for the coarser mesh (Mesh 1) is larger than that for the finer mesh (Mesh 2). Figure 36 represent the dissipated energy curves for the 2D crack branching example. In Fig. 36, the dissipated energy increases as the time increases. Meanwhile, the results for both meshes and all the time steps are in quite good agreement.

Figure 37 presents the curves of crack tip velocity achieved from the post-processing. As depicted in previous example, the crack tip is found in the iso-curve of the phase-field $\phi=0.75$. All the curves in Fig. 37 have the similar trend for the crack branching case. As we notice in the simulation, all the velocities are smaller than $0.5v_{R}$. This finding is in good agreement with the results of Borden et al. [24]. In addition, the crack widening occurs when the time $t=28$ $\mu$s $\sim$ 30 $\mu$s. The crack starts to branch when the time goes to 48 $\mu$s $\sim$ 51 $\mu$s, which is later than the results of Borden et al. [24]. This also proves that the staggered scheme can delay the time for the crack branching.

To end the example of 2D crack branching, we test the influence of the critical energy release rate $G_{c}$ on the crack pattern and the crack-tip velocity as shown in Figs. 38 and 39. The results are now presented only using Mesh 1 and $\Delta t=0.1$ $\mu$s. When $G_{c}=0.5$ J/m2, multiple crack branching can be seen in Fig. 38. Figure 38 also shows that, for a larger $G_{c}$, the crack will propagate at a larger angle with the horizontal after the first branching. The crack is more hard to propagate as well. Figure 39 depicts the crack-tip velocity for different $G_{c}$. The maximum crack-tip velocity decreases with the increase in $G_{c}$. Larger $G_{c}$ will cause crack propagation at a larger velocity, which is in good agreement with the 2D dynamic shear example.

5 Conclusions

In this work, we present the implementation of the numerical phase-field modeling for crack propagation in commercial finite element software COMSOL. The previous existing phase-field models for quasi-static and dynamic crack propagation problems are reviewed and implemented in proper forms in COMSOL. The elastic strain energy density is decomposed into two individual parts resulting from compression and tension, respectively. Thus, only the tension induced cracks are obtained. In COMSOL, the simulation is facilitated by establishing four modules and placing the coupling terms appropriately. The coupled system is solved using a staggered scheme. Here, an implicit time integration scheme is used and employ the Newton-Raphson iteration is adopted to compute the individual fields.

A number of 2D and 3D benchmark examples for quasi-static and dynamic crack propagations are tested. We also check the effects of length scale parameter, critical energy release rate, mesh size, and time step size on the results. All the simulations have the correct crack patterns and satisfactory accuracy, showing the feasibility of implementing phase-field model for crack propagation by COMSOL, even in 3D spaces. For future work, COMSOL will be more helpful and effective for implementing and extending the phase-field model to problems with more fields.

Acknowledgement

The financial support provided by the Sino-German (CSC-DAAD) Postdoc Scholarship Program 2016, the Natural Science Foundation of China (51474157), and RISE-project BESTOFRAC (734370) is gratefully acknowledged.