VOLUME 36, NUMBER 4

PHYSICAL REVIEW B

1

AUGUST 198'7

Cauchy relations for central-force random networks

E. J. Garboczi World Industries, Inc. , Research and Development, 2500 Columbia Avenue, P. O. Box 3511, Lancaster, Pennsylvania 17604 (Received 24 November 1986)

Armstrong

For many years it has been known that the elastic moduli of networks whose only interatomic forces are central pair potentials and whose atoms occupy centers of inversion at equilibrium obey the Cauchy relations. Recently it has been shown analytically that such networks with one bond missing (thus eliminating all centers of inversion at atomic sites) still obey the Cauchy relations exactly, while similar networks with one site missing do not. The usual hypothesis that must be satisfied in order for the Cauchy relations to hold must be generalized. This paper presents a possible hypothesis and random netexplores its implications via computer simulation of three different two-dimensional works.

I.

INTRODUCTION

The elastic stiffness tensor of a solid is a fourth-rank tensor, written as C;jI,I. In three dimensions, there are 3 =81 elements. Hooke's law is written as

~jiCij kl ~kl with summation over k, l implied. o.;~ is the stress tensor and cI,I is the strain tensor. Both the stress and the strain tensors can be shown to be symmetric in their respective k~1. This indices, imparting to C~~I the symmetry reduces the number of possible independent elements of density C;~1,.~ to 36. The existence of the strain-energy adds the symmetry (ij)~(kl), further reducing the possible number of independent elements of the elastic stiffness " are defined by the tensor to 21.' The Cauchy relations additional symmetry making C~ I,I completely symmetric in all four indices and thus having at most 15 independent elements. This last symmetry implies the existence of six Cauchy relations:

i~j,

j~k,

—C1212, C1123 = 12 3 , C1122

C1133 1233

C1313

=

~

C2233

1323 ,

1322

C2323

=

1223

or, written in the Voigt notation, C13 = Css C14 — C6s

~

C63

= Cs4,

C23

= C44

Cs2

C64

The usual hypothesis that is sufficient for the Cauchy relations to be obeyed by the elastic moduli of a network is that each atom must be at a center of inversion and the only interatomic forces must be central pair potentials. In the elastic limit, each pair potential can be of specified as a force constant. It is the arrangement atoms plus force constants that must be invariant under inversion as stated in the above hypothesis. Recent results, however, indicate that this hypothesis can be somewhat relaxed. It was shown that the elastic moduli of a Bravais lattice with one bond's central-force 36

constant different from the rest still obeyed the Cauchy relations exactly. Numerical results implied that the Cauchy relations were still being obeyed (within numerical error) as a finite fraction of bonds were randomly changed. Of course, all atomic centers of inversion are lost when even one bond is different from the rest. One might be tempted to generalize the hypothesis for the Cauchy relations to allow for random arrangements of force constants on a centrosymmetric lattice, thus preserving some sort of average or macroscopic inversion symmetry. However, the situation is more complicated than that. A recent result' for site percolation on central-force networks showed that the elastic moduli of a central-force Bravais lattice rigorously do not obey the Cauchy relations when the bonds belonging to a single site are different from the rest. Thus in a generalized hypothesis for the Cauchy relations, one must not allow simply for random arrangements of force constants; one must also specify what kinds of random distributions are allowed. This paper attempts to construct a generalized hypothesis for the Cauchy relations that is consistent with all previous results. This hypothesis is stated in Sec. II. Sections III —V present the results of computer simulations of three different models that give further support to the sufficiency of this hypothesis. Section VI concludes this paper with a discussion of results and suggestions for further work.

II.

GENERALIZED HYPOTHESIS FOR THE CAUCHY RELATIONS

The elastic moduli of networks whose bonds were removed (or had force constants changed) according to bond-percolation statistics were found to obey the Cauchy relations within numerical error for all concentrations of bond defects, and exactly in the limit of a single-bond dedefects). This result was fect (or a few nonoverlapping In found to be true in two and three dimensions. bond percolation, the probability that a given bond is present is p, the fraction of bonds present (0 &p & 1). This

1987

The American Physical Society

E. J. GARBOCZI

2116

probability does not depend on the presence or absence of any other bond. By way of contrast, networks whose bonds were resites had elastic moved by the deletion of individual moduli that did not obey the Cauchy relations at any value of p„ the fraction of sites present. In site percolation, each site has a probability p, of being present, which does not depend on the presence or absence of any other site. However, the probability of a given bond being present is p, and is dependent on the status of nearby bonds. For example, if the bond connecting sites 3 and B is missing, then all the bonds connected to site 3 must be missing or all the bonds connected to site must be missing or both, because site 3 or site B or both must have been removed. These differences between site and bond-percolation bond probabilities was the motivation for the following which is conjectured to be generalized hypothesis, sufhcient for the elastic moduli of a network to obey the Cauchy relations. The hypothesis is in two parts: (1) The sites of the network must be arranged so that, neglecting the values of the force constants, there is a center of inversion symmetry at every site, and (2) only central pair potentials are allowed as interatomic forces; but the value of each central-force constant in the elastic limit may be randomly selected, as long as its value is selected independently of the value of any other force constant. If requirements (1) and (2) are met, then it is proposed that the elastic moduli of such a network will obey the Cauchy relations. It is clear that this conjecture includes the former hypothesis as a special case. It should be noted at this point that in two dimensions, choosing z to be the out-of-plane coordinate, the Cauchy relations (3) reduce to the single relation C~2 = C66 (xx =1, yy =2, xy =6). All the models studied in the following sections are two dimensional.

8

III.

RANDOM GEOMETRY MODEL

different, deviating from unity by order s. The connectivity of the unperturbed lattice is preserved. The lattice potential (4) becomes

U= 'a g [(u; —u, ).5,'J ]', —,

