Material frame representation of equivalent stress tensor for discrete solids Текст научной статьи по специальности «Физика»

УДК 539.2 + 539.3

Материальное представление эквивалентного тензора напряжений

для дискретных твердых тел

В.А. Кузькин1,2, A.M. Кривцов1,2, R.E. Jones3, J.A. Zimmerman3

1 Институт проблем машиноведения РАН, Санкт-Петербург, 199178, Россия

2 Санкт-Петербургский государственный политехнический университет, Санкт-Петербург, 195251, Россия

3 Сандийская национальная лаборатория, Ливермор, СА 94551-0969, США

В работе получены1 выражения для эквивалентных тензоров напряжений Коши и Пиолы для дискретных твердых тел. В случае однородной деформации данные выражения являются точными. Для вывода использованы материальное описание, длинноволновое приближение и разделение движений на континуальные и тепловые. Полученные выражения для эквивалентных тензоров напряжений применимы в случае произвольных многочастичных взаимодействий. Рассмотрено несколько примеров вычисления напряжений в дискретных твердых телах с парными и многочастичными взаимодействиями. Проведено сравнение с известными выражениями для тензора напряжений в дискретной среде (Люси, Харди и Хайнца-Пола-Биндера) и континуальным определением напряжения. Показано, что в случае однородной деформации и конечной температуры напряжения, рассчитанные на основе формулы, предложенной в данной работе, в точности совпадают с напряжениями, получаемыми на основе классического континуального определения. Тензоры напряжений Люси и Харди дают аналогичные результаты только в случае осреднения по достаточно большому объему. Таким образом, при выгшслении напряжений в дискретныж твердыж телах выражение, предложенное в данной работе, может быть использовано в качестве эталона.

Ключевые слова: тензор напряжений Коши, материальное описание, молекулярная динамика

Material frame representation of equivalent stress tensor for discrete solids

V.A. Kuzkin1,2, A.M. Krivtsov1,2, R.E. Jones3, and J.A. Zimmerman3

1 Institute for Problems in Mechanical Engineering RAS, St. Petersburg, 199178, Russia

2 Saint Petersburg State Polytechnical University, St. Petersburg, 195251, Russia

3 Sandia National Laboratories, Livermore, СА 94551-0969, USA

In this paper, we derive expressions for equivalent Cauchy and Piola stress tensors that can be applied to discrete solids and are exact for the case of homogeneous deformation. The main principles used for this derivation are material frame formulation, long wave approximation and decomposition of particle motion into continuum and thermal parts. Equivalent Cauchy and Piola stress tensors for discrete solids are expressed in terms of averaged interparticle distances and forces. No assumptions about interparticle forces are used in the derivation, thereby ensuring our expressions are valid irrespective of the choice of interatomic potential used to model the discrete solid. The derived expressions are used for calculation of the local Cauchy stress in several test problems. The results are compared with prediction of the classical continuum definition (force per unit area) as well as existing discrete formulations (Hardy, Lucy, and HeinzPaul-Binder stress tensors). It is shown that in the case of homogeneous deformations and finite temperatures the proposed expression leads to exactly the same values of stresses as classical continuum definition. Hardy and Lucy stress tensors give the same result only if the stress is averaged over a sufficiently large volume. Thus, given the lack of sensitivity to averaging volume size, the derived expressions can be used as benchmarks for calculation of stresses in discrete solids.

Keywords: Cauchy stress tensor, material frame formulation, molecular dynamics

1. Introduction

Many phenomena in materials science and engineering at different length scales can be considered from both discrete [1-3] and continuum [4, 5] points of view. In particular, joint application of discrete and continuum techniques

is important for investigation of novel nanomaterials [6, 7], such as graphene [8, 9]. However, given the differences between discrete and continuum descriptions of matter, comparison and coupling [10] of the results obtained using these two approaches is difficult. The comparison requires

© Кузькин B.A., Кривцов A.M., Джонс P., Циммерман Д., 2014

the estimation of equivalent continuum fields, such as the stress tensor, for discrete systems. Several approaches for calculation of continuum fields from discrete systems have been proposed in the literature. One of the first expressions for a stress tensor derived for a system of particles was derived by Clausius [11]. In equilibrium, this virial stress is equivalent to an average Cauchy stress in some finite domain containing particles. Calculation of highly localized pointwise stresses is a more challenging problem. Statistical physics, and a particular averaging kernel — the Dirac delta, was employed by Irving and Kirkwood [12] to calculate stresses that are consistent with the macroscopic balance laws. Though the obtained expression for the stress tensor has great importance from theoretical point of view, it does not result in a continuous field. Additional modifications are required in order to calculate local stresses using Irving-Kirkwood approach. In particular, Irving-Kirk-wood expression has been adjusted for calculation of local stresses and used in molecular dynamics simulations in the paper [13].

Other approaches for estimating a local continuum field include the method [14], which uses Fourier transformations, and the one by Hardy [15-19], which uses kernels with finite support suitable for computer simulations. With regard to Hardy's approach, the dependence of stress on the form of localization functions and the radius of localization volume has been investigated in the paper [18], and the generalization of Hardy's approach for calculation of stress in the material formulation has been carried out [19]. A simpler formalism proposed by Lucy [20] has been used by Hoover [21] for calculation of stresses in shock waves. This approach avoids integration of the kernel along the bond and is therefore computationally more efficient.

The Lucy and Hardy formalisms are based on the spatial frame formulation, i.e. continuum properties are computed at the points fixed in space. These and other spatial frame formulations are commonly used in fluid dynamics. In solid mechanics theories, such as thermoelasticity, plasticity, and fracture mechanics, the material frame formulation is more widespread. In the material frame formulation, stresses are calculated at material points. This approach is more efficient for dealing with problems featuring complexities such as free surfaces, interfaces, cracks, and inclusions. Accurate formulation of equivalent continuum parameters for discrete systems in the material frame is important, for example, for development of hybrid discrete-continuum methods for solution of solid mechanics problems [19, 22].

The choice between spatial and material frame formulations strongly influences the expressions for equivalent continuum properties of discrete systems. This difference has been first emphasized by Hoover [2]. It has been shown that the virial theorem for solids has two different formulations. In the material frame formulation, the virial stress has no kinetic part while the potential part depends on in-terparticle forces and distances that are averaged separately.

In the spatial frame formulation, the expression for the virial stress contains both kinetic and potential parts.

