Mechanical Behavior of Nano Structures Using Atomic-Scale Finite Element Method (AFEM)

Abstract

This work presents a detailed description of the formulation and implementation of the Atomistic Finite Element Method AFEM, exemplified in the analysis of one- and two-dimensional atomic domains governed by the Lennard Jones interatomic potential. The methodology to synthesize element stiffness matrices and load vectors, the potential energy modification of the atomistic finite elements (AFE) to account for boundary edge effects, the inclusion of boundary conditions is carefully described. The conceptual relation between the cut-off radius of interatomic potentials and the number of nodes in the AFE is addressed and exemplified for the 1D case. For the 1D case elements with 3, 5 and 7 nodes were addressed. The AFEM has been used to describe the mechanical behavior of one-dimensional atomic arrays as well as two-dimensional lattices of atoms. The examples also included the analysis of pristine domains, as well as domains with missing atoms, defects, or vacancies. Results are compared with classical molecular dynamic simulations (MD) performed using a commercial package. The results have been very encouraging in terms of accuracy and in the computational effort necessary to execute both methodologies, AFEM and MD. The methodology can be expanded to model any domain described by an interatomic energy potential.

Keywords:
AFEM; Interatomic potentials; Lennard Jones potential; molecular dynamics

1 INTRODUCTION

The last decades have witnessed an enormous increase in the development of nanotechnologies, which have brought a large impact on areas such as engineering, medicine, biology, chemistry, as well as in aerospace, construction and manufacturing industries. The scientific community has been working to analyze and develop mechanical, electronic and optoelectronic devices at nano-scale, which can be applied to many areas which range from sensors engineering to medical therapies and drug delivery systems (The Royal Society and The Royal Academy of Engineering, 2004; Novoselov et al., 2012; Chen et al., 2013). Nanomechanics of materials is a new branch of mechanics which studies the properties and behavior of nanoscale materials and structures in response to applied loads and prescribed boundary conditions. From the dimensional perspective, a structure with at least one dimension smaller than 100 nm is considered to be a nanostructure (The Royal Society and The Royal Academy of Engineering, 2004).

The mechanical behaviour of nanoscale systems can be analyzed by using Ab-initio methods (Dirac, 1981), semi-empirical quantum methods (Tadmor and Miller, 2011) and also modified classical continuum mechanics. Ab-initio calculations are most accurate compared to semi-empirical quantum methods, but computationally very expensive and the analyses are still limited to a few hundred or thousand atoms. The molecular dynamics (MD) (Alder et al., 1959) has been one of the most commonly method used to analyze the behaviour of nanomaterials. Although MD methods can be applied to analyze larger systems, than those that can be treated by Ab-initio methods, it is still not applicable to systems consisting of many millions of atoms, which are frequently necessary to model practical devices.

For nanostructures greater than a few nanometres, multi-scale continuum methods such as Gurtin-Murdoch theory (GM-T) (Gurtin and Murdoch, 1975), that bridges surface and bulk energies, is quite accurate and computationally inexpensive (Dingerville et al., 2005; Miller et al., 2000; Liu et al., 2010; Sapsathiarn et al., 2012). However, GM-T and other continuum methods are not strictly applicable to model discrete atomic structures. The computational cost in modelling nanomaterials is a fundamental challenge. In order to overcome this challenge, the Atomic-scale Finite Element Method (AFEM) was developed and applied for multiscale analysis of carbon nanotubes (CNTs) by Liu and his co-authors (Liu et al., 2004; Liu et al., 2005). Presently, the AFEM can model and analyze the mechanical behavior of many materials at the nanoscale (Malakouti, 2006; Tao et al., 2016; Damasceno et al., 2015).

The AFEM is based on an energy approach. It requires an interatomic energy potential, describing local or non-local bonding forces of an atom interacting with a chosen set of other surrounding atoms. The choice of the interacting surrounding atoms and their array will, ultimately, lead to the formulation of the specific Atomic Finite Element. In principle, the method can be applied to all atomic systems that may be described by an interatomic energy potential.

The purpose of this work is to introduce and discuss fundamental issues about the formulation and application of the Atomic-Scale Finite Element Method (AFEM). This methodology will be applied, exemplary, to analyze the mechanical behavior of one-dimensional (1D) and two-dimensional (2D) atomic structures and lattices. The Lennard-Jones interatomic potential will be used in the formulation and analysis. In particular, the article describes how to relate the number of nodes of a 1D Atomic Finite Element with the cut-off radius of the Lennard-Jones potential. It also addresses, in detail, the modifications at the element level that have to be performed to model bounded atomic domains and to introduce proper boundary conditions. The obtained numerical results are compared with classical MD simulations. Accuracy and computational costs involved in both methodologies are addressed.

2 ATOMIC SCALE FINITE ELEMENT METHOD FORMULATION

The AFEM, as described by Liu et al. (2004), is based on an energy approach. It requires an interatomic energy potential, U, describing local or non-local bonding forces of an atom interacting with a chosen set of other surrounding atoms. The interatomic energy potentials are described in terms of the coordinates of the individual atoms, x_i, and some constitutive parameters of the considered atomic elements and arrays.

For a system with N atoms, the total interatomic energy considers the contribution of all interacting atoms:

${\text{U}}_{\text{tot}}\text{=}{\displaystyle \sum _{\text{i < j}}^{\text{N}}{\text{U}}_{ij}\left({\underset{\_}{\text{x}}}_{\text{j}}\text{-}{\underset{\_}{\text{x}}}_{\text{i}}\right)}$ (1)

Consider the energy W_f, related to the sum of the work of external forces, ${\underline{F}}_i^{ext}$ acting on individual atoms indicated by ‘i':

$W_f=\sum _{i=1}^N{\underline{F}}_i^{ext}·{\underline{x}}_i$ (2)

