In this work, we present a concurrent atomistic-continuum (CAC) method for modeling and simulation of crystalline materials. The CAC formulation extends the Irving-Kirkwood procedure for deriving transport equations and fluxes for homogenized molecular systems to that for polyatomic crystalline materials by employing a concurrent two-level description of the structure and dynamics of crystals. A multiscale representation of conservation laws is formulated, as a direct consequence of Newton's second law, in terms of instantaneous expressions of unit cell-averaged quantities using the mathematical theory of distributions. Finite element (FE) solutions to the conservation equations, as well as fluxes and temperature in the FE representation, are introduced, followed by numerical examples of the atomic-scale structure of interfaces, dynamics of fracture and dislocations, and phonon thermal transport across grain boundaries. In addition to providing a methodology for concurrent multiscale simulation of transport processes under a single theoretical framework, the CAC formulation can also be used to compute fluxes (stress and heat flux) in atomistic and coarse-grained atomistic simulations.

Solid state physics describes the structure of all crystals in terms of primitive unit cells. These primitive unit cells, each containing an identical group of atoms, are the building blocks of crystals. When added continuously in space, they form the structure of a crystal and completely fill the space the crystal occupies,1 cf. Fig. 1.

FIG. 1.

A two-level structural description of crystals: a continuously distributed lattice + a discrete basis.

FIG. 1.

A two-level structural description of crystals: a continuously distributed lattice + a discrete basis.

Close modal

Although the motion of an atom in a crystal is an irregularly fluctuating function of time, it turns out that the dynamics of a crystal can be most readily described, not in terms of individual atoms, but in terms of traveling waves, named lattice vibrations by Max Born. These vibrational modes are called phonons and are of two types: acoustic and optical, cf. Fig. 2. With the former, atoms in a unit cell move in the same phase, resulting in the deformation of the lattice; with the latter, atoms undergo relative motion within the lattice cells, leaving the lattice unchanged. Thus, the displacements of atoms in a crystal can be described as the sum of continuous lattice displacements and the internal displacements of the atoms relative to the lattice.2 

FIG. 2.

Transverse atomic motion in (a) an acoustic mode and (b) an optical mode.

FIG. 2.

Transverse atomic motion in (a) an acoustic mode and (b) an optical mode.

Close modal

It is seen that the solid state physics description of the structure and dynamics of crystals is based on a concurrent atomistic-continuum (CAC) approach: the structure is continuous at the lattice level but discrete at the atomic level; the lattice deformation is continuous until a defect is generated, cf. Fig. 3, while the internal motion is discrete.

FIG. 3.

Schematics illustrating that (a) the lattice deformation is continuous until (b) structural discontinuities generated by slip and/or twinning.

FIG. 3.

Schematics illustrating that (a) the lattice deformation is continuous until (b) structural discontinuities generated by slip and/or twinning.

Close modal

Historically, atomistic and continuum descriptions offer two fundamentally different approaches to our understanding of crystalline matter. From the atomistic viewpoint, matter consists of discrete particles; while from the continuum viewpoint, matter is infinitely divisible. The continuum view leads to formulations of field equations of conservation laws. Supplemented by constitutive relations, such as Hooke's law and Fourier's law, these conservation equations serve as the governing equations in continuum modeling of materials.

This work introduces the concurrent atomistic-continuum (CAC) method that links and unifies the atomistic and continuum descriptions of crystalline materials. The Irving-Kirkwood procedure for linking the molecular description to hydrodynamical equations is reviewed in Sec. II; the CAC formalism is introduced in Sec. III; the finite element (FE) formulation of CAC as well as the fluxes and temperature in the FE representation is outlined in Sec. IV; numerical examples are presented in Sec. V, followed by a summary in Sec. VI. Detailed derivations of the conservation equations using the mathematical theory of distributions for monatomic crystals are included in  Appendix A and that for polyatomic crystals in  Appendix B.

The statistical mechanics formulation to link discrete and continuum descriptions of conservation laws was pioneered by Irving and Kirkwood (IK) in 1950.3 In their landmark paper, IK defined local densities in terms of molecular variables using the Dirac delta for “localized density,”4 “since mass or momentum of any molecule may be considered as localized at that molecule.”3 A local density of a dynamic phase function A(rk,vk), where rk andvk are the position and velocity of the kth particle, is defined as an ensemble-averaged point function a¯(x,t) as

a¯(x,t)k=1NA(rk,vk)δ(rkx),
(1)

where the symbol “≡” means “equal by definition.”

Equation (1) is usually interpreted as mapping the phase functions in the 6N-dimensional phase space to the local densities in the 4-dimensions of physical space and time. Based on the local density definition in Eq. (1), IK derived the rate of change for the densities of mass, momentum, and energy in the form of partial differential equations that relate the density of a conserved quantity a¯(x,t) and its flux J¯(x,t) as

ta¯(x,t)=xJ¯(x,t).
(2)

Formulas for fluxes (stress and heat flux) are then obtained through comparing the rate equations with the differential equations of hydrodynamics.

The IK formalism has inspired numerous efforts to link molecular variables to field quantities in classical continuum mechanics.5–8 To understand the motivation of a new formalism, we recall the following noted by Irving and Kirkwood3 and by Kirkwood.9 

  1. Local densities and fluxes in the IK formulation are point functions. “The point functions, although averaged neither over space nor time, satisfy equations that are identical in form to the equations of hydrodynamics for a single component, single phase system.”3 These equations were called “the hydrodynamical-like equations” in the IK paper.3 It was noted that “to obtain the hydrodynamical equations themselves it is merely necessary to perform the appropriate space and time averages.”3 

  2. The densities for mass, momentum, and energy were defined without restriction to a single component single phase system. For fluxes (stress and heat flux), IK stated that “we must impose this restriction” and the flux formulas obtained with their formalism are “only valid for a single component, single phase system.”3 

  3. The statistical mechanical theory of transport processes published in a series of 14 papers by Kirkwood and co-workers only considers rigid molecules, ignoring the internal degrees of freedom (DOFs) of the molecules. In the first one of the series, Kirkwood envisioned extension of the formulation to molecules possessing internal degrees of freedom (DOFs).9 

The CAC formalism is an extension of the IK procedure for homogeneous hydrodynamic systems to a two-level structural description of crystals that includes the atomic DOFs within each lattice cell.10,11 With such an extension, the CAC formulation is applicable to polyatomic crystalline materials and can be used to solve atomic trajectories or quantify instantaneous fluxes, conferring full characteristics of a concurrent atomistic-continuum method.

For a polyatomic crystal of Nl unit cells, with each unit cell containing Nα atoms, the density of a phase function A(rkξ,vkξ) at point r in the physical space can be expressed as

a¯(r,t)=k=1Nlξ=1NaA(rkξ,vkξ)δ(rkξr),
(3)

where rkξ and vkξ are the position vector and velocity vector of the ξth atom in the kth unit cell, and the symbol ⟨ ⟩ denotes ensemble averaging that involves NlNa fold integrals of the phase function A and a probability distribution function in the phase space.

Recall that the structure and deformation of a crystal are continuous in space at the lattice level but are discrete at the atomic level. It is, therefore, consistent with solid state physics to describe the crystal structure in terms of the positions of lattice cells and the relative positions of atoms inside the lattice cells. For this purpose, we resolve the sum in the definition of the local density in Eq. (3) into contributions from groups of atoms within the same lattice cell and express the locations relative to the unit cell for those atoms, i.e.,

a¯(x,y,t)k=1Nlξ=1Na1A(rkξ,vkξ)δ(xrk)δ(yΔrkξ),
(4)

where rk is the position of kth unit cell and Δrkξ is the relative position of ξth atom in the unit cell k.

The density in Eq. (4) is expressed in terms of the direct product of two deltas, with the two field variables x and y representing a distinction of large- and small-scale variations of the variable r in Eq. (3). Note that the ensemble average in Eq. (4) also involves NlNa fold integrals and, although expressed in terms of different phase variables, the values of a local density defined in Eqs. (3) and (4) can be mathematically equivalent. However, different from that in Eq. (2), a conservation equation for the density defined in Eq. (4) relates a conserved quantity to two fluxes, J¯1 and J¯2, owing to the use of large- and small-scale variations in the latter, i.e.,

ta¯(x,y,t)=xJ¯1(x,t)+yJ¯2(x,y,t).
(5)

Evans and Morriss argued that, since conservation laws hold instantaneously, conservation equations and fluxes should be definable without ensemble averaging.12 This has been demonstrated recently for monatomic crystals.13–16 To derive the conservation equations and flux expressions that hold instantaneously based on the two-level structural description of crystals, we consider a unit cell of volume V located at point x in the physical space; within the unit cell, there are Nα volume elements, each containing one atom and having volume Vα satisfying αVα=V. Since the distribution of the unit cells is continuous and homogeneous in space, we may define a local density per unit-cell volume as

a(x,t)k=1Nl(ξ=1NaA(rkξ,vkξ))δV(xrk)=α=1Naaα(x,t),
(6)

where aα(x,t)=kA(rkα,vkα)δV(xrk) is the contribution of the αth atom to a(x,t) and can be expressed as

aα(x,t)=k=1Nlξ=1NaA(rkξ,vkξ)δV(xrk)δVα(yΔrkξ)=a(x,y,t)
(7)

and