In the present paper, we generalize the expression for stress tensor in the material frame formulation obtained by Hoover [2]. The long wave approximation [23-27] is used for transition from discrete description to equivalent continuum. It is shown that the expression derived by Hoover [2], with minor modifications, can be used for calculation of local stresses in static and dynamic problems. In the case of homogeneous deformations and finite temperatures, this expression leads to exactly the same values of stresses as classical definition used in continuum mechanics (force per unit area). In addition, comparison with commonly used equivalent stress tensors is carried out for several examples with uniform and nonuniform stress fields. The advantages of the derived expression for the stress tensor are demonstrated.

2. Hypotheses

Let us consider a discrete system consisting of particles arranged in an infinite ideal crystal lattice. Herein, only crystals with simple structure are investigated (i.e. crystals that are invariant to translation on any vector connecting lattice nodes).

We use two main principles for transition from discrete system to an equivalent continuum: (i) a decomposition of particle motions into slow continuum and fast thermal modes [25-27] and, (ii) the long wave approximation [23, 24]. First let us focus on the particle motion decomposition, as this approach has been intensely examined by many researchers in recent years. In the literature, it is effected using different types of averaging such as: spatial averaging, time averaging, averaging over phase space or over frequency spectrum, etc. In [28] it was noted that unique decomposition is impossible because rules for a choice of averaging parameters like averaging time, representative volume, etc. do not exist. The only restriction for these parameters is that they should depend on time and spatial scales of the problem being solved, and some work to quantify these scales has been done in [29].

Let us denote average (f ) and oscillating (thermal) f components of physical quantity f as:

f = if) + f, f=f-(f)- (1) Different expressions for the averaging operator (•) are proposed in the literature.

For the case of a one-dimensional chain of N particles, the following operator was defined for the velocity v [25]:

1 t+TI2 n+A/ 2

iv)n =— J E v(T)dT . (2)

T A t-Tl2 k=n-A/2

Here iv) n denotes the averaged value of physical quantity v over a spatial region centered at n and comprised of the A particles surrounding n (inclusive of n), and over a temporal domain of size T that surrounds the current time t.

As an alternative approach, the direct and inverse Fourier transformations were used for decomposition (1) in [28]. The Fourier transformation of v is given by:

V (v) = <V> (v) + V (v) = J ve^dx,

V2n _

where

<vxv) -Pf"" V (v) =

0, v>v„„f,

0 ^^ V(v), v>v c

(3)

(4)

Here vcut is the cut-off frequency, which is on the order of THz [23] and depends on the material. The quantities (v) and v are obtained through the inverse Fourier transformation:

<v> = J Ve_vtdv, v = J Ve_'vtdv. (5)

_L

The approaches developed by Hardy [16-19] and Lucy [20, 21] employ the kernel density estimation; in particular, the velocity is given by

E mk v kw (x - rk)

<v>(x, t) =-

(6)

E mw (x - rk)

k

Here vk, mk, rk are the velocity, mass and position of particle k; x is the coordinate of the spatial point, where the velocity is calculated; and w is the localization/kernel function that determines the weights in this weighted average. Applying decomposition (1), the thermal component of the velocity for particle n is

v„(x, t) = v„(t) -(v)(x, t) (7)

and is associated with both the particle n itself and the spatial point x. This type of decomposition (as well as all those discussed) is not unique. It strongly depends on the choice of localization function. In particular, if localization volume surrounds a single atom, then the thermal component of the velocity is equal to zero.

Since a unique averaging operator does not exist, a formulation for estimating continuum fields should not be dependent on the specific choice of the operator, nor should the results depend qualitatively on this choice. In the following derivations, properties of the particular averaging methods are not used unless otherwise stated.

Another principle used in the present paper (independent of the decomposition) for transition from a discrete system to an equivalent continuum is the long wave approximation [23, 24]. The approximation assumes that if an average component of any physical quantity is slowly changing in space over distances of order of the interatomic spacing, then the average component can be considered as a continuum function of a space variable and can be expanded into power series with respect to the interatomic distance. Moreover, the resulting series should converge rapidly. Development and use of this approach is examined in further detail in the next section.

3. The approach based on the long wave approximation

3.1. Kinematics

We consider two states of a discrete solid and its equivalent continuum: the reference configuration and the current configuration. It is assumed that the mapping between these configurations exists. This assumption is usual for solid mechanics theories based on the material frame formulation. For the sake of simplicity let us take the undeformed configuration of the crystal lattice as the reference. Using a local numbering convention, similar to one established in papers [26, 27], we focus on a reference particle "0" and mark all of its neighbors by the index a. The vector that connects 0 with a in the reference configuration is a0a. By the definition, vectors a 0a have the following property

a0(-a) =-a0a • (8)

Note that no thermal motion is associated with the reference configuration, i.e. a0a = 0. In a current configuration vector connecting particle 0 to its neighbor a is denoted

A0a :

def ~

A 0a = A a- A 0 = <A 0a) + A 0a , (9)

where A0, Aa are radius-vectors of particles 0 and a. The use of lower-case letters to denote reference configuration and upper-case letters to denote current configuration was done in [26, 27], and we maintain this notational style here for consistency.

Let us consider kinematics of the discrete system in the long wave approximation. It is assumed that average values of particle positions are identical to positions of corresponding points of continuum media. Positions for points of equivalent continuum in the reference and current configurations are denoted as r and R respectively. Thus, the position of particle 0 is represented by the vector r0 in the reference configuration of the equivalent continuum, and is represented by the vector R(r0 ) = (A0 ) in the current configuration. Note that R(r0 ) is a mapping that implicitly depends on time as well as the reference position r0. Here, we have dropped this dependence on time for brevity. The average position of particle 0's neighbor a is determined by vector R(r0 + a0a). Then vectors A0a and a0a connecting the particles are related by the following formula

< A 0a) =< A a)-< A 0 ) =

= R(r0 + a0a ) _ R(r0 ) = a0a ' R(r0 X

(10)

° def

where V = d/dr is the spatial gradient operator in the reference configuration [5].

This expression is the Taylor series expansion of the equivalent continuum position for particle a relative to the position for particle 0, truncated to first order. One can see that expression (10) is similar to the formula used in continuum mechanics that connects material vectors dr and dR. In the literature formula (10) is usually called Cauchy-Born

rule (see, e.g., [30]). If particle positions are known, then formula (10) can be used for calculation of deformation gradient V R(r0) for equivalent continuum:

V R(r0) = 1 E a0aa0a I -E a0a<A0a >

-1

(11)

Note that in the long wave approximation the function R(r) is assumed to be slowly changing on the distance of order of | a0a |. Therefore defining the function R(r) at particle positions is sufficient for calculation of gradients.

