Continuum limit for three-dimensional mass-spring networks and discrete Korn's inequality

Continuum limit for three-dimensional mass-spring networks and discrete Korn's inequality

ARTICLE IN PRESS Journal of the Mechanics and Physics of Solids 54 (2006) 635–669 www.elsevier.com/locate/jmps Continuum limit for three-dimensional...

383KB Sizes 1 Downloads 46 Views

ARTICLE IN PRESS

Journal of the Mechanics and Physics of Solids 54 (2006) 635–669 www.elsevier.com/locate/jmps

Continuum limit for three-dimensional mass-spring networks and discrete Korn’s inequality M. Berezhnyya, L. Berlyandb, a

Institute of Low Temperature and Engineering, Ukrainian Academy of Science, Lenin Ave 47, Kharkiv 61164, Ukraine b Department of Mathematics and Materials Research Institute, Penn State University, University Park, PA-16802-6401, 218 Mc Allister Building, USA Received 15 July 2004; received in revised form 27 July 2005; accepted 17 September 2005

Abstract In a bounded domain O  R3 we consider a discrete network of a large number of concentrated masses (particles) connected by elastic springs. We provide sufficient conditions on the geometry of the array of particles, under which the network admits a rigorous continuum limit. Our proof is based on the discrete Korn’s inequality. Proof of this inequality is the key point of our consideration. In particular, we derive an explicit upper bound on the Korn’s constant. For generic non-periodic arrays of particles we describe the continuum limit in terms of the local energy characteristic on the mesoscale (intermediate scale between the interparticle distances (small scale) and the domain sizes (large scale)), which represents local energy in the neighborhood of a point. For a periodic array of particles we compute coefficients in the limiting continuum problems in terms of the elastic constants of the springs. r 2005 Elsevier Ltd. All rights reserved. Keywords: Homogenization; Elasticity; Mesoscale; Triangulized network; Discrete Korn’s inequality

1. Formulation of the problem and main result Consider a system of N point particles (N is a large number) located in a smooth bounded domain O 2 R3 . Introduce a small parameter  ¼ p31ffiffiffi ! 0 and denote by xðÞ i the locations of N Corresponding author. Tel.: +1 814 863 9683; fax: +1 814 865 3735.

E-mail addresses: [email protected] (M. Berezhnyy), [email protected] (L. Berlyand). 0022-5096/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jmps.2005.09.006

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

636

particles (we will omit the superscript  where it will not cause any confusion). We say that ðÞ the particles xðÞ i and xj are neighboring particles if the distance between them is of order . More precisely, there exist constants c1 Xc2 40 which do not depend on  such that ðÞ c2 pjxðÞ i  xj jpc1 .

(1.1)

We now introduce the mass-spring network. Assume that each particle has a finite mass 3 mðÞ i ¼ M i   ð0om1 pM i pm2 o1Þ, where m1 and m2 do not depend on . Particles are assumed to be of infinitely small volume; that is, we consider a system of point masses. Some neighboring particles (distances between neighbors are of order ) are connected by elastic springs, so that the system of particles and springs form a three-dimensional graph (network) whose edges correspond to the springs and the vertices correspond to particles. The only condition on these edges is that the triangulization condition in the sense of Definition 1 (see below) is satisfied. The properties of this graph which are essential in our consideration will be specified later (e.g., the degree of a vertex, which is the number of edges adjacent to a given vertex). Throughout this paper we assume that each particle has a finite number of neighbors which does not depend on . We define the interparticle interaction as follows: denote the displacement of the ith particle by uðÞ i ðtÞ; i 2 1; N. The displacements of the particles which are connected by elastic springs are small in the following sense: ðÞ juðÞ i  uj jpc.

(1.2)

This assumption is crucial for our consideration, and the applicability of the resulting homogenized formulas depends on whether this assumption is actually satisfied. Violation of (1.2) may lead to instabilities (Friesecke and Theil, 2002) where a lattice model with quadratic springs contains a geometric non-linearity. Assume that the springs provide a linear elastic interaction. Then the elastic energy of two neighboring particles is given by ðÞ ðÞ ðÞ hC ij ðuðÞ i  uj Þ; ui  uj i,

where h; i stands for a dot product in R3 and the matrix of elastic constants C ij for ith and jth particles (C ij  0 if these particles do not interact) is determined by a scalar spring constant kij : * ðÞ + ðÞ ui  uðÞ xðÞ j i  xj ðÞ ðÞ ðÞ ðÞ ij ij ðÞ ; e ¼ C  ðui  uj Þ ¼ k e ; e , (1.3) ij ij ij ðÞ ðÞ jxðÞ jxðÞ i  xj j i  xj j which means that we consider the central interaction between neighboring particles. Then the potential energy of the network is given by 1 X ij ðÞ ðÞ ðÞ ðÞ hC  ðui  uðÞ (1.4) HðuðÞ 1 ; . . . ; uN Þ ¼ H 0 þ j Þ; ui  uj i, 2 i;j where the summation is taken over all pairs of interacting particles, and H 0 is an arbitrary constant. We consider the elastic constants of the form kij ¼ kij 2 ;

k1 pkij pk2 ,

(1.5) ui

where the constants k2 Xk1 40 do not depend on . Since the displacements are of order  and the number of terms in the sum (1.4) is of order 3 , the scaling factor 2 in (1.5)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

637

ðÞ implies that HðuðÞ 1 ; . . . ; uN Þ in (1.4) is finite. Note that the condition k 1 40 is crucial in our consideration. If k1 ¼ 0 or k1 ! 0 as  ! 0, then an instability (due to the blow up of the Korn’s constant in (3.5)) may appear; for an example, see (Friesecke and Theil, 2002). Since the energy (1.4) is invariant under translations of the set of the particles as a whole, the energy minimization defines many equilibrium states. The unique minimum is defined by the condition that all particles that are at a distance smaller than c1  (see (1.1)) from the clamped boundary qO (the corresponding displacements uðÞ i are equal to zero). Thus, in N our consideration the system of particles has the unique equilibrium state fxðÞ i gi¼1 . Denote by M the number of particles which are located in the boundary layer UðqO; c1 Þ, M ¼ Oð2 Þ5N ¼ 3 . In what follows we will call such particles boundary particles. Then the motion of the network is described by a system of 3N ODEs: ðÞ ðÞ € ðÞ mðÞ i u i ¼ ruðÞ Hðu1 ; . . . ; uN Þ;

i 2 1; N,

(1.6)

i

subject to the following initial and boundary conditions, respectively: uðÞ i ð0Þ ¼ 0;

ðÞ u_ ðÞ i ð0Þ ¼ ai ;

uðÞ i ðtÞ  0;

t40

i 2 1; N,

ði 2 1; MÞ,

(1.7) (1.8)

where aðÞ i ði 2 1; NÞ is a set of given constants such that the kinetic energy N X

ðÞ 2 mðÞ i jai j pC

(1.9)

i¼1

is bounded uniformly in . Condition (1.7) means that the initial displacements of particles are zero and the system is driven by the initial velocities. Condition (1.8) means that the boundary particles are clamped. Here homogeneous Dirichlet boundary conditions are assumed for technical simplicity only. Generalization for non-homogeneous boundary conditions, as well as for other type of boundary conditions (Neumann, Robin), is straightforward as is usually the case in linear homogenization problems when the homogenization limit (the PDE or the constitutive equation) does not depend on the type of external boundary conditions. The goal of our study is to derive the continuum limit for the problem (1.6)–(1.8) as  ! 0. The understanding of the relationship between the atomistic (discrete) and the continuum moduli of solids dates back to the work of Cauchy–Born (Born and Huang, 1954), where linear elasticity equations were derived based on physical considerations for periodic mass-spring lattices. However, the compactness of the family of solutions of discrete problems was not established there. Such compactness is necessary for a rigorous justification of the continuum limit, which, in particular, provides the limits of validity for the formulas for macroscopic elastic moduli obtained previously by the Cauchy–Born rule. More recently, the applicability of the Cauchy–Born rule under various assumptions was investigated by a number of authors. In the work of Blanc et al. (2002) the continuum limit was derived based on the hypothesis that the microscopic displacements are equal to the macroscopic ones (see also the work of E. and Ming, 2004). The validity and failure of the Cauchy–Born rule for a two-dimensional mass-spring lattice was studied by Friesecke and Theil (2002). The examples of failure of this rule underline the necessity of a mathematical justification of its limits of validity.

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

638

Several one-dimensional non-linear problems were studied by Braides et al. (1999), as well as Truskinovsky (1996). Friesecke and James (2000) studied the passage from atomistic to continuum for thin films. Our work addresses the following novel aspects. First, we rigorously establish the compactness of the solutions of the discrete problems and therefore justify the continuum limit for a non-periodic three-dimensional spacial array of point masses. Second, our proof is based on the discrete Korn’s inequality (3.5), where u ðxÞ is a vector-function defined in the domain O as a linear interpolation of the Laplace transform in t of the discrete vectorfunction ui ðtÞ; t40 which solves (1.6)–(1.8). We prove this inequality for generic arrays of particles (satisfying the so-called non-periodic triangulization condition introduced in Section 2) and obtain an upper bound on the Korn’s constant in terms of the geometric characteristics of the network c1 ; c2 and the physical parameter k1 which are defined in (1.1)–(1.2) and (1.5) respectively. In particular, this bound allows us to see when this constant goes to infinity, which may lead to instabilities analogous to those discussed in the work of Friesecke and Theil (2002). We derive and justify the following continuum limit: 8 3 X > q2 uðx; tÞ q > > > rðxÞ  anpqr ðxÞnp ½uðx; tÞer ¼ 0; x 2 O; t40; > 2 > qt qx q > n;p;q;r¼1 < uðx; tÞ ¼ 0; x 2 qO; tX0; >  > > > quðx; tÞ > > uðx; 0Þ ¼ 0; ¼ aðxÞ; > : qt 

x 2 O;

t¼0