The total energy in this system with N atoms, E_tot, is composed of the interatomic bonding energy and the (eventual) work of external forces.

$E_{tot}=\sum _{i<j}^N{U}_{ij}\left({\underline{x}}_j-{\underline{x}}_i\right)-\sum _{i=1}^N{\underline{F}}_i^{ext}·{\underline{x}}_i$ (3)

As usual, the desired equilibrium condition is obtained by determining the stationary value of the total energy in the system with respect to the position x_i of the individual atoms:

$\frac{d{E}_{tot}}{d\underline{x}}=0$ (4)

The total energy E_tot can be expanded in a Taylor series around the equilibrium position of the atoms, x⁽⁰⁾:

${\text{E}}_{\text{tot}}\left(\underset{\_}{\text{x}}\right)\approx {\text{E}}_{\text{tot}}\left({\underset{\_}{\text{x}}}^{\left(\text{0}\right)}\right)\text{+}{\frac{{\text{dE}}_{\text{tot}}}{\text{d}\underset{\_}{\text{x}}}|}_{\underset{\_}{\text{x}}\text{=}{\underset{\_}{\text{x}}}^{\left(\text{0}\right)}}\left(\underset{\_}{\text{x}}\text{-}{\underset{\_}{\text{x}}}^{\left(\text{0}\right)}\right)\text{+}\frac{\text{1}}{\text{2}}{\left(\underset{\_}{\text{x}}\text{-}{\underset{\_}{\text{x}}}^{\left(\text{0}\right)}\right)}^{\text{T}}{\frac{{\text{d}}^{\text{2}}{\text{E}}_{\text{tot}}}{\text{d}\underset{\_}{\text{x}}\text{d}\underset{\_}{\text{x}}}|}_{\underset{\_}{\text{x}}\text{=}{\underset{\_}{\text{x}}}^{\left(\text{0}\right)}}\left(\underset{\_}{\text{x}}\text{-}{\underset{\_}{\text{x}}}^{\left(\text{0}\right)}\right)$ (5)

Defining the displacement u as the difference between the actual, x, and the equilibrium position, x⁽⁰⁾

$u=\left(\underline{x}-{\underline{x}}^{\left(0\right)}\right)$ (6)

And substituting the Equation (5) into Equation (4) leads to the AFEM equation system:

$\left[\text{K}\left(u\right)\right]\left\{u\right\}\text{=}\left\{\text{P}\right\}$ (7)

Where [K(u)] corresponds to the nonlinear stiffness matrix, {u}, the displacement increment vector and {P} the non-equilibrium load vector, respectively, given by:

$\left[K\left(u\right)\right]=\left[\frac{d^2{E}_{tot}}{d\underline{x}d\underline{x}}\overline{)\underline{x}}={\underline{x}}^{\left(0\right)}\right]=\left[\frac{d^2{U}_{tot}}{d\underline{x}d\underline{x}}\overline{)\underline{x}={\underline{x}}^{\left(0\right)}}\right]$ (8)

And

$\left\{\text{P}\right\}\text{=}\left\{{\text{-}\frac{{\text{dE}}_{\text{tot}}}{\text{d}\underset{\_}{\text{x}}}|}_{\underset{\_}{\text{x}}\text{=}{\underset{\_}{\text{x}}}^{\left(\text{0}\right)}}\right\}\text{=}\left\{\text{F}{\text{-}\frac{{\text{dU}}_{\text{tot}}}{\text{d}\underset{\_}{\text{x}}}|}_{\underset{\_}{\text{x}}\text{=}{\underset{\_}{\text{x}}}^{\left(\text{0}\right)}}\right\}$ (9)

3 LENNARD JONES POTENTIAL

In this section, a brief review of the Lennard-Jones (LJ) potential proposed by Jones (1924) is introduced. It is used to calculate the potential energy of the interaction between two non-bonding atoms based on their distance of separation, r_ij:

$U\left(r_{ij}\right)=\{\begin{array}{l}4{\epsilon }_{LJ}\left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6\right]for{r}_{ij}<r_c\hfill \\ 0\begin{array}{cccccc}& & & & & \end{array}for{r}_{ij}\ge r_c\hfill \end{array}$ (10)