Using Eq. (10) one can derive relations between vectors i A 0a), a 0a and measures of deformation used in nonlinear theory of elasticity [5]. For example, the follo-

def

wing identity is satisfied for Cauchy-Green measure G =

def o o t

= (V R) • (V R )T:

1 iA0a) |2 = a0a • G • a0a - (12)

Thus, the long wave approximation allows us to connect deformations of the discrete system and the equivalent continuum.

3.2. Equivalent stress tensor

Consider dynamics of infinite perfect crystal of simple structure. Let us assume that the total force F0 acting on the reference particle 0 is represented in the form:

Fo =E Foa ({A 0a }a.A )- (13)

a

Here F0a is the force acting on the reference particle due to its neighbor a, {A 0a }aeA is the set of all vectors A 0a from particle 0 to its neighbor a lying within the set A. It is assumed that the decomposition (13) is carried out in such a way that the forces F0a satisfy Newton's third law with respect to the neighbor a. The example of decomposition (13) in the case of three-body forces is considered in Sect. 4.2.

Note that decomposition (13) used in the present paper, as well as many other papers dealing with calculation of stresses [15, 19] for systems with multibody interactions is not the only possibility. Alternative ways of decomposing atomic forces and constructing a corresponding expression for the stress tensor are derived in [31]. The application of this approach is demonstrated in next section using simple example with three-body forces.

Now let us derive the equation of motion for the equivalent continuum using the equation of motion for the reference particle and the decomposition of motions:

mi v 0 X =E F0a ({A 0a} oeA

a (14)

mv 0 =E F0a ({A0a }aeA X

a

where m is the mass of particle 0. The first equation in (14) describes slow motion of the system that corresponds to average value of particles velocity (v0). The motion (v0) can be considered as the motion of a continuum media. The second equation describes the thermal oscillation charac-

terized by thermal component of displacement v0. One can see that both iF0a) and F0a depend on the total particle-neighbor distances, A 0a=iA 0a) + A 0a, and thus couple two equations in Eq. (14).

We assume that average interparticle forces can be approximated by continuum functions fa (r), i.e.

fa M = <F0a), f-a M = <Ffl(-a) ) (15)

Evidently this assumption is satisfied exactly in the case of homogeneous deformations, since in this case fa are constant. The functions fa (r) satisfy Newton's third law:

fa (r0 ) = -f-a (r0 + a 0a ),

f-a (r0) = -fa (r0 - a0a)'

By their physical meaning, the functions fa (r0) are similar to traction vectors at material point r in directions a0a. Using the introduced functions fa (r), one obtains

iF0a + F0(-a)) = fa (r0 ) + f-a (r0 ) =

= fa (r0) - fa (r0 - a0a ) ~ a0a •V fa (r0) =

= V7^(a0afa(r0)). (17)

Given that the functions fa are slowly changing between atoms and, hence, higher order terms in the expansion can be neglected via the long wave approximation. Let us carry out the following transformations in the equation of motion of particle 0:

mv0 = E F0a = E F0(-a) = 1/2 E (F0a + F0(-a) ). (18)

Substituting formula (17) into averaged equation (18) and dividing both parts by elementary cell volume in the reference configuration V* one obtains

„ ( 1 >

- < V0 > =V-

T1" E a0afa (r0)

2V a

(19)

Let us compare formula (19) with the equation of motion for continuum in Piola's form [5]:

p*v = V • P, (20)

where P is the Piola stress tensor, p* = m/V* is a density in the reference configuration. Comparing equations (19) and (20) one can deduce that Piola stress tensor at the point r0 has the form

P(r0) = E a0aiF0a), (21)

2V* a

where the identity (15) was used. Strictly speaking the formula (21) determines the stress tensor to within a field of zero divergence. This tensor field corresponds to some equilibrium stress in the crystal. Note that the only assumption made during the derivation is represented by formula (17). In the case of homogenous deformation of the crystal formula (17) is satisfied exactly. Thus in this case the expression (21) for Piola stress tensor is also exact. Note that the expression for Piola stress tensor is used, in particular, for investigation of stability of crystals [32]. In this case it is

more convenient for analytical derivations than the Cauchy stress tensor.

