A new multiscale methodology for modeling of single and multi-body solid structures

Computational Materials Science 63 (2012) 1–11

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A new multiscale methodology for modeling of single and multi-body solid structures M.H. Korayem ⇑, S. Sadeghzadeh, V. Rahneshin Robotic Research Laboratory, Center of Excellence in Experimental Solid Mechanics and Dynamics, School of Mechanical Engineering, Iran University of Science and Technology, Tehran, Iran

a r t i c l e

i n f o

Article history: Received 5 April 2012 Received in revised form 12 May 2012 Accepted 24 May 2012 Available online 26 June 2012 Keywords: Multi-scale modeling Open systems Solid structures Macro/nano mechanics FIMM (Fixed Interfacial Multiscale Method)

a b s t r a c t Fixed Interfacial Multiscale Method (FIMM), a concurrent multiscale model is developed for mechanical analysis of solid nanostructures at the finite temperatures. A direct link between the nano field atoms and macro field nodes using the local atomic volume displacements associated with every macro field node in their common zone has been replaced with the previous methods. FIMM leaves negligible relative motion of atoms in every atomic volume by moving the interface into the macro part. A nano plate in extension and a nanorobotic manipulator have been simulated for comparison. The responses have been compared with the results of semi continuum approach, Molecular Dynamics (MDs), a developed simple Coarse Grained MD (CGMD) and two recently developed multiscale methods. Comparisons show that FIMM is at least 450% better than other methods in numerical cost, and it is also more accurate than other multiscale methods. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Multiscale problems may be observed in a lot of physical, chemical, electrical, mechanical, thermal, and even coupled field systems. Nature reveals this problem as multi size and multi time scales; where in the size scale problems, some ordinary (relatively large) or small (relatively less) dimensions may be included. This leads to various processing time scales for materials and related systems. As listed in Table 1, four characteristic workspaces can be introduced to cover all sizes and time dimensions: (a) The sub-atomic scale (i.e. sub nanometers and sub nanoseconds); in which the quantum dots have the critical role. For this class, the quantum based methods can be employed [1]. (b) The nano scale (i.e. Angstrom to several hundred nanometers and femto to several hundred nanoseconds); where atoms or even electrons play crucial roles and their interactions can be described by an appropriate classical interatomic potential. For material properties at this scale, Molecular Dynamics (MDs) and Monte Carlo (MC) simulations are usually performed [2]. The classical simulations are able to provide valuable insights into atomic processes involving considerably larger systems which can reach up to about 1000 billion atoms [3]. The time scale in MD simulations is up to several hundred nanoseconds. ⇑ Corresponding author. E-mail addresses: [email protected] (M.H. Korayem), [email protected] (S. Sadeghzadeh), [email protected] (V. Rahneshin). 0927-0256/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2012.05.059

(c) The micro scale (i.e. several hundred nano to several hundred micro meters and nano to micro seconds); which includes mesoscale behaviors, too. Although, the larger dimensions in this domain may be simulated with some advanced macro models, the smaller parts still need the atomic based approaches [4]. Protein-based systems, clusters, and related lattice defects such as dislocations, grain boundaries, and other microstructural elements play an important role. Coarse Grained Molecular Dynamics (CGMDs) approach has been developed recently to reach larger sizes and time scales in proteins [5]. This approach has been utilized for general systems with some corrections [6,7]. (d) The macro scale (i.e. 1 mm and more); where the continuum fields such as density, velocity, temperature, displacement, and stress fields are the players. Continuity is the most principal and constitutive law which represents the behavior of the system. The derived continuum equations can be solved analytically or can be discretized and solved numerically by Finite Element (FE) [8], Finite Difference (FD) and other advanced numerical techniques [9]. Nevertheless, FE method (FEM) is the most common approach which has been employed in this regime. When a process occurs only in a specific scale (for example, the motion of a ball on a surface), the aforementioned methods can be used. For a lot of advanced technologies, some processes initiate from one scale and move to another one. Crack initiation and propagation can be mentioned as an example, where it initiates from nano scale and moves to micro and then macro fields. This forces

2

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11 Table 1 A complete span of various models in material modeling.

Length Scale Time Scale Theory→ → Features↓ Structure unit

Sub atomic Sub nano Sec

Nano Nano Sec

Small scale approaches Quantum Mechanics

MD

Phonons

Atom

Micro Micro Sec

Macro Sec Large scale Multiscale approaches approaches Continuum CM-MD CM-SCGMD SCGMD Mechanics coupling coupling (CM) Site Element-Atom Element-Site Element

Unit Shape System Shape Lattice Deformation

Discrete

Time Scale

Below 1pSec

Length Scale

1fm-1A

Applicable to

Sub atomic studies

ContinuousContinuous Discrete 1mSec1fSec-1nSec 1fSec-1μSec 1nSec- 1 mSec 1nSec-1Sec several hours Discrete

1A

Discrete

- 10nm 1fm-100nm

Molecular sciences

Nano to micro Medicine

ContinuousDiscrete

1nm-10μm

1μm-several meters

NEMS, Nanorobotics, Nanorobotics, Macro scale MEMs, manipulation, systems NEMS drug delivery

Macro Super micro Molecular Super nano processes in processes sciences and Not applicable to metric systems long with large sub atomic duration time scales studies