Where ${\text{r}}_{\text{ij}}=\sqrt{{\left({\text{x}}_{\text{i}}{\text{- x}}_{\text{j}}\right)}^2+{\left({\text{y}}_{\text{i}}{\text{- y}}_{\text{j}}\right)}^2}$ is the bond length, and the r_c is some prescribed cut-off radius, which switches off the atomic interaction when the bond length exceeds its value. Only two Lennard-Jones parameters are sufficient to describe the intermolecular interactions, σ and ( _LJ , as shown in Figure 1. The parameter ( _LJ is the depth of the potential well and σ is the inter-particle distance at which the potential is zero. The ${\left(\raisebox{1ex}{$\sigma $}\!\left/ \!\raisebox{-1ex}{${\text{r}}_{\text{ij}}$}\right.\right)}^{12}$ term represents the strong repulsion and the high energy on the bond, as the distance between the pair of atoms decreases. The ${\left(\raisebox{1ex}{$\sigma $}\!\left/ \!\raisebox{-1ex}{${\text{r}}_{\text{ij}}$}\right.\right)}^6$ term represents the attractive term, which describes the weak attraction on the bond, as the distance between the pair of atoms increase.

The internal force, F(r_ij), between the two atoms is obtained by differentiating the potential with respect to the intermolecular distance, r_ij:

$\text{F}\left({\text{r}}_{\text{ij}}\right)\text{= -}\frac{\text{dU}\left({\text{r}}_{\text{ij}}\right)}{{\text{dr}}_{\text{ij}}}\text{= 4}{\epsilon }_{LJ}\left[\text{12}\frac{{\sigma }^{\text{12}}}{{\text{r}}_{\text{ij}}{}^{\text{13}}}\text{- 6}\frac{{\sigma }^{\text{6}}}{{\text{r}}_{\text{ij}}{}^{\text{7}}}\right]$ (11)

Figure 1
Lennard-Jones potential U(r), and internal force F(r).

Figure 1 shows the LJ potential, U(r), and the interatomic force, F(r), as a function of the distance, r_ij. By setting the force F(r) =0, the equilibrium distance between the atoms, r_eq is found to be:

${\text{r}}_{\text{eq}}{\text{= 2}}^{\frac{\text{1}}{\text{6}}}\sigma \text{= 1}\text{.1225}\sigma$ (12)

If the distance between two atoms is less than r_eq, the Lennard-Jones force is repulsive. If the distance is greater than r_eq, the force is attractive. The maximal force is obtained at the intermolecular distance,

${\text{r}}_{\text{fmax}}\text{=}\frac{\text{1}}{\sqrt[\text{6}]{\frac{\text{7}}{\text{26}}}}\sigma \text{= 1}\text{.2445}\sigma$ (13)

It should be noted that for intermolecular distances larger than r_fmax the interatomic bonding force decreases continuously. When the distance between atoms reaches the equilibrium distance between the atoms, r_eq, the bonding force is assumed to vanish.

The cut-off radius r_c is an important parameter to be chosen. It determines with how many neighboring atoms, a given atom is interacting. The cut-off radius r_c also determines whether the atomic interaction has a local or non-local character. In this article the cut-off radius of the interatomic bonding will be related to the number of nodes within a specific Atomic Finite Element. This will be illustrated in the next paragraphs within the context of 1D elements.

Consider an atomic chain with 5 atoms in equilibrium position, as shown in Figure 2a. Every considered atom interacts with one or more of its neighbors, according to the relation between the interatomic equilibrium distance r_eq and the chosen cut-off radius r_c. Three distinct conditions will be analyzed. First, consider that the cut-off radius obeys the relation (r_eq < r_c < 2r_eq). This means that one atom only interacts with the first nearest neighboring atoms, as shown in Figure 2b. For the second condition, (2 r_eq < r_c < 3 r_eq), every atom interacts with the first and the second nearest neighboring atoms, as shown in Figure 2c. The third condition, (3 r_eq < r_c < 4 r_eq), is illustrated in Figure 2d.

Figure 2
An atomic chain with 5 atoms and distinct relations between the equilibrium distance r_eq and the chosen cut-off radius of the LJ potentials, r_c.

Two aspects should be stressed. First, when the cut-off radius is smaller than two times the equilibrium distance, r_c <2 r_eq, the atomic bondings have a local character. For larger cut-off radius, r_c >2 r_eq, the bonding extends to the second or more distant neighbors and gives a non-local character to the interatomic interactions. The second comment is related to the finitude of the atomic chain considered. The atoms 1 and 5, indicated in Figure 2a, represent the physical limit of the considered chain. The lack of atoms to the left of atom number 1 and to the right of atom number 5 impacts directly the energy bonds of the chain. If r_c >2 r_eq, non-local bonds exist and not only the bond energy of the most external atoms (1 and 5, at the present case) if affected. This must be carefully considered when modelling a bounded atomic domain with the AFEM.

4 ATOMIC FINTE ELEMENTS

In this section the formulation for 1D and 2D Atomic Finite Elements are, exemplarily, presented for the case the Lennard-Jones interatomic potential. The idea is to focus on the fundamental issues of the AFEM formulation and not on material properties of the domains to be analyzed.

4.1 The 1D Atomic Element

The notion of the atomic finite element will be introduced based on three examples. Consider Figure 3, showing an atomic finite element with three nodes, designated as i-1, i and i+1. The definition of the atomic finite element and all the derived calculations, such as internal forces and stiffness matrix, are related to the central node, in this case, node i. It is assumed that the relation between the equilibrium distance and the potential cut-off radius satisfies (r_eq < r_c < 2 r_eq). The arrows in Figure 3 a shows that the central atom i, is only connected to its first neighbor atom to the left (i-1) and to the right (i+1). The total bonding energy for this 3 node atomic finite element i is composed of the bonding (i-1, i) and (i, i+1):

${\text{U}}_{\text{tot}}^{\text{i}}{\text{= U}}_{\text{i}\text{,i-1}}{\text{+ U}}_{\text{i}\text{,i+1}}$ (14)

Figure 3
A three-node atomic-scale finite element.

The expressions for element stiffness matrix, ${\text{K}}_{\text{i}}^{\text{3e}}$ and non-equilibrium force vector, ${\text{P}}_{\text{i}}^{\text{3e}}$ and for this atomic three node finite element are obtained by substituting the expression for the bonding energy (14) and the external forces in equations (8) and (9), resulting delivering, respectively:

${\text{K}}_{\text{i}}^{\text{3e}}\text{=}{\left[\begin{array}{ccc}\text{0}& \frac{\text{1}}{\text{2}}\frac{{\text{d}}^{\text{2}}{\text{U}}_{\text{tot}}}{{\text{d x}}_{\text{i-1}}{\text{d x}}_{\text{i}}}& \text{0}\\ \frac{\text{1}}{\text{2}}\frac{{\text{d}}^{\text{2}}{\text{U}}_{\text{tot}}}{{\text{d x}}_{\text{i-1}}{\text{d x}}_{\text{i}}}& \frac{{\text{d}}^{\text{2}}{\text{U}}_{\text{tot}}}{{\text{d x}}_{\text{i}}{\text{d x}}_{\text{i}}}& \frac{\text{1}}{\text{2}}\frac{{\text{d}}^{\text{2}}{\text{U}}_{\text{tot}}}{{\text{d x}}_{\text{i+1}}{\text{d x}}_{\text{i}}}\\ \text{0}& \frac{\text{1}}{\text{2}}\frac{{\text{d}}^{\text{2}}{\text{U}}_{\text{tot}}}{{\text{d x}}_{\text{i+1}}{\text{d x}}_{\text{i}}}& \text{0}\end{array}\right]}_{\text{3x3}}$ (15)

$P_i^{3e}={\left[{\underline{F}}_i^{ext}-\begin{array}{c}0\\ \frac{d{U}_{tot}}{d{x}_i}\\ 0\end{array}\right]}_{3x1}$ (16)

The Lennard-Jones potential is a pairwise potential, describing the interaction between two atoms, that is the interaction of atom (i) over atom (i+1) and the reciprocal interaction of atom (i+1) on atom (i). In the present formulation the atomic element is represented by node (i) and only the interaction of node (i) with (i-1) and (i) with (i+1) will be considered. This is implied by the directions of the arrows in Figure 3. The reciprocal interaction, of the neighboring atomic element centered in the node (i-1) on the node (i), of the present element, will be added later at the mesh assembly procedure. So only half of the bonding energy between (i) and (i-1), and (i) and (i+1) has to be considered. This is expressed by the 1/2 factor in the element stiffness matrix.

Figure 4a shows a 5-node atomic finite element. The central node is denoted by (i) and it interacts with the two nearest nodes, (i-1) and (i-2) to the left and (i+1) and (i+2) to the right. For this element the condition (2r_eq < r_c < 3r_eq) is considered. Analogously, for the 7-node element, shown in Figure 4b, the relation between equilibrium distance and cut-off radius is (3r_eq < r_c < 4r_eq).

Figure 4
(a) A five-node atomic-scale finite element, (b) a seven-node atomic-scale finite element.

The total bonding energy of these elements are, respectively, given by:

${\text{U}}_{\text{tot}}^{\text{i-5e}}{\text{= U}}_{\text{i}\text{,i-1}}{\text{+ U}}_{\text{i}\text{,i-2}}{\text{+ U}}_{\text{i}\text{,i+1}}{\text{+ U}}_{\text{i}\text{,i+2}}$ (17)

${\text{U}}_{\text{tot}}^{\text{i-7e}}{\text{= U}}_{\text{i}\text{,i-1}}{\text{+ U}}_{\text{i}\text{,i-2}}{\text{+ U}}_{\text{i}\text{,i-3}}{\text{+ U}}_{\text{i}\text{,i+1}}{\text{+ U}}_{\text{i}\text{,i+2}}{\text{+ U}}_{\text{i}\text{,i+3}}$ (18)

4.2 Atomic Element in Two Dimensions

In this section two-dimensional atomic finite elements are presented. The starting point is the 7-node 2D element shown in Figure 5. As it will be shown, it may be used to model the geometry of hexagonal Bravais type lattices (Liu et al., 2006). The central atom, (i), has connection with all of surrounding atoms, (i+1), (i+2), (i+3), (i+4), (i+5) and (i+6). The total energy for this basic 2D 7-node atomic finite element is the sum of all pairwise interacting nodes:

${\text{U}}_{\text{i}}^{\text{2D\_tot}}{\text{= U}}_{\text{i}\text{,i+1}}{\text{+ U}}_{\text{i}\text{,i+2}}{\text{+ U}}_{\text{i}\text{,i+3}}{\text{+ U}}_{\text{i}\text{,i+4}}{\text{+ U}}_{\text{i}\text{,i+5}}{\text{+ U}}_{\text{i}\text{,i+6}}$ (19)

Figure 5
Atomic finite element with seven nodes.

The element stiffness matrix and the element non-equilibrium force vector for the atomic finite element with seven nodes are given by, respectively:

$\left[{\underline{K}}_e^{2D}\right]=\left[\begin{array}{cccccccc}0& \cdots & 0& \frac{d^2{U}_{tot}}{d{x}_{\left(i+3\right)}d{x}_i}& \frac{d^{{}_2}{U}_{tot}}{d{y}_{\left(i+3\right)}d{x}_i}& 0& \cdots & 0\\ ⋮& \ddots & ⋮& \frac{d^2{U}_{tot}}{d{x}_{\left(i+3\right)}d{y}_i}& \frac{d^2{U}_{tot}}{d{y}_{\left(i+3\right)}d{y}_i}& ⋮& \ddots & ⋮\\ 0& \cdots & 0& ⋮& ⋮& 0& \cdots & 0\\ \frac{d^2{U}_{tot}}{d{x}_{\left(i+3\right)}d{x}_i}& \frac{d^2{U}_{tot}}{d{y}_{\left(i+3\right)}d{x}_i}& \cdots & \frac{d^2{U}_{tot}}{d{x}_{i}d{x}_i}& \frac{d^2{U}_{tot}}{d{y}_{i}d{x}_i}& \cdots & \frac{d^2{U}_{tot}}{d{x}_{\left(i+6\right)}d{x}_i}& \frac{d^2{U}_{tot}}{d{y}_{\left(i+6\right)}d{x}_i}\\ \frac{d^2{U}_{tot}}{d{x}_{\left(i+3\right)}d{y}_i}& \frac{d^2{U}_{tot}}{d{y}_{\left(i+3\right)}d{y}_i}& \cdots & \frac{d^2{U}_{tot}}{d{x}_{i}d{y}_i}& \frac{d^2{U}_{tot}}{d{y}_{i}d{y}_i}& \cdots & \frac{d^2{U}_{tot}}{d{x}_{\left(i+6\right)}d{y}_i}& \frac{d^2{U}_{tot}}{d{y}_{\left(i+6\right)}d{y}_i}\\ 0& \cdots & 0& ⋮& ⋮& 0& \cdots & 0\\ ⋮& \ddots & ⋮& \frac{d^2{U}_{tot}}{d{x}_{\left(i+6\right)}d{x}_i}& \frac{d^2{U}_{tot}}{d{y}_{\left(i+6\right)}d{x}_i}& ⋮& \ddots & ⋮\\ 0& \cdots & 0& \frac{d^2{U}_{tot}}{d{x}_{\left(i+6\right)}d{y}_i}& \frac{d^2{U}_{tot}}{d{y}_{\left(i+6\right)}d{y}_i}& 0& \cdots & 0\end{array}\right]$ (20)

${\underline{P}}_e=\left[\begin{array}{c}{\underline{F}}_i^{ext}-\\ {\underline{F}}_i^{ext}-\end{array}\begin{array}{c}\begin{array}{c}{\left(0\right)}_{6x1}\\ \frac{d{U}_{tot}}{d{x}_i}\\ \frac{d{U}_{tot}}{d{y}_i}\end{array}\\ {\left(0\right)}_{6x1}\end{array}\right]$ (21)

Figures 6a and 6b illustrate two other 2D atomic finite elements. These atomic finite elements are defined based on the types of Bravais lattices (Liu et al., 2006). shown in Figure 7a and 7b respectively. The extension of the AFEM to these elements is straight forward.

Figure 6
Atomic finite elements for different two-dimensional Bravais lattices.

Figure 7
Square and Trigonal f two-dimensional Bravais lattices.

5 AFEM IMPLEMENTATION

The assemblage procedure in one and two dimensions, the concept of modified atomic finite element and the application of boundary conditions will be discussed in this section. Figure 8 illustrates the AFEM assemblage procedure for a 1D atomic chain containing 4 atoms. The atomic finite element with three nodes is considered. The atomic chain is bounded so that to the left side of atom 1 and to the right side of atom 4 there are no atoms. Because of this lack of neighboring atoms, the atomic finite elements on the edges of the mesh must be modified. So even before the imposition of the boundary conditions, the potential energy of the edge element must be modified to account for the non-existence of one or more neighboring atoms.

Figure 8
Example of assemblage scheme for an atomic chain of 4 atoms.

Using a 3-node atomic finite element, the assemblage procedure for the 4 element chain, shown in Figure 8, requires 4 atomic finite elements, called, respectively, El=1, El=2, El=3, El=4. Initially, all the atomic elements have three local nodes, designated by (-1), (0) and (+1). In the sequence, the atomic elements, El=1 and El=4, must be modified. For the atomic El = 1, the contribution of energy from the central node (0) onto left node (-1) does not exist and, therefore, must be eliminated from the calculations of the stiffness matrix and global non-equilibrium force vector. Also, the degree of freedom of node (-1) of EI=1 must be eliminated from the assemblage procedure. This is graphically indicated by the dashed arrow between node (0) and node (-1) for EI=1. Analogously, for EI=4, the energy contribution of node (0) to the node (+1) must be eliminated from the energy calculation procedure. Again, the dashed arrow in EI=4, between node (0) and node (+1) illustrates the required modification.

The expression for the total energy of the 4 atom chain, Utotal, is given by the sum of the contributions of the energy of individual elements, U^EI, as shown in equation (22):

${\text{U}}_{\text{total}}\text{=}\frac{\text{1}}{\text{2}}\left({\text{U}}^{\text{El=1}}{\text{+U}}^{\text{El=2}}{\text{+U}}^{\text{El=3}}{\text{+U}}^{\text{El=4}}\right)$ (22)

The the expressions for the energy of each of the 4 atomic elements used to model the 4-atom chain of Figure 8 are given in the expressions (23) to (26):

${\text{U}}^{\text{El=1}}\text{= 4}{\epsilon }_{LJ}\left[{\left(\frac{\sigma }{{\text{r}}_{\text{12}}}\right)}^{\text{12}}\text{-}{\left(\frac{\sigma }{{\text{r}}_{\text{12}}}\right)}^{\text{6}}\right]$ (23)

${\text{U}}^{\text{El=2}}\text{= 4}{\epsilon }_{LJ}\left[{\left(\frac{\sigma }{{\text{r}}_{\text{21}}}\right)}^{\text{12}}\text{-}{\left(\frac{\sigma }{{\text{r}}_{\text{21}}}\right)}^{\text{6}}\right]\text{+4}\epsilon \left[{\left(\frac{\sigma }{{\text{r}}_{\text{23}}}\right)}^{\text{12}}\text{-}{\left(\frac{\sigma }{{\text{r}}_{\text{23}}}\right)}^{\text{6}}\right]$ (24)

${\text{U}}^{\text{El=3}}\text{= 4}{\epsilon }_{LJ}\left[{\left(\frac{\sigma }{{\text{r}}_{\text{32}}}\right)}^{\text{12}}\text{-}{\left(\frac{\sigma }{{\text{r}}_{\text{32}}}\right)}^{\text{6}}\right]\text{+4}\epsilon \left[{\left(\frac{\sigma }{{\text{r}}_{\text{34}}}\right)}^{\text{12}}\text{-}{\left(\frac{\sigma }{{\text{r}}_{\text{34}}}\right)}^{\text{6}}\right]$ (25)

${\text{U}}^{\text{El=4}}\text{= 4}{\epsilon }_{LJ}\left[{\left(\frac{\sigma }{{\text{r}}_{\text{43}}}\right)}^{\text{12}}\text{-}{\left(\frac{\sigma }{{\text{r}}_{\text{43}}}\right)}^{\text{6}}\right]$ (26)

It should be noticed that the energy expressions for the elements in equations (23) to (26) are given in terms of the distance between the atoms, considering the global numbering scheme. Elements EI=1 and EI=4 only present the energy contribution from one bond. On the other hand, elements EI=2 and EI=3 present two energy contributions, because these elements present neighboring atoms to the right and to the left.

The third important stage of the assemblage produce in AFEM is the application of the boundary conditions. Figure 9 illustrates the same atomic chain of Figure 8 but including a set of boundary conditions. In this case, the displacement of atom 1 is blocked, so, the degree of freedom of this atom must be excluded from the global stiffness matrix, and from of the global non-equilibrium force vector.

Figure 9
Introduction the boundary conditions at the 4 atoms chain.

The assemblage procedure of AFEM in two dimensions system is the same used for the one dimensional system. Figure 10 illustrates an example, which consists of a triangular (hexagonal) mesh with 38 atoms. The basic atomic finite element considered has 7 nodes, and was shown in Figure 5. In order to apply boundary conditions or to analyze meshes with defects and missing atoms, the basic atomic elements must also be modified to account for the missing neighboring atoms. This is illustrated in Figure 10. Element number 23 is an example of a complete atomic element. Elements 1, 6 and 25 are incomplete, or modified elements. Just like the 1D case, the bond energy of the missing neighboring atoms must be excluded from the total energy of the element.

Figure 10
Atomic structure with complete and modified elements.

The total energy expression for the mentioned elements 23, 1, 6, and 25 are:

${\text{U}}_{\text{tot}}^{23}{\text{= U}}_{\mathrm{23,16}}{\text{+ U}}_{\text{23}\text{,19}}{\text{+ U}}_{\text{23}\text{,20}}{\text{+ U}}_{\text{23}\text{,26}}{\text{+ U}}_{\text{23}\text{,27}}{\text{+ U}}_{\text{23}\text{,30}}$ (27)

${\text{U}}_{\text{tot}}^1{\text{= U}}_{\mathrm{1,4}}{\text{+ U}}_{\mathrm{1,5}}{\text{+ U}}_{\mathrm{1,8}}$ (28)

${\text{U}}_{\text{tot}}^6{\text{= U}}_{\mathrm{6,2}}{\text{+ U}}_{\mathrm{6,3}}{\text{+ U}}_{\mathrm{6,9}}{\text{+ U}}_{\mathrm{6,10}}{\text{+ U}}_{\mathrm{6,13}}$ (29)

${\text{U}}_{\text{tot}}^{25}{\text{= U}}_{\mathrm{25,18}}{\text{+ U}}_{\mathrm{25,22}}{\text{+ U}}_{\mathrm{25,29}}{\text{+ U}}_{\mathrm{25,32}}$ (30)

The total bond energy of the complete element 23 and of the modified elements 1, 6 and 25 is given in Table 1. The element energy determined within the atomistic finite element procedure and the corresponding energy obtained by the Molecular Dynamics (MD) package LAMMPS (Plimpton, 1995) are also indicated. It can be seen that the energy for each AFE is distinct, according to the number of missing surrounding atoms. Table 1 also shows that the energy values calculated by AFEM and MD are in accordance.

Table 1
Comparison of total element energy obtained from AFEM and MD.

	Total Energy
	$U_{tot}^{23}$ (eV)	$U_{tot}^1$ (eV)	$U_{tot}^6$ (eV)	$U_{tot}^{25}$ (eV)
AFEM	-3.0	-1.50	-2.50	-2.0
MD	-3.0	-1.50	-2.50	-2.0

The matrices and load vectors for the complete or modified elements are assembled and, in the sequence, the boundary conditions are imposed. It is stressed again, that modifying the atomic elements to account for missing bonds and the imposition of boundary conditions are two distinct solution steps within the AFEM. In the present implementation the resulting non-linear system, equation (7) was solved iteratively by the Newton-Raphson method.

6 RESULTS

6.1 Comparison of AFEM and MD in One Dimension

In this section, the 1D AFEM is investigated. Initially, the influence of the number of element nodes in the atomic elements is addressed. As already discussed, the number of nodes in the AFE is related to the chosen cut-off radius for the Lennard Jones potential. The first example is shown in Figure 11. It consists of a 7 atom chain, fixed at atom 1 and loaded by a force F at node 7. This example will be solve using atomic elements with 3, 5 and 7 nodes. Before imposing the boundary conditions, the distance between the atoms is set to be equal to the equilibrium position, r_eq = 1.122 Å. After applying the boundary conditions, the load vector is increased until the atomic chain fails by bond breaking. At every load step the non-linear equation (7) is solved iteratively.

Figure 11
A schematic diagram of a one-dimensional atomic chain.

Figure 12 shows the resulting force-strain relations. The strain measure ( is obtained by dividing the displacement of atom 7, u₇ by the original chain length L_o, ( = u₇/L_o. The curves show that there is very little difference between the results obtained considering the atomic elements with 5-nodes and 7-nodes. It suggests that 1D AFEM simulations using atomic finite element with 5 nodes is almost as accurate as atomic finite element with 7-nodes.

Figure 12
External force by strain (3, 5 and 7 nodes elements).

In the sequence, the accuracy of the obtained AFEM results will be validated by comparison with the results obtained from MD package LAMMPS. To simplify, the values of parameters of the LJ potential were assigned unit values, σ = 1 and ε_LJ =1. The AFEM results were determined for the 5 node element. In the MD simulations, the Lennard Jones potential was used at a temperature of 1 K. Non-periodic boundary conditions were used. For the AFEM and MD simulations the cut-off radius is set equal to r_c = 2.5σ (Smit, 1992). The time integration step for the MD simulations is 0.05 fs.

Figure 13a shows force-strain curves obtained by both methodologies. Figure 13c depicts the relative error between strains obtained by both solutions. There is a good agreement between both solutions, with a relative error smaller that 10^-5, except at the initial equilibrium position. The MD simulation needed 1000 iterations to achieve the equilibrium position for every load step. On the other hand, the number of iterations for the AFEM to reach a non-equilibrium force vector smaller than 10^-11 for distinct values of the load force is shown in Figure 13b. The number of iterative steps for every force in the AFEM varied from 4 to 8 according to the loading stage.

Figure 13
1D AFEM results - a) force-strain relation; b) number of iterations for AFEM, c) MD-AFEM relative error, d) convergence measure.