δV(xrk)1V{1ifxrkV,0ifxrkV,
(8)
δVα(yΔrkξ){1ifyΔrkξVαorξ=α,1ifyΔrkξVαorξα.
(9)

Since α=1NαδVα(yΔrkξ)=1, one can see that Eq. (7) satisfies Eq. (6), and a(x,y,t)=aα(x,t) holds for all yVα.

Note that the box functions defined in Eqs. (8) and (9) have a jump discontinuity at the enclosing surfaces V and Vα, respectively. These functions do not have derivatives in the classical sense. It thus requires a generalization of the concept of a function. The corresponding mathematical theory is known as the theory of distributions.17 This theory will be used throughout this work as the mathematical tool to derive the conservation laws for physical quantities described by nondifferentiable functions.

In atomistic simulations, the equation of motion is solved in discrete time steps. Accordingly, local densities should be further averaged over a time-interval. The contribution of the αth atom to the mass density ρα, momentum density pα(pα=ραvα), and energy density Eα(Eα=ραeα) per unit volume can then be defined as

ρα(x,t)k=1Nlξ=1Namξδ¯V(xrk)δ¯Vα(yΔrkξ),
(10)
pα(x,t)k=1Nlξ=1Namξvkξδ¯V(xrk)δ¯Vα(yΔrkξ),
(11)
Eα(x,t)k=1Nlξ=1Na1(12mξvkξ2+Φkξ)δ¯V(xrk)δ¯Vα(yΔrkξ),
(12)

where Φkξ is the potential energy of atom ; δ¯V denotes the average of δV over a time-interval Δt, i.e., δ¯V(rkx)=1Δt0ΔtδV[rk(t+τ)x]dτ. The time evolution of a conserved quantity, i.e., the conservation equation, can then be expressed in terms of a lattice-level flux J as

ta(x,t)=1VVJ(x+x,t)ndS
(13)

or in terms of atomic-level fluxes J1 and J2 as

taα(x,t)=1VVJ1(x+x,t)ndS1VVαJ2(x,y+y,t)ndSα,
(14)

where n is the outward unit normal vector to the surface element, J1 is the part of flux that flows across the enclosing surface ∂V, and J2 is the part of flux that flows back and forth inside V and hence only cross the bounding surface ∂Vα, as shown in Fig. 4.

FIG. 4.

2D schematic for J1, the flux across the bounding surface of the unit cell (the rectangle) ∂V, and J2, the flux that does not cross ∂V but rather ∂Vα (the triangle) within the unit cell.

FIG. 4.

2D schematic for J1, the flux across the bounding surface of the unit cell (the rectangle) ∂V, and J2, the flux that does not cross ∂V but rather ∂Vα (the triangle) within the unit cell.

Close modal

It is noted that by integrating Eq. (5) over V and invoking the divergence theorem, one obtains an integral equation in the form of Eq. (14). However, Eq. (5) is not valid instantaneously and it requires the fluxes to be differentiable. By contrast, the conservation equation in Eq. (14) holds instantaneously in the distributional sense and the fluxes are not required to be continuous functions.

Averaging a local density over a unit cell of volume V leads to the integral form of conservation laws that equate the rate of the change of conserved quantities in the unit cell to their fluxes across the enclosing surface of the unit cell. Following from Eq. (10), the equation of conservation of mass can be derived in the form of Eq. (14) as

ραt=k=1Nlξ=1Namkξtδ¯V(xrk)δ¯Vα(yΔrkξ)+k=1Nlξ=1Namkξδ¯V(xrk)tδ¯Vα(yΔrkξ)=1VVραvndS1VVαραΔvαndSα,
(15)

where v=x˙, vα=pα/ρα, and Δvα=vαv.

The momentum conservation law can also be derived, using the momentum density defined in Eq. (11), from Newton's second law, as

(ραvα)t=k=1Nlξ=1Namξv˙kξδ¯V(xrk)δ¯Vα(yΔrkξ)+k=1Nlξ=1Namξvkξtδ¯V(xrk)δ¯Vα(yΔrkξ)+k=1Nlξ=1Namξvkξtδ¯V(xrk)tδ¯Vα(yΔrkξ)=1VV(tαραvαvn)dS+1VVα(ταραvαΔvαn)dSα,
(16)

where tα and τα are the momentum flux (stress) that cross ∂V and ∂Vα, respectively, as a result of atomic interaction and motion involving atom α. The total stress vector t on a surface element An of area A and normal n, located at position y in the unit cell at x, from all the atoms can then be expressed as

t(x,y,n)=α=1Na(tα+τα)=tpot(x,y,n)+tkin(x,y,n),
(17)

where tpot is the potential part of t that measures the interaction force per unit area transmitted across An, i.e.,

tpot(x,y,n)=k,l(k)Nlξ,ηNaΦkξrlηrlηrkξδ¯An(φx)dφ+kNlξ,ηξNaVδ¯V(xrk)Φkξrkηrkηrkξδ¯An(φy)dφ,
(18)

and tkinis the kinetic part of t that measures the flow of momentum per unit area and time across An, i.e.,

tkin(x,y,n)=k=1Nlmv~kv~kδ¯An(xrk)k=1NlVδ¯V(xrk)α=1NamαΔv~kαΔv~kαδ¯An(yΔrkα).
(19)

Here, v~k=vkv, Δv~kα=ΔvkαΔvα, rlηrkξδ¯An(φx)dφ, defined in Eq. (A3), and δ¯An(yΔrkα), defined in Eq. (A4), represent the potential and kinetic flux as a line-plane intersection problem in space and time, respectively.

Similarly, denoting qα and jα as the parts of heat flux that cross ∂V and ∂Vα, respectively, the time rate of change of the energy density Eα, can be derived as

(ραeα)t=tk=1Nlξ=1NaEkξδ¯V(xrk)δ¯Vα(yΔrkξ)=1VV(qα+tαvαvραeα)ndS+1VVα(jα+ταvαΔvαραeα)ndSα.
(20)

The total heat flux is then given by

q(x,y)=α=1Na(qα+jα)=qpot(x,y)+qkin(x,y)=q(x,y,ei)ei,
(21)

where ei (i = 1, 2, 3) are the orthonormal basis at x,

qpot(x,y,ei)=k,lkNlα,ηNaΦkαrlηv~kαrlηrkαδ¯Ai(φx)dφ+kNlVδ¯V(rkx)α,ηαNaΦkαrkηv~kαrlηrkαδ¯Ai(φy)dφ
(22)

and

qkin(x,y,ei)=k=1Nlα=1Nα(12mαv~kα2+Φkα)v~kδ¯Ai(xrk)k=1NlVδ¯V(rkx)αNa(12mαv~kα2+Φkα)Δv~kαδ¯Ai(yΔrkα).
(23)

The link between atomistic and continuum descriptions at this point has fundamentally departed from the IK formalism. As distinct from classical continuum mechanics, this new formulation describes the structure and properties of a crystalline system as a continuous function in x at the lattice level, cf. Figs. 3 and 5, while discrete in α or y at the atomic level; fluxes (stress and heat flux) are composed of components resulting from the motion and deformation of lattice cells (denoted as t and q, respectively) and the rearrangement of atoms within the cells (denoted as τ and j, respectively). Most importantly, instead of linking molecular variables to classical continuum mechanics equations, the CAC formalism leads to a new representation of conservation laws and fluxes.

FIG. 5.

Schematic of the mass density distribution in space for a crystal with two different atoms in the primitive unit cell: (a) using Eq. (3), [(b) and (c)] using Eq. (7), showing that the mass density is a continuous and homogeneous function at the lattice level, while it is discrete and inhomogeneous at the atomic scale.

FIG. 5.

Schematic of the mass density distribution in space for a crystal with two different atoms in the primitive unit cell: (a) using Eq. (3), [(b) and (c)] using Eq. (7), showing that the mass density is a continuous and homogeneous function at the lattice level, while it is discrete and inhomogeneous at the atomic scale.

Close modal

An important difference between atomistic and classical continuum theories is the definition of temperature; in the former, temperature is a derived quantity, while in the latter, temperature is a fundamental state variable. Specifically, the abiding definition of temperature in classical molecular simulations is the kinetic temperature defined by the kinetic theory in terms of the mean kinetic energy of the random motion of atoms, while that in quantum calculations is based on the mean energy of phonons.

The kinetic theory provides a microscopic description of temperature through the classical mechanics of particles. The equipartition theorem of the kinetic theory asserts that the average kinetic energy of a particle in an equilibrium gaseous system with temperature T is 12mvkα2=32kBT. Away from thermal equilibrium, the particle velocity is usually replaced by v~kα, the difference between particle velocity and the velocity field. Such a temperature definition is consistent with the physical picture of an ideal-gas thermometer in which the velocity is measured relative to the co-moving frame of the kinetic thermometer.18 The nonequilibrium temperature is expressed19 as

32kBTN(x)=k=1Nl12mαv~kα2Vδ¯V(xrk).
(24)

This temperature definition, however, has unambiguous meaning only in thermal equilibrium or an approximate local thermal equilibrium condition such as that in a steady state. To elucidate this point, we decompose the kinetic energy density into two parts,

Eαkin(x,t)=12k=1Nlξ=1Namξ(vkξ)2δ¯V(xrk)δ¯Vα(yΔrkξ)=12ρα(vα)2+12k=1Nlmα(v~kα)2δ¯V(xrk)=kα1+kα2,
(25)

where kα1 is the kinetic energy density represented by the velocity field and kα2 is the kinetic energy due to the velocity difference between particle velocity and velocity field, usually interpreted as the random motion of particles. Although the definition of total kinetic energy is unambiguous, the velocity field can only be uniquely defined when it is zero (for thermal equilibrium) or a constant (for steady state). This is because the velocity field depends on length and time scales, i.e., the averaging, of the linear momentum. Thus, regardless of whether it is in an atomistic or a multiscale simulation, temperature can only be well defined at thermal equilibrium or steady state. It is a concept describing the state of a system rather than the dynamics of an atom.

Differing from the classical concept, the quantum description of temperature is expressed in terms of phonons. There are more phonons at high temperature and fewer phonons at low temperature. The average phonon number is determined by the equilibrium temperature T through the Planck constant ħ and the Boltzmann constant kB. This relationship is called the Bose-Einstein distribution, i.e.,

n(κ,ν)=[exp(ω(κ,ν)kBT)1]1,
(26)

where n(κ,ν) is the phonon number for branch ν and wave vector κ in thermal equilibrium. The total kinetic energy of an equilibrium system can then be expressed in terms of all the available phonon modes in the systems, i.e.,

K=12κNlν3Naω(k,ν)[12+n(k,ν)].
(27)

Recall that the number of phonon modes is equal to the number of the particle degrees of freedom of the system and that the total kinetic energy per phonon mode in the classical limit is equal to kBT/2. The difference between the quantum and classical temperatures can thus be quantified through equating the phonon and classical descriptions of the total kinetic energies,20 i.e.,

K=12κNlν3Naω(k,ν)[12+n(k,ν)]=12κNlν3Naω(k,ν){12+[exp(ω(k,ν)kBTquantum)1]1}=32NlNakBTclassical.
(28)

For a given quantum temperature, the corresponding classical temperature can be calculated using Eq. (28). Figure 6 presents such calculations for LiF and SrTiO3. As can be seen from Fig. 6, the two temperatures are different at low temperature, but the difference becomes negligible as the temperature increases. At 0 K quantum temperature, LiF and SrTiO3 have a zero-point energy equivalent to 237.1 K and 279.7 K classical temperature, respectively, and hence there are corresponding motions at 0 K, i.e., the zero-point vibrations.

FIG. 6.

Quantum temperature vs classical temperature (blue solid curves) for (a) LiF and (b) SrTiO3, showing that the two temperatures converge at high temperature, indicated as the solid blue curves approaching the dotted green lines.

FIG. 6.

Quantum temperature vs classical temperature (blue solid curves) for (a) LiF and (b) SrTiO3, showing that the two temperatures converge at high temperature, indicated as the solid blue curves approaching the dotted green lines.

Close modal

As shown in Fig. 5, a local density such as the mass density in a polyatomic system is not a continuous function at the atomic scale. By decomposing the atomic position into the position of the lattice cell, x, and its internal position within the unit cell at x, we are able to express a density aα(x,t) or a(x, y, t) at an equilibrium or steady state system as a homogeneous and continuous function in x. Consequently, the field equations of the conservation laws can be discretized in x and solved using continuum-based numerical methods such as the finite element (FE) method.

As shown in Eq. (25), the decomposition of the total kinetic energy into kα1 and kα2 depends on the length scale of the velocity field. In particular, in a FE model kα1 is represented by FE nodal velocities, in which the displacement field in an element is interpolated with FE shape functions Sξ(x), i.e.,

u^α(x,t)=Sξ(x)Uξα(t),
(29)

where Uξα(t) is the displacement vector of αth atom embedded within the ξth node of the element; the contribution of the αth atom to the kinetic energy in an element of volume Ve is thus given by

Vekα1dV=Ve12ρα(u˙α)2dV=Ve12ρα(Sξ(x)U˙ξα)2dV.
(30)

The kinetic energy of an element can also be expressed in terms of phonons using the phonon dispersion relations of the FE model.21 However, a FE model employing the usual trilinear shape functions cuts off phonons whose wavelengths are smaller than the element size. The motion of the FEs thus only represents a subset of the phonons of the system. This means that the FE mesh, i.e., the FE displacement approximation, uniquely determines the decomposition of total kinetic energy into kα1 and kα2, where kα1 is represented by FE nodal velocities, while kα2 is related to kα1 through the temperature-dependent phonon number n(κ,ν).

A special case is for systems at high-temperature thermal equilibrium, in which every phonon mode or every particle DOF has the same energy kBT/2. The total kinetic energy of an 8-node element containing nl=Ve/V unit cells is thus given by

Vekα1dV+Vekα2dV=nl×32kBT,
(31)

with

Vekα1(x,t)dV=Ve12ρα[Sξ(x)U˙ξα]2dV=8(32kBT)
(32)

and

Vekα2(x,t)dV=Ve12k=1Nlmα(v~kα)2δ¯V(xrk)dV=(nl8)(32kBT).
(33)

Equations (32) and (33) can be used to prescribe FE nodal velocities for a given initial equilibrium temperature T or to find kα2 from kα1.

Substituting the conservation equations of mass into the linear momentum equation, we obtain

ραu¨α=fαint(x)+1VVtαkindS+1VVαταkindSα.
(34)

Denoting the sum of the two surface integrals of kinetic stresses as fαT(x,t), noting that the averaged fαT(x,t) over a long-time duration can be linked to kα2 and hence the local temperature, we have

1t0tfαT(x,τ)dτ1t0t1VVtαkindS+1t0t1VVαταkindSα=λxT,
(35)

where xT and λ can be determined from x(kα1+kα2), in which xkα1 can be calculated based on FE nodal velocities, and xkα2 can be determined from xkα1 if we assume local thermal equilibrium and employ the assumption used in the linearized phonon Boltzmann equation, i.e.,xnk(x,t)=xnk0(x,t),22 where nk0 and nk are the equilibrium and nonequilibrium phonon mode distribution for wave vector k, respectively.

Equation (34) can be solved using the finite element method, with the weighted residual of Eq. (34) over an element being given by

VeSη(x)(ραu¨αfαintfαextfαT)dV=0.
(36)

This is the weak (Galerkin) form of Eq. (34). The integral can be numerically evaluated using Gaussian quadrature, and the equation can be solved for equilibrium or nonequilibrium processes.

For systems at thermal equilibrium,xT=0, fαT is a fluctuating function of time with a zero mean everywhere except at the system's boundary if thermal expansion is constrained. It can thus be modeled as a constant force applied at the boundary. When atomic-scale fluctuations are important to the simulated phenomena such as thermal activation of dislocation motion or phase transition, fαT can be modeled as a periodic or randomly fluctuating force in time with the frequencies of the phonons that are cut off by the FE shape functions.

For nonequilibrium processes, xT is a constant within an element if the usual linear interpolation is used. fαT can then be modeled as a body force with a constant mean in each element but with periodic or random fluctuation in t. It is noted that, to model it as a random force, fαT needs to be partitioned into a fluctuating part and a damping part, with the latter to be determined by and to balance the fluctuation force, according to the fluctuation-dissipation theorem.23 

Similar to that in MD,13–15,24 fluxes across a surface element in a CAC model can be computed through calculations of the interaction forces and atomic motions that cross the surface element. Not only can the stress and heat flux at any given surface element be calculated, their spatial distributions can also be plotted as contours for the entire model.

With the usual trilinear shape functions to approximate the displacement field in FE, maximum fluxes occur at the nodes of the elements. Thus, only fluxes at the nodes and a few other locations on the boundaries of the elements are needed for the flux contours, such as the unit cells shown in Fig. 7. In regions of structural disorders or discontinuities, such as those near a dislocation core or crack tip, atomic-scale fluxes can be calculated using Eqs. (18) and (19) and plotted as a contour at atomic resolution.

FIG. 7.

(a) 2D schematic of an element containing many unit cells, and each unit cell contains two atoms. The shaded unit cells are those that can be used to calculate fluxes for interpolation of the flux fields in the element so as to generate contours of fluxes within the element; (b) shear stress contour near a dislocation core.

FIG. 7.

(a) 2D schematic of an element containing many unit cells, and each unit cell contains two atoms. The shaded unit cells are those that can be used to calculate fluxes for interpolation of the flux fields in the element so as to generate contours of fluxes within the element; (b) shear stress contour near a dislocation core.

Close modal

Currently, there are two public versions of CAC codes, the CAC package in LAMMPS and PyCAC. The general procedure for a CAC simulation can be summarized as follows:

  • Select an interatomic potential(s) for the system;

  • Discretize the system based on its morphological length scales, using atomic resolution for atomic-scale structural disorders such as grain boundaries (GBs), while using FEs for single crystal regions with the mesh size determined by, e.g., the smallest possible dislocation spacing;

  • Relax the system with initial temperature. One may need to systematically change the initial configurations to obtain the dislocation or GB structures that have the minimal energy;

  • Apply loading and boundary conditions;

  • Run simulation to solve for and to visualize atomic trajectories, which are either atomically resolved or interpolated through FE nodal displacements, for observations of dynamic processes or emerging phenomena, using software such as Ovito for atomic visualization or Tecplot for FE visualization; and

  • Calculate fluxes (stress or heat flux) or other transport properties.

The accuracy and efficiency of the CAC method have been tested through one-to-one comparisons with MD in space- and time-resolved simulations of phenomena such as crack initiation and branching,26,27 phase transitions,28 dislocation nucleation,29–32 loop formations33–36 interactions with interfaces and other defects,37–43 as well as phonon-dislocation interaction44,45 and phonon-GB interaction.46 In general, for an element containing N by N by N unit cells, the ratio of the number of mathematical operations in CAC to that in MD scales as O(N−3).42 This means that CAC can be very efficient by using large elements, e.g., for polycrystalline materials with micrometer-sized grains or for crystals with low dislocation density, i.e., large dislocation spacing.

To illustrate the applicability of CAC, in this section, we present three sets of prior simulation results. These simulations were performed at low temperature without considering the fluctuations from phonons whose wavelengths are smaller than the element size.

CAC has been employed to study the structure and energy of grain boundaries (GBs) and their dynamic interaction with cracks and dislocations in SrTiO3 polycrystals.42,43 In these studies, GBs are atomically resolved, while regions away from GBs are coarse-grained. Simulation results show that both the energy and atomic displacements of SrTiO3 GBs from CAC agree well with those that from MD if the width of the atomically resolved region for the GBs is larger than 12 Å.42 The ghost forces reported at the atomic-continuum interface in many multiscale models are not present in the CAC models.42 The atomic-scale GB structures from CAC42 are shown in Fig. 8 to compare well with those from experimental results.47 

FIG. 8.

CAC results of (a) ∑5(310), (b) ∑13(510), and (c) ∑5(210) GBs in SrTiO3, and (d) STEM images of ∑5(210) GB;47 the unit structures around each GB plane are outlined with solid lines.42 Reproduced with permission from Yang and Chen, Proc. R. Soc. A Math. Phys. Eng. Sci. 471, 20140758 (2015). Copyright 2015 The Royal Society.

FIG. 8.

CAC results of (a) ∑5(310), (b) ∑13(510), and (c) ∑5(210) GBs in SrTiO3, and (d) STEM images of ∑5(210) GB;47 the unit structures around each GB plane are outlined with solid lines.42 Reproduced with permission from Yang and Chen, Proc. R. Soc. A Math. Phys. Eng. Sci. 471, 20140758 (2015). Copyright 2015 The Royal Society.

Close modal

Figure 9 presents CAC simulation results of the misfit dislocation networks in the phase interfaces in heteroepitaxial thin films. CAC is shown to have reproduced the experimentally observed dislocation networks in PbTe/PbSe (001) interfaces.48 The dislocation spacing in the system as a function of the epilayer thickness obtained by CAC is also in good agreement with the experimental measurements. In addition, the critical thickness above which misfit dislocations are observed to form in the heteroepitaxial thin film compares very well with experimental observations.49 

FIG. 9.

CAC simulation results of (a) the dislocation network in a PbTe/PbSe (001) interface and (b) a comparison between CAC simulation results and experimental measurements of the dislocation spacing as a function of PbTe epilayer thickness; the inset in (a) is a STM image of the dislocation network in PbTe/PbSe (001).48 Reproduced with permission from Springholz and Wiesauer, Phys. Rev. Lett. 88, 015507 (2001). Copyright 2001 American Physical Society.

FIG. 9.

CAC simulation results of (a) the dislocation network in a PbTe/PbSe (001) interface and (b) a comparison between CAC simulation results and experimental measurements of the dislocation spacing as a function of PbTe epilayer thickness; the inset in (a) is a STM image of the dislocation network in PbTe/PbSe (001).48 Reproduced with permission from Springholz and Wiesauer, Phys. Rev. Lett. 88, 015507 (2001). Copyright 2001 American Physical Society.

Close modal

Since the only constitutive relation in CAC is the nonlocal force field, and CAC is of an integral form, continuity between FEs is not required. Consequently, nucleation and propagation of dislocations and cracks can be simulated in the FE region via sliding and separation between elements. Simulation results can be output in terms of elements or in terms of atomic trajectories. This makes CAC a convenient tool for the simulation and visualization of the dynamic processes of defect nucleation, interaction, and propagation.

Figure 10 presents CAC simulation results for the loading of a notched solid normal to the notch depth. The simulations have captured the dynamics processes of crack initiation, propagation, and branching. They also reveal that the interactions of the stress waves that emitted from the propagating crack with those that are reflected back by the boundaries of the specimen lead to a sharp increase of the crack speed, which in turn triggers the dynamic instability and hence the branching of cracks in micrometer-scale specimens.26 In addition, we have measured the critical stress intensity factor and observed a crack speed near the speed of Rayleigh wave before crack branching.

FIG. 10.

Time sequence of CAC simulations of a brittle crystal showing (a) stress waves emitting from a propagating crack and [(b)–(d)] crack branching as a result of the interaction between waves propagating from the crack tips and those reflected from the specimen boundaries.50 Reproduced with permission from Deng et al., Int. J. Plast. 26, 1402–1414 (2010). Copyright 2010 Elsevier Ltd.

FIG. 10.

Time sequence of CAC simulations of a brittle crystal showing (a) stress waves emitting from a propagating crack and [(b)–(d)] crack branching as a result of the interaction between waves propagating from the crack tips and those reflected from the specimen boundaries.50 Reproduced with permission from Deng et al., Int. J. Plast. 26, 1402–1414 (2010). Copyright 2010 Elsevier Ltd.

Close modal

CAC has also been used to characterize the complex dynamics of fast moving dislocations, including the energy intensities and the wavelengths of phonons emitted from sonic dislocations, as well as the velocity-dependent stress fluctuations around the core of moving dislocations. Figure 11 presents snapshots comparing dislocation motion simulated using CAC with a uniform mesh and that using MD, with the number of DOFs of the CAC model being 1% of that of MD. The dislocation velocity is measured to be ∼2900 m/s in both CAC and MD simulations. Mach cones in a V-shaped pattern of the phonon wave-fronts are observed in the wake of the sonic dislocations. Analysis of simulation results based on a wavelet transform shows that the faster a dislocation is moving, the longer the emitted phonon wavelength.29–32 

FIG. 11.

Snapshots of dislocation motion in a V-notched specimen; only the left part of the specimen in MD and the right part of the specimen in CAC are displayed. The atoms in MD and the elements in CAC are color coded in displaying stress component σyy. The color, red or blue, respectively, indicates tensile or compressive stress. Reproduced with permission from Xiong et al., Acta Mater. 104, 143–155 (2016). Copyright 2015 Acta Materialia Inc.

FIG. 11.

Snapshots of dislocation motion in a V-notched specimen; only the left part of the specimen in MD and the right part of the specimen in CAC are displayed. The atoms in MD and the elements in CAC are color coded in displaying stress component σyy. The color, red or blue, respectively, indicates tensile or compressive stress. Reproduced with permission from Xiong et al., Acta Mater. 104, 143–155 (2016). Copyright 2015 Acta Materialia Inc.

Close modal

Figure 12 presents CAC simulation results of the sequential slip transfer of a0 ̸2 [110](111) dislocations across a Σ3 (111) coherent twin boundary (CTB) in Ni56 under ostensibly quasistatic strain rate conditions. Simulation results shows that CAC can serve as an effective tool for the study of slip transfer reactions of sequences of incoming and outgoing long-range dislocation pileups of mixed dislocation characters. Five embedded atom method (EAM) potentials51–55 were compared; results demonstrate the uncertainty in computed dislocation-interface reactions associated with the use of a variety of interatomic potentials and suggest that the applicability of dislocation/GB interaction criteria in the literature for slip transfer may have limitations.

FIG. 12.

Snapshots of dislocation pileup with dominant leading screw character impinging against a CTB. Atoms are colored by adaptive common neighbor analysis: red are of hexagonal-close packed local structures, blue are BCC atoms, and all FCC atoms are deleted. In (a), incoming dislocations approach the CTB. In (b), the leading dislocation is constricted at the CTB, where two Shockley partial dislocations are recombined into a full dislocation. All five potentials produced the same results in (a) and (b). In (c), with Mishin-EAM51 and Voter-EAM potentials,52 the dislocation cross-slips into the outgoing twinned grain via redissociation into two partials. In (d) with Angelo-EAM,53 Foiles-EAM,54 and Zhou-EAM,55 the redissociated dislocation is absorbed by the CTB, with two partials gliding on the twin plane in opposite directions.56 Reproduced with permission from Xu et al., JOM 69, 814–821 (2017). Copyright 2017 The Minerals, Metals & Materials Society.

FIG. 12.

Snapshots of dislocation pileup with dominant leading screw character impinging against a CTB. Atoms are colored by adaptive common neighbor analysis: red are of hexagonal-close packed local structures, blue are BCC atoms, and all FCC atoms are deleted. In (a), incoming dislocations approach the CTB. In (b), the leading dislocation is constricted at the CTB, where two Shockley partial dislocations are recombined into a full dislocation. All five potentials produced the same results in (a) and (b). In (c), with Mishin-EAM51 and Voter-EAM potentials,52 the dislocation cross-slips into the outgoing twinned grain via redissociation into two partials. In (d) with Angelo-EAM,53 Foiles-EAM,54 and Zhou-EAM,55 the redissociated dislocation is absorbed by the CTB, with two partials gliding on the twin plane in opposite directions.56 Reproduced with permission from Xu et al., JOM 69, 814–821 (2017). Copyright 2017 The Minerals, Metals & Materials Society.

Close modal

CAC solves for approximate atomic trajectories. Collective motion of atoms, i.e., phonons, and the nucleation and propagation of defects, as well as phonon scattering by defects and other phonons, naturally emerges in the simulation. It thus provides a unified treatment for phonon transport and defect dynamics, as both are represented by the same atomic motion, without the need to assume the nature of phonon transport or scattering, or the mechanism for defect nucleation, propagation, and interactions.

Figure 13 presents CAC simulation results of a heat pulse propagating across GBs.46 A phonon representation of heat pulses composed of spatio-temporal Gaussian wave packets obeying the Bose-Einstein distribution57 is used to mimic the coherent lattice excitation by ultrashort lasers. As can be seen from Fig. 13, the simulation has captured the following:

  1. The phenomenon of phonon focusing “caustics,”58 manifested as the concentration of energy flux in the sixfold symmetric focusing pattern shown in Fig. 13(a);

  2. The Kapitza resistance, shown in Fig. 13(d) as the discontinuities in the kinetic energy at the GBs;

  3. Simultaneous ballistic and diffusive phonon transport across GBs, and a co-existence of coherent-incoherent phonon scattering by the GBs. These are clearly shown in Figs. 13(b)13(d) as both wavelike transport of the propagating heat pulse and the breaking of the 6-fold symmetry of the phonon focusing pattern due to diffusive transport and incoherent scattering by the GBs; and

  4. Local reconstruction of the GB structures, as shown in Fig. 13(e), as a result of the phonon-GB interaction during the heat pulse propagation.

FIG. 13.

Snapshots of kinetic energy distribution in space during the propagation a heat pulse across GBs. GB locations are indicated by white arrows. Reproduced with permission from Chen et al., Acta Mater. 136, 355–365 (2017). Copyright 2017 Acta Materialia Inc.

FIG. 13.

Snapshots of kinetic energy distribution in space during the propagation a heat pulse across GBs. GB locations are indicated by white arrows. Reproduced with permission from Chen et al., Acta Mater. 136, 355–365 (2017). Copyright 2017 Acta Materialia Inc.

Close modal

We have presented the CAC method for simulation of transport processes in crystalline materials. The formalism extends the Irving-Kirkwood statistical mechanical theory of transport processes for single-phase, single-component molecular systems to a concurrently coupled atomic-continuum description of polyatomic crystals. Conservation equations are derived in terms of instantaneous expressions of cell-averaged quantities using the mathematical theory of distributions.59 Fluxes are then obtained that quantify the flow of conserved quantities across the lattice cells as well as those that flow back and forth within the cells, as direct consequences of the local density definition and Newton's second law.

Basic features of the CAC formulation, in comparison with classical continuum mechanics and molecular dynamics, may be summarized as follows:

  1. As in classical continuum mechanics, the CAC governing equations are field equations, and hence can be discretized and solved using the finite element method. However, unlike the situation with continuum mechanics, each element node in CAC represents a unit cell that further contains a group of atoms; the atomic interaction and atomic-level crystal structure are thus fully built into the formulation, with no additional constitutive equations other than an interatomic potential being needed.

  2. By including the internal DOFs of each lattice cell associated with atomic movement in the formulation, CAC can resolve the atomic structure of crystalline materials and reproduce both acoustic and optical phonon modes. It also can measure fluxes at the atomic scale, such as stress near a dislocation core or heat flux across an interface or a defect.

  3. Since CAC is derived in a bottom-up manner from the atomistic, it becomes identical to the underlying atomistic model when discretized at the lattice level. CAC thus has the full applicability of MD in predicting materials behavior without the need for a priori assumption for mechanisms. It can also go beyond MD to scale up for mesoscale simulations by utilizing the continuous structure and deformation in regions where structure disorders or other discontinuities are absent or where atomic resolution is not needed.

  4. CAC is thus naturally a concurrent multiscale tool that resolves details to full atomic resolution at regions of interest, while coarse-graining to thousands of atoms per finite element elsewhere; it admits the propagation of dislocations, addresses stacking fault structures in the FE regions, and provides a unified treatment for material microstructures, defects, and phonon transport.

A comparison of representative coarse grained (CG) methods can be found in Ref. 60, and a detailed review of dynamic multiscale methods can be found in Ref. 61. CAC is currently the only concurrently coupled atomistic and continuum method whose continuum description is not based on the classical continuum mechanics equations but a reformulated field equation of conservation laws. As distinct from other CG and multiscale methods, the CAC formulation links and unifies atomistic and continuum description of transport equations. There are two unique capabilities of the CAC formulation:

  1. It can serve as a multiscale theory for concurrent atomistic-continuum simulation of polyatomic crystalline materials under a single set of governing equations, thereby expanding the current atomic-interaction-based simulation capabilities from the nanoscale to the mesoscale crystalline materials.

  2. It provides consistent formulas for calculations of continuum quantities such as stress and heat flux in atomistic or CG atomistic simulations of heterogeneous materials, noting that the majority of flux formulas in the literature or currently implemented in MD software are only applicable to homogenized, single-phase, single-component materials.

Irving and Kirkwood recognized that the classical hydrodynamics is a single-scale theory as a result of a single-scale materials description.3 Beyond the single-scale continuum mechanics or hydrodynamics, we have the generalized continuum mechanics (GCM) that advocate and employ a concurrent two-level structural description of materials.62–64 Well-established GCM theories include the Micromorphic theory,65 the Micropolar theory,66 and the Microstructural theory.67 Unlike CAC, these GCMs were formulated via a top down approach.68 Although the link between molecular dynamics and GCM has been attempted,69–72 the popular GCMs are two-level continuum-continuum formulations,60 lacking a description for the discrete structure or motion at the atomic scale.

We have also presented three sets of simulation results in this tutorial to illustrate the applicability of CAC in reproducing atomic-scale structures of interfaces, dynamic processes of dislocations and fracture, and phonon transport across GBs. These are low-temperature simulations, in which thermal fluctuations with wavelength smaller than the size of the elements in the FE regions are ignored.

For simulation of defect dynamics at finite temperature, the FE mesh size in CAC is determined by the density of defects as well as by the dominant wavelengths of phonons that interact with the defects. In regions where atomic-scale accuracy is necessary to capture interface and defect reconstruction, full atomic resolution is needed.

For low-temperature transport processes in which phonon-phonon scattering has a negligible effect, the usual FE trilinear shape functions may be replaced by those that combine the trilinear shape functions and Bloch wave functions consisting of phonons whose wavelengths are cut off by the shape function. This allows short-wavelength phonons from an atomically resolved region to traverse a coarse-grained FE region and enter another atomically resolved region73 so as to model phonon scattering by dislocations, impurities, GB or phase interfaces with all possible phonons.

For finite-temperature phonon transport when high-frequency short-wavelength phonons are dominated by thermally resistive74 or diffusive processes, the surface integrals of the kinetic momentum fluxes, i.e., fαT, in the momentum equation can be modeled as a body force with a constant mean but fluctuating in time in terms of the frequencies of the missing phonons. This is intended to reproduce the phonon density of state of the underlying atomistic model, thereby approximating the effect of the thermal fluctuations by the missing phonons. Phonon transport represented by the motion and deformation of the elements, whether it is diffusive, ballistic, or simultaneously ballistic and diffusive, will naturally emerge in the simulation as the consequence of the conservation laws, the interatomic potential, and the initial temperature.

Finally, we would like to note that, while the conservation equations and expressions of fluxes and temperature in the CAC formulation are identical to those in an atomistic model of particles, the FE solutions to the conservation equations are approximate numerical solutions. Consequently, simulated phenomena and calculated fluxes using the FE approximations feature a trade-off of efficiency and accuracy. There are thus coarse-graining errors that may increase with the FE mesh size, depending on the nature of the phenomena and the morphological length scales of the materials. Advanced algorithms that can improve the accuracy and efficiency of the simulation, such as adaptive mesh refinement and advanced multiple time step algorithm, are needed. Also, the current CAC formulation is based on the traditional interatomic potentials with well-defined site energies. Future work will be needed to extend the formulation to quantum-accurate, machine-learned potentials.

This material is based upon research supported by the U.S. Department of Energy (DOE), Basic Energy Sciences under Award No. DE-SC0006539. The contributions of D.L.M. and Y.C. are partially supported by the National Science Foundation (NSF) under Grant Nos. CMMI 1761553 and 1761512, respectively. We thank Zexi Zheng for generating the data plotted in Fig. 6.

The classical formulations of boundary value problems are based on the assumption that their solutions are smooth (i.e., with continuous derivatives) and that the equation is satisfied at every point in the domain of the problems. For problems that involve point mass, point charge, or other discreteness or discontinuities, the smoothness assumption breaks down. The mathematical tool for such problems is the theory of distributions, also known as the theory of generalized functions, developed by Laurent Schwartz in the 1950s.17 One generalized function frequently used in science and engineering is the Dirac delta-function introduced by Paul Dirac in his research in quantum mechanics and is employed in most statistical mechanics formulations to link a particle description to a field representation. The theory of distributions enables a rigorous solution to the mathematical difficulties in linking discrete and continuum descriptions that have been bypassed by ad hoc constructions or by heuristic arguments.

 Appendixes A and B present formal derivations of conservation laws based on the mathematical theory of distributions. A more detailed discussion of the integral form of conservation laws and atomistic formulas for fluxes for monatomic crystals can be found in Refs. 13–16. Due to additional quantities such as temperature discussed in this work, there are notational changes in this paper.

1. Local density definition and the balance equation of linear momentum for monatomic crystals

We define the volumetric density of mass, linear momentum, and energy over a unit cell and a time-interval as

ρ(x,t)kmkδ¯V(xrk),ρ(x,t)v(x,t)kmkvkδ¯V(xrk),ρ(x,t)e(x,t)k(12mkvk2+Φk)δ¯V(xrk),
(A1)

where δ¯V is the time step-averaged box function δV defined in (8) and is used to serve the role of averaging the IK point-function local density over a unit cell V and a time-interval Δt.

The derivative of δV does not exist everywhere as a classical function, but it exists as a distribution. A distribution f(x) is not required to have a value at any point x, it is defined by its action on a test function φ(x) that infinitely many times differentiable and has a bounded support. For example, the Dirac delta-function δ(x) is defined as (δ, φ) = φ(0). In this sense, the derivative of δV is the surface distribution that acts on a test function φ as (tδV(xrk,φ(x,t))=VφvkndS. Any smooth function g(x) can be viewed as a distribution whose action on a test function φ(x) is defined as the integral of the product g(x) φ(x) over the whole range of the variable x. Furthermore, any distribution f(x) is proved to be the limit of a sequence of smooth functions fɛ in the sense that the numerical sequence (fɛ, φ) converges to the number (f, φ) for any test function as the regularization parameter ɛ approaches 0. In this sense, the derivative of δV can formally be written as a definite integral of a δ-function, that is, as the limit of the corresponding integral in which the δ-function is replaced by a regularized distribution δɛ that is a smooth function and converges to δ as ε0 in the distributional sense,59 i.e.,

tδV(xrk)=limε01VVvknδε(x+xrk)dS1VVvknδ(x+xrk)dS.
(A2)

For the rest of the appendixes, any formal integral of a δ-function, like the one in the right-hand side of [Eq. (A2)], is understood as a distribution defined by such a limit. Any relation between distributions obtained in this way is proved to be independent of the choice of regularization.59 

The internal force density is also a volume density. An important fact is that the total internal force in a volume element V is equal to the total interaction force transmitted across the enclosing surface ∂V, which can be described as a line-plane intersection problem. We denote an oriented surface element that is centered at x with normal n and area A as An; the intersection of the line segment Lkl and surface element An can be expressed as

LklδAn(xφ)dφ1AAnLklδ(x+xφ)dφndS=1A{1ifLklintersectsAnandnLkl0,1ifLklintersectsAnandnLkl<0,0ifLkldoesnotintersectAn.
(A3)

It is noted that the time-interval average of LklδAn(xφ)dφ over a time step is itself.

The crossing of an atom through a surface element during a time step Δt can also be expressed as the intersection of the path of the atom rk and a surface element An with area A and normal n as

δ¯An(xrk)v~k1AΔtAn0Δtδ(xrk)v~kndτdS=1AΔt{1ifrkintersectsAninΔtandnv~k0,1ifrkintersectsAninΔtandnv~k<0,0ifrkldoesnotintersectAn.
(A4)

Expressing the total internal force in the volume element V as the total interaction forces crossing the bounding surface ∂V, we obtain the internal force density in terms of site energy Φk, according to Eq. (A3), as

kFkδV(xrk)=k,lΦkrlδV(xrk)=k,lk(ΦkrlΦlrk)δV(xrk)=1VVkV,lV(ΦkrlΦlrk)Lkldφδ(x+xφ)dS=1VVk,l(k)ΦkrlLkldφδ(x+xφ)dS.
(A5)

The time rate of change of the linear momentum within V can be derived, as a result of Newton's second law applying to the particles in V, using Eqs. (A5) and (A2), as

tkmkvkδ¯V(xrk)=kmkv˙kδ¯V(xrk)+kmkvktδ¯V(xrk)=kFkδ¯V(xrk)1Δt0Δtdτ1VVkmkvkvkδ(x+xrk)ndS=1Δt0Δtdτ1VVk,l(k)ΦkrlLkldφδ(x+xφ)dS1Δt0Δtdτ1VV{kmkv~kv~kδ(x+xrk)+ρvv}ndS=Vtd2xVρvvndS,
(A6)

where the stress vector t on the surface element of normal n is expressed as a line-plane intersection problem, cf. Eqs. (A3) and (A4), as

t(x,n)=k,l(k)ΦkrlLklδ¯An(φx)dφkmkv~kv~kδ¯An(xrk)=k<lFklLklδ¯An(φx)dφkmkv~kv~kδ¯An(xrk),
(A7)

where Fkl(Φl/rkΦk/rl) is a measure of the strength of the interaction between atom k and atom l.

2. Conservation equation of energy for monatomic crystals

Similarly, the time rate of change of the energy density can also be expressed in terms of atomic variables as

(ρe)t=tkEkδ¯V(xrk)=kE˙kδ¯V(xrk)+kEktδ¯V(xrk).
(A8)

The first term in the right-hand side of Eq. (A8) can be expressed as

kE˙kδ¯V(xrk)=kddt(12mkvk2+Φk)δ¯V(xrk)=k(mkvkv˙k+Φ˙k)δ¯V(xrk)=k(Fkvk+lΦkrlr˙l)δ¯V(xrk)=k,l(Φlrkvk+Φkrlvl)δ¯V(xrk)=1Δt0Δtdτ1VVk,l(k)ΦkrlvlLkldφδ(x+xφ)dS,
(A9)

and the second term, according to Eq. (A2), as

kEktδ¯V(xrk)=1Δt0ΔtdτVk1VEkvkδ(xxrk)ndS=1Δt0ΔtdτVk1V(12mkv~k2+mkv~kv+12mkv2+Φk)v~kδ(x+xrk)ndS1VVρevndS.
(A10)

Substituting Eqs. (A9) and (A10) into Eq. (A8), we obtain the integral form of the energy conservation law as

(ρe)t=V(qpot+σpotv)ndS+V(qkin+σkinvρev)ndS
(A11)

and the heat flux vector in the indicial notation as

qi(x,t)=k(12mkv~k2+Φk)v~kδ¯Ai(xrk)+k,l(k)Φkrlv~lLklδ¯Ai(φx)dφ.
(A12)

1. Balance equation of linear momentum for polyatomic crystals

The rate of change of the density of linear momentum defined in Eq. (11) can be expressed as

(ραvα)t=k=1Nlξ=1Namξv˙kξδ¯V(rkx)δ¯Vα(Δrkξy)+k=1Nlξ=1Namξvkξtδ¯V(rkx)δ¯Vα(Δrkξy)+k=1Nlξ=1Namξvkξδ¯V(rkx)tδ¯Vα(Δrkξy).
(B1)

The first term in the right-hand side of Eq. (B1) can be derived according to Newton's second law and, in the absence of a body force, it is equal to the internal force density fαint(x). Note that the translational symmetry of potential energies requires that

Φkξrkξ=l(k)NlηNaΦkξrlηη(ξ)NaΦkξrkη.
(B2)

This leads to the expression of the total force acting on atom due to all the other atoms as

Φrkξ=lNlηNaΦlηrkξ=l(k)NlηNaΦlηrkξ+ηξNaΦkηrkξ+Φkξrkξ=l(k)NlηNaΦlηrkξ+ηξNaΦkηrkξl(k)NlηNaΦkξrlηη(ξ)NaΦkξrkη=l(k)NlηNa(ΦlηrkξΦkξrlη)+ηξNa(ΦkηrkξΦkξrkη).
(B3)

Figure 14 shows that the total interaction force on atom α in the unit cell at point x can be decomposed into the interaction forces on atom α due to the interactions between atoms in different unit cells and that between atoms in the same unit cell. This then leads to the definitions of surface tractions, tαpot and ταpot. The internal force density fαint(x) can then be written, according to Eq. (B3), using the line-plane intersection defined in Eq. (A3), as

fαint(x)=1VVkV,lVNlηNa(ΦlηrkαΦkαrlη)rlηrkαδ(x+xφ)dφdS+VαkNlδ¯V(xrk)ηαNa(ΦkηrkαΦkαrkη)rkηrkαδ(y+yφ)dφdSα1VVtαpotdS+1VVαταpotdSα,
(B4)

where

tαpot(x,n)=kV,lVNlηNa(ΦlηrkαΦkαrlη)rlηrkαδ¯An(xφ)dφ,
(B5)
ταpot(x,y,n)=kNlVδ¯V(xrk)ηαNa(ΦkηrkαΦkαrkη)rkηrkαδ¯An(yφ)dφ=kVNlηαNa(ΦkηrkαΦkαrkη)rkηrkαδ¯An(yφ)dφ.
(B6)
FIG. 14.

(a) The volume element V of the unit cell at point x; (b) two contributions to fαint(x): (1) the interaction forces that cross ∂V and (2) the interaction forces that only cross Vα; (c) the volume element Vα.

FIG. 14.

(a) The volume element V of the unit cell at point x; (b) two contributions to fαint(x): (1) the interaction forces that cross ∂V and (2) the interaction forces that only cross Vα; (c) the volume element Vα.

Close modal

The second term in Eq. (B1) can be expressed according to Eq. (A2), noting that all the cross terms involving the velocity differences vanish over a time-interval average and α=1NαδVα(yΔrkξ)=1, as

k=1Nlξ=1Namξvkξtδ¯V(xrk)δ¯Vα(yΔrkξ)=1VVk=1Nlξ=1Namξvkξvknδ¯(x+xrk)δ¯Vα(yΔrkξ)dS=1VVk=1Nlξ=1Namξ(v~k+Δv~kξ+vα)(v~k+v)nδ¯(x+xrk)δ¯Vα(yΔrkξ)dS=1VVk=1Nlmαv~kv~knδ¯(x+xrk)dS1VVραvαvndS1VV[tαkinραvαvn]dS
(B7)

and the third term as

k=1Nlξ=1Namξvkξδ¯V(xrk)tδ¯Vα(yΔrkξ)=Vαk=1NlδV(xrk)ξ=1NamξvkξΔvkξδ(y+yΔrkξ)ndSα=k=1NlδV(xrk)Vαξ=1NamξΔv~kξΔv~kξδ(y+yΔrkξ)ndSα+1VVαρvαΔvαndSα1VVα[ταkinρvαΔvαn]dSα,
(B8)

with

tαkin(x,n)=k=1Nlmαv~kv~kδ¯An(xrk),
(B9)
ταkin(x,y,n)=kNlVδ¯V(xrk)mαΔv~kαΔv~kαδ¯An(yΔrkα)=kVNlmαΔv~kαΔv~kαδ¯An(yΔrkα).
(B10)

Substituting Eqs. (B4), (B7), and (B8) into Eq. (B11), we obtain the conservation equation of linear momentum as

(ραvα)t=1VV(tαραvαvn)dS+1VVα(ταραvαΔvαn)dSα,
(B11)

where

tα=tαpot(x,n)+tαkin(x,n),andτα=ταpot(x,y,n)+ταkin(x,y,n).
(B12)

The total stress vector on a surface element at point y in the unit cell located at x due to the interaction and motion of all the atoms is then expressed as a sum of potential and kinetic parts as

t(x,y,n)=α=1Na(tα+τα)=α=1Na(tαpot+tαkin+ταpot+ταkin)=tpot(x,y,n)+tkin(x,y,n),
(B13)

where

tpot(x,y,n)=α=1Na(tαpot+ταpot)=kV,lVNlα,ηNa(ΦlηrkαΦkαrlη)rlηrkαδ¯An(xφ)dφ+kVNlα,ηαNa(ΦkηrkαΦkαrkη)rkηrkαδ¯An(yφ)dφ=k,l(k)Nlα,ηNarlηrkαΦkαrlηδ¯An(xφ)dφ+kNlα,ηαNaVδ¯V(xrk)rkηrkαΦkαrkηδ¯An(yφ)dφ
(B14)

and

tkin(x,y,n)=α=1Na(tαkin+ταkin)=α=1Nak=1Nlmαv~kv~kδ¯An(xrk)α=1Nak=1NlmαΔv~kαΔv~kαVδ¯V(xrk)δ¯An(yΔrkα).
(B15)

The stress tensor can then be expressed in terms of the stress vector as

σij(x,y)=ti(x,y,ej)orσαij(x,y)=tαi(x,y,ej),i,j=1,2,3.
(B16)

2. Conservation equation of energy for polyatomic crystals

The rate of change of the energy density can be written as

(ραeα)t=k=1Nlξ=1NaE˙kξδ¯V(xrk)δ¯Vα(yΔrkξ)+k=1Nlξ=1NaEkξtδ¯V(xrk)δ¯Vα(yΔrkξ)+k=1Nlξ=1NaEkξδ¯V(xrk)tδ¯Vα(yΔrkξ).
(B17)

With

E˙kξ=mkξvkξv˙kξ+Φ˙kξ=lNlηNa(ΦlηrkξvkξΦkξrlηvlη)=l(k)NlηNaΦlηrkξvkξηξNaΦkηrkξvkξΦkξrkξvkξ+l(k)NlηNaΦkξrlηvlη+ηξNaΦkξrkηvkη+Φkξrkξvkξ=l(k)NlηNa(ΦlηrkξvkξΦkξrlηvlη)ηξNa(ΦkηrkξvkξΦkξrkηvkη),
(B18)

the first term in the right-hand side of Eq. (B17) can be expressed, using Eq. (B18) and the line-plane intersection theorem, as

k=1Nlξ=1NaE˙kξδ¯V(xrk)δ¯Vα(yΔrkξ)=1VVkV,lVNlηNa(ΦlηrkαvkαΦkαrlηvlη)rlηrkαδ(x+xφ)dφdS+kNlδ¯V(xrk)VαηαNa(ΦkηrkαvkαΦkαrkηvkη)rkηrkαδ(y+yφ)dφdSα.
(B19)

The second term can be derived according to the distributional derivative of the box function expressed in Eq. (A2) as

k=1Nlξ=1NaEkξtδ¯V(xrk)δ¯Vα(yΔrkξ)=1VVk=1Nlξ=1NaEkξvknδ(x+xrk)δ¯Vα(yΔrkξ)dS=1VVk=1Nl(12mαv~kα2+v~kαvα+(vα)2+Φkα)v~knδ(x+xrk)dSVραeαvndS=1VVk=1Nl(12mαv~kα2+Φkα)v~kδ(x+xrk)ndS1V1Δt0ΔtdtVk=1Nlv~kαvαv~kδ¯(x+xrk)ndSVραeαvndS1VV(qαkin+tαkinvαραeαv)ndS
(B20)

and the third term as

k=1Nlξ=1NaEkξδ¯V(xrk)tδ¯Vα(yΔrkξ)=k=1Nlδ¯V(xrk)Vαξ=1NaEkξΔvkξnδ(y+yΔrkξ)dSα=k=1Nlδ¯V(xrk)Vαξ=1Na(12mξv~kξ2+Φkξ+v~kξvα+(vα)2)Δv~kξnδ(yΔrkξ)dSαVVαραeαΔvαndS1VVα(jαkin+ταkinvαραeαΔvα)ndSα.
(B21)

This leads to

(ραeα)t=1VV(qα+tαvαραeαv)ndS+1VVα(jα+ταvαραeαΔvα)ndSα.
(B22)

The total heat flux is thus expressed as

q(x,y)=α=1Na(qα+jα)=qiei,
(B23)

with

qi(x,y)=k,lkNlα,ηNaΦkαrlηv~kαrlηrkαδ¯Ai(xφ)dφk=1Nlα=1Nα(12mαv~kα2+Φkα)v~kδ¯Ai(xrk)
(B24)

and

ji(x,y)=kNlVδ¯V(xrk)α,ηαNaΦkαrkηv~kαrkηrkαδ¯Ai(yφ)dφk=1NlVδ¯V(rkx)αNa(12mαv~kα2+Φkα)×Δv~kαδ¯Ai(yΔrkα).
(B25)

Irving and Kirkwood noted that a definition for the site energies of particles is required for definition of energy density as well as for the derivation of a field representation of energy conservation law in terms of particle variables.3 The conservation equations and the expressions for fluxes derived in the appendixes are valid for any analytical many-body potentials with well-defined site energies. A detailed comparison and discussion of the fluxes defined as a surface density presented in  Appendix A and the fluxes formulas expressed as a volume density can be found in Ref. 13. Specifically, analytical and MD simulation results presented in Ref. 13 show that replacing the Dirac δ with a finite-sized volume weighting function changes the fundamental nature of fluxes as a surface density. Consequently, the dynamic balance between the rate of change of the total momentum or energy in a volume element and the fluxes across the surface boundary cannot be established, leading to the failure of the volume averaged flux formulas in satisfying the momentum and energy conservation laws as well as typical transport boundary conditions. The equivalence between flux formulas expressed as a pair density derived using the kinetic theory75,76 and that expressed as a single density derived using statistical mechanics,3 as ensemble averaged point functions or as surface-averages, can be found in Ref. 16.

1.
C.
Kittel
,
Introduction to Solid State Physics
(
John Wiley & Sons, Inc.
,
New York
,
1956
).
2.
M.
Born
and
K.
Huang
,
Dynamical Theory of Crystal Lattices
(
Clarendon Press
,
Oxford
,
1954
).
3.
J.
Irving
and
J. G.
Kirkwood
, “
The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics
,”
J. Chem. Phys.
18
,
817
829
(
1950
).
4.
P.
Dirac
,
Principles of Quantum Mechanics
(Clarendon Press, Oxford,
1930
).
5.
B.
Todd
,
P. J.
Daivis
, and
D. J.
Evans
, “
Heat flux vector in highly inhomogeneous nonequilibrium fluids
,”
Phys. Rev. E
51
,
4362
(
1995
).
6.
B. D.
Todd
,
D. J.
Evans
, and
P. J.
Daivis
, “
Pressure tensor for inhomogeneous fluids
,”
Phys. Rev. E
52
,
1627
1638
(
1995
).
7.
Y.
Chen
, “
Local stress and heat flux in atomistic systems involving three-body forces
,”
J. Chem. Phys.
124
,
054113
(
2006
).
8.
E.
Smith
,
D.
Heyes
,
D.
Dini
, and
T.
Zaki
, “
Control-volume representation of molecular dynamics
,”
Phys. Rev. E
85
,
056705
(
2012
).
9.
J. G.
Kirkwood
, “
The statistical mechanical theory of transport processes. I. General theory
,”
J. Chem. Phys.
14
,
180
201
(
1946
).
10.
Y.
Chen
, “
Reformulation of microscopic balance equations for multiscale materials modeling
,”
J. Chem. Phys.
130
,
134706
(
2009
).
11.
Y.
Chen
and
J.
Lee
, “
Atomistic formulation of a multiscale field theory for nano/micro solids
,”
Philos. Mag.
85
,
4095
4126
(
2005
).
12.
D. J.
Evans
and
G.
Morriss
,
Statistical Mechanics of Nonequilibrium Liquids
(
Cambridge University Press
,
Cambridge
,
2008
).
13.
Y.
Chen
and
A.
Diaz
, “
Physical foundation and consistent formulation of atomic-level fluxes in transport processes
,”
Phys. Rev. E
98
,
052113
(
2018
).
14.
Y.
Chen
and
A.
Diaz
, “
Local momentum and heat fluxes in transient transport processes and inhomogeneous systems
,”
Phys. Rev. E
94
,
053309
(
2016
).
15.
Y.
Chen
, “
The origin of the distinction between microscopic formulas for stress and Cauchy stress
,”
Europhys. Lett.
116
,
34003
(
2016
).
16.
A.
Diaz
,
D.
Davydov
, and
Y.
Chen
, “
On the equivalence of the two foundational formulations for atomistic flux in inhomogeneous transport processes
,”
Proc. R. Soc. A
475
,
20180688
(
2019
).
17.
L.
Schwartz
,
Théorie des Distributions
(
Hermann
,
Paris
,
1950
).
18.
W. G.
Hoover
,
Computational Statistical Mechanics
(
Elsevier
,
New York
,
1991
).
19.
S.
Root
,
R. J.
Hardy
, and
D. R.
Swanson
, “
Continuum predictions from molecular dynamics simulations: Shock waves
,”
J. Chem. Phys.
118
,
3161
3165
(
2003
).
20.
J.
Turney
,
A.
McGaughey
, and
C.
Amon
, “
Assessing the applicability of quantum corrections to classical thermal conductivity predictions
,”
Phys. Rev. B
79
,
224305
(
2009
).
21.
Y.
Li
 et al, “
Phonon spectrum and phonon focusing in coarse-grained atomistic simulations
,”
Comput. Mater. Sci.
162
,
21
32
(
2019
).
22.
A.
Cantarero
and
F. X.
Àlvarez
,
Nanoscale Thermoelectrics
(
Springer
,
New York
,
2014
), pp.
1
39
.
23.
R.
Kubo
, “
The fluctuation-dissipation theorem
,”
Rep. Prog. Phys.
29
,
255
(
1966
).
24.
J.
Rigelesaiyin
,
A.
Diaz
,
W.
Li
,
L.
Xiong
, and
Y.
Chen
, “
Asymmetry of the atomic-level stress tensor in homogeneous and inhomogeneous materials
,”
Proc. R. Soc. A Math. Phys. Eng. Sci.
474
,
20180155
(
2018
).
25.
S.
Xu
 et al, “
PyCAC: The concurrent atomistic-continuum simulation environment
,”
J. Mater. Res.
33
,
857
871
(
2018
).
26.
Q.
Deng
,
L.
Xiong
, and
Y.
Chen
, “
Coarse-graining atomistic dynamics of brittle fracture by finite element method
,”
Int. J. Plast.
26
,
1402
1414
(
2010
).
27.
Q.
Deng
and
Y.
Chen
, “
A coarse-grained atomistic method for 3D dynamic fracture simulation
,”
J. Multiscale Comput. Eng.
11
,
227
237
(
2013
).
28.
L.
Xiong
and
Y.
Chen
, “
Coarse-grained simulations of single-crystal silicon
,”
Model. Simul. Mater. Sci. Eng.
17
,
035002
(
2009
).
29.
L.
Xiong
,
Q.
Deng
,
G.
Tucker
,
D. L.
McDowell
, and
Y.
Chen
, “
A concurrent scheme for passing dislocations from atomistic to continuum domains
,”
Acta Mater.
60
,
899
913
(
2012
).
30.
L.
Xiong
,
G.
Tucker
,
D. L.
McDowell
, and
Y.
Chen
, “
Coarse-grained atomistic simulation of dislocations
,”
J. Mech. Phys. Solids
59
,
160
177
(
2011
).
31.
S.
Yang
,
L.
Xiong
,
Q.
Deng
, and
Y.
Chen
, “
Concurrent atomistic and continuum simulation of strontium titanate
,”
Acta Mater.
61
,
89
102
(
2013
).
32.
L.
Xiong
 et al, “
Coarse-grained elastodynamics of fast moving dislocations
,”
Acta Mater.
104
,
143
155
(
2016
).
33.
L.
Xiong
,
Q.
Deng
,
G. J.
Tucker
,
D. L.
McDowell
, and
Y.
Chen
, “
Coarse-grained atomistic simulations of dislocations in Al, Ni and Cu crystals
,”
Int. J. Plast.
38
,
86
101
(
2012
).
34.
L.
Xiong
,
D. L.
McDowell
, and
Y.
Chen
, “
Nucleation and growth of dislocation loops in Cu, Al and Si by a concurrent atomistic-continuum method
,”
Scr. Mater.
67
,
633
636
(
2012
).
35.
S.
Xu
,
L.
Xiong
,
Y.
Chen
, and
D. L.
McDowell
, “
Edge dislocations bowing out from a row of collinear obstacles in Al
,”
Scr. Mater.
123
,
135
139
(
2016
).
36.
S.
Xu
,
L.
Xiong
,
Y.
Chen
, and
D. L.
McDowell
, “
An analysis of key characteristics of the Frank-Read source process in FCC metals
,”
J. Mech. Phys. Solids
96
,
460
476
(
2016
).
37.
S.
Xu
,
L.
Xiong
,
Y.
Chen
, and
D. L.
McDowell
, “
Sequential slip transfer of mixed-character dislocations across Σ3 coherent twin boundary in FCC metals: A concurrent atomistic-continuum study
,”
npj Comput. Mater.
2
,
15016
(
2016
).
38.
L.
Xiong
,
S.
Xu
,
D. L.
McDowell
, and
Y.
Chen
, “
Concurrent atomistic–continuum simulations of dislocation–void interactions in fcc crystals
,”
Int. J. Plast.
65
,
33
42
(
2015
).
39.
S.
Xu
,
L.
Xiong
,
Y.
Chen
, and
D. L.
McDowell
, “
Shear stress- and line length-dependent screw dislocation cross-slip in FCC Ni
,”
Acta Mater.
122
,
412
419
(
2017
).
40.
S.
Xu
,
L.
Xiong
,
Y.
Chen
, and
D.
McDowell
, “
Validation of the concurrent atomistic-continuum method on screw dislocation/stacking fault interactions
,”
Crystals
7
,
120
(
2017
).
41.
S.
Yang
and
Y.
Chen
, “Concurrent atomistic-continuum simulation of defects in polyatomic ionic materials,” in Multiscale Materials Modeling for Nanomechanics, edited by C. Weinberger and G. Tucker (Springer, Cham, Switzerland) Vol. 245, pp. 261–296.
42.
S.
Yang
and
Y.
Chen
, “
Concurrent atomistic and continuum simulation of bi-crystal strontium titanate with tilt grain boundary
,”
Proc. R. Soc. A Math. Phys. Eng. Sci.
471
,
20140758
(
2015
).
43.
S.
Yang
,
N.
Zhang
, and
Y.
Chen
, “
Concurrent atomistic–continuum simulation of polycrystalline strontium titanate
,”
Philos. Mag.
95
,
2697
2716
(
2015
).
44.
X.
Chen
,
L.
Xiong
,
D. L.
McDowell
, and
Y.
Chen
, “
Effects of phonons on mobility of dislocations and dislocation arrays
,”
Scr. Mater.
137
,
22
26
(
2017
).
45.
L.
Xiong
,
D. L.
McDowell
, and
Y.
Chen
, “
Sub-THz phonon drag on dislocations by coarse-grained atomistic simulations
,”
Int. J. Plast.
55
,
268
278
(
2014
).
46.
X.
Chen
 et al, “
Ballistic-diffusive phonon heat transport across grain boundaries
,”
Acta Mater.
136
,
355
365
(
2017
).
47.
H.-S.
Lee
 et al, “
Defect energetics in SrTiO3 symmetric tilt grain boundaries
,”
Phys. Rev. B
83
,
104110
(
2011
).
48.
G.
Springholz
and
K.
Wiesauer
, “
Nanoscale dislocation patterning in PbTe/PbSe (001) lattice-mismatched heteroepitaxy
,”
Phys. Rev. Lett.
88
,
015507
(
2001
).
49.
Y.
Li
,
Z.
Fan
,
W.
Li
,
D. L.
McDowell
, and
Y.
Chen
, “
A multiscale study of misfit dislocations in PbTe/PbSe(001) heteroepitaxy
,”
J. Mater. Res.
,
34
,
2306
2314
(
2019
).
50.
Q.
Deng
, “Coarse-graining atomistic dynamics of fracture by finite element method formulation, parallelization and applications,” Ph.D. dissertation (
University of Florida
,
2011
).
51.
Y.
Mishin
,
D.
Farkas
,
M.
Mehl
, and
D.
Papaconstantopoulos
, “
Interatomic potentials for monoatomic metals from experimental data and ab initio calculations
,”
Phys. Rev. B
59
,
3393
(
1999
).
52.
A. F.
Voter
and
S. P.
Chen
, “
Accurate interatomic potentials for Ni, Al and Ni3Al
.”
MRS Online Proc. Libr. Arch.
82,
175 (
1986
).
53.
J. E.
Angelo
,
N. R.
Moody
, and
M. I.
Baskes
, “
Trapping of hydrogen to lattice defects in nickel
,”
Model. Simul. Mater. Sci. Eng.
3
,
289
(
1995
).
54.
S. M.
Foiles
and
J.
Hoyt
, “
Computation of grain boundary stiffness and mobility from boundary fluctuations
,”
Acta Mater.
54
,
3351
3357
(
2006
).
55.
X.
Zhou
,
R.
Johnson
, and
H.
Wadley
, “
Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers
,”
Phys. Rev. B
69
,
144113
(
2004
).
56.
S.
Xu
,
L.
Xiong
,
Y.
Chen
, and
D. L.
McDowell
, “
Comparing EAM potentials to model slip transfer of sequential mixed character dislocations across Two symmetric tilt grain boundaries in Ni
,”
JOM
69
,
814
821
(
2017
).
57.
X.
Chen
,
A.
Chernatynskiy
,
L.
Xiong
, and
Y.
Chen
, “
A coherent phonon pulse model for transient phonon thermal transport
,”
Comput. Phys. Commun.
195
,
112
116
(
2015
).
58.
J. P.
Wolfe
,
Imaging Phonons: Acoustic Wave Propagation in Solids
(
Cambridge University Press
,
Cambridge
,
1998
).
59.
V. S.
Vladimirov
,
A.
Jeffrey
, and
F. E.
Schroeck
, “
Equations of mathematical physics
,”
Am. J. Phys.
39
,
1548
1548
(
1971
).
60.
Y.
Chen
,
J.
Zimmerman
,
A.
Krivtsov
, and
D.
McDowell
, “
Assessment of atomistic coarse-graining methods
,”
Int. J. Eng. Sci.
49
,
1337
1349
(
2011
).
61.
A.
Diaz
,
D.
McDowell
, and
Y.
Chen
,
Generalized Models and Non-Classical Approaches in Complex Materials 2
(
Springer
,
New York
,
2018
), pp.
55
77
.
62.
A. C.
Eringen
,
Microcontinuum Field Theories: Foundations and Solids
(
Springer
,
New York
,
1999
), Vol. 487.
63.
Y.
Chen
,
J. D.
Lee
, and
A.
Eskandarian
, “
Atomistic viewpoint of the applicability of microcontinuum theories
,”
Int. J. Solids Struct.
41
,
2085
2097
(
2004
).
64.
Y.
Chen
,
J. D.
Lee
, and
A.
Eskandarian
, “
Examining the physical foundation of continuum theories from the viewpoint of phonon dispersion relation
,”
Int. J. Eng. Sci.
41
,
61
83
(
2003
).
65.
A. C.
Eringen
,
Mechanics of Micromorphic Continua
(
Springer
,
New York
,
1968
).
66.
A. C.
Eringen
,
Microcontinuum Field Theories
(
Springer
,
New York
,
1999
), pp.
101
248
.
67.
R. D.
Mindlin
, “
Micro-structure in linear elasticity
,”
Arch. Rational Mech. Anal.
16
,
51
78
(
1964
).
68.
J. D.
Lee
and
Y.
Chen
, “
Constitutive relations of micromorphic thermoplasticity
,”
Int. J. Eng. Sci.
41
,
387
399
(
2003
).
69.
Y.
Chen
,
J.
Lee
, and
A.
Eskandarian
, “
Atomistic counterpart of micromorphic theory
,”
Acta Mech.
161
,
81
102
(
2003
).
70.
Y.
Chen
and
J. D.
Lee
, “
Connecting molecular dynamics to micromorphic theory. (I). Instantaneous and averaged mechanical variables
,”
Phys. A Stat. Mech. Appl.
322
,
359
376
(
2003
).
71.
Y.
Chen
and
J. D.
Lee
, “
Connecting molecular dynamics to micromorphic theory. (II). Balance laws
,”
Phys. A Stat. Mech. Appl.
322
,
377
392
(
2003
).
72.
Y.
Chen
and
J. D.
Lee
, “
Determining material constants in micromorphic theory through phonon dispersion relations
,”
Int. J. Eng. Sci.
41
,
871
886
(
2003
).
73.
X.
Chen
,
A.
Diaz
,
L.
Xiong
,
D. L.
McDowell
, and
Y.
Chen
, “
Passing waves from atomistic to continuum
,”
J. Comput. Phys.
354
,
393
402
(
2018
).
74.
B. H.
Armstrong
, “
N processes, the relaxation-time approximation, and lattice thermal conductivity
,”
Phys. Rev. B
32
,
3381
(
1985
).
75.
J. G.
Kirkwood
,
F. P.
Buff
, and
M. S.
Green
, “
The statistical mechanical theory of transport processes. III. The coefficients of shear and bulk viscosity of liquids
,”
J. Chem. Phys.
17
,
988
994
(
1949
).
76.
M.
Born
and
H.
Green
, “
A general kinetic theory of liquids. III. Dynamical properties
,”
Proc. R. Soc. Lond. A
190
,
455
474
(
1947
).