where er ðr ¼ 1; 3Þ form an orthonormal basis in R3 and np ½u ¼ 1=2ðqun =qxp þ qup =qxn Þ is the strain tensor. The functions anpqr ðxÞ; rðxÞ and the vector-function aðxÞ are determined in Sections 2 and 4 (formulas ((2.7), (4.1) and (4.2)). For generic non-periodic arrays of particles the elastic coefficients anpqr ðxÞ are determined via mesocharacteristics introduced in Section 2. For periodic lattices, coefficients anpqr ðxÞ are constants and are calculated explicitly in terms of the spring constant for a cubic lattice (Section 6). Finally, we mention here the work of Vogelius (1991), where the continuum limit was rigorously derived for planar electrical non-periodic networks. The compactness of solutions in this problem does not required Korn’s inequality since this problem is scalar. These results were subsequently extended by Krasniansky (1997) for more general scalar problems. 2. Mesocharacteristic We now introduce a mesocharacteristic which describes the elastic interaction on a mesoscale h, where 5h5diam O. Here we assume h does not depend on , and the consecutive limit is zero: lim lim

h!0 !0

 ¼ 0. h

These assumptions reflect the physics of the problem. Indeed, while the mesoscale (the size of a representative piece of an inhomogeneous medium) does not depend on the

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

639

microscale (size of an inhomogeneity, in our case the length of the springs), the microscale is an order of magnitude smaller than the mesoscale. To this end, we assume that the array of particles is such that the points xðÞ i partition the domain O into simplices (pyramids) Pa whose edges are the springs between interacting particles. The angles of the simplices are uniformly bounded from below by a strictly positive constant which does not depend on . The edges of these simplices are not necessarily equal and the simplices are not necessarily identical. The only condition on xðÞ i is (1.1). Definition 1. The network which satisfies this condition is called a triangulized network. We will show that this partition’s condition (non-periodic triangulization) is sufficient for the discrete Korn’s inequality (3.5) to hold. We introduce the following notations: jnp ðxÞ ¼ 12ðen xp þ ep xn Þ;

cnp ðxÞ ¼ 12ðen xp  ep xn Þ.

(2.1)

Consider a point y 2 O and introduce a cube K yh of side length h centered at y, where h is the mesoscale defined above. Introduce the following functional (called the mesocharacteristic): 2   3 X X  X 1  ðÞ E K y ðvÞ ¼ hC ij ðvi  vj Þ; vi  vj i þ h2g 3 T jk jjk ðxi  yÞ ,  vi  h   2 i; j y y i j;k¼1 Kh

Kh

(2.2) P

where i; j K y means the summation over all pairs of interacting particles which are located h inside the cube K yh (the total number of particles in K yh is denoted by p). Denote by v the collection of displacement vectors v1 ; . . . ; vp so that v ¼ ðv1 ; . . . ; vp Þ and let fT jk g3j;k¼1 be an arbitrary second rank tensor with constant components such that T jk ¼ T kj . The first sum in (2.2) represents the elastic energy in K yh . The second sumPis a penalty term which represents the deviation of the vectors vi from the linear part 3j;k¼1 T jk jjk ðxðÞ i  yÞ. The third sum in (2.2) can be viewed as a linear part (differential) of a homogenized vectorfunction uðxÞ (see (3.13)), and 0ogo2 is a technical parameter. We seek the minimum of y this functional among all vectors vi which correspond to particles xðÞ i 2 K h ; i ¼ 1; . . . ; p. p The minimizing displacement vectors are denoted by fwi gi¼1 ; they exist since C ij is a positive definite operator due to (1.3): min E K y ðvÞ ¼ E K y ðwÞ; v

h

h

w ¼ ðw1 ; . . . ; wp Þ.

Next, we consider a specific set of tensors defined via the Kronecker delta notation: 3 1 T ðnpÞ ¼ 12ðen  ep þ ep  en Þ ¼ fT ðnpÞ jk ¼ 2ðdjn dkp þ djp dkn Þgj;k21 ;

n; p ¼ 1; 3.

p ðnpÞ ðnpÞ If the vectors fwnp i gi¼1 minimize the functional E K y ðvÞ in which T jk ¼ T jk , then the sets of h np vectors fwi g and fwi g satisfy the following algebraic systems, respectively (the Euler–Lagrange equations): Xi np 2g 3 np C ij ðwnp  wi ¼ h2g 3 jnp ðxðÞ (2.3) i  wj Þ þ h i  yÞ, j Ky h

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

640

Xi

C ij ðwi  wj Þ þ h2g 3 wi ¼ h2g 3

j Ky h

3 X

T jk jjk ðxðÞ i  yÞ

(2.4)

j;k¼1

P i y for all particles xðÞ j K y stands for the summation over all particles which i 2 K h , where h are located in K yh and interact (that is, are connected by the springs) with a given ðÞ particle xi . P It follows from (2.3)–(2.4) that wi ¼ 3n;p¼1 wnp i T np . Substitute this formula into (2.2). Then, min E K y ðvÞ ¼ E K y ðwÞ ¼ v

h

h

3 X

anpqr ðy; ; h; gÞT np T qr ,

(2.5)

n;p;q;r¼1

where anpqr ðy; ; h; gÞ ¼

1X qr qr hC ij ðwnp  wnp j Þ; wi  wj i 2 i; j y  i Kh X qr 2g 3 np ðÞ qr ðÞ  hwnp þh i  j ðxi  yÞ; wi  j ðxi  yÞi.

ð2:6Þ

i Ky h

These coefficients play an important role in our consideration since they define the effective elastic moduli in the continuum limit. Namely, we define a fourth rank tensor anpqr ðyÞ; y 2 O, as follows: anpqr ðy; ; h; gÞ anpqr ðy; ; h; gÞ ¼ lim lim . 3 !0 h!0 !0 h h3

anpqr ðyÞ ¼ lim lim h!0

(2.7)

We consider spatial arrays of particles such that the limits (2.7) exist for any g 2 ð0; 2Þ. In Section 6, we provide an example where the limits (2.7) are calculated explicitly. 3. Compactness of the family of solutions and Korn’s inequality We begin by converting the evolutionary problem (1.6)–(1.8) into a stationary one by ðÞ using the Laplace transform in t. We denote the Laplace transform of uðÞ i ðtÞ by ui ðlÞ, and ðÞ ðÞ by a slight abuse of notation simply write ui instead of ui ðlÞ: Z 1 ðÞ ðÞ lt ui ðlÞ ¼ ðLui ÞðlÞ ¼ uðÞ dt. i ðtÞe 0

Applying the Laplace transform to (1.6)–(1.8) we obtain ðÞ ðÞ ðÞ ðÞ 2 ðÞ mðÞ i l ui þ ruðÞ Hðu1 ; . . . ; uN Þ  mi ai ¼ 0; i

uðÞ i ðlÞ  0;

i 2 1; M.

i 2 1; N,

(3.1)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

641

Fix l40. Taking into account (1.4), we observe that the solutions of (3.1) minimize the following quadratic functional: N N X X 1 X ij ðÞ ðÞ ðÞ ðÞ ðÞ 2 ðÞ ðÞ 2 F ðu Þ ¼ hC  ðui  uðÞ Þ; u  u i þ l m ju j  2 mðÞ i i hai ; ui i, j i j i 2 i; j i¼1 i¼1 (3.2) where u ¼ In order to establish the compactness of u ðxÞ, we need a uniform in  bound on the interaction energy: X ðÞ ðÞ ðÞ hC ij ðuðÞ i  uj Þ; ui  uj ipC, 

ðÞ ðuðÞ 1 ; . . . ; uN Þ.

i; j

where C does not depend on . It follows from (3.2) that F ðu ÞpF ð0Þ ¼ 0. Next, using the Cauchy–Schwartz inequality, we obtain N N X X 1 X ij ðÞ ðÞ ðÞ ðÞ ðÞ 2 ðÞ ðÞ 2 hC  ðui  uðÞ Þ; u  u i þ l m ju j p2 mðÞ i i jhai ; ui ij j i j i 2 i; j i¼1 i¼1 !1=2 !1=2 !1=2 N N N N X X X X ðÞ ðÞ 2 ðÞ ðÞ 2 ðÞ ðÞ 2 ðÞ 2 p2 mi jai j mi jui j p2 mi jai j mðÞ i jui j i¼1

þ



i¼1

i¼1

!1=2

N 2 X mðÞ jaðÞ j2 l i¼1 i i !1=2 N X 1 X ij ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ 2 2 hC  ðui  uj Þ; ui  uj i þ l mi jui j . 2 i; j i¼1

1 X ij ðÞ ðÞ ðÞ hC  ðui  uðÞ j Þ; ui  uj i 2 2l i; j

!1=2

i¼1

¼

Hence (using (1.9)), N N X 1 X ij ðÞ 4X 4 ðÞ ðÞ ðÞ ðÞ 2 ðÞ 2 2 hC  ðui  uðÞ Þ; u  u i þ l m ju j p mðÞ i i jai j p 2 C ¼ const. j i j i 2 2 i; j l i¼1 l i¼1 (3.3) From this, we conclude that the sum   is bounded uniformly in . We now construct the linear spline u ðxÞ which corresponds to the discrete vectorfunction uðÞ i : P

u ðxÞ ¼

N X

ij ðÞ i;j hC  ðui

i uðÞ i L ðxÞ,

ðÞ uðÞ j Þ; ui

uðÞ j i

(3.4)

i¼1

where Li ðxÞ; i 2 1; N, are constructed as follows. Fix a site xðÞ i and define a linear function ðÞ Li ðxÞ such that Li ðxðÞ Þ ¼ 1. Consider all sites x that are connected to xðÞ i j i by an edge of our triangulized network and set Li ðxðÞ j Þ ¼ 0 for all such neighboring sites. Then extend the function Li ðxÞ by linearity into all simplices adjacent to the site xðÞ i . Finally, set i i ðÞ i L ðxÞ  0 outside these simplices. Clearly, L ðxk Þ ¼ dik and L ðxÞ is linear. Next, define u ðxÞ by (3.4). Then, X ðÞ X0 ðÞ 2 3 jui  uðÞ jui j2 pC 2  ku ðxÞk2H 1 ðOÞ , C 1  ku ðxÞk2H 1 ðOÞ p j j þ i; j

i

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

642

P0 where the summation i; j is taken over all pairs of sites connected by the edges of our network, and the constants C 1 and C 2 do not depend on . We now prove the following discrete Korn’s inequality. We formulate it here for the Laplace transform uðÞ i ðlÞ defined above; however, the same proof works for any vectorfunction defined on the set xðÞ i and satisfying the boundary condition (1.8). Theorem 1. Suppose that the network satisfies the triangulization condition from Definition 1, and let uðÞ i be Laplace transforms of solutions of (1.6)–(1.8). Then the following inequality holds: X ðÞ ðÞ ðÞ hC ij ðuðÞ (3.5) ku ðxÞk2H 1 ðOÞ pC  i  uj Þ; ui  uj i, i; j

where the constant C ¼ 4c31 =3k1 c2 does not depend on  (see (1.1) and (1.5)). Proof. Consider a simplex Pa of our triangulization and denote by xa its center of mass. Since the spline u ðxÞ is a linear vector-function in the interior of this simplex, we can write u ðxÞ ¼ u ðxa Þ þ ðru ðxa Þ; x  xa Þ;

x 2 Pa .

Next we separate the symmetric and antisymmetric parts in the linear term: u ðxÞ ¼ u ðxa Þ þ

3 X

fnp ½u ðxa Þjnp ðx  xa Þ þ wnp ½u ðxa Þcnp ðx  xa Þg;

x 2 Pa ,

n;p¼1

where np ½u ¼ 1=2ðqun =qxp þ qup =qxn Þ; wnp ½u ¼ 1=2ðqun =qxp  qup =qxn Þ. Denote by

(3.6) Pa i;j

ðÞ the summation over all edges xðÞ i  xj of the simplex Pa . Then, using Eqs. (3.4) and (3.6), we get a a X X ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ hC ij ðuðÞ  u Þ; u  u i ¼ hC ij ðu ðxðÞ j i i j i Þ  u ðxj ÞÞ; u ðxi Þ  u ðxj Þi i; j

¼

a X i; j 3 X

i; j

* C ij

3 X

ðÞ ðÞ np ðÞ a fnp ½u ðxa Þjnp ðxðÞ i  xj Þ þ wnp ½u ðx Þc ðxi  xj Þg;

n;p¼1

fqr ½u ðxa Þjqr ðxðÞ i

+ 

xðÞ j Þ

þ

wqr ½u ðxa Þcqr ðxðÞ i



xðÞ j Þg

q;r¼1

X

3 a X c X ðÞ ðÞ ðÞ np ½u ðxa Þqr ½u ðxa Þ ðxðÞ in  xj n Þðxip  xj p Þ  n;p;q;r¼1 i; j ðÞ ðÞ ðÞ ðxðÞ iq  xj q Þðxir  xj r Þ,

ð3:7Þ

where the constant c ¼ k1 is independent of . Hereafter we denote all constants independent of  by c, and the values xðÞ in ðn 2 1; 3Þ are the coordinates of the vector ðÞ ðÞ ðÞ xðÞ i ¼ ðxi1 ; xi2 ; xi3 Þ. The last inequality in (3.7) follows from (1.1) and (1.3). Finally, we P ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ denote by X npqr the sum ai;j ðxðÞ in  xj n Þðxip  xj p Þðxiq  xj q Þðxir  xj r Þ. Consider now a

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

set of real numbers t ¼ ftnp g3n;p¼1 such that tnp ¼ tpn ; jtj2 ¼ 3 X

X npqr tnp tqr ¼

n;p;q;r¼1

a X

3 X

i; j

n;p¼1

643

P3

2 n;p¼1 tnp 40.Then

!2 ðÞ ðÞ ðÞ ðxðÞ in  xj n Þðxip  xj p Þtnp

X0.

(3.8)

We now show that this quadratic form is positive definite by a contradiction argument. Assume that there exists a non-zero set ftnp g3n;p¼1 such that the form (3.8) is equal to zero on this set. Then 3 X

ðÞ ðÞ ðÞ ðxðÞ in  xj n Þðxip  xj p Þtnp ¼ 0

(3.9)

n;p¼1 ðÞ for each pair of indices ði; jÞ that corresponds to the edge xðÞ of simplex Pa . i  xj Equations (3.9) form a system of 6 linear equations for 6 unknowns ðt11 ; t12 ; t13 ; t22 ; t23 ; t33 Þ. Simple but tedious calculations show that the absolute value of the determinant of this system is equal to 23  64  jPa j4 , where jPa j is the volume of simplex Pa . Due to the non-degeneracy condition in Definition 1, jPa ja0. Thus the system (3.9) has only the trivial solution, which establishes a contradiction. Next we show that !2 a 3 X X ðÞ ðÞ ðÞ ðÞ F ðtÞ ¼ ðxin  xj n Þðxip  xjp Þtnp XcðÞjtj2 . (3.10) i; j

n;p¼1

Indeed, if (3.10) does not hold, then for every d40 there exists a set of numbers ftnp g3n;p¼1 t t such that 0pF ðjtj Þod. This implies that F ðjtj Þ ¼ 0, which contradicts the fact that the t system (3.9) has only the trivial solution, since j jtj j ¼ 1a0. Taking into account (1.1), we conclude that cðÞ ¼ c2 4 ; c2 40. Choosing tnp ¼ np ½u ðxa Þ in (3.10) and using (3.7), we obtain a 3 X X ðÞ ðÞ ðÞ hC ij ðuðÞ  u Þ; u  u iXk  c 2np ½u ðxa Þ3 . 1 2 j i j i i; j

n;p¼1

Since the vector-function u ðxÞ is linear inside each Pa , np ½u ðxÞ is constant in Pa and thus Z X a 3 X ðÞ ðÞ ðÞ hC ij ðuðÞ  u Þ; u  u iXc 2np ½u ðxÞ dx, i j i j Pa n;p¼1

i; j

where c ¼ 6k1 c2 =c31 . Summing up over all simplices of the triangulization of the network in ðÞ O, and taking into account that each edge xðÞ i  xj enters this sum no more than 4 times, we have X0 X ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ hC ij ðuðÞ  u Þ; u  u iX4 hC ij ðuðÞ 4 i j i j i  uj Þ; ui  uj i i; j

Xc

i; j

Z X 3 O n;p¼1

2np ½u ðxÞ dx.

ð3:11Þ

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

644

ðÞ ðÞ 1 Since u ðxðÞ i Þ ¼ ui ¼ 0 for xi 2 qO, the spline u ðxÞ 2 H 0 ðOÞ. Now we can use the standard continuum Korn’s inequality Z X 3 2 ku ðxÞkH 1 ðOÞ p2 2np ½u ðxÞ dx; u ðxÞ 2 H 10 ðOÞ. (3.12) 0

O n;p¼1

Inequalities (3.11) and (3.12) imply (3.5). Theorem 1 is proved.

&

Finally, we establish the compactness of the family fu ðxÞg. Since u ðxÞ 2 H 10 ðOÞ, inequalities (3.3), (3.5) imply that the sequence of splines u ðxÞ is bounded uniformly in  in the space H 10 ðOÞ: ku ðxÞk2H 1 ðOÞ pconst. 0

Hence it is weakly compact, and due to the standard embedding theorem there exists a subsequence of the sequence u ðxÞ (which we denote for convenience u ðxÞ) such that u ðxÞ*uðxÞ weakly in H 10 ðOÞ;

u ðxÞ ! uðxÞ strongly in L2 ðOÞ;

uðxÞ 2 H 10 ðOÞ. (3.13)

As we shall show in the next section, the limiting vector-function uðxÞ will be a solution of the homogenized problem. 4. Minimizer of the homogenized problem In order to describe the homogenized problem, we introduce functions rðxÞ40 and aðxÞ defined via the Voronoi cell U i . Recall that for a given set of sites xi in the domain O, the Voronoi tessellation is defined via the Voronoi cells of the points xi . The Voronoi cell U i of a point xi is defined as U i ¼ fx 2 O : jxi  xjpjxj  xj; jaig. It is known that the cells of the Voronoi tesselation are polyhedra, and thus we obtain a partition of the domain O into polyhedra U i . If wU i ðxÞ is the characteristic function of U i , P mðÞ i we introduce the distributed density function r ðxÞ ¼ N i¼1 jU i j wU i ðxÞ, where jU i j is the volume of U i . We suppose that the sequence of functions r ðxÞ weak-star converges in L1 ðOÞ to a function rðxÞ40 as  ! 0: r ðxÞ*rðxÞ ðweak-star in L1 ðOÞÞ.

(4.1) PN

ðÞ i¼1 ai wU i ðxÞ

Next, we introduce the distributed velocity vector-function a ðxÞ ¼ by assuming that the sequence of vector-functions a ðxÞ converges strongly in L2 ðOÞ to a vector-function aðxÞ as  ! 0: L2 ðOÞ

a ðxÞ ! aðxÞ.

(4.2)

We further assume that rðxÞ and aðxÞ are smooth; otherwise, a standard approximation by smooth functions can be used. Note that the existence of limits (4.1)–(4.2) is a rather general restriction on the spatial distributions of the locations of point masses and their initial velocities. Since we do not consider any spatial periodicity, it is necessary to impose some conditions on these spatial distributions.

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

645

Theorem 2. Suppose that uðxÞ is the limiting vector-function defined in (3.13), the array of particles satisfies the non-periodic triangulization condition in the sense of Definition 1, and the limits (2.7), (4.1), (4.2) exist. Then uðxÞ minimizes the limiting (homogenized) functional ) Z ( X 3 2 2 F0 ðvÞ ¼ anpqr ðxÞnp ½vqr ½v þ l rðxÞjvðxÞj  2rðxÞhaðxÞ; vðxÞi dx (4.3) O

n;p;q;r¼1

in the class of vector-functions from H 10 ðOÞ. We first briefly outline the scheme of the proof. Our goal is to establish the following inequality: F0 ðuÞpF0 ðwÞ;

8wðxÞ 2 H 10 ðOÞ.

(4.4)

We do so first for smooth functions w 2 C 20 ðOÞ, and then we use a standard smoothing procedure. Since we use the variational duality approach, the proof consists of two steps: (i) Establishing an upper bound, and (ii) establishing a lower bound. In the first step, we prove the following inequality: lim F ½u pF0 ½w;

!0

8w 2 C 20 ðOÞ,

(4.5)

where the minimizer u ¼ u ðxðÞ i Þ is defined by (3.4). The key point of this step is the construction of a discrete vector-function wh for any wðxÞ 2 C 20 ðOÞ. This function is called the quasiminimizer because if w ¼ u (where u is defined by (3.13)), then the resulting function uh provides the almost-minimum energy to the functional F . This function is constructed by a mesopartition of the domain O into mesocubes of size h ð5h51Þ. The fact that u minimizes the functional F among all discrete vector-functions implies that F ½u pF ½wh . The vector-function wh has the following two properties. First, the function wh ¼

N X

wh ðxðÞ i ÞwU i ðxÞ

i¼1

approximates wðxÞ in the following sense: lim lim kwh ðxÞ  wðxÞkL2 ðOÞ ¼ 0.

h!0 !0

Furthermore, the explicit construction of the vector-function wh allows us to pass to the limit in F ½wh  and obtain the following inequality: lim lim F ½wh pF0 ½w,

h!0 !0

(4.6)

where F0 is the homogenized functional (4.3). Since the construction of wh is the key point of the proof, we now explain the idea behind this construction. In each cube of the mesopartition, this discrete vector-function

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

646

has the form wxha ðxðÞ i Þ ¼ wðxa Þ þ

3 X

ðÞ np ½wðxa Þvnp a ðxi Þ þ

n;p¼1 ðÞ vectors vnp a ðxi Þ

3 X

wnp ½wðxa Þcnp ðxðÞ i  xa Þ,

(4.7)

n;p¼1

where the constant ði 2 1; pÞ minimize the functional E ðnpÞ ðvÞ (see (2.2) with K xh a ðÞ np np ðÞ T ¼ T ). To explain (4.7), we note that if vnp ðx Þ is replaced by j ðx  x Þ, then it becomes a a i i the linear part of the Taylor expansion of wðxÞ in the mesocube K xh a centered at point xa . The vector-function (4.7) possesses the following properties. First, it approximates (in some sense) the vector-function wðxÞ in the cube K xh a . Second, it almost minimizes the local energy functional (mesocharacteristic) (2.2) when T jk ¼ jk ½wðxa Þ. Indeed, the first and the second terms in (4.7) do not contribute to the elastic energy. The latter property allows us to pass to the limit as  ! 0 and compute the limiting functional F0 via mesocharacteristics. If K xh a is a partition of the domain O, then a global quasiminimizer can be constructed, roughly speaking, by gluing together the quasiminimizers of (2.2) with the help of an appropriate partition of unity. In the second step we establish the following inequality: (4.8) F0 ½up lim F ½u . !0

Clearly, (4.8) and (4.5) imply (4.4). Proof of Theorem 2. (i) The upper bound. Take any vector-function wðxÞ 2 C 20 ðOÞ. Choose a mesoscale parameter h such that there S exists a cover of the domain O by cubes K xh a centered at points xa : O  a2L K xh a . These g

points form a cubic lattice of period h  h1þ2 ; 0ogo0 such that the cubes overlap. Due to the overlap of the cubes, we can further select smaller cubes K x0a with the same center xa h g 0 and edges of length h ¼ h  2h1þ2 . It is well known that there exists a set of functions ffa ðxÞ 2 C 1 0 ðOÞga2L (called a partition of unity) such that ( X 1; x 2 K xh a ; c fa ðxÞ ¼ fa ðxÞ  1; x 2 O. 0pfa ðxÞp1; jrfa ðxÞjp 1þg ; xa 0; xeK h ; h 2 a2L We now construct a discrete vector-function ( 3 X X ðÞ ðÞ wh ðxi Þ ¼ wðxa Þ þ np ½wðxa Þvnp a ðxi Þ a2L

þ

3 X

n;p¼1

wnp ½wðxa Þcnp ðxðÞ i  xa Þ

(4.9)

)  fa ðxðÞ i Þ,

ð4:10Þ

n;p¼1 ðÞ ðnpÞ where vnp a ðxi Þ are the minimizers of the functional E K xh a ðvÞ defined by (2.2). Taking into account (2.6) and (2.7), we estimate the interaction energy in cube K xh a : X ðÞ 3 np ðÞ qr ðÞ qr ðÞ hC ij ðvnp a ðxi Þ  va ðxj ÞÞ; va ðxi Þ  va ðxj Þi ¼ Oðh Þ, i; j K xa h

3

X i

K xh a

ðÞ 5þg np ðÞ qr ðÞ qr ðÞ hvnp Þ. a ðxi Þ  j ðxi  xa Þ; va ðxi Þ  j ðxi  xa Þi ¼ Oðh

(4.11)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

647

Next, we estimate the energy in the overlapping parts K xh a nK x0a as h ! 0. To this end, we h subtract the energy of the particles in the smaller cube K x0a from the energy of the particles h in the larger cube K xh a . This leads to an upper bound of the energy functional because some particles in the cube K x0a may interact with the particles in the cube K xh a (i.e., the h corresponding springs cross the boundary qK x0a ): h X ðÞ np ðÞ qr ðÞ qr ðÞ hC ij ðvnp ðx Þ  v ðx ÞÞ; v ðx a a a i j i Þ  va ðxj Þi i; j K xa nK xa 0 h h

þ h2g 3

X i

p

8
ðÞ np ðÞ qr ðÞ qr ðÞ hC ij ðvnp a ðxi Þ  va ðxj ÞÞ; va ðxi Þ  va ðxj Þi

i; j K xa h

þh2g 3



ðÞ np ðÞ qr ðÞ qr ðÞ hvnp a ðxi Þ  j ðxi  xa Þ; va ðxi Þ  j ðxi  xa Þi

K xh a nK x0a h

X



jnp ðxðÞ i



ðÞ xa Þ; vqr a ðxi Þ



jqr ðxðÞ i

i K xa h

8 > :

ðÞ hvnp a ðxi Þ

9 =  xa Þi ;

ðÞ np ðÞ qr ðÞ qr ðÞ hC ij ðvnp a ðxi Þ  va ðxj ÞÞ; va ðxi Þ  va ðxj Þi

i; j K xa 0 h

0

þ ðh Þ2g 3

X i

ðÞ np ðÞ qr ðÞ qr ðÞ hvnp a ðxi Þ  j ðxi  xa Þ; va ðxi Þ  j ðxi  xa Þi

K x0a h

0

þ ððh Þ2g  h2g Þ3

X

9 = ;

ðÞ np ðÞ qr ðÞ hvnp a ðxi Þ  j ðxi  xa Þ; va ðxi Þ

i K x0a h

 jqr ðxðÞ i  xa Þig

 1 1  ¼ anpqr ðxa ; ; h; gÞ  anpqr ðxa ; ; h ; gÞ þ Oððh Þ Þ 0 ðh Þ2þg h2þg 0 1 0 !2þg !3 1 0 0 0 0 h A ¼ Oðh3 Þ@1  h A ¼ Oðh3  ðh Þ3 Þ þ Oððh Þ3 Þ@1  h h 0 1 !2þg 0 g 0 h 3 @ A ¼ Oðh3 Þð1  ðð1  2h2 Þ3 Þ þ Oððh Þ Þ 1  h 0

0

g

h!0

0

5þg



g

þ Oððh Þ3 Þð1  ð1  2h2 Þ2þg Þ p Oðh3 Þh2 ¼ oðh3 Þ. In the first equality we used (2.7) and (4.11); in the second equality we used (2.7). Thus X ðÞ 3 np ðÞ qr ðÞ qr ðÞ hC ij ðvnp a ðxi Þ  va ðxj ÞÞ; va ðxi Þ  va ðxj i ¼ oðh Þ, i; j K xa nK xa 0 h



3

X

h

ðÞ 5þg np ðÞ qr ðÞ qr ðÞ hvnp Þ. a ðxi Þ  j ðxi  xa Þ; va ðxi Þ  j ðxi  xa Þi ¼ oðh

i K xa nK x0a h h

(4.12)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

648

Since w 2 C 20 ðOÞ, it can be written in the form 3 X

wðxÞ ¼ wðxa Þ þ

fnp ½wðxa Þjnp ðx  xa Þ þ wnp ½u ðxa Þcnp ðx  xa Þg þ ga ðxÞ,

(4.13)

n;p¼1

where ga ðxÞ ¼ Oðh2 Þ. Substituting (4.13) into (4.10) we obtain ( 3 X X ðÞ ðÞ np ðÞ wh ðxi Þ ¼ wðxðÞ Þ þ np ½wðxa Þ  ðvnp a ðxi Þ  j ðxi  xa ÞÞ i a2L

n;p¼1

)

ga ðxðÞ i Þ

 fa ðxðÞ i Þ.

ð4:14Þ

This allows us to evaluate F ðwh Þ: F ðwh Þ ¼

1 X ij ðÞ ðÞ ðÞ hC  ðwh ðxðÞ i Þ  wh ðxj ÞÞ; wh ðxi Þ  wh ðxj Þi 2 i; j þ l2

N X

ðÞ 2 mðÞ i jwh ðxi Þj  2

i¼1

N X

ðÞ ðÞ mðÞ i hai ; wh ðxi Þi.

ð4:15Þ

i¼1

We now use the estimates (4.11) and (4.12) in order to obtain the upper bound (4.6). Consider the first term of (4.15): 1 X ij ðÞ ðÞ ðÞ hC  ðwh ðxðÞ i Þ  wh ðxj ÞÞ; wh ðxi Þ  wh ðxj Þi 2 i; j 1XX ðÞ ðÞ ðÞ hC ij ðwh ðxðÞ ¼ i Þ  wh ðxj ÞÞ; wh ðxi Þ  wh ðxj Þi þ Ah 2 a2L i; j xa  K

h

0

¼ H h þ Ah ,

ð4:16Þ

ðÞ xa where Ah contains terms that correspond to neighboring particles xðÞ i 2 K h0 ; xj 2 ðÞ ðÞ xa xa xa xa K h nK 0 or xi ; xj 2 K h nK 0 . Substituting (4.10) in H h , we get h

h

1XX 2 a2L i; j

*

( C ij

K x0a

3 X

ðÞ np ðÞ np ½wðxa Þ  ðvnp a ðxi Þ  va ðxj ÞÞ

n;p¼1

h

þ

3 X

) wnp ½wðxa Þ  c

np

ðxðÞ i



xðÞ j Þ

,

n;p¼1 3 X q;r¼1

qr ½wðxa Þ 

ðÞ ðvqr a ðxi Þ



ðÞ vqr a ðxj ÞÞ

þ

3 X q;r¼1

+ np

wqr ½wðxa Þ  c

ðxðÞ i



xðÞ j Þ

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

¼

3 X X

np ½wðxa Þ  qr ½wðxa Þ

a2L n;p;q;r¼1

0

B1 X @ 2 i; j

649

K x0a

1 ðÞ np ðÞ qr ðÞ qr ðÞ C hC ij ðvnp a ðxi Þ  va ðxj ÞÞ; va ðxi Þ  va ðxj ÞiA

h

3 1X X

np ½wðxa Þ  wqr ½wðxa Þ 2 a2L n;p;q;r¼1 X ðÞ ðÞ qr ðÞ np ðÞ hC ij ðvnp  a ðxi Þ  va ðxj ÞÞ; c ðxi  xj Þi þ

i; j K xa 0 h

3 1X X

wnp ½wðxa Þ  qr ½wðxa Þ 2 a2L n;p;q;r¼1 X ðÞ qr ðÞ qr ðÞ hC ij cnp ðxðÞ  i  xj Þ; va ðxi Þ  va ðxj Þi þ

i; j K xa 0 h

3 1X X þ wnp ½wðxa Þ  wqr ½wðxa Þ 2 a2L n;p;q;r¼1 X ðÞ ðÞ qr ðÞ hC ij cnp ðxðÞ  i  xj Þ; c ðxi  xj Þi. i; j K xa 0 h

Taking into account (1.4) and the definition of cnp ðxÞ (see (2.1)), we see that ðÞ ðÞ qr ðÞ np ðÞ hC ij ðvnp a ðxi Þ  va ðxj ÞÞ; c ðxi  xj Þi ðÞ qr ðÞ qr ðÞ ¼ hC ij cnp ðxðÞ i  xj Þ; va ðxi Þ  va ðxj Þi ðÞ ðÞ qr ðÞ ¼ hC ij cnp ðxðÞ i  xj Þ; c ðxi  xj Þi ¼ 0.

The expression H h in (4.16) can then be bounded from above by 8 > 3 <1 X X X ðÞ np ðÞ qr ðÞ np ½wðxa Þ  qr ½wðxa Þ  hC ij ðvnp a ðxi Þ  va ðxj ÞÞ; va ðxi Þ > 2 : i; j xa a2L n;p;q;r¼1 K

ðÞ 2g 3  vqr  a ðxj Þi þ h

X

h

0

ðÞ np ðÞ qr ðÞ qr ðÞ hvnp a ðxi Þ  j ðxi  xa Þ; va ðxi Þ  j ðxj  xa Þi

i K x0a h

9 > = > ;

which, due to (2.6) and (3.11), gives in the limit as  ! 0 and h ! 0 Z

3 X

np ½wðxÞ  qr ½wðxÞ  anpqr ðxÞ dx.

(4.17)

O n;p;q;r¼1

We now prove that Ah vanishes as  ! 0 and h ! 0. To this end we consider the case ðÞ xa xa when the particles xðÞ i and xj are located in K h nK 0 . Then, making use of (4.9) and (4.14), h

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

650

we have XX

ðÞ ðÞ ðÞ hC ij ðwh ðxðÞ i Þ  wh ðxj ÞÞ; wh ðxi Þ  wh ðxj Þi

a2L i; j K xa nK xa 0 h h

¼

*

XX

" ðÞ wðxðÞ i Þ  wðxj Þ

C ij

a2L i; j K xa nK xa 0 h h

þ

3 X X

ðÞ np ðÞ np ½wðxb Þ  fðvnp b ðxi Þ  j ðxi  xb Þ

b2L n;p¼1

gb ðxðÞ i ÞÞ



# fb ðxðÞ i Þ

ðÞ wðxðÞ i Þ  wðxj Þ þ



ðÞ ðvnp b ðxj Þ

3 XX

np

j

ðxðÞ j

 xb Þ 

gb ðxðÞ j ÞÞ



fb ðxðÞ j Þg

ðÞ np ðÞ np ½wðxg Þ  fðvnp g ðxi Þ  j ðxi  xg Þ

g2L n;p¼1

gg ðxðÞ i ÞÞ



fg ðxðÞ i Þ



,

+

ðÞ ðvnp g ðxj Þ

np

j

ðxðÞ j

 xg Þ 

gg ðxðÞ j ÞÞ



fg ðxðÞ j Þg

.

ð4:18Þ

Next, using (1.4) and the fact that for a smooth vector-function wðxÞ the inequality kwðxÞ  wðyÞkR3 pckx  ykR3 holds, we get      X X   ðÞ ðÞ ðÞ ðÞ hC ij ðwðxi Þ  wðxj ÞÞ; wðxi Þ  wðxj Þi    a2L i; j xa xa   K h nK 0 h XX ðÞ ðÞ pc kij   jwðxi Þ  wðxj Þj2 a2L i; j K xa nK xa 0 h h

pc

XX

3 pc

a2L i; j K xa nK xa 0 h

pc

X a2L

X h2 ðh  h0 Þ a2L

jOj

 N  3

h

g

h3þ2 pc

g h!0 jOj 3þ2g  h pch2 ! 0. 3 h

Straightforward calculations show that ðÞ ðÞ ðÞ np ðÞ np ðÞ np ðÞ ðvnp g ðxi Þ  j ðxi  xg ÞÞ  fg ðxi Þ  ðvg ðxj Þ  j ðxj  xg ÞÞ  fg ðxj Þ ðÞ ðÞ ðÞ np ðÞ np ðÞ ¼ ðvnp g ðxi Þ  vg ðxj Þ  j ðxi  xj ÞÞ  fg ðxi Þ ðÞ ðÞ ðÞ np ðÞ þ ðvnp g ðxj Þ  j ðxj  xg ÞÞ  ðfg ðxi Þ  fg ðxj ÞÞ.

Next, we use this equality to estimate some of the terms in (4.18):   * X X a X 3 X  ðÞ ðÞ np ðÞ C ij ðwðxðÞ Þ  wðx ÞÞ; np ½wðxg Þ  ðvnp  g ðxi Þ  vg ðxj Þ i j  a2L i; j xa xa g2L n;p¼1  K h nK 0 h

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

651

 +  3 XX X  ðÞ ðÞ ðÞ pc jnp ðxðÞ  x ÞÞ  f ðx Þ fhC ij ðwðxðÞ  g i i j i Þ  wðxj ÞÞ,  a2L i; j K xa nK xa n;p¼1  0 h h

ðÞ ðÞ ij 2 np ðÞ np ðÞ vnp g ðxi Þ  vg ðxj Þi þ k   jj ðxi  xj Þjg 80 11=2 > 3 < XX X B ðÞ ðÞ ðÞ C pc @ hC ij ðwðxðÞ i Þ  wðxj ÞÞ; wðxi Þ  wðxj ÞiA > : a2L i; j K xa nK xa n;p¼1 h

h

0

BX X @

0

3 X

11=2 ðÞ np ðÞ np ðÞ np ðÞ C hC ij ðvnp g ðxi Þ  vg ðxj ÞÞ; vg ðxi Þ  vg ðxj ÞiA

a2L i; j K xa nK xa n;p¼1 0 h

9 > =

h

XX

8 > <

X g kij 3 pc h4  oðh3 Þ > > ; : a2L i; j K xa nK xa a2L 0 h h ( )  1=2 g g jOj h!0 pc h4   oðh3 Þ þ h2 ! 0. 3 h þ

!1=2 þ

XX

kij 3

a2L i; j K xa nK xa 0 h h

9 > = > ;

P Here we use the fact that the symbol ag2L stands for the summation over the cubes that contain the slab K xh a nK x0a . The number of such cubes is bounded by a constant c which does h not depend on h. In the first inequality we used (1.4); in the second we used (4.12). Next, we estimate the remaining terms in (4.18) using the Cauchy–Schwartz inequality, (1.4), (4.9) and (4.12):   * X X a X 3 X  ðÞ ðÞ np ðÞ C ij ðwðxðÞ Þ  wðx ÞÞ; np ½wðxg Þ  ðvnp  g ðxj Þ  j ðxj  xg ÞÞ i j  a2L i; j xa xa g2L n;p¼1  K h nK 0 h  + ðÞ  X 3 X jxðÞ  i  xj j ðÞ ðÞ ij 2 np ðÞ np ðÞ ðfg ðxi Þ  fg ðxj ÞÞ p k  jvg ðxj Þ  j ðxj  xa Þj  g  i; j xa xa h1þ2 n;p¼1  K h nK 0 h 0 11=2 0 3 X c X BX C BX ðÞ p 1þg kij 3 A  @ kij 3  jvnp @ a ðxj Þ 2 h a2L i; j K xa nK xa i; j K xa nK xa n;p¼1 h

h

0

h

11=2 2C  jnp ðxðÞ j  xa Þj A

¼c

X a2L

g

oðh3þ4 Þpc

p

c X h

1þ2g

a2L

jOj 3þ4g h!0  h ! 0. h3

3

g

h

0

5

g

h2þ4  oðh2þ2 Þ

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

652

Similarly,   * ! X X a X 3 X  np ðÞ np ðÞ ðÞ ij C np ½wðxb Þ  ðvb ðxi Þ  vb ðxj ÞÞ  fb ðxi Þ ;   a2L i; j xa xa n;p¼1 b2L  K h nK 0 h  +  a X 3 X  ðÞ ðÞ ðÞ qr ðÞ qr ½wðxg Þ  ðvqr g ðxj Þ  j ðxj  xg ÞÞ  ðfg ðxi Þ  fg ðxj ÞÞ   g2L q;r¼1  * 3 XX X ðÞ np ðÞ pc C ij ð ðvnp b ðxi Þ  vb ðxj ÞÞ; a2L i; j K xa nK xa 0 h

n;p¼1

h

a X 3 X

+

ðÞ ðvqr g ðxj Þ



jqr ðxðÞ j

 xg ÞÞ



g2L q;r¼1

 h

1þ2g

0

11=2 3 X X X c B ðÞ np ðÞ np ðÞ np ðÞ C p 1þg @ hC ij ðvnp g ðxi Þ  vg ðxj ÞÞ; vg ðxi Þ  vg ðxj ÞiA 2 h a2L i; j K xa nK xa n;p¼1 h

0

BX X @

h

0

3

a2L i; j K xa nK xa 0 h

11=2

3 X q;r¼1

ðÞ qr ðÞ 2C jvqr a ðxj Þ  j ðxj  xa Þj A

h

p

X

c g

h1þ2

!1=2 3

oðh Þ

a2L



X

!1=2 oðh

5þg

Þ

a2L



1=2  5þg 1=2  3 1=2  5þg 1=2 oðh3 Þ oðh Þ oðh Þ oðh Þ h!0 p 1þg   ¼c  ! 0 5þg 3 3 3 2 h h h h h c

and   * X X 3 XX  ðÞ ðÞ np ðÞ C ij np ½wðxb Þ  ðvnp  b ðxj Þ  j ðxj  xb ÞÞ  ðfb ðxi Þ  a2L i; j xa xa n;p¼1 b2L  K h nK 0

 +  X  ðÞ ðÞ qr ðÞ qr ðÞ fb ðxðÞ ÞÞ ;  ½wðx Þ  ðv ðx Þ  j ðx  x ÞÞ  ðf ðx Þ  f ðx ÞÞ  qr g g g i g j g j j j  g2L q;r¼1  h

!

pc

3 X

XX



a2L i; j K xa nK xa 0 h

3 XX

ðÞ np ðÞ jvnp b ðxj Þ  j ðxj  xb Þj

b2L n;p¼1

h



3 XX g2L q;r¼1

ðÞ qr ðÞ jvqr g ðxj Þ  j ðxj  xg Þj 

2 h2þg

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

0

c BX X p 2þg @ h a2L i; j

3 

BX X @

3 

a2L i; j K xa nK xa 0 h

"

c h

ðÞ np ðÞ 2C jvnp a ðxj Þ  j ðxj  xa Þj A

h

0

p

11=2

3 X n;p¼1

K xh a nK x0a

h

X

2þg

q;r¼1

#

oðh5þg Þ p

a2L

3 X

c h

2þg

653

11=2 ðÞ qr ðÞ 2C jvqr a ðxj Þ  j ðxj  xa Þj A



jOj oðh5þg Þ h!0  oðh5þg Þpc  5þg ! 0. 3 h h

ð4:19Þ

We now rewrite the terms in (4.18) which contain the vector-function ga ðxÞ: ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ ga ðxðÞ j Þfa ðxj Þ  ga ðxi Þfa ðxj Þ ¼ ðga ðxj Þ  ga ðxi ÞÞ  fa ðxj Þ ðÞ ðÞ þ ga ðxðÞ j Þ  ðfa ðxj Þ  fa ðxj ÞÞ.

Then, analogous to the derivation of (4.19), we get XX ðÞ ðÞ ðÞ lim lim hC ij ðwh ðxðÞ i Þ  wh ðxj ÞÞ; wh ðxi Þ  wh ðxj Þi ¼ 0. h!0 !0

a2L i; j K xa nK xa 0 h h

ðÞ xa xa xa The case when xðÞ i 2 K h0 and xj 2 K h nK h0 can be considered similarly. Thus we have shown that the second term in the RHS of (4.18) vanishes:

lim lim Ah ¼ 0.

(4.20)

h!0 !0

We now show that the vector-functions wh ðxÞ ¼ vector-function wðxÞ in L2 ðOÞ. Indeed,

PN

ðÞ i¼1 wh ðxi Þ

 wU i ðxÞ converge to the

2   X N   ðÞ kwh ðxÞ  ¼ wh ðxi Þ  wU i ðxÞ  wðxÞ   i¼1 L2 ðOÞ  " ( X N 3 X X  ðÞ ¼ wðxðÞ Þþ np ½wðxa Þ  ðvnp a ðxi Þ i  i¼1 a2L n;p¼1 2 ) #   ðÞ ðÞ np ðÞ j ðxi  xa ÞÞ  ga ðxi Þ  fa ðxi Þ  wU i ðxÞ  wðxÞ  L2 ðOÞ 2   " (  X X N N X X 3    ðÞ p2 wðxðÞ þ 2  ½wðxa Þ  ðvnp a ðxi Þ i Þ  wU i ðxÞ  wðxÞ   i¼1  i¼1 a2L n;p¼1 np L2 ðOÞ 2 ) #   ðÞ ðÞ np ðÞ j ðxi  xa ÞÞ  ga ðxi Þ  fa ðxi Þ  wU i ðxÞ  wðxÞk2L2 ðOÞ

L2 ðOÞ

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

654

 (  3 X X X ðÞ ðÞ 2  p2 kwðxi Þ  wðxÞkL2 ðU i Þ þ 2 np ½wðxa Þ  ðvnp a ðxi Þ x  a2L i K a n;p¼1 i¼1 h 2 )  N X  ðÞ 2 jnp ðxðÞ  x ÞÞ  g ðx Þ  w ðxÞ p2c max jx  xðÞ  a U a i i i j  jU i j i  x2U i i¼1 L2 ðOÞ  ( )  N X X 3 X X ðÞ ðÞ ðÞ np  þ2 np ½wðxa Þ  ðvnp a ðxi Þ  j ðxi  xa ÞÞ  ga ðxi Þ  j¼1  a2L i K xa n;p¼1 h 2   N X 3 X 2c 2 3    wU i ðxÞ p 3     þ 2c np ½wðxa Þ     j¼1 n;p¼1 N X

L2 ðU j Þ

ðÞ ðvnp aj ðxj Þ



np

j

ðxðÞ j

 xaj ÞÞ 

( N 3 X X

p2c   þ 2c 2

j¼1

ðÞ kvnp aj ðxj Þ

2 

 gaj ðxðÞ j Þ  j

np

L2 ðU j Þ

ðxðÞ j



)

xaj Þk2L2 ðU j Þ

þ

2 kgaj ðxðÞ j ÞkL2 ðU j Þ

, ð4:21Þ

n;p¼1

where aj is the index of a cube which contains the jth particle. Clearly, the term 2c2 vanishes as  ! 0. We next show that the second term in the RHS of the last inequality in (4.21) also vanishes as  ! 0: ( ) N 3 X X ðÞ 2 np ðÞ np ðÞ 2 kvaj ðxj Þ  j ðxj  xaj ÞkL2 ðU j Þ þ kgaj ðxj ÞkL2 ðU j Þ j¼1

¼

n;p¼1

( N 3 X X j¼1

pc

3 X n;p¼1

np

j

ðxðÞ j

2

 xaj Þj þ

2 jgaj ðxðÞ j Þj

n;p¼1

( N 3 X X j¼1

pc

) ðÞ jvnp aj ðxj Þ

 jU i j

) ðÞ jvnp aj ðxj Þ

j

np

ðxðÞ j

2

4

 xaj Þj þ Oðh Þ

 3

n;p¼1

3

N X

ðÞ 4 np ðÞ 2 3 jvnp aj ðxj Þ  j ðxj  xaj Þj þ c  N  Oðh Þ   .

j¼1

Since N ¼ 3 , the term c  N  Oðh4 Þ  3 vanishes as h ! 0. Next we observe that (4.11) implies 3

N X

ðÞ np ðÞ 2 3 jvnp aj ðxj Þ  j ðxj  xaj Þj p

XX a2L

j¼1

¼

X a2L

ðÞ np ðÞ 2 jvnp a ðxj Þ  j ðxj  xa Þj

j K xa h

Oðh5þg Þpc

jOj h!0  Oðh5þg Þ ! 0. 3 h

Thus we have shown that lim lim kwh ðxÞ  wðxÞkL2 ðOÞ ¼ 0.

h!0 !0

(4.22)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

655

We now pass to the limit as  ! 0 and h ! 0 in the second term of (4.15). We have Z N N X X ðÞ ðÞ ðÞ ðÞ ðÞ mðÞ ha ; w ðx Þi ¼ r ðx Þha ðx Þ; w ðx Þi  jU j ¼ r ðxÞha ðxÞ; wh ðxÞi dx. h  h i  i i i i i i i¼1

O

i¼1

Now using (4.1), (4.2), (4.22) and the last equality, we obtain Z N X ðÞ ðÞ lim lim mðÞ ha ; w ðx Þi ¼ rðxÞhaðxÞ; wðxÞi dx. h i i i h!0 !0

Analogously, Z N X ðÞ 2 mðÞ jw ðx Þj ¼ rðxÞjwðxÞj2 dx. lim lim h i i h!0 !0

(4.23)

O

i¼1

(4.24)

O

i¼1

Combining (4.15)–(4.17), (4.20), (4.23) and (4.24) we obtain the upper bound (4.6). (ii) We now establish the lower bound (4.8). Assume for simplicity that the vector-function uðxÞ defined in (3.13) belongs to the class C 20 ðOÞ. We partition the domain O by cubes K xh a whose centers form a cubic lattice of period h, and thus the cubes do not overlap: [ x O K xh a ; K xh a \ K h b ¼ ; aab. a2L

Note that these cubes partition the domain O, as opposed to covering it as in (i). We now construct another quasiminimizer wa ðxÞ in each cube K xh a which almost minimizes the interaction (first) term in (3.2) as follows: wa ðxÞ ¼ u ðxÞ  uðxa Þ 

3 X

wnp ½uðxa Þcnp ðx  xa Þ.

n;p¼1

Then 3 X

min E K xa ðvÞ ¼ v

h

anpqr ðxa ; ; h; gÞT np T qr pE K xa ðwa Þ. h

n;p;q;r¼1

Choosing T np ¼ np ½uðxa Þ, we get 3 X

anpqr ðxa ; ; h; gÞnp ½uðxa Þqr ½uðxa Þ

n;p;q;r¼1

p

1X 2 i; j

a ðÞ a ðÞ a ðÞ hC ij ðwa ðxðÞ i Þ  w ðxj ÞÞ; w ðxi Þ  w ðxj Þi K xh a

þ h2g 3

X

2    3 X   a ðÞ ðÞ np ½uðxa Þjnp ðxi  xa Þ . w ðxi Þ    xa n;p¼1

i K h

Using (1.3) and the obvious equality ðÞ ðÞ a ðÞ wa ðxðÞ i Þ  w ðxj Þ ¼ u ðxi Þ  u ðxj Þ 

3 X n;p¼1

ðÞ wnp ½uðxa Þcnp ðxðÞ i  xj Þ,

ð4:25Þ

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

656

we obtain ðÞ ðÞ ðÞ ðÞ ij a ðÞ a ðÞ a ðÞ hC ij ðwa ðxðÞ i Þ  w ðxj ÞÞ; w ðxi Þ  w ðxj Þi ¼ hC  ðu ðxi Þ  u ðxj ÞÞ; u ðxi Þ  u ðxj Þi.

We now evaluate the second term of the RHS of inequality (4.25): 2    3 X   a ðÞ np ðÞ np ½uðxa Þj ðxi  xa Þ w ðxi Þ    n;p¼1   3 X  ¼ u ðxðÞ Þ  uðx Þ  wnp ½uðxa Þcnp ðxðÞ a i i  xa Þ  n;p¼1 2  3 X  np ðÞ  np ½uðxa Þj ðxi  xa Þ  n;p¼1    ðÞ ðÞ 2 p2ju ðxðÞ i Þ  uðxi Þj þ 2uðxi Þ  uðxa Þ  3 X



np

np ½uðxa Þj

ðxðÞ i

 xa Þ 

n;p¼1

3 X

np

wnp ½uðxa Þc

ðxðÞ i

n;p¼1

2    x a Þ 

ðÞ 2 4 ¼ 2ju ðxðÞ i Þ  uðxi Þj þ Oðh Þ.

Summing up inequality (4.25) over all cubes of our partition and passing to the limit as  ! 0 and h ! 0, we obtain lim lim

h!0 !0

3 X X

anpqr ðxa ; ; h; gÞnp ½uðxa Þqr ½uðxa Þ

a2L n;p;q;r¼1

Z ¼

3 X

anpqr ðxÞnp ½uðxÞqr ½uðxÞ dx

O n;p;q;r¼1

(

p lim lim

h!0 !0

þ2h

1 X ij ðÞ ðÞ ðÞ hC  ðu ðxðÞ i Þ  u ðxj ÞÞ; u ðxi Þ  u ðxj Þi 2 i; j

2g 3

 

N X i¼1

ju ðxðÞ i Þ



2 uðxðÞ i Þj

þh

2g 3

 

N X

)

4

Oðh Þ .

i¼1

Obviously, the last term is of order h2g  Oðh4 Þ, and it vanishes in the limit as  ! 0, since 0ogo2. The second term is equal to zero in the limit as  ! 0, due to (3.13) and the following lemma: Lemma 1. If the sequence of splines u ðxÞ defined by (3.4) is bounded uniformly on  in P ðÞ H 1 ðOÞ, then the sequence of the piecewise-constant vector-functions u ðxÞ ¼ N i¼1 ui  wU i ðxÞ converges in L2 ðOÞ to the vector-function uðxÞ defined by (3.13). Proof of this lemma is presented in Appendix A.

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

657

Using Lemma 1, we now obtain 3

N X

ðÞ 2 3 ju ðxðÞ i Þ  uðxi Þj ¼ 

i¼1

pc

N X

ðÞ 2 ju ðxðÞ i Þ  u ðxi Þj

i¼1 N X

ðÞ 2 ju ðxðÞ i Þ  u ðxi Þj  jU i j

i¼1

¼c

N X

ku  u k2L2 ðU i Þ ¼ cku  u k2L2 ðOÞ p2cðku  uk2L2 ðOÞ þ ku  u k2L2 ðOÞ Þ,

i¼1

P ðÞ 2 where u ðxÞ ¼ N i¼1 uðxi Þ  wU i ðxÞ. Clearly ku  ukL2 ðOÞ ! 0 as  ! 0, due to Lemma 1. 2 For the term ku  u kL2 ðOÞ , we have N N Z X X 2 2 ku  u kL2 ðU i Þ ¼ juðxÞ  uðxðÞ i Þj dx i¼1

i¼1

pc

N X i¼1

Ui

2 max jx  xðÞ i j  jU i jpc x2U i

N X

!0

5 ¼ c2 ! 0.

i¼1

Thus we have shown that Z X 3 anpqr ðxÞnp ½uðxÞqr ½uðxÞ dx O n;p;q;r¼1

1 X ij ðÞ ðÞ ðÞ hC  ðu ðxðÞ i Þ  u ðxj ÞÞ; u ðxi Þ  u ðxj Þi. !0 2 i; j

p lim

ð4:26Þ

We now estimate the second and the third terms of functional F . We have Z N N X X ðÞ 2 ðÞ ðÞ ðÞ l2 mðÞ ju j  2 m ha ; u i ¼ ½l2 r ðxÞju ðxÞj2  2r ðxÞha ðxÞ; u ðxÞi dx. i i i i i i¼1

O

i¼1

Then, taking into consideration (4.1),(4.2) and using Lemma 1, we conclude that ( ) N N X X ðÞ ðÞ ðÞ ðÞ ðÞ lim l2 mi jui j2  2 mi hai ; ui i !0

i¼1

Z

i¼1

½l2 rðxÞjuðxÞj2  2rðxÞhaðxÞ; uðxÞi dx.

¼

ð4:27Þ

O

Combining (4.26) and (4.27), we get F0 ½up lim F ½u ,

(4.28)

!0

provided that the limiting vector-function uðxÞ is smooth ðC 20 ðOÞÞ. If uðxÞ 2 H 10 ðOÞ, we use a standard approximation procedure: since the domain O has a 1 C -boundary, for any vector-function uðxÞ 2 H 10 ðOÞ there exists a sequence of vectorfunctions ud ðxÞ 2 C 20 ðOÞ (Mizohata, 1973) such that lim kud  ukH 1 ðOÞ ¼ 0.

d!0

We now use the following:

(4.29)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

658

Lemma 2. For any vector-function ud ðxÞ 2 C 20 ðOÞ, there exists a sequence of vector-functions ud ðxÞ 2 H 10 ðOÞ of the form (3.4) such that ð1Þ lim kud  ud kL2 ðOÞ ¼ 0; !0

ð2Þ lim kud  u kH 1 ðOÞ ¼ 0; d!0

ð3Þ lim kud  ud kL2 ðOÞ ¼ 0. !0

Proof of this lemma is presented in Appendix B. It follows from (4.28) and Lemma 2 that F0 ½ud p lim F ½ud . !0

Let d ! 0. Then, taking into account (4.29) and (2) from Lemma 2, we obtain the lower bound (4.8), since the functionals F0 ½u and F ½u  are continuous in the space H 1 ðOÞ. Finally, we combine (4.6) and (4.8) to obtain F0 ½up lim F ½u p lim lim F ½wh pF0 ½w; !0

h!0 !0

8w 2 C 20 ðOÞ.

Using the smoothing approximation (4.29), we conclude that F0 ½upF0 ½w;

8w 2 H 10 ðOÞ.

5. The homogenization theorem Clearly, the limiting vector-function uðxÞ defined in (3.13) satisfies the Euler–Lagrange equations 8 3 > <  P q fanpqr ðxÞnp ½uðxÞger þ l2 rðxÞuðxÞ ¼ rðxÞaðxÞ; x 2 O; (5.1) n;p;q;r¼1 qxq > : uðxÞ ¼ 0; x 2 qO: It is not difficult to show that this problem has a unique solution. The consequence of this is the fact that the whole sequence u ðxÞ converges to the uniquely defined vectorfunction uðxÞ in L2 ðOÞ when l40. By taking the inverse Laplace transform we show that uðx; tÞ ¼ L1 ðuðx; lÞÞ satisfies the following initial boundary value problem: 8 3 X > q2 uðx; tÞ q > > > rðxÞ  anpqr ðxÞnp ½uðx; tÞer ¼ 0; x 2 O; t40; > 2 > qt qx q > n;p;q;r¼1 < uðx; tÞ ¼ 0; x 2 qO; tX0; >  > > > quðx; tÞ > > > : uðx; 0Þ ¼ 0; qt  ¼ aðxÞ;

(5.2)

x 2 O:

t¼0

Note that the functions u ðx; lÞ and uðx; lÞ were defined for real l only. In order to apply L1 we need to extend these functions into Re l40 and verify that the extended functions have the appropriate rate of decay as jlj ! 1 (see Appendix C for details). Thus we prove the following:

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

659

Theorem 3. Let u be a solution of the problem (1.6)–(1.9). Suppose that the array of the network (particles xðÞ i ) satisfies the non-periodic triangulization condition in the sense of ðÞ Definition 1, and the limits (2.7), (4.1), (4.2) exist. Then the solution u ðxðÞ i ; tÞ ¼ ui ðtÞ ði 2 1; NÞ of the discrete problem (1.6)–(1.8) converges to the solution uðx; tÞ of the homogenized continuum problem (5.2) in P the following sense: the linear interpolation of the solution of the ðÞ i discrete problem u ðx; tÞ ¼ N i¼1 ui ðtÞ  L ðxÞ (see (3.4))converges weakly in L2 ðO  ½0; TÞ (for any T40) to the vector-function uðx; tÞ. 6. Explicit formulae for the elastic moduli of a periodic array of particles We now show the existence of the limit (2.7) for a particular example of a periodic cubic lattice. This periodic example recovers previously known formulae (obtained by the Cauchy–Born rule), but in addition it clearly states the limits of validity of these formulae which were not obtained in previous derivations by the Cauchy–Born rule. As the recent analysis by Friesecke and Theil shows, this rule may fail and thus a rigorous mathematical derivation of the limits of validity is necessary. We consider a cubic lattice in which each vertex is connected by a spring to its nearest neighbors NN (the edges of the unit periodicity cube), to its next nearest neighbors NNN (the diagonals of the faces of the cube) and to its next-to-next nearest neighbors NNNN (the diagonals of the cube). For simplicity, we assume that the elastic constants of all these springs are the same. It is straightforward to carry out our calculation for a more realistic case when there are three different spring constants (macroelasticity parameters) for NN, NNN and NNNN. Our basic periodicity cell is depicted on Fig. 1, where it is shown that each vertex is connected to 33  1 ¼ 26 vertices in the lattice. On this figure a fixed point xi ¼ ð0; 0; 0Þ is shown as a large dark ball and all its neighbors xj are shown as smaller balls. The coordinates of these neighbors are xkj ¼ , 0 or . We prove the following: Theorem 4. For the cubic lattice described above (see also Fig. 1), the elastic moduli in the homogenized equation (5.2) are constants and are given by the following formulae:   pffiffiffi pffiffiffi pffiffiffi pffiffiffi annnn ¼ k 1 þ 2 þ 49 3 ; annpp ¼ anpnp ¼ k 12 2 þ 49 3 ; n; p 2 1; 3, anpqr ¼ 0 in all other cases; k ¼ kij (see (1.5)). Remark. If we introduce the notations l ¼ annnn , m ¼ annpp ¼ anpnp , then the elasticity equations (5.2) have the form 8 3 X > q2 uðx; tÞ q2 ur ðx; tÞ > > rðxÞ  ðl  3mÞ er  mDuðx; tÞ  2mrdiv uðx; tÞ ¼ 0; > > 2 > qt qx2r > r¼1 > > < t40; > uðx; tÞ ¼ 0; x 2 qO; tX0; > >  > >  > quðx; tÞ >  > > : uðx; 0Þ ¼ 0; qt  ¼ aðxÞ; x 2 O: t¼0

x 2 O;

ARTICLE IN PRESS 660

M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

Fig. 1. The basic periodicity cell.

Proof of Theorem 4. We compute the coefficients anpqr ðyÞ using formulae (2.6)–(2.7). Taking into consideration (1.3) and the definition of jnp ðxÞ, we see that for particles which np ðÞ are not located near the faces of the cube K yh , the vector wnp i coincides with j ðxi  yÞ (see ðÞ (2.3)). Thus for an arbitrary particle xi we look for a solution of the Eq. (2.3) of the form np ðÞ wnp i ¼ j ðxi  yÞ þ vi ,

(6.1)

where vi ¼ 0 for internal particles. We call a particle xðÞ i internal if all 26 particles connected y ðÞ y to xðÞ lie inside the cube K . For a particle x 2 K i i h h which is not internal (hereafter a boundary particle), vi a0. Our strategy is to substitute (6.1) into (2.6) and show that contribution of vi vanishes in the limit (2.7). Thus the calculation of anpqr ðyÞ amounts to establishing the following relation:

anpqr ðyÞ ¼ lim lim

1 2

P

ij np ðÞ i;j K y hC  j ðxi h

ðÞ qr ðÞ  xðÞ j Þ; j ðxi  xj Þi

h3

h!0 !0

Substituting (6.1) into (2.6), we obtain

anpqr ðyÞ ¼ lim lim

h!0 !0

1 2

P

ij np ðÞ i;j K y hC  j ðxi h

ðÞ qr ðÞ  xðÞ j Þ; j ðxi  xj Þi

h3

.

(6.2)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

661

8
1 ðÞ þ lim lim hC ij ðvi  vj Þ; jqr ðxðÞ i  xj Þi 2 h!0 !0: i; j y  Kh X ij np ðÞ hC  j ðxi  xðÞ þ j Þ; vi  vj i i; j K y h

9 = 1 X X þ hC ij ðvi  vj Þ; vi  vj i þ h2g 3 jvi j2  3 . ; h y i; j y i

ð6:3Þ

Kh

Kh

We next show that 8
ð6:4Þ

h

(the first term in (6.3) is a finite number). We first estimate the following quantity: A¼ where

P

X

hC ij ðvi  vj Þ; vi  vj i þ h2g 3

0

0 jvi j2 pc  h    @

i Ky h

i; j K y h iKy

X

X0

11=2 jvi j2 A

,

(6.5)

i Ky h

stands for the summation over the boundary particles. Using the inequality

h

ðnpÞ np ðÞ E ðnpÞ ðwnp i ÞpE K y ðj ðxi K yh h

 yÞÞ together with (2.2) and (6.1), we obtain         X   X  ðÞ ðÞ  ij np ðÞ ij np ðÞ    hC  j ðxi  xj Þ; vi  vj i þ  hC  ðvi  vj Þ; j ðxi  xj Þi. Ap   i; j K y   i; j K y h

Then

(6.6)

h

   X   ðÞ ij np ðÞ  Ap2 hC  j ðxi  xj Þ; vi  vj i   i; j K y h     X0  X0 ðÞ ij np ðÞ  ¼ 4 hC  j ðxi  xj Þ; vi ip4c  2  jvi j,   i; j K y i Ky

ð6:7Þ

h

h

where the first inequality follows from the symmetry of the operators C ij , and the equality follows from the fact that C ij ¼ C ji and vi ¼ 0 for the internal particles. Finally, the second inequality in (6.7) follows from (1.3). Next, 11=2 0 11=2 0 11=2 0 X0 X0 X0 X0 h 2 2 jvi jp@ 1A  @ jvi j A pc   @ jvi j A .  y y y y i i i i Kh

Kh

Thus, we establish (6.5).

Kh

Kh

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

662

Second, we prove the estimate X0 g jvi j2 pA  h1þ2 . 2 

(6.8)

i Kz h

We begin with the obvious equality jvi j2 ¼ jvki j2 þ

iki X ðjvki þl j2  jvki þl1 j2 Þ,

(6.9)

l¼1 y ðÞ ðÞ where xðÞ i is a particle located near a cube’s face F h , and the vector xi  xki and the ðÞ ðÞ vectors xki þl  xki þl1 ðl ¼ 1; i  ki Þ are perpendicular to this face. Now we evaluate jvki þl j2  jvki þl1 j2 using Young’s inequality:

jvki þl j2  jvki þl1 j2 ¼ ðjvki þl j  jvki þl1 jÞ  ðjvki þl j þ jvki þl1 jÞ 1 pmðjvki þl j  jvki þl1 jÞ2 þ ðjvki þl j þ jvki þl1 jÞ2 , 4m

ð6:10Þ

where m is any strictly positive number. ðÞ y ðÞ ðÞ Let xðÞ ni  xi be the segment perpendicular to the face F h so that xni and xni lie near the ðÞ face F yh and the opposite face, respectively. The distance between xðÞ i and xni is thus of order h. ðÞ ðÞ Recall that xðÞ ki is an internal particle lying on the segment xni  xi . Note that in order y to sum up over all particles in the cube K h , we first fix a boundary particle xðÞ i and sum ðÞ over all particles lying in the segment xðÞ  x , and then sum over all boundary particles ni i y y xðÞ lying near the face F (that is, summing over all segments perpendicular to F ). i h h Making such summations in (6.9) and applying (6.10), we get 8 9
Kh

Kh

Kh

where the constant c does not depend on  or h. This implies that 8 9 <3 X = X X0 X 1 2 jvi j2 pc jvi j2 þ 2  m  jvi  vj j2 þ 2   jvi j2 . :h i y ; m i y y i i; j y Kh

Kh

(6.12)

Kh

Kh

We now use the standard Korn’s inequality for vector-functions from H 1 ðOÞ (see, for example, Duvaut and Lions (1972)): (Z ) Z 3 X 2 2 2 np ½vðxÞ dx þ jvðxÞj dx . (6.13) kvkH 1 ðOÞ pc O n;p¼1

O

Following the derivation of (3.5) with (6.13) in place of (3.12), we obtain

c 

X0 i; j

2

jvi  vj j þ 

3

X i

! jvi j

2

p

X i; j

hC ij ðvi  vj Þ; vi  vj i þ 3

X i

jvi j2 .

(6.14)

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

663

g

Now we choose m ¼ ðh1þ2 Þ= in (6.10). Since ho1, we obtain 8 < 3 X X0 g X jvi j2 pc 1þg jvi j2 þ h1þ2 hC ij ðvi  vj Þ; vi  vj i 2 : 2 y y h i K i K i; j h h 9 = 3 X X g  þh1þ2 3 jvi j2 þ 1þg jvi j2 ; h 2 i Ky i Ky h h 8 9
which establishes (6.8). Combining (6.8) and (6.5), we get 0 Apc  h  @2

11=2

X0

jvi j2 A

pc  h 

pffiffiffiffi 1þ2g Ah 2 ,

i Ky h

which implies that 0 Apch

3þ2g

;

h@2

X0

11=2 jvi j2 A

g

pch3þ2 .

(6.15)

i Ky h

Now (6.4) follows from (6.7) and (6.16). Clearly, (6.4) and (6.3) imply (6.2). To complete the proof of Theorem 4, we fix all values kij ¼ k and use (1.3) to obtain i X ðÞ ðÞ qr ðÞ hC ij jnp ðxðÞ i  xj Þ; j ðxi  xj Þi ¼ 0, j

pffiffiffi  i X pffiffiffi 4 3 ðÞ ðÞ ij nn ðÞ nn ðÞ 3 hC  j ðxi  xj Þ; j ðxi  xj Þi ¼ 2k 1 þ 2 þ , 9 j pffiffiffi  i X pffiffiffi 8 3 ðÞ ðÞ ij np ðÞ np ðÞ 3 2þ hC  j ðxi  xj Þ; j ðxi  xj Þi ¼ k , 9 j pffiffiffi  i X pffiffiffi 8 3 ðÞ ðÞ ij nn ðÞ pp ðÞ 3 2þ hC  j ðxi  xj Þ; j ðxi  xj Þi ¼ k , 9 j where

Pi j

ðÞ means the summation over all particles xðÞ j connected with the particle xi . Then

anpqr ðyÞ ¼ 0,

ARTICLE IN PRESS 664

M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

annnn ðyÞ ¼ lim lim

pffiffiffi  pffiffiffi 4 3 1X 3 2k 1 þ 2 þ 2 i y 9 Kh

h3 pffiffiffi  3  pffiffiffi 4 3 h  OðÞ 3  k 1 þ 2 þ  9 ¼ lim lim h!0 !0 h3 pffiffiffi   pffiffiffi 4 3 h3 3    k 1 þ 2 þ pffiffiffi 4 pffiffiffi 9 3 3 . ¼ lim lim ¼k 1þ 2þ h!0 !0 9 h3 h!0 !0

ð6:16Þ

Note that in (6.16) we estimated the number of -cells in the mesocube K yh as ððh  OðÞÞ=Þ3 , which is asymptotically equivalent to h3 =3 . Clearly, in this calculation we only need such asymptotic estimates rather than the exact number of such cells. Analogously   1 pffiffiffi 4 pffiffiffi 2þ 3 . annpp ðyÞ ¼ anpnp ðyÞ ¼ k 2 9 So the limits (2.7) exist and the coefficients anpqr ðyÞ are constants as claimed. The theorem is proved. Acknowledgements The authors thank E. Khruslov and L.Truskinovsky for very helpful discussions, suggestions, and for providing useful references. The authors are grateful to both referees for comments and suggestions which led to an improvement of the manuscript. The work of L.B. was supported by NSF Grant DMS-0204637. Appendix A. Proof of Lemma 1 Lemma 1. If the sequence of splines u ðxÞ, defined by (3.4), is bounded uniformly in  in P ðÞ H 1 ðOÞ, then the sequence of piecewise-constant vector-functions u ðxÞ ¼ N i¼1 ui wU i ðxÞ converges to the vector-function uðxÞ defined by (3.13). !0

Proof. It is sufficient to show that ku  u kL2 ðOÞ ! 0. From the (3.4), any coordinate of ð2Þ ð3Þ the spline u ðxÞ ¼ ðuð1Þ  ðxÞ; u ðxÞ; u ðxÞÞ can be written in the form Z x ðkÞ ðÞ uðkÞ ðxÞ ¼ hruðkÞ x 2 U i ; k 2 1; 3,   ðrÞ; dri þ u ðxi Þ; xðÞ i

inside of each Voronoi cell, where the integration is taken along the segment x  xðÞ i . Then   2 N Z Z x X   ðkÞ 2 ðkÞ ðkÞ ku  u kL2 ðOÞ ¼ hru ðrÞ; dri dx.  ðÞ   x i¼1 U i i

Since the segment x 

xðÞ i

can be parameterized as

rðt; xÞ ¼ ð1  tÞxðÞ i þ tx;

t 2 ½0; 1,

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

665

the integral along this segment can be written and bounded as follows: Z x Z 1 Z 1 ðÞ ðkÞ ðkÞ hru ðrÞ; dri ¼ hrr u ðrðt; xÞÞ; x  xi i dtpc   jrr uðkÞ  ðrðt; xÞÞj dt. xðÞ i

0

0

From this we get 2  Z Z 2 Z 1   x 1     ðkÞ 2  pc2 hru ðrÞ; dri pc jrr uðkÞ ðrðt; xÞÞj dt jrr uðkÞ    ðrðt; xÞÞj dt.    xðÞ 0 0 i

In the last inequality we used the Cauchy–Schwartz inequality. Next, we integrate over Voronoi cell U i and sum up over all cells: Z 1 Z 1 Z N Z  N X X 2 2 ðkÞ 2 c2 jrr uðkÞ ðrðt; xÞÞj dt dx ¼ c jr u ðrðt; xÞÞj dx dt r   i¼1

Ui

0

i¼1

# Z 1" Z N X 1 2 ðkÞ 2 ¼ c jrr u ðrÞj dr dt. t3 ð1tÞxðÞ 0 þtU i i¼1 i

0

Ui

In the last equality we used the Fubini theorem to change the order of integration. Since rr uðkÞ  ðrÞ is constant inside of simplices Pa and the integration is taken over the Voronoi cell rescaled by factor t, we obtain the following inequality: Z Z ðkÞ 2 3 2 jrr u ðrÞj dr ¼ t jrr uðkÞ  ðrÞj dr. ð1tÞxðÞ þtU i i

Ui

Hence, N X i¼1

# Z 1" Z N Z X 1 ðkÞ 2 2 2 c jr u ðrÞj dr dt ¼ c jrr uðkÞ r   ðrÞj dr 3 ðÞ t 0 ð1tÞxi þtU i i¼1 U i Z 2 2 ðkÞ 2 ¼ c2 jrr uðkÞ  ðrÞj drpc ku kH 1 ðOÞ 2

0

O 2 !0

pc ! 0. The lemma is proved.

&

Appendix B. Proof of Lemma 2 Lemma 2. For any vector-function ud ðxÞ 2 C 20 ðOÞ (see (4.29)), there exists a sequence of vector-functions ud ðxÞ 2 H 10 ðOÞ of the form (3.4) such that ð1Þ lim kud  ud kL2 ðOÞ ¼ 0; !0

ð2Þ lim kud  u kH 1 ðOÞ ¼ 0; d!0

ð3Þ lim kud  ud kL2 ðOÞ ¼ 0. !0

Proof. Clearly, for any vector-function wðxÞ 2 C 20 ðOÞ there exists a sequence of splines wb ðxÞ ¼

N X i¼1

i wðxðÞ i Þ  L ðxÞ

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

666

such that lim kwb  wkH 1 ðOÞ ¼ 0.

(B.1)

!0

From this, it follows that kwb kH 1 ðOÞ p2kwkH 1 ðOÞ for obðwÞ. Next, using (4.29), we see that for any vector-function wðxÞ 2 H 10 ðOÞ there exists a sequence of vector-functions wb ðxÞ 2 H 10 ðOÞ of the form (3.4) such that ðaÞ lim kwb  wkH 1 ðOÞ ¼ 0; !0

ðbÞ kwb kH 1 ðOÞ p2kwkH 1 ðOÞ

for obðwÞ. So for every vector-function wd ¼ ud  u there exists a sequence of vectorc functions wc d with properties (a) and (b). Now set ud ¼ u þ w d . Then cd  wd kL2 ðOÞ lim kud  ud kL2 ðOÞ ¼ lim ku  u þ w

!0

!0

cd  wd kL2 ðOÞ ¼ 0, p lim ku  ukL2 ðOÞ þ lim kw !0

!0

cd kH 1 ðOÞ p2  lim kwd kH 1 ðOÞ ¼ 2  lim kud  ukH 1 ðOÞ ¼ 0. lim kud  u kH 1 ðOÞ ¼ lim kw

d!0

d!0

d!0

d!0

Thus, properties (1) and (2) hold. Property (3) follows from Lemma 1. Lemma 2 is proved. & Appendix C. Inverse Laplace transform Return to problem (3.1): l2 u þ

1 ru Hðu Þ ¼ a , m 

where N X 3 i X 1 1 X ðÞ ru Hðu Þ ¼ hC ij ðuðÞ i  uj Þ; ek ie3ði1Þþk . ðÞ m m j i¼1 k¼1 i

Here ðÞ u ¼ ðuðÞ 1 ; . . . ; uN Þ;

ðÞ a ¼ ðaðÞ 1 ; . . . ; aN Þ;

1ðÞ 2ðÞ 3ðÞ uðÞ i ¼ ðui ; ui ; ui Þ,

1ðÞ 2ðÞ 3ðÞ aðÞ i ¼ ðai ; ai ; ai Þ

and the vectors e3ði1Þþk ði 2 1; N; k 2 1; 3Þ form an orthonormal basis in R3N . As we have shown, the solution u ðx; lÞ of this problem converges to the solution uðx; lÞ of problem (5.1) in L2 ðOÞ for fixed l40. Define now two operators which map C3N into C3N and H 10 ðOÞ \ H 2 ðOÞ into L2 ðOÞ, respectively: Au ¼

1 ru Hðu Þ; m 

Bu ¼ 

3 X 1 q  fanpqr ðxÞnp ½uðxÞger . rðxÞ n;p;q;r¼1 qxq

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

667

It is easy to see that operator A and the Friedrich’s extension of operator B are both selfadjoint and non-negative if we set Z N X ðÞ ðAu ; v Þ ¼ mðÞ hðAu Þ ; v i 3 ; ðBu; vÞ ¼ rðxÞhBuðxÞ; vðxÞiC3 dx.  i i C i O

i¼1

Rewrite Eqs. (3.1) and (5.1) in the form Cu  ðl2 Þu ¼ f,

(C.1)

where C is a self-adjoint operator. Next, since the eigenvalues of a non-negative operator are non-negative, we see that for any RHS in (C.1) ða 2 C3N ; aðxÞ 2 L2 ðOÞÞ, there exist unique solutions of these problems for Re l40 which are analytic in l (Ahiezer and Glazman, 1963). Denote them as before by u ðx; lÞ and uðx; lÞ, respectively. It is easy to show (since the operators A and B are non-negative) that in the domain Re l4l1 40 ku kL2 ðOÞ p

c ; jlj

kukL2 ðOÞ p

c . jlj

(C.2)

Moreover, the sequence of vector-functions u ðx; lÞ (which are analytic in the domain Re l4l1 ) is compact on this domain and converges in L2 ðOÞ to the analytic vector-function uðx; lÞ for any real lX2l1 . Note that the limiting point 2l1 of the set lX2l1 belongs to Re l4l1 , which allows us to apply the Vitaly theorem (Marcushevich, 1978). It follows from this theorem that the sequence u ðx; lÞ converges uniformly in l on the whole domain Re l4l1 to an analytic vector-function u1 ðx; lÞ. Moreover, u1 ðx; lÞ ¼ uðx; lÞ for real lX2l1 . Then, due to the uniqueness theorem for analytic functions, the vector-functions u1 ðx; lÞ and uðx; lÞ coincide on the whole domain Re l4l1 . Thus, we show that lim ku ðx; lÞ  uðx; lÞkL2 ðOÞ ¼ 0,

!0

(C.3)

uniformly in l for Re l4l1 . We now introduce the notation l ¼ s þ is. Since the solutions u ðx; lÞ and uðx; lÞ of the problems (3.1) and (5.1) are the Laplace transforms of the solutions u ðx; tÞ and uðx; tÞ of the problems (1.6)–(1.8) and (5.2), respectively, we have Z 2l1 þi1 Z 2l1 þi1 1 1 lt u ðx; tÞ ¼ e u ðx; lÞ dl; uðx; tÞ ¼ elt uðx; lÞ dl. (C.4) 2pi 2l1 i1 2pi 2l1 i1 Let cðxÞ be an arbitrary vector-function of the class L2 ðOÞ, and let jðtÞ be an arbitrary function of the class C 10 ð0; TÞ. We now show that Z TZ hu ðx; tÞ  uðx; tÞ; cðxÞjðtÞi dx dt ¼ 0; 8T40. (C.5) lim !0

0

O

Indeed, (C.2), (C.4) and the condition jðtÞ 2 C 10 ð0; TÞ (from this condition it follows R T using lt 1 that 0 e jðtÞ dt ¼ Oðjlj Þ), we have Z T Z     hu ðx; tÞ  uðx; tÞ; cðxÞjðtÞi dx dt lim !0

0

O

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

668

Z  ¼ lim !0

0

T

Z

O

1 2pi

Z

2l1 þi1 2l1 i1

  elt ðu ðx; lÞ  uðx; lÞÞ dl; cðxÞjðtÞ dx dt

Z  Z Z T  1  2l1 þi1 lt ¼ lim  dl hðu ðx; lÞ  uðx; lÞÞ; cðxÞi dx  e jðtÞ dt !0 2p 2l i1 0 O 1 Z þ1 ds pc lim ku ðx; 2l1 þ isÞ  uðx; 2l1 þ isÞkL2 ðOÞ !0 1 s Z p ds ¼ c lim ku ðx; 2l1 þ isÞ  uðx; 2l1 þ isÞkL2 ðOÞ !0 s p Z p ds þ ku ðx; 2l1 þ isÞ  uðx; 2l1 þ isÞkL2 ðOÞ s 1  Z þ1 ds þ ku ðx; 2l1 þ isÞ  uðx; 2l1 þ isÞkL2 ðOÞ . s p

ðC:6Þ

Due to (C.2), the last two integrals in s in (C.6) are uniformly small in  for sufficiently large p. The first integral in s vanishes as  ! 0 due to (C.3). Thus (C.5) follows. Next, P since the set of vector-functions f i ci ðxÞji ðtÞ; ci ðxÞ 2 L2 ðOÞ; ji ðtÞ 2 C 10 ð0; TÞg is dense in L2 ðO  ½0; TÞ, Theorem 3 follows. Due to standard properties of the Laplace transform and (5.1), we conclude that the vector-function uðx; tÞ is a solution of the following initial boundary value problem: 8 3 2 X > q > > rðxÞ q uðx; tÞ  > anpqr ðxÞnp ½uðx; tÞer ¼ 0; x 2 O; t40; > 2 > qt qx q > n;p;q;r¼1 < (C.7) uðx; tÞ ¼ 0; x 2 qO; tX0; >  > >  > quðx; tÞ > > > : uðx; 0Þ ¼ 0; qt  ¼ aðxÞ; x 2 O: t¼0

References Ahiezer, N.I., Glazman, I.M., 1963. Theory of Linear Operators in Hilbert Space. Frederic Ungar, New York. Blanc, X., Le Bris, C., Lions, P.-L., 2002. From molecular models to continuum mechanics. Arch. Rational Mech. Anal. 164, 341–381. Born, M., Huang, K., 1954. Dynamical Theory of Crystal Lattices. Oxford University Press, Oxford. Braides, A., Dal-Maso, G., Garroni, A., 1999. Variational formulation of softening phenomena in fracture mechanics: the one-dimensional case. Arch. Rational Mech. Anal. 146, 23–58. Duvaut, G., Lions, J.-L., 1972. Inequalities in Mechanics and Physics. Dunod, Paris (French). E, W., Ming, P., 2004. Analysis of multiscale problems. J. Comput. Math. 22, 210–219. Friesecke, G., James, R.D., 2000. A scheme for the passage from atomistic to continuum theory for thin film, nanotubes and nanorods. J. Mech. Phys. Solids 48, 1519–1540. Friesecke, G., Theil, F., 2002. Validity and failure of the Cauchy–Born rule in 2D mass-spring lattice. J. Nonlinear Sci. 12, 445–478. Krasniansky, M., 1997. Vector Homogenized Model for Microinhomogeneous Finite-Difference Equation, GAKUTO International Series Mathematical Sciences and Applications, vol. 9. Homogenization and Applications to Material Sciences, pp. 257–270. Marcushevich, A.I., 1978. The Short Course of Theory of Analytic Functions. Nauka, Moscow (Russian). Mizohata, S., 1973. The Theory of Partial Differential Equations. Cambridge University Press, XI, London.

ARTICLE IN PRESS M. Berezhnyy, L. Berlyand / J. Mech. Phys. Solids 54 (2006) 635–669

669

Truskinovsky, L., 1996. Fracture as a phase transitions. In: Batra, R.C., Beatty, M.F. (Eds.), Contemporary Research in the Mechanics and Mathematics of Materials. CIMNE, Barcelona, pp. 322–332. Vogelius, M., 1991. A homogenization result for planar, polygonal networks. Math. Modelling Numer. Anal. 25 (4), 483–514.

Further Reading Braides, A., Gelli, M.S., 2002. Continuum limits of discrete systems without convexity hypotheses. Math. Mech. Solids 7, 41–46. Khruslov, E.Ya., 1979. Asymptotic behavior of solutions of the second boundary value problem under fragmentation of the boundary of the domain. Math. USSR-Sb. 35. Khruslov, E.Ya., 1981. On convergence of solutions of the second boundary value problem in weakly connected domains, Theory of Operators and Functional Spaces and its Application, Naukova dumka, Kyiv (Russian). Kosevich, A.M., 1988. The theory of crystal lattice (the physical mechanics of crystal), Publishing House at Kharkov State University of Publishing Union Vyscha Shkola (Russian). Kozlov, S.M., 1987. Averaging of difference schemes, Math. USSR, Sb. Krasniansky, M.B., 1995. Homogenization of random walks on lattices with weak connections. Math. Phys. Anal. Geometry, N 1, 51–67. Marchenko, V.A., Khruslov, E.Ya., 1974. Boundary Value Problems in Domains with Fine Grained Boundary. Kyiv, Naukova dumka (Russian). Martynenko, A.I., 1990. Operating Calculation. Vyscha Shcola, Kyiv (Russian). Oleinic, O.A., Shamaev, A.S., Iosif’yan, G.A., 1982. Mathematical problems in elasticity and homogenization. Studies in Mathematics and its Applications, vol. 26, North-Holland Publishing Co., Amsterdam. Timoshenko, S.P., 1957. History of Science about Resistance of Materials. Moscow, State Publishing House of Technical and Theoretical Literature (Russian). Timoshenko, S.P., Goodier, J.N., 1970. Theory of Elasticity, third ed. Engineering Societies Monographs, International Student Edition. New York, Tokyo. Zhikov, V.V., Kozlov, S.M., Oleinik, O.A., 1994. Homogenization of Differential Operators. Springer, New York.