&ij

= a g I(u; —u, ).5;, ]', '

j

j

K = —,'(C„+C„).

The simulations were performed by redefining the lattice vectors that define the supercell to reflect a macroscopic strain E. All sites were then allowed to relax to their new equilibrium positions while periodic boundary conditions were maintained. After every 50 moves per site, the effective modulus was computed via the total energy, —,'C;JE . When the last change after 50 moves was less than 10 of the previous value of the modulus, the relaxation procedure was judged to be complete. Figure 1 shows Ci~, C66, and Ci2 plotted against s . the real-space matrix dynamical By constructing D =Dp+ V from the potential (5), where V is the "defect" matrix and Dp is the dynamical matrix for potential (4), and examining V in first-order perturbation theory, one can check that the lowest-order correction to the elastic moduli goes as not s. The straight lines in Fig. 1 are the extrapolation of the computed initial slope for the elastic moduli versus s . This slope was computed from the unperturbed elastic moduli and the values of the moduli at s = 10 . The agreement with the data points is quite good out to about s =0. 16 (s =0.4). Figure 2

s,

&.

5

1. 0

(4)

—,

&

where 5,'~ is a unit vector from i to that is different for each i, nearest-neighbor pair and u, , uz are measured from the new equilibrium state. Computer simulations were performed on this model to compute C„, C«, and C» as a function of s. A supercell with 40X46=1840 sites was used under periodic boundTen independent were configurations ary conditions. averaged to determine the elastic moduli at each value of s, with C„being averaged over the x and y directions as well. C» was calculated from direct computations of Cll bulk and the modulus the relation K, using

The random geometry model is studied to investigate whether part (1) of the generalized Cauchy hypothesis, the stipulation that the geometric arrangement of the sites of an elastic network must be centrosymmetric, is needed. The model consists of a two-dimensional triangular net with nearest-neighbor Hooke-spring forces that all have force constants a. Before any changes are made, the harmonic lattice potential is Up

36

0. 5

0

(ij )

where u, , uj are the displacements from equilibrium of the ith and jth sites, the angular brackets under the summation sign restrict the summation to nearest-neighbor pairs, and 6,~ is a unit vector from i to It should be noted that there are only six different values of 6,z in the crystalline triangular net. The nearest-neighbor spacing is 1. This network is then perturbed by displacing each site by the vector s (cosO, sinO), where O is a random variable evenly distributed between 0 and 2m, and s & 1. This new atomic arrangement is theo defined as the new equilibrium position. Each "nearest-neighbor" bond length is now

66

j.

O. 5

0 00 ~

0. 05

12

0. 10

I

I

0. 15

0. 20

0. 25

S2

FIG. 1. Showing Cll, C66, and C» vs s' for the randomgeometry model. The solid lines are extrapolations of the computed initial slope.

CAUCHY RELATIONS FOR CENTRAL-FORCE RANDOM NETWORKS

36

2117

1. 5

TABLE I. Elastic moduli results from simulations of a 100 && 116= 11 600 site central-force triangular network. The nearest-neighbor force constants were evenly distributed as indicated. Each modulus is an average over 20 independent

1. 0

configurations, with CI2 averaged over the x and y directions as well. I is defined as I =(o.2 — ai)/o. i, with o. i &a &a2.

1&a&1.2

0. 5

1

0. 0

=0.2

Cee Cl~

&o. &2.0 I =1.0

C12

I

0. 00

Cee

0. 25

Cia 1

Showing C=(Cee —CI2)/Cee vs s for the randomgeometry model. The solid line is an extrapolation of the computed initial slope.

EMT

Modulus

Distribution

&a &200

Cee

0.475 438

0.475 433+0.000 176 0.475 435+0.000 673

0.633 12

0.633 11+0.000 88 0.633 12+0.000 93

3.519

3.522+0. 019 3.531+0.078 30.91+0.22

30.87

=199

I

Simulation

31.09+0. 31

FIG. 2.

=1.

shows C=(C66 — C~q)/C66 plotted versus s . It is clear that C ~2& C«even for the smallest values of s . The straight line is the extrapolated initial slope. The agreement with the data points is worse than in Fig. 1, because the division of two uncertain numbers always introduces more error.

IV. FLAT DISTRIBUTION OF FORCE-CONSTANTS MODEL The model studied in this section is a triangular network with the harmonic lattice potential given in Eq. (4), but with nearest-neighbor central-force constants e distributed such that a& &n &aq, with o; equally likely to take on any value in the interval. If the generalized hypothesis is correct, then this network should obey the Cauchy relations. Computational results were averaged over 20 configurations of a 100&& 116=11600 site supercell with periodic boundary conditions. Results for CI2 were also over the x and y directions for each averaged The relaxation procedure of Sec. III was configuration. used. Table I contains the results for the model for ai —1, and o.'2 —1.2, 2.0, 20, and 200. Ci~, C«, and the bulk modulus K were computed independently [approximately 850 h of central-processing-unit (CPU) time on a VAX 11/785 computer, double precision], and C~2 was calculated via the relation K= —,'(C~~+C~q). Only C~2 and C66 results are presented in the table. The standard deviation in each modulus was calculated over the 20 The elfective-medium-theory configurations. (EMT) results were obtained using the single-bond theory derived in Ref. 7, with the probability distribution P(a) = 1/ 'a2 — cx)), a) &a &cz2. The parameter I is defined as I =(a2 —a~)/a~, and so is a measure of the disorder present in the force constant distribution. When I =0, there are only a] force constants present, and so C&z must equal C66. As can be seen in Table I, for I =0. 2, C~2 and C66

