On a \Gamma-convergence analysis of a quasicontinuum methodThis work is supported by a grant of the Deutsche Forschungsgemeinschaft (DFG) SCHL 1706/2-1.

On a -convergence analysis of a quasicontinuum methodthanks: This work is supported by a grant of the Deutsche Forschungsgemeinschaft (DFG) SCHL 1706/2-1.

Mathias Schäffner and Anja Schlömerkemper University of Würzburg, Department of Mathematics, Emil-Fischer-Straße 40, 97074 Würzburg, Germany. (mathias.schaeffner@mathematik.uni-wuerzburg.de, anja.schloemerkemper@mathematik.uni-wuerzburg.de)

In this article, we investigate a quasicontinuum method by means of analytical tools. More precisely, we compare a discrete-to-continuum analysis of an atomistic one-dimensional model problem with a corresponding quasicontinuum model. We consider next and next-to-nearest neighbour interactions of Lennard-Jones type and focus on the so-called quasinonlocal quasicontinuum approximation. Our analysis, which applies -convergence techniques, shows that, in an elastic setting, minimizers and the minimal energies of the fully atomistic problem and its related quasicontinuum approximation have the same limiting behaviour as the number of atoms tends to infinity. In case of fracture this is in general not true. It turns out that the choice of representative atoms in the quasicontinuum approximation has an impact on the fracture energy and on the location of fracture. We give sufficient conditions for the choice of representative atoms such that, also in case of fracture, the minimal energies of the fully atomistic energy and its quasicontinuum approximation coincide in the limit and such that the crack is located in the atomistic region of the quasicontinuum model as desired.

Key words. Quasicontinuum method; atomistic-to-continuum; –convergence; fracture.

AMS subject classifications. 49J45, 74R10, 74G15, 74G10, 74G65, 70C20

1 Introduction

The quasicontinuum (QC) method was introduced by Tadmor, Ortiz and Phillips [35] as a computational tool for atomistic simulations of crystalline solids at zero temperature. The key idea is to split the computational domain into regions where a very detailed (atomistic, nonlocal) description is needed and regions where a coarser (continuum, local) description is sufficient. The QC-method and improvements of it are successfully used to study crystal defects such as dislocations, nanoindentations or cracks and their impact on the overall behaviour of the material, see e.g. [25].

There are various types of QC-methods: Some are formulated in an energy based framework, some in a force based framework; further, different couplings between the atomistic and continuum parts and different models in the continuum region are considered. In the previous decade, many articles related to the numerical analysis of such coupling methods were published. We refer to [15, 23] for recent overviews, in particular on the large literature including work on error analysis.

In this article, we consider a one-dimensional problem and focus on the so-called quasinonlocal quasicontinuum (QNL) method, first proposed in [33]. The QNL-method and further generalizations of it (see e.g. [16, 30]) are energy-based QC-methods and are constructed to overcome asymmetries (so called ghost-forces) at the atomistic/continuum interface which arise in the classical energy based QC-method.
We are interested in an analytical approach in order to verify the QNL-method as an appropriate mechanical model by means of a discrete-to-continuum limit. This is embedded into the general aim of deriving continuum theories from atomistic models, see e.g. [3, Section 4.1], where also the need of a rigorous justification of QC-methods is addressed.
Our approach, announced in [34], is based on -convergence, which is a notion for the convergence of variational problems, see e.g. [6]. We start with a one-dimensional fully atomistic model problem which takes nearest and next-to-nearest neighbour interactions into account. The limiting behaviour of the corresponding discrete model was analyzed by means of -convergence techniques in [31] for a large number of atoms. In particular the -limit and the first order -limit are derived there, which take into account boundary layer effects.
From the fully atomistic model problem, we construct an approximation based on the QNL-method. In particular, we keep the nearest and next-to-nearest neighbour interactions in the atomistic (nonlocal) region and approximate the next-to-nearest neighbour interactions in the continuum (local) region by certain nearest neighbour interactions as outlined below. Furthermore, we reduce the degree of freedom of the energy by fixing certain representative atoms and let the deformation of all atoms depend only on the deformation of these representative atoms.
It turns out that the choice of the representative atoms has a considerable impact on the validity of the QC-method, see Theorem LABEL:5:cor, which is the main result of this work. This theorem asserts that the QC-method is valid if the representative atoms are chosen in such a way that there is at least one non-representative atom between two neighbouring representative atoms in the local region and in particular at the interface between the local and nonlocal regions. In Proposition LABEL:lemma:fracture, we prove that the mentioned sufficient condition on the choice of the representative atoms is indeed sharp by showing that in cases where the condition is not satisfied the limiting energy functional of the QC-method does not have the same minima as the limiting energy of the fully atomistic model and thus should not be considered an appropriate approximation. This implies by means of analytical tools that in numerical simulations of fracture one has to make sure to pick a sufficiently large mesh in the continuum region and at the interface.

The outline of this article is as follows. In Section 2 we present the two discrete models, namely the fully atomistic and the quasicontinuum model, in detail. In Sections 3 and 4 we investigate the limiting behaviour of the quasicontinuum energy functional by deriving the -limits of zeroth and first order. It turns out that the -limit of zeroth order of the fully atomistic and the quasicontinuum model coincide (Theorem LABEL:theorem:zero). If the boundary conditions are such that the specimen behaves elastically, we prove that both models also have the same -limit of first order (Theorem LABEL:theorem:elastic).
If the boundary conditions are such that fracture occurs, the quasicontinuum approximation leads to a -limit of first order (Theorem LABEL:theorem:fracture) that is in general different from the one obtained earlier for the fully atomistic model ([31], cf. Theorem LABEL:th:fracturennn). To compare the fully atomistic and the quasicontinuum model also in this regime, we analyze the -limits of first order further in Section 5. As mentioned above, it turns out that if we use a sufficiently coarse mesh in the continuum region, the minimal energies of the two first order -limits coincide (Theorem LABEL:5:cor). In fact we are able to show that in our particular model problem it is sufficient that the mesh size in the continuum region is at least twice the atomistic lattice distance. With this choice, fracture occurs always in the atomistic region as desired.
Furthermore, the -convergence results imply, under suitable assumptions, a rate of convergence of the minimal energy of the quasicontinuum model to the minimal energy of the fully atomistic model (Theorem LABEL:th:limmin). Finally, we show that the condition on the mesh size is sharp. In Proposition LABEL:lemma:fracture, we provide examples where the corresponding -limit has a different minimal energy and different minimizers than the fully atomistic system, which is due to poorly chosen meshes. This yields an analytical understanding of why meshes have to be chosen coarse enough in the continuum region.

Similar models as the one we consider here, were investigated previously in terms of numerical analysis. We refer especially to [14, 21, 26, 28, 29] where the QNL method is studied in one dimension. By proving notions of consistency and stability, those authors perform an error analysis in terms of the lattice spacing. To our knowledge, most of the results do not hold for “fractured” deformations. However, in [27] a Galerkin approximation of a discrete system is considered and error bounds are proven also for states with a single crack of which the position is prescribed. Recently, a different approach based on bifurcation theory is used in [22] to study the QC-approximation in the context of crack growth.

In [4], a different one-dimensional atomistic-continuum coupling method is investigated. Similar as in the QC-method the domain is splitted in a discrete and a continuum region. In the discrete part the energy is given by nearest neighbour Lennard-Jones interaction and in the continuum part by an integral functional with Lennard-Jones energy density. It is shown that fracture is more favourable in the continuum than in the discrete region. To overcome this, the energy density of the continuum model is modified by introducing an additional term which depends on the lattice distance in the discrete region. Furthermore, in [5, p. 420] it is remarked that if the continuum model is replaced by a typical discretized version, the fracture is favourable in the discrete region. As mentioned above, we here treat a similar issue in the QNL-method, see in particular Theorem LABEL:5:cor, Proposition LABEL:lemma:fracture.

The techniques of our analysis of the QNL method are related to earlier approaches based on -convergence for the passage from discrete to continuum models in one dimension, see [8, 9, 10, 11, 12, 31, 32]; see also [18, 19] for a treatment of two dimensional models. Recently, -convergence was used in [17] to study a QC approximation. In [17] a different atomistic model, namely a harmonic and defect-free crystal, is considered. Under general conditions it is shown that a quasicontinuum approximation based on summation rules has the same continuum limit as the fully atomistic system.

Common in all those works based on -convergence is that primarily information about the global minimum and minimizers are obtained. Since atomistic solutions are not necessary global minimizers, it would be of interest to obtain also results for local minimizers, for instance in the lines of [7, 9]. In this article, we treat systems with nearest and next-to-nearest neighbour interaction. A natural question is how the sufficient conditions on the choice of representative atoms change if we consider also interacting neighbours, . Therefore the corresponding fully atomistic model has first to be studied, which is part of ongoing research.

2 Setting of the Problem

First we describe our atomistic model problem which is the same as in [31]. We consider a one-dimensional lattice given by with and interpret this as a chain of atoms. We denote by the deformation of the atoms from the reference configuration and write as shorthand. We identify such functions with their piecewise affine interpolations and define

The energy of a deformation is given by

where and are potentials of Lennard-Jones type which will be specified in [LJ1]–[LJ4] below. Moreover, we impose boundary conditions on the first and last two atoms. For given we set


To consider only deformations which satisfy (LABEL:def:boundarycond), we define the functional


The goal is to solve the minimization problem

which we consider as our atomistic problem.

The idea of energy based quasicontinuum approximations is to replace the above minimization problem by a simpler one of which minimizers and minimal energies are good approximations of the ones for . Typically this new problem is obtained in two steps:

  • Define an energy where the long range (in our case next-to-nearest neighbour) interactions are replaced by certain nearest neighbour interactions in some regions.

  • Reduce the degree of freedom by choosing a smaller set of admissible functions.

To obtain (a), the next-to-nearest neighbour interactions are approximated as

see e.g. [28]. While this approximation turns out to be appropriate in the bulk, this is not the case close to surfaces, where the second neighbour interactions create boundary layers. This motivates to construct a quasinonlocal quasicontinuum model accordingly: For given let with . For we define the energy by using the above approximation for , cf. Fig. LABEL:fig:qcchain and keeping the atomistic descriptions elsewhere

Analogously to we define the functional

A crucial step for the following analysis is to rewrite the energy in a proper way. By defining


and , sometimes called Cauchy-Born energy density (see [28]), we can write


for satisfying (LABEL:def:boundarycond). To emphasize the local structure of the continuum approximation, we rewrite the summation over the terms with in (LABEL:hnlqc) as an integral. To this end we use the fact that is constant on for and thus



for satisfying (LABEL:def:boundarycond).

Figure 1: Illustration of the quasicontinuum approximation. Here denotes the scaled distance between the corresponding atoms in the deformed configuration and the two dotted lines stand for . Moreover, the black balls symbolise the repatoms.

To obtain (b) we consider instead of the deformation of all atoms just the deformation of a possibly much smaller set of so called representative atoms (repatoms). We denote the set of repatoms by with and define


Since we are interested in the energy for deformations , we define


In the following sections we study as tends to infinity. Therefore, we will assume that is such that


Hence, in particular . The above assumption corresponds to the case that the size of the atomistic region becomes unbounded on a microscopic scale (i), but shrinks to a point on a macroscopic scale (ii). While assumption (i) is crucial, see also Remark LABEL:rem:elastic, the assumption (ii) can be easily replaced by , and . In this case the analysis is essentially the same, but in the case of fracture, see Theorem LABEL:theorem:fracture, one has to distinguish more cases. We assume (LABEL:ass:kn) (ii) here because it is the canonical case from a conceptual point of view. Otherwise the atomistic region and continuum region would be on the same macroscopic scale.

3 Zero-Order -Limit

In this section we derive the -limit of the discrete energy (LABEL:def:hnlqct), which we refer to as zero-order -limit. This limit involves the convex and lower semicontinuous envelope of the effective potential energy which is already introduced in [11] defined by


We state the assumptions on the functions , and under which the following results are obtained.

  1. ] (strict convexity)

  2. ] (uniqueness of minimal energy configurations) For every such that we have where is defined as


    This implies

  3. ] (regularity and behaviour at , ). be in on their domains such that on its domain. Let and . Moreover, we assume the following limiting behaviour

  4. ] (structure of , and ). , are such that there exists a convex function


    and there exist constants such that


    Further, there exist such that


    and is strictly convex in on its domain for . Moreover, it holds and for all .

Remark 3.1

(a) The main examples we think of are Lennard-Jones interactions, defined classically as


and . The calculations in [31, Remark 4.1] show that defined as above satisfy [LJ1]–[LJ4]. Another example of interatomic potentials which satisfy the above assumptions, see [31, Remark 4.1], are Morse-potentials, defined for as


(b) The assumptions [LJ1]–[LJ4] imply that . In particular, we have


(c) Note that [LJ4] and (LABEL:lim:j012) imply that either or that there exists such that or for . In [LJ3], we assume for simplicity. However, this could be dropped making suitable assumptions on in the following statements.

To define appropriate function spaces, we use a similar notation as in [8] and [31]. Let be a function with bounded variation. Then we say that if satisfies the Dirichlet boundary conditions and . To allow jumps in respectively , the boundary conditions are replaced by respectively in this case. Analogously, we define for special functions with bounded variations and the above boundary conditions. Let (or in ), then we denote by the jump set of in , and for we set . Moreover we denote by the singular part of the measure with respect to the Lebesgue measure.

Let us now state and prove the zeroth-order -limit of the functional . It turns out that the limiting functional is equal to the -limit of the functional , cf. [31].

Theorem 3.2

Suppose [LJ1]–[LJ4] are satisfied and let . Let satisfy (LABEL:ass:kn) and let with be such that


Then the -limit of defined in (LABEL:def:hnl) and of defined in (LABEL:def:hnlqct) with respect to the –topology is the functional defined by

on .

Proof. The result for follows from [31, Theorem 3.1]. Thus we prove the result for . The following compactness property and lower bound follow from [10, Theorem 3.7] and [11, Theorem 3.1]. For the readers convenience, we present direct proofs here.

Compactness. Let be a sequence with equibounded energy . The definition of and the properties of , imply that . Define the set . Next, we make use of the fact that are bounded from below and that the energy is equibounded. Moreover, we apply (LABEL:ass:psi2) and Jensen’s inequality to obtain

for some independent of . By (LABEL:ass:psi1), we have that for some constant independent of . Moreover, by using the boundary conditions, we obtain

Since , we obtain by the Poincaré-inequality that is equibounded. Thus, we can extract a subsequence of which converges weakly to some , see [2, Theorem 3.23]. As argued in [31, Theorem 3.1], we have .

Liminf inequality. Let and be a sequence with equibounded energy which converges to in . The above compactness property and [2, Proposition 3.13] imply that converges to weakly in . By using [LJ3], [LJ4], we obtain for the recession function

with arbitrary. For every there exists such that for every . For large enough, we deduce from (LABEL:hnlqc:integral) by the definition of and [LJ4]

Note that by it follows for all , thus there exists such that

The last inequality is a direct implication of [2, Theorem 2.34], using that weakly converges to . By using that the right-hand side above is finite only if , we obtain the liminf inequality from the arbitrariness of .

Limsup inequality. To show the existence of a recovery sequence, we first do not take the boundary conditions into account. Therefore, we define the functional by

For every we show existence of a sequence converging to in such that


As outlined in the proof of [10, Theorem 3.5] it is enough to show the above inequality for linear and for with a single jump: by density, this proves the statement for and the general estimate follows by relaxation arguments. Firstly, we consider functions with a single jump. Let with , and . By (LABEL:T_n) there exists with and such that for . We define now a sequence by


Obviously we have in . The functions are defined such that for and for all with . Using , (LABEL:T_n), [LJ3] and [LJ4] this implies

Now let for some . For every sequence satisfying (LABEL:T_n) we find a sequence of natural numbers such that

We define for every a set with , where such that there exist which satisfy

From we deduce and thus . Let us define such that and

By using and , we obtain

and thus in . Indeed, by , and , the last term tends to zero as . For the limsup inequality we argue similarly as in the case of a jump before. By definition, we have