Int J Mech Mater Des (2009) 5:79–110 DOI 10.1007/s10999-008-9087-x
Coupling atomistics and continuum in solids: status, prospects, and challenges J. M. Wernik Æ S. A. Meguid
Received: 19 July 2008 / Accepted: 18 August 2008 / Published online: 14 September 2008 Ó Springer Science+Business Media, B.V. 2008
Abstract The recent focus of the scientific community on multiscale computer modeling techniques of nano-engineered materials stems from the desire to develop more realistic methodologies that are capable of accurately describing the varied time and length scales associated with this class of materials. Of importance is the ability to model the atomistic region using the appropriate techniques such as quantum mechanics/molecular dynamics, and the continuum region using homogenized properties. The continuity of atomistic and continuum regions in a solid necessitates a seamless coupling between these two regions. This is carried out using a transition region. In view of the large discrepancy between length and time scales in atomistic and continuum regions, the development of the transition region has been the main concern of the research community. It is the purpose of this review to critically discuss the issues concerning the transition region and the efforts made by the scientific community in treating them. In particular, this review addresses issues concerning the coupling of molecular dynamics to finite element modeling techniques.
J. M. Wernik S. A. Meguid (&) Engineering Mechanics and Aerospace Design Lab, Department of Mechanical and Industrial Engineering, University of Toronto, 5 King’s College Road, Toronto, ON, Canada M5S 3G8 e-mail:
[email protected] J. M. Wernik e-mail:
[email protected]
Three aspects of this review are accordingly considered. The first is concerned with the current state of atomistic–continuum coupling techniques in computational mechanics. The second is concerned with present the research conducted in the Engineering Mechanics and Design Laboratory at the University of Toronto in the field of nano-reinforced interfaces. Finally, we present the limitations of the current techniques and suggestions for improvements. Keywords Seamless coupling Nano-reinforced composites Sequential Concurrent Multiscale modeling Finite element Molecular dynamics Transition region Atomistics Continuum Carbon nanotubes Nano-reinforced interfaces Nomenclature * ai Acceleration of ith atom * ri Position of ith atom [B] Matrix containing derivatives of shape functions rij Distance between atom i and atom j Ci Set of atoms around representative atom i Rf(t) Bridging scale random force term {d} FE nodal displacement vector t Time [D] Elasticity matrix {T} FE surface traction vector {e} FE strain vector * u Interpolated displacement field
123
80
Ea U E*a Ui Ec * vi Ee * Vðr Þ Ei Vi EQC Vij fa Vijk fi W f(t) u, v, w f imp(t) q* Fi h(t-s) f= g P gc {r} mi l MA w M na [N] Ne Nn Nrep
J. M. Wernik, S. A. Meguid
Atomistic energy FE internal energy Representative atom energy Displacement on node i Continuum energy Velocity of ith atom Elemental energy interatomic potential energy Energy of atom i Individual atom potential Quasicontinuum energy Two-body potential Force on representative atom a Three-body potential EAM electron density dependent embedding energy FE external work Bridging scale interatomic force x, y, z displacements in FE Bridging scale impedance force Mass density Force acting on ith atom Bridging scale time history kernel FE combined body and applied force vector FE potential energy Atomic level force experienced by cluster atom c FE stress vector Mass of ith atom Lennard-Jones potential well depth parameter Atomic mass matrix Lennard-Jones hard sphere radius parameter Finite element mass matrix Representative atom weight function Shape function matrix Number of atoms in element e Number of nodes Number of representative atoms
Abbreviations AFEM Atomic scale finite element method CADD Combined atomistic discrete dislocation method CB Cauchy-Born rule CFRP Carbon fibre reinforced plastic CGMD Coarse grained molecular dynamics
123
CNT DFT EA EAM EANP EANT FE FEAt GFRP GLARE LJ MD MEMS MM MPM NRPC QC QM PDE RVE SWCNT SW MWCNT TB
Carbon nanotube Density functional theory Pure epoxy adhesive Embedded atom method potential Epoxy adhesive reinforced with nanopowder Epoxy adhesive reinforced with carbon nanotubes Finite element Finite element atomistic method Glass fibre reinforced plastic Glass reinforced aluminum composites Lennard-Jones potential Molecular dynamics Micro-electrical-mechanical devices Micromechanics Material point method Nano-reinforced polymer composites Quasicontinuum method Quantum mechanics Partial differential equation Representative volume element Single-walled carbon nanotube Stillinger-Weber potential Multi-walled carbon nanotube Tight binding method
1 Introduction There are a variety of modeling methods currently in use by the research community. They aim not only to simulate material behavior at a particular scale of interest but also to assist in developing new materials with highly desirable properties. These scales can range from the basic atomistic to the much coarser continuum levels. The hierarchy of modeling methods consists of quantum mechanics, molecular dynamics, micromechanics, and finally continuum mechanics. Quantum mechanics is the study of mechanical systems whose dimensions fall on the atomic scale such as that of atoms and electrons. It provides the ability to predict the total energy and the atomic structure of a system of electrons and nuclei. Quantum effects can be best described using either the density functional theory or the tight binding method. Both methods are capable of capturing the physics of the problem at the Angstrom level. However, it is unrealistic to extrapolate the results
Coupling atomistics and continuum in solids
of these methods to the continuum level for engineers to use in system/component design. Similarly, molecular dynamics is a computational method fit for scales in the nanometer range. All of the physics in the molecular dynamics method are contained in the forces acting on each atom in the system. These forces are determined from interatomic potentials. On the other hand, micromechanics modeling techniques employ continuum mechanics techniques to develop field variables for solids containing defects and inhomogeneities at the micrometer level. Finally, continuum mechanics is best suited for scales in the millimeter range where homogenization and averaging techniques are accurate enough to describe the material behavior. Individually, each of these methods is accurate and best suited for its own length scale and modeling inaccuracies may arise from improper enforcement of a specific technique associated with a particular length scale on another. Furthermore, it would not only be impractical but also very costly to attempt to model an entire continuum using molecular dynamics because each individual atom would have to be simulated in the model. The largest atomistic simulations performed to date by the fastest supercomputer are on the order of 1 billion atoms, which only represent a volume of 1 lm3 (Abraham et al. 2002). Many engineering problems are characterized in terms of multiple scales and require a novel approach to describe their behaviors. This is carried out using multiscale computational modeling techniques to describe the behaviors of materials on scales ranging from atomistic to continuum. The mechanical deformation and failure of many engineering materials are in fact multiscale phenomena and the observed macro-scale behavior is governed by processes that occur on many different length and time scales. These processes are often dependant on each other and effect the overall deformation. It is therefore necessary to model these systems using a variety of length scales, which accurately represent the governing physics. The importance of multiscale modeling can be manifested in the following example shown in Fig. 1. The figure demonstrates a typical heterogeneous solid containing microdefects, inhomogeneities, microcracks, and a main crack. The modeling of this heterogeneous solid involves the four simulation techniques mentioned above. As the main crack propagates into the process zone, the atoms at the front of the crack tip
81
Fig. 1 Multiscale example
break their bonds with neighboring atoms. This behavior is accurately represented with quantum mechanics (QM) using the tight-binding method (TB). The accuracy of this approach stems from the fact that ad hoc phenomenological failure models have been avoided. Further away from the crack tip, the material is being deformed but not to the point where bonds are breaking. In this case, the atom displacements and energies are captured using molecular dynamics (MD) (Broughton and Abraham 1999; Rudd and Broughton 2000). The interaction between the main crack and the neighboring microcracks, microdefects, and inhomogeneities can also be modeled using QM and MD. However, this is only the case if the distance between them and the size of the defect or inhomogeniety is on the corresponding scale of interest. Otherwise, the more appropriate method of micromechanics (Meguid and Wang 1934, 1994; Gong and Meguid 1993; Wang and Meguid 1999) can be used to capture the long-range interactions. In the far field, the stress/strain fields are of lower magnitudes and averaged properties across the length scales are sufficient to accurately represent the behavior of the continuum using continuum methods such as the finite element methods (FEM). It is worth noting that the system shown in Fig. 1 is decomposed into seven regions of differing dynamics; as follows: (i) tight binding (TB), (ii) TB-MD interface region,
123
82
(iii) molecular dynamics (MD), (iv) MD-MM interface region, (v) micromechanics (MM), (vi) MM-FE interface region, and (vii) finite element (FE). The problem is further complicated by the fact that the continuum stress and strain fields interact with the other scales of interest. It is crucial that this two-way interaction between the different scales is accounted for in any viable multiscale simulation. Capturing the interaction between the scales of interest is of key concern to multiscale modeling methods and is entirely dependent on the coupling scheme used. There are two main techniques currently adopted in the treatment of this class of problems. The first is the atomistic-based continuum technique, in which the atomistic region is smeared and included in a continuum model that describes the entire field of atomistic and continuum. In this case, the problem can be subdivided into both atomistic and continuum domains although the overall model is averaged into one continuum representation. Both sub-domains contain their own constitutive laws but the governing relations for the problem as a whole are obtained through the homogenization of the resulting continuum. The atomistic-based continuum approach will not be further considered in this review. The second, which is discrete in nature, relies on coupling atomistic solutions, using for example MD, with discrete models using the finite element method. The discrete methods can be further subdivided into sequential and concurrent coupling methods. Sequential methods, sometimes referred to as hierarchical methods, pass information (displacement or force fields) from the finer scale as boundary conditions to the coarser one. The sequential approach assumes that the problem considered can be easily separated into processes that are governed by different length and time scales. Therefore, the coarse-scale physics are completely determined by the fine scale along with the corresponding boundary conditions. Many physical phenomena can be readily modeled using single scale methods where the physics of the coarser scale are condensed into few parameters. Further details regarding sequential coupling will be presented in Sect. 3. Concurrent methods, on the other hand, perform the entire multiscale simulation simultaneously and continually feed information from one length scale to the other in a seamless fashion. Concurrent methods are better suited in representing scales with a heavy dependence on each other
123
J. M. Wernik, S. A. Meguid
because of the continuous transfer of information between the different scales. In essence, the passing of information ensures consistency amongst the field variables between the two simulation methods. Researchers have developed a wide variety of different concurrent methods that have been applied to problems ranging from the modeling of MEMS devices (Rudd and Broughton 1999; Rudd 2000) to nanoindentation material response (Tadmor et al. 1996a, 1999). In this review, our focus is on discrete coupling methods which include both sequential and concurrent methods. Several excellent reviews have been published which address the methods used to couple the atomistics to the continuum. These reviews include those by Liu et al. (2004a, b), Park and Liu (2004), Rudd and Broughton (2000), Vvedensky (2004), and Curtin and Miller (2003). An excellent review has also been published on single scale modeling methods by Ghoniem et al. (2003). The objectives of this review, however, are threefold. In the first, we describe the current state of literature in the field of atomistic/continuum modeling, while in the second, we present the experimental and theoretical research currently being conducted by the authors in the field of nano-reinforced interfaces. In the third, we discuss the limitations of these existing techniques. The layout of this review is as follows. Following this introduction section, we will describe the fundamental formulations pertaining to molecular dynamics and continuum mechanics in Sect. 2. In Sect. 3, we will discuss the sequential and concurrent coupling approaches with a predominant focus on concurrent coupling methods. In Sect. 4, we will describe the research currently being conducted at the Engineering Mechanics and Aircraft Design Laboratory at the University of Toronto in the area of nano-reinforced interfaces. Finally, in Sect. 5, we will outline some of the limitations of the coupling efforts as well as suggest possible enhancements.
2 Fundamental formulations 2.1 Molecular dynamics Molecular dynamics is a computational method used to simulate the time dependent behavior of molecular
Coupling atomistics and continuum in solids
83
systems. MD simulations provide a means of investigating the structure, dynamics, and thermodynamics of individual molecules. The simulations generate trajectories that describe the atomic positions, velocities, and accelerations of individual particles as they vary with time. These trajectories are then used to obtain average values of system properties such as energy, pressure and temperature. At the core of the molecular dynamics method is the process of determining the forces acting on individual atoms. These forces are then used to evolve the simulations in time by incorporating them into Newton’s second law, as follows: *
*
F i ¼ m i ai
ð1Þ
*
where F i is the force acting on particle i, mi is the * mass of particle i, and ai is its acceleration. The atomic forces are determined by taking the gradient * of the systems total potential energy, Vðr Þ; as given by *
*
F i ¼ ri Vðr Þ
ð2Þ
where V(r) is the system’s position-dependent interatomic potential. It represents the potential energy of the system when the atoms are arranged in a specific configuration. Combining Eqs. 1 and 2 *
*
ri Vðr Þ ¼ mi ai
ð3Þ
illustrates that the derivative of the potential energy with respect to the particle position can be related to the acceleration of that particle. Only the initial positions, velocities, and an expression for the total system potential energy need to be provided to conduct a molecular dynamic simulation. In most MD simulations, the initial atomic positions are determined by the material and the problem being simulated, whereby they are chosen to coincide with the atomic lattice locations. The initial velocities are generated randomly based on statistical mechanics at a specified temperature. The other quantity of interest needed to conduct an MD simulation is an expression for the total potential energy. The total potential energy of an atomic system is described using classical interatomic potentials that can be expressed in a wide variety of forms. The general expression for total atomistic potential energy can be expressed as
*
*
*
Vðr 1 ; r 2 ; . . .; r N Þ ¼
X i
þ
*
Vi ðr i Þ þ X
X
*
*
Vij ðr i ; r j Þ
i;j;j [ i *
*
*
Vijk ðr i ; r j ; r k Þ þ
i;j;k;k [ j [ i
ð4Þ where N is the number of atoms in the system. The functions Vm are called m-body potentials and describe the interaction of one to m pairs of atoms. The first term represents the energy of the isolated atom i. The second accounts for two-body interactions of the atoms i and j, and the third for the threebody interactions, and so on. In the field of multiscale modeling, many different interatomic potentials have been used. This is primarily because of the large variance in the materials being studied. The choice of potential for a molecular dynamics simulation is determined by factors such as the material being modeled, the bond type, the desired accuracy, and available computational resources. Here we restrict our attention to the Lennard-Jones (LJ), embedded-atom (EAM), and Stillinger-Weber (SW) potentials for the purpose of demonstrating their applicability and differences. The Lennard-Jones potential is termed a pair potential. Pair potentials only account for the energies associated with the relative displacement of two atoms and can be an oversimplification in some cases. The Lennard-Jones (Jones 1924) potential is defined as " # 1 1 Vij ¼ 4l ð5Þ 12 6 rij =w rij =w where l is the potential well depth, w is the hard sphere radius of the atom or the distance at which Vij is zero, and rij is the distance between the two atoms i and j. The 12th power term is meant to model the repulsion between two atoms as they approach each other. The 6th power term adds cohesion to the system, and is intended to mimic van der Waals type forces. The van der Waals forces are weaker than the repulsion term hence its lower order. The LJ potential is used for large-scale simulations where computational efficiency is highly desirable and the focus is on fundamental issues rather than the study of specific material properties. The LJ potential can be considered accurate for materials that have no electrons
123
84
J. M. Wernik, S. A. Meguid
available for bonding and only interact through van der Waals forces. The EAM and SW types of interatomic potentials are termed multi-body potentials. These types of potentials not only account for the relative pair-wise displacement energies, but can also include energies associated with electron densities, bond bending, torsion, and other effects which are neglected in pair potentials. The EAM potential (Daw and Baskes 1984) is a common choice for metals. The EAM constitutive law involves an energy associated with embedding the atom into a local electron density of the surrounding atoms. The EAM energy of an atom is given by 1X Ei ¼ fi ð qi Þ þ Vij ðrij Þ ð6Þ 2 j6¼i where fi is an electron density dependent embedding energy and Vij is a pair potential between atom i and its neighbor j. The electron density at atom i, qi, is the superposition of density contributions from each of the neighbors such that X i ¼ q qj ðrij Þ ð7Þ j6¼i
The SW type (Stillinger and Weber 1985) energies involve both two- and three-body interaction terms. The SW potentials are normally used for covalent materials, semiconductors and when multi-body interactions are important. The main difference between the EAM and SW type atomistic laws is that instead of an embedding energy, SW includes a term involvingthree body interactions to account for the directional bonding in covalent materials. The SW energy of an atom i is 1X 1 X * * * Ei ¼ Vij ðrij Þ þ Vijk ðr i ; r j ; r k Þ ð8Þ 2 j6¼i 6 i;j;k;k [ j [ i where Vijk is the three-body interaction potential. Nearly all interatomic potentials used in MD simulations incorporate a cut-off distance whereby interactions between atoms that extend beyond this distance are ignored. The purpose of this truncation is to reduce the number of degrees of freedom that need to be included in determining the energy associated with a particular atom. This also lends to the fact that the energy in an MD simulation is non-local. This effectively means that the energy in a bond between a pair of atoms depends not only on the position of the
123
immediate neighbors but also with atoms within a specified radius. This radius is always chosen as the cut-off distance of the interatomic potential being used. The total potential energy is a function of the atomic positions of all the atoms in the system. Due to the complexity of such a function, the resulting equations of motion are integrated numerically. Different numerical algorithms have been employed in MD simulations for integrating the equations of motion. The most popular of which is the Verlet algorithm (Verlet 1967). All the algorithms assume that the positions, velocities and accelerations can be approximated by a Taylor series expansion as given by 1 rðt þ DtÞ ¼ rðtÞ þ vðtÞDt þ aðtÞDt2 þ 2 vðt þ DtÞ ¼ vðtÞ þ aðtÞDt þ
1 d3 rðtÞ 2 Dt þ 2 dt3
ð9Þ ð10Þ
and a similar expression for the acceleration. In the above Dt is the timestep used throughout the simulation. The verlet algorithm proceeds by only writing two third-order Taylor series expansions of the position, one forward in time and the other backward in time as follows: 1 rðt þ DtÞ ¼ rðtÞ þ vðtÞDt þ aðtÞDt2 2 1 d 3 rðtÞ 3 þ Dt 6 dt3 1 rðt DtÞ ¼ rðtÞ vðtÞDt þ aðtÞDt2 3 2 1 d rðtÞ 3 Dt 6 dt3
ð11Þ
ð12Þ
Then by combining the above two equations, we obtain rðt þ DtÞ ¼ 2rðtÞ rðt DtÞ þ aðtÞDt2
ð13Þ
where a(t) is expressed as the force divided by the mass, and the force itself is a function of the positions as given by Eq. 2. This is the general form of the verlet algorithm. This method is time-reversible, has good energy conservation properties and has low memory storage because the velocities are not included in the time integration. However, removing the velocity from the integration introduces numerical inaccuracies. This has lead to the development of other algorithms that are based off of the Verlet
Coupling atomistics and continuum in solids
85
method, namely the velocity Verlet (Allen and Tildesley 1987) and Verlet leap-frog (Hockney 1970) algorithms which explicitly involve the velocity in the time evolution of the atomic coordinates. We refer the reader to three excellent reviews that address the molecular dynamics method in its entirety that have been published by Rafii-Tabar (2000), Carlsson (1990), and Ercolessi (1997) for further reading. 2.2 Finite element method The traditional framework in mechanics has always been the continuum. Under this framework, materials are assumed to be composed of a divisible continuous medium, with a constitutive relation that remains the same for a wide range of system sizes. Continuum equations are typically in the form of deterministic or stochastic partial differential equations (PDE’s). The underlying atomic structure of matter is neglected altogether and is replaced with a continuous and differentiable mass density. Similar replacements are made for other physical quantities such as energy and momentum. Differential equations are then formulated from basic physical principles, such as the conservation of energy or momentum. There are a large variety of numerical methods that can be used for solving continuum partial differential equations, the most popular being the finite element method (FEM). The finite element method is a numerical method for approximating a solution to a system of partial differential equations. The FEM proceeds by dividing the continuum into a number of elements, each connected to the next by nodes. This discretization process converts the PDE’s into a set of coupled ordinary equations that are solved at the nodes of the FE mesh and interpolated throughout the interior of the elements using shape functions. The main advantages of the FEM are its flexibility in geometry, refinement, and loading conditions. It should be noted that the FEM is local, which means that the energy within a body does not change throughout each element and only depends on the energy of the nodes of that element. The total potential energy under the FE framework consists of two parts, the internal energy U and external work W P¼UW
ð14Þ
The internal energy is the strain energy caused by deformation of the body and can be written as
U¼
1 2
Z
frgT fegdX ¼
X
1 2
Z
fegT ½DfegdX
ð15Þ
X
T where frg ¼ rxx ryyrzz sxy syz sxz denotes the T stress vector, feg ¼ exx eyy ezz cxy cyz cxz denotes the strain vector, [D] is the elasticity matrix, and X indicates that integration is over the entire domain. The external work can be written as 8 9 Z < =x = W ¼ fu v wg =y dX : ; = X 8z 9 Z < Tx = ð16Þ þ fu v wg Ty dC : ; T z C where u, v, and w represent the displacements in x, y, z directions, respectively, {=} is the force vector which contains both applied and body forces, {T} is the surface traction vector, and C indicates that the integration of the traction occurs only over the surface of the body. After discretization of the region into a number of elements, point-wise discretization of the displacements u, v, and w for the respective x, y, and z directions, is achieved using shape functions [N] for each element, such that the total potential energy becomes Z 1 P g ¼ fd gT ½BT ½D½BfdgdX 2 Xg 8 9 Z < =x = T T fd g ½N = dX : y; = Xg 8 z9 Z < Tx = fd gT ð17Þ ½N T Ty dC : ; T g z C where [B] is the matrix containing the derivatives of the shape functions, and {d} is a vector containing the displacements. A full description of the finite element method can be found in either of the texts by Hughes (1987) or Zienkiewicz (1991).
3 Coupling methods in nanomechanics In this section we will introduce both sequential and concurrent coupling methods with emphasis on the latter. It is important to distinguish between the two approaches and identify how they proceed in coupling the scales of interest.
123
86
3.1 Sequential coupling Sequential coupling methods, as their name suggests, perform the simulations in step-by-step fashion in a bottom-up approach. In doing so, the simulations are running independently of each other and a complete separation of both length and time scales is achieved. Sequential methods employ parameter passing between the different simulation methods. The outputs of the finer scale simulation are used as boundary conditions for the coarser scale. These boundary conditions can consist of displacement or force fields. Due to the fact that the simulations are not simultaneously, it is important to feedback the force fields after each step and compare them with the previous ones in order to ensure that the simulations are converging. However, it should be noted that the uniqueness of the solution is not guaranteed. The number of scales that can be coupled is entirely dependent on the problem that is being modeled. A schematic of the sequential approach is given in Fig. 2. The problem is divided into two scales: MD modeling and FE modeling. The figure also illustrates the complete separation of the scales and the process of parameter passing between each simulation. The necessity of feeding back the force fields for comparison is also illustrated in the figure. The complete separation of the simulations may lead to the loss of information transfer between the fine and coarser regions. This problem may become more pronounced as more scales are coupled together. The sequential
Fig. 2 Schematic of the sequential approach
123
J. M. Wernik, S. A. Meguid
coupling approach is best suited for coupling scales with weak dependence. As the scales of interest begin to show a stronger dependence, the one time passing of boundary conditions may be insufficient to ensure accuracy or realistic results. This will necessitate an iterative passing of information between scales inherent in concurrent multiscale methods. Some examples of sequential coupling studies available in literature include the study of tidal basins starting from the Angstrom scale of water molecules by Clementi and Reddaway (Clementi 1988). In his research, an interatomic potential was first derived from first principle quantum mechanical calculations and passed to a molecular dynamic simulation used to obtain the densities and velocities of the individual water molecules. These parameters were then employed in a fluid dynamics calculation to model the tidal circulation. Hao et al. (2003) developed another prime example of a sequential model in the design of a ductile fracture simulator of high toughness steels. This model transcended scales from the quantum to continuum with laws developed by each scale that govern the behavior of the coarser scale. First principle calculations derive the force separation laws between atoms while coarser scale simulations, such as finite element methods, provide a constitutive model for crack propagation. The analysis of carbon nanotube functionalization and deformation has also been an area of active research for sequential modeling approaches (Odegard and Frankland 2005). Namilae and Chandra (2005) developed a sequential model to study the atomic scale interface effects on composite behavior. In their study, they used molecular dynamics to simulate fiber pullout tests of carbon nanotubes embedded in a matrix. The MD results were used to evaluate the cohesive zone parameters which were then used as boundary conditions for the FE simulations to obtain the macroscopic response of the composite. It has, however, been identified that under some circumstances the characterization of carbon nanotube deformation can accurately be modeled using purely continuum concepts (Yakobson et al. 1996). Carbon nanotubes, which have diameters on the order of 1 nm, behave much like slender hollow cylinders under bending and seem to display continuum characteristics such as the formation of kinks as strain energy release mechanisms and a constitutive behavior similar to Hooke’s law. In this special
Coupling atomistics and continuum in solids
case, quantum (Zhao and Ding 2008; Galano and Francisco-Marquez 2008; An et al. 2007; Troya et al. 2003; Wang et al. 2006; Ding 2005), molecular (Cao and Chen 2006a, b, 2007; Esfarjani et al. 2006; Zhao and Cummings 2006; Gong et al. 2008; Buehler et al. 2004; Srivastava and Wei 2003; Mylvaganam et al. 2006), and continuum (Yakobson et al. 1996; Cao and Chen 2006b; Pantano et al. 2004) models can all provide comparable and accurate results. Additional sequential coupling method studies can be found in Cavallotti et al. (2005), Haslam et al. (2002), Almeida et al. (2008). 3.2 Concurrent coupling Of importance to this section is the coupling of the atomistic and continuum regions to develop a unified and concurrent treatment of a physical system. In this case, we must introduce a transition region to couple the different length scales associated with molecular dynamic and finite element simulations as well as allow two way transfer of the field variables of interest between these two regions. The literature contains numerous methods of concurrent coupling. In this review, we limit our attention to the following; (i) the combined finite element atomistic method (FEAt), (ii) the material point method (MPM), (iii) the local quasicontinuum method (QC), (iv) the non-local quasicontinuum method, (v) the bridging scale method, (vi) combined atomistic discrete dislocation method (CADD), (vii) the atomic-scale finite element method (AFEM), and (viii) coarse grained
87
molecular dynamics (CGMD). There are other methods which will be omitted from discussion either because the focus was on scales of different interest or the general approach is very similar to the methods discussed above. The additional coupling methods can be viewed in the literature (Broughton and Abraham 1999; Abraham and Broughton 1998; Abraham 2000; Fish et al. 2007). 3.2.1 The MD/FE transition region Most MD/FE coupling methods incorporate a region where the MD and FE domains interact and transfer boundary conditions to ensure equilibrium and compatibility of these domains. This region is normally referred to as the transition region or ‘handshake’ zone, and its general form is depicted in Fig. 3. The treatment of the material in the transition region between the atomistic and continuum regions is a critical aspect of the problem. It is in this region that most of the approximations are made due to the inherent incompatibility of the two domains; namely, the non-local formulation of the molecular dynamics region and the local formulation of the finite element region. Several features are common amongst all concurrent coupling methods which employ this type of transition region. First, one side of the transition region is neighbored by the MD region and the other by the FE region. The transition region consists of an overlap of the two domains where both atoms and finite element nodes are present. It has been a standard practice to reduce the finite element node spacing to
Fig. 3 General depiction of the transition region used in coupling methods
123
88
coincide with the atomic spacing so as to make a smooth and seamless transition between the two domains. This results in a one-to-one correspondence between the atoms and nodes directly along the interface between these two regions; this effectively means that the nodes directly overlap the atomic positions. This eliminates the need to interpolate the atom positions throughout the elements at the interface as well as to ensure continuity of the displacements. The atoms and nodes share the same position. As one moves away from the interface to the continuum direction, the finite elements become larger and the nodes become increasingly sparse. This reduces the number of degrees of freedom of the problem by representing the areas of little interest with a coarser continuum representation. As in standard finite element methods, there is no need to limit the choice of elements used in the model to one specific type. In the case of Fig. 3, four-noded quadrilateral elements are used in the immediate vicinity of the interface. As one moves further into the continuum domain, eight-noded quadrilateral elements are used. The choice of element type, as with many other factors, is dictated by the desire to balance between accuracy and computational cost. The overlap atoms in the transition region serve the purpose of providing the MD atoms along the interface with a complete set of neighbors for them to interact with according to the non-local formulation of the MD region. They typically extend out a distance equal to the cutoff distance of the interatomic potential used. The overlap atoms are not true atoms in the fact that they do not contribute their energies to the system but only contribute their effects on the MD atoms. Their positions are determined through interpolating the displacements from the nodes of the elements in which they reside. Weight functions are introduced at the interface to ensure continuity of the field variables (forces and displacements). These weight functions vary according to the type of elements being used, the number of nodes of each element that fall on the interface, and in the case of the forces, the nodal separation distance. 3.2.2 Finite element-atomistic method (FEAt) One of the earliest concurrent coupling methods developed was the combined finite element and atomistic method (FEAt) by Kohlhoff et al. (1991).
123
J. M. Wernik, S. A. Meguid
This method was applied to the study of crack propagation on cleavage and non-cleavage planes in b.c.c. crystals as well as brittle fracture (Kohlhoff et al. 1991; Gumbsch and Beltz 1995; Gumbsch 1995). This coupling method is capable of treating plane static and dynamic problems. It uses non-local elasticity theory to describe the transition region between the MD and FE domains. It attempts to address the non-local/local mismatch between the two regions by using a non-local continuum formulation. This approach will automatically correct for unrealistic forces called ‘ghost forces’ that occur due to the locality mismatch, which will be discussed in the following section. In their efforts to couple the MD and FE domains, Kohlhoff and coworkers established a one-to-one correspondence between atoms and nodes throughout the entire transition region. In this case, the transition region extends a distance of twice the cutoff distance used by the interatomic potential. This arrangement ensures compatibility between the two domains at the interface. However, equilibrium between the lattice and the continuum is not guaranteed a priori, and additional conditions are necessary. These conditions are imposed on the elastic constants by expanding the elastic energy in terms of Taylor series with respect to the strains, viz; EðeÞ ¼ Eð0Þ þ oE=oeij 0 eij 1 þ o2 E= oeij oekl eij ekl 2 0 1 3 þ o E= oeij oekl oemn 0 eij ekl emn þ 6 ð18Þ Since the strains are assumed equal at the interface, equilibrium of the stresses implies that the elastic constants in both the atomistic and continuum domains should be equal. It is therefore assumed that the coefficients in the series are equal to the elastic constants; viz: cij ¼ oE=oeij j0 ð19Þ 2 ð20Þ cijkl ¼ o E= oeij oekl 0 3 cijklmn ¼ o E= oeij oekl oemn 0 ð21Þ Terms up to and including the third order are matched to ensure equilibrium. The fact that the Taylor expansion has to be cut off at some point introduces one of the approximations of this method.
Coupling atomistics and continuum in solids
3.2.3 Material point method (MPM) The material point method is an alternative to FE for representing the continuum domain. This method was developed by Sulsky et al. (1995) as a particle-in-cell method. It was predominantly used to model fluid flow but has been recently extended to the multiscale modeling of silicon under tension by Lu et al. (2006a) and to high energy cluster impacts by Guo and Yang (2006). The MPM is relatively new to this field, and can be used either in two- or three-dimensional simulations. The material point method is a finite-temperature dynamic MD coupling method that uses two representations of the continuum; one based on a computational mesh and the other based on a collection of material ‘points’ or particles. These two representations form the Eulerian and Lagrangian descriptions of the material, respectively. The structure of the method is depicted in Fig. 4. The MPM method makes use of the benefits of each description while avoiding the associated problems with each. Traditionally, the finite element method is used as a representation of the continuum. However, if the material is seriously distorted the FEM will typically result in severe mesh distortion and the consequence is ill conditioning of the element stiffness matrix leading to mesh lockup or
89
entanglement. This is avoided in the MPM in the way it handles the deformations. In the MPM, the material continuum is discretized into a collection of material points whereby each material point is assigned a mass based on the volume it covers and the initial mass density. It also carries all relevant physical characteristics, such as position, velocity, acceleration, stress, strain and constitutive parameters. These material points form the Langrangian description of the material. In the continuum framework, motion is governed by the conservation equations including mass, momentum and energy. In the MPM, the conservation of momentum is solved on a background Eulerian mesh to avoid the meshlockup associated with the Langrangian description. The FE shape functions are adopted to map information between material points and background grids. The MPM incorporates a hierarchical mesh refinement technique that allows the FE mesh to be reduced to the atomistic scale thereby forming a oneto-one correspondence between material points and MD atoms. A computational grid is constructed from 4- and 5-noded isoparametric quadrilateral elements for two-dimensional simulations. These elements are then used to define the shape functions associated with the spatial nodes. The coupling scheme used employs a similar transition region as the one depicted in Fig. 3. The only difference is that the MD atoms overlap the material points as apposed to the FE mesh nodes. It ensures compatibility of the strain fields of the system, just as in the FEAt method, but it also does not guarantee equilibrium of the lattice and continuum. Therefore, the MPM method also follows the same approach of matching the elastic constants to the coefficients of the third order Taylor series expansion of the interatomic potential. Thus, compatibility is achieved as long as the elastic constants of the continuum are equal to those defined by the interatomic potential in the MD region. 3.2.4 The local quasicontinuum method
Fig. 4 MPM structure
The quasicontinuum method was first conceived by Tadmor et al. (1996b) at Brown University and has since had a great deal of development and has been applied to a number of different applications. The method was first applied to modeling defects in single crystal FCC metals and used to study nanoindentation
123
90
in thin films (Tadmor et al. 1996a). Shenoy et al. (1998, 1999) generalized the method and extended it to treat the interactions of dislocations with grain boundaries in polycrystalline solids. Miller et al. (1998a, b) used the method to study fracture at the atomic scale and investigated the effect of grain orientation on fracture and the interaction of cracks with grain boundaries. Knap and Ortiz (2001) have developed an entirely non-local three-dimensional version of the model and applied it to the study of nanoindentation. More recently, a method of coupling the atomistics with discrete dislocation dynamics using the QC method was developed by Qu et al. (2005). A recent review of the QC method can be found in Miller and Tadmor (2002). The QC method is a static, zero-temperature coupling method. Efforts are being made to extend it to account for finitetemperature (Dupuy et al. 2005), and to develop a dynamic version. The most recent study employing the quasicontinuum method was performed by Chandraseker and Mukherjee (2007) in which he used it to estimate the elastic properties of singlewalled carbon nanotubes (SWNT) subject to coupled extension and twist deformations. An important feature of the QC method is that it allows the fully atomistic regions to evolve with the deformation during simulation using an adaptive mesh refinement procedure similar to the ones employed in standard finite element methods. The QC method approaches the issue of coupling by determining a total energy functional for the entire system and enforcing its equilibrium by minimizing this energy with respect to the atomic positions. As with all other coupling methods, the goal of the QC method is to drastically reduce the degrees of freedom and computational cost of the problem, while maintaining full atomistic details in regions where it is required. This is achieved through the implementation of two approximations: (i) the use of representative atoms, and (ii) the implementation of the Cauchy-Born (CB) rule to map continuum scale deformation fields to the atomistic scale. In (i), the formulations introduce a measure of the magnitude of deformation using the deformation gradient F. This stipulates that if the deformation gradient changes gradually at the atomic scale, then it is not necessary to track the displacement of every atom in the system. In this case, a small subset of atoms is chosen to capture this gradual change. This particular subset of atoms is referred to as
123
J. M. Wernik, S. A. Meguid
representative atoms. The representative atoms act as atoms in the atomistic region and also serve the purpose of acting as nodes in the continuum region. The positions of the representative atoms can be determined from FE methods and the positions of the atoms located within the elements can be determined from using interpolating functions. The implementation of the representative atoms in this manner will reduce the degrees of freedom of the problem to only those of the representative atoms. The selection of representative atoms is determined from the magnitude of the change in the deformation gradient. If the gradient changes drastically, this will indicate that there is a need for greater details in this region in which case more of the atomic positions are chosen as representative atom locations. However, if the deformation gradient changes slowly, this indicates that there is little or no deformation occurring in this region and that fewer representative atoms may be selected to describe this region. In essence, the representative atom density should be high in regions of large deformation. Figure 5 illustrates this process for three cases of varying deformation gradient with application to an edge dislocation example. In the first two cases the deformation gradient is either constant or changing slowly over a prescribed distance of Dx. In these cases it is sufficient to capture the deformation using few atoms chosen as representative atoms. The third case, however, shows a large variance in the deformation gradient over the same distance in the immediate
Fig. 5 Representative atom selection based on deformation gradient
Coupling atomistics and continuum in solids
91
vicinity of the dislocation core. In order to accurately capture this change more atoms will be required. Once the representative atoms have been identified the FE mesh is constructed by using them to define the nodes. In regions where all atoms are chosen as representative atoms, the dislocation core, the simulation reduces to a purely atomistic model while the other areas maintain a coarser FE representation. For the purpose of energy minimization, the energy of each atom must still be determined to construct the total energy functional. This brings us to item (ii) stated above. The use of linear interpolating functions to determine displacements within the elements effectively means that the deformation gradient will be uniform throughout the elements. The QC method makes use of the Cauchy-Born rule which stipulates that a uniform deformation gradient at the continuum scale can be mapped to the same uniform deformation in the atomistic scale as depicted in Fig. 6. This means that every atom in an element will have the same energy. Therefore, the energy of any element in the FE region can be determined by multiplying the energy of one atom by the number of atoms in the element; as follows: *
*
Ee ðuÞ ¼ Ne Ei ðuÞ
EQC ¼ Ea þ Ec ¼
X i2ðA;IÞ
*
Ei ðuÞ þ
NX element
*
Xe Ee ðuÞ
ð23Þ
e¼1
where Ea is the total energy in the atomistic region and Ec is the total energy in the continuum region. The atomistic energy is simply the sum of all the individual atom energies in the model. The summation includes both the atoms in the MD region and those that lie along the interface between the MD and FE regions. The continuum energy is a summation of all the element energies where Nelement is the number of elements in the model. To ensure appropriate force contributions from both regions for the interface atoms/nodes, the energies of the continuum elements directly adjacent to the interface are weighted differently in the total potential energy using the weight functions xe. The QC method then employs the conjugate gradient method to drive the system towards equilibrium. The objective here is to determine the zero force positions of the atoms in the model by minimizing the total potential energy function with respect to the atomic positions.
ð22Þ
where Ne is the number of atoms in element e, Ei is * the energy of one atom in element e and u is the interpolated displacement field. Now that the energy formulation, for both the atomistic and continuum
Fig. 6 Cauchy Born rule depiction
regions, has been developed, the total energy for the entire coupled system can be described as follows:
3.2.5 The non-local quasicontinuum method When applying the Cauchy-Born rule in the above formulation of the QC method it is assumed that the deformation is uniform throughout the elements. The uniform deformation constraint, by the way it has been applied, is only valid within the elements and not along the interfaces or free surfaces. To correctly account for the energy of the interfaces between elements and surfaces, the non-uniform surroundings of the atoms along the interface must be correctly represented. For example, a representative atom that resides along the interface between two elements will be subject to a deformation gradient that will be different in both elements. Therefore, the CB rule cannot be applied to this atom. In order to accurately account for this variation in the deformation gradient, the continuum must be represented in a non-local framework (Knap and Ortiz 2001). There are two ways of implementing non-locality in the continuum region; the energy-based and force-based formulations. The energy-based formulation begins by modifying the continuum energy functional to include only
123
92
J. M. Wernik, S. A. Meguid
representative atom energies. Each representative atom can be considered as a cluster of atoms, the size of which depends on the location of the representative atom in the system. In this case, the energy functional will include a weighting scheme that corresponds to the number of atoms in each cluster. Under this formulation the continuum energy is defined as the sum of all the representative atom energies as given by Ec ¼
Nrep X
*
na Ea ðuÞ
ð24Þ
a¼1
where Nrep is the number of representative atoms, na and Ea are the weight function and energy for representative atom a, respectively. The representative atom energy is then determined from its deformed environment. The simulation then proceeds with the same energy minimization as in the previous description. The force-based formulation takes on a different approach in which the need for an energy functional for the entire system is abandoned. Instead, it recognizes that the energy minimization process is equivalent to finding the zero-force position of each degree of freedom. Therefore, this form of the QC method adopts the approach of working directly from an expression for the forces on each representative atom. The force on a node can be defined as *
*
Fi ¼
oEa *
¼
oU j *
*
Nn X oEi ðuÞ ou *
i¼1
*
ou oU j
ð25Þ
*
where u and U j are the interpolated displacement field and the displacement of node j, respectively. The summation of the force expression can be reduced to a summation over a cluster of atoms around a representative atom to reduce the computational burden and still maintain a certain level of accuracy. Therefore, " # Nn X* X * * Fi ni gc Ni ðX c Þ ð26Þ c2Ci
where Ci refers to the set of atoms around representative atom i, gc is the atomic-level force experienced * by cluster atom c in the displacement field u; and ni is a weight function for a representative atom i. The disadvantages of both non-local QC formulations lie in the increase in computational cost when
123
compared with the local QC method. In determining the force or energy of each cluster of atoms, one is required to determine the updated position of each atom that is represented in its respective cluster. The benefit of these non-local formulations is the smooth transition region. However, the reduction of the degrees of freedom in regions of low-density representative atoms leads to some errors. The non-local formulations constitute some of the major developments of the QC method. However, the extension towards both dynamic and finite-temperatures remain a challenge. Efforts have also been made to extend the original QC method to couple with quantum mechanical calculations based on density functional theory in which the method can extend over three separate length scales and used to study the effect of chemical interactions on macroscopic material behavior (Lu et al. 2006a, b). This model was successfully applied to the study of the effects hydrogen impurities impose in the core of an edge dislocation in aluminum as well as on the long-range effects of the dislocation stress field. 3.2.6 The bridging scale method The bridging scale is a dynamic finite-temperature coupling method developed by Wagner and Liu (2003) that can be used for either two- or threedimensional simulations. This method has been applied to quasi-static nanotube bending, dynamic crack propagation, dynamic shear banding, and wave propagation (Park et al. 2005; Liu et al. 2006; Qian et al. 2004). Most recently, the bridging scale method has been used as the basis for developing a new multiscale procedure to couple a mesoscale discrete particle model and a macroscale continuum model of incompressible fluid flow. This method is called the mesoscopic bridging scale (MBS) method and was developed by Kojic et al. (2008) but the details of this method will be omitted because it is somewhat outside the scope of this review. The basic idea behind the bridging scale method is to combine the MD and FE systems together into one unified system. In so doing, some redundancy will exist between the overlapping systems. In order to remove the redundant MD degrees of freedom from the combined system, an impedance force must be added to the MD atoms at the boundary between MD and FE regions. These impedance forces are
Coupling atomistics and continuum in solids
introduced to compensate for the effects of the removed MD atoms. In this way, the FE region will exist everywhere, including those regions where MD is present, whereas the MD region will only exist in areas of interest or areas requiring high accuracy. The bridging scale method provides two major advantages when compared to other methods. The first is that the MD and FE equations of motion do not require the same time step because the continuum representation exists everywhere and the displacement field is decomposed into linearly independent coarse and fine scales. These coarse and fine scales correspond to the FE and MD displacement fields, respectively. This also implies that the FE mesh does not need to be refined and may maintain a coarse representation. The second major benefit stems from the addition of the impedance force. This force is responsible for allowing high-frequency waves to naturally dissipate out from the MD without propagating back and corrupting both the MD and FE simulations with oscillatory wave reflections. The schematic approach to the bridging scale method is illustrated in Fig. 7. At the top of the figure
Fig. 7 Bridging scale schematic
93
both MD and FE systems are shown as separate, distinct systems. When the two systems are combined into a unified system, the coupled multi-scale equations of motion can be written as, * € * MA u ¼ F
ð27Þ
* * € Md ¼ N T F
ð28Þ *
atomic mass matrix, u are the MD where MA is the * displacements, F is the total MD force vector, M is the finite element mass matrix, d are the FE nodal displacements, and N is a matrix containing the finite element shape functions. Equation 27 is the MD equation of motion and the displacements q can be determined by using any MD solver, while the MD forces can be determined from any valid interatomic potential energy function. In addition, one can use standard finite element methods to solve Eq. 28. It is important to note that the MD equation of motion is only solved in the MD region, which for the time being exists everywhere, but the FE equation of motion is solved across the entire system. The coupling of these equations is accomplished * T through the coarse scale internal force, N F ; which is * a function of the MD internal force, F : In the MD region the coarse scale internal force is determined by extrapolating the MD internal force to the nodes of the mesh by the use of FE shape functions. In this way, the MD internal forces can be thought of as defining the constitutive relation for the FE internal forces. We should also note that the FE equation of motion is redundant in the regions where both FE and MD exist. The FE equation of motion is only an approximation to the MD equation of motion. Therefore, if the FE nodal spacing is reduced to coincide with atom positions on a one-to-one correspondence, the coarse scale internal forces will equal the MD internal forces. Removing the redundancy in the coarse scale equation of motion will lead to the system depicted in the final stage of Fig. 7. This system is governed by the following coupled MD and FE equations of motion €* imp MA u ¼ f ðtÞ þ f0;m;n ðtÞ þ Rf0;m;n ðtÞ
ð29Þ
* € * Md ¼ N T F
ð30Þ
where
123
94
imp f0;m;n ðt Þ ¼
J. M. Wernik, S. A. Meguid
XZ m0 ;n0
t
hmm0 ;nn0 ðt sÞ
0
q0;m0 ;n0 ðsÞ u0;m0 ;n0 ðsÞ Rd0;m0 ;n0 ðsÞ ds ð31Þ where (l,m,n) and (l0 ,m0 ,n0 ) denote a system used to define unit cells. Equation 29 is the MD equation of motion. The first term on the right-hand side is the interatomic forces, which can be obtained from interatomic potential energy functions. The second term is the impedance force; it acts to replicate the unwanted MD degrees of freedom in their absence. This force only acts on the atoms along the interface between the MD and FE regions. This force also contains the time history kernel h(t - s), and acts to disperse the fine scale energy from the MD region into the surrounding FE region. The impedance force or equivalently, time history kernel, distinguishes this method from the others because it successfully allows high-frequency waves that cannot be captured by the coarse scale to be dissipated out of the MD region. The third term in the MD equation of motion is the random force term. It accounts for the thermal forces due to the removed fine scale degrees of freedom in their absence. For further details regarding this method and the derivation of the equations of motion the reader is referred to Wagner and Liu (2003). 3.2.7 Coupled atomistic dislocation dynamics method (CADD) The CADD method was developed by Shilkrot et al. (2004) and has been applied to the modeling of the two-dimensional nanoindentation process in single crystal films (Miller et al. 2004). The CADD method is claimed to offer two advantages over other methods; the ability to accommodate discrete dislocations in the continuum region and an algorithm for automatically detecting dislocations as they move from the atomistic region to the continuum region and then correctly ‘converting’ the atomistic dislocations into discrete dislocations or vice-versa. Other methods also attempt to follow deformations as they evolve, such as the QC method, but they do not have a means of transforming the dislocation into a continuum representation. Therefore, every dislocation in the model will require a full atomistic description which imposes a limit on
123
how many can be included in the simulation before the computational cost becomes that of a fully atomistic simulation. The CADD method is similar to all other MD/FE coupling methods in that it too begins from a separation of the spatial regions which are modeled either by a fully atomistic model or continuum finite elements. CADD extends these methods by addressing the above problem by allowing for the presence and movement of continuum discrete dislocations that interact with each other and with the atoms in the atomistic region via their elastic fields. The CADD method has been formulated specifically for 2D problems. The method has not been successfully extended to three dimensions although this is a current area of research. The CADD method was also primarily intended only for static zero-temperature simulations but extensions to both static and dynamic finite-temperature frameworks have been achieved for modeling of the nanoindentation process (Shiari et al. 2005; Shlarl et al. 2008). The CADD method couples to a discrete dislocation framework using similar coupling schemes presented above. In this way, this method can be thought of as an extension to the above coupling methods. The main distinguishing feature of the CADD method, the multiscale treatment of dislocations, is independent of the coupling scheme used, such that CADD can use any of the techniques previously developed. Therefore, the method of detecting the dislocations as they move across the atomistic/continuum interface is the primary contribution to literature. Just as in the non-local formulation of the QC method, the CADD method does not have a total energy functional defined for the entire system. Instead, this method enforces equilibrium through finding the zero-force position of each degree of freedom. Due to this methods ability to use previously developed coupling schemes, there is no reason to repeat the details so the reader is referred to Shilkrot et al. (2004) for a detailed illustration of the methods ability to detect and pass dislocations from the atomistic region to the continuum region and vice-versa. 3.2.8 The atomic-scale finite element method (AFEM) The atomic-scale finite element method was developed by Liu et al. (2004b) as an alternative to
Coupling atomistics and continuum in solids
traditional molecular dynamics. It is similar to the MPM method in which the traditional FE continuum representation was discarded for the material point approach. In this case, the atomistic representation is being replaced by the AFEM. This method is capable of modeling in both two-and three-dimensions, but so far is limited to static problems only. The AFEM is quoted as being capable of achieving larger and longer atomistic simulations with accuracy comparable to that of molecular dynamics. The AFEM has been widely used in the static study of carbon nanotube deformation processes (Liu et al. 2004b, 2005; Guo et al. 2007; Leung et al. 2006) and crack propagation (Adelzadeh et al. 2008), but shows promise for solving other physics related optimization problems. The increase in attainable system sizes and shorter simulation times associated with the AFEM method lies in the fact that it is an N-order method whereas typical MD simulations employing the conjugate gradient method are of order N2. This is because the AFEM uses both first and second order derivatives of the system energy in the minimization processes, whereas the conjugate gradient method only uses the first. By doing so, the AFEM is capable of reaching a stable state of energy minimization within one step (for linear systems), whereas typical conjugate gradient methods require multiple iterations. It is important to note that for nonlinear systems with complex interatomic potentials, the processes of employing the second derivative for energy minimization does not provide any major benefit, and the AFEM is no longer an N-order method. The reduction in computational cost shows great potential for increased system complexity and longer time scales but only for linear systems. An additional benefit is that the AFEM has exactly the same structure as the continuum finite element method. This means it works within the same theoretical framework of employing a system of linear equations (Ku = F), and therefore enables a seamless coupling between the two regions by establishing a unified system of governing equations. One of the major differences between the AFEM and FEM and reasons for increase in accuracy is that the former does not involve any traditional approximations of the conventional FE method, such as shape functions. The atomic-scale finite element method relies on developing material-specific atomistic elements
95
similar to the elements found in standard finite element methods. The elements are then used to represent the material under study. This presents possibly the largest limitation for the AFEM. For every material under study, a new atomistic element must be analytically developed to match the atomic structure and atomistic interactions depending on the interatomic potential used. This can be a very difficult and tedious process for classes of materials which have complicated long-range atomic interactions. An important characteristic, which must be implemented into the design of the AFEM elements, is the nonlocality of the discrete atoms it is representing, or equivalently, the multi-body interactions amongst atoms. It is also necessary to develop transitional elements between atomistic and continuum regions when coupling AFEM to FE. The transitional elements ensure smooth transition between the two types of elements and possess both local and nonlocal characteristics. The transitional elements are also responsible for removing the effects of ghost forces which plague most coupling methods. In this way, the transitional elements are a form of the transition region common in coupling methods and encompass the majority of the approximations. 3.2.9 Coarse grained molecular dynamics (CGMD) Refining the continuum region to atomic scales can be problematic and should be avoided. Continuum concepts rely on a continuous representation of the material and use phenomenological formulations and homogenization techniques to determine the field variables. Alternatively, the atomic scale consists of discrete particles and cannot be represented using a continuum approach. If the FE mesh is refined to atomic dimensions, the fundamental relations no longer hold true and the framework of continuum mechanics no longer applies. This problem will be discussed in detail in the coming section, but it is primarily for this reason, in addition to the reduced degrees of freedom, that a coarser representation of molecular dynamics called coarse-grained molecular dynamics (CGMD) (Rudd 2005; Rudd and Broughton 1998) has been introduced to the field of multiscale modeling. CGMD provides a means of attaining larger system sizes while avoiding the problem of scaling the continuum to fine scales. Traditionally, CGMD has been widely used in the field of
123
96
computational chemistry for the purpose of reducing the computational burden associated with the modeling of protein structures, lipid bilayers (Aoyagi et al. 2002), and polymer chains (Izvekov and Voth 2005), to name just a few. However, this coarser representation of MD also shows great potential for the field of atomistic and continuum coupling. Rudd and Broughton (2000, 2005) coupled CGMD to standard MD and avoided the standard FE altogether. This is because he was primarily interested in modeling MEMS devices and quantifying phonon spectra, where large continuum scales are generally avoidable. However, it is also possible to use CGMD as the finer scale in a coupled simulation to FE. This would in essence bring the atomistics region closer to the stable limits of continuum mechanics. CGMD has specifically been formulated to address several of the issues that limit the successfulness of other coupling methods. It is designed for dynamic and finite-temperature simulations thereby making it applicable to a larger variety of problems. It also avoids the effects of ghost forces by implementing a well-behaved physical response to stationary strain fields. CGMD avoids the issues related to high frequency waves propagating out of the MD region and reflecting back causing unphysical effects. This method can also be implemented with a wide variety of different MD models, including manybody interatomic potentials that extend beyond nearest neighbors. The general idea behind CGMD is to group the atomistic degrees of freedom into a cluster or bead of atoms. Each bead or cluster carries all the material properties associated with the atoms it is representing. In this case, the beads account for all the reduced degrees of freedom of the problem. Figure 8 demonstrates the general idea behind the CGMD method. The figure illustrates the coarse-graining of a low molecular weight polymer with epoxide groups at each end (Bockenheimer et al. 2002). At the top of the figure the chemical structure of the polymer is
Fig. 8 CGMD bead representation
123
J. M. Wernik, S. A. Meguid
shown. Directly beneath the structure is the corresponding coarse-grained representation whereby each bead embodies all the properties of the molecules it is representing. The coarse graining of molecular dynamics proceeds by obtaining a set of finite element equations from an underlying atomistic Hamiltonian. In this way the CGMD model is derived completely from MD principles and contains no continuum parameters. Although it employs a meshing scheme similar to that of FE, it avoids the continuum concepts needed in FE representations. As a result, this method provides a seamless coupling scheme between the two domains and controls the errors that normally arise in such coupling methods. A few key characteristics of the CGMD method merit discussion. First, the mesh can be refined, just as the other methods, to have the nodes correspond with atomic positions in regions where accuracy is of importance. When the entire mesh is refined in this manner, the coarse grained equations of motion reduce to exactly that of the MD equations. Secondly, the same timestep is used throughout the entire simulation. Although the coarse grained region can be run with a longer timestep, the computational efficiency of this method provides little incentive to do so. This is because one might run into problems with the use of a spatially varying timestep. Lastly, the way CGMD treats and manipulates the element stiffness matrix throughout the course of the simulation drastically reduces the computational burden. In this way, the simulation of billions of atoms can be achieved on a simple desktop workstation.
4 Emdl research activities in nanomechanics The Engineering Mechanics and Aerospace Design Laboratory (EMDL) at the University of Toronto is devoting attention to three main areas in nanoengineering/nanoscience. These are: (i) the mechanical behavior of single-walled (SWCNT) and multi-walled carbon nanotubes (MWCNT), (ii) continuum based modeling of the thermal effects on the mechanical behavior of carbon nanotubes, (iii) the use of nanofillers to reinforce bi-material interfaces. Further details regarding (i) and (ii) can be viewed in Wong et al. (2004), Liew et al. (2004), Meguid and Chen (2007). In this topical review, we provide a summary of our recent work on nano-reinforced interfaces.
Coupling atomistics and continuum in solids
The outstanding mechanical properties of carbon nanotubes make them ideal candidates for use in a wide variety of applications (Endo et al. 2004) including as reinforcement agents in composite materials. However, excellent fiber properties do not necessarily translate into the same properties for the bulk composite. Several issues pertaining to the alignment, dispersion, length, size, chirality, orientation and load transfer need to be optimized in order to achieve the best properties of the composite. Since experimentation at the nanoscale is still developing, the only way to quantify the effects of such parameters is predominantly through modeling. The large discrepancy in the scales involved in quantifying the macroscopic properties of nanocomposites from the addition of CNTs shows a clear need for integrated modeling methods. To date, there have been a number of multiscale models developed for the characterization of nano-reinforced polymer composites (NRPC). However, most of these efforts have adopted the atomistic-based continuum technique briefly described in the introductory section and outlined in references Li (2003), Odegard et al. (2002, 2003), Leung et al. (2005), Li and Chou (2006), Hu et al. (2005), Shi et al. (2005), Tserpes et al. (2008), Gao and Li (2005), Namilae and Chandra (2005), Li and Chou (2003), Zeng et al. (2008), Gates et al. (2005). The current state of literature indicates that comparatively limited efforts are being made in the development of concurrent multiscale models for these classes of problems. The current state of literature shows that a need exists for developing multiscale models for nanoreinforced interfaces. In this case, nano-fillers such as nanotubes and nanopowders are embedded in epoxy adhesives for improved bonding of either similar or dissimilar adherents. Our nano-reinforced interface research is concerned with the determination of the influence of nanofillers on the strength and toughness of thermoset epoxy adhesives used in the aerospace industry. This is due to the recent shift in aircraft design from the use of aluminum-based to compositebased airframes. Composite-based airframes exhibit excellent mechanical properties and lighter weight when compared to their aluminum counterparts. They are also capable of being constructed from adhesive joining processes as apposed to traditional mechanical fastening. Advantages of adhesive bonding over conventional mechanical fastening processes are
97
several. These include a better distribution of stresses over the entire bonded area, ability to optimize the geometry of the bonded joint, the ability to bond dissimilar materials, good fatigue resistance and damping, and reduced weight of the assembly. However, when compared with the adherend materials, an adhesive bonding layer generally offers a low level of mechanical strength and fracture toughness, especially at elevated temperatures and long term loading. The use of multi-layered laminates in the aircraft industry relies heavily on a strong and tough interface between the different layers. Unfortunately, the price of high strength epoxy is reduced toughness. This could be in part because of the complexity in the bonding process, the presence of porosities, micro-defects, and second-phase particles in the adhesive layer. The stringent toughness requirements of critical loading-bearing members necessitate the reinforcement of the adhesive layers and interfaces. The influence of nano-reinforcement in adhesively bonded joints on tensile strength, interfacial shear strength, and fracture toughness was investigated experimentally using SEM, TEM, and mechanical testing (Sun 2007), and will subsequently be studied theoretically using multiscale modeling. The experimental results reveal that the addition of nanofillers into the thermoset epoxy increases its strength and toughness (Meguid and Sun 2004). However, the increase in the amount of nanofillers beyond a certain critical value, could lead to deterioration in the strength and toughness of the resulting interface. Given below is a summary of the experimental work conducted at EMDL, the corresponding findings, and a description of the approach that will be adopted in the development of a concurrent multiscale model. 4.1 Substrates and nanofillers The new Airbus A380 is made of approximately 30% of glass reinforced aluminum composites (GLARE). This composite is composed of alternate glass fiber epoxies (GFRP) and aluminium layers as depicted in Fig. 9. A competitor of the A380, the Boeing 787, predominantly uses carbon fiber composites (CFRP) (Baker et al. 2004). In the research conducted by our group, Meguid and Sun (2004) used adherents made of these same materials; aluminum, CFRP, GFRP, and the epoxy was aerospace grade with low temperature
123
98
J. M. Wernik, S. A. Meguid
SEM images of the SWCNTs and TEM micrographs of the alumina nanofiber powder, respectively. A comparison of the nanofiller mechanical properties is also shown in the following table.
Property
Carbon nanotube
Alumina nanopowder
Size Young’s modulus
1–2 nm in diameter 1 TPa (in-plane)
2–4 nm in diameter 300 GPa
200 GPa (in-plane)
2000 MPa
Fig. 9 Glare pre-preg (After Baker et al. 2004)
Tensile strength
0.1–0.5 MPa (out-of-plane) Aspect ratio
10–30
20–80
Surface area
300–600 m2/g
300–700 m2/g
4.2 Experimental results
Fig. 10 SEM micrograph of carbon nanotubes (After Meguid and Sun 2004)
Fig. 11 TEM micrographs of alumina nanofiber powder (After Meguid and Sun 2004)
capabilities. The nano-reinforcing fillers consisted of single-walled carbon nanotubes (SWCNTs) and alumina nanofiber powder. These fillers were varied in weight percentage and there effects on the adhesive properties were studied. Figures 10 and 11 show the
123
Tensile and shear tests were conducted to study the effect of type and composition of nanofiller on the mechanical properties of the nano-reinforced interface. Figure 12a shows the tensile load against the displacement data for three different cases: epoxy adhesive with 2.5% alumina nanopowder (EANP), epoxy adhesive with 2.5% carbon nanotubes (EANT) and pure epoxy adhesive (EA). The figure indicates that the bonding strength increases as a result of the presence of uniformly dispersed carbon nanotubes and alumina nanoparticles. Figure 12b and c show improvement in the Young’s modulus as well as the ultimate tensile strength for the cases involving different weight fractions of homogeneously dispersed nanofillers used in the study. However, for weight percentages above 10%, the properties degrade to below the pure epoxy adhesive, indicative of the sensitivity of the epoxy to the strength of the concentration of the dispersed nanofillers. The shear test results are depicted in Fig. 13. Figure 13a shows load–displacement data resulting from the shear lap tests for the same three different cases. The filler concentration was 2.5% by weight for both EANP and EANT samples. The figure reveals that there is also an increase in the shear resistance as a result of dispersed carbon nanotubes and alumina nanoparticles. Figure 13b and c show the improvement in the shear modulus as well as the shear strength for the cases involving different concentrations of
Coupling atomistics and continuum in solids
Fig. 12 Tensile test results for nano-reinforced epoxy adhesives: (a) Load–displacement for EANP, EANT and EA, (b) Ultimate tensile strengths for both EANT and EANP specimens, (c) Young’s modulus for both EANT and EANP specimens (After Meguid and Sun 2004)
dispersed alumina nanopowders and carbon nanotubes. Analogues to the tensile tests, the results also reveal the sensitivity of the shear properties to the concentration of the nanofillers. An increase in the weight fraction of the nanofillers beyond 7–8% results in a reduction in the shear properties of the adhesive. Figure 14 shows three SEM micrographs of the failure surface of the carbon CFRP composite. They
99
Fig. 13 Shear test results for nano-reinforced epoxy adhesives: (a) Load–displacement for EANP, EANT and EA, (b) Ultimate shear strengths for both EANT and EANP specimens, (c) Shear modulus for both EANT and EANP specimens (After Meguid and Sun 2004)
provide useful information about the strength of the interface as a result of the dispersion of the alumina nanopowder. For example, Fig. 14a shows the surface failure at the interface from the CFRP side due to shear loading for a weight fraction of 2.5% dispersed alumina nanopowder. Two distinct regions exist: the first is the dark region showing exposed fibres of the
123
100
J. M. Wernik, S. A. Meguid
Fig. 14 SEM micrograph of shear test failure surfaces for nanopowder specimens: (a) 2.5% dispersed alumina nanopowder, (b) 7.5% dispersed alumina nanopowder, (c) 15% dispersed alumina nanopowder (After Meguid and Sun 2004)
CFRP substrate, while the second contains remnants of the nano-reinforced adhesive. The rough areas of the CFRP suggest a stronger bond exists at the interface as a result of the strengthening effect resulting from the dispersed nanofillers. Figure 14b and c, with their respective nanopowder weight concentrations of 7.5% and 15%, show reduced areas of the exposed fibres of CFRP. Figure 14b shows clearly fibre pull out at the interface of a section of the bare CFRP. However, larger areas are covered by the epoxy adhesive indicative of the reduction in the strength in comparison with Fig. 14a for the 2.5% case. For carbon nanotubes reinforced interfaces, Fig. 15a–c show typical failures of three different interfaces corresponding to nanotube weight concentrations of 2.5%, 7.5% and 15%, respectively. Again, it is observed that for the 2.5% case, large areas of carbon fibres of the substrate are exposed due to the strong adhesive bonding. We attribute that dependence of the tensile and shear properties of nano-reinforced interface to the
123
following. The nanosized fillers are characterized by large surface areas per unit gram. As the number of adhesively joined points increases, the cohesive strength of the epoxy increases leading to a higher mechanical strength of the interface. Our experimental results show that there is a limit to the number of dispersed nanofillers beyond which a drop in the properties is observed. Once the nanoparticles fully fill the gaps and porosities and all contact points are established, the additional particles cannot interact effectively within the epoxy adhesive and consequently, poor matrix infiltration occurs. It is also believed that agglomeration of the nanoparticles could act as failure initiation sites, which could result in lowering the strength and stiffness of the adhesive. 4.3 Development of a multiscale model In order to support the experimental work discussed above, the development of a multiscale model is required. In the case of NRPC, it has been standard
Coupling atomistics and continuum in solids
101
Fig. 15 SEM micrograph of failure surfaces for carbon nanotube specimens: (a) 2.5% dispersed carbon nanotubes, (b) 7.5% dispersed carbon nanotubes, (c) 15% dispersed carbon nanotubes (After Meguid and Sun 2004)
practice to develop representative volume elements (RVE) that can be used as the building block to study the macroscopic material response to different loading conditions. In this case, only the individual nanotube and surrounding polymer matrix are considered. Then, through the use of the atomistic-based continuum technique, the effective properties of the RVE are obtained. However, for the study of nanoreinforced interfaces, the problem necessitates a model capable of representing the system in its entirety. Here we must capture not only the interaction between nanotube and matrix but also the interaction between neighboring nanotubes, epoxy, and substrates as they may significantly affect the bulk properties. This will call for a reduction in the number of degrees of freedom of this complex system size. This can be achieved through concurrent coupling of CGMD and FE (Meguid and Sun 2005). In using a CGMD/FE coupled model to
Fig. 16 Multiscale model of nano-reinforced interface
represent the nano-reinforced interface system, we must correctly decompose the problem into CGMD, FE, and transition regions as depicted in Fig. 16. Here we have chosen to model the entire epoxy matrix and nanofiller subsystem in a CGMD representation, while the substrates are modeled in a FE framework.
123
102
The development of a concurrent CGMD/FE multiscale model presents a number of challenges. Most importantly, the interaction and load transfer between nanotube and epoxy matrix must correctly be accounted for in a CGMD system. Due to the similarity between this class of problems and those of NRPC, the same considerations pertaining to the load transfer between reinforcing agent and polymer matrix need to be included. Different approaches can be adopted in characterizing the mechanisms and magnitudes of load transfer between a nanotube and polymer matrix. The interfacial characteristics between the CNTs and polymer matrix remain unclear and researchers have reported a large range of interfacial shear stresses. Due to the lack of consistency, three approaches are possible. First, it can be assumed that no chemical bonding exists between the nanotubes and polymer matrix. In this case, only non-bonded potentials and van der Waals forces are considered. To avoid weak interfacial strength, some researchers propose that the chain of the polymer wrap around the nanotube in a helical fashion to enhance the non-bonded nanotube-polymer interaction, which has been observed experimentally (Lordi and Yao 2000). Figure 17 illustrates this process for a polystyrene polymer chain wrapping
Fig. 17 Polystyrene polymer chain wrapped around a SWCNT (a) side view (b) top view (After Hu et al. 2005)
123
J. M. Wernik, S. A. Meguid
around a SWCNT. The second approach is to assume that there exist strong chemical bonds. In this case, C–C covalent bonds are included between the nanotube and polymer, which increases the interfacial strength significantly. The third consideration assumes that covalent cross-links form between the nanotube and polymer matrix. In this case, only a small percentage of covalent bonds form from the introduction of multifunctional amines which act as intermediary bonding sites between the nanotube and polymer chains. However, it is possible that the chemical bonding in the form of functionalization may compromise the properties of the nanotube by introducing structural changes in the graphitic layers of the nanotube (Fiedler et al. 2006). Frankland et al. (2002) studied the effect of functionalization, in the form of chemical cross-links, on the mechanical properties of a polymer composite. They found that by introducing chemical cross-links, involving as little as 1% of the carbon atoms in the nanotube, the properties of the nanotube-polymer interface can be increased significantly. A depiction of the model used in one of their studies is shown in Fig. 18. When developing the concurrent multiscale model it will be important to clearly identify which of the above processes will be considered.
Coupling atomistics and continuum in solids
103
Fig. 18 Functionalization of CNT through chemical cross-linking: (a) Nonfuctionalized polyethylene nanotube RVE, (b) crosslink between nanotube and polymer chain, (c) arrangement of cross-links, (d) functionalized polyethylene nanotube RVE (After Frankland et al. 2002)
Next, we must also accurately represent the individual CNTs and epoxy polymer chains without having to explicitly model every atom. It will be difficult to establish the optimal magnitude of coarse graining in the molecular region. This is because if the nanotubes and polymer chains are coarse grained to the extent that each individual chain or tube is represented as a single bead, the interaction mechanisms mentioned above will be extremely difficult to implement and the simulation may loose a certain level of accuracy. On the other hand, if the level of coarse graining is too low the number of degrees of freedom will be excessively large and the simulation may take an extremely long time to converge. As with many other factors, there is a delicate balance between computational cost and level of accuracy. Of course, as with all of the other concurrent coupling methods presented in the review, the exact details of the coupling scheme will have to be decided. Since no efforts have been made in coupling CGMD to FE, serious attention will have to been given to this area of the problem. Furthermore, the coupling scheme will have to avoid the limitations presented in the following section, which alone can be a great challenge. These are just some of the issues that need to be addressed in order to accurately model the nano-reinforced interface problem. Further details regarding the concurrent CGMD/FE multiscale model will be the subject of our forthcoming publication.
5 Limitations of current coupling methods In spite of the excellent contributions made by the referenced scientific community, these contributions
suffer from a number of challenges. These are summarized below. 5.1 Finite element reduction Many of the methods discussed above adopt a similar approach of refining the FE mesh down to the size of the atomic lattice dimensions. The purpose of this refinement is to ensure seamless transition from the FE to the MD regions. This approach makes for a much easier computation in the transition region between the two domains by allowing the interface atoms to occupy the same location as a node thereby defining the nodal positions by the MD atom positions. As mentioned earlier, the fundamental assumption that continuum mechanics relies on is the homogeneous and continuous description of the material in question. Essentially, continuum theory requires that variables such as temperature, density, pressure, and stress can be defined by an averaging process and that these variables are assumed to be smoothly varying continuous functions of position. As the length scales of interest become equal to microns and nanometers, quantities such as stress, strain, density and temperature begin to lose their traditional meaning. Since molecular dynamic simulations occur on the nanometer scale or less, the material must be described using discrete points such as atoms or molecules in order to accurately capture the physics that occur at this scale. Therefore, by scaling the continuum region down to the MD scale the fundamental continuum assumptions become invalidated and no longer hold true. It is difficult to establish at exactly what point of mesh refinement the continuum region loses its validity but in order to do so we must keep the
123
104
underlying assumptions in mind. The continuum formulations remain intact provided that the length scales of interest are much larger than the length scales of molecular variation in the system under consideration. This is so that the concepts of stress and strain remain valid and applicable. When modeling systems of molecular scale, we must ask ourselves if the molecular details of the system can be ignored. The methods presented above have adopted the approach of fine-scaling the continuum region when the logical alternative would have been to coarsen the atomistic region. 5.2 Wave reflection The problem in this case stems from the fact that the wavelength emitted by the MD region is much smaller than what can be captured by the FE region. Since energy conservation is enforced, the wave is reflected back into the MD region leading to oscillatory solutions. When reducing the FE mesh to atomic dimensions, it is assumed that the FE mesh would be able to capture these high-frequency waves and allow them to pass through unimpeded, although other effects such as differences in compliance between the two regions will still cause some wave reflection. The most successful way of permitting high-frequency wave travel through the interface is through the incorporation of an impedance force acting on the interface atoms. This was done in the bridging scale method which showed that if the impedance force was not incorporated only 30% of the MD energy was transferred into the continuum whereas if the force was incorporated nearly 99.9% of the energy was transferred (Liu et al. 2004a, b). Further details can be viewed in references Karpov et al. (2005), Wagner et al. (2004). 5.3 Ghost forces The respective local and non-local mismatch between the FE and MD regions introduces spurious effects across the interface of the transition region. The molecular dynamics non-local formulation implies that the energies of individual atoms interact not only with their immediate neighbours but also with others a distance rcut-off away. The local formulation of the continuum implies that energy is uniformly defined
123
J. M. Wernik, S. A. Meguid
throughout an element and is only dependant on the energy defined at the nodes of that element. Obtaining the nodal and atomic forces in the system, through minimizing the total potential energy, leads to non-accurate forces within the transition region. This can be attributed to the introduction of the overlap atoms as a means of maintaining the proper neighboring for the non-local MD atoms. The overlap atom energies are omitted from the total energy functional because they introduce redundant degrees of freedom. The basic idea behind these ‘ghostforces’ is that some atoms in the MD region have missing force contributions because the overlap atom energies are not included in the calculation. Furthermore, some FE nodes have extra force contributions acting on them. Therefore, some atoms are not being treated as true atoms and some nodes are not being treated as true nodes. Since the interfacial forces are in error, the atom and nodal displacements at the interface are also incorrect. An example of how ghost forces arise can be found in the review by Curtin and Miller (2003). A means of correcting ghost forces in the QC method has been developed by the Rodney and coworkers (Shenoy et al. 1999). He also extended the method to three dimensions and used it to study junction formation and destruction between interacting dislocations (Rodney and Phillips 1999). In an attempt to overcome this difficulty, they introduced the missing forces as dead loads. Ghost forces are calculated from the initial reference state of the material system and are then subtracted from the relevant degrees of freedom throughout the entire course of the simulation. With these forces acting as dead loads, the QC energy functional is modified to include the work done by these forces during the deformation. It should be noted that this method of correcting for ghost forces could be applied to other methods as well. An alternative approach to eliminating ghost forces is to abandon the requirement for a welldefined energy functional for the system and to drive the system to equilibrium by seeking a configuration for which the force on all the atoms is zero. In this way, the continuum region is made non-local and the locality mismatch is avoided. In using this approach, the forces do not need to be determined from the derivatives of the energy functional. Instead, expressions for the forces are approximated. This approach
Coupling atomistics and continuum in solids
is adopted in the force-based formulation of the nonlocal QC method as well as in the CADD method. 5.4 Uniform deformation gradient When employing the Cauchy-Born rule in certain coupling methods, the continuum deformation is restricted to being uniform. This implies that the underlying atomistic deformation is also uniform and that each atom bond exhibits the same deformation as its neighbor throughout the entire element. In regions employing a coarse FE mesh, atomistic details will not be captured due to the use of this approximation. In engineering problems, this is rarely the case and this assumption imposes unrealistic conditions on the validity of these methods. For example, plasticity and contact problems can cause sharp gradients in the deformation field and this assumption would not maintain its validity for these types of problems. In order to capture the atomistic details the FE mesh would have to be refined to the point where the number of degrees of freedom would be excessively large. The uniform deformation gradient also implies that the only type of element that can be used throughout the FE region is restricted to three-noded triangular elements. These elements are known for their computational efficiency but lack accuracy. 5.5 Unrealistic interatomic potentials A simulation is realistic if it can replicate the behavior of a real system. For molecular dynamic simulations this can be achieved only to the extent that the interatomic forces are similar to those that real atoms would experience in the same configuration. The forces are determined as gradients of the interatomic potential energy function. The accuracy of the simulation therefore depends on the ability of the potential to replicate the behavior of the material. In order for these simulations to be viable, it is important to incorporate interatomic potentials that accurately represent the interaction between atoms for a particular material. These can depend on; attractive and repulsive forces, bonding angles, bond types, electron densities, vacancies, cohesive energies, and bond types, just to mention a few. The complex nature of these interactions leads to difficulties in developing accurate and valid expressions for these potentials. Furthermore, it is very difficult to
105
justify or validate the accuracy of these potentials. It is also important to note that interatomic potentials have incorporated a cutoff distance in their formulations to disregard the interactions between far neighbouring atoms. This results in simpler calculations but introduces a separate problem as well. When a potential is truncated or ‘cut short’, a particle pair that passes the cutoff distance will introduce a jump or step function into the energy. Many of these events can corrupt the energy conservation of an MD simulation. It is for the above reasons that the scientific community has adopted the use of the most cost effective expressions, which can be drastically over simplified. There is a direct trade-off between accuracy and cost when determining the choice of potential for the system. 5.6 Solution algorithms and timestep The timestep in FE simulation is governed by the smallest element in the mesh. Therefore, if the FE mesh is refined down to the atomic scale, the simulation will evolve very slowly in the sense that many timesteps will be used to simulate the dynamics in the areas of little interest. It also seems inefficient that the FE simulation should evolve at the same time scales as the atomistic level. This problem will directly affect the computational cost of the simulation. The bridging scale method was developed as a means of avoiding this issue by completely separating the two length scales whereby the FE region exists everywhere and the MD region only exists in areas of interest. This way, there is no restriction on the FE mesh to be refined to atomic dimensions, thereby the continuum simulation can evolve on a larger timestep. Time integration algorithms are based on finite difference methods, where the time is discretized on a finite grid. The timestep is the distance between successive points on the grid. These algorithms are only approximate and have two errors associated with them. First, truncation errors arise because these algorithms involve Taylor expansions which need to be truncated at some point. Secondly, round-off errors arise from the implementation of the algorithms on operating systems which use a finite number of digits in their arithmetic operations. Both these errors can lead to a divergence of the solution. As very small timesteps need to be used in the atomistic domain, the
123
106
highly iterative nature of these multiscale simulations can cause error amplification. 5.7 Degrees of freedom/relevance to design One of the main goals of multiscale modeling is to drastically reduce the number of degrees of freedom of engineering problems while maintaining accuracy in regions of interest so that the modeling of largescale problems could be realized. The methods developed thus far have not achieved the goal of attaining large realistic system sizes. Instead, these methods have focused on very limited problems and only managed to minimally reduce the computational cost of these problems. Even with the use of parallel computing these methods still require a long time to simulate the problems of interest. All single scale modeling methods provide results that lead to an understanding of how a material or component will behave under specific conditions. This information is then passed on to designers to learn from and incorporate into future designs. The multiscale methods developed thus far have not demonstrated the usefulness for engineering design. They do not encompass the applicability or relevance to a wide variety of engineering problems. They are problem specific and only demonstrate the minor ability to reduce the degrees of freedom of the problem rather than produce results that are meaningful to engineers. Unless the efforts are directed to real engineering problems, the efforts become mere numerical gymnastics. Clearly, the results can be used to assist in developing new material systems. However, this is also limited by the prohibitive cost associated with nanotechnology.
6 Conclusion and future work The contents of this review have demonstrated the significant efforts made by the research community in developing a seamless coupling scheme between atomistic and continuum scales. The research being conducted in the field of atomistic/continuum coupling is bringing us ever closer to a truly seamless coupling scheme capable of modeling a vast variety of problems. We are now able to model problems of ever increasing size, complexity and accuracy. However, despite these efforts, there still exist major challenges
123
J. M. Wernik, S. A. Meguid
and limitations in achieving such a goal. The methods presented above have their own advantages, but each lacks where another prevails. An ideal coupling method would incorporate all the positive aspects of each method into one sole approach. This method would be capable of both static and dynamic simulations at finite-temperatures. It would substantially reduce the degrees of freedom of any problem thereby reducing the computational burden and extending its abilities to large-scale problems. Other features of the method would include the ability to model in both two- and three-dimensional space, be capable of longer timescales, and applicable to a larger set of engineering problems. It would also avoid the issue of disparate timescales, the major bottleneck in most simulations. Few methods have achieved some of these goals, while others have developed ways of accelerating the simulations for very specific cases. The development of this type of method seems more realistic than the advancement of computers with processing abilities capable of modeling entire macroscale problems using purely atomistic approaches. For this to be realized, it is crucial that the scientific community acknowledges the advancements being made in computational material science and combine efforts to achieve this goal. The field of nanoengineering shows great potential for application of multiscale modeling methods. This is because of the intimate dependence between the atomistic and continuum scales in such systems as well as it is very difficult to accurately model nano systems using individual methods either because of the computational cost or complexity of the system. With the development of a truly seamless coupling method realistic system sizes can be achieved and used to validate designs and experimental processes. The study of nano-reinforced interfaces for one can use multiscale models to optimize filler type and size for improved toughness in epoxy adhesives used in the aerospace industry. Without the use of these models, one might be at liberty to make unjust recommendations which could ultimately lead to flawed or unsafe designs. The field of multiscale modeling has great potential and is ultimately reliant on the development of a seamless coupling scheme between the atomistic and continuum scales in a computationally efficient manner. This review has addressed the main atomistic/ continuum coupling methods used in computational
Coupling atomistics and continuum in solids
nanomechanics. There are other methods which have been omitted from discussion either because the focus was outside the scales interest or the general approach is very similar to the methods discussed above. The key limitations and challenges of the above methods have also been addressed for the reader to take into consideration. The review addressed also the recent developments in the field of nano-reinforced composites. This included a discussion on the research conducted by the authors in nano-reinforced interfaces. The issues of concern in nanomechanics are the subjects of our ongoing research at EMDL. We will continue to focus our efforts towards developing a perfectly seamless transition region that can overcome many of the difficulties that currently exist in the field of multiscale modeling as mentioned above. Specifically, efforts will be applied to studying the effects of nanofillers on the bond strength and toughness of interfaces. We will attempt to further this research through the development of multiscale models that can accurately simulate the nano-interface problem. These models are currently being developed. Details of our efforts will be the subject of our forthcoming publication. Acknowledgements The authors wish to acknowledge the financial support provided by Natural Sciences and Engineering Research Council (NSERC) of Canada.
References Abraham, F.F.: Madly spanning the length scales in dynamic fracture. Comp. Model. Eng. Sci. 1(4), 63–69 (2000) Abraham, F.F., Broughton, J.Q.: Spanning the continuum to quantum length scales in a dynamic simulation of brittle fracture. Europhys. Lett. 44(6), 783–787 (1998). doi: 10.1209/epl/i1998-00536-9 Abraham, F.F., Walkup, R., Gao, H., Duchaineau, M., De la Rubia, T.D., Seager, M.: Simulating materials failure by using up to one billion atoms and the world’s fastest computer: brittle fracture. Proc. Natl. Acad. Sci. USA 99(9), 5777–5782 (2002). doi:10.1073/pnas.062012699 Adelzadeh, M., Shodja, H.M., Rafii-Tabar, H.: Computational modeling of the interaction of two edge cracks, and two edge cracks interacting with a nanovoid, via an atomistic finite element method. Comput. Mater. Sci. 42, 186–193 (2008). doi:10.1016/j.commatsci.2007.07.012 Allen, M.P., Tildesley, D.J.: Computer Simulation of Liquids. Claredon Press, Oxford (1987) Almeida, T.S., Coutinho, K., Cabral, B.J.C., Canuto, S.: Electronic properties of liquid ammonia: a sequential molecular
107 dynamics/quantum mechanics approach. J. Chem. Phys. 128, 014506 (2008) An, W., Wu, X., Yang, J.L., Zeng, X.C.: Adsorption and surface reactivity on single-walled boron nitride nanotubes containing stone-wales defects. J. Phys. Chem. C 111, 14105–14112 (2007). doi:10.1021/jp072443w Aoyagi, T., et al.: A general-purpose coarse-grained molecular dynamics program. Comput. Phys. Commun. 145, 267– 279 (2002). doi:10.1016/S0010-4655(02)00271-0 Baker, A., Dutton, S., Kelly, D.: Composite Materials for Aircraft Structures, 2nd edn. American Institute of Aeronautics and Astronautics, USA (2004) Bockenheimer, C., Valeske, B., Possart, W.: Network structure in epoxy aluminium bonds after mechanical treatment. Int. J. Adhes. Adhes. 22, 349–356 (2002). doi:10.1016/ S0143-7496(02)00014-3 Broughton, J., Abraham, F.: Concurrent coupling of length scales methodology and application. Phys. Rev. B 60, 4 (1999). doi:10.1103/PhysRevB.60.2391 Buehler, M., Kong, Y., Gao, H.: Deformation mechanisms of very long single-wall carbon nanotubes subject to compressive loading. J. Eng. Mater. Technol. 126, 245–249 (2004) Cao, G., Chen, X.: Buckling behavior of single-walled carbon nanotubes and a targeted molecular mechanics approach. Phys. Rev. B, 74, 10pp (2006a) Cao, G., Chen, X.: Buckling of single-walled carbon nanotubes upon bending: molecular dynamics simulations and the finite element method. Phys. Rev. B. 73, 11pp (2006b) Cao, G., Chen, X.: The effects of chirality and boundary conditions on the mechanical properties of single-walled carbon nanotubes. Int. J. Solids Struct. 44, 5447–5465 (2007). doi:10.1016/j.ijsolstr.2007.01.005 Carlsson, A.: Beyond pair potentials in elemental transition metals and semiconductors. Solid. State. Phys. 43, 1–91 (1990). doi:10.1016/S0081-1947(08)60323-9 Cavallotti, C., Di Stanislao, M., Moscatelli, D., Veneroni, A.: Materials computation towards technological impact: the multiscale approach to thin films deposition. Electrochim. Acta 50, 4566–4575 (2005). doi:10.1016/j.electacta.2004. 10.092 Chandraseker, K.I., Mukherjee, S.: Atomistic-continuum and ab initio estimation of the elastic moduli of single-walled carbon nanotubes. Comput. Mater. Sci. 40, 147–158 (2007). doi:10.1016/j.commatsci.2006.11.014 Clementi, E.: Global scientific and engineering simulations on scalar, vector and parallel LCAP-type supercomputers. Philos. Trans. R. Soc. Lond. A 326, 445–470 (1988). doi:10.1098/rsta.1988.0097 Curtin, W., Miller, R.: Atomistic/continuum coupling in computational materials science. Model Simul. Mater. Sci. Eng. 11, R33–R68 (2003). doi:10.1088/0965-0393/ 11/3/201 Daw, M., Baskes, M.: Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals. Phys. Rev. B. 29, 12 (1984). doi:10.1103/Phys RevB.29.6443 Ding, F.: Theoretical study of the stability of defects in singlewalled carbon nanotubes as a function of their distance from the nanotube end. Phys. Rev. B. 72, 245–409 (2005)
123
108 Dupuy L., Tadmor, E., Miller, R., Phillips, R.: Finite-temperature quasicontinuum: molecular dynamics without all the atoms. PRL 95 (2005) Endo, M., Hayashi, T., Kim, Y.A., Terrones, M., Dresselhaus, M.S.: Applications of carbon nanotubes in the twenty-first century. Philos. Trans. R. Soc. Lond. A 362, 2223–2238 (2004). doi:10.1098/rsta.2004.1437 Ercolessi, F.: A Molecular Dynamics Primer. Spring College in Computational Physics, Trieste (1997) Esfarjani, K., Gorjizadeh, N., Nasrollahi, Z.: Molecular dynamics of single wall carbon nanotube growth on nickel surface. Comput. Mater. Sci. 36, 117–120 (2006). doi: 10.1016/j.commatsci.2005.02.022 Fiedler, B., Gojny, F.H., Wichmann, M.H.G., Nolte, M.C.M., Schulte, K.: Fundamental aspects of nano-reinforced composites. Compos. Sci. Technol. 66, 3115–3125 (2006). doi:10.1016/j.compscitech.2005.01.014 Fish, J., Nuggehally, M.A., Shepard, M.S., Picu, C.R., Badia, S., Parks, M.L., et al.: Concurrent AtC coupling based on a blend of the continuum stress and the atomistic force. Comput. Methods Appl. Mech. Eng. 196, 4548–4560 (2007). doi:10.1016/j.cma.2007.05.020 Frankland, S.J.V., Caglar, A., Brenner, D.W., Griebel, M.: Molecular simulation of the influence of chemical crosslinks on the shear strength of carbon nanotube-polymer interfaces. J. Phys. Chem. B. 106, 12 (2002). doi:10.1021/ jp015591? Galano, A., Francisco-Marquez, M.: Reactivity of silicon and germanium doped CNTs toward aromatic sulfur compounds: a theoretical approach. Chem. Phys. 345, 87–94 (2008). doi:10.1016/j.chemphys.2008.01.040 Gao, X.L., Li, K.: A shear-lag model for carbon nanotubereinforced polymer composites. Int. J. Solids Struct. 42, 1649–1667 (2005). doi:10.1016/j.ijsolstr.2004.08.020 Gates, T.S., Odegard, G.M., Frankland, S.J.V., Clancy, T.C.: Computational materials: multi-scale modeling and simulation of nanostructured materials. Compos. Sci. Technol. 65, 2416–2434 (2005). doi:10.1016/j.compsci tech.2005.06.009 Ghoniem, N., Busso, E., Kioussis, N., Huang, H.: Multiscale modeling of nanomechanics and micromechanics an overview. Philos. Mag. V83(31–34), 3475–3528 (2003). doi:10.1080/14786430310001607388 Gong, S.X., Meguid, S.A.: On the elastic fields of an elliptical inhomogeneity under plane deformation. Proc. R. Soc. Lond. A. 443(1919), 457–471 (1993) Gong, N., Liang, Y., Yao, Y., Liu, B.: Static and dynamic analysis of carbon nanotube cantilever based on molecular dynamics simulation. Key Eng. Mater. 375–376, 631–635 (2008) Gumbsch, P.: An atomistic study of brittle fracture: toward explicit failure criteria from atomistic modeling. J. Mater. Res. 10(11), 2897–2907 (1995). doi:10.1557/JMR.1995. 2897 Gumbsch, P., Beltz, G.: On the continuum versus atomistic descriptions of dislocation nucleation and cleavage in nickel. Model. Simul. Mater. Sci. Eng. 3, 597–613 (1995). doi:10.1088/0965-0393/3/5/002 Guo, Z., Yang, W.: MPM/MD handshaking method for multiscale simulation and its application to high energy
123
J. M. Wernik, S. A. Meguid cluster impacts. Int. J. Mech. Sci. 48, 145–159 (2006). doi:10.1016/j.ijmecsci.2005.08.007 Guo, X., Leung, A.Y.T., Jiang, H., He, X.Q., Huang, Y.: Critical strain of carbon nanotubes: an atomic-scale finite element study. J. Appl. Mech. 74, 347–351 (2007) Hao, S., Moran, B., Liu, W., Olson, G.: A hierarchical multiphysics model for design of high toughness steels. J. Comput. Aided Mater. Des. 10, 99–142 (2003). doi:10.1023/B:JCAD.0000036813.66891.41 Haslam, A.J., Moldovan, D., Phillpot, S.R., Wolf, D., Gleiter, H.: Combined atomistic and mesoscale simulation of grain growth in nanocrystalline thin films. Comput. Mater. Sci. 23, 15–32 (2002). doi:10.1016/S0927-0256(01)002 18-X Hockney, R.W.: The potential calculation and some applications. Methods Comput. Phys. 9, 135–211 (1970) Hu, N., Fukunaga, H., Lu, C., Kameyama, M., Yan, B.: Prediction of elastic properties of carbon nanotube reinforced composites. Proc. R. Soc. 461, 1685–1710 (2005). doi: 10.1098/rspa.2004.1422 Hughes, T.R.: The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Englewood Cliffs, NJ Prentice-Hall (1987) Izvekov, S., Voth, G.: A multiscale coarse-graining method for biomolecular systems. Phys. Chem. B. 109, 2469–2473 (2005) Jones, J.: On the determination of molecular fields: II from the equation of state of a gas. Proc. R. Soc. A106, 463–477 (1924). doi:10.1098/rspa.1924.0082 Karpov, E., Wagner, G., Liu, W.: A green’s function approach to deriving non-reflecting boundary conditions in molecular dynamics simulations. Int. J. Numer. Methods Eng. 62, 1250–1262 (2005). doi:10.1002/nme.1234 Knap, J., Ortiz, M.: An analysis of the quasicontinuum method. J. Mech. Phys. Solids 49, 1899–1923 (2001). doi:10.1016/ S0022-5096(01)00034-5 Kohlhoff, S., Gumbsch, P., Fischmeister, H.: Crack propagation in b.c.c crystals studied with a combined finiteelement and atomistic model. Philos. Mag. A 64(4), 851– 878 (1991). doi:10.1080/01418619108213953 Kojic, M., Filipovic, N., Tsuda, A.: A mesoscopic bridging scale method for fluids and coupling dissipative particle dynamics with continuum finite element method. Comput. Methods Appl. Eng. 197, 821–833 (2008). doi:10.1016/ j.cma.2007.09.011 Leung, A.Y.T., Guo, X., He, X.Q., Kitipornchai, S.: A continuum model for zigzag single-walled carbon nanotubes. Appl. Phys. Lett. 86, 083110 (2005) Leung, A.Y.T., Guo, X., He, X.Q., Jiang, H., Hunag, Y.: Postbuckling of carbon nanotubes by atomic-scale finite element. J. Appl. Phys. 99, 124308 (2006) Li, C., Chou, T.: Elastic moduli of multi-walled carbon nanotubes and the effect of van der waals forces. Compos. Sci. Technol. 63, 1517–1524 (2003a). doi:10.1016/S02663538(03)00072-1 Li, C., Chou, T.: A structural mechanics approach for the analysis of carbon nanotubes. Int. J. Solids Struct. 40, 2487–2499 (2003b). doi:10.1016/S0020-7683(03)00056-8 Li, C., Chou, T.: Multiscale modeling of compressive behavior of carbon nanotube/polymer composites. Compos. Sci.
Coupling atomistics and continuum in solids Technol. 66, 2409–2414 (2006). doi:10.1016/j.compsci tech.2006.01.013 Liew, K.M., Wong, C.H., He, X.Q., Tan, M.J., Meguid, S.A.: Nanomechanics of single and multiwalled carbon nanotubes. Phys. Rev. B. 69, 11 (2004). doi:10.1103/Phys RevB.69.115429 Liu, W., Karpov, E., Zhang, S., Park, H.: An introduction to computational nanomechanics and materials. Comput. Methods Appl. Mech. Eng. 193, 1529–1578 (2004a). doi: 10.1016/j.cma.2003.12.008 Liu, B., Huang, Y., Jiang, H., Qu, S., Hwang, K.C.: The atomic-scale finite element method. Comput. Methods Appl. Mech. Eng. 193, 1849–1864 (2004b). doi:10.1016/ j.cma.2003.12.037 Liu, B., Jiang, H., Huang, Y., Qu, S., Yu, M.-F., Hwang, K.C.: Atomic-scale finite element method in multiscale computation with applications to carbon nanotubes. Phys. Rev. B. 72, 035435 (2005) Liu, W.K., Park, H.S., Qian, D., Karpov, E.G., Kadowaki, H., Wagner, G.J.: Bridging scale methods for nanomechanics and materials. Comput. Methods Appl. Mech. Eng. 195, 1407–1421 (2006). doi:10.1016/j.cma.2005.05.042 Lordi, V., Yao, N.: Molecular mechanics of binding in carbonnanotube-polymer composites. J. Mater. Res. 15, 12 (2000) Lu, G., Tadmor, E.B., Kaxiras, E.: From electrons to finite elements: a concurrent multiscale approach for metals. Phys. Rev. B. 73, 124108 (2006a) Lu, H., Daphalapurkar, N., Wang, B., Roy, S., Komanduri, R.: Multiscale simulation from atomistic to continuum–coupling molecular dynamics (MD) with the material point method (MPM). Philos. Mag 86(20), 2971–2994 (2006b). doi:10.1080/14786430600625578 Meguid, S.A., Chen, B.J.: Modelling temperature-dependent fracture nucleation of SWCNTs using atomistic-based continuum theory. Int. J. Sol. Struct. 44(11–12), 3828– 3839 (2007) Meguid, S.A., Sun, Y.: On the tensile and shear strength of nano-reinforced composite interfaces. Mater. Des. 325, 289–296 (2004) Meguid, S.A., Sun, Y.: Intelligent condition monitoring of aerospace composites: part I—nano reinforced surfaces & interfaces. Int. J. Mech. Mater. Des. 37–52 (2005) Meguid, S.A., Wang, X.D.: On the dynamic interaction between a microdefect and a main crack. Proc. R. Soc. Lond. A. 448(1934), 449–464 (1995) Meguid, S.A., Wang, X.D.: The dynamic interaction of a microcrack with a main crack under antiplane loading. Int. J. Solids Struct. 31(8), 1085–1097 (1994). doi:10.1016/ 0020-7683(94)90165-1 Micheal Lai, W., Rubin, D., Krempl, E.: Introduction to Continuum Mechanics, Revised Edition. Pergamon Press, NY (1978) Miller, R., Tadmor, E.: The quasicontinuum method: overview, applications and current directions. J. Comput. Aided. Mater. Des. 9, 203–239 (2002). doi:10.1023/A:102609801 0127 Miller, R., Tadmor, E., Phillips, R., Ortiz, M.: Quasicontinuum simulation of fracture at the atomic scale. Model. Simul. Mater. Sci. Eng. 6, 607–638 (1998a). doi:10.1088/09650393/6/5/008
109 Miller, R., Ortiz, M., Phillips, R., Shenoy, V., Tadmor, E.: Quasicontinuum models of fracture and plasticity. Eng. Fract. Mech. 61, 427–444 (1998b). doi:10.1016/S00137944(98)00047-2 Miller, R.E., Shilkrot, L.E., Curtin, W.A.: A coupled atomistics and discrete dislocation plasticity simulation of nanoindentation into single crystal films. Acta Mater. 52, 271– 284 (2004). doi:10.1016/j.actamat.2003.09.011 Mylvaganam, K., Vodenitcharova, T., Zhang, L.C.: The bending-kinking analysis of a single-walled carbon nanotube–a combined molecular dynamics and continuum mechanics technique. J. Mater. Sci. 41, 3341–3347 (2006). doi:10.1007/s10853-005-5389-7 Namilae, S., Chandra, N.: Multiscale model to study the effect of interfaces in carbon nanotube-based composites. J. Eng. Mater. Technol. 127 222–232 (2005) Odegard, G.M., Frankland, S.J.V.: Effect of nanotube functionalisation on the elastic properties of polyethylene nanotube composites. AIAA J. 43(8), 1828–1835 (2005). doi:10.2514/1.9468 Odegard, G.M., Gates, T.S., Nicholson, L.M., Wise, K.E.: Equivalent-continuum modeling of nano-structured materials. Compos. Sci. Technol. 62, 1869–1880 (2002). doi: 10.1016/S0266-3538(02)00113-6 Odegard, M., Gates, T.S., Wise, K.E., Park, C., Siochi, E.J.: Constitutive modeling of nanotube-reinforced polymer composites. Compos. Sci. Technol. 63, 1671–1687 (2003). doi:10.1016/S0266-3538(03)00063-0 Pantano, A., Parks, D.M., Boyce, M.C.: Mechanics of deformation of single- and multi-wall carbon nanotubes. J. Mech. Phys Solids 52, 789–821 (2004). doi:10.1016/ j.jmps.2003.08.004 Park, H., Liu, W.: An introduction and tutorial on multiplescale analysis in solids. Comput. Methods Appl. Mech. Eng. 193, 1733–1772 (2004). doi:10.1016/j.cma.2003. 12.054 Park, H.S., Karpov, E.G., Liu, W.K., Klein, P.A.: The bridging scale for two-dimensional atomistic/continuum coupling. Philos. Mag. 85(1), 79–113 (2005). doi:10.1080/14786430 412331300163 Qian, D., Wagner, G., Liu, W.: A multiscale projection method for the analysis of carbon nanotubes. Comput. Methods Appl. Mech. Eng. 193(17–20), 1603–1632 (2004). doi: 10.1016/j.cma.2003.12.016 Qu, S., Shastry, V., Curtin, W., Miller, R.: A finite-temperature dynamic coupled atomistic/discrete dislocation method. Model. Simul. Mater. Sci. Eng. 13, 1101–1118 (2005). doi:10.1088/0965-0393/13/7/007 Rafii-Tabar, H.: Modelling the nano-scale phenomena in condensed matter physics via computer-based numerical simulations. Phys. Rep. 325, 239–310 (2000). doi:10.1016/ S0370-1573(99)00087-3 Rodney, D., Phillips, R.: Structure and strength of dislocation junctions: an atomic level analysis. PRL 82, 1704–1707 (1999). doi:10.1103/PhysRevLett.82.1704 Rudd, R.E.: Coupling of length scales in MEMS modelling: the atomic limit of finite elements. Int. Soc. Opt. Eng. 4019, 16–25 (2000) Rudd, R.E., Broughton, J.Q.: Coarse-grained molecular dynamics and the atomic limit of finite elements. Phys. Rev. B. 58, 10 (1998). doi:10.1103/PhysRevB.58.R5893
123
110 Rudd, R.E., Broughton, J.Q.: Coupling of length scales and atomistic simulation of MEMS resonators. Int. Soc. Opt. Eng. 3680, 104–113 (1999) Rudd, R., Broughton, J.: Concurrent coupling of length scales in solid state systems. Phys. Status Solidif. 217, 251 (2000). doi:10.1002/(SICI)1521-3951(200001)217:1\251 ::AID-PSSB251[3.0.CO;2-A Rudd, R.E., Broughton, J.Q.: Coarse-grained molecular dynamics: nonlinear finite elements and finite temperature. Phys. Rev. B. 72, 144104 (2005) Shenoy, V., Miller, R., Tadmor, E., Phillips, R., Ortiz, M.: Quasicontinuum models of interfacial structure and deformation. Phys. Rev. Lett. 80, 742–745 (1998). doi: 10.1103/PhysRevLett.80.742 Shenoy, V., Miller, R., Tadmor, E., Rodney, D., Phillips, R., Ortiz, M.: An adaptive finite element approach to atomicscale mechanics–the quasicontinuum method. J. Mech. Phys Solids 47, 611–642 (1999a). doi:10.1016/S00225096(98)00051-9 Shenoy, V., Miller, R., Tadmor, E., Rodney, D., Phillips, R., Ortiz, M.: An adaptive finite element approach to atomicscale mechanics—the quasicontinuum method. J. Mech. Phys. Solids 47, 611–642 (1999b). doi:10.1016/S00225096(98)00051-9 Shi, D., Feng, X., Hang, H., Huang, Y., Hwang, K.: Multiscale analysis of fracture of carbon nanotubes embedded in composites. Int. J. Fract. 134, 369–386 (2005). doi: 10.1007/s10704-005-3073-1 Shiari, B., Miller, R., Curtin, W.: Coupled atomistic/discrete dislocation simulations of nanoindentation at finite temperatures, J. Eng. Mater. Technol. 127, 358–368 (2005) Shilkrot, L.E., Miller, R.E., Curtin, W.A.: Multiscale plasticity modeling: coupled atomistics and discrete dislocation mechanics. J. Mech. Phys. Solids. 52, 755–787 (2004). doi:10.1016/j.jmps.2003.09.023 Shlarl, B., Miller, R.E., Klug, D.D.: Multiscale modeling of solids at the nanoscale: dynamic approach. Can. J. Phys. 86, 391–400 (2008). doi:10.1139/P07-145 Spencer, A.J.M.: Continuum Mechanics. Dover Publications, Dover edition, NY (2004) Srivastava, D., Wei, C.: Nanomechanics of carbon nanotubes and composites. Appl. Mech. Rev 56, 2 (2003). doi: 10.1115/1.1538625 Stillinger, F., Weber, T.: Computer simulation of local order in condensed phases of silicon. Phys. Rev. B 31, 8 (1985). doi:10.1103/PhysRevB.31.5262 Sulsky, D., Zhou, S., Schreyer, H.: Application of a particle-incell-method to solid mechanics. Comput. Phys. Commun. 87, 236–252 (1995). doi:10.1016/0010-4655(94)00170-7 Sun, Y.: Influence of nanofillers on bond strength and toughness, University of Toronto, Ph.D. Thesis, 2007 Tadmor, E., Phillips, R., Ortiz, M.: Mixed atomistic and continuum models of deformation in solids. Langmuir 12, 4529–4534 (1996a). doi:10.1021/la9508912 Tadmor, E., Ortiz, M., Phillips, R.: Quasicontinuum analysis of defects in solids. Philos. Mag. A 73(6), 1529–1563 (1996b). doi:10.1080/01418619608243000
123
J. M. Wernik, S. A. Meguid Tadmor, E., Miller, R., Phillips, R., Ortiz, M.: Nanoindentation and incipient plasticity. J. Mater. Res. 14, 2233–2250 (1999). doi:10.1557/JMR.1999.0300 Troya, D., Mielke, S.L., Schatz, G.C.: Carbon nanotube fracturedifferences between quantum mechanical mechanisms and those of empirical potentials. Chem. Phys. Lett. 382, 133– 141 (2003). doi:10.1016/j.cplett.2003.10.068 Tserpes, K.I., Papanikos, P., Labeas, G.: Sp.G. Pantelakis, Multi-scale modeling of tensile behavior of carbon nanotube-reinforced composites. Theor. Appl. Fract. Mech. 49, 51–60 (2008). doi:10.1016/j.tafmec.2007.10.004 Verlet, L.: Computer experiments on classical fluids. I. Thermodynamical properties of lennard-jones molecules. Phys. Rev. 159, 98–103 (1967). doi:10.1103/PhysRev.159.98 Vvedensky, D.: Multiscale modeling of nanostructures. J. Phys. Condens. Matter 16, R1537–R1576 (2004). doi:10.1088/ 0953-8984/16/50/R01 Wagner, G.J., Liu, W.K.: Coupling of atomistic and continuum simulations using a bridging scale decomposition. J. Comput. Phys. 190, 249–274 (2003). doi:10.1016/S00219991(03)00273-0 Wagner, G., Karpov, E., Liu, W.: Molecular dynamics boundary conditions for regular crystal lattices. Comput. Methods Appl. Mech. Eng. 193, 1579–1601 (2004). doi:10.1016/j.cma.2003.12.012 Wang, X.D., Meguid, S.A.: Dynamic interaction between a matrix crack and a circular inhomogeneity with a distinct interphase. Int. J. Solids Struct. 36, 517–531 (1999). doi:10.1016/S0020-7683(98)00039-0 Wang, C., Zhou, G., Liu, H., Wu, J., Qiu, Y., Gu, B., et al.: Chemical functionalisation of carbon nanotubes by carboxyl groups on stone-wales defects: a density functional theory study. J. Phys. Chem. B 110, 10266–10271 (2006). doi:10.1021/jp060412f Wong, C.H., Liew, K.M., He, X.Q., Tan, M.J., Meguid, S.A.: Modeling and simulation of multi-walled carbon nanotubes using molecular dynamics simulation. NSTI. Nanotech. 3, 2004 (2004) Yakobson, B.I., Brabec, C.J., Bernhole, J.: Nanomechanics of carbon tubes: instabilities beyond linear response. Phys. Rev. Lett. 76, 14 (1996). doi:10.1103/PhysRevLett.76. 2511 Zeng, Q.H., Yu, A.B., Lu, G.Q.: Multiscale modeling and simulation of polymer nanocomposites. Prog. Polym. Sci. 33, 191–269 (2008). doi:10.1016/j.progpolymsci.2007.09. 002 Zhao, X., Cummings, P.T.: Molecular dynamics study of carbon nanotube oscillators revisited. J. Chem. Phys. 124, 134–705 (2006) Zhao, J., Ding, Y.: Silicon carbide nanotubes functionalized by transition metal atoms: a density-functional study. J. Phys. Chem. C 112, 2558–2564 (2008). doi:10.1021/jp073722m Zienkiewicz, O.C.: The Finite Element Method, vol 1–2, 4th edn. McGraw-Hill, London (1991)