0 the differ by less than 1 part in 200000, while for I two moduli differ by about 1 part in 60000. In both cases, C~q —C66 is much less than o. tq or o 66, the standard deviations in the two moduli. In the case of the larger I values, I =19 and 199, the disagreement between C~2 and C66 is rather larger (1 part in 350 and 150, respectively), but C~2 —C66 is still less than cr~q or cr66 It is important to note that for a2 »a~ (I &&1), the fiucbetween simulations of tuations independent configurations were larger both in an absolute and relative sense than for smaller I . For example, o. )2/C)2 -—0.001 for I =0.2, while for I =199, o. ]2/C~q-— 0. 01, an order of magnitude increase. However, for all four values of I considered, one can state that C ~2 and C«are equal within "experimental" error. ~

V. CORRELATED BOND-PERCOLATION MODEL

The correlated bond-percolation model was originally described in Ref. 11. A slightly modified version was used in this work. Qne assigns a random number x, to ' ' each site of a lattice, where — —, &x, & —, and x, is evenly distributed. Then for each bond connecting sites i and j, one variable defines the random x„=—,' ( x, + x, + x, —x, ), where 0&x„& 1. One then chooses x, 0&x & 1, and removes all bonds with x,-, ~x. It is not difficult to see that the actual fraction of bonds remaining after this operation is carried out is p (x), the distribution function of x treated as a random variable. Using elementary techniques one can show that -

~

~

p='

~

~

1

—(1 —x)

for

—,

3x

for

0&x &

'

&x

&

1, '

—,

I,

Figure 3 shows p versus

dp/dx at x

= '.

x. Note the discontinuity

in

—,

networks Simulations were performed on triangular with values of x ranging between 0. 1 and 0.9. However,

E. J. GARBOCZI

2118

0. 12

1. 00

0. 75

o. oa

N

Q

0. 50

0. 06

E

(3 C3

0. 25

0. 00

0. 03

CL

0. 00

0. 75

0. 50

0. 25

0 00

1. 00

0. 0

0. 2

0. 4

0. 6

FIG. 3. Showing p vs x for a correlated bond-percolation process, where x is the parameter defined in Sec. V and p =p(x) is the fraction of bonds present.

FIG. 5. Showing the parameters C=(Cl~ —C«)/Cf, t, (circle) C l2 —2C«) /C«(square) vs p. and = ( C l — and C were calculated from independent computations of Cll, K, and C«.

instead of removing bonds with x,,-)x, such a bond's force constant was changed from 1 to 0.2. This procedure avoided the numerical difriculties of the rigidity phase transition at p==', (x =0.47), but still allowed a fairly large range of elastic moduli to be explored. Fifteen configurations of a 60& 70=4200 site supercell were aver750 h of CPU aged for each value of x. Approximately time on a VAX 11/785, in double precision, were required to complete the correlated bond-percolation simulations. The relaxation procedure of Sec. III was used. Figure 4 shows K and C«plotted as a function of the bond fraction p. In Figs. 4 —6, p is the fraction of bonds with force constant o, =1.0. The solid lines are the generalized EMT derived in Ref. 9, with p (x) used instead of the bond-percolation probability. It is seen in Fig. 4 that

this EMT is correct as p~1 and p~0, as both of these cases are in the few-defect limit. When x = 1 —b (or x =b)), for any network one can choose 6 sufficiently small so that only one bond is diff'erent from the rest. The EMT is exactly correct in the one-bond-defect limit. For 0.2 &p &0. 8, the EMT is clearly incorrect, though not by any great amount. Similar behavior was also seen in the conductivity case. Figure 5 shows C=(C~q — C66)IC66 and I=(C~~ —Ci2 —2C66)/C66 versus p. Cauchy's relation is clearly violated for intermediate values of p and is obeyed as The parameter is a measure of the anand p p isotropy of the network. The infinite triangular net with all bonds the same is isotropic and should remain so as bonds are randomly changed. The small nonzero values

"

I

I

&

~1.

~0

I

2. 2

1. 0

2. 1

0. 6 D

0. 4

2. O

0. 2 0. 0

FIG. 4.

0. 0

l

l

0. 2

0. 4

Showing the bulk modulus

0. 6

! 0. 8

l

l

o. 4

0. 6

0. 8

K/C«vs

p for the data plotted in

I

0. 0

0. 2

K and shear modulus C«

p, where p is the fraction of bonds with force constant o.'=1. A fraction 1 — p of the bonds have n =0.2. vs

1. a

1. 0

FIG. 6. Fig. 4.

Showing the ratio

CAUCHY RELATIONS FOR CENTRAL-FORCE RANDOM NET%'ORKS

36

of I observed are due to the effect of a network with finite size. I is much smaller than C, so that the violation of Cauchy's relation is not a finite-size effect. Figure 6 is a plot of the ratio K/C66 versus p. For an isotropic two-dimensional network that satisfies the Cauare nonzero, we chy relations, K/C66 ——2. When C and can write

