0
Review Articles

Finite-Temperature Quasi-Continuum

[+] Author and Article Information
E. B. Tadmor

Professor
Department of Aerospace Engineering and Mechanics,
The University of Minnesota, Twin Cities,
Minneapolis, MN 55455
e-mail: tadmor@aem.umn.edu

F. Legoll

Professor
Laboratoire Navier,
Université Paris-Est,
Ecole des Ponts ParisTech,
77455 Marne La Vallée Cedex 2, France;
INRIA Rocquencourt,
MICMAC Project-Team,
78153 Le Chesnay Cedex, France
e-mail: legoll@lami.enpc.fr

W. K. Kim

Department of Aerospace Engineering and Mechanics,
The University of Minnesota, Twin Cities,
Minneapolis, MN 55455
e-mail: kwk4814@gmail.com

L. M. Dupuy

CEA, DEN, DMN, SRMA, LC2M,
91191, Gif-sur-Yvette, France
e-mail: laurent.dupuy@cea.fr

R. E. Miller

Professor
Department of Mechanical and Aerospace Engineering,
Carleton University,
Ottawa, ON, K1S 5B6, Canada
e-mail: ron_miller@carleton.ca

It is worth noting that the same problem arises in conventional dynamic finite element calculations, where the discreteness of the mesh leads to anisotropic wave dispersion and reflection in nonuniform grids [11,12].

See; for example, the multiscale model of [22], which was modeled on the “stadium boundary conditions” of [23] and the Langevin thermostat employed by Holland and Marder [24] in earlier, purely atomistic simulations. Also, Shiari et al. [25,26] have proposed a simpler type of damping band whereby a conventional Nosé–Hoover thermostat is applied in the damping region. Damping-band methods are related to so-called “partial thermostating” as discussed in [27].

Note that Kimmer and Jones [45] have recently done a critical comparison of the computational cost and accuracy of both the quasi-harmonic and local harmonic approximations over a range of temperatures. This work sheds further light on the level of approximation inherent in our formulation, and is largely consistent with our observations discussed in Sec. 2. Khoei et al. [46] have also examined the stability (as a function of strain) of the temperature-dependent Cauchy–Born rule, showing that care must be taken at large strains where the eigenvalues of the crystal's dynamical matrix can become negative.

We distinguish between vibrational displacements from the changing mean position (w=q-q¯) and absolute displacement from a fixed reference position (u=q-Q).

Note that we do not make explicit the dependence of Z or A on q¯ since we do not consider q¯ to be a free parameter. As pointed out above, a specific value for q¯ is determined later for a given temperature and boundary condition.

It is important to note that, in general, rigid-body modes must be removed from Φc when computing the free energy in Eq. (24)—otherwise the determinant is simply zero [6]. For the coarse-grained system considered here, this is not an issue since the presence of the atomistic region atoms precludes this possibility.

The reader is referred to Ref. [56] for the study of a related question, in the one-dimensional setting, under a different regime, namely when Nc>>Nat.

In standard QC parlance, the repatoms in the atomistic region are called nonlocal repatoms and those in the continuum region are local repatoms.

There is no need to use interpolation to obtain the positions of atoms in the atomistic region, since all atoms are represented there. In fact, the finite element mesh is not needed in the atomistic region (except for the elements spanning the interface with the continuum region). It is retained to facilitate meshing and to enable easy expansion and shrinking (i.e., adaption) of the atomistic region.

For a two-dimensional QC model, the volume of an element is equal to its area multiplied by the length of the periodic cell in the out-of-plane direction.

We note that there is no contradiction between allowing macroscopic dynamics and the assumption of (metastable) thermodynamic equilibrium inherent in the free energy calculations used in the continuum region. Continuum mechanics is a dynamical theory which enables phenomena such as stress waves. However, it is couched in terms of field variables such as “mean positions” and “temperature” that are only defined under conditions of thermodynamic equilibrium. This apparent contradiction is addressed by the assumption of “local thermodynamic equilibrium” that ensures that the variation of continuum fields in space and time are sufficiently small and slow to allow the system locally to sample its restricted phase space. Naturally, this places certain limitations on the sort of dynamical processes that can be explored in the continuum region.

As pointed out above, we begin the derivation of the hot-QC-dynamics method with the quasi-harmonic expression for the free energy (as opposed to the local harmonic free energy) since this makes it easier to estimate the accuracy of the derivation. Eventually, to obtain a practical formulation, we will revert to the local harmonic approximation.

