A Posteriori Error Estimation and Adaptive Algorithm for the Atomistic/Continuum Coupling in 2D

# A Posteriori Error Estimation and Adaptive Algorithm for the Atomistic/Continuum Coupling in 2D

## Abstract

Atomistic/continuum coupling methods aim to achieve optimal balance between accuracy and efficiency. Adaptivity is the key for the efficient implementation of such methods. In this paper, we carry out a rigorous a posterior error analysis which includes the residual estimate, the stability constant estimate and the error bound, for a consistent atomistic/continuum coupling method in 2D. Adaptive mesh refinement algorithm can be designed correspondingly and the convergence rate with respect to degrees of freedom is optimal compare with a priori error estimates.

wangh@scu.edu.cn

mliao@xs.ustb.edu.cn

plin@maths.dundee.ac.uk

lzhang2012@sjtu.edu.cn

HW was partially supported by NSFC grant 11501389, 11471214 and Sichuan University Starting Up Research Funding No. 2082204194117. PL and ML were partially supported by NSFC grant 91430106 and Fundamental Research Funds for Central Universities Nos. 06108038 and FRF-BR-13-023. PL and HW were partially supported by EMS RSF grant. LZ was partially supported by NSFC grant 11471214, 11571314 and the One Thousand Plan of China for young scientists.

## 1Introduction

Atomistic/continuum (a/c) coupling methods are a class of computational multiscale methods that aim to combine the accuracy of the atomistic model and the efficiency of the continuum model for crystalline solids with defects [25]. To be precise, the atomistic model can be applied in a small neighborhood of the localized defects such as vacancies, dislocations, and cracks, while the continuum model (e.g., Cauchy-Born rule) can be employed away from the defect cores where elastic deformation occurs. The construction and analysis of different a/c coupling methods have attracted considerable attention in the research community in recent years [15]. We refer the readers to [22] for a review of such methods.

The goal of the mathematical analysis is to find the optimal relation of accuracy vs. degrees of freedom. The a priori analysis has been given for several typical a/c coupling methods, for example the QNL method (quasi-nonlocal quasicontinuum method) [23], the BQCE method (blended energy-based quasi-continuum method) [14], the BQCF method (blended force-based quasi-continuum method) [17], the GRAC method (geometric reconstruction based atomistic/continuum coupling) [36] and the BGFC method (atomistic/continuum blending with ghost force correction) [38]. In contrast, only few research articles are concerned with the a posteriori error control of the these methods. Prudhomme et al. [39] uses the goal-oriented approach to provide a posteriori error control for a three dimensional nanoindentation problem with the quantity of interest being the force acting on the indenter. The error estimator, whose effectiveness is validated only numerically, is given as a modification of a rigorous derived residual functional. Arndt and Luskin [1] use the same approach to analytically derive the a posteriori error estimates for a one dimensional Frenkel-Kontorova model. The goal oriented a posteriori error estimators are used to optimize the choice of the atomistic region as well as the finite element mesh in the continuum region. All of these work employ the original energy-based quasicontinuum method as the underlying model which is later shown to be inconsistent [6]. Recently, Kochmann et al. [45] proposed an adaptivity strategy for the so-call “fully-nonlocal quasi-continuum” method which apply a discrete model in the entire computational domain without coupling of different models. Such method is in a slightly different spirit from the problem we consider in the present work.

The a posteriori error bounds in the energy norms are first derived in [32] through a similar approach as the present work for two different types of coupling schemes which are (later characterized in [16]) the QCP method and the QNL method. However, the QCP method analyzed in [32] only contains finite element discretization without the coupling of different models while the extension to higher dimensions of the scheme in [26] is not clear. A recent advance in this direction [35] is the a posteriori error analysis of a consistent energy-based coupling method developed in [40], where the a posteriori error estimators in the energy norm and in energy itself were proposed.