I

K /C66

—2+ C+

I,

' —,

so that the deviation of K/C66 from 2 in Fig. 6 is due almost entirely to the violation of Cauchy's relation and not finite-size anisotropy.

VI. DISCUSSION The generalized hypothesis stated in Sec. II is essentialto classify what kinds of random centralforce systems will obey the Cauchy relations. It should be emphasized that by "random" is meant random central-force constants superimposed on a centrosymmetric lattice. The three models studied above have given support to the suSciency of this hypothesis. The first model studied, the random-geometry model, showed that the condition of having atoms at the sites of a centrosymmetric (neglecting the force constants) lattice must be retained in any generalized Cauchy hypothesis. The model described in Sec. IV had nearest-neighbor force constants that were randomly and independently distributed on a centrosymmetric lattice. One should take note that every nearest-neighbor bond had a different force constant. For I =0.2 and 1.0, the Cauchy relation was found to be obeyed to 1 part in 200000 and 1 part in 60000, respectively, thus confirming the prediction of the generalized Cauchy hypothesis. Since I is a measure of the force-constant disorder, it should be emphasized that these two values of I represent a substantial perturbation of the network. It is also interesting to note that in these two cases the EMT and simulation results agreed within the same small amount. One might speculate that this EMT could be exact for this model. It is known that this EMT does not predict the correct critical threshold and exponent for models studied previously that undergo a riIt is possible that the close agreement gidity transition. with EMT in the I =0.2 and 1.0 cases was due to the fairly small range over which the force constants were distributed. The I =19 and 199 simulations were carried out to check whether for o. 2»o. the Cauchy relation would still be obeyed but with increasing disagreement between EMT and simulation. It was found that EMT and simulation results still only differed by less than a standard deviation, with this difference being of the same orly an attempt

C„—

der as C« . Limitations on computer time prohibited using more configurations on larger networks to reduce the standard deviations at the two higher I values. The computer time necessary to completely relax a network increases rapidly with I . Of the 850 CPU hours required for the simulation results given in Sec. IV, 700 h were spent just on the I =19 and 199 computations. The question of the exactness of this EMT for this model must remain open for the present. It should be noted, however, that if the Cauchy relations are obeyed by this or any random system, then some form of renormalized forceconstant description must be exact. The correlation bond-percolation model in Sec. V was chosen because it is intermediate between ordinary bond and site percolation, in the following sense: As x decreases from 1 (or increases from zero), only one bond is altered initially, so that the single-bond EMT is exact for p very close to 1 (or 0), but correlations quickly build up as more bonds are altered that cause the distribution of unaltered bonds to differ significantly from simple bond percolation. The conjecture would predict that the Cauchy relation would be obeyed for p close to 1 (or zero) with the largest disagreement near p=0. 5. This was exactly the behavior observed. In Fig. 5, the parameter ( C~2 C66 )/C66 increased as p decreased from 1, and had its maximum value at about p =0. 5. C then decreased as p decreased further, as the +=0.2 bonds became dominant and the o:=1.0 bonds became the "de~

~

C:

"

fects. The next step in this work would seem to be an attempt to formulate an analytical proof that the generalized Cauchy hypothesis given here does indeed lead to the Cauchy relations. This is not an easy task. Lacking a general proof, it would be useful to take the flat distribution of force constants model described in Sec. IV and attempt to show analytically that the Cauchy relations hold exactly for this specific case. Such a proof might be simply generalized. Of course, even if one could formulate a proof using the full hypothesis presented in this paper, this would still not necessarily eliminate the possibility of an even more general hypothesis being su%cient for the CauAddition careful numerical work on chy relations. different random systems might find a clear counterexample to support such further generality. A. CENO%I. EDGMENTS

~

In the Voigt notation, the six independent combinations of x, y, and z are written as xx=1, yy=2, zz=3, yz=4, zx =5, and xy— = 6. In this notation C;,kl is written as C „, m, n =1,6. The symmetry (ij)~(kl) is equivalent to saying that C „ is a symmetric second-rank tensor. A. Cauchy, Exercises Math. 3, 213 (1828). B. Saint-Venant, Memoires des Savants e'trangers 14, (1855).

2119

I would like to thank M. F. Thorpe for many helpful conversations about the Cauchy relations, and E. W. Behrens for a helpful conversation about elasticity theory.

~A. E. H. Love, A Treatise on the Mathematical Theory of Elasti city (Cambridge University Press, Cambridge, 1934). sM. Born and K. Huang, Dynamical Theory of Crystal Lattices (Clarendon, Oxford, 1966), Chap. 3. 6J. H. Weiner, Statistical Mechanics of Elasticity (Wiley, New York, 1983), Chap. 9; I. Stakgold, Q. Appl. Math. 8, 169

(1950-1951).

2120 7S. Feng, M. F. Thorpe, and E. J. Garboczi, 276 (1985). 8E. J. Garboczi and M. F. Thorpe, Phys. Rev. E. J. Garboczi and M. F. Thorpe, Phys. Rev. ' M. F. Thorpe and E. J. Garboczi, Phys.

E. J. GARBOCZI Phys. Rev. B

3I,

B 31, 7276 (1985). B 33, 3289 (1986). Rev. B 35, 8579

36

(1987). ''Scott Kirkpatrick, Rev. Mod. Phys. 45, 574 (1973).

' This transition

smaller systems.

was seen in earlier unpublished

