A phase-field method for modeling cracks with frictional contact

TL;DR

This paper presents a novel phase-field method for modeling cracks with frictional contact that simplifies the representation of complex crack geometries and contact constraints without explicit algorithms.

Contribution

It introduces a new phase-field approach that captures arbitrary crack geometries and contact conditions without requiring explicit contact algorithms.

Findings

Successfully models crack growth with frictional contact.

Verifies the method against existing discrete models.

Demonstrates capability for complex crack and contact scenarios.

Abstract

We introduce a phase-field method for continuous modeling of cracks with frictional contacts. Compared with standard discrete methods for frictional contacts, the phase-field method has two attractive features: (1) it can represent arbitrary crack geometry without an explicit function or basis enrichment, and (2) it does not require an algorithm for imposing contact constraints. The first feature, which is common in phase-field models of fracture, is attained by regularizing a sharp interface geometry using a surface density functional. The second feature, which is a unique advantage for contact problems, is achieved by a new approach that calculates the stress tensor in the regularized interface region depending on the contact condition of the interface. Particularly, under a slip condition, this approach updates stress components in the slip direction using a standard contact constitutive law, while making other stress components compatible with stress in the bulk region to ensure non-penetrating deformation in other directions. We verify the proposed phase-field method using stationary interface problems simulated by discrete methods in the literature. Subsequently, by allowing the phase field to evolve according to brittle fracture theory, we demonstrate the proposed method's capability for modeling crack growth with frictional contact.

Equations

$$\gamma (d,\nabla d):=\frac{1}{2}\left(\frac{d^2}{L}+L\nabla d\cdot \nabla d\right).$$

$$\nabla \cdot \mathbf{\sigma }(\mathbf{\epsilon })=0\text{in}\Omega ,$$

$$\mathbf{u}=\hat{\mathbf{u}}$$

$$\mathbf{\upsilon }\cdot \mathbf{\sigma }=\hat{\mathbf{t}}$$

$$g(d)\in [0,1],g(0)=1,g(1)=0,g^{\prime }(d)<0.$$

$$g(d)=(1-d{)}^2,$$

$$\mathbf{\sigma }=g(d){\mathbf{\sigma }}_{\text{bulk}}+[1-g(d)]{\mathbf{\sigma }}_{\text{interface}}.$$

$$g^{\prime }(d)\mathcal{H}+G_c\left(\frac{d}{L}-L\nabla \cdot \nabla d\right)=0\text{in}\Omega ,$$

$${\mathbf{\sigma }}_{\text{bulk}}={\mathbb{C}}_{\text{bulk}}:\mathbf{\epsilon },$$

$${\mathbb{C}}_{\text{bulk}}={\mathbb{C}}^{\mathrm{e}}:=\lambda 1\otimes 1+2G\mathbb{I}.$$

$${\epsilon }_{\mathrm{N}}$$

$$\mathbf{\sigma }=g(d){\mathbf{\sigma }}_{\text{bulk}}.$$

$$f:=|\tau |-\mu p_{\mathrm{N}}\le 0.$$

$$\tau \equiv {\sigma }_{nm}=\mathbf{\sigma }:(\mathbf{n}\otimes \mathbf{m}).$$

$$\mathbf{\sigma }={\mathbf{\sigma }}_{\text{bulk}}.$$

$${\mathbf{\sigma }}_{\text{interface}}={\mathbf{\sigma }}_{\text{friction}}+{\mathbf{\sigma }}_{\text{no-penetration}}.$$

$${\mathbf{\sigma }}_{\text{no-penetration}}={\mathbf{\sigma }}_{\text{bulk}}-{\tau }_{\text{bulk}}(\mathbf{n}\otimes \mathbf{m}+\mathbf{m}\otimes \mathbf{n}),$$

$${\mathbf{\sigma }}_{\text{friction}}=\mu p_{\mathrm{N},\text{bulk}}\mathrm{sign}({\tau }_{\text{bulk}})(\mathbf{n}\otimes \mathbf{m}+\mathbf{m}\otimes \mathbf{n}).$$

$${\mathbf{\sigma }}_{\text{interface}}={\mathbf{\sigma }}_{\text{bulk}}+(\mu p_{\mathrm{N},\text{bulk}}\mathrm{sign}({\tau }_{\text{bulk}})-{\tau }_{\text{bulk}})(\mathbf{n}\otimes \mathbf{m}+\mathbf{m}\otimes \mathbf{n}).$$

$$\mathbf{\sigma }={\mathbf{\sigma }}_{\text{bulk}}+[1-g(d)](\mu p_{\mathrm{N},\text{bulk}}\mathrm{sign}({\tau }_{\text{bulk}})-{\tau }_{\text{bulk}})(\mathbf{n}\otimes \mathbf{m}+\mathbf{m}\otimes \mathbf{n}).$$

$${\mathcal{S}}_u$$

$${\mathcal{S}}_d$$

$${\mathcal{V}}_u$$

$${\mathcal{V}}_d$$

$$-{\int }_{\Omega }{\nabla }^{\mathrm{s}}\mathbf{\eta }:\mathbf{\sigma }\mathrm{d}V+{\int }_{{\partial }_t\Omega }\mathbf{\eta }\cdot \hat{\mathbf{t}}\mathrm{d}A=0,$$

$${\int }_{\Omega }\varphi g^{\prime }(d)\mathcal{H}\mathrm{d}V+{\int }_{\Omega }{G}_c\left(\varphi \frac{d}{L}+L\nabla \varphi \cdot \nabla d\right)\mathrm{d}V=0.$$

$$\mathcal{R}=-\mathcal{J}\mathrm{\Delta }\mathbf{U},$$