Despite all those developments, the rigorous mathematical justification of a posteriori error estimate beyond 1D is still missing. In this paper, we present a rigorous a posteriori error estimate for a consistent energy-based a/c method in two dimension which is of physical significance, and has not been considered so far to the best knowledge of the authors. We use the residual-based approach [46] to establish the estimate in negative Sobolev norms following [35]. Two features distinguish our problem from the classic residual-based estimate for finite element approximation of the elliptic equations. The first one is the existence of the modeling error which is in origin different from the applications of quadrature rules. The second one is that the mesh may not be further refined when it almost coincides with the reference lattice but a model adaptation should be imposed. The 2D results rely on the analysis and implementation of divergence free field, which characterize the essential difference of 2D compared with 1D results in [26], where the analysis can be carried out by explicit calculations.

Similar to the a priori analysis of GRAC in [36], we constrain ourselves to the case of nearest-neighbor interactions. Although the analysis can be extended to finite range interactions and to other a/c coupling methods, we decide not to include these so that the main ideas and steps are clearly presented without the distraction from the unnecessary complexity of the presentation.

The paper is organized as follows. In § Section 2 we set up the atomistic, continuum and coupling models for point defects. In § Section 3.1 we derive the residual for the coupling scheme. In § Section 3.2 we obtain the stability bound for the coupling method, and then, in § Section 3.3 we present rigorous a posteriori error estimates. We present the numerical algorithms and results in § Section 4 and we give the conclusion and the outlook of our work in § Section 5. Some auxiliary results are given in § Section 6.

## 2Formulation

We will first give a brief review of a model for crystal defects in an infinite lattice in the spirit of [9] in § Section 2.1 and the Cauchy-Born continuum model in § Section 2.2. We then present a generic form of a/c coupling schemes in § Section 2.3, the consistent scheme GRAC will be introduced specifically in § Section 2.4.

### 2.1Atomistic model

#### Atomistic lattice and defects

Given , non-singular, is the homogeneous reference lattice which represents a perfect single lattice crystal formed by identical atoms and possessing no defects. is the reference lattice with some local defects. The mismatch between and represents possible defects , which are contained in some localized defect cores such that the atoms in do not interact with defects (see § ? and § ? regarding interaction neighbourhood). For example, for a crystal with a single point defect at , and one can choose a proper radius such that , where . For different types of point defects, we have

• for a vacancy at ;

• for an interstitial at but ;

• for an impurity at , the difference of the impurity atom with other atoms can be characterized by the inhomogeneity of interaction potentials (§ ?)

This characterization of localized defects can be straightforwardly generalized to multiple point defects and micro-cracks, for example, see the setup of the model problem in § Section 4.2. Straight screw dislocations can be enforced through the appropriate choice of boundary conditions [9].

#### Lattice function and lattice function space

Given , , denote the set of vector-valued lattice functions by

A deformed configuration is a lattice function . Let be the identity map, the displacement is defined by for any .

For each , prescribe an interaction neighbourhood with some cut-off radius . The interaction range is defined as the union of lattice vectors defined by the finite difference of lattice points in and .

To measure the error for lattice functions we need to introduce function norms and function spaces on the lattice. Define the “finite difference stencil” . Higher-order finite differences, e.g., and can be defined in a canonical way. A lattice function norm can hence be defined using those notations. For , let the lattice energy-norm (a discrete -semi-norm) be

The associated lattice function space is defined by

We choose

to be the collection of all the nearest neighbour bonds in the reference lattice, and for , denote . Then the energy norm can be reformulated as

The homogeneous lattice naturally induces a simplicial micro-triangulation . In 2D, , where . Let be the P1 nodal basis function associated with the origin; namely, is piecewise linear with resepct to , and and for and . The nodal interpolant of can be written as

We can introduce the discrete homogeneeous Sobolev spaces

with semi-norm . It is know from [30] that and are equivalent.

#### Interaction potential

For each , let denote the site energy associated with the lattice site . In this paper, we consider the general multibody interaction potential of the generic pair functional form [44]. Namely, the potential is a function of the distances between atoms within interaction range and there is no angular dependence. Accordingly, we have the following equivalent forms of interaction potentials,

A great number of practical potentials are in the form , including the widely used Embedded Atom Model (EAM) [5] and Finnis-Sinclair model [11]. For example, assuming a finite interaction neighborhood and an interaction range for , EAM potential reads

for a pair potential , an electron density function and an embedding function . The potential is homogeneous outside the defect region , namely, and for . and have the following symmetry: , and .