work done on

PHYSICAL REVIEW B

1

AUGUST 198'7

Cauchy relations for central-force random networks

E. J. Garboczi World Industries, Inc. , Research and Development, 2500 Columbia Avenue, P. O. Box 3511, Lancaster, Pennsylvania 17604 (Received 24 November 1986)

Armstrong

For many years it has been known that the elastic moduli of networks whose only interatomic forces are central pair potentials and whose atoms occupy centers of inversion at equilibrium obey the Cauchy relations. Recently it has been shown analytically that such networks with one bond missing (thus eliminating all centers of inversion at atomic sites) still obey the Cauchy relations exactly, while similar networks with one site missing do not. The usual hypothesis that must be satisfied in order for the Cauchy relations to hold must be generalized. This paper presents a possible hypothesis and random netexplores its implications via computer simulation of three different two-dimensional works.

I.

INTRODUCTION

The elastic stiffness tensor of a solid is a fourth-rank tensor, written as C;jI,I. In three dimensions, there are 3 =81 elements. Hooke's law is written as

~jiCij kl ~kl with summation over k, l implied. o.;~ is the stress tensor and cI,I is the strain tensor. Both the stress and the strain tensors can be shown to be symmetric in their respective k~1. This indices, imparting to C~~I the symmetry reduces the number of possible independent elements of density C;~1,.~ to 36. The existence of the strain-energy adds the symmetry (ij)~(kl), further reducing the possible number of independent elements of the elastic stiffness " are defined by the tensor to 21.' The Cauchy relations additional symmetry making C~ I,I completely symmetric in all four indices and thus having at most 15 independent elements. This last symmetry implies the existence of six Cauchy relations:

i~j,

j~k,

—C1212, C1123 = 12 3 , C1122

C1133 1233

C1313

=

~

C2233

1323 ,

1322

C2323

=

1223

or, written in the Voigt notation, C13 = Css C14 — C6s

~

C63

= Cs4,

C23

= C44

Cs2

C64

The usual hypothesis that is sufficient for the Cauchy relations to be obeyed by the elastic moduli of a network is that each atom must be at a center of inversion and the only interatomic forces must be central pair potentials. In the elastic limit, each pair potential can be of specified as a force constant. It is the arrangement atoms plus force constants that must be invariant under inversion as stated in the above hypothesis. Recent results, however, indicate that this hypothesis can be somewhat relaxed. It was shown that the elastic moduli of a Bravais lattice with one bond's central-force 36

constant different from the rest still obeyed the Cauchy relations exactly. Numerical results implied that the Cauchy relations were still being obeyed (within numerical error) as a finite fraction of bonds were randomly changed. Of course, all atomic centers of inversion are lost when even one bond is different from the rest. One might be tempted to generalize the hypothesis for the Cauchy relations to allow for random arrangements of force constants on a centrosymmetric lattice, thus preserving some sort of average or macroscopic inversion symmetry. However, the situation is more complicated than that. A recent result' for site percolation on central-force networks showed that the elastic moduli of a central-force Bravais lattice rigorously do not obey the Cauchy relations when the bonds belonging to a single site are different from the rest. Thus in a generalized hypothesis for the Cauchy relations, one must not allow simply for random arrangements of force constants; one must also specify what kinds of random distributions are allowed. This paper attempts to construct a generalized hypothesis for the Cauchy relations that is consistent with all previous results. This hypothesis is stated in Sec. II. Sections III —V present the results of computer simulations of three different models that give further support to the sufficiency of this hypothesis. Section VI concludes this paper with a discussion of results and suggestions for further work.

II.

GENERALIZED HYPOTHESIS FOR THE CAUCHY RELATIONS

The elastic moduli of networks whose bonds were removed (or had force constants changed) according to bond-percolation statistics were found to obey the Cauchy relations within numerical error for all concentrations of bond defects, and exactly in the limit of a single-bond dedefects). This result was fect (or a few nonoverlapping In found to be true in two and three dimensions. bond percolation, the probability that a given bond is present is p, the fraction of bonds present (0 &p & 1). This

1987

The American Physical Society

E. J. GARBOCZI

2116

probability does not depend on the presence or absence of any other bond. By way of contrast, networks whose bonds were resites had elastic moved by the deletion of individual moduli that did not obey the Cauchy relations at any value of p„ the fraction of sites present. In site percolation, each site has a probability p, of being present, which does not depend on the presence or absence of any other site. However, the probability of a given bond being present is p, and is dependent on the status of nearby bonds. For example, if the bond connecting sites 3 and B is missing, then all the bonds connected to site 3 must be missing or all the bonds connected to site must be missing or both, because site 3 or site B or both must have been removed. These differences between site and bond-percolation bond probabilities was the motivation for the following which is conjectured to be generalized hypothesis, sufhcient for the elastic moduli of a network to obey the Cauchy relations. The hypothesis is in two parts: (1) The sites of the network must be arranged so that, neglecting the values of the force constants, there is a center of inversion symmetry at every site, and (2) only central pair potentials are allowed as interatomic forces; but the value of each central-force constant in the elastic limit may be randomly selected, as long as its value is selected independently of the value of any other force constant. If requirements (1) and (2) are met, then it is proposed that the elastic moduli of such a network will obey the Cauchy relations. It is clear that this conjecture includes the former hypothesis as a special case. It should be noted at this point that in two dimensions, choosing z to be the out-of-plane coordinate, the Cauchy relations (3) reduce to the single relation C~2 = C66 (xx =1, yy =2, xy =6). All the models studied in the following sections are two dimensional.

8

III.

RANDOM GEOMETRY MODEL

different, deviating from unity by order s. The connectivity of the unperturbed lattice is preserved. The lattice potential (4) becomes

U= 'a g [(u; —u, ).5,'J ]', —,