Let us carry our the same derivations in a current configuration using the following transformations: f-a (Ro) = f (Ro +<Ao(-a) » = = -f„ (Ro) + <Ao«>-Vf„ (Ro), (22)

where V = d/dR is a nabla operator in the current configuration. In formula (22) it was used that in the long wave approximation ( A o(-a) > = -( A oa>. In the case of homogenous deformation of the crystal (22) is exact. Using formula (22) one can rewrite (17) in the current configuration:

(Foa + Fo(-a) > = fa (Ro ) + f-a (Ro ) ~

= <Ao«>-Vf„ (Ro). (23)

Substituting expression (23) into equation (18) and dividing both parts by the volume of elementary cell in the current configuration V one obtains

m < vo >'=^ E< Ao«>-Vfa (Ro) =

V 2V a

= V-|^E<Aoa>fa l-EV

2V

1

<Aoa> Ifa. (24)

2V

The argument R0 of functions fa is omitted in the right side of the given equation for brevity. The second term in the right side is equal to zero, since

2V <A0a

>! = ^ v.

1 2V

(

I

V o T V*(VR )T

V

4-,-^

=o

0a

= 0, (25)

where Eq. (10) and Piola's identity were used (see, for example, [5]). Then the equation of motion (24) has the following form

1

-< v o >' = V

E< A oa>fa ( R o)

(26)

V ' { 2V'

The requirement of equivalence of discrete and continuum systems leads to expressions for the Cauchy stress tensor

T(R o) = E< A oa><Foa>

2V a

(27)

and density p = m/V in the current configuration. The expressions (21), (27) for the Piola and Cauchy stress tensors satisfy the following well-known relation

t = V*(V R )T • P/V

(see, for example, [5]).

Note that the expression for the stress tensor (27) does not explicitly depend on particle velocities. Also, unlike usual expression for virial stress tensor, in formula (27) the forces and vectors connecting particles are averaged separately. However for solids it does not lead to any contradictions. In [2] it is shown that for solids there are two equivalent formulations of virial theorem. According to the theorem the stress tensor can be represented either using (AoaFoa> and kinetic part or using (Aoa> (Foa>. In the latter case there is no kinetic part. The expression for volume-averaged stress calculated using formula (27) coin-

cides with the expression from [2]. The derivation in [2] is based on the assumption (A0a) = 0, valid only in the case of statics. The derivation given in the present paper shows that expression (27) has much wider range of applicability. It satisfies the law of momentum balance (26) for equivalent continuum and therefore can be used for calculation of local stresses in both static and dynamic problems.

In the next paragraph, we will demonstrate that in the case of homogenous deformations and finite temperature, formula (27) becomes exact (see Figs. 1, 2). Note that the only assumption about the interparticle forces used in the present paragraph is given by formula (13).

4. Comparison of different expressions for Cauchy stress tensor

Let us compare the expression for stress tensor (27) with analogous expressions used in the literature. In the following examples the operator (•) in formula (27) corresponds to time averaging. The expression for the Cauchy stress tensor at spatial point x and time t using the Hardy formalism [16] is

1

th (xt ) = 2 EE A jFjBj(x) - E m v iv i-w(x - Ai \ 2 i j i

Bj (x) = J w( A j -ÀAj - x)dA,

(28)

where A j = A j - A t; A t is the position of particle i and vt is the thermal velocity as defined in Eq. (7).

Another approach developed by Lucy [20] is used in [21]. The following expression for Cauchy stress tensor at spatial point x and time t is proposed:

Tl(x, t) = E

(

E Aij Fij - mi Vi Vi

I

w(x - Ai ),

(29)

where thermal velocities vt are calculated using formulas (6), (7). The given approach avoids integration over bond length. In the present work, the localization function proposed by Lucy [20] is used with both expressions (28) and (29):

f r,Y r, f

w(r ) =

_5_ nR

o, r > Rc

r 1 + 3 — 1 --1

Rc c . Rc,

r < R

(30)

where Rc is a range of localization function. This function has two smooth derivatives everywhere and therefore is appropriate for calculation of continuum fields. Note that in contrast to formula (27), expressions (28) and (29) explicitly depend on particle velocities.

Using molecular dynamics simulations, let us examine the consistency of the expressions just described with continuum mechanics. In all the examples considered below particles interact via the spline potential [33]. Corresponding interatomic force Fj acting on particle i due to particle j has the form

F _ ^ k ( )

a

ГГ Л 8

a

V J J

Г V44 a

V J J

A j '

(31)

k ( а ) =

1,0 < Aj < b,

- b2

2 l2 vacut - b ,

0, Aij ^ «cut,

, ь < а < «cut '

where b = (13/7)1/6a, a is an equilibrium distance, e0 is a bond energy. The force (31) coincides with Lennard-Jones force for Aj < b and smoothly goes to zero as Aj goes to

acut-

4.1. Example 1: Cold pressure in a crystal

First, consider a two-dimensional triangular lattice system compressed volumetrically by 9.75 % (bonds are compressed by 5 %). Let us calculate the pressure at points where the particles are located using formulas (27)-(29). In the given calculations the cut-off radius is equal to 1.5a. The pressure corresponding to formula (27) is calculated analytically:

P_

Ps

6 MA2

Vd

a

14

V _

V5—

(32)

where A = 0.95a; d = 2 is a dimensionality; M = 6 is a coordination number; ps =e0 /a2 is usual Lennard-Jones scaling. The resulting value is p/ps ~ 11.29033. Molecular dynamics simulations are carried out in order to calculate pressure using formulas (28), (29), and classical continuum definition. In the latter case, the average normal component of the force acting on one layer of particles per unit length is calculated. The resulting value coincides with prediction of formula (32). Pressures calculated using the Hardy (28) and Lucy (29) stress tensors converge to the same value with increasing range of the localization function Rc. Trapezoidal rule with 50 points per bond length is used for calculation of integrals in Hardy's expression (28). The dependence of relative error in pressure on Rc is shown in Fig. 1. Note that for Rc > 1.9a the difference between the Lucy and Hardy stresses is less than one percent and is decreasing with Rc.

One can see from Fig. 1 that in contrast to formula (27), formulas (28) and (29) overestimate the pressure, at least for Rc < 2a. At the same time, formula (27) gives the exact value of pressure.

4.2. Example 2: Shear of a square lattice with three-body forces at finite temperature

Let us demonstrate the possibility of using formula (27) for calculation of stresses in systems with multibody interactions. Consider the simplest problem with three-body interactions, notably pure shear of perfect square lattice at

finite temperature. In this problem the different definitions of stress can be compared analytically. Let particles interact via angular springs connecting two neighboring bonds. The average length of all bonds are equal to A. Only nearest neighbors are taken into account. Average positions of the particle number 0 and its neighbors are shown in Fig. 3.

Assume that potential energy of the spring connecting bonds A01 and A02 is given by the following expressions

U012 = U012(A)1, A)2, A12)' A12 = A02 - A01. (33)

In order to calculate stresses using classical definition (force per unit area) let us introduce interparticle forces Fj using the following definition [15, 19, 34]

Aij

. , (34)

AJ

where Utot is the total potential energy of the lattice. Forces

F _ dUtot

ij ЭА,

e , e _ ij ij

F- defined by formula (34) are central and satisfy the relations F = ^ F-, F- = -F-, where Fi is the total force acting on particle i. Using formula (33) and definition (34) one can calculate the forces F10, F12, F^-2^:

F10 _ -2

dU

012 + dU01(-2) Л

dA01 dA01

c01'

F12 _ 2-

dU

012

dA

"^12' F1(-2) = 2 el(-2). a12 dA1( -2)

The remaining forces in the system are calculated using symmetry of the problem and the third Newton's law. Let us calculate the stress vector tcl acting on crystal cross-section with normal n orthogonal to vector e02 (see Fig. 3). According to classical continuum mechanics definition, the stress vector in two dimensions is equal to the force acting on the cross-section per unit length. Then t cl has the following form

dU

01( -2)

(35)

t cl _ — ( F10 + F12

A

+ F1( -2) >,

(36)

Fig. 1. Relative error in pressure calculated using Lucy and Hardy definitions for cold, two-dimensional crystal as a function of localization radius Rc. Here pcl/ p s = 11.29033 is a pressure calculated using classical continuum mechanics definition (force per unit length). Triangles corresponds to pressure calculated using formula (27)

where A is the nearest-neighbor distance. Now let us use formula (27). Using symmetry of the problem one can rewrite formula (27) as follows:

T = V1« A 01X F01 > + < A 02 X F02 ) +

+ <A 03 ><F03) + < A 04 ><F04)), (37)

where the identity <A0iXF0i)=<A0(-i)XF0(-i)), i = 1, 2, 3, 4 is used. Substituting the relations <A03) = <A01 +

+A02), < A04 ) = < A02 - A01), <F03 ) = -<F1(-2)), <F04 ) =