the use of multiscale modeling as a solution for advanced technologies. Furthermore, some other reasons for the urgent need for multiscale modeling can be addressed. Strong dependence of some macroscopic parameters such as temperature and pressure on the small-scale behaviors, characterization of phenomena that have different manifestations at various scales, synchronous study of small and large scales for manipulating the small-scale objects, and direct imaging and strolling in small scales can be mentioned. As an alternative approach, recent efforts have focused on combining aforementioned methodologies together, and some multiscale frameworks have been obtained [10]. Upon introducing the multiscale issue, material and physics scientists found it completely proper as a solution of many unknown problems. In total there are various frameworks to analysis the multiscale problems. For instance, the metacommunity concept in ecology [11] can be addressed. Remarkably, however, much of formal multiscale methods in physics and material science are focused on two categories: sequential and concurrent methods. In the sequential models, the properties are calculated at one scale and then passed on to another scale. In other words, the information obtained from one model is enriched by another model. These approaches include two groups, in which, the information from micro-scale systems is transferred to the continuum mechanics model. The first group are the models based on the Cauchy–Born Hypothesis (CBH) [12,13], in which, the information acquired from the atomic structure of the object clearly reveal themselves in the calculation of the elastic stress and tensor of the material’s properties. The second group of these models is based on the Virtual Atom Cluster (VAC) [14], in which, the structure of the material is enriched on the basis of the information obtained from quantum mechanics. This method has been used for the study of carbon

1nm-1mm

Molecular Small, huge sciences and and non sub atomic homogeneou studies s systems

nano tubes. The characteristic of the systems that are suited for a sequential approach is that the large-scale variations decouple from the small-scale physics, or the large-scale variations appear homogeneous and quasi static from the small-scale point of view. Based on this category, some developments have been directed. For example, Chen and Lee [15] developed a methodology to link atomic to microcontinuum theories through phonon dispersion relations and Zeng et al. [16] extended the applications to nonlocal micromorphic theory. The second category of multiscale simulations is concurrent, parallel or explicit approaches. In the concurrent models, several models exist simultaneously and information is exchanged between them concurrently. This approach is necessary for systems that are inherently multiscale; that is, a system whose behavior at each scale depends strongly on what happens at the other scales. This group of multi-scale models has been used for about a decade, and seven major instances could be cited in the development of these approaches: Quasi Continuum (QC) method [17,18], Macroscopic Atomistic Ab initio Dynamics (MAAD) method [19,20], Coarse-Grained Molecular Dynamics (CGMDs) [5,6], Heterogeneous Multiscale Method (HMM) [21], Bridging-scale Methods (BSCMs) [22,23], Multiscale Field Theory (MFT) [24] and Embedded Statistical Coupling Method (ESCM) [25]. Table 2 lists the general specs of the sequential and all concurrent based methods. In this paper, a Fixed Interfacial Multiscale Method (FIMM) is proposed for computationally and mathematically efficient modeling of solid structures. The approach is applicable to multi-body mechanical systems. In FIMM, a direct link between the nano field atoms and macro field nodes by the local atomic volume displacements associated with every macro field node in their common zone has been replaced with the previous methods. Moreover, con-

Table 2 General specs of the multiscale (sequential and all concurrent based) methods, state of the present method. Methods

Specs QC

MAAD

CGMD

HMM

BSCM

MFT

ESCM

Offline Coupling between scales

Same as Continuum, including small scale energy

Bridges various scales

Piecing a hierarchy of computational approaches None

Continuum Mechanics, especially FEM

Enlarging the unit of small scale (atoms) to several atoms (site) Small scale theories

An incomplete macroscale model, with microscale model as a supplement

Based on

Linking the small and large scales by an interfacial medium (tight-binding (TB) approximation) Crack propagation

Framework of HMM

Compatibility in common zone

Statistical averaging of atomistic displacements in local atomic volumes Restatement of the standard boundary value problem

Elemental selection in small scale and relating to large scale DOFs

Continuity between interfacial medium and both large and small scales

Matching condition at the small-large scale interface

Enforcing DOFs in common nodes

CBH, VAC, . . .

MD/TB/FE, MD/TB, TB/FE

Equality of potential energy of a site and included atoms None

Field description of conservation laws at atomic scale ha Continuous local densities of fundamental quantities Conservation laws

None

Direct Coupling method, VAC

None

FIMM (present)

Mechanical specification study

Crack propagation, . . .

Fast Computing of MD simulations

Complex fluids, solids, interface problems, . . .

Particular class of problems with linear coarse grains in solids

Any type of crystalline material.

General domains, Very large or small systems, . . .

Non Newtonian systems, Plasticity, Multi body dynamics system (may be developed), . . .

Multi body dynamics system, ...

Very large or very small cases

Determining macroscale properties from parametrization of microstructural behavior

Systems with vast common interface, relatively big area for small scale, . . .

Non conservative systems, . . .

Constraints

Included sub methods Application

Not Applicable to

Atomic to microcontinuum link, CBH, VAC Weak dependence between various scales Strongly dependent scales, Complex dynamic systems, . . .

Statistical averages of the atomistic quantities

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

Sequential Definition

3

4

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