&ij

= a g I(u; —u, ).5;, ]', '

j

j

K = —,'(C„+C„).

The simulations were performed by redefining the lattice vectors that define the supercell to reflect a macroscopic strain E. All sites were then allowed to relax to their new equilibrium positions while periodic boundary conditions were maintained. After every 50 moves per site, the effective modulus was computed via the total energy, —,'C;JE . When the last change after 50 moves was less than 10 of the previous value of the modulus, the relaxation procedure was judged to be complete. Figure 1 shows Ci~, C66, and Ci2 plotted against s . the real-space matrix dynamical By constructing D =Dp+ V from the potential (5), where V is the "defect" matrix and Dp is the dynamical matrix for potential (4), and examining V in first-order perturbation theory, one can check that the lowest-order correction to the elastic moduli goes as not s. The straight lines in Fig. 1 are the extrapolation of the computed initial slope for the elastic moduli versus s . This slope was computed from the unperturbed elastic moduli and the values of the moduli at s = 10 . The agreement with the data points is quite good out to about s =0. 16 (s =0.4). Figure 2

s,

&.

5

1. 0

(4)

—,

&

where 5,'~ is a unit vector from i to that is different for each i, nearest-neighbor pair and u, , uz are measured from the new equilibrium state. Computer simulations were performed on this model to compute C„, C«, and C» as a function of s. A supercell with 40X46=1840 sites was used under periodic boundTen independent were configurations ary conditions. averaged to determine the elastic moduli at each value of s, with C„being averaged over the x and y directions as well. C» was calculated from direct computations of Cll bulk and the modulus the relation K, using

The random geometry model is studied to investigate whether part (1) of the generalized Cauchy hypothesis, the stipulation that the geometric arrangement of the sites of an elastic network must be centrosymmetric, is needed. The model consists of a two-dimensional triangular net with nearest-neighbor Hooke-spring forces that all have force constants a. Before any changes are made, the harmonic lattice potential is Up

36

0. 5

0

(ij )

where u, , uj are the displacements from equilibrium of the ith and jth sites, the angular brackets under the summation sign restrict the summation to nearest-neighbor pairs, and 6,~ is a unit vector from i to It should be noted that there are only six different values of 6,z in the crystalline triangular net. The nearest-neighbor spacing is 1. This network is then perturbed by displacing each site by the vector s (cosO, sinO), where O is a random variable evenly distributed between 0 and 2m, and s & 1. This new atomic arrangement is theo defined as the new equilibrium position. Each "nearest-neighbor" bond length is now

66

j.

O. 5

0 00 ~

0. 05

12

0. 10

I

I

0. 15

0. 20

0. 25

S2

FIG. 1. Showing Cll, C66, and C» vs s' for the randomgeometry model. The solid lines are extrapolations of the computed initial slope.

CAUCHY RELATIONS FOR CENTRAL-FORCE RANDOM NETWORKS

36

2117

1. 5

TABLE I. Elastic moduli results from simulations of a 100 && 116= 11 600 site central-force triangular network. The nearest-neighbor force constants were evenly distributed as indicated. Each modulus is an average over 20 independent

1. 0

configurations, with CI2 averaged over the x and y directions as well. I is defined as I =(o.2 — ai)/o. i, with o. i &a &a2.

1&a&1.2

0. 5

1

0. 0

=0.2

Cee Cl~

&o. &2.0 I =1.0

C12

I

0. 00

Cee

0. 25

Cia 1

Showing C=(Cee —CI2)/Cee vs s for the randomgeometry model. The solid line is an extrapolation of the computed initial slope.

EMT

Modulus

Distribution

&a &200

Cee

0.475 438

0.475 433+0.000 176 0.475 435+0.000 673

0.633 12

0.633 11+0.000 88 0.633 12+0.000 93

3.519

3.522+0. 019 3.531+0.078 30.91+0.22

30.87

=199

I

Simulation

31.09+0. 31

FIG. 2.

=1.

shows C=(C66 — C~q)/C66 plotted versus s . It is clear that C ~2& C«even for the smallest values of s . The straight line is the extrapolated initial slope. The agreement with the data points is worse than in Fig. 1, because the division of two uncertain numbers always introduces more error.

IV. FLAT DISTRIBUTION OF FORCE-CONSTANTS MODEL The model studied in this section is a triangular network with the harmonic lattice potential given in Eq. (4), but with nearest-neighbor central-force constants e distributed such that a& &n &aq, with o; equally likely to take on any value in the interval. If the generalized hypothesis is correct, then this network should obey the Cauchy relations. Computational results were averaged over 20 configurations of a 100&& 116=11600 site supercell with periodic boundary conditions. Results for CI2 were also over the x and y directions for each averaged The relaxation procedure of Sec. III was configuration. used. Table I contains the results for the model for ai —1, and o.'2 —1.2, 2.0, 20, and 200. Ci~, C«, and the bulk modulus K were computed independently [approximately 850 h of central-processing-unit (CPU) time on a VAX 11/785 computer, double precision], and C~2 was calculated via the relation K= —,'(C~~+C~q). Only C~2 and C66 results are presented in the table. The standard deviation in each modulus was calculated over the 20 The elfective-medium-theory configurations. (EMT) results were obtained using the single-bond theory derived in Ref. 7, with the probability distribution P(a) = 1/ 'a2 — cx)), a) &a &cz2. The parameter I is defined as I =(a2 —a~)/a~, and so is a measure of the disorder present in the force constant distribution. When I =0, there are only a] force constants present, and so C&z must equal C66. As can be seen in Table I, for I =0. 2, C~2 and C66