= <F12) into formula (37), one obtains the expression for the stress tensor: 1

* = -■V <A oi><Fio + f12 +Fi( -2) >-- V < A 2o >< F2o + F21 + F2(-1) >.

(38)

Then the stress vector acting on the cross-section with normal n orthogonal to vector e02 (see Fig. 3) is calculated using the Cauchy formula t = n • t. As a result

t = n T = -

n-<A o1 >

<F1o + F12 + F1(-2) >= t

(39)

Here the expression V = -An • <A01) for the volume of the elementary cell was used. From formula (39) it follows that the stress vector calculated using formula (27) exactly coincides with classical continuum definition (36).

Alternative approach for calculation of the stress vector in systems with multibody interactions is suggested in [31]. It is stated that every angular spring causes three forces with zero sum. For example, for the spring connecting particles 0, 1 and 2:

T-^12 .j?o12 .j?o12 л тo12def dUo12

Fo + F1 + F2 = o, Fi =-"

9A,

(40)

Here i = 0, 1, 2. The forces F°12, i = 0, 1, 2 contribute to the potential part of the stress vector only if the corresponding angular spring is dissected by the cross section. The sign of

Fig. 2. Relative error in local thermal pressure calculated using formulae (27)-(29) as a function of localization radius Rc. Here pcl is the pressure calculated using classical continuum mechanics definition (force per unit length)

contribution depends on position of the particle i with respect to the cross section. Then the potential part of the stress vector is equal to sum of contributions from all crossed angular springs. In the present example it has the form:

tPH°KB=A <F012 + Fi1(-2) - Fo°12 -

F o1( -2) Fo

-f2°12

_F01( -2),

(41)

Here it is used that angular springs connecting particles 0, 1, 2 and 0, 1, -2 cause the same average forces as angular springs connecting particles 0, 1, -4 and 0, 1, 3. The forces in formula (41) are as follows:

F12 =

dU

o12

dU

o12

Hi

coi

ЭД)2

co2'

F012 __ dUoi2

dU,

oi2

F2o12 =-

dAoi dU

-oi

dA

12

oi2

dU

dAo2

F oi(-2) = dUoi(-2) o dAoi

eo2 "

oi2

dA

c12,

12,

12

dU

e oi

oi(-2)

(42)

