High-order phase-field model with the local and second-order max-entropy approximants

Abstract

We approximate the fracture surface energy functional based on phase-field method with smooth local maximum entropy (LME) and second-order maximum entropy (SME) approximants. The higher-order continuity of the meshfree methods such as LME and SME approximants allows to directly solve the fourth-order phase-field equations without splitting the fourth-order differential equation into two second-order differential equations. We will first show that the crack surface functional can be captured more accurately in the fourth-order model with smooth approximants such as LME, SME and B-spline. Furthermore, smaller length scale parameter is needed for the fourth-order model to approximate the energy functional. We also study SME approximants and drive the formulations. The proposed meshfree fourth-order phase-field formulation show more stable results for SME compared to LME meshfree methods.

Keywords

second-order maximum entropy / local maximum entropy / second- and fourth-order phase-field models / B-spline

Introduction

During the last few decades the numerical simulation of scientific problems has gained importance and often plays a key role in design decisions [¹,²]. This has been mainly motivated by the impossibility to have analytical solutions in most practical situations and the costs of obtaining meaningful and detailed information from experiments. Numerical methods, such as the finite element method have been used with some degrees of success, but often they are unable to capture some physical properties of the phenomenon. For instance, modeling of moving discontinuities with classical finite elements is difficult to automate because of the requirement that the mesh must conform to the surfaces of discontinuity. It also usually requires local refinement near the fracture zone, in particular, near the crack tips where singularities in the stress field occur [³–⁶]. The numerical methods to model fracture have centralized on two main approaches: Discrete and continuum damage descriptions. Examples of the discrete description of fracture include the extended finite element method (XFEM) [⁷,⁸] or the generalized finite element method (GFEM) [⁹] and numerical manifold method [¹⁰]. These methods allow arbitrary propagating discontinuities without remeshing. One key challenge of such methods is describing the crack geometry and tracking the paths of the cracks as the fracture progresses. This becomes increasingly challenging for complex fracture patterns. In general, these numerical approximations track the evolution of the fracture during the simulations, but they have shown to be inefficient regarding, for example, crack branching in three-dimensional applications. XIGA formulations (extended isogeometric analysis) for fracture [¹¹–¹³] aim to combine the advantages of isogeometric analysis and the extended finite element method. However, they also do not resolve the issue of complex track cracking procedures. An alternative to model complex fracture are meshfree methods [¹,¹⁴–¹⁸]. While those contributions also rely on the representation of the crack surface, some meshfree methods such as the cracking particles method [¹⁹,²⁰] model fracture as a set of crack segments and therefore can capture—on cost of the accuracy in the crack kinematics—complex fracture patterns quite naturally. Besides discrete crack models, continuous descriptions of fracture in solids have been presented. Among the most popular approaches are gradient models [²¹–²⁶] and non-local models [²⁷–³⁰]. They introduce an intrinsic length scale and diffuse fracture over a certain width. Recently, the screened-Poisson equation method was introduced in Refs. [³¹,³²] that can be classified as either discrete or continuum-based. However, this method uses standard shape functions which is not appropriate for problems that require smooth shape functions.

In recent years, phase-field methods are widely used in science and engineering to model a variety of physics [³³,³⁴]. In the phase-field approach, the problem is reformulated in terms of a coupled system of partial differential equations (PDEs). A continuous scalar-valued phase-field is introduced into the model to indicate the material phases. The evolution of the phase-field is governed by a PDE that leads to a coupling between the bulk and the surface energies. In this model, the proposed energy functionals closely similar to the potential functional presented by Mumford and Shah [³⁵], which has been used in image segmentation. The existence of solution to the Mumford-Shah functional minimization was proven by Ambrosio in Ref. [³⁶]. In Ref. [³⁷], an approximation by an elliptic functional defined on Sobolev spaces was developed, based on the theory of Г-convergence. In general, the phase-field approach to model systems with sharp interfaces consists in incorporating a continuous field variable, the so-called phase-field variable, which differentiates between multiple physical phases within a given system through a smooth transition. Also, the phase-field model is used in phase transition problems such as Cahn-Hilliard equation [³⁸], image information or image inpainting [³⁹], multi-phase systems and fracture mechanics [⁴⁰]. The predecessors of phase-field approaches to fracture can be traced back to 1998 in Refs. [⁴¹,⁴²], where the brittle crack propagation problem was regularized and recast as a minimization problem. Phase-field approaches for fracture [⁴³] bear certain similarities to gradient models but they converge to a discrete crack model when the characteristic length tends to zero. The main drawback of the phase-field method is the higher computational cost of solving a coupled PDEs system.