The energy of an infinite configuration is typically ill-defined. However, if we redefine the potential as the difference , which is equivalent to assuming , the energy functional

is a meaningful object. By [9], given some natural technical assumptions for the site potentials , is well-defined for , and if is in its variables, is times Fréchet differentiable. In particular, we define as the Lipschitz constant of .

Under the above conditions, the goal of the atomistic problem is to find a strongly stable equilibrium , such that

We call strongly stable if there exists such that

.

It is proven in [9] that, if the homogeneous lattice is stable, for a critical point of , , exhibit the following generic decay,

where .

### 2.2Continuum model

Let be the homogeneous site potential outside the defect core. To formulate atomistic to continuum coupling schemes we require a continuum model compatible with defined through a strain energy density function . A typical choice in the multi-scale context is the Cauchy–Born model [8], the energy density in the homogeneous crystal is defined by

### 2.3A/C coupling

We begin by giving a generic formulation of an a/c coupling method and employ concepts and notation from various earlier works, such as [25], but adapt the formulation to the settings in this paper.

First, we choose a computational domain : let be a simply connected, polygonal and closed set, such that for some . We decompose , where the atomistic region is again simply connected and polygonal, and contains the defect core: . Let be the radius of . Take as a regular simplicial partition (triangles for or tetrahedra for ) of the continuum region .

Next, we decompose the set of atoms into a core atomistic set and an interface set (typically a few “layers” of atoms surrounding ) such that . Let be the triangulation of induced by . is the associated triangulation of . Sometimes, it is also useful to define . Please see Figure 1 for an illustration of the computational mesh.

We may use instead of by dropping the subscript if there is no confusion. Let be the set of nodes in , and be the set of edges in .

Hence, the space of coarse-grained displacements is

Suppose is the Voronoi cell associated with a lattice point . For each , we can choose the associated effective Voronoi cell such that for , but for the effective Voronoi cell depends on the geometry of the interface (see [36]). Let be the volume ratio of with resepect to . For each element we define the effective volume of by

We note that only if .

Now we are ready to define the generic a/c coupling energy functional ,

where is a modified interface site potential which satisfies certain consistency conditions, and are suitable coefficients. The details will be discussed immediately in Section § Section 2.4 and references therein.

The goal of the approximate variational problem for a/c coupling is to find

in the subscript of and can be omitted if there is no confusion.

### 2.4Consistent Atomistic/Continuum Formulation

The choice of the interface potential in is the key in the formulation of atomistic/continuum coupling methods. In order to demonstrate the a posteriori error estimate for the generic a/c coupling methods, we shall restrict ourselves to the GRAC type methods [36].

#### The patch tests

A key condition that has been widely discussed in the a/c coupling literature is that should exhibit no “ghost forces”. We call this condition the force patch test, namely,

In addition, to guarantee that approximates the atomistic energy , it is reasonable to require that the interface potentials satisfy an energy patch test

If an a/c method satisfies the patch test and , it is called a consistent a/c method.

#### GRAC: Geometric reconstruction based consistent a/c method

To complete the definition of the consistent a/c coupling energy and of the associated variational problem , we must specify the interface region and the interface site potential. The geometric reconstruction approach was pioneered by Shimokawa et al [43], and then modified and extended in [7]. We refer to [37] for details of the implementation of geometric reconstruction based consistent atomistic/continuum (GRAC) energy for multibody potentials with general interaction range and arbitrary interfaces,. The extension of GRAC to 3D is a work in progress [10].

Our prototype implementation of GRAC is for the 2D triangular lattice with

Let , then , , are the nearest neighbour directions in .

Given the homogeneous site potential , we can represent in terms of . For each , let be free parameters, and define

A convenient short-hand notation is

we name the parameters as the reconstruction parameters.

The parameters are chosen so that the resulting energy functional satisfies the energy and force patch tests and . A sufficient (and likely necessary) condition for the energy patch test is that for all and . This is equivalent to

In addition, optimal condition and stabilisation mechanism were proposed in [37] and [28] to improve the accuracy and stability of GRAC scheme.

#### Stress formulation

The stress tensor based formulation can be obtained from the first variation of the energy. For any , and , there exist piecewise constant tensor fields , and which satisfy the following identities

here is the micro-triangulation induced by the reference lattice . We call an atomistic stress tensor, a continuum stress tensor, and an a/c stress tensor. For the nearest neighbour interactions, we can choose the following atomistic stress tensor, continuum stress tensor, and a/c stress tensor following respectively from the first variations -,

We call piecewise constant tensor field divergence free if

By definitions , it is easy to know that the force patch test condition is equivalent to that is divergence free for any constant deformation gradient . The discrete divergence free tensor fields over the triangulation can be characterized by the non-conforming Crouzeix-Raviart finite elements [36]. The Crouzeix-Raviart finite element space over is defined as

The following lemma in [36] characterizes the discrete divergence-free tensor field.

We also have the following corollary,

#### A Priori Error Estimate

In the analytical framework proposed in [19], the numerical error can be split into 3 parts: the modeling error due to the discrepancy between the atomistic model and the continuum model at the interface and the finite element edges, the coarsening error due to finite element discretization of the solution space in the continuum region, and the truncation error due to the finite size of the computational domain. It is proven in [9] that there exists a strongly stable solution to and a constant for GRAC such that,

With the generic decay property , and the following quasi-optimal conditions:

• the radius of the atomistic region satisfies,

• the finite element mesh is graded so that the mesh size function for satisfies,

It holds that there exists a constant , depending on , , , , and such that for sufficiently large,

In particular, when , and when P1 finite elements are used in the continuum region, we have,

where is the overall degrees of freedom.

## 3Error Analysis

We present the a posteriori error analysis in this section. In § Section 3.1, we derive the residual estimate for the consistent GRAC a/c coupling scheme introduced in § Section 2.4. Then, in § Section 3.2, we give a lower bound for the stability constant which is computable from the a/c solution . Finally, we put forward the main results Theorem ? and Theorem ? in § Section 3.3.

### 3.1Residual Estimate

To be more precise, we restrict ourselves to the case of nearest neighbour multibody interactions, namely, we use the so-called “grac23” method introduced in [36] as the a/c coupling mechanism. We will extend the formulation to general short-range multibody interactions in future works.

For lattice function , we denote its continuous and piecewise affine interpolant with respect to the micro-triangulation by . Identifying , we can define the (piecewise constant) gradient and the spaces of compact and finite energy displacements, respectively, by

It can be shown that that is dense in [9].

The first variation of the atomistic variational problem is to find such that

The first variation of the a/c coupling variational problem is to find such that

We introduce the truncation operator as in [9] by first choosing a cut-off function for and for . Define for by

where is defined by

The residual is defined as an operator on which is given by

By , denote , and take , where is the modified Clément operator [4] whose definition will be made clear in the following sections. By we can separate the residual into three groups

The first group , represents the truncation error; the second group represents the modeling error; and the third group represents the coarsening error. Note that the first variation of the coupling energy is given by

We will deal with the contribution from those three groups separately in the following subsections.

#### Truncation error

To analyze the truncation error , we need the Lemma 7.3 for the truncation operator in [9], namely, if the radius of the computational domain is sufficiently large (for example, in the nearest neighbour case, we only need ), the following estimates hold

where are independent of .

For any , the stress-based formulation of the first variation , the fact that for , the equivalence of and , and Cauchy-Schwarz inequality lead to,

Thus, the truncation error estimator is given by

#### Modeling error

In the analysis of the modeling error , the stress based formulation leads to,

As a result, we define the modeling error estimator by,

#### Coarsening error

For the coarsening error , we first observe that

Here, we choose to be the modified Clément interpolation [4] of . For any node in the triangulation , let be the piecewise linear function with respect to , and be the support of . The Clément interpolation operator can be defined as

By definition, satisfies the Dirichlet boundary condition. The Clement interpolation possesses the following property [3], for any element , interior edge ,

where is the diameter of , is the length of , the element patch , and the edge patch . The constants and depend only on the shape regularity of .

We also have the stability of the Clément interpolation,

where the constant only depends on the shape regularity of .

For notational convenience, we suppose that each interior edge has a prescribed orientation. and are the triangles on the left hand side and right hand side of the edge , and are the associated unit norm vector. The integration by parts of leads to,