dAo(-

e o(-2),

F01( -2)

F-21( -2)

dU

oi( -2)

dU

oi( -2)

dAoi

eoi

dA

e1(-2)'

l1( -2)

dU

oi( -2)

dU

eo( -2)

oi(-2)

dAo( -2)

Substituting (42) into formula (41) and using formulas (35)

4_2) 1("2)'

for the forces Fio, F12, F1(-2) one obtains:

,pot _

1

lHKB " A+ F12 + F1(-2) > = tcl"

(43)

Thus the potential part of the stress vector calculated using the Heinz-Paul-Binder formalism exactly coincides with the classical definition (36). In [31] it is suggested that the stress vector should additionally have kinetic part. The present example shows that at least for solids the kinetic part should not be added.

The expression for the stress tensor (27) gives the exact value of the stress vector. Therefore it can be used for calculation of stresses in systems with multibody interactions.

4.3. Example 3: Thermal pressure in a crystal

Consider calculation of local Cauchy stress in two-dimensional crystal at a finite temperature. Assume that ini-

Fig. 3. Shear of the square lattice. Averaged positions of the particle number 0 and its nearest neighbors

tial configuration of the lattice is undeformed. Periodic boundary conditions in both directions are used. Initial velocities of particles are chosen so that the sample has the following temperature: kBT/e0 = 0.05, where kB is the Boltzmann constant. Temperature is calculated as the mean kinetic energy of thermal motion [2]. Consequently, the pressure is due to thermal expansion only. The pressure is calculated using formulas (27), (28), (29). Lucy and Hardy expressions (28) and (29) are additionally time averaged over 2 -104 T0, where T0 is the period corresponding to the Einstein frequency. In the case of formula (27), time averaging is carried out using operator (•). In all three cases standard error of the mean is approximately 0.2 %. As in the previous example, the Lucy and Hardy stresses are calculated for 1 < Rc/a < 3. The results are compared with classical continuum mechanics definition of stress (force per unit length), which is used as a benchmark. The deviations from the benchmark are shown in Fig. 2. One can see from Fig. 2 that the error of the Lucy and Hardy stresses is decreasing with increasing Rc. The error for Rc > 3a is less than 1 %. Additionally let us note that the Lucy and Hardy stresses are indistinguishable for Rc > 2.1a. At the same time, the pressure obtained using definition (27) coincides with the classical continuum definition within the standard error. Thus at homogeneous deformations and finite temperature, the Eq. (27) can be used as a benchmark.

4.4. Example 4: Stresses around a pair of dislocations

Let us consider the case of an inhomogeneous stress distribution and compare the results with the prediction of linear elasticity theory. Consider the stress field generated by a pair of edge dislocations as modeled by a two-dimensional crystal lattice with pair interactions. An analogous problem in three dimensions is considered in [17]. It is well-known that triangular lattice is isotropic in the case of small deformations. Thus, one can compare the stress field induced in the triangular lattice with analogous stress field in two-dimensional isotropic linear elastic continuum. First, let us derive expressions for stresses using well-known results for three-dimensional continuum under plane strain conditions [35].

Consider an infinite three-dimensional continuum under the plane strain conditions containing single dislocation. Assume that x axis is directed along the Burgers vector and dislocation line coincides with z axis. Then the normal stresses have the form [35]:

_ _ bG* y(3x2 + y2)

" 2n(1 -v*) (x2 + y2)2 '

2 2 (44)

__ bG* y(x2 - y2)

Tyy _ 2n(1-v*)(x2 + y2)2 ' where G*, v* are the shear modulus and Poisson's ratio of the continuum, b is the modulus of the Burgers vector. Stresses in a two dimensional continuum with a dislocation

can be obtained replacing G* with the shear modulus of a two dimensional continuum G and v* by v/(v +1), where v is the Poisson's ratio for a two dimensional continuum. For a triangular lattice with spline interactions,

v = 1/3, G = 18V3ps. (45)

Let us exploit the principle of superposition to determine the complete stress field for a system of two dislocations with equal and opposite Burgers vectors separated in space. This field is

Txx (x y) = Txx (x y, b) + Txx (x - x0, y - yo, - bl

t (46)

Tyy(x y) =Tyy(x y, b) +Tyy(x - Xo, y - y,, - b),

where the first dislocation is at the coordinate system origin and the second dislocation has coordinates {x0, y0}.

Now let us compare analytical stresses (46) with results of molecular dynamics simulations. In the given example the cut-off radius is equal to 1.75a. The dislocations are created using the following procedure. Rectangular sample containing 40 000 particles is considered; 24 particles are removed as it is shown in Fig. 4 (see [36]).

Interatomic distances in the initial configuration are set to be equal 0.995a in order to collapse the set of vacancies shown in Fig. 4. Additionally the crystal is compressed by 0.25 % in x direction in order to favor formation of dislocations with Burgers vectors in the same direction. Periodic boundary conditions are used during the formation of dislocation. Finally boundaries are released and crystal is equilibrated at zero temperature. The process of formation of the dislocations is shown in Fig. 4. Coordinates of the second dislocation are xo = 12a, y0 = 12V3a. The absolute values of the Burgers vector for both dislocations is equal to a. The distribution of stresses along axis y at x = 0 (see Fig. 4) is calculated using formulae (27)-(29). Localization radius R = 2.1a is used in formulae (28), (29). The given value corresponds to minimum error in pressure in the case of homogeneous deformations (see Fig. 1). In the case of formula (27) stresses are calculated at atomic positions. Stresses in between atoms are obtained using linear interpolation. For the sake of simplicity it is assumed that V = V* =Sa 2 /2. The resulting distributions of stresses Txx, T along y axis are shown in Fig. 5.

One can see from the figure that in contrast to the continuum solution, the results of molecular dynamics simulation do not contain singularity at y = 0. Obviously the singularity is a mathematical consequence of the assumptions of linear elasticity. The stresses calculated using formula (27) proposed in the present paper are closer to the prediction of elasticity theory (44), (46). Formulae (28) and (29) give analogous results only for |y| > 5a. The average difference between the Lucy and Hardy stresses is less than 1 %. Better agreement between formula (27) and continuum theory is caused by the fact that, in contrast to the Hardy and Lucy formulas, formula (27) does not contain spatial averaging.

jjjjjjjjjjjjjejjjsljjjjjjo iJJJJJJJJJJJviJJJJJJJJJJJJ. JOJJJJJ^JOJJJJJJJJJOJJJJJ IJJJOJJJJJJJOJJJJJJJJJJJJ. JjJJJJJJJjJJJJJJJJJjJJJJJ

JJOJJJJJJJJJJJJJJJJJJJJJJ iOJJJ^JJJ^JJJOJJJJJJJOOJOL

iJJJJJJJOJJJJJOJO^l JJJO^O, jJJJJJjjjjjjjjjjj jjjJOJj „, JJJJJJ^JJJJ JJJJJJJ. JJJ; У JJJJJJJJJJJ JJJJJJJJ iJJO ^ J JJ JJJO JJ J JJJJJOVXIL

--------IJJJJJJOJ JJJJJJJJJ

J J J J J J J j jjJJJJJJj-. UJJ.JJJJ JJJJJJJOJJ JJJJJJJ JJJJJJJJJJ..

»«••«•9 JJJJJJJJJJJ JJJJJJ jj jjj jjj jjj..

^JJJJJ JJJJJJJJJJJJ

J j j J _ _

' j « • • • ol jjjjjjj

' j J J J J j

■Jjjjjj^jjjj jjjjjjjjjjjjj

'JJJJJJJJJJ j J J J J J J J J J J J j . ••••«•О«** JJJJJJJJJJJJJJ Ijjjjjjjjj J j J J J J J J J J J J J J, JJJJJJ^JJ JJJJJJJJjjjJjjj

-------1 J j J j JJJ J J J J ■ J OJ О.

JJjJjJJJJjJJJJJJ JJJJJJJJJJJJJJJj. JJJJJJJJJJJJJJJJJ JJJ JJJ JJJJ JJJJ JJJ..

------ , j j j j j j j j j J»0 J j J j j

i jjjJjjjJjjjjjjjjJC ,.-)»«*«<

Jjjjjjjjjjjjjjjjjj X jjjjj

' J J J J J J J J J J J j J J J J Ji. Ijjj9<

' jjjjjjjjJJJJJJjJjjjjjjjj, JJJjJJJJJJJJJJJJJJJjJJJJJ

JJJJJJJJJJJJJjjJJJJJJJJJJ JjJJjJJJJJJjjJJJJJJJJJJJ.. JJJJJJJJJJJJJJJJJJjJjJJJJ jjjjjjjjjjjjjjjjjjjjjjjj. JJJJJJJJJJJJJJJJJJJJJJJJJ JjjjjjjjjjjjjjjjjjjJJJjJ* JJJJJJJJJJJJJJJJJJJJJJJJJ jJJJJJJJJJJJJJJJJJ^OOJJJ-JjjjjjjJJJJJJJJJJJ О»JJjj jJjjjjjjjjjjjjjJJ }JJJJJJJJJJJJJJJ

J j j j j j j j j j j j j j j »Jtf • • • °

J J j j j j j j j j j j j j l^j* * -

jjjjjjjjjjjjjj • ® ® • e

i j j j j j j j j j j j j j • ••♦•Л*Л*Л4

jjjjjjjjjjjjjjjjjjjjjjjj

"••#VMV»VMV»vavav *g

•••ЭjjjjJOjjj J J J j j j j J j j j^ I jJJJJjJjJjJJJJJJJJJJJJJJ jjjjjjjJ J J J j J J j J j j

jjjj •"•VAt*VAVAVAfA*

JJJJjjjjjjjJJJ^JJJJJJJJJ ■ J j j J J j j j j j J J J J J J J J J J J J J ijjj jjjjjjjj J JJ j j J J J J J J J VmViVmV1

J j J j J J jjjJjjjjjjjjjjjjjj

• * JJjjjj JJJJJ j J J j j j j J JJ с JjjjJjjjJj j j j j j j j Jj j j J J J jjj Лаа» J J J J J J J J J J J J J j j > • * j jjjj JJJJJJJ Jjjjjjjjj JJJJJJ JJJJJJjjJjJjJJJJJ, »••a«« »••»99»»j«»99 jjjJJjJjJJjjJJJJjjjJjJJJ-JjjjjjjjjjjjjjjJjjjJJJJJJ jjjjJJjjJjjjjJJjjJjJjjjj, JJJjJJjJJJJJJJJJJjJJJJJJJ jjjjjjjjjJJjJjjjJJjJJJJJ^ JJJJJJJJJJJJJJJJJJJJJJJJJ