Note that in the development of the hot-QC method we implicitly assumed that detΦ=detΦc is always positive, which is valid at low or moderate temperatures. However, in hot-QC-dynamic, we allow for the thermal vibration of repatoms in the continuum region. As a result, the instantaneous deformation gradients in some elements can be far from their equilibrium values causing det Φ0. In order to avoid the singularity that occurs when this happens, we introduce a small positive number δ2 into the argument of the logarithm in Eq. (103). In practice δ2 can be computed by calculating det Φ in the undeformed configuration and dividing it by a large number.

Note that Q is dimensionless, P has units of energy × time, and M has units of energy × time 2, therefore strictly speaking Q, P, and M are not actually variables of position, momentum and mass. Yet they play a similar role in the equations.

The QC potential in Eq. (86) involves the free energy given by Eq. (42), which is the sum of a 0 K contribution and of a temperature-dependent contribution. The former is exactly the 0 K QC potential whose gradient may be corrected to remove so-called ghost forces [3] which introduce spurious effects at the interface between the local and the nonlocal regions.

The masses for the Nosé–Hoover chain are selected based on the natural frequency of the system and through trial and error (see Ref. [79] for details). In the current simulations, the same value of 8.62 × 10−5eV·ps2 was used for all links in the Nosé–Hoover chain. Note that similarly to the friction coefficient in Langevin thermostat, the masses of the Nosé–Hoover chain thermostat only affect the equilibration time scale. Thus the simulation results eventually converge to the equilibrium values regardless of the mass values.

The temperature uniformity problem described in Sec. 4.1.2 for the Nosé–Poincaré method does not occur in the thermal expansion simulations because the meshes are uniform. We therefore use this approach to take advantage of its faster convergence.

Ideally, one would like the computational models to exhibit no mesh dependence and to predict the same thermal expansion along all symmetry-related directions regardless of the mesh structure. In reality, due to the approximations inherent to these methods, some degree of mesh dependence becomes apparent at increasing temperature.

The MD simulation was performed using 4000 atoms in a cubic periodic simulation box with the sides aligned with the [100], [010], and [001] directions. The temperature was controlled using a Nosé–Hoover thermostat and zero pressure conditions were imposed using a Parrinello–Rahman barostat (see Ref. [6] for details on the thermostat and barostat). The system was equilibrated for 10 ps (10,000 time steps) and then averages were taken for the three periodic lengths of the box over the next 20 ps (20,000 time steps).

Recall that both the static and dynamic variants of hot-QC involve the assumption that the nodes in the continuum region occupy stationary positions corresponding to a minimum of the effective free energy of the system. The vibrations of the nodes mean that they incorrectly carry kinetic energy (see Eq. (66)) and leads to the spurious entropy contribution which the entropy correction is attempting to remove. The fact that the vibrations of the nodes are reduced by the mesh entropy correction supports this interpretation.

It is well-known that the time step in integrating the Langevin dynamics must be reduced when the friction coefficient is increased. It is standard to require that the time step Δt is much smaller than 1/ξ0 (for ξi=ξ0mi) [83]. In our simulations Δt = 10−3 ps, while the smallest value of 1/ξ0 is 0.1 ps, so this condition is satisfied.

1Corresponding author.

Manuscript received May 6, 2012; final manuscript received September 7, 2012; published online xx xx, xxxx. Editor: Harry Dankowicz.

Appl. Mech. Rev 65(1), 010803 (Mar 21, 2013) (27 pages) Paper No: AMR-12-1028; doi: 10.1115/1.4023013 History: Received May 06, 2012; Revised September 07, 2012

A generalization of the quasi-continuum (QC) method to finite temperature is presented. The resulting “hot-QC” formulation is a partitioned domain multiscale method in which atomistic regions modeled via molecular dynamics coexist with surrounding continuum regions. Hot-QC can be used to study equilibrium properties of systems under constant or quasistatic loading conditions. Two variants of the method are presented which differ in how continuum regions are evolved. In “hot-QC-static” the free energy of the continuum is minimized at each step as the atomistic region evolves dynamically. In “hot-QC-dynamic” both the atomistic and continuum regions evolve dynamically in tandem. The latter approach is computationally more efficient, but introduces an anomalous “mesh entropy” which must be corrected. Following a brief review of related finite-temperature methods, this review article provides the theoretical background for hot-QC (including new results), discusses the implementational details, and demonstrates the utility of the method via example test cases including nanoindentation at finite temperature.

Copyright © 2013 by ASME
Your Session has timed out. Please sign back in to continue.

References

Figures

Grahic Jump Location
Fig. 1