Figure 13d shows a typical behavior of the quantity u•P, that is the dot product between displacement vector and the global non-equilibrium force vector. It also represents a measure of convergence of the iterative solution (Kim, 2006). The figure shows that convergence obtained for the present results is very fast.

The results of Figures 13a to 13d suggest that the 1D AFEM formulation and implementation is correct, that the results are accurate and the number of iterative steps, required for the solution, is significant smaller than those of the typical MD simulations.

6.2 Comparison of AFEM and MD in Two Dimensions

This section is dedicated to the study of a 2D AFE. A two-dimensional hexagonal (or triangular) lattice, shown in Figure 14, presenting 655 atoms is analyzed. The 2D atomic finite element (AFE) with 7 nodes, given in Figure 5, is considered. Zero displacement boundary conditions are imposed on all atoms on the left edge, shown inside of the red dashed box. External force with the same amplitude are applied on all atoms on the right edge in x direction, as shown by the arrows in Figure 14. The Lennard Jones potential was used in MD at a temperature of 1 K. The Lennard-Jones potential parameters are, σ = 1= and ( _LJ = 1. The cut-off radius for the MD simulations is set to r_c = 1.5Å. Non-periodic boundary conditions were used.

Figure 14
A schematic diagram of a two-dimensional lattice.

Figure 15
(a) Force-strain relation; (b) force-displacement relation; (c) convergence of AFEM; (d) number of iterations for distinc force levels.

Figure 15a shows a force-strain diagram obtained by AFEM and MD. The results are obtained for the central atom at the left border of the mesh. A similar result, showing the change in total length of the atomic mesh, ∆L, is given if Figure 15b. The AFEM and MD solutions agree very well. The residual (u•P) is shown in Figure 15c and the number of iterations required to solve the non-linear system for every load step can be found in Figure 15d. The MD simulation required 50000 iterative steps to achieve a final configuration with 4 significant digits, whereas the AFEM never needed more than 7 steps to achieve the same accuracy.

As a second 2D example, consider the same mesh of the previous example, shown in Figure 16. Zero displacement boundary condition imposed to the atoms within the red dashed box at the left border. An external force F_ext is applied in the vertical direction upon atom number 739. The black lattice in Figure 16 represents the original atomic mesh. The blue mesh shows the deformed mesh at its last loading step before bond breaking. The atom 739 is connected to atoms 728 and 729 as illustrated in the detail shown in Figure 16. It is simple to see that the bonding force between atoms 739 and 729 is in compression and that the force between atoms 739 and 728 is traction.

Figure 16
Initial and final configuration of a 2D atomic mesh.

Figure 17 shows in dashed lines the interatomic bonding force developed by the Lennard-Jones potential as a function of the interatomic distance r _ij , as anticipated by equation (11). This figure also show the bonding forces between atoms 739-729 (compression) and between atoms 739-728 (traction) for distinct amplitude stages of the external force, F_ext. Traction is marked in blue color and compression in red color.

Figure 17
internal force by r.

As expected, the interatomic bonding forces start with a zero value at the equilibrium position and increase continuously in traction and compression for increasing values of the external load. The bonding forces follow exactly the path predicted in equation (11). It should be noticed that this is a highly non-linear behavior, in which interatomic displacements for compressed bonds is much smaller than those of the traction bonds for the same external loading level. Figure 17 also shows the break bonding situation. After the condition, in which the traction bond, in blue, reaches the maximum bonding force at r_fmax = 1.2445 σ, given by equation (13), the bonding forces decrease for larger distances, the bond can no longer withstand the imposed external loads and breaks. The iterative system will no longer converge, unless a new mesh without the lose atom is considered.

In the third example, the effect of missing atoms or vacancies on the mechanical behavior of an atomic mesh is investigated. Figure 18a shows a pristine mesh, with no vacancies. On the other hand, Figure 18b depicts the same mesh but with 4 missing atoms. In both cases, zero displacement boundary conditions are prescribed for the atoms in the right edge within the red dashed box. A force of equal, increasing amplitude is applied to the 4 atoms at the left border. The external load is increased incrementally until one or more atomic bonds within the mesh break. Figure 18c shows the stress (force)-strain behavior of the pristine and of the defect mesh. The same strain measure of the previous examples is used. It can be seen that vacancies do have a large effect on the force-strain behavior. For validation purposes, Figure 18d depicts the deformed configuration of the defect mesh (Figure 16b), determined by AFEM and MD. There is a good visual agreement between the two methodologies.

Figure 18
(a) Pristine mesh; (b) Defected mesh, (c) stress-strain curves, and (d) final equilibrium configuration from AFEM and MD.

As a last example, consider the square Bravais lattice shown in Figure 19a. The square 2D AFE with 5 nodes, shown in Figure 6a, is used to model the atomic domain. The displacement of all atoms with the box at the left edge of the mesh are blocked. Upon all atoms within the dashed box at the right side, an equal and increasing force is applied. Figure 19b shows the resulting force-strain curves obtained by both methodologies. Both curves show the same trend, but there is no perfect agreement. This should be further investigated, because all previous and even more complex meshes had presented a much better agreement between both methodologies.

Figure 19
(a) Square bravais lattice; (b) Force-strain curves obtained from AFEM and MD.

7 CONCLUSIONS

This work presents a detailed description of the formulation and implementation of the Atomistic Finite Element Method (AFEM), exemplified in the analysis of one- and two-dimensional atomic domains governed by the Lennard-Jones interatomic potential. The methodology to synthesize element stiffness matrices and load vectors, the potential energy modification of the atomistic finite elements (AFE) to account for boundary edge effects, the inclusion of boundary conditions were carefully described, in a way that the authors had not previously found in the AFEM literature. The conceptual relation between the cut-off radius of interatomic potentials and the number of nodes in the AFE is addressed and exemplified for the 1D case. For the 1D case elements with 3, 5 and 7 nodes were addressed. The AFEM has been used to describe the mechanical behavior or one-dimensional atomic arrays as well as two-dimensional lattices of atoms. The reported studies also included the analysis of pristine domains, as well as domains with missing atoms, defects, or vacancies. Almost all results were compared with classical molecular dynamic simulations (MD) performed using a commercial package. The results have been very encouraging in terms of accuracy and in the computational effort necessary to execute both methodologies, AFEM and MD. The methodology presented can be extended to other potential, like Tersoff (Tersoff, 1987), REBO (Brenner, 1990), AIREBO (Stuart, 2000), which have been applied to model nano structures composed of graphene and materials alike.

Acknowledgments

The research leading to this article has been funded by the São Paulo Research Foundation (Fapesp) through grants 2012/17948-4, 2013/23085-1, 2015/00209-2 and 2013/08293-7 (CEPID). Support from CAPES and CNPQ has been fundamental to the development of the present research. This is gratefully acknowledged.

References

• Alder, B.J., Wainwright, T.E.J. (1959). Studies in Molecular Dynamics. I. General Method, The Journal of Chemical Physics 31: 459-466.
• Brenner, D.W. (1990), Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films, Physical Review B 42: 9458.
• Chen, C., Hone, J. (2013). Graphene Nanoelectromechanical Systems. Proceedings of the IEEE 101(7).
• Damasceno, D.A., Rajapakse, R.K.N.D., Mesquita, E. (2015). Mechanical behavior of atoms using atomic-scale finite element method. Proceedings of the XXXVI Ibero-Latin American Congress on Computational Methods in Engineering.
• Dingerville, R., Qu, J., Cherakoui, M. (2005). Surface free energy and its effect on the elastic behavior of nano-sized particles wires and films, Journal of the Mechanics and Physics of Solids 53: 1827-1854.
• Dirac, P.A.M. (1981). The Principles of Quantum Mechanics, Clarendon Press, Oxford.
• Gurtin, M., Ian Murdoch, A. (1975). A continuum theory of elastic material surfaces, Archive for Rational Mechanics and Analysis 57: 291-323.
• Jones, J.E. (1924). On the determination of molecular fields-II. From the equation of state of a gas. Proceedings of the Royal Society A 106: 463-477.
• Kim, K. (2006). The Atomic-scale Finite Element Method for Analyzing Mechanical Behavior of Carbon Nanotube and Quartz, M.Sc. Dissertation, Virginia Polytechnic Institute, Virginia.
• Liu, B., Huang, Y., Jiang, H., Qu, S., Hwang, K.C. (2004). The atomic-scale finite element method, Computer Methods in Applied Mechanics and Engineering 193: 1849-1864.
• Liu, B., Jiang, H., Huang, Y., Qu, S., Yu, M.F. (2005). Atomic-scale finite element method in multiscale computation with applications to carbon nanotubes, Physical Review B, 72: 1098-0121.
• Liu, C., Rajapakse, R.K.N.D. (2010). Continuum models incorporating surface energy for static and dynamic response of nanoscale beams, IEEE Transactions on Nanotechnology 9: 422-431.
• Liu, W.K., Karpov, E.G., Park, H.S. (2006). Nano Mechanics and Materials: Theory, Multiscale Methods and Applications, John Wiley & Sons.
• Malakouti, M., Montazeri, A. (2016). Nanomechanics analysis of perfect and defected graphene sheets via a novel atomicscale finite element method, Superlattices and Microstructures 94: 1 -12.
• Miller, R.E., (2000). Size-dependent elastic properties of nanosized structural elements, Nanotechnology, 11:139-147.
• Nanoscience and nanotechnologies: opportunities and uncertainties (2004). The Royal Society and The Royal Academy of Engineering, London.
• Novoselov, K.S., Fal′Ko, V.I., Colombo, L., Gellert, P.R., Schwab, M.G., Kim, K. (2012). A roadmap for graphene, Nature 490: 192.
• Plimpton, S. (1995), Fast Parallel Algorithms for Short-Range Molecular Dynamics, Journal of Computational Physics 117: 1-19.
• Sapsathiarn, Y., Rajapakse, R.K.N.D. (2012). A model for large deflections of nanobeams and experimental comparison, IEEE Transactions on Nanotechnology 11, 247-254.
• Smit, B. (1992), Phase diagrams of Lennard-Jones fluids, The Journal of Chemical Physics 96: 8639-8640.
• Stuart, S.J.; Tutein, A.B.; Harrison, J.A. (2000), A reactive potential for hydrocarbons with intermolecular interactions, The Journal of Chemical Physics 112: 6472 -6486.
• Tadmor, E.B., Miller, R.E. (2011). Modeling materials: continuum, atomistic, and multiscale techniques. Cambridge University Press, New York.
• Tao, J., Sun, Y., (2016). The Elastic Property of Bulk Silicon Nanomaterials through an Atomic Simulation Method. Journal of Nanomaterials 2016: 6.
• Tersoff, J. (1987), New empirical approach for the structure and energy of covalent systems, Physical Review B 37: 6991.

Publication Dates

• Publication in this collection
  Dec 2017

History

• Received
  21 May 2017
• Accepted
  03 Aug 2017