J J J J j j j J j j J j J j j J JO j J J J J jv JjjJjjjjJjjjjjJJjj'JJjjJJj

J j J J j j j j j j J J J j j j j^ J j j j j js

JJJJJJJJjjjJOjJjjJOjjjjjj IjjjJjJjJjjjJjjJJJJOjJJJj..

JJJJJJJJJJJJJJJJJJ jjjjjj

j j j JO J j J J J j j j J J j jjjjjjjjjjjjjj

ijjjjjjjjj ••••••♦•AVA*A<

JjjjjjJJJJjjJJJJJOjJJJJJ iJJOJJJJQJJJJJJJdJJJ j J J J ■ J jjjJjjjjjJJJjJJJJOJJJJJJ ijjjjjjjjjjjjjjjjjjjjjjj^; JJJJjJjjjjOjJjjJJJJjjjJO-

IjjjjjJjjjJjjJJJJJJJJJJJJ jjjjjjjjjjjjjjJJjjjJjjjj

UJJJJJJjjjJJ

jjjjjajjjjjjj

tjjjojjjjjjjjjj•V >jjjjjjjjjjjjjjjjjjjjjjjj

jjjjjjjjjjjjjjj«JJJJJJJ

iJJjjJJJjjjjOjjJJJJJJJJO^ joojjjjjjjjjJjjjjjjjjjjj

IJJJJJJJOJJJJJOjjjjjjjjjj JJJJJJJJJJJOOJJOJJJJJJJJ >jjjjja^jjjji>JjJJJJjJJJJ^ JjjjJJ^jJJjjjJjJJJjjjjjj jjjjjjjjjjjjjjjjjjJJJJJJ

J j J J j 'JJ jjjjjjjjjjjjj Jjjj ............

JjjjjjjjjjjJJjjjJJJjjjjjj

j j j j j »Л jjjjjjjjjj •••••••

jjjjjjjjjjjjjjjjjjjjjjjj^

JJJJ'JJjjjjJjJJJJjejJjJIjJJ

Fig. 4. Formation of the pair of dislocations. Only atoms near the dislocations are shown

Finally, let us note that comparison of stress distributions obtained using formulas (27)-(29) for |y| < 5a is not so straightforward. On the one hand, the difference between stresses at some points in this interval is relatively high. On the other hand, given the limitations of linear elasticity, we would need to appeal to a more complex model of the continuum (e.g [37]) to determine which formulation is more accurate. Such an analysis was performed for the Hardy stress in [17] in three dimensions.

5. Results and discussions

An approach for the calculation of equivalent continuum parameters for discrete solids in the material frame formulation was presented. Two main principles were used for the transformation: the decomposition of particle motions into continuum and thermal parts, and the long wave approximation [23, 24]. The relation between kinematics of discrete system and kinematics of equivalent continuum was established. Equivalent Cauchy-Green measure of deformation for discrete system was introduced. The transition from a single particle equation of motion to equation of motion for equivalent continuum was carried out using the long wave approximation. No assumptions about interpar-ticle forces were used. Expressions connecting the Cauchy and Piola stress tensors with averaged interparticle forces and distances were derived for the case of multibody inter-

4

2

£ 0 H

-2 -4

-10 -5 0 5 10

y/a

Fig. 5. The dependence of t (a) and tw (b) ony at x = 0. The lack

actions. In the case of pair interactions the volume-averaged expression for the stress tensor exactly coincides with expression derived by Hoover [2].

Four test problems with homogeneous and nonhomo-geneous stress fields and finite thermal motion were considered. In the case of pair interactions the expression for Cauchy stress tensor (27) was compared with the Hardy [16] and Lucy [20, 21] expression. Additionally, it was shown that in all considered examples the difference between the Hardy and Lucy stresses is of order of 1 %, noting that the Lucy expressions is more efficient from the computational point of view. The calculation of stresses in systems with multibody interactions was demonstrated using the simple example, notably pure shear of square lattice at finite thermal motion. The stress vector was calculated analytically using formula (27) and Heinz-Paul-Binder approach [31]. The stress vector calculated using formula (27) coincide with classical continuum definition and the potential part of the Heinz-Paul-Binder stress vector.

It was shown that in the case of homogenous deformations and finite temperatures the stresses calculated by definition (27) exactly coincide with classical continuum stresses (force per unit length). The Hardy (28) and Lucy (29) expressions give the same result only if the averaging over a sufficiently large volume is used. In the case of non-homogeneous deformations the stress calculated using formula

4

2

-2 -4

-10 -5 0 5 10

-y/a

of symmetry is caused by the presence of the second dislocation

(27) is closer to prediction of elasticity theory. Better agreement is caused by the fact that, in contract to the Hardy and Lucy expressions, formula (27) does not contain spatial averaging.

The authors are deeply grateful to Prof. William Graham Hoover for the fruitful discussions of the present paper. This work was financially supported by Sandia National Laboratories and European Research Agency (FP7-PE0PLE-2009-IAPP Marie Curie IAPP transfer of knowledge program, Project Reference No. 251475).

References

1. Mechanics — From Discrete to Continuous / Ed. by V.M. Fomin, A.N. Andreev, et al. - Novosibirsk: Izd. SO RAN, 2008. - 243 p.

2. Hoover W.G. Molecular Dynamics: Lecture Notes in Physics. - Berlin: Springer, 1986. -142 p.

3. Psakhie S. G., Korostelev S.Yu., Smolin A.Yu., Shilko E.V., Dmitriev A.I., Horie Y., Ostermeyer G.P., Pegel M., Blatnik S., Zavsek S. Movable cellular automata method for simulating materials with mesostructure // Theor. Appl. Fract. Mech. - 2001. - V. 37. - No. 1-3. - P. 311-334.

4. Love A.E.H. A Treatise on the Mathematical Theory of Elasticity. -New York: Dover, 1944. - 643 p.

5. Lurie A.I. Nonlinear Theory of Elasticity. - Amsterdam: North-Holland, 1990. - 617 p.

6. Goldstein R.V., Morozov N.F. Mechanics of deformation and fracture of nanomaterials and nanotechnology // Phys. Mesomech. - 2007. -V. 10. - No. 5-6. - P. 235-246.

7. Panin V.E., Egorushkin VE. Nanostructural states in solids // Phys. Met. Metallography. - 2010. - V. 110. - No. 5. - P. 464-473.