sidering the mechanics of the problem, and using a system of equations in matrix form, a dynamic algorithm has been presented in order to solve the multi scale dynamic problem. The macro and nano field’s computational systems are independent of each other and only relate through an iterative update of their boundary conditions. This method presents an improved coupling approach which is inherently applicable to three-dimensional domains. In addition, it prevents the resolving of the continuum model into atomic resolution, and allows finite temperature cases to be applied. One of the prominent features of the present work is the presentation of reliable solutions for problems that include natural, forced, body, and interfacial degrees of freedom. Since solids are fairly rigid at the macro zone, the interfacial volume of the considered system has been moved to macro part and then, it assumed to be rigid. Thus, FIMM leaves negligible relative motion of atoms in each atomic volume by moving the interface into the macro part. Now, previous nano field and a bit of macro part form the new nano field. This leads to larger dimensions for the nano field, in comparison to the previous methods. In Table 2, we placed the FIMM as a sub method of ESCM. One major difference between ESCM and FIMM is the constraint of coupling, where ESCM uses a local average of atoms included in an interfacial volume. In the following parts, both the macro (continuum model) and nano (atomistic model) field theories are discussed briefly. Then, FIMM has been presented in details and validated by comparison with MD, CGMD, ESCM and MAAD approaches for a clamped silicon wafer in plane strain condition. At the end, several useful results and discussions have been introduced. 2. Macro and nano field (MF and NF) theories The FEM is proposed as a solution of MF. In FEM, displacements are expressed as the product of some weighting functions (named shape functions) into the nodal displacements. Substitution of the displacement filed to constitutive relations and then implementation of the Hamilton principle, leads to the governing equation of motion for a finite element as 1 €e þ ½C uu e q_ e þ ½K uu  K u/ K 1 ½Muu e q // K /u e qe ¼ F qe  ½K u/ K // e F /e

ð1Þ