0 the differ by less than 1 part in 200000, while for I two moduli differ by about 1 part in 60000. In both cases, C~q —C66 is much less than o. tq or o 66, the standard deviations in the two moduli. In the case of the larger I values, I =19 and 199, the disagreement between C~2 and C66 is rather larger (1 part in 350 and 150, respectively), but C~2 —C66 is still less than cr~q or cr66 It is important to note that for a2 »a~ (I &&1), the fiucbetween simulations of tuations independent configurations were larger both in an absolute and relative sense than for smaller I . For example, o. )2/C)2 -—0.001 for I =0.2, while for I =199, o. ]2/C~q-— 0. 01, an order of magnitude increase. However, for all four values of I considered, one can state that C ~2 and C«are equal within "experimental" error. ~

V. CORRELATED BOND-PERCOLATION MODEL

The correlated bond-percolation model was originally described in Ref. 11. A slightly modified version was used in this work. Qne assigns a random number x, to ' ' each site of a lattice, where — —, &x, & —, and x, is evenly distributed. Then for each bond connecting sites i and j, one variable defines the random x„=—,' ( x, + x, + x, —x, ), where 0&x„& 1. One then chooses x, 0&x & 1, and removes all bonds with x,-, ~x. It is not difficult to see that the actual fraction of bonds remaining after this operation is carried out is p (x), the distribution function of x treated as a random variable. Using elementary techniques one can show that -

~

~

p='

~

~

1

—(1 —x)

for

—,

3x

for

0&x &

'

&x

&

1, '

—,

I,

Figure 3 shows p versus

dp/dx at x

= '.

x. Note the discontinuity

in

—,

networks Simulations were performed on triangular with values of x ranging between 0. 1 and 0.9. However,

E. J. GARBOCZI

2118

0. 12

1. 00

0. 75

o. oa

N

Q

0. 50

0. 06

E

(3 C3

0. 25

0. 00

0. 03

CL

0. 00

0. 75

0. 50

0. 25

0 00

1. 00

0. 0

0. 2

0. 4

0. 6

FIG. 3. Showing p vs x for a correlated bond-percolation process, where x is the parameter defined in Sec. V and p =p(x) is the fraction of bonds present.

FIG. 5. Showing the parameters C=(Cl~ —C«)/Cf, t, (circle) C l2 —2C«) /C«(square) vs p. and = ( C l — and C were calculated from independent computations of Cll, K, and C«.

instead of removing bonds with x,,-)x, such a bond's force constant was changed from 1 to 0.2. This procedure avoided the numerical difriculties of the rigidity phase transition at p==', (x =0.47), but still allowed a fairly large range of elastic moduli to be explored. Fifteen configurations of a 60& 70=4200 site supercell were aver750 h of CPU aged for each value of x. Approximately time on a VAX 11/785, in double precision, were required to complete the correlated bond-percolation simulations. The relaxation procedure of Sec. III was used. Figure 4 shows K and C«plotted as a function of the bond fraction p. In Figs. 4 —6, p is the fraction of bonds with force constant o, =1.0. The solid lines are the generalized EMT derived in Ref. 9, with p (x) used instead of the bond-percolation probability. It is seen in Fig. 4 that

this EMT is correct as p~1 and p~0, as both of these cases are in the few-defect limit. When x = 1 —b (or x =b)), for any network one can choose 6 sufficiently small so that only one bond is diff'erent from the rest. The EMT is exactly correct in the one-bond-defect limit. For 0.2 &p &0. 8, the EMT is clearly incorrect, though not by any great amount. Similar behavior was also seen in the conductivity case. Figure 5 shows C=(C~q — C66)IC66 and I=(C~~ —Ci2 —2C66)/C66 versus p. Cauchy's relation is clearly violated for intermediate values of p and is obeyed as The parameter is a measure of the anand p p isotropy of the network. The infinite triangular net with all bonds the same is isotropic and should remain so as bonds are randomly changed. The small nonzero values

"

I

I

&

~1.

~0

I

2. 2

1. 0

2. 1

0. 6 D

0. 4

2. O

0. 2 0. 0

FIG. 4.

0. 0

l

l

0. 2

0. 4

Showing the bulk modulus

0. 6

! 0. 8

l

l

o. 4

0. 6

0. 8

K/C«vs

p for the data plotted in

I

0. 0

0. 2

K and shear modulus C«

p, where p is the fraction of bonds with force constant o.'=1. A fraction 1 — p of the bonds have n =0.2. vs

1. a

1. 0

FIG. 6. Fig. 4.

Showing the ratio

CAUCHY RELATIONS FOR CENTRAL-FORCE RANDOM NET%'ORKS

36

of I observed are due to the effect of a network with finite size. I is much smaller than C, so that the violation of Cauchy's relation is not a finite-size effect. Figure 6 is a plot of the ratio K/C66 versus p. For an isotropic two-dimensional network that satisfies the Cauare nonzero, we chy relations, K/C66 ——2. When C and can write