8. Alyokhin V.V., Annin B.D., Babichev A.V., Korobeynikov S.N. Free vibrations and buckling of graphene sheets // Doklady Physics. - 2013. -V. 58. - No. 11. - P. 487-490.

9. Bunch J.S., van der Zande A.M., Verbridge S.S., FrankI.W., Tanenba-um D.M., Parpia J.M., CraigheadH.G., McEuen P.L. Electromechanical resonators from graphene sheets // Science. - 2007. - P. 315-490.

10. Liu W.K., QianD., GonellaS, LiS, Chen W, ChirputkarS. Multiscale methods for mechanical science of complex materials: Bridging from quantum to stochastic multiresolution continuum // Int. J. Numer. Meth. Engng. - 2010. - V. 83. - P. 961080.

11. Clausius R.J.E. On a mechanical theorem applicable to heat // Phil. Mag. - 1870. - V. 40. - P. 122-127.

12. Irving J.H., Kirkwood J.G. The statistical mechanical theory oftrans-port processes. IV. The equations of hydrodynamics // J. Chem. Phys. -1950. - V. 18. - P. 817-829.

13. Lutsko J.F. Stress and elastic constants in anisotropic solids: molecular dynamics techniques // J. Appl. Phys. - 1988. - V. 64. - P. 11521154.

14. Kunin I.A. Theory of Elastic Media with Microstructure. - Springer, 1982. - 291 p.

15. Admal N.C., Tadmor E.B. A unified interpretation of stress in molecular systems // J. Elast. - 2011. - V. 100. - No. 1-2. - P. 63-143.

16. Hardy R.J. Formulae for determining local properties in molecular-dynamics simulations: Shock waves // J. Chem. Phys. - 1982. - V. 76. -P. 622-628.

17. Webb E.B., Zimmerman J.A., Seel S.C. Reconsideration of continuum thermomechanical quantities in atomic scale simulations // J. Math. Mech. Sol. - 2008. - V. 13. - P. 221-266.

18. Zimmerman J.A., Webb E.B., Hoyt J.J., Jones R.E., Klein P.A., Bammann D.J. Calculation of stress in atomistic simulation // Model. Simul. Mater. Sci. Eng. - 2004. - V. 12. - P. S319-S332.

19. Zimmerman J.A., Jones R.E., Templeton J.A. A material frame approach for evaluating continuum variables in atomistic simulations // J. Comp. Phys. - 2010. - V. 229. - P. 2364-2389.

20. Lucy L.B. A numerical approach to the testing of the fission hypothesis // Astronomical J. - 1977. - V. 82. - P. 1013-1024.

21. Hoover W.G., Hoover C.G. New Trends in Statistical Physics: Festschrift in Honor of Professor Dr. Leopoldo Garcia-Colin's 80th Birthday / Ed. by A. Macias, L. Dagdug. - Singapore: World Scientific, 2010. - 360 p.

22. Xiao S.P., Belytschko T. A bridging domain method for coupling continua with molecular dynamics // Comp. Meth. Appl. Mech. Eng. -2004. - V. 193. - P. 1645-1669.

23. Born M., Huang K. Dynamical Theory of Crystal Lattices. - Oxford: Clarendon Press, 1988. - 420 p.

24. Krivtsov A.M. Influence of velocities dispersion on spall strength of material // Z. Angew. Math. Mech. - 1999. - V. 79. - P. 511-512.

25. Krivtsov A.M. From nonlinear oscillations to equation of state in simple discrete systems // Chaos, Solitons & Fractals. - 2003. - V. 17. -P. 79-87.

26. Krivtsov A.M., Kuzkin V.A. Microscopic Derivation of the Equation of State for Perfect Crystals // Proc. 6th Int. Conf. on Engineering Computational Technology / Ed. by M. Papadrakakis, B.H.V. Topping. - Stirlingshire: Civil-Comp Press, 2008. - Paper 145.

27. Krivtsov A.M., Kuzkin V.A. Derivation of equations of state for perfect crystals with simple structure // Mech. Solids. - 2011. - V. 46. -No. 3. - P. 387-399.

28. Zhou M. Thermomechanical continuum representation of atomistic deformation at arbitrary size scales // Proc. R. Soc. A. - 2005. -V.461.- P. 3447-3472.

29. Ulz M.H., Mandadapu K.K., Papadopoulos P. On the estimation of spatial averaging volume for determining stress using atomistic methods // Model. Simul. Mater. Sci. Eng. - 2013. - V 21. - P. 1501015015.

30. Arroyo M., Belytschko T. Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy-Born rule // Phys. Rev. B. - 2004. -V. 69. - P. 115415.

31. Heinz H., Paul W., Binder K. Calculation of local pressure tensors in systems with many-body interactions // Phys. Rev. E. - 2005. - V. 72. -P. 066704.

32. Podolskaya E.A., Krivtsov A.M., Panchenko A.Yu., Tkachev P.V. Stability of ideal infinite two dimensional crystal lattice // Dokl. Phys. -2012. - V. 57. - No. 2. - P. 92-95.

33. Le-Zakharov A.A., Krivtsov A.M. Molecular dynamics investigation of heat conduction in crystals with defects // Dokl. Phys. - 2008. -V. 53. - No. 5. - P. 261-264.

34. Kuzkin V.A. Interatomic force in systems with multibody interactions // Phys. Rev. E. - 2010. - V. 82. - P. 016704.

35. Kosevich A.M. The Crystal Lattice: Phonons, Solitons, Dislocations. -Berlin: Wiley, 1999. - 324 p.

36. Kaski K., Kuronen A., Robles M. Computer Simulation Studies in Condensed Matter Physics XIV / Ed. by D.P. Landau, S.P. Levis, H.B. Schüttler. - Berlin: Springer, 2001. - P. 12-26.

37. Eringen A.C. Edge dislocation in nonlocal elasticity // Int. J. Eng. Sci. - 1977. - V. 15. - No. 3. - P. 177-183.

Поступила в редакцию 04.03.2014 г.

Сведения об авторах

Кузькин Виталий Андреевич, к.ф.-м.н., снс ИПМаш РАН, зам. зав. каф. СПбГПУ, [email protected] Кривцов Антон Мирославович, д.ф.-м.н., проф., зав. лаб. ИПМаш РАН, зав. каф. СПбГПУ, [email protected] Jones Reese E., д.ф.-м.н., иссл. Сандийской национальной лаборатории, [email protected] Zimmerman Jonathan A., д.ф.-м.н., иссл. Сандийской национальной лаборатории, [email protected]