A QC model for an atomically-sharp crack tip in a bcc iron crystal (the horizontal axis is oriented along the [100] direction and the vertical axis along [010]). Frame (a) shows the atoms in the vicinity of the crack tip. The dark gray square is the “atomistic region” of QC. The continuum region encompasses all atoms outside this square and extends beyond boundaries of the figure. Frame (b) shows the corresponding QC model. Atoms retained in the model (so called “representative atoms” or “repatoms”) are shown as filled circles. These repatoms serve as the nodes of a finite element mesh.

Grahic Jump Location
Fig. 2

Hot-QC mesh used to study thermalization and temperature distribution. In the outside zone, the continuum approximation is applied and the repatoms are local. In the central square, the repatoms are either local (simulation 1) or nonlocal (simulations 2 and 3).

Grahic Jump Location
Fig. 3

For Nosé–Poincaré dynamics, evolution as a function of t of (a) Arep(t), (b) Brep(t), and (c) C(t) for different choices of M (note the different scales for the horizontal axes in the three figures)

Grahic Jump Location
Fig. 4

For Langevin dynamics, evolution as a function of t of (a) Arep(t), (b) Brep(t), and (c) C(t) for different values of ξ0 (with the choice ξi=ξ0mi). Note the different scales for the horizontal axes in the three figures.

Grahic Jump Location
Fig. 5

For Langevin dynamics, evolution as a function of t of err(t) defined by Eq. (162) for different choices of ξ0 (with the choice ξi=ξ0mi)

Grahic Jump Location
Fig. 6

Evolution as a function of t of (a) Arep(t), (b) Brep(t), and (c) C(t) for different thermalizing models

Grahic Jump Location
Fig. 7

Distribution of Ai(tmax) for four thermalizing schemes at the end of the simulation: (a) Nosé–Poincaré method with M = 0.01kBT = 2.585 × 10−4eV·ps2, (b) Langevin dynamics with ξi = 10mi, (c) Langevin dynamics with ξi=mi, and (d) Langevin dynamics with ξi = 1. The legend in frame (a) applies to all frames.

Grahic Jump Location
Fig. 8

Distribution of Bi(tmax) for four thermalizing schemes at the end of the simulation: (a) Nosé–Poincaré method with M = 0.01kBT = 2.585 × 10−4eV·ps2, (b) Langevin dynamics with ξi = 10mi, (c) Langevin dynamics with ξi=mi, and (d) Langevin dynamics with ξi = 1. The legend in frame (a) applies to all frames.

Grahic Jump Location
Fig. 9

Average temperature in the continuum and atomistic regions for (a) Langevin thermostat, (b) Nosé–Poincaré thermostat, (c) Nosé–Hoover chain method, (d) Nosé–Hoover chain method with separate temperature control for the continuum and atomistic regions. The nominal set temperature is T = 100 K.

Grahic Jump Location
Fig. 10

Meshes of the model: (a) 10 × 10 mesh with 121 nodes and 200 elements, (b) 20 × 20 mesh with 441 nodes and 800 elements, (c) 40 × 40 mesh with 1681 nodes and 3200 elements, and (d) 80 × 80 mesh with 6561 nodes and 12,800 elements

Grahic Jump Location
Fig. 11

Lattice parameters as function of temperature: (a) 10 × 10 mesh, (b) 20 × 20 mesh, (c) 40 × 40 mesh, and (d) 80 × 80 mesh

Grahic Jump Location
Fig. 12

Thermal fluctuation in the x-displacement of the bottom left node of the hot-QC-dynamic mesh

Grahic Jump Location
Fig. 13

Mesh of the nanoindentation model: (a) entire system, and (b) zoomed view of the fully-refined mesh under the indenter

Grahic Jump Location
Fig. 14

Force as a function of indentation depth for (a) the Langevin thermostat, and (b) the Nosé–Poincaré thermostat

Grahic Jump Location
Fig. 15

Mesh after dislocation nucleation at T = 100 K. The figure shows a closeup of the region under the indenter (which is embedded in a larger continuum region as shown in the inset). The nucleated dislocation dissociates into a pair of Shockley partial dislocations. See text for discussion.

Grahic Jump Location
Fig. 16

Force versus indentation depth curves for nanoindentation at different temperatures. (a) All curves superposed to demonstrate the elastic softening. (b) Curves shifted horizontally for visibility. Arrows indicate the points on each curve where dislocation nucleation first occurs.

Grahic Jump Location
Fig. 17

Dislocation nucleation load as a function of temperature. The bars indicate the standard error in the measurement obtained from 10 simulations at the same temperature. The dashed line is an analytical fit based on a modified Tomlinson model. See text for details.

Tables

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In