I

K /C66

—2+ C+

I,

' —,

so that the deviation of K/C66 from 2 in Fig. 6 is due almost entirely to the violation of Cauchy's relation and not finite-size anisotropy.

VI. DISCUSSION The generalized hypothesis stated in Sec. II is essentialto classify what kinds of random centralforce systems will obey the Cauchy relations. It should be emphasized that by "random" is meant random central-force constants superimposed on a centrosymmetric lattice. The three models studied above have given support to the suSciency of this hypothesis. The first model studied, the random-geometry model, showed that the condition of having atoms at the sites of a centrosymmetric (neglecting the force constants) lattice must be retained in any generalized Cauchy hypothesis. The model described in Sec. IV had nearest-neighbor force constants that were randomly and independently distributed on a centrosymmetric lattice. One should take note that every nearest-neighbor bond had a different force constant. For I =0.2 and 1.0, the Cauchy relation was found to be obeyed to 1 part in 200000 and 1 part in 60000, respectively, thus confirming the prediction of the generalized Cauchy hypothesis. Since I is a measure of the force-constant disorder, it should be emphasized that these two values of I represent a substantial perturbation of the network. It is also interesting to note that in these two cases the EMT and simulation results agreed within the same small amount. One might speculate that this EMT could be exact for this model. It is known that this EMT does not predict the correct critical threshold and exponent for models studied previously that undergo a riIt is possible that the close agreement gidity transition. with EMT in the I =0.2 and 1.0 cases was due to the fairly small range over which the force constants were distributed. The I =19 and 199 simulations were carried out to check whether for o. 2»o. the Cauchy relation would still be obeyed but with increasing disagreement between EMT and simulation. It was found that EMT and simulation results still only differed by less than a standard deviation, with this difference being of the same orly an attempt

C„—

der as C« . Limitations on computer time prohibited using more configurations on larger networks to reduce the standard deviations at the two higher I values. The computer time necessary to completely relax a network increases rapidly with I . Of the 850 CPU hours required for the simulation results given in Sec. IV, 700 h were spent just on the I =19 and 199 computations. The question of the exactness of this EMT for this model must remain open for the present. It should be noted, however, that if the Cauchy relations are obeyed by this or any random system, then some form of renormalized forceconstant description must be exact. The correlation bond-percolation model in Sec. V was chosen because it is intermediate between ordinary bond and site percolation, in the following sense: As x decreases from 1 (or increases from zero), only one bond is altered initially, so that the single-bond EMT is exact for p very close to 1 (or 0), but correlations quickly build up as more bonds are altered that cause the distribution of unaltered bonds to differ significantly from simple bond percolation. The conjecture would predict that the Cauchy relation would be obeyed for p close to 1 (or zero) with the largest disagreement near p=0. 5. This was exactly the behavior observed. In Fig. 5, the parameter ( C~2 C66 )/C66 increased as p decreased from 1, and had its maximum value at about p =0. 5. C then decreased as p decreased further, as the +=0.2 bonds became dominant and the o:=1.0 bonds became the "de~

~

C:

"

fects. The next step in this work would seem to be an attempt to formulate an analytical proof that the generalized Cauchy hypothesis given here does indeed lead to the Cauchy relations. This is not an easy task. Lacking a general proof, it would be useful to take the flat distribution of force constants model described in Sec. IV and attempt to show analytically that the Cauchy relations hold exactly for this specific case. Such a proof might be simply generalized. Of course, even if one could formulate a proof using the full hypothesis presented in this paper, this would still not necessarily eliminate the possibility of an even more general hypothesis being su%cient for the CauAddition careful numerical work on chy relations. different random systems might find a clear counterexample to support such further generality. A. CENO%I. EDGMENTS

~

In the Voigt notation, the six independent combinations of x, y, and z are written as xx=1, yy=2, zz=3, yz=4, zx =5, and xy— = 6. In this notation C;,kl is written as C „, m, n =1,6. The symmetry (ij)~(kl) is equivalent to saying that C „ is a symmetric second-rank tensor. A. Cauchy, Exercises Math. 3, 213 (1828). B. Saint-Venant, Memoires des Savants e'trangers 14, (1855).

2119

I would like to thank M. F. Thorpe for many helpful conversations about the Cauchy relations, and E. W. Behrens for a helpful conversation about elasticity theory.

~A. E. H. Love, A Treatise on the Mathematical Theory of Elasti city (Cambridge University Press, Cambridge, 1934). sM. Born and K. Huang, Dynamical Theory of Crystal Lattices (Clarendon, Oxford, 1966), Chap. 3. 6J. H. Weiner, Statistical Mechanics of Elasticity (Wiley, New York, 1983), Chap. 9; I. Stakgold, Q. Appl. Math. 8, 169

(1950-1951).

2120 7S. Feng, M. F. Thorpe, and E. J. Garboczi, 276 (1985). 8E. J. Garboczi and M. F. Thorpe, Phys. Rev. E. J. Garboczi and M. F. Thorpe, Phys. Rev. ' M. F. Thorpe and E. J. Garboczi, Phys.

E. J. GARBOCZI Phys. Rev. B

3I,

B 31, 7276 (1985). B 33, 3289 (1986). Rev. B 35, 8579

36

(1987). ''Scott Kirkpatrick, Rev. Mod. Phys. 45, 574 (1973).

' This transition

smaller systems.

was seen in earlier unpublished

work done on