where [Muu]e, [Kuu]e, [Ku/]e, Fqe, [K//]e, F/e, [Cuu]e, and qe are respectively, the element’s mass matrix, stiffness matrix, multi field coupling hardness matrix, mechanical load, secondary field hardness matrix, secondary field force vector, structural damping matrix, and the vector of change of degrees of freedom in the considered system. Usually, the fist field is mechanical and the second may be magnetic, electric and so on. A Simple Coarse Grained Molecular Dynamics (SCGMDs) has been introduced for a solution of NF problem. Consider the general form of MD model as

€ i ¼ ri U mi u

qn1 ;

um1 ;

ðm ¼ anÞ

ð3Þ

where a is the coarsening ratio and T is the transformation matrix. Decomposition of u into the remained (a) and unfolded (b) degrees of freedoms may leads to

 an1 ; bða1Þn1

a ¼ R1 ðq  PbÞ

q ¼ Ra þ Pb;

ð4Þ

The transformation matrix, T, makes the sub matrices, R and P. Considering kinetic (T) and potential (V) energies, generalized forces (Q) and implementing the Lagrange equations, we have

  d @T @T @V   ¼Q dt @ q_ @q @q

ð5Þ

Neglecting the nonlinear terms and spanning the b vector with related part of q vector, it can be equalized to q. Then, the Lagrange equation leads to 1

€  ðR1 kR ðI  PÞ þ kÞq ¼ aQ ðR1 mR1 ðI  PÞ þ mÞq

ð6Þ

‘m’ and ‘k’ are derived from MD model. ‘I’ is the identity matrix and Q is the generalized force in MD model. Above relation has been solved in a discrete form for matrix singularity avoidance. In same thermodynamic conditions, the length and mass scale ratios between MD and SCGMD are equal to ||T|| (norm of q/u) and ||R1R1(I  P) + I|| (norm of M/m). If these two norms are equal to g and f, d SCGMD can be introduced as

rffiffiffiffiffiffi m and Dtmax jCGMD L kT sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi kR1 R1 ðI  PÞ þ Ikm ’g f  kTkL kT

Dtmax jMD

ð7Þ

Practically, values of g and f can be found by the order of coarse graining. Base on the lattice configuration, it may be varied and of course, they depend only on the structure of MD and considered CGMD models. For MD model, Rafii–Tabar–Sutton (RTS) multi-body long-range potential is used in the current study which is an extended form of Sutton–Chen potential with the capability of modeling the unlike material’s interactions. The general form of RTS potential for binary A–B unlike materials is [26–27]:

X 1 XX AA ^i Vðr ij Þ  d p 2 i ji i

HRTS ¼ I

qffiffiffiffiffiffi

qAi  dBB

qffiffiffiffiffiffi X ^i Þ qBi ð1  p

ð8Þ

i

With

^j V AA ðrij Þ þ ð1  p ^i p ^i Þð1  p ^j ÞV BB ðrij Þ þ ½p ^i ð1  p ^j Þ Vðrij Þ ¼ p ^i ÞV AB ðr ij Þ ^j ð1  p þp

qAi ¼

X

X ^j UAA ðr ij Þ þ ð1  p ^j ÞUAB ðr ij Þ UA ðr ij Þ ¼ ½p

j–i

qBi ¼

ð9Þ ð10Þ

j–i

X

UB ðr ij Þ ¼

X ^j ÞUBB ðr ij Þ þ p ^j UAB ðr ij Þ ½ð1  p

j–i

ð11Þ

j–i

^i is the site occupancy operator and defined as: where p

ð2Þ

where m, u and U(u1:uN) are the mass of each particle, displacements, and interatomic potentials. N is the total number of atoms. SCGMD approach is based on the notion that, if instead of one atom, several numbers of atoms are taken as a unit (site), a larger volume of material and also more simulation time can be achieved. Some approaches have been introduced recently [6,7]. The only crucial issue in these models is predicting the interactive potential. Let u and q be the displacements of atoms and sites in MD and SCGMD, respectively. Then a projection may be written as

q ¼ Tu;

 u¼

^i ¼ p



1 if site i is occupied by an A atom 0 if site i is occupied by an B atom

ð12Þ

The functions Vxy(r) and Uxy(r) are defined as

V xy ðrÞ ¼ exy ½axy =rn

Uxy ðrÞ ¼ ½axy =rm

xy

xy

ð13Þ ð14Þ

And the constants are defined by AA

d

¼ eAA C AA

mAB ¼

d

BB

¼ eBB C BB

  1  AA 1 m þ mBB nAB ¼ nAA þ nBB 2 2

ð15Þ

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

aAB

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ aAA aBB

5

pffiffiffiffiffiffiffiffiffiffiffiffiffi

eAB ¼ eAA eBB

where e is a parameter with the dimensions of energy, ‘a’ is a parameter with the dimensions of length and is normally taken to be the equilibrium lattice constant, ‘m’ and ‘n’ are positive constants with n > m. This potential has the advantage that all the parameters can be easily obtained from the Sutton–Chen elemental parameters of metals [27]. MD and CGMD have been solved with a VelocityVerlet (VV) algorithm. The Verlet and Velocity-Verlet schemes require only one force calculation per integration scheme. Verlet algorithms are simplistic and give very good energy conservation characteristics albeit at modest integration step sizes [2]. 3. Macro-micro coupling model In the present method, a direct link between the single NF atoms and MF nodes by means of the local atomic volume displacements associated with every MF node in their common zone has been replaced with the previous methods. The MF and NF computational systems are independent of each other and only relate through an iterative update of their boundary conditions. To generalize the problem, the macro-nano-related problems are divided into two groups of closed and open systems. The group of problems, where all the side boundaries of the nano domain overlap the interfacial degrees of freedom, are called ‘‘closed systems’’; and the group of problems, where the side boundaries of the nano region, in addition to the interfacial degrees of freedom, possess limited (and in some cases, unlimited) degrees of freedom, are called ‘‘open systems’’. Figs. 1 and 2 illustrate the general cases of the closed and open systems, respectively. In the closed system, usually one nano field and one macro field exist. For example, in the crack propagation problem, the fine region of crack propagation is designated as the nano field, and the coarse region in which the crack is growing is designated as the macro field. In the open system, several macro fields could be interacting with several nano fields. If the numbers of macro and nano fields are equal to M and N, respectively, and the area of each field is indicated by X, then for the closed system, we can write:

8 i; j 2 NF > < Xi \ Xj ¼ Ø; Q Xi \ Xj ¼ or Ø; i; j 2 MF > : Xi \ Xj –Ø; i 2 NF; j 2 MF

ð16Þ

Fig. 2. General case definition of open mechanical systems.

Q In the above relations, £ and denote the empty and non-empty spaces, respectively. With this notation, it can be easily proved (considering the presented definitions) that a closed system is a special case of an open system. Therefore, in an open system, there may be more than one nano field and each of the nano fields may be in contact with one another in different ways. These contacts (for example in the nanomanipulation process using nanorobots) may not occur during a certain time range, and after that duration, these contacts may be established. Here, exclusively mechanical systems are taken into account, and therefore, the mentioned contacts are of the second order only, and volumetric sharing is not considered. Through the use of coupling models for closed environments, the size limitation of the atomic model could be minimized, such that an inner region (with complex dynamic processes and large deformation gradients) could exist inside an outer region (with small deformation gradients). It is not like this in open systems, where the effect of size will be considerable. For the coupling of MF and NF in closed systems, four regions are considered throughout the system shown in Fig. 3. These four regions, in the order of proceeding from micro to nano environ-

And for the open system:

8 Q > < Xi \ Xj ¼ Q or Ø; i; j 2 NF Xi \ Xj ¼ or Ø; i; j 2 MF > Q : Xi \ Xj ¼ or Ø; i 2 NF; j 2 MF

Fig. 1. General case definition of closed mechanical systems.

ð17Þ

Fig. 3. Common region between NF and MF in the closed system.

6

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

ments, consist of: Macro Field (MF), Unfolded Volume (UV), Interfacial Volume (IV), and Nano Field (NF). The IV region is in fact, a region where the terminal atoms of a NF model have surrounded a MF node in the model. The IU region is the region between the end nodes of MF and the end of the NF model. The two regions of MF and NF need no further explanation. In view of the presented cases, IVs estimate the mean displacements of NF in the center of mass of these displacements. These averages are later used as the initial conditions of displacements in the relevant interfacial nodes. It should be mentioned that, the IV need not match the macro element that surrounds it, with respect to the size and shape. For the analysis of open systems, in addition to the aforementioned regions, Free Boundaries of nano field (FB) and Common Boundaries of nano field (CB) are also defined. Regardless of the type of initial state, these two regions may have the potential of undergoing different changes during the analysis time range. The CB region is usually affected by the concept of ‘‘cut off radius’’ (Fig. 4). In order for the FB region to behave freely (at the surface), changes should be made to the model. This is the philosophy behind the establishment of the FB region. The existence of free surface creates unwanted effects in the NF system. In comparison with the cases in which the boundary is affected by an external load, this occurrence in FB is not so critical. In addition to unwanted effects, since atoms at the free surface or close to it don’t have a complete set of neighboring atoms, the coordination between the atoms falls apart. To remedy this lack of coordination, and to make the atoms stable in the interfacial NF region, two approaches can be adopted. The first approach is to offer an additional volume of atoms away from the center, which forms the surface NF region. The second approach is to consider a number of the same NF system atoms as an unfolded volume. In case of using the first approach, although the surface NF region eliminates the effects of the free surface, it applies an unwanted virtual stiffness to the system, which elastically constrains the deformation of the inner NF region. To counteract this effect, the unwanted virtual hardness should be compensated. Since the effects of surface in solids are controllable, to a large extent, by the inner volume, in this article, it is suggested to use the second approach. Of course, in places where the limitation of size exists (like the tip of a cone-shaped region), the use of the first approach is inevitable.

tween these two regions. In these schemes, the iterations begin with the displacements of the MF and NF interface. These displacements are obtained as the statistical average from the atomic positions of every IV, and by averaging in the time duration of NF. These average displacements are then applied in the MF region, as displacement boundary conditions (~ dI ). Then, the obtained MF boundary value problem is solved to yield the new interfacial reaction forces, i.e. ~ RI . Then, these forces are applied to the atoms located in IVs; and so, the fixed-reaction boundary conditions are defined in the NF system. During the iterations in which MF is solved, the reaction boundary conditions are fixed, and they are applied to the NF region to guarantee the correct application of the elastic field from the MF domain. In solving the problems of statics, the iteration cycle of NF and MF continues until the system reaches a lasting equilibrium of displacements and forces between the continuum and atomic fields. While for the problems of dynamics, after the establishment of static equilibrium (using the aforementioned method), the system should be solved dynamically (generally via numerical methods). This issue constitutes one of the substantial complexities of the present work.

3.1. Algorithm for establishment of coupling



In general, the coupling of MF and NF is accomplished through schemes based on the establishment of iterative equilibrium be-

3.2. Problems of damped dynamics of MF The general equations of motion (for instance; relation (1)) of the macro part’s displacement, when no external forces exist, are as follows:

LQ¼R

ð18Þ

That, L is motion operator and Q is the vector of displacements and R is external forces. For example, in relation (1), they are as

8 1 d2 d > < L ¼ ½Muu e dt2 þ ½C uu e dt þ ½K uu  K u/ K // K /u e Q ¼ qe > : R ¼ F qe  ½K u/ K 1 // e F /e

ð19Þ

Partitioning the vector of displacements to forced degrees of freedom (F), natural boundary displacements (B) and inner points of the domain (V) and finally interfacial displacements (I), the interfacial forces and displacements can be more achievable. Then, L is broken down as La , in which a = V, F, I, B. The general state of a multi-scale problem includes the initial value and the external force problems. In the first case, the problem is:

La  Q ¼ 0 Qðt 0 Þ ¼ gðtÞ

;

a ¼ V; F; I; B

ð20Þ

The general dynamic displacement vector could be expressed as:

QðtÞ ¼ uðtÞ þ gðtÞ

ð21Þ

where gjðtÞ and u(t) are the vector of initial degrees of freedom and the elastic vector of the whole system, respectively. By substituting the general dynamic displacement vector in the equation of motion, and using the principle of superposition, we have:

La  Q ¼ La  gðtÞ

ð22Þ

And thus, the initial value problem is converted into the external force problem. By arranging the total displacement vector, the state form of the equations can be expressed as:

D  P ¼ AP þ BL:gðtÞ

ð23Þ

where D = d/dt and state matrices are:

" A¼

Fig. 4. Common region between NF and MF in the open system.

0

I

C1 2 C0

C1 2 C1

# ;

" B¼

L0 ¼ C 0 ; L1 ¼ C 1  D; L2 ¼ C 2  D2

0 C1 2

# ; ð24Þ

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

_ T is the state vector. The solution of unknown In which, P ¼ fQQg elastic displacements within the MF region, i.e. u, can be obtained by solving the above state equation and by applying the initial condition IC = 0 in the first step, and applying the condition IC = ICi in the ith step of the macro solution. Then, the interfacial forces are obtained through:

RI ðt n Þ ¼ LI  Q

ð25Þ

3.3. Problem of damped dynamics of NF The dynamics of the atom i with mass m(i), at the position r(i), and in the NF regions are described through the Newton’s equations of motion. Thus, for different regions belonging to the general NF we have:

mi~ f i þ fiD r€i ¼ ~

In the inner NF region

ð26aÞ

~ Rk f i þ Ik þ fiD ; ð~ ri ;~ r_i Þjfirst step mi~ r€i ¼ ~ NI _ ~ CðiÞ ; Q ~ CðiÞ ¼ ðQ Þjit’s last step Interfacial NFðkth IVÞ mi~ fi þ r€i ¼ ~

~ f kcuv NkUV

þ

~ f kAuv NkUV

þ fiD

ð26bÞ

þ

PAtomRC ~k ðf c Þi i¼1

fiD

NkCB In the common NF region ðfor the kth CBÞ

ð26eÞ

The atoms in the inner NF region only experience the atomic force P ~ f ij and the frictional forces of fiD , which result from their f i ¼ j~ neighboring atoms. The atoms existing in the interfacial NF region (which belong to the kth IV) also experience an additional force (~ RkI ) which is distributed among N kI atoms. In addition, the continuity of the fields of the zero- and first-order degrees of freedom should be guaranteed in it. The atoms existing in the surface NF region (which belong to the kth UV) also experience an opposing force (~ f kcuv ) which is distributed among N kS atoms. Moreover, they tolerate the force of ~ f kAuv due to the crushing of the system, which itself should be divided by N kS atoms. In the common boundaries of the NF region, some magnitude of force, which arises from the forces of atoms inside the cut-off radius of the contact surface of two NF fields, should be considered. It should be mentioned that, except in the areas of direct contact between surfaces (where the effect of friction is modeled with the inclusion of some impacts), in the above equations, the viscous friction force fiD is uniformly applied to the atoms inside the IVs and UVs. The complete algorithm for the simulation of coupling has been given in Fig. 5.

Surface NF region ðfor the kth UVÞ ð26cÞ

~ f kcf r€i ¼ fi þ k þ fiD mi~ NFB

mi~ r€i ¼ fi þ

7

In the free NF region ðfor the kth FBÞ ð26dÞ

4. Results and discussions Note that all the coefficients in the SC potential and motion Eqs. (8)–(15) are determined in recent works for molecular problems of metals [26,27]. They are modified for implementation into a simple coarse-grained procedure (3–4, 6–7) to reach more time and length scales. This is one major difference between the SCGMD represen-

Fig. 5. Dynamic analysis algorithm of MF and NF coupling.

8

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

(a)

(b) MD particles Fixed MD particles CGMD sites Fixed CGMD sites FEM nodes FEM elements sides

↑y

→

x

L x = 1nm, L y = 0.6 nm, Thickness = 0.05 nm E = 155.8Gpa, ρ = 0.002329 kg/m3 , Fx = 5 nN

(c)

(d)

Fig. 6. Considered configurations for simulation of tensional deformation of a nano plate: (a) FEM model, (b) MD (blue (small) circles) and CGMD (red (large) circles) models, (c) FIM with MD model and FIM with CGMD model, (d) coordinates and legends of figures. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

FEM

MD

SCGMD

1.5

Normalized Deflection

tation of the problem and MD simulation. The SCGMD coefficient (7) should be determined, in addition to the common MD coefficients. However, the quantity of SCGMD coefficient is limited to a specific range. The SCGMD coefficient may be chosen so that the site size should be no less than the lattice spacing. In the case that it equals lattice spacing, the SCGMD has the same degrees of freedom as MD. For a mechanical system, both MD simulation and SCGMD techniques can be utilized. As a result, the SCGMD based simulation shall be computationally much more efficient for statistical, finite size and finite temperature problems, for simultaneously large length and time scales phenomena, and especially, for dynamic, time-dependent and non-equilibrium phenomena and systems. Leaving the Quantum Mechanics approach, where the subatomic behaviors are studied, other theories are compared together. For simplicity, as depicted in Fig. 6, extension mode of a 1  0.6  0.05 nm silicon wafer has been simulated that is clamped at the left. It has been modeled by FEM (as a CM approach), MD, SCGMD, FIM–MD and FIM–SCGMD approaches. They show the (five) considered configurations and the comparison between the results of these theories, respectively. Before anything else, it should be stressed that the nature of coupling of various domains in multiscale methods makes them very intractable. Therefore, so many simulations should be performed to reach moderately good results. Some parameters should be modified so that they are compatible with existent boundary and ambient conditions. Anyway, organizing a FEM code, a semi macro solution has been obtained and it is normalized using the MD response in amplitude. It is not valid for small scales; so, we used it only for validation of frequency of MD response. An exact solution for the problem of extension may be introduced yet. Fig. 7 shows the results of MD model, the presented CGMD and FEM models for the mentioned nano plate under extension. Deflection and time simulation have been normalized with respect to the t and d, that

1

0.5

0

-0.5 0

0.5

1

1.5

2

Normalized Time Fig. 7. Comparison of the results of the Molecular Dynamics (MDs) model, presented Coarse Grained MD (CGMD) and FEM model for a nano plate under extension.

t ¼

rffiffiffiffiffiffi qL ; E

d ¼

FxL AE

ð27Þ

‘q’, ‘L’, ‘E’ and ‘A’ are the bulk mechanical properties of considered silicon plate and ‘Fx’ is applied force on the right end. Fig. 8 adds the results of FIM model and compares the results of the FIM (FIM–MD and FIM–CGMD) simulation and the continuum model (FEM) and molecular model (MD and CGMD) for a nano plate under extension. Table 3 lists details of presented diagrams. Utilized computational device is an Intel Core i7-2600 K CPU (8 MB L3 Cache, 3.40 GHz). In this table, RK and VV are Rung–Kutta and Velocity-Verlet methods, respectively. Fig. 7 clearly shows that MD and CGMD models, which represent the motion of the atoms at low rate loading, are completely compatible with the semi continuum approach. However, the result of CGMD has more discrepancy compared to the result of

9

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

FEM

MD

SCGMD

FIM-SCGMD

FIM-MD

Normalized Deflection

1.5

1

0.5

0

-0.5 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

Normalized Time Fig. 8. Comparison of the results of the FIM (FIM–MD and FIM–CGMD) simulation and the continuum model (FEM) and molecular model (MD and CGMD) for a nano plate under extension.

FEM. In total, as depicted also in Fig. 8, utilization of CGMD loses some accuracy because it neglects some degrees of freedom in the molecular model to reduce the computational costs. Comparison of the elapsed times in Table 3 vividly shows the effects of this reduction. However, the accuracy should be in attention for less error. The MD and FEM results are approximately the same. However, the result of the coupled model (FIM–MD) shows some discrepancies. This arises from the interfacial bridging algorithm and the resting period from the initiation of loading. As it may be seen, the responses have been more alike after a transient time. 4.1. Application to multi-body nanorobotic systems Many of previous studies in multiscale subject were limited to a specified single dominant body. For instance, crack initiation and propagation of a bulk mechanical system has been considered in

a crack problem or melting point has been searched in a specific material by a thermal study. Based on their multi-body manner, proteins were studied from a multi-body point of view. Perhaps CGMD was a successive resultant of these studies. As mentioned earlier, CGMD has been largely limited to sub micrometric dimensions. In the case of larger systems, multiscale methods are the best solution. For example, in [28], a multiscale simulation study was carried out in order to gain in-depth understanding of machining mechanism of nanometric cutting of single crystal copper. In [29], based on the EAM (embedded atom method), multiscale simulations were performed to study the effects of a nanocavity on nanoindentation of copper film with a hard frictionless indenter by employing the quasicontinuum method (QCM). Here, as a solution in the case of simultaneous intervention of super micro and nano metric domains, the results of FIM have been compared with two recently developed multiscale methods; Coupling Length Scale (CLS) method [20] and the Embedded Statistical Coupling Method (ESCM) [25]. In CLS, surfaces and other regions in which atomistic effects are critical can be modeled accurately with an atomistic simulation. Finite elements offer an adequate and efficient model of micro regions. Then the two, FE and atomistics, are melded together to run concurrently through consistent boundary conditions at the handshaking interface. The main problem is the shrinking of FE mesh until it synchronizes with MD. Nevertheless, this would lead to more elapsed time. On the other hand, ESCM tries to overcome this limitation. It replaces a direct linkage of individual MD atoms and finite element (FE) nodes with a statistical averaging of atomistic displacements in local atomic volumes associated with each FE node in an interface region. The present method may be addressed as a modification of ESCM for fewer computational costs and more accuracy in the case of implementation on solid structures. It is interesting to compare the performance of FIM with aforementioned methods for a multibody system. Scanning Probe Microscopes (SPMs) and especially Atomic Force Microscopes (AFMs) are good instances, where various scales, from nano to

Table 3 Some details for presented comparisons. Models

FEM MD CGMD FIM–MD FIM–CGMD

Specs Averaged response

Frequency (stabilized)

Rising time

Number of DoFs (Large + Small)

Solution method

Elapsed time (s)

Averaged error (toward the MD)

Total rank

0.9911 1.0599 1.2260 1.1607 1.1455

1.22 1.316 1.132 1.105 1.376

3.052 2.959 3.058 2.942 2.604

500 + 0 0 + 1342 0 + 168 250 + 672 250 + 96

RK VV VV RK–VV RK–VV

50 15,720 1820 6003 746

6.94178 – 13.54812 8.684415 7.472719

FEM MD 3 2 1

Fig. 9. General dimensions and arrangement of AFM nano manipulator.

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

macro scales, are in a positive correlation together in order to identify and manipulate the nano objects. Fig. 9 depicts the general dimensions and arrangement of the AFM nano manipulator. A silicon nitride tip has been attached at the end of a micro cantilever that is actuated by a piezo tube. Since the piezotubes (XY piezotube and Z piezotube) dimensions (length ffi 12 mm) are much more than other parts, it is assumed to be rigid compared to them. The larger and smaller parts of the nano manipulator are labeled as MF and NF, respectively. A nanoparticle with 20 nm in diameter has been considered for manipulation. As an example of the affected part of the substrate, a small domain of it has been considered as NF. Fig. 10 depicts the configuration for implementation of FIM on the AFM nano manipulator. Figs. 11 and 12 show the horizontal and vertical travelled distance of nano particle on the substrate. By the aforementioned computational device, the elapsed time of these methods has been compared in Table 3. The total number of sites in CGMD part for nano particle, tip and substrate are 200, 325 and 1047, respectively. The FEM model (of FIM–CGMD, CLS and ESCM) includes 551 nodes and 392 quadratic elements. The total number of atoms in MD part for nano particle, tip and substrate are 1024, 1650, and 5250, respectively. The elapsed time for CLS, ESCM and FIM–CGMD are 156,712, 66054.5 and 12004.7 s. These show that the numerical performance of FIM is 450% more than ESCM, where it has also responsed 1205% better than CLS. This proves the cost efficiency of FIM compared to the other methods, as its responses are also comparable with other results. Figs. 11 and 12 show that FIM optimistically leads to a better response.

Fig. 10. Schematic arrangement for implementation of FIM–CGMD on the AFM nano manipulator.

Horizontal Travelled Distance of Particle (m)

x 10

-8

18 16 14 12 CLS ESCM FIM Desired

10 8 6 4 2 0 0

2

4

6

8

Time (Sec)

10

12

14 -9

x 10

Fig. 11. Comparison of the results of CLS, ESCM and FIM–CGMD; horizontal travelled distance of nano particle on the substrate.

x 10

-8

0

Vertical Travelled Distance of Particle (m)

10

CLS ESCM FIM Desired

-0.5

-1

-1.5

-2

-2.5

0

0.2

0.4

0.6

0.8

Time (Sec)

1

1.2

1.4 x 10

-8

Fig. 12. Comparison of the results of CLS, ESCM and FIM–CGMD; vertical travelled distance of nano particle on the substrate.

5. Conclusion The Fixed Interfacial Multiscale (FIM) method developed in this paper is used to simulate the extension of a thin nano plate (as a closed system) for validation, and then to simulate the AFM nano manipulator (as an open system). The presented approach is proposed for computationally and mathematically efficient modeling of solid structures. As it was completely emphasized, the approach is applicable to multi-body mechanical systems. Since solids are fairly rigid at the macro zone, the interfacial volume of the considered system has been moved to macro part and then, it is assumed to be rigid. Thus, FIMM leaves negligible relative motion of atoms in each atomic volume by moving the interface into the macro part. After some description and model developments, extension mode of a 1  0.6  0.05 nm silicon wafer, which is clamped at the left end, has been simulated. Comparisons show that MD and CGMD models are completely compatible with the semi continuum approach. However, there are more discrepancy between the result of CGMD and FEM. In total, utilization of CGMD loses some accuracy because it neglects some degrees of freedom in the molecular model in order to reduce the computational costs. Comparison of the elapsed times vividly shows the effect of this reduction. The MD and FEM results are approximately the same. However, the result of the coupled model (FIM–MD) shows some discrepancies. This arises from the interfacial bridging algorithm and the resting time since the initiation of loading. As it may be seen, the responses are more alike after a transient time. At the end, as an open system, the AFM nano manipulator has been modeled by FIM and the results have been compared with two other multiscale approaches recently presented. The elapsed time for the implementation of CLS, ESCM and FIM–CGMD to the nano manipulation operation are 156,712, 66,054.5 and 12,004.7 s. These show that the numerical performance of FIM is 450% more than ESCM, where it has also responsed 1205% better than CLS. This proves the cost efficiency of FIM compared to the other methods, as its responses are also comparable with other results. The presented diagrams show that FIM optimistically leads to a better response. Multiscale modeling of materials is in its very early stage of development, and there are many scientific and mathematical problems to be addressed in the future. It is a rich field of physical, numerical, computational, and mathematical challenges. It is also going to play a key role in the simulation and design methodology for the newly emerging field of nanotechnology. FIM may be fo-

M.H. Korayem et al. / Computational Materials Science 63 (2012) 1–11

cused on various existent multiscale problems in the future. It can be readily perform by considering the dynamics of proposed system. References [1] A. Szabo, N.S. Ostlund, Modern Quantum Chemistry, McGraw-Hill, New York, 1989. [2] P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, 1987. [3] M. Kunaseth, R. K. Kalia, A. Nakano, P. Vashishta, Performance modeling, analysis, and optimization of cell-list based molecular dynamics, in: Proceeding of International Conference on Scientific Computing, CSC, 2010, pp. 209–215. [4] T. Zohdi, P. Wrig-Ger, Introduction to computational micromechanics, in: Series: Lecture Notes in Applied and Computational Mechanics, vol. 20, Springer-Verlag, 2005. [5] M. Levitt, Journal of Molecular Biology 59 (1976) 104–107. [6] R.E. Rudd, J.Q. Broughton, Physical Review B 58 (10) (1998) 5893–5896. [7] M.H. Korayem, V. Rahneshin, S. Sadeghzadeh, Computational Materials Science 60 (2012) (2012) 201–211. [8] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1987. [9] M. Bernd (Ed.), Advances in extended and multifield theories for continua, in: Series: Lecture Notes in Applied and Computational Mechanics, vol. 59, Springer-Verleg, 2011. [10] D.D. Vvedensky, Journal of Physics: Condensed Matter 16 (2004) 1537–1576. [11] M.A. Leibold et al., Ecology Letters 7 (2004) 601–613. [12] M.J. Abdolhosseini Qomi, A. Aghaei, A.R. Khoei, International Journal for Numerical Methods in Engineering 85 (7) (2011) 827–846.

11

[13] M.J. Aghaei, International Journal of Solids and Structures 46 (9) (2009) 1925– 1936. [14] D. Qian, R.H. Gondhalekar, International Journal for Multiscale Computational Engineering 2 (2004) 277–289. [15] Y. Chen, J.D. Lee, International Journal of Engineering Science 41 (2003) 871– 886. [16] X. Zeng, Y. Chen, J.D. Lee, International Journal of Engineering Science 44 (2006) 1334–1345. [17] E.B. Tadmor, M. Ortiz, R. Phillips, Philosophical Magazine A 73 (1996) 1529– 1563. [18] V.B. Shenoy, R. Miller, E.B. Tadmor, R. Phillips, M. Ortiz, Physical Review Letters 80 (1998) (1998) 742–745. [19] F. Abraham, J. Broughton, N. Bernstein, E. Kaxiras, Computers in Physics 12 (1998) 538–546. [20] J. Broughton, N. Bernstein, E. Kaxiras, F. Abraham, Physical Review 60 (1999) 2391–2403. [21] B. Engquist, Communication in Mathematical Science 1 (1) (2003) 87–133. [22] G.J. Wagner, G.K. Eduard, W.K. Liu, Computer Methods in Applied Mechanics and Engineering 193 (17–20) (2003) 1579–1601. [23] S.P. Xiao, T. Belytschko, Computer Methods in Applied Mechanics and Engineering 193 (2004) 1645–1669. [24] Y. Chen, J.D. Lee, Philosophical Magazine 85 (2005) 4095–4126. [25] E. Saether, V. Yamakov, E.H. Glaessgen, International Journal for Numerical Methods in Engineering 78 (2009) 1292–1319. [26] A.P. Sutton, J. Chen, Philosophical Magazine 61 (1990) (1990) 139–146. [27] H. Rafii-Tabar, Physics Reports 325 (2000) (2000) 239–310. [28] H.M. Pen, Y.C. Liang, X.C. Luo, Q.S. Bai, S. Goel, J.M. Ritchie, Computational Materials Science 50 (12) (2011) 3431–3441. [29] W. Yu, Sh. Shen, Computational Materials Science 46 (2) (2009) 425–430.