Dealing computationally with some physical and engineering problems such as Kirchhoff-Love thin shell theory, is challenging. This is due to of appearing the high-order derivatives in the weak form, and therefore the Galerkin method requires at least C¹-continuous basis functions. To overcome this difficulty, the smooth basis functions such as isogeometric analysis [⁴⁴–⁴⁶] and maximum entropy [⁴⁷], were introduced. Maximum entropy basis functions are a relatively new class of approximation functions, as they were first introduced in the context of polygonal interpolation. The idea of these functions is to maximize the Shannon entropy [⁴⁸] of the basis functions, which gives a measure of the uncertainty in the approximation scheme. The principle of maximum entropy maximum entropy (max-ent) was developed by Jaynes [⁴⁹,⁵⁰], who showed that there is a natural correspondence between statistical mechanics and information theory. In particular, max-ent offers the least-biased statistical inference when the shape functions are viewed as probability distributions subject to the approximation constraints (such as linear reproducing properties). However, without additional constraints, the basis functions are non-local, which due to increased overlapping makes them unsuitable for analysis using Galerkin methods. The large overlapping of the basis functions, generally leads to more expensive numerical integration due to large number of evaluation points. It also, produces non-sparse stiffness matrix which causes more expensive linear system to solve. The local maximum-entropy (LME) approximation schemes were developed in Ref. [⁵¹] using a framework similar to meshfree methods. Here, the support of the basis functions is introduced as a thermalization (or penalty) parameter b in the constraint equations. When b= 0, then the max-ent principle is fully satisfied and the basis functions will be least-biased. For example, if only zeroth-order consistency is required, the shape functions are Shepard approximants [⁵²] with Gaussian weight function. When b is large, then the shape functions have minimal support. In particular, they become the usual linear finite element functions defined on a Delaunay triangulation of the domain associated with the given node set. In Ref. [⁵¹] it was shown that for some values of b, the approximation properties of the max-ent basis functions are greatly superior to those of the finite element linear functions, even when the added computational cost due to larger support is taken into account. The second-order maximum entropy (SME) method was developed in Refs. [⁵¹,⁵³] satisfying the second-order consistency condition. Though, as shown in Refs. [⁵¹,⁵³] that the second-order consistency condition makes the constraints unfeasible in general. Some methods were suggested to deal with this obstacle [⁵⁴,⁵⁵]. However, the presented methods made some other difficulties by themselves. Here, we follow Ref. [⁵³] which considers extra non-negative parameters in the second-order consistency condition to overcome this problem. Then, the locality of the SME basis functions is controlled by parameter a. When a is large, the support size of the SME basis functions is large, and when the parameter a is small the support is small as well. Subsequent studies [¹,¹⁴,⁴⁰,⁵⁵–⁵⁸], show that max-ent shape functions are suitable for solving a variety of problems such as linear and geometrically nonlinear thin shell analysis, compressible and nearly-incompressible elasticity and incompressible media problems.

The LME and SME approximations have several advantages over other meshfree methods such as the element-free Galerkin method [⁵⁹] or reproducing kernel particle method as demonstrated by numerical examples in Refs. [⁴⁷,⁵¹,⁵⁵,⁵⁸,⁶⁰]. For example, in contrast to above mentioned methods,

● LME and SME approximations fulfill the weak Kronecker-Delta properties facilitating the imposition of Dirichlet boundary conditions;

● Their shape functions are always positive which leads to non-negative values in the off-diagonals in the mass matrix;

● They require less integration points for the same accuracy.

In this paper, we study the performance of a fourth-order phase-field model introduced in Ref. [⁴⁰]. We show it is more efficient to use fourth-order phase-field model with LME and SME basis functions, due to the smoothness of these shape functions.

The paper is organized as follows. The brief reviews on LME and SME basis functions are presented in Sections 2 and 3, respectively. In Section 4, we present the general theory and motivation for the second and fourth-order phase-field models. In Section 5, we demonstrate the capabilities of the method through some numerical studies. In Section 5, the convergence rate and the error of modeling of the fourth-order phase-field model are compared to the second-order phase-field model, for LME, SME and B-spline. Some concluding remarks are given in Section 6.

Local maximum entropy

Max-ent shape functions are a relatively new class of approximation functions, as they were first introduced in Ref. [⁴⁷] in the context of polygonal interpolation. The idea of these functions is to maximize the Shannon entropy [⁶¹] of the basis functions, which gives a measure of the uncertainty in the approximation scheme. The principle of max-ent was developed by Jaynes [⁴⁹,⁵⁰], who showed that there was a natural correspondence between statistical mechanics and information theory. In particular, max-ent offers the least-biased statistical inference when the shape functions are viewed as probability distributions subject to the approximation constraints (such as linear reproducing properties). However, without additional constraints, the basis functions are non-local, which due to increased overlapping makes them unsuitable for analysis using Galerkin methods. The large overlapping of the basis functions generally leads to more expensive numerical integration schemes due to large number of evaluation points. It also produces non-sparse stiffness matrix which require significantly more computational resources to solve. The LME meshfree approximants, introduced in Ref. [⁵¹], are related to other convex approximation schemes, such as natural neighbor approximants [⁶²], subdivision approximants [⁶³], or B-spline and non-uniform rational B-splines (NURBS) basis functions [⁴⁶]. These basis functions are non-negative and satisfy the zeroth-order and first-order consistency conditions.

(1)p_a(x)\ge \mathrm{0},

(2){\displaystyle \sum _{a=\mathrm{1}}^N{p}_a(x)}=\mathrm{1},

(3){\displaystyle \sum _{a=\mathrm{1}}^N{p}_a(x)}{x}_a=x\mathrm{.}

The main idea of max-ent is to maximize the Shannon’s entropy, H(p_{\mathrm{1}},p_{\mathrm{2}},\hspace{0.17em}\mathrm{...},p_N), subject to the consistency constraints as follows:

(4)\{\begin{array}{l}\mathrm{(max-ent)}\hspace{0.17em}\mathrm{For}\hspace{0.17em}a\hspace{0.17em}\mathrm{fixed}\hspace{0.17em}x,\hspace{0.17em}\mathrm{maximize}\\ H(p_{\mathrm{1}},p_{\mathrm{2}},\mathrm{...},p_N)=-{\displaystyle \sum _{a=\mathrm{1}}^N{p}_a\log p_a}\\ \mathrm{subject}\hspace{0.17em}\mathrm{to}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{p}_a\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}a=\mathrm{1},\mathrm{2},\mathrm{...},N\\ {\displaystyle \sum _{a=\mathrm{1}}^N{p}_a}=\mathrm{1}\\ {\displaystyle \sum _{a=\mathrm{1}}^N{p}_a{x}_a}=x\end{array}\mathrm{.}

Solving the (max-ent) problem produces the set of basis functions, a=1,\mathrm{2},\mathrm{...},N. However, these basis functions are non-local, i.e., they have support in all of convX, and are not suitable for use in a Galerkin approximation because it would lead to a full, non-banded matrix. Nevertheless, they have been used in Ref. [⁴⁷] as basis functions for polygonal elements.

Another optimization problem which takes into account the locality of the shape functions is Rajan’s form of the Delaunay triangulation [⁶⁴]. This can be stated as the following linear program:

(5)\{\begin{array}{l}(\mathrm{Rajan’s\; form})\mathrm{For}\hspace{0.17em}a\hspace{0.17em}\mathrm{fixed}\hspace{0.17em}x\hspace{0.17em}\mathrm{minimize}\\ \mathrm{subject\; to}\hspace{0.17em}{p}_a\ge \mathrm{0},a=\mathrm{1},\mathrm{2},\mathrm{...},N\\ {\displaystyle \sum _{a=\mathrm{1}}^N{p}_a}=\mathrm{1}\\ {\displaystyle \sum _{a=\mathrm{1}}^N{p}_a{x}_a}=x\end{array}\mathrm{.}

It is easy to see that U(x,p_{\mathrm{1}},p_{\mathrm{2}},\mathrm{...},p_N) is minimized when the shape functions p_{\mathrm{1}},p_{\mathrm{2}},\mathrm{...},p_N decay rapidly as the distance from the corresponding nodes x_a increases. There, the shape functions that satisfy (Rajan’s form) problem will have small supports, where the support can be defined up to a small tolerance \epsilon by \mathrm{supp}(p_a)=\{x:p_a(x)>\epsilon \}.

The main idea of LME approximants is to compromise between the max-ent problem and the (Rajan’s form) problem by introducing parameters {\beta }_a that control the support of the p_a. Therefore,

(6)\{\begin{array}{l}\mathrm{For}\hspace{0.17em}a\hspace{0.17em}\mathrm{fixed}\hspace{0.17em}x\hspace{0.17em}\mathrm{minimize}\\ {\displaystyle \sum _{a=\mathrm{1}}^N{\beta }_a{p}_a{\left|x-x_a\right|}^{\mathrm{2}}}+{\displaystyle \sum _{a=\mathrm{1}}^N{p}_a\log p_a}\\ \mathrm{subject}\hspace{0.17em}\mathrm{to}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{p}_a\ge \mathrm{0},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}a=\mathrm{1},\mathrm{2},\mathrm{...},N\\ {\displaystyle \sum _{a=\mathrm{1}}^N{p}_a}=\mathrm{1}\\ {\displaystyle \sum _{a=\mathrm{1}}^N{p}_a{x}_a}=x\end{array}\mathrm{.}

Finally, the local max-ent basis functions are written in the form:

p_a(x)={\displaystyle \frac{\mathrm{1}}{Z(x,\lambda *(x))}}\exp [-{\beta }_a|x-x_a{|}^{\mathrm{2}}+\lambda *(x)\cdot (x-x_a)],

where Z(x,\lambda )={\displaystyle \sum _{b=\mathrm{1}}^N\exp [-{\beta }_b{\left|x-x_b\right|}^{\mathrm{2}}+\lambda (x-x_b)]} is a function associated with a set of nodes X={\{x_a\}}_{a=\mathrm{1},\mathrm{2},\mathrm{...},N} and \lambda *(x) is defined by {\lambda }^{*}\left(x\right)=\arg \underset{\lambda \in {\mathbb{R}}^d}{\min }\log Z(x,\lambda ).

This optimization problem is efficiency solved with Newton Raphson’s method.

The LME basis functions are smooth (C^{\infty }) [⁵¹]. In this paper, we choose \beta ={\displaystyle \frac{\gamma }{h^{\mathrm{2}}}}, where h is a measure of the nodal spacing and \gamma is constant over the domain. In this case, the basis functions are smooth and their degree of locality is controlled by the parameter \gamma. As it was studied in Refs. [¹,⁵¹], the optimal support size of the basis functions or \gamma value is problem dependent. It is also shown that choosing \gamma \ge \mathrm{4.8} gives almost the same results as standard finite element method. While decreasing \gamma gives smoother LME basis functions. The LME basis functions satisfy a weak Kronecker Delta property at the boundary of the convex hull of the nodes. Therefore, the basis functions that correspond to interior nodes vanish on the boundary.

Second-order maximum entropy

The max-ent basis functions can be designed to satisfy quadratic consistency as well,

(7)s_a(x)\ge \mathrm{0},

(8){\displaystyle \sum _{a=\mathrm{1}}^N{s}_a(x)}=\mathrm{1},

(9){\displaystyle \sum _{a=\mathrm{1}}^N{s}_a(x)x_a}=x,

(10){\displaystyle \sum _{a=\mathrm{1}}^N{s}_a(x)x_a^{\mathrm{2}}}=x^{\mathrm{2}}\mathrm{.}

However, it was shown in Ref. [⁵³], there is no such set of convex approximants to fulfill Eq. (10) except forx=x_a. According to Ref. [⁵³], the second order consistency constraint can be replaced by

{\displaystyle \sum _{a=\mathrm{1}}^N{s}_a(x)(x_a^{\mathrm{2}}-d_a)}=x^{\mathrm{2}},

where d_a is a non-negative parameter or in multi-dimensional are symmetric semi-positive de finite matrices for each node. Therefore, the optimization problem to obtain SME approximants is written as follows:

(11)\{\begin{array}{l}\mathrm{(SME)}\hspace{0.17em}\mathrm{For}\hspace{0.17em}a\hspace{0.17em}\mathrm{fixed}\hspace{0.17em}x\hspace{0.17em}\mathrm{maximize}\\ H(s_{\mathrm{1}},s_{\mathrm{2}},\mathrm{...},s_N)=-{\displaystyle \sum _{a=\mathrm{1}}^N{s}_a\log s_a}\\ \mathrm{subject\; to}\hspace{0.17em}{s}_a\ge \mathrm{0},a=\mathrm{1},\mathrm{2},\mathrm{...},N\\ {\displaystyle \sum _{a=\mathrm{1}}^N{s}_a}=\mathrm{1}\\ {\displaystyle \sum _{a=\mathrm{1}}^N{s}_a(x_a-x)}=\mathrm{0}\\ {\displaystyle \sum _{a=\mathrm{1}}^N{s}_a\{{(x_a-x)}^{\mathrm{2}}-d_a\}}=\mathrm{0}\end{array}\mathrm{.}

Following the Lagrange multipliers method as in Ref. [⁵³], the second-order max-ent basis functions are written in the form:

s_a(x)={\displaystyle \frac{\mathrm{1}}{Z(x,\lambda *(x),\mu *(x))}}\exp [\lambda *(x)\cdot (x-x_a)-\mu *(x):{\mathit{D}}_a],

where

Z(x,\lambda ,\mu )={\displaystyle \sum _{b=\mathrm{1}}^N\exp [\lambda \cdot (x-x_b)-\mu :{\mathit{D}}_a]},

(12){\mathit{D}}_a={(x_a-x)}^2-d_a,

and \lambda *(x), \mu *(x) are the optimal Lagrange multipliers obtained by minimizing the reduced dual function

(13)g(\lambda ,\mu )=\mathrm{log\; Z}(x,\lambda ,\mu )\mathrm{.}

The locality of the SME basis functions, is controlled by slack parameter \alpha, which d_a={\displaystyle \frac{\alpha }{\mathrm{4}}}\hspace{0.17em}\max \{h_{a-\mathrm{1}}^{\mathrm{2}},h_a^{\mathrm{2}}\} for a=\mathrm{3},\mathrm{4},\mathrm{...},N-\mathrm{1} and for the nodes located next to the boundary is: d_{\mathrm{2}}=\max \{h_{\mathrm{1}}^{\mathrm{2}},{\displaystyle \frac{\alpha }{\mathrm{4}}}\hspace{0.17em}{h}_{\mathrm{2}}^{\mathrm{2}}\} and d_{N-\mathrm{1}}=\max \{h_{N-\mathrm{1}}^{\mathrm{2}},{\displaystyle \frac{\alpha }{\mathrm{4}}}\hspace{0.17em}{h}_{N-\mathrm{2}}^{\mathrm{2}}\}. When \alpha is large, the basis functions have large support size on the domain and for the small values of \alpha, the support size is small. The range of values \mathrm{1.2}\le \alpha \le \mathrm{3} are recommended by Ref. [⁵³]. Also, the SME basis functions are smooth (C^{\infty }), non-negative and satisfy the weak Kronecker Delta property.

Second- and fourth-order phase-field model

Consider an in finite bar with cross-section \Gamma, occupying a region \Omega =(-\infty ,+\infty )\times \Gamma. According to Fig. 1(a), a crack is placed at x = 0. The phase-field variable v(x)\in [\mathrm{0},\mathrm{1}] with

(14)v(x)=\{\begin{array}{ll}\mathrm{0}& x=\mathrm{0}\\ \mathrm{1}& \mathrm{otherwise}\end{array},

is introduced to describe the crack topology, where v = 0 indicates the crack (total damage) while v = 1 refers to intact state. This phase-field is discontinuous at x = 0 and satisfies the following conditions

(15)\{\begin{array}{l}v(\mathrm{0})=\mathrm{0}\\ v(\pm \infty )=\mathrm{1}\end{array}\mathrm{.}

A function that fulfils the criterion Eq. (15) and \tilde{v}(x)\to v(x), when the diffusivity l_{\mathrm{0}}\to \mathrm{0}, is

(16)\tilde{v}(x)=\mathrm{1}-e^{-\frac{\left|x\right|}{\mathrm{2}{l}_{\mathrm{0}}}}\mathrm{.}

Furthermore, Eq. (16) is the solution of a second-order differential equation:

(17)\tilde{v}(x)-\mathrm{4}{l}_{\mathrm{0}}^{\mathrm{2}}{\tilde{v}}^{″}-\mathrm{1}=\mathrm{0},\mathrm{i}\mathrm{n}\mathrm{\Omega },

subject to the essential boundary conditions given in Eq. (15). This second-order differential equation is the Euler equation of the variational principle

(18)\tilde{v}=\arg \left\{\underset{\tilde{\tilde{v}}\in W}{\inf }I(\tilde{v})\right\}

where W=\{\tilde{v}|\tilde{v}(\mathrm{0})=\mathrm{0},\tilde{v}(\pm \infty )=\mathrm{1}\} and

(19)I\tilde{v}=\frac{1}{4}{\int }_{-\infty }^{+\infty }[{(1-\tilde{v})}^2+{4l}_0^2{\tilde{v}}^{\prime 2}]\mathrm{d}x,

The value of Eq. (16) is

(20)I\left(\tilde{v}(x)=\mathrm{1}-e^{-\frac{\left|x\right|}{\mathrm{2}{l}_{\mathrm{0}}}}\right)=l_{\mathrm{0}}\Gamma ,

which computes the area of the crack by

(21){\mathrm{\Gamma }}_{l_0}=\frac{1}{{4l}_0}{\int }_{-\infty }^{+\infty }[{(1-\tilde{v})}^2+{4l}_0^2{\tilde{v}}^{\prime 2}]\mathrm{d}x,

with dV=\Gamma dx. Hence, the second-order differential equation leads to second-order phase-field method to compute the crack surface. In order to approximate function \tilde{v}, a C^{\mathrm{0}} basis functions such as the one provided by the standard finite element method (FEM) is sufficient. In meshfree methods such as LME [⁵¹] and SME [⁵³] approximants and B-spline [⁴⁶] interpolants, the basis functions are smooth. Hence, the kink (at x = 0) in the phase-field might not be well-captured with these methods. Let us assume another approximation function [⁶⁵]

(22)\tilde{v}(x)=\mathrm{1}-e^{-\frac{\left|x\right|}{l_{\mathrm{0}}}}\left(\mathrm{1}+{\displaystyle \frac{\left|x\right|}{l_{\mathrm{0}}}}\right)\mathrm{.}

This function satisfies the conditions Eq. (15) and introduces smoother approximation as is demonstrated in Fig. 1(c). We find that Eq. (22) is the solution of the fourth-order differential equation

(23)\tilde{v}-\mathrm{1}-\mathrm{2}{l}_{\mathrm{0}}^{\mathrm{2}}{\tilde{v}}^{\prime \prime }+l_{\mathrm{0}}^{\mathrm{4}}{\tilde{v}}^{(\mathrm{4})}=\mathrm{0},\mathrm{i}\mathrm{n}\mathrm{\Omega },

under the assumption that \tilde{v}(\mathrm{0})=\mathrm{0}, {\tilde{v}}^{\prime }(\mathrm{0})=\mathrm{0} and {\tilde{v}}^{\alpha }(x)\to \mathrm{0} as x\to \pm \infty for all \alpha \ge \mathrm{0} [⁶⁵]. The functional behind these ordinary differential equations is

(24)I\tilde{v}-\frac{1}{2}{\int }_{-\infty }^{+\infty }[{(1-\tilde{v})}^2+\mathrm{2}{l}_{\mathrm{0}}^{\mathrm{2}}{\tilde{v}}^{\prime 2}+l_{\mathrm{0}}^{\mathrm{4}}{\tilde{v}}^{\prime \prime 2}]\mathrm{d}x\mathrm{.}

The value of this functional for Eq. (22) is \mathrm{2}{l}_{\mathrm{0}}\Gamma, which gives

(25){\mathrm{\Gamma }}_{l_0}=\frac{1}{{4l}_0}{\int }_{-\infty }^{+\infty }[{(1-\tilde{v})}^2+\mathrm{2}{l}_{\mathrm{0}}^{\mathrm{2}}{\tilde{v}}^{\prime 2}+l_{\mathrm{0}}^{\mathrm{4}}{\tilde{v}}^{\prime \prime 2}]\mathrm{d}x,

which is the fourth-order phase-field model to compute crack surface.

Finite elements based on Lagrange polynomials are not well suited for the fourth-order phase-field model which requires C^{\mathrm{1}} continuity in the phase-field. A common approach is to split the fourth-order differential equation into two second-order differential equations as proposed e.g., in Ref. [⁶⁶]. Due to the higher order continuity of the LME and SME approximants, they are ideally suited to directly solve the fourth-order phase-field equations. Such an approach has been proven to be successful in isogeometric analysis (IGA) [⁶⁵]. However, in contrast to IGA which requires at least a quadratic basis, only linear completeness is needed in the LME approximants.

Numerical examples

Solution of phase-field equation

First, we determine how well the proposed second and fourth order models approximate the “crack topology”. For simplicity, we consider the 1D model explained in Section 4, with L = 1. We set the Dirichlet boundary condition v = 0, for x = 0. We perform calculations with the second- and fourth-order models for a fixed l₀= 0.1, and vary grid spacing h. We also consider several values of \gamma and for LME and SME, respectively. In Figs. 2 and 3, Eqs. (19) and (24) are considered as energy for the second- and fourth-order phase-field models, respectively. We expect, I(\tilde{v})\approx l_{\mathrm{0}}\Gamma for the second-order phase-field formulation. Here \Gamma=1. In Fig. 2(a), for LME with \gamma \ge \mathrm{5.8} the results are acceptable. For \gamma= 5.8, the LME basis functions are much sharper and rigorously approach the C^{\mathrm{0}} finite element shape functions for \gamma \to \infty. As indicated in Fig. 2(b), the second-order phase-field formulation is not proper to be used with SME approximants. This is due to the smooth basis functions used to approximate a non-smooth solution. In Fig. 3 with fourth-order phase-field functional, for both LME with small values of \gamma and SME approximants, \mathrm{2}h\le l_{\mathrm{0}} resolves the regularized crack surface such that I(\tilde{v})\approx \mathrm{2}{l}_{\mathrm{0}}\Gamma.

Next, we compute the error in the L^{\mathrm{2}} norm and H^{\mathrm{1}} semi-norm for the second-order model by considering Eq. (16) as exact solutions. As it is obvious from Fig. 4, the best results for LME, are obtained for \gamma= 5.8. The rate of convergence in the L^{\mathrm{2}} norm in this case is about 1.7, while the rate of convergence for \gamma<2.8 is about 1.0. The convergence rate in the H^{\mathrm{1}} semi-norm is optimal, i.e., equal 1, with \gamma= 5.8 but sub-optimal (order 0.4) for \gamma \le2.8. The results for SME, are shown in Fig. 5. In this case, the convergence rate for L^{\mathrm{2}} norm is optimal and for the H^{\mathrm{1}} semi-norm is sub-optimal. When \alpha decreases, the support size of the SME basis functions also decreases, but unlike LME basis function, SME basis functions do not approach FEM. In view of the fact, even when \alpha is small the SME basis functions are smooth on peak and satisfy the second order consistency. This study shows that the second-order phase-field model is not suitable to be used with max-ent approximants.

Finally, we compute the error in the L^{\mathrm{2}} norm, H^{\mathrm{1}} semi-norm and H^{\mathrm{2}} semi-norm for the fourth-order phase-field model with the LME and SME approximants and B-spline interpolants. For the LME basis functions, the best results are obtained for \gamma= 0.6, and the rate of convergence of the L^{\mathrm{2}} norm error is about 4, which is super convergent. As can be seen from Fig. 6, for values of \gamma \ge0.8 no convergence in the L^{\mathrm{2}} norm is obtained. Note that, in contrast to the at least quadratic NURBS-formulation, we employ basis functions that possess up to first-order reproducing conditions. Strictly speaking, to solve the fourth-order PDE, functions that reproduce quadratic polynomials are required [⁶⁷], and hence we should use quadratic max-ent approximants. However, it has already been shown that, in practice, linear max-ent approximants with low value of \gamma effectively behave like quadratic consistent approximants, and the convergence rate in fourth-order PDE (thin shells and phase-field models of membranes) is only degraded for very fine meshes and very small errors [¹⁴,⁵⁷]. Figure 6 agrees with this behavior, in that convergence is degraded only for high values of \gamma and for very fine grids. However, we also consider quadratic or second-order max-ent approximants. Figure 7 indicates the results for SME with different values of \alpha. In this case, due to use of the quadratic basis functions to solve the fourth-order PDE, as it should be, the degraded does not occur. As illustrated in Fig. 7, the accuracy of the results increases as the locality of the SME basis functions increases, the value of \alpha decreases. The results for the LME basis functions with \gamma= 0.6, closely follow that of SME approximants with \alpha= 2.5.

Figure 8 illustrates the results for the B-spline interpolants. B-spline and T-spline interpolants are widely used in modeling engineering and science problems which have been described in detail in Refs. [⁴⁶,⁶⁸,⁶⁹]. These interpolants are generalization of NURBS and most commonly used in isogeometric analysis [⁶⁹,⁷⁰,⁷¹]. A B-spline basis function of order p + 1 is a piecewise polynomial function of degree p. In Fig. 8, the convergence rate for B-spline of degree p = 4 in the L^{\mathrm{2}} norm is around 3.6 and in the H^{\mathrm{1}} semi-norm and H^{\mathrm{2}} semi-norm is 3 and 2, respectively. Figures 7 and 8 show the super convergent results with the fourth order phase-filed model for the basis functions that reproduce quadratic polynomials.

Conclusions

We have applied a fourth-order phase-field model in combination with smooth LME and SME approximants and B-spline interpolants. The smoothness and higher-order continuity of the LME and SME approximants, i.e., C^{\infty }, allows to directly solve the fourth-order phase-field equation without splitting it into two second-order differential equations. We compared the second- and the fourth-order phase-field model with smooth approximants. The fourth-order phase-filed model gives more accurate and efficient results with smooth approximants. We conclude that for the LME approximants there is an optimal value of \gamma \le0.8 for the fourth-order phase-field model and \gamma \ge4.8 for the second-order phase-field model that maximizes the accuracy. However, the second-order phase-field model is not appropriate to be used with the SME approximants. While, the fourth-order phase-field model with SME approximants, gives an excellent result. We also conclude that for the SME approximants there is an optimal value of \alpha= 2.5 for the fourth-order phase-field model that maximizes the accuracy. We highlight the good behavior of the first-order LME approximants with \gamma \le0.8 to approximate functionals depending on second-order derivatives. We conclude that the super convergent results for the fourth-order phase-filed model, is obtained with LME, SME and B-spline basis functions. The proposed smooth approximations also show potential for other problems which will be examined in the future, such as cracks in thin shell bodies with complex geometry and topology and Cahn-Hilliard equation.

References

[1] Amiri F, Anitescu C, Arroyo M, Bordas S P A, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57doi:10.1007/BF00052492

[2] Winne D H, Wundt B M. Application of the Griffith-Irwin theory of crack propagation to the bursting behavior of disks, including analytical and experimental studies. Transactions of the American Society of Mechanical Engineers, 1958, 80: 1643–1655

[3] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotations. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[4] Areias P, Rabczuk T, Camanho P P. Finite strain fracture of 2D problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72: 50–63

[5] Areias P, Rabczuk T, Dias da Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[6] Nguyen-Xuan H, Liu G R, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 252–273

[7] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[8] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150

[9] Strouboulis T, Babuska I, Copps K. The design and analysis of the generalized finite element method. Computer Methods in Applied Mechanics and Engineering, 2000, 181(1–3): 43–69

[10] .Ma G, An X, He L. The numerical manifold method: A review. International Journal of Computational Methods, 2010, 7(1): 1–32

[11] Bhardwaj G, Singh I V, Mishra B K, Bui T Q. Numerical simulation of functionally graded cracked plates using {NURBS} based {XIGA} under different loads and boundary conditions. Composite Structures, 2015, 126: 347–359

[12] Ghorashi S Sh, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based {XIGA} for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

[13] Nguyen-Thanh N, Valizadeh N, Nguyen M N, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on Kirchhoff-Love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[14] Amiri F, Millan D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69(2): 102–109

[15] Li S, Liu W K. Meshfree and particle methods and their applications. Applied Mechanics Reviews, 2002, 55(1): 1–34

[16] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23–24): 1391–1411

[17] Wang S, Zhang H. Partition of unity-based thermomechanical meshfree method for two-dimensional crack problems. Archive of Applied Mechanics, 2011, 81(10): 1351–1363

[18] Zhuang X, Cai Y, Augarde C. A meshless sub-region radial point interpolation method for accurate calculation of crack tip fields. Theoretical and Applied Fracture Mechanics, 2014, 69: 118–125

[19] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[20] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[21] de Borst R. Damage, Material Instabilities, and Failure. Hoboken: John Wiley & Sons, 2004

[22] de Borst R, Benallal A, Heeres O M. A gradient-enhanced damage approach to fracture. Journal de Physique Archives IV France, 1996, 6(C6): 491–502

[23] Liang J, Zhang Z, Prevost J H, Suo Z. Time-dependent crack behavior in an integrated structure. International Journal of Fracture, 2004, 125(3–4): 335–348

[24] Peerlings R H J, de Borst R, Brekelmans W A M, De Vree J H P. Gradient enhanced damage for quasi-brittle materials. International Journal for Numerical Methods in Engineering, 1996, 39(19): 3391–3403

[25] Peerlings R H J, de Borst R, Brekelmans W A M, de Vree J H P, Spee I. Some observations on localization in non-local and gradient damage models. European Journal of Mechanics. A, Solids, 1996, 15(6): 937–953

[26] Sluys L J, de Borst R. Dispersive properties of gradient-dependent and rate-dependent media. Mechanics of Materials, 1994, 18(2): 131–149

[27] Bažant Z. Imbricate continuum and its variational derivation. Journal of Engineering Mechanics, 1984, 110(12): 1693–1712

[28] Bažant Z, Chang T P. Instability of nonlocal continuum and strain averaging. Journal of Engineering Mechanics, 1984, 110(10): 1441–1450

[29] Bažant Z, Jirasek M. Nonlocal integral formulations of plasticity and damage: Survey of progress. Journal of Engineering Mechanics, 2002, 128(11): 1119–1149

[30] Bažant Z, Lin F B. Nonlocal smeared cracking model for concrete fracture. Journal of Structural Engineering, 1988, 114(11): 2493–2510

[31] Areias P, Rabczuk T, de Sá J C. A novel two-stage discrete crack method based on the screened Poisson equation and local mesh refinement. Computational Mechanics, 2016, 58(6): 1003–1018

[32] Areias P, Msekh M A, Rabczuk T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158(Suppl C): 116–143

[33] Hofacker M, Miehe C. A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns. International Journal for Numerical Methods in Engineering, 2013, 93(3): 276–301

[34] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273–1311

[35] Mumford D, Shah J. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on Pure and Applied Mathematics, 1989, 42(5): 577–685

[36] Ambrosio L. Variational problems in SBV and image segmentation. Acta Applicandae Mathematicae, 1989, 17(1): 1–40

[37] Ambrosio L, Tortorelli V M. Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence. Communications on Pure and Applied Mathematics, 1990, 43(8): 999–1036

[38] Gómez H, Calo V M, Bazilevs Y, Hughes T J R. Isogeometric analysis of the Cahn-Hilliard phase-field model. Computer Methods in Applied Mechanics and Engineering, 2008, 197(49–50): 4333–4352

[39] Esedoglu S, Shen J. Digital inpainting based on the Mumford-Shah-Euler image model. European Journal of Applied Mathematics, 2002, 13(4): 353–370

[40] Amiri F, Millan D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ent approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 254–275

[41] Bourdin B, Francfort G A, Marigo J J. The variational approach to fracture. Journal of Elasticity, 2008, 91(1–3): 5–148

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

[43] Borden M J, Verhoosel C V, Scott M A, Hughes T J R, Landis C M. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 2012, 217–220: 77–95

[44] Areias P, Rabczuk T, Msekh M A. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312(Suppl C): 322–350

[45] Kiendl J, Ambati M, De Lorenzis L, Gomez H, Reali A. Phase-field description of brittle fracture in plates and shells. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 374–394

[46] Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

[47] Sukumar N. Construction of polygonal interpolants: A maximum entropy approach. International Journal for Numerical Methods in Engineering, 2004, 61(12): 2159–2181

[48] Shannon C E. A mathematical theory of communication. Sigmobile Mobile Computing & Communications Review, 2001, 5(1): 3–55

[49] Jaynes E T. Information theory and statistical mechanics. Physical Review, 1957, 106(4): 620–630

[50] Jaynes E T. Information theory and statistical mechanics. II. Physical Review, 1957, 108(2): 171–190

[51] Arroyo M, Ortiz M. Local maximum-entropy approximation schemes: A seamless bridge between finite elements and meshfree methods. International Journal for Numerical Methods in Engineering, 2006, 65(13): 2167–2202

[52] Shepard D. A two-dimensional interpolation function for irregularly-spaced data. In: Proceedings of the 1968 23rd ACM National Conference. New York: ACM, 1968, 517–524

[53] Rosolen A, Millan D, Arroyo M. Second-order convex maximum entropy approximants with applications to high-order PDE. International Journal for Numerical Methods in Engineering, 2013, 94(2): 150–182

[54] Cyron C J, Arroyo M, Ortiz M. Smooth, second order, non-negative meshfree approximants selected by maximum entropy. International Journal for Numerical Methods in Engineering, 2009, 79(13): 1605–1632

[55] Ortiz A, Puso M A, Sukumar N. Maximum-entropy meshfree method for incompressible media problems. Finite Elements in Analysis and Design, 2011, 47(6): 572–585

[56] Millán D, Rosolen A, Arroyo M. Thin shell analysis from scattered points with maximum-entropy approximants. International Journal for Numerical Methods in Engineering, 2011, 85(6): 723–751

[57] Millán D, Rosolen A, Arroyo M. Nonlinear manifold learning for meshfree finite deformations thin shell analysis. International Journal for Numerical Methods in Engineering, 2013, 93(7): 685–713

[58] Ortiz A, Puso M A, Sukumar N. Maximum-entropy meshfree method for compressible and near-incompressible elasticity. Computer Methods in Applied Mechanics and Engineering, 2010, 199(25–28): 1859–1871

[59] Belytschko T, Lu Y Y, Gu L. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[60] Rosolen A, Millan D, Arroyo M. On the optimum support size in meshfree methods: A variational adaptivity approach with maximum entropy approximants. International Journal for Numerical Methods in Engineering, 2010, 82(7): 868–895

[61] Shannon C. A mathematical theory of communication. Bell System Technical Journal, 1948, 27(4): 623–656

[62] Sukumar N, Moran B, Belytschko T. The natural element method in solid mechanics. International Journal for Numerical Methods in Engineering, 1998, 43(5): 839–887

[63] Cirak F, Ortiz M, Schröder P. Subdivision surfaces: A new paradigm for thin-shell finite-element analysis. International Journal for Numerical Methods in Engineering, 2000, 47(12): 2039–2072

[64] Rajan V T. Optimality of the Delaunay triangulation in R^d. Discrete & Computational Geometry, 1994, 12(2): 189–202

[65] Borden M J, Hughes T J R, Landis C M, Verhoosel C V. A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering, 2014, 273(5): 100–118

[66] Apostolatos A, Schmidt R, Wüchner R, Bletzinger K U. A Nitsche-type formulation and comparison of the most common domain decomposition methods in isogeometric analysis. International Journal for Numerical Methods in Engineering, 2014, 97(7): 473–504

[67] Hughes T J R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. New York: Dover Publications, 2000

[68] Cottrell J A, Hughes T J R, Bazilevs Y. Isogeometric Analysis: Toward Integration of CAD and FEA. Chichester: Wiley, 2009

[69] Piegl L, Tiller W. The NURBS Book. 2nd ed. New York: Springer, 1997

[70] Bazilevs Y, Calo V M, Cottrell J A, Evans J A, Hughes T J R, Lipton S, Scott M A, Sederberg T W. Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5–8): 229–263

[71] Rogers D F. An Introduction to NURBS With Historical Perspective. San Diego: Academic Press, 2001

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature