Analysis of Energy-Based Blended Quasicontinuum Approximations

We begin by presenting a model of a one-dimensional lattice of atoms. Our reference configuration is the integer lattice $\epsilon\mathbb{Z}$ with lattice spacing $\epsilon>0$. The admissible deformations are N-periodic, mean-zero displacements of uniformly strained states of $\epsilon\mathbb{Z}$.

Precisely, for $N\in\mathbb{N}$, we let

be the set of N-periodic, mean-zero displacements. For $F\in(0,\infty)$ we let

be the uniformly deformed state with macroscopic strain $F$, and we let

be the set of admissible deformations with fixed macroscopic strain $F$. We will fix $\epsilon=\frac{1}{N}$ throughout the remainder of the paper, so that the reference length of a period is one, independently of $\epsilon$. For notational convenience, we have let the deformations $y\in\mathcal{Y}_{F}$ be functions of $\mathbb{Z}$ instead of $\epsilon\mathbb{Z}$ even though we still think of $\epsilon\mathbb{Z}$ as the reference configuration. We also remark that the “displacements” $u$ which appear in the definition of $\mathcal{Y}_{F}$ are not displacements measured relative to the reference lattice $\epsilon\mathbb{Z}$, but are instead displacements relative to the uniformly deformed state $y^{F}$.

In the language of mechanics, one would say that our choice of the deformation space $\mathcal{Y}_{F}$ means that we have imposed “periodic boundary conditions with macroscopic strain $F$.” Periodic boundary conditions are mathematically similar to Dirichlet boundary conditions in that constraints are applied to the space of admissible deformations. In our case, the constraint is given by the macroscopic strain $F$ which models the macroscopic scale stretching or contraction of the chain. Periodic boundary conditions are useful because they eliminate surface effects near the boundaries of the chain.

For a deformation $y\in\mathcal{Y}_{F},$ we define a stored energy per period which we call $\Phi(y)$. Our stored energy $\Phi(y)$ sums the energies of all the interactions between first nearest neighbors and second nearest neighbors in the lattice. We compute the energies of interaction using a potential $\phi:(0,\infty]\rightarrow\mathbb{R}$ which we assume satisfies

1. (1)

   $\phi\in C^{4}((0,\infty],\mathbb{R});$

2. (2)

   there exists $r^{*}$ so that $\phi^{\prime\prime}(r)>0$ for all $r\in(0,r^{*})$, and $\phi^{\prime\prime}(r)<0$ for all $r\in(r^{*},\infty);$

3. (3)

   $\phi$ and its derivatives decay at infinity.

Examples of commonly used interaction potentials which satisfy the above conditions include the Lennard-Jones potential and the Morse potential. Specifically, we define

In (2.1), we call the terms of the form $\phi(y^{\prime}_{\xi})$ first nearest neighbor interactions, and we call terms of the form $\phi(y^{\prime}_{\xi}+y^{\prime}_{\xi+1})$ second nearest neighbor interactions.

The reader will observe that we have scaled the energy by $\epsilon$ and the interatomic distances by $\epsilon^{-1}$. This is done for two related reasons. First, quasicontinuum methods are only needed when when $N$ is large (equivalently, when $\epsilon$ is small). By introducing the scaling above, we are able to distinguish those parts of the error which are most significant in the practically relevant limit of large $N$. The reader should bear this in mind when interpreting our error estimates in Sections 3, 4, and 5. Second, we desire that in the limit as $\epsilon$ tends to zero (equivalently, as $N$ tends to infinity), the atomistic energy converges to the fully continuum Cauchy-Born energy (2.5). This is explained in detail in Section 2.2 below. Roughly speaking, this ensures that the atomistic energy is compatible with the continuum energy to which it is coupled in the BQCE and BQNL schemes. The compatibility of the two energies is reflected in the fact that the modeling error of the Cauchy-Born energy (2.4) is of order $\epsilon^{2}$ (See Remark 3.2).

Now suppose we want to model the response of our lattice to a dead load. We let $\mathcal{U}^{*}$ denote the algebraic dual of $\mathcal{U}$, and we equip the space $\mathcal{U}$ with the $\ell^{2}_{\epsilon}$ inner product given by

Generally speaking, a dead load $f$ should be thought of as an element of $\mathcal{U}^{*}$. However, since $\mathcal{U}$ is a Hilbert space with the inner product (2.2), any $g\in\mathcal{U}^{*}$ is represented by some $f\in\mathcal{U}$. That is, for any $g\in\mathcal{U}^{*}$ there exists $f\in\mathcal{U}$ so that for all $v\in\mathcal{U}$, $g(v)=<f,v>$. Thus, we will think of dead loads as elements of $\mathcal{U}$.

The total energy for an atomic chain subject to a dead load $f\in\mathcal{U}$ is then given by

and the equilibrium deformation of the atoms solves the minimization problem:

Let $\delta\Phi$ denote the first variation of the energy $\Phi$, and let $\delta^{2}\Phi$ denote the second variation. We will say that a solution $y\in\mathcal{Y}_{F}$ of the minimization problem (2.3) is strongly stable if $\delta^{2}\Phi(y)[u,u]>0$ for all $u\in\mathcal{U}\setminus\{0\}$, and we will call a deformation $y\in\mathcal{Y}_{F}$ an equilibrium if it solves

We call an energy $\Psi(y)$ local if it can be written in the form

where $\psi:(0,\infty]\rightarrow\mathbb{R}$ is called a strain energy density. Otherwise, we say that an energy is nonlocal. Observe that energy (2.1) is nonlocal. The development and use of a strain energy density $\phi$ allows coarse-graining in the continuum region by utilizing the full range of algorithms, codes, and analysis developed for the continuum finite element method. (An alternative quasicontinuum approach uses quadrature to approximate the nonlocal coarse-grained energy without the construction of a strain energy density [21]. Peridynamics offer another promising fully nonlocal approach to coarse-graining [37].) Thus, it is useful to devise local energies which approximate (2.1). The key to developing such energies is to observe that as long as $y^{\prime}_{\xi}$ varies slowly between neighboring lattice points, we have

If we replace the second nearest neighbor terms $\phi(y^{\prime}_{\xi}+y^{\prime}_{\xi+1})$ by $\frac{1}{2}\{\phi(2y^{\prime}_{\xi})+\phi(2y^{\prime}_{\xi+1})\}$ in (2.1), we then obtain the Cauchy-Born energy

which is a local energy commonly used to approximate (2.1).

An alternative derivation of the Cauchy-Born energy is as follows. Suppose that $Y:[0,1]\rightarrow\mathbb{R}$ is a $C^{\infty}$ deformation of $[0,1]$, and define $y^{\epsilon}_{\xi}=Y(\epsilon\xi)$. It can be shown that

where $\psi^{cb}(r):=\phi(r)+\phi(2r)$ is called the Cauchy-Born strain energy density. We call the energy $\Psi^{cb}$ the fully continuum Cauchy-Born energy. We refer the reader to [8] for a detailed discussion of the convergence of $\Phi(y^{\epsilon})$ to $\Psi^{cb}(Y)$ as $\epsilon$ goes to zero.

Now suppose that $Y$ is in the Lagrange $P^{1}$ finite element space on $[0,1]$ with nodes at the reference position of each atom; so the nodes are $\epsilon\xi$ for $\xi=1,\dots,N$. In that case, it is easy to show that

Thus, we see that the Cauchy-Born energy $\Phi^{cb}$ is a finite element approximation of the fully continuum Cauchy-Born energy $\Psi^{cb}$. Consequently, we will refer to $\Phi^{cb}$ as a “continuum” energy. We say that the energy $\Phi^{cb}$ is not coarse-grained since $\Phi^{cb}$ has the same number of degrees of freedom as the atomistic energy $\Phi$: both depend on the deformed position of each atom.

Suppose we want to approximate a solution to the minimization problem (2.3) which has defects in localized regions but which is smooth throughout most of the lattice. Efficiency requires that most of the lattice be modeled using the Cauchy-Born energy (2.4), while accuracy requires that the defects be treated using the nonlocal energy (2.1). Thus, we desire a coupling of the two models.

We distinguish two approaches to deriving coupled energies. We call these the strain-energy based approach and bond based approach. In the strain-energy based approach, one first defines an atomistic energy per atom. For our one-dimensional chain, an appropriate energy per atom is

We call $\Phi^{a}_{\xi}(y)$ the atomistic energy at atom $\xi$. Observe that $\Phi^{a}_{\xi}(y)$ is half the total energy of all bonds involving atom $\xi$, so $\Phi(y)=\epsilon\sum_{\xi=1}^{N}\Phi^{a}_{\xi}(y)$. The energy per atom should be interpreted as the energy per length $\epsilon$ at atom $\xi$; thus, $\Phi_{\xi}(y)$ is analogous to a continuum strain-energy density.

To define a coupled energy we now blend the atomistic energy per atom with the continuum energy. Schematically, we choose some blending functions $\alpha:[0,1]\rightarrow[0,1]$ and $\beta:\epsilon\mathbb{Z}\rightarrow[0,1]$. Then for a deformation $Y:[0,1]\rightarrow\mathbb{R}$ with corresponding atomistic deformation $y^{\epsilon}_{\xi}:=Y(\epsilon\xi)$, we define

In equation (2.7), one should imagine that the function $\alpha$ is zero in a small region containing any defects, and one throughout the bulk of the lattice. One should imagine that $\beta$ is one near defects, and zero throughout the bulk of the lattice. For an abruptly coupled energy, one should imagine that $\alpha$ is the characteristic function of the continuum region, and $\beta$ is the characteristic function of the atomistic region. By contrast, for a blended energy, one should imagine that $\alpha$ and $\beta$ are more smooth. The energy-based quasicontinuum (QCE) method [28], the bridging domain method [6], and the Arlequin method [5] all take essentially the form (2.7). Moreover, although the AtC method [3] was derived in variational form, the equilibrium condition is the Euler-Lagrange equation of an energy similar to (2.7).

We will now explain the QCE energy in detail. Let $[0,1]=\mathcal{A}\cup\mathcal{C}$ be a partition of the domain $[0,1]$ into an atomistic region $\mathcal{A}$ and a continuum region $\mathcal{C}$. Suppose that $Y$ is in the $P^{1}$ Lagrange finite element space with nodes at every lattice site, and let $y_{\xi}:=Y(\epsilon\xi)$ be the corresponding deformation of the atomistic lattice $\epsilon\mathbb{Z}$. In this special case, the QCE energy reduces to

We call $\Phi^{c}_{\xi}$ the continuum energy at $\xi$.

Our blended energy-based quasicontinuum energy (BQCE) is based on a blend of the atomistic energy at $\xi$ (2.6) and the continuum energy at $\xi$ (2.8). Let $\gamma:\mathbb{Z}\rightarrow[0,1]$ be a blending function. For $y\in\mathcal{Y}_{F}$ we define

We observe that the QCE energy is simply the BQCE energy whose blending function $\gamma$ is the characteristic function of the continuum region. We will show in Theorem 3.1 that BQCE does not pass the patch test where we recall that a quasicontinuum energy passes the patch test if there are no forces under uniform strain; for our one-dimensional problem an energy $\Phi^{coupled}$ passes that patch test if $\delta\Phi^{coupled}(y^{F})=0$ for all $F\in(0,\infty)$. We will call a method which passes the patch test patch test consistent.

We will now discuss the QNL and BQNL energies in detail. The atomistic and continuum energies of the bond between the nearest neighbors at reference positions $\epsilon\xi$ and $\epsilon(\xi-1)$ are taken to be $\phi(y^{\prime}_{\xi})$. This is the same as the energy of that bond in the fully atomistic model (2.1). The atomistic energy of the bond between the second nearest neighbors at reference positions $\epsilon(\xi+1)$ and $\epsilon(\xi-1)$ is taken to be $\phi(y^{\prime}_{\xi+1}+y^{\prime}_{\xi})$, and the continuum energy of that bond is taken to be $\frac{1}{2}\{\phi(y^{\prime}_{\xi+1})+\phi(y^{\prime}_{\xi})\}$. As discussed in Section 2.2, this choice of continuum bond energy is related to the Cauchy-Born rule. The coupled energy is then given by

where $[0,1]=\mathcal{A}\cup\mathcal{C}$ is a partition of $[0,1]$.

Our blended quasinonlocal QC (BQNL) energy is based on a smooth blending of the atomistic and continuum bond energies used to define the QNL energy (2.10). Let $\eta:\mathbb{Z}\rightarrow[0,1]$ be a blending function. We define

We remark that the QNL energy is the BQNL energy whose blending function is the characteristic function of the atomistic region $\mathcal{A}$.

We will now give a simple proof that BQNL passes the patch test. We discuss the modeling error of BQNL in more detail in Section 3. The patch test consistency of BQNL is a consequence of the following result.

In light Lemma 2.1, one expects that BQNL will inherit many of the properties of the QNL energy. In particular, since any BQNL energy is an affine combination of QNL energies, and since the QNL energy is patch test consistent [39], the BQNL energy passes the patch test. Moreover, we show in Remark 4.2 that the BQNL energy predicts the critical strain of the atomistic energy as accurately as the QNL method.

For our analysis, it will be convenient to define a single energy which incorporates both BQCE and BQNL as special cases. We will use the blended quasicontinuum (BQC) energy

where $\alpha,\beta:\mathbb{Z}\rightarrow[0,1]$ are blending functions. We observe that the BQCE energy with blending function $\gamma$ is the same as the BQC energy with blending functions

The BQNL energy with blending function $\eta$ is the BQC energy with blending functions

Here we collect all the notation and auxiliary theorems which we will use below. The reader may feel free to skim this section at first and return only as necessary.

For differences, we write

For means, we write

For $y\in\mathcal{Y}_{F}$, we will think of $y^{\prime}$, $y^{\prime\prime}$, $y^{\prime\prime\prime}$, and $\overline{y}$ as N-periodic functions defined on the reference configuration.

We will also use certain bounds on the potential function $\phi$ and its derivatives:

We will denote the convex hull of the set $A$ by

We define the following norms

Correspondingly, we let ${\ell^{p}_{\epsilon}}$ denote the space $\mathcal{U}$ equipped with the norm $\|\cdot\|_{{\ell^{p}_{\epsilon}}}$, and we let $\mathcal{U}^{{1},{p}}$ denote the space $\mathcal{U}$ equipped with the norm $\|\cdot\|_{\mathcal{U}^{{1},{p}}}$. Additionally, let $Y^{*}$ denote the topological dual of the Banach space $Y$, and let $\mathcal{U}^{{-1},{p}}:=(\mathcal{U}^{{1},{q}})^{*}$ where $p,q\in[1,\infty]$ with $\frac{1}{p}+\frac{1}{q}=1$. We let ${\mathcal{A}}$ denote the atomistic region, ${\mathcal{C}}$ denote the continuum region, and ${\mathcal{I}}$ denote the interface. We let $\|\cdot\|_{{\ell^{p}_{\epsilon}}(\mathcal{P})}$ be the $\|\cdot\|_{{\ell^{p}_{\epsilon}}}$ norm taken over the set $\mathcal{P}$ for $\mathcal{P}={\mathcal{A}},{\mathcal{C}},\mbox{ or }{\mathcal{I}}$. We denote the closed ball of radius $r$ at $x$ in $X$ by

for $X$ one of the spaces ${\ell^{p}_{\epsilon}}$ or $\mathcal{U}^{{1},{p}}$.

We will write $\delta\Psi$ for the first variation of a differentiable energy functional $\Psi$. The first variation is a map from $\mathcal{Y}_{F}$ into $\mathcal{U}^{*}$; we will let $\delta\Psi(y)[u]$ denote the first variation of $\Psi$ at the deformation $y\in\mathcal{Y}_{F}$ evaluated on the test function $u\in\mathcal{U}$. We will use the letter $u$ to denote a test function belonging to $\mathcal{U}$ throughout the remainder of the paper. The letter $y$ will be used to denote a deformation belonging to $\mathcal{Y}_{F}$. We warn the reader that in expressions such as $\delta\Psi(y)[u]$, $u$ denotes an arbitrary test function, not the displacement corresponding to the deformation $y$. Similarly, we will let $\delta^{2}\Psi$ denote the second variation of $\Psi$. The second variation can be interpreted either as a map from $\mathcal{Y}_{F}$ into the space of bilinear forms on $\mathcal{U}$, or as a map from $\mathcal{Y}_{F}$ into $L(\mathcal{U},\mathcal{U}^{*})$. We will let $\delta^{2}\Psi(y)[u,v]$ denote the second variation of $\Psi$ at $y\in\mathcal{Y}_{F}$ evaluated on the test functions $u,v\in\mathcal{U}$.

We will need the following version of the Inverse Function Theorem which appears as Lemma 1 in [31].

In Section 5, we will use Theorem 2.1 to prove a priori existence with error estimates for the BQCE and BQNL energies. We will apply Theorem 2.1 using bounds on $\eta$ derived from the modeling estimates in Section 3 and bounds on $\sigma$ derived from the stability estimates given in Section 4.