$$[\mathcal{R}{]}_e^i:=-{\int }_{{\Omega }_e}{\nabla }^{\mathrm{s}}{\mathbf{\eta }}^i:\mathbf{\sigma }\mathrm{d}V+{\int }_{{\partial }_t{\Omega }_e}{\mathbf{\eta }}^i\cdot \hat{\mathbf{t}}\mathrm{d}A,$$

$$[\mathcal{J}{]}_e^{i,j}:=-{\int }_{{\Omega }_e}{\nabla }^{\mathrm{s}}{\mathbf{\eta }}^i:\mathbb{C}:{\nabla }^{\mathrm{s}}{\mathbf{\eta }}^j\mathrm{d}V,$$

$$\displaystyle\mathbb{C}=\left\{\begin{array}[]{ll}g(d)\mathbb{C}_{\text{bulk}}&\text{for a no-contact condition,}\\ \mathbb{C}_{\text{bulk}}&\text{for a stick condition,}\\ \mathbb{C}_{\text{bulk}}+[1-g(d)](\mathbb{C}_{f}-\mathbb{C}_{\tau})&\text{for a slip condition.}\end{array}\right.$$

$${\mathbb{C}}_f$$

$${\mathbb{C}}_{\tau }$$

Full text

††journal:  

1 Introduction

Frictional cracks are ubiquitous in natural and manufactured systems. For example, in the Earth’s crust, frictional cracks appear over a wide range of scales . They are also widespread in many branches of science and engineering including material sciences and civil and mechanical engineering. Accordingly, the numerical modeling of motion and friction in crack surfaces has long been an important subject, and there is a large body of literature on computational contact mechanics (see, for example, Laursen [10], Wriggers [11], and references therein).

At present, standard numerical methods treat frictional cracks as discrete discontinuities subjected to constraints on contact conditions. The discontinuities should be aligned with element boundaries in classical finite element methods, while they can be embedded inside elements in modern methods Irrespective of their alignment with elements, the contact surfaces must satisfy a set of constraints including the no-penetration constraint under compression. Imposing these constraints on discrete interfaces, however, is an outstanding challenge in computational contact mechanics. A large number of studies have addressed this challenge by employing various types of algorithms such as Nevertheless, the optimal way to treat these contact constraints is yet an unresolved issue. Also importantly, numerical methods for frictional contacts require significant effort for implementation, especially when one wants to accommodate complex interface geometry. For these reasons, a simple numerical method is desired that can model frictional contact problems with low implementation cost.

In this paper, we propose a phase-field method for frictional crack problems that can efficiently handle complex crack geometry and contact conditions. Phase-field modeling is a continuous (as opposed to discrete) approach to interface problems that approximates a sharp interface as a region where the phase field attains a certain value. Its upshot is that one can represent an interface without any function or algorithm for describing its geometry, which is highly advantageous when the geometry is complex and may evolve with time. For this reason, phase-field modeling has found widespread applications in a variety of scientific and engineering problems. Recently, it has enjoyed considerable success in computational modeling of a wide range of fracturing processes Nevertheless, to our knowledge, the present work is the first attempt to apply a phase-field approach to cracks with frictional contact.

The key idea of the proposed phase-field method is to incorporate contact behaviors and constraints through suitable calculation of the stress tensor in the regularized interface region. In existing phase-field models of fracture, the stress tensor in the interface region is either degraded or maintained according to the sign of some part of the strain tensor. This way roughly considers a no-contact condition and a stick contact condition, but these two conditions are usually not distinguished in a manner consistent with contact mechanics. Furthermore, and perhaps more importantly, phase-field fracture models have not incorporated a slip contact condition in which relative motion between interacting surfaces takes places according to friction. Note that a slip condition is a major challenge in computational contact mechanics because it requires one to model slip behavior while imposing the no-penetration contact in non-slip directions. In this work, we propose a new approach that incorporates and distinguishes between all contact conditions by a proper calculation of the stress tensor in the interface region. Our use of stress tensor is consistent with other types of smeared crack formulations for frictional cracks (e.g. Borja [45]), but our way to calculate stress is completely different from existing smeared methods as it builds on a stress calculation procedure in phase-field modeling of fracture.

Compared with standard discrete methods for frictional contacts, the proposed phase-field method has two attractive features. First, it can represent arbitrary interface geometry without an explicit function or enriched basis functions, which is indeed a hallmark of all phase-field methods. Second, it can accommodate contact constraints without a dedicated algorithm, which is a unique advantage for contact problems. This new feature is attained by making the components of the stress tensor in the non-slip directions compatible between the interface and bulk regions. Notably, this way is a modification of the volumetric–deviatoric decomposition approach proposed by Amor et al. [46] for considering the no-penetration constraint in unilateral frictionless contact, The rest components of the stress tensor, which are relevant to stick/slip behavior, are determined using a standard constitutive law for frictional cracks. In this way, the proposed phase-field method translates a discrete problem with constraints into a continuous problem with multiple constitutive responses. As the continuous problem can be solved by the standard finite element method, it would require significantly less effort for implementation as compared with other methods that can simulate frictional cracks passing through the interior of elements.

The paper is organized as follows. In Section 2, we develop a phase-field formulation for a boundary-value problem that involves frictional contact. Following a standard diffuse approximation of crack geometry in phase-field modeling of fracture, we introduce a new approach that explicitly considers and calculates the stress tensor in the interface region according to the contact condition of the interface. In Section 3, we discretize the phase-field contact formulation using the standard finite element method. In doing so, we present algorithms for calculating the stress, stress–strain tangent, and unit normal and slip vectors at quadrature points in the regularized interface region. In Section 4, we verify the proposed phase-field method using stationary interface problems that have been simulated by discrete methods in the literature. We then combine the proposed method with an evolution equation for brittle fracture, and demonstrate the method’s capability for modeling crack growth with frictional contact. In Section 5, we conclude the work.

2 Phase-field formulation for cracks with frictional contact

In this section, we develop a phase-field formulation for continuous modeling of frictional cracks in solids. The formulation builds on methods originally developed for phase-field modeling of crack propagation, but it can be useful for modeling general frictional interfaces in solids. For this reason, we will focus on the use of a phase-field approach to geometric approximation of frictional contact problems, without delving into aspects of fracture mechanics. Also, to make the following presentation simple, we will assume without loss of generality that the material is free of inertial and body forces, isotropic, elastic, and geometrically and materially linear. If necessary, these assumptions may be relaxed in standard ways in solid mechanics.

2.1 Problem statement and phase-field approximation

Consider a domain $\Omega\in\mathbb{R}^{\mathrm{dim}}$ in a $\mathrm{dim}$-dimensional space with external boundary $\partial\Omega$. The boundary is partitioned into the displacement (Dirichlet) boundary, $\partial_{u}\Omega$, and the traction (Neumann) boundary, $\partial_{t}\Omega$, such that $\overline{\partial_{u}\Omega\cap\partial_{t}\Omega}=\emptyset$ and $\overline{\partial_{u}\Omega\cup\partial_{t}\Omega}=\partial\Omega$. The domain has a set of internal discontinuities, which is denoted by $\Gamma$. A discontinuity has two surfaces that are either separated or in contact. When the two surfaces are in contact, relative motion may or may not exist between them depending on the magnitudes of tractions and frictions therein.

We begin our formulation by approximating the discontinuities’ geometry using a standard approach in phase-field modeling of fracture. Let us define the phase-field variable, $d\in[0,1]$, such that it denotes a fully discontinuous (interface) region by $d=1$ and a fully continuous (bulk) region by $d=0$. We then introduce a surface density functional for length regularization of the sharp geometry of $\Gamma$. Among several forms of the functional proposed in the literature, here we adopt the most popular one, given by

[TABLE]

Here, $L$ is the length parameter introduced for regularization of sharp geometry, which determines the size of the diffuse approximation zone. Figure 1 illustrates how this phase-field approach approximates the original domain with sharp discontinuity. Note that the diffuse approximation naturally gives rise to regions in which $0<d<1$.

Once discontinuities have been diffusely approximated in the way described above, a continuous version of the problem can be stated as follows. Find the displacement field in this domain, $\boldsymbol{u}$, that satisfies the balance of linear momentum

[TABLE]

where $\boldsymbol{\sigma}$ is the (Cauchy) stress tensor and $\boldsymbol{\varepsilon}$ is the infinitesimal strain tensor defined as the symmetric gradient of $\boldsymbol{u}$. The boundary conditions of this problem are given by

[TABLE]

where $\hat{\boldsymbol{u}}$ and $\hat{\boldsymbol{t}}$ are prescribed boundary conditions of displacement and traction vectors, respectively, and $\boldsymbol{\upsilon}$ is the outward unit normal vector at the domain boundary. Note that no boundary condition is imposed on $\Gamma$ as the discontinuities have already been smeared in the domain $\Omega$ through the phase-field approximation described above.

As the whole domain is now regarded as a continuum, stress tensors in the bulk and interface systems should be continuously interpolated. For this interpolation, we introduce a function of the phase field, $g(d)$, that satisfies

[TABLE]

This function is usually called the degradation function in phase-field modeling of fracture. In this work, we use the most common form of $g(d)$ in the literature, given by

[TABLE]

Using this function, we can express the stress tensor in the domain as

[TABLE]

where $\boldsymbol{\sigma}_{\text{bulk}}$ and $\boldsymbol{\sigma}_{\text{interface}}$ are stress tensors in the bulk and interface systems, respectively. It is noted that Eq. (7) is a generalization of the way in which stress is calculated in phase-field models of fracture.

In addition, the phase-field variable is postulated to satisfy the following partial differential equation

[TABLE]

which is also adopted from phase-field modeling of fracture. Here, $\mathcal{H}$ and $G_{c}$ are positive parameters corresponding to the crack driving force and the critical fracture energy in the context of phase-field modeling of fracture. As such, for simulation of growing cracks, they may be calculated according to a phase-field formulation for fracture. However, if the crack interface is assumed to be stationary, one may take any positive values for these two parameters for initialization of the phase field. In this case, Eq. (8) is solved only once in the beginning of the problem to initialize the phase-field variable, and the phase field remains constant throughout the course of loading.

2.2 Calculation of stress tensors according to contact conditions

So far, the only difference between our formulation and the most standard phase-field formulation for fracture is that here we have explicitly considered the stress tensor in the interface system, $\boldsymbol{\sigma}_{\text{interface}}$. This modification has been made to incorporate contact-dependent mechanical responses of the interface system into the phase-field formulation. In the following, we propose a specific procedure for calculating the bulk and interface stress tensors in Eq. (7) according to the contact condition of the interface.

First, we calculate the bulk stress through a standard stress–strain relationship in continuum mechanics, namely

[TABLE]

where $\mathbb{C}_{\text{bulk}}$ is the fourth-order stress–strain tangent tensor of the bulk region. For a linear elastic material, the stress–strain tangent is given by

[TABLE]

Here, $\lambda$ and $G$ are the Lamé parameters which can be converted into Young’s modulus $E$ and Poisson’s ratio $\nu$, $\boldsymbol{1}$ is the second-order identity tensor, and $\mathbb{I}$ is the fourth-order symmetric identity tensor. In short, the bulk stress tensor is evaluated as usual.

Next, we propose a new way to calculate the interface stress tensor depending on the contact condition of the interface system. To identify the contact condition, we introduce a coordinate system that is oriented with respect to the interface normal and tangential directions. Figure 2 depicts this interface-oriented coordinate system in a 2D domain. Hereafter, we denote by $\boldsymbol{n}$ the unit vector in the interface normal direction and by $\boldsymbol{m}$ the unit vector in the slip direction. These unit vectors are assumed to be known for now.

In the interface-oriented coordinate system, the normal strain along the interface normal direction is calculated as

[TABLE]

This strain can be used to distinguish between contact and no-contact conditions. By definition, this strain plays the role of the gap function in classical contact mechanics.

When $\varepsilon_{\mathrm{N}}>0$, the interface has a gap between its two surfaces, which corresponds to an open (non-contacting) crack in phase-field modeling of fracture. In this case, the interface system is stress-free, i.e. $\boldsymbol{\sigma}_{\text{interface}}=\boldsymbol{0}$. We thus evaluate the stress tensor under this no-contact condition as

[TABLE]

It is noted that the above expression is the same as the stress equation for an open crack in phase-field models of fracture.

By contrast, when $\varepsilon_{\mathrm{N}}\leq 0$, the interface is considered being in contact. The contact condition of a cohesive–frictional interface is either a stick condition or a slip condition. The distinction between stick and slip conditions can be made by the following yield function

[TABLE]

Here, $\mu$ is the friction coefficient of the interface, and $p_{\mathrm{N}}$ is the contact normal pressure, defined as $p_{\mathrm{N}}\equiv-\sigma_{nn}=-\boldsymbol{\sigma}:(\boldsymbol{n}\operatorname{\otimes}\boldsymbol{n})$ in accordance to Eq. (11). Note that $p_{\mathrm{N}}\geq 0$ whenever $\varepsilon_{\mathrm{N}}\leq 0$. Lastly, $\tau$ is the resolved shear stress in the interface, which can be calculated in this continuum formulation as (see Fig. 2)

[TABLE]

Note that $\sigma_{nm}=\sigma_{mn}$ by the symmetry of the stress tensor. A contacting interface is under a stick condition when $f<0$, whereas it is under a slip condition when $f=0$. Therefore, stick and slip conditions are distinguished in the same way in classical contact mechanics.

Under a stick condition, no relative motion exists between the bulk and interface systems, so $\boldsymbol{\sigma}_{\text{bulk}}=\boldsymbol{\sigma}_{\text{interface}}$. Therefore, the stress tensor under a stick condition can be calculated as

[TABLE]

We note that this corresponds to the standard way how a phase-field model of fracture treats a closed crack under compression. However, existing phase-field models do not allow slip motion along the interface. This limitation is the major motivation of this work and tackled in the following.

When the interface is under a slip condition, the stress tensor in the interface system is non-zero and different from the bulk stress. We thus have to evaluate the interface stress such that it incorporates the frictional contact behavior and the no-penetration constraint simultaneously. For this purpose, we decompose the interface stress tensor into a friction part, $\boldsymbol{\sigma}_{\text{friction}}$, and a no-penetration part, $\boldsymbol{\sigma}_{\text{no-penetration}}$, as

[TABLE]

Of these two, the no-penetration part is determined to make the deformation along the interface and bulk systems compatible in all directions except the slip direction. This compatibility can be attained by defining the no-penetration part as follows:

[TABLE]

where $\tau_{\text{bulk}}:=\boldsymbol{\sigma}_{\text{bulk}}:(\boldsymbol{n}\operatorname{\otimes}\boldsymbol{m})$. In words, the no-penetration part is computed by fully degrading the slip direction components of the bulk stress tensor. Then, the friction part is calculated to substitute the degraded components, according to a prescribed contact constitutive law. Under a slip condition, $|\tau|=\mu p_{\mathrm{N}}$ because the yield function (13) is zero, and here $p_{\mathrm{N}}$ must be equal to $p_{\mathrm{N},\text{bulk}}:=-\boldsymbol{\sigma}_{\text{bulk}}:(\boldsymbol{n}\operatorname{\otimes}\boldsymbol{n})$ because the bulk stress is not degraded in the contact normal direction. Also, it is likely that the signs of $\tau$ and $\tau_{\text{bulk}}$ are identical. Therefore, we can express the friction part of the stress tensor as

[TABLE]

Inserting Eqs. (18) and (17) into Eq. (16), we obtain the interface stress tensor as

[TABLE]

One can see that this stress tensor is obtained by replacing the slip-relevant part of the bulk stress tensor with $\boldsymbol{\sigma}_{\text{friction}}$. Substituting the above equation into Eq. (7) gives the overall stress tensor under a slip condition as

[TABLE]

Note that for all contact conditions, we get $\boldsymbol{\sigma}=\boldsymbol{\sigma}_{\text{interface}}$ when $d=1$ and $\boldsymbol{\sigma}=\boldsymbol{\sigma}_{\text{bulk}}$ when $d=0$.

To summarize, we have proposed an approach that calculates the stress tensor in the interface system according to the contact condition of the interface. In this approach, the interface stress is null under a no-contact condition and equal to the bulk stress under a stick contact condition. The interface stress under a slip condition is calculated as a combination of the friction part and the no-penetration part so that the contact constitutive behavior and the no-penetration constraints are incorporated into the phase-field formulation.

Remark 1*.*

The foregoing expressions for the interface stress tensor can be re-interpreted based on the decomposition of $\boldsymbol{\sigma}_{\text{interface}}=\boldsymbol{\sigma}_{\text{friction}}+\boldsymbol{\sigma}_{\text{no-penetration}}$. The friction part, $\boldsymbol{\sigma}_{\text{friction}}$, is zero under a no-contact condition, compatible with the bulk stress under a stick condition, and calculated from a friction constitutive law under a slip condition. The no-penetration part, $\boldsymbol{\sigma}_{\text{no-penetration}}$, is zero under a no-contact condition while it is compatible with the bulk stress under stick and slip contact conditions.

3 Discretization and algorithms

This section describes discretization methods and algorithms for numerical solution of the proposed formulation using the standard finite element method.

3.1 Finite element discretization

The proposed phase-field formulation can be readily solved by the standard finite element method. For finite element discretization, we first define trial solution spaces for the displacement field and the phase field as

[TABLE]

where $H^{1}$ denotes a Sobolev space of order one. Weighting function spaces for the two fields are accordingly defined as

[TABLE]

Applying the standard weighted residual procedure, we obtain the following two variational equations:

[TABLE]

Here, Eq. (25) is the linear momentum balance equation and Eq. (26) is the phase-field equation. Both of them are solved in each load step if the interface system is subjected to growth during the course of loading. However, for a stationary interface problem, the phase-field equation (26) only needs to be solved once in the initialization stage of the problem. The Galerkin and matrix forms of these equations can be developed in a standard manner, so they are omitted for brevity.

As we are considering linear elasticity in this work, the momentum balance equation (25) is linear if the contact condition inside the domain is fixed. This is because the phase-field method has formulated a frictional crack problem as a continuum problem with stiffness that may spatially vary according to the phase field. However, when the contact condition of a point is subject to change after loading, the problem is incrementally nonlinear. Therefore, it is necessary to apply a nonlinear solution method for the momentum equation. Note that the solution method will converge quickly whenever the contact condition remains unchanged from an initial guess. The phase-field equation (26) is always linear so it can be solved easily.

In this work, we use Newton’s method to solve the discretized momentum balance equation. During Newton iterations, the increment of the nodal displacement vector, denoted by $\mathrm{\Delta}{\boldsymbol{U}}$, can be obtained by solving

[TABLE]

where $\boldsymbol{\mathcal{R}}$ is the residual vector, of which element-wise contribution can be calculated as

[TABLE]

and $\boldsymbol{\mathcal{J}}$ is the Jacobian matrix, of which element-wise contribution can be calculated as

[TABLE]

with $e$ denoting an element index and $i,j$ denoting shape function indices. $\mathbb{C}$ is the stress–strain tangent that is the same as $\mathbb{C}_{\text{bulk}}$ in the bulk region but may be different from it otherwise.

Remark 2*.*

3.2 Update of stress and tangent tensors

To evaluate Eqs. (28) and (29) during finite element assembly, we need to calculate the stress tensor, $\boldsymbol{\sigma}$, and the stress–strain tangent tensor, $\mathbb{C}$, at every quadrature point. Consider a typical Newton update step for which the strain tensor, $\boldsymbol{\varepsilon}$, and the phase-field variable, $d$, are given at the quadrature point of interest. It is also assumed that the unit vector in the interface normal direction, $\boldsymbol{n}$ and the unit vector along the slip direction, $\boldsymbol{m}$, are also known at this point.

The stress tensor at a quadrature point can be updated as described in Algorithm 1. The algorithm first checks whether the current quadrature point belongs to a bulk region where $d=0$. If not, the algorithm identifies the contact condition at the quadrature point and then updates the stress tensor following the approach proposed in Section 2.2. To distinguish between stick and slip conditions, the algorithm uses a standard predictor–corrector approach employing the bulk stress as a trial stress. Therefore, the yield function, $f$, is calculated using $\tau$ and $p_{\mathrm{N}}$ of the bulk stress.

The stress–strain tangent tensor should also be evaluated to assemble the Jacobian matrix, see Eq. (29). This calculation is trivial for a bulk region. For an interface region where $0<d\leq 1$, it is given by

[TABLE]

Here, for a slip condition, $\mathbb{C}_{f}$ is defined as

[TABLE]

and $\mathbb{C}_{\tau}$ is defined as

[TABLE]

Note that all these tensors can be calculated in a straightforward manner as long as the vectors $\boldsymbol{n}$ and $\boldsymbol{m}$ are given.

3.3 Calculation of unit normal and slip vectors

As described above, the unit vector in the interface normal direction, $\boldsymbol{n}$, and the unit vector in the slip direction, $\boldsymbol{m}$, are crucial for the proposed phase-field formulation. Usually, in phase-field modeling, the interface normal vector is approximated as $\boldsymbol{n}\approx\operatorname{\nabla}{d}/\|\operatorname{\nabla}{d}\|$. The accuracy of this approximation, however, seems to be insufficient for our purpose, for two main reasons: (1) because a crack tip has been approximated bluntly, $\operatorname{\nabla}{d}/\|\operatorname{\nabla}{d}\|$ calculated around the crack tip region is indeed nearly orthogonal to the desired $\boldsymbol{n}$, and (2) unless $d$ is very close to 1, $\operatorname{\nabla}{d}/\|\operatorname{\nabla}{d}\|$ may not be close enough to the desired $\boldsymbol{n}$.

Therefore, for a more accurate calculation of $\boldsymbol{n}$ and $\boldsymbol{m}$, we devise a new algorithm that first identifies a lower-dimensional crack path from phase-field values and then estimates the unit vectors from the identified crack path. Hereafter, we will restrict our attention to 2D problems because identifying a crack path in a 3D phase field is an outstanding challenge. In a 2D phase-field problem, a crack path would be a 1D line that connects points where $d=1$. Drawing on this idea, we construct such a line through the procedure described in Algorithm 2. In essence, this algorithm finds nodes where $d\approx 1$ that are distant at least the length parameter $L$ and connect them in a piecewise linear manner to approximate the crack path. From this piecewise linearly approximated crack path, we can compute $\boldsymbol{n}$ and $\boldsymbol{m}$ and then assign them to nearby quadrature points.

The proposed algorithm for calculating $\boldsymbol{n}$ and $\boldsymbol{m}$ is simple and appears to be sufficiently accurate for phase-field modeling of frictional interfaces. We note, however, that any other algorithm can be used for the same purpose as long as it gives reliable results for the unit vectors. For example, Ziaei-Rad et al. [49] have proposed a variational method for identifying a crack path in phase-field modeling. The use of such an advanced method will likely improve the accuracy of the overall numerical solution, although it requires significantly more effort for implementation. Furthermore, due to lack of a good algorithm for estimating $\boldsymbol{n}$ and $\boldsymbol{m}$ in 3D, applications of the proposed phase-field formulation will be limited to 2D problems in the next section. Overcoming this limitation for 3D problems will be a future research topic.

4 Numerical examples

This section has two objectives: (1) to verify the proposed phase-field formulation for frictional interfaces, and (2) to demonstrate the capability of the proposed method for modeling crack growth with frictional contact. For the first objective, we adopt three numerical examples of frictional interfaces that have been simulated by discrete methods in the literature. Yet the interfaces in these benchmark examples are stationary (i.e. not allowed to advance). Therefore, for the second objective, we introduce a fourth example whereby a preexisting crack propagates according to the phase-field equation (26).

Results in this section have been obtained using the deal.II finite element library [50, 51]. Bilinear quadrilateral elements have been used for all numerical examples. Plane strain conditions are assumed throughout.

4.1 Square domain with an internal crack

Our first example is compression of an internally cracked domain depicted in Fig. 3. This problem was initially presented in Dolbow et al. [15] and then used by other studies including Liu and Borja [16] and Annavarapu et al. [31]. Here we also consider this problem for verification of the phase-field formulation. The domain is a 1 m wide square and possesses a crack whose tips are located at coordinates $(0.3,0.33)$ m and $(0.7,0.68)$ m. Note that these tip locations are adopted from Annavarapu et al. [31] and slightly different from those in Liu and Borja [16]. We assign the elasticity parameters of the material as $E=10000$ MPa and $\nu=0.3$, and the friction coefficient of the crack as $\mu=0.1$. The top boundary of the domain is subjected to a prescribed displacement of $-0.1$ m (downward).

To investigate the convergence of numerical solution with the length parameter and the element size, we consider three values of $L$, namely $0.008$ m, $0.004$ m, and $0.002$ m and three values of $L/h$, namely 4, 8, and 16. The values of $L/h$ are selected based on results in Borden et al. [33] that indicate $L/h\geq 4$ gives reasonably accurate solutions when the same type of surface density functional is employed. Because such fine discretization is necessary only for the interface region and its nearby, we locally refine elements around a node where the phase-field variable is greater than a threshold, until their size reaches a prescribed $L/h$ value. To determine the threshold value, we recall that the spatial variation of the phase field in the chosen surface density functional is given by $d=\exp(-|x|/L)$, with $x$ denoting the distance from the point where $d=1$ [32]. As $d=\exp(-1)\approx 0.378$ when $x=L$, we set the threshold as 0.1 to make the locally refined region sufficiently wide. The same mesh refinement scheme will be used throughout this section.

For initialization of the phase-field variable, we adopt a standard way in phase-field fracture modeling that prescribes $\mathcal{H}$ at quadrature points around a preexisting crack (see Appendix A of Borden et al. [33] for example). With $\mathcal{H}$ values prescribed to make $d=1$ at the initial crack, we solve the phase-field equation (26) once to obtain a phase-field distribution that will be used throughout the problem. Figure 4 shows phase-field distributions in the $L=0.008$ m, $0.004$ m, and $0.002$ m cases when $L/h=8$. It is clear that the diffuse approximation zone becomes narrower as $L$ decreases. After initializing the phase field in this way, we simulate the problem through 10 load steps with a uniform displacement increment of -0.01 m on the top boundary.

We begin by checking to see whether numerical solutions converge with mesh refinement, Although not presented, results of mesh refinement studies conducted with other length parameters showed more or less the same patterns. Therefore, we have found that the numerical model converges with mesh refinement, and that refinement levels of $L/h\geq 4$ provide reasonable accuracy.

Next, we examine how numerical solutions are sensitive to the length parameter for phase-field approximation.

One can see that the results show little sensitivity to $L$. Although a smaller $L$ leads to a sharper displacement jump across the crack and thereby larger contact stresses, the $L=0.008$ m case also shows fairly good results. This observation indicates that a rather diffuse approximation still can provide reasonably good solutions. Therefore, it can be concluded that the phase-field method does not require a very small $L$ for practical purposes while allowing us to obtain a more accurate solution by reducing $L$.

Having confirmed that the proposed method gives consistent solutions, we now verify it with results in the literature. We particularly compare our results from the $L=0.002$ m and $L/h=8$ case with results in Annavarapu et al. [31] obtained by the combination of XFEM and a weighted Nitsche method. Figure 9 shows this comparison. As can be seen, the phase-field and XFEM results are nearly identical in both qualitative and quantitative aspects. Therefore, we have verified that the proposed phase-field method can provide numerical solutions comparable to those obtained by advanced discrete methods for frictional cracks.

Lastly, we remark that the phase-field formulation for this problem is linear in all the load steps, requiring only a single Newton update. This is because the contact condition along the crack is identified as a slip condition from the initial stress-free condition (as zero stress makes $f=0$), and it is indeed a slip condition throughout the problem. So this example was nothing but a linear elasticity problem with heterogeneous stiffness. This means that, although the phase-field method requires a quite fine mesh, its fast convergence can counterbalance the overall computational cost.

4.2 Square domain with an inclined interface

The purpose of our second example is to investigate the ability of the proposed phase-field method to distinguish between stick and slip conditions. For this purpose, we adopt the problem of a square domain with an inclined interface, which was also first used in Dolbow et al. [15] and later revisited by Annavarapu et al. [31], among others. The setup of this problem is illustrated in Fig. 10. Similar to the previous example, a 1 m wide square domain is compressed from the top, but here the discontinuous interface is extended to the side boundaries of the domain. The interface is inclined from the horizontal with an angle $\theta=\tan^{-1}(0.2)$. Therefore, when the friction coefficient $\mu$ is smaller than 0.2, the upper block should slip along the interface; otherwise, the upper and lower blocks should be sticked together and behave as a whole. Accordingly, this problem serves as a good benchmark example for examining the capability for distinguishing between stick and slip behaviors.

Because the domain size remains the same as the previous example, we consider the same three length parameters, $L=0.008$ m, 0.004 m, and 0.002 m. We discretize the domain and initialize the phase field using the same way in the previous example. The refinement level is now fixed as $L/h=8$. Following Annavarapu et al. [31], we consider two cases of friction coefficients, namely $\mu=0.19$ and $\mu=0.21$. The elasticity parameters are set as $E=1000$ MPa and $\nu=0.3$ for both the upper and lower blocks. We again use 10 load steps with a uniform displacement increment of $-0.01$ m.

Figure 11 comparatively shows the results of $\mu=0.19$ and $\mu=0.21$ cases in terms of the $x$-displacement field. We can find that the domain is under a slip condition when $\mu=0.19<\tan\theta$ and under a stick condition when $\mu=0.21>\tan\theta$. Therefore, we have confirmed that the phase-field method can also distinguish between stick and slip conditions appropriately.

Also for this problem, we check the sensitivity to the length parameter by repeating the same problem with $L=0.008$ m, 0.004 m, and 0.002 m. We have found that stick and slip conditions are correctly distinguished with all the three length parameters. Because the stick case is a standard linear elasticity problem, we only present the displacement and contact stress results of the slip case ($\mu=0.19$) in Figs. 12 and 13. It can again be seen that the numerical solutions show little sensitivity to $L$, at least for the length parameters considered which are reasonably small compared with the domain size.

Figure 14 shows the Newton convergence behaviors of the stick and slip cases when $L=0.002$ m and $L/h=8$. We can see that except the first load step of the stick case, all load steps converged after a single update, which evinces the linearity of the formulation. The first step of the stick case required multiple iterations because a slip condition is initially assumed for the stress-free initial condition and it has to be corrected by a Newton iteration. From the second step, as a stick condition is identified from the last converged step, the problem remained linear. To confirm this statement, we have also repeated the same problems by changing the initial contact condition to a stick condition. Then, as shown in Fig. 15, the first load step of the slip case required two iterations for convergence, and all other load steps in the stick and slip cases converged after a single Newton update. Therefore, we can conclude that the formulation is linear if the initial guess of the contact condition is correct, and that an incorrect guess of the contact condition can be rectified during a Newton iteration.

4.3 Sliding of a block

In our third example, we simulate a problem whereby stick and slip conditions coexist along an interface. The problem is an elastic block sliding on a rigid foundation, which was introduced by Oden and Pires [52] and later used by other works such as Wriggers et al. [53] and Simo and Laursen [25]. A slightly modified but essentially the same problem was also presented in Annavarapu et al. [31]. As depicted in Fig. 16, this problem considers a rectangular elastic block on a rigid foundation and applies tractions on its top and right boundaries. The elasticity parameters of the block are $E=1000$ kPa and $\nu=0.3$. The rigid foundation is approximated with a $10^{9}$ times larger Young’s modulus and a zero Poisson’s ratio, as done in Annavarapu et al. [31]. Emulating the setup of the original problem, the interface is frictional with $\mu=0.5$ in the 3.6 m-long middle part, but it is frictionless elsewhere. Under the given condition, the frictional part will mostly be sticked to the foundation but the frictionless parts will slip. We test two cases of length parameters, $L=0.016$ m and 0.008 m, with meshes of $L/h=8$. The problem is solved in a single step as in previous works.

Figure 17 compares deformed geometries obtained by our phase-field formulation with numerical result in Simo and Laursen [25] obtained by classical finite elements with an augmented Lagrangian method. It can be seen that the two results are fairly similar and that the interface is partially slipped. For a more direct comparison, in Fig. 18 the classical result is overlapped to phase-field results obtained with $L=0.016$ m and 0.008 m. We observe that the two results are matched remarkably well, for both the $L=0.016$ m and 0.008 m cases. This agreement again demonstrates that the phase-field formulation can correctly identify and reproduce stick and slip behaviors.

Previous studies have commonly used the contact normal pressures and tangential stresses of this problem to study the performance of contact algorithms. Here we also use these quantities to fully verify that the phase-field formulation can treat contact constraints well without an algorithm. Figure 19 presents the contact normal pressures and tangential stresses at the interface (when $L=0.008$ m), comparing them with data digitized from Simo and Laursen [25] and Oden and Pires [52]. One can see that the contact pressures and stresses of the phase-field and classical solutions are also in excellent agreement. This comparison has verified that the contact stresses and pressures of the phase-field method are accurate even when slip and stick conditions are mixed. Remarkably, embedded discontinuity methods have sometimes produced oscillatory solutions to the contact tractions of the same problem (see Fig. 8(d) of Annavarapu et al. [31] for example). However, as demonstrated throughout this section, the phase-field method is always free of such oscillation because it reformulates a contact problem as a continuum problem for which no algorithm is necessary for contact constraints. Therefore, an appealing feature of the phase-field method is that the method can address arbitrary crack geometry as embedded discontinuity methods, but without oscillations in contact pressures and stresses.

Lastly, we plot the Newton convergence behaviors of the $L=0.016$ m and 0.008 m cases in Fig. 20. As shown, this problem requires more iterations than prior examples because the interface here involves both stick and slip conditions (the initial contact condition was slip). Therefore, Newton’s method did not converge well initially. However, once the contact conditions of all points became correctly identified, the residual decreased at a rate for a linear problem. This behavior agrees well with our observation in the previous example.

4.4 Propagation of an inclined frictional crack

Following verification with stationary interface problems, we simulate propagation of a frictional crack to demonstrate the capability of the phase-field method for modeling growth of a crack with frictional contact. To this end, we now allow a crack to evolve according to fracture mechanics theory by solving the phase-field equation (26) in every load step. The phase-field equation and the momentum balance equation (25) are solved sequentially as proposed by Miehe et al. [54]. Because this sequential solution method is now fairly standard in the literature, its details are omitted for brevity. Also, because the phase-field equation has a physical meaning now, $W$ and $G_{c}$ in the phase-field equation should be calculated from physical quantities, rather than being arbitrarily assigned as before. Considering brittle shear fracture, we regard $W$ as the deviatoric part of strain energy and $G_{c}$ as the mode II fracture energy. In other words, we have modified a standard phase-field model for brittle fracture to accommodate frictional contact.

The setup of our particular problem is illustrated in Fig. 21. The domain is a 2 m wide and 4 m tall rectangle that possesses a 45*∘* inclined crack from coordinates (0.0,0.7) m to (1.3,2.0) m. The material parameters of the domain are: $E=10000$ MPa, $\nu=0.3$, and $G_{c}=50$ kJ/m2. To investigate the effect of friction on this problem, we consider three values of friction coefficients, $\mu=0.01$, 0.10, and 0.30. For phase-field modeling, we use $L=0.016$ m and locally refine the mesh until $L/h$ reaches 8 along the existing and expected crack path. Note that the mesh is structured and so the elements are not aligned with the crack direction. Once the preexisting crack is initialized as before, we vertically compress the domain with a constant displacement rate of $2\times 10^{-4}$ m per load step. Because the contact condition of this problem is rather simple, in most load steps Newton’s method converged after a single update.

Figure 22 shows how the phase-field variable and the vertical displacement field evolve during the course of loading when $\mu=0.3$. As shown, the phase-field model well simulates propagation of the preexisting crack along the 45*∘* direction until it reaches the upper right size of the domain. During the propagation stage, the displacement field is discontinuous across the crack but still continuous through the non-fractured region. After the crack has fully developed, however, the upper and lower parts of the domain are completely disconnected, and the upper part slips along the crack. Note that this post-fracture process is essentially the same as stationary interface problems simulated earlier in this section. We also note that other friction coefficient cases show qualitatively identical responses in terms of the crack path and the displacement pattern.

In Fig. 23 we plot the load–displacement curves of the three friction coefficient cases. As expected, the peak load and displacement increase with the friction efficient. We also see that the two cases show more or less the same pattern, in which the material fails in a brittle manner and then exhibits a residual strength. The residual strength also increases with the friction coefficient, which evinces the contribution from the frictional resistance along the crack.

Before closing this section, we would like to demonstrate the critical role of contact treatment in phase-field modeling of crack propagation under compression. For this purpose, we simulate the same problem with the model of Amor et al. [46], whereby the contact condition is treated by the volumetric–deviatoric decomposition of the stress tensor. This stress decomposition scheme is the only difference from our phase-field formulation used above. We note that when our formulation attempted to simulate this problem without friction ($\mu=0$), it did not converge from the very first load step because all nodes along the preexisting crack slip immediately. However, although the volumetric–deviatoric decomposition assumes frictionless contact, it can still simulate this problem until the crack fully develops. This indicates that the inexact contact treatment of the volumetric–deviatoric decomposition provides non-physical frictional resistance along the interface. On a related note, we have also found that the volumetric–deviatoric decomposition is unable to distinguish stick and slip conditions correctly for the second example of this section.

Figure 24 compares simulation results from the two phase-field formulations when $\hat{u}_{y}=-0.020$ m. To minimize the influence of friction on the comparison, the $\mu=0.01$ case is shown in this figure. One can see that when the volumetric–deviatoric decomposition is used for this problem, the crack path becomes kinked, giving rise to rather unrealistic deformation responses. This difference demonstrates that inappropriate estimation of contact stresses can impact the crack driving force to the extent that alters the crack path direction. Therefore, it can be concluded that accurate treatment of contact condition is critical to the application of phase-field modeling to compression-induced fracture propagation, which is a classic problem in geomechanics [55, 56, 57, 58].

5 Closure

A phase-field method has been proposed for modeling cracks with frictional contact. Built on standard approaches in phase-field modeling of fracture, the proposed method calculates the stress tensor in a regularized interface region by identifying the contact condition in an interface-oriented coordinate system. By doing so, the phase-field method accommodates contact behavior in the interface direction while imposing no-penetration constraints in other directions. Using benchmark examples in the literature, we have verified that the proposed method can provide numerical solutions very close to those obtained by a discrete method, showing little sensitivity to the length parameter for phase-field regularization. Moreover, by allowing the crack to evolve according to brittle fracture theory, we have demonstrated that the proposed phase-field method can also simulate propagation of frictional cracks.

The proposed phase-field method has two key features that make it an appealing alternative to standard discrete methods. First, it can model a crack passing through the interior of elements without an explicit representation of geometry or enrichment of basis functions. Second, it does not require a sophisticated algorithm for imposing contact constraints on crack surfaces. Thanks to these two features, the phase-field method can be implemented far more easily than most of existing methods for frictional cracks.

Therefore, it is believed that the phase-field method can be an attractive option even for modeling frictional interfaces that are stationary, let alone those that evolve.

Acknowledgments

This work was supported by the Research Grants Council of Hong Kong (Project 27205918). The first author also acknowledges financial support from the Hong Kong PhD Fellowship.