Combining self-consistent and numerical methods for the calculation of elastic fields and effective properties of 3D-matrix composites with periodic and random microstructures

International Journal of Engineering Science 49 (2011) 420–442

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Combining self-consistent and numerical methods for the calculation of elastic fields and effective properties of 3D-matrix composites with periodic and random microstructures S. Kanaun a,⇑, E. Pervago b a b

ITESM-CEM, Atizapan, Edo de México, México 52926, Mexico Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas 152, Mexico

a r t i c l e

i n f o

Article history: Received 9 November 2010 Received in revised form 7 January 2011 Accepted 23 January 2011 Available online 18 February 2011 Keywords: Heterogeneous medium Integral equations Gaussian approximating functions Toeplitz matrix Fast fourier transform Effective elastic properties

a b s t r a c t The work is devoted to the calculation of static elastic fields in 3D-composite materials consisting of a homogeneous host medium (matrix) and an array of isolated heterogeneous inclusions. A self-consistent effective field method allows reducing this problem to the problem for a typical cell of the composite that contains a finite number of the inclusions. The volume integral equations for strain and stress fields in a heterogeneous medium are used. Discretization of these equations is performed by the radial Gaussian functions centered at a system of approximating nodes. Such functions allow calculating the elements of the matrix of the discretized problem in explicit analytical form. For a regular grid of approximating nodes, the matrix of the discretized problem has the Toeplitz properties, and matrix–vector products with such matrices may be calculated by the fast fourier transform technique. The latter accelerates significantly the iterative procedure. First, the method is applied to the calculation of elastic fields in a homogeneous medium with a spherical heterogeneous inclusion and then, to composites with periodic and random sets of spherical inclusions. Simple cubic and FCC lattices of the inclusions which material is stiffer or softer than the material of the matrix are considered. The calculations are performed for cells that contain various numbers of the inclusions, and the predicted effective constants of the composites are compared with the numerical solutions of other authors. Finally, a composite material with a random set of spherical inclusions is considered. It is shown that the consideration of a composite cell that contains a dozen of randomly distributed inclusions allows predicting the composite effective elastic constants with sufficient accuracy. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction We consider elastic fields in heterogeneous and composite media and their effective elastic properties. Conventional numerical methods applied to the solutions of these problems are the Finite Element Method (FEM) (see, e.g., Gusev, 1999; Segurado & Llorca, 2002; Zohdi & Wriggers, 2005), the fast fourier transform (FFT) technique applied to the solution of the volume integral equations of elasticity for heterogeneous media (Moulinec & Suquet, 1994, 1998), and the multipole expansion method (Davis, Hatch, Yener, & Chompooming, 1993; Kushch, 1997). In all these methods, the so-called representative volume element (RVE) of a composite material with periodic boundary conditions in its sides is considered. This ⇑ Corresponding author. E-mail addresses: [email protected] (S. Kanaun), [email protected] (E. Pervago). 0020-7225/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2011.01.001

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

421

volume contains a sufficiently large number of inclusions in order to be macroscopically homogeneous. The number of the inclusions that should be contained in RVE is not fully clear, but numerical results of Gusev (1999) and Segurado and Llorca (2002) show that this number is not very large: ‘‘a couple of dozens’’ as stated in Gusev (1999). On the other hand, several self-consistent methods have been developed for the calculation of effective elastic constants of heterogeneous materials (see Kanaun & Levin, 2008 for overview). These methods reduce the problem of interaction between many inhomogeneities to a problem for one inhomogeneity (the one-particle problem) using heuristic hypotheses about stress–strain state of each inclusion in the composite. In the effective field method (EFM), it is assumed that every inclusion is subjected to an effective external field that is the same for all the inclusions. In the effective medium method (EMM), every inclusion is supposed to be embedded in the medium with the effective properties of the composite. If the one-particle problem may be solved analytically, the self-consistent methods yield analytical equations for the effective elastic constants. These equations are in good agreement with experimental data and numerical solutions provided the volume fraction of the inclusions is not very large. Unfortunately, analytical solutions of the one-particle problems exist only for the ellipsoidal inclusions. The application of self-consistent methods may be substantially extended to inclusions of arbitrary shapes if numerical methods are used for the solution of the one-particle problems. Also, one-particle problems may then be formulated as the problems for a typical cell of the composite that is similar to RVE. This cell contains a finite number of inclusions and is treated as one particle. As a result, the homogenization problem may be solved with much higher accuracy compared to the conventional self-consistent methods. In addition, detailed elastic fields in a typical cell of the composites can be found. In the present work, a self-consistent effective field method (EFM) is developed for the calculation of the elastic fields in 3D-matrix composites and prediction of their effective elastic properties. The ‘‘one-particle problem’’ is formulated as a problem for a cell of a composite that contains several isolated inclusions. Such a cell is subjected to an effective external field that takes into account the presence of the surrounding inclusions. In solving the 3D-elasticity problem for a cell, a numerical method based on integral equations for elastic fields in heterogeneous materials, radial Gaussian approximating functions, and the fast fourier transform technique are used. The results for composites with periodic and random lattices of spherical inclusions are compared with numerical solutions of other authors available in the literature. It is shown that the consideration of a cell that contains about a dozen of inclusions yields the effective elastic constants of regular and random composites with sufficient accuracy.

2. Integral equations for elastic fields in a homogeneous medium with a set of isolated inclusions Let an infinite homogeneous medium with the stiffness tensor C0 contain a number of isolated inhomogeneities that occupy regions Vk (k = 1, 2, . . .), with stiffness tensors Ck(x), where x(x1, x2, x3) is a point of 3D-space. The medium is subjected to an external strain field e0(x) (or stress field r0(x)), and the objective is to calculate elastic fields in the medium with the inclusions. This problem may be reduced to volume integral equations for the strain e(x) or stress r(x) tensors in the heterogeS neous medium. Let V = kVk, and V(x) be the characteristic function of the region V occupied by the inclusions: V(x) = 1 if x 2 V, V(x) = 0 if x R V. The strain field e(x) in the medium with the inclusions satisfies the following integral equation (see, e.g., Kunin, 1983; Kanaun & Levin, 2008):

eij ðxÞ þ

Z

0

K ijkl ðx  x0 ÞC 1klmn ðx0 Þemn ðx0 ÞVðx0 Þdx ¼ e0ij ðxÞ:

ð1Þ

Here C1(x) = C(x)  C0, C(x) = C(k)(x) when x 2 Vk, and C(x) = C0 when x R V. The kernel K(x) of the integral operator in this equation is the second derivative of the Green function G(x) of the homogeneous host medium C0

  K ijkl ðxÞ ¼  @ i @ k Gjl ðxÞ ðijÞðklÞ ;

ð2Þ

where G(x) is the solution of the equation:

@ i C 0ijkl @ k Glm ðxÞ ¼ djm dðxÞ;

ð3Þ

where d(x) is Dirac’s delta-function. The subscripts in parentheses mean symmetrization. A similar equation may be written for the stress field r(x):

rij ðxÞ 

Z

0

Sijkl ðx  x0 ÞB1klmn ðx0 Þrmn ðx0 ÞVðx0 Þdx ¼ r0ij ðxÞ:

ð4Þ

Here B1(x) = B(x)  B0, B(x) = C1(x), and B0 = (C0)1, r0(x) is an external stress field applied to the medium, and

Sijkl ðxÞ ¼ C 0ijmn K mnrs ðxÞC 0rskl  C 0ijkl dðxÞ:

ð5Þ

Note that the functions K(x) and S(x) in Eqs. (2) and (5) behave as jxj3 when jxj ? 0, and therefore, the integral operators K and S in Eqs. (1) and (4) are singular. Actions of these operators on a continuous tensor function fij(x) with a finite support are defined by the equations (see Kunin, 1983; Mikhlin, 1965):

422

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Fig. 1. A cell W0 of the composite material; cuboid W is the region where the numerical solution is constructed.

ðKf Þij ðxÞ ¼ Aijkl fkl ðxÞ þ v :p: ðSf Þij ðxÞ ¼ Dijkl fkl ðxÞ þ v :p:

Z Z

0

ð6Þ

Sijkl ðx  x0 Þfkl ðx0 Þdx :

ð7Þ

K ijkl ðx  x0 Þfkl ðx0 Þdx ; 0

R Here v :p: . . . dx is the Cauchy principal value of the integral that exists for K(x) and S(x) in the forms (2) and (5). The constant tensors A and D have the forms:

A¼

1 4p

Z

e KðkÞd X;

X1

D¼

1 4p

Z

e SðkÞdX:

ð8Þ

X1

e Here KðkÞ and e SðkÞ are the fourier transforms of the functions K(x) and S(x), and X1 is the surface of a unit sphere in the k-space of the fourier transforms. Eqs. (6) and (7) define regularizations of the singular integrals in Eqs. (1) and (4). If the e medium is isotropic, the function KðkÞ takes the form:

e ijkl ðkÞ ¼ K

Z

K ijkl ðxÞ expðik  xÞdx ¼

1 h

l0

i E5ijkl ðmÞ  j0 E6ijkl ðmÞ ;

m¼

k ; jkj

ð9Þ

where k is the vector parameter of the fourier transform, k  x is the scalar product of the vectors k and x,

j0 ¼

k0 þ l0 k0 þ 2l0

ð10Þ

and k0 and l0 are the Lame constants of the host medium. E5(m) and E6(m) are the elements of the tensor basis Ek(m) (k = 1, 2, 3, 4, 5, 6) introduced in Kunin (1983) for presentation of fourth-rank tensors:

E1ijkl ¼ dik djl jðijÞðklÞ ; E4ijkl ¼ mi mj dkl ;

E2ijkl ¼ dij dkl ;

E3ijkl ¼ dij mk ml ;

E5ijkl ¼ mi mk djl jðijÞðklÞ ;

E6ijkl ¼ mi mj mk ml :

ð11Þ

The fourier transform e SðkÞ of the kernel S(x) in Eq. (5) has the form

h i e S ijkl ðkÞ ¼ 2l0 P 1ijkl ðmÞ þ ð2j0  1ÞP2ijkl ðmÞ ;

ð12Þ

P1 ðmÞ ¼ E1  2E5 ðmÞ þ E6 ðmÞ;

ð13Þ

P2 ðmÞ ¼ E2  E3 ðmÞ  E4 ðmÞ þ E6 ðmÞ:

Let fij(x) be a smooth tensor-function with fourier transform e f ij ðkÞ. If e f ij ðkÞ is bounded and tends to zero at infinity faster 3 than jkj , the actions of the convolution operators K and S on such a function are defined by the equations:

ðKf Þij ðxÞ ¼

Z

1 3

ð2pÞ

e ijkl ðkÞef kl ðkÞ expðik  xÞdk; K

ð14Þ

423

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

ðSf Þij ðxÞ ¼

Z

1 ð2pÞ3

e S ijkl ðkÞef kl ðkÞ expðik  xÞdk:

ð15Þ

The integrals in the right-hand sides of these equations exist as ordinary integrals. Since the strain e(x) and stress r(x) tensors under the integrals in Eqs. (1) and (4) are multiplied by the function V(x), the elastic fields inside the inclusions (in the region V) are the principal unknowns of the problem. The strain and stress fields in the medium are reconstructed from the same Eqs. (1) and (4) if the fields inside the inclusions are known. Another important implication of Eqs. (1) and (4) is that for the numerical solution of these equations, any region W that includes V may be considered. In particular, by the solution of the problem for a medium with a finite number of inclusions, a cuboid W that contains all these inclusions may be considered (see Fig. 1). Unique solutions of Eqs. (1) and (4) exist if the tensor of the elastic constants C(x) and the inverse tensor C1(x) = B(x) do not degenerate inside V (see Mikhlin, 1965). 3. Numerical solution of the integral Eqs. (1) and (4) 3.1. Discretization of Eq. (1) by the Gaussian approximating functions For the numerical solution of the integral Eqs. (1) and (4), they should be firstly discretized using an appropriate class of approximating functions. In this work, the radial Gaussian functions are used for this purpose. According to Maz’ya and Schmidt (2007), a bounded smooth function u(x) defined on 3D-space can be presented in the form of the following series:

uðxÞ  uh ðxÞ ¼

X

ðrÞ

ðrÞ

u uðx  x Þ;

uðxÞ ¼

r

1 ðpHÞ3=2

exp 

jxj2 2

Hh

! ð16Þ

:

Here x(r) (r = 1, 2, . . .) are the nodes of a regular grid, h is the grid step, u(r) = u(x(r)) is the value of function u(x) at the node x = x(r), and H is a non-dimensional parameter of the order 1. In Maz’ya and Schmidt (2007), Eq. (16) is called the ‘‘approximate approximation’’ because its error does not vanish when h ? 0. However, the non-vanishing part of the error has the order of exp (p2H) and may be neglected in practical calculations. Let us start with Eq. (1) for strains and find its solution in the form similar to (16):

eij ðxÞ 

N X

eijðsÞ uðx  xðsÞ Þ:

ð17Þ

s¼1

Here x(s) (s = 1, 2, . . . , N) are the nodes of a regular grid that covers a cuboid W (Fig. 1) that contains the region V occupied by the inclusions, N is the total number of the nodes in W, e(s) are unknown coefficients of the approximation. After substituting (17) into the integral in Eq. (1) we obtain the following equation:

eij ðxÞ þ

N X

ðsÞ 0 Pijkl ðx  xðsÞ ÞC 1ðsÞ klmn emn ¼ eij ðxÞ;

1ðsÞ

C ijkl ¼ C 1ijkl ðxðsÞ Þ;

ð18Þ

s¼1

where the tensor Pijkl(x) has the form

Pijkl ðxÞ ¼ ¼

Z

0

K ijkl ðx  x0 Þuðx0 Þdx ¼ h

Z h

3

ð2pÞ3 l0

1 ð2pÞ3

Z

e ijkl ðkÞ u e ðkÞ exp ðik  xÞdk K

! 2 i Hh jkj2 E5ijkl ðmÞ  j0 E6ijkl ðmÞ exp   ik  x dk: 4

ð19Þ

e Here the definition (14) of the operator K and Eq. (9) for KðkÞ are used,

! 2 Hh jkj2 e ðkÞ ¼ h exp  : u 4 3

ð20Þ

1ðsÞ

Note that in Eq. (18), C ijkl ¼ 0 if x(s) R V. After introducing a spherical coordinate system in the k-space and integrating firstly over the unite sphere, and then, over the radius jkj, the integral in the right hand side of Eq. (19) is calculated explicitly, and the tensor P(x) takes the form

P¼

1 n

l0

h io w1 E1 þ U0 E5 ðnÞ  j0 w2 ðE2 þ 2E1 Þ þ U1 ðE3 ðnÞ þ E4 ðnÞ þ 4E5 ðnÞÞ þ U2 E6 ðnÞ ;

U0 ¼ w0  3w1 ;

U1 ¼ w1  5w2 ;

U2 ¼ w0  10w1 þ 35w2 :

n¼

x ; jxj

ð21Þ ð22Þ

Here Ek(n) (k = 1, . . . , 6) are the elements of the tensor basis (11). The three scalar functions w0, w1, w2 in Eqs. (21) and (22) depend on jxj and have the forms:

424

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

wa ¼ wa

w0 ðzÞ ¼

  jxj ; h 1 ðpHÞ3=2

a ¼ 0; 1; 2;

ð23Þ

 2 z ; exp  H

ð24Þ

w1 ðzÞ ¼

 2  pffiffiffiffiffiffiffi    z z pffiffiffiffi 2z exp  þ pHErf pffiffiffiffi ; H 4p3=2 z3 H H

ð25Þ

w2 ðzÞ ¼

 2  pffiffiffiffiffiffiffi   z z p ffiffiffiffi þ 6 p H z exp  p ð 3H þ 2z ÞErf : 16p2 z5 H H

ð26Þ

1

1

Here Erf(z) is the probability integral:

2 Erf ðzÞ ¼ pffiffiffiffi

Z

p

z

2

et dt:

ð27Þ

0

The system of linear algebraic equations for the unknowns e(s) in Eq. (17) follows from Eq. (18) if the latter is satisfied at all the nodes (the collocation method). Note that if the nodes of the approximation (17) compose a cubic grid, the coefficients of the approximation e(s) coincide with the values of the function e(x) in the corresponding nodes (e(s) = e(x(s))) (Mikhlin, 1965). As a result, we obtain the system of linear algebraic equations for the coefficients e(s) in the form

eðrÞ ij þ

N X

1ðsÞ ðsÞ 0ðrÞ PðrsÞ ijkl C lmn emn ¼ eij ;

r ¼ 1; 2; . . . ; N;

ð28Þ

s¼1 ðrÞ Pðr;sÞ  xðsÞ Þ; ijkl ¼ Pijkl ðx

0ðrÞ 1ðrÞ C ijkl ¼ C 1ijkl xðrÞ ; eij ¼ e0ij ðxðrÞ Þ:

ð29Þ

This system may be written in the canonical form as follows:

ðI þ BÞX ¼ F;

ð30Þ

where I is the unit matrix of the dimensions 6N  6N, and the vectors of the unknowns X and the right-hand side F are

X ¼ jX ð1Þ ; X ð2Þ ; . . . ; X ð6NÞ jT ;

X ðrÞ ¼

8 ðrÞ e11 ; > > > > > > > eðrNÞ ; > 22 > > > > > < ðr2NÞ

e33

F ¼ jF ð1Þ ; F ð2Þ ; . . . ; F ð6NÞ jT ;

r6N N < r 6 2N; ; 2N < r 6 3N;

> > eðr3NÞ ; 3N < r 6 4N; > 12 > > > > ðr4NÞ > > e13 ; 4N < r 6 5N; > > > > : ðr5NÞ e23 ; 5N < r 6 6N;

F ðrÞ ¼

8 0ðrÞ e11 ; > > > > > > > e0ðrNÞ ; > 22 > > > > > < 0ðr2NÞ

e22

ð31Þ r 6 N; N < r 6 2N;

; 2N < r 6 3N;

0ðr3NÞ > > e12 ; 3N < r 6 4N; > > > > > 0ðr4NÞ > > e13 ; 4N < r 6 5N; > > > > : 0ðr5NÞ e23 ; 5N < r 6 6N:

ð32Þ

Here j  jT is the transposition operation. The matrix B in Eq. (30) has the dimensions 6N  6N and consists of 36 sub-matrices bpq of the dimensions N  N,

B ¼ jbpq j;

p; q ¼ 1; 2; . . . ; 6;

ð33Þ

ðmkÞ

¼ Pppij C ijqq ;

ðmkÞ

1ðkÞ

p; q ¼ 1; 2; 3;

ðmkÞ

ðmkÞ

1ðkÞ

b5q ¼ P13ij C ijqq ;

ðmkÞ

ðmkÞ

1ðkÞ

bq5 ¼ Pqqij C ij13 ;

ðmkÞ

ðmkÞ

1ðkÞ

b45 ¼ P12ij C ij13 ;

ðmkÞ

ðmkÞ

1ðkÞ

b55 ¼ P13ij C ij13 ;

ðmkÞ

ðmkÞ

1ðkÞ

b65 ¼ P23ij C ij13 ;

bpq

b4q ¼ P12ij C ijqq ; bq4 ¼ Pqqij C ij12 ; b44 ¼ P12ij C ij12 ; b54 ¼ P13ij C ij12 ; b64 ¼ P23ij C ij12 ;

ð34Þ

ðmkÞ

ðmkÞ

1ðkÞ

b6q ¼ P23ij C ijqq ;

ðmkÞ

ðmkÞ

1ðkÞ

ðmkÞ

ðmkÞ

1ðkÞ

bq6 ¼ Pqqij C ij23 ;

ðmkÞ

ðmkÞ

1ðkÞ

ðmkÞ

ðmkÞ

1ðkÞ

b46 ¼ P12ij C ij23 ;

ðmkÞ

ðmkÞ

1ðkÞ

ð37Þ

ðmkÞ

ðmkÞ

1ðkÞ

b56 ¼ P13ij C ij23 ;

ðmkÞ

ðmkÞ

1ðkÞ

ð38Þ

ðmkÞ

ðmkÞ

1ðkÞ

b66 ¼ P23ij C ij23 :

ðmkÞ

ðmkÞ

1ðkÞ

ð39Þ

ð35Þ q ¼ 1; 2; 3;

ð36Þ

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

425

In these equations, m, k = 1, 2, . . . , N; summation with respect to repeating indices i, j from 1 to 3 is implied. These equations follow from (28) and (30). As it is seen from Eqs. (21)–(26) the elements of the matrix B in Eqs. (33)–(39) have simple analytical forms and are calculated fast. 3.2. Discretization of Eq. (4) by the Gaussian approximating functions Let us consider Eq. (4) and find its solution in the form similar to (17):

rij ðxÞ 

N X

ðsÞ ij

1

ðsÞ

r uðx  x Þ; uðxÞ ¼

s¼1

ðpHÞ3=2

exp 

jxj2 2

Hh

! ð40Þ

:

Substitution of this approximation in Eq. (4) and application of the collocation methods lead to the following linear algebraic ðsÞ system for the coefficients rij of the approximation

rðrÞ ij 

N X

ðr;sÞ 1ðsÞ 0ðrÞ Cijkl Bklmn rðsÞ mn ¼ rij ;

r ¼ 1; 2; . . . ; N;

ð41Þ

s¼1 1ðsÞ 0ðrÞ ðrÞ rðrÞ Bijkl ¼ B1ijkl ðxðsÞ Þ; rij ¼ r0ij ðxðrÞ Þ; ij ¼ rij ðx Þ;

ðrÞ Cðr;sÞ  xðsÞ Þ; ijkl ¼ Cijkl ðx

Cijkl ðxÞ ¼

Z

ð42Þ 0

Sijkl ðx  x0 Þuðx0 Þdx :

ð43Þ

The last integral is calculated explicitly, similar to the integral in Eq. (19), and the tensor C(x) takes the form:

CðxÞ ¼ 2l0

6 X m¼1

gm

  jxj m E ðnÞ; h

g 1 ¼ ðw0  2w1 Þ þ 4j0 w2 ;

n¼

x ; jxj

ð44Þ

g 2 ¼ ð1  2j0 Þðw0  2w1 Þ þ 2j0 w2 ; 1 g 5 ¼  ðU0  16j0 U1 Þ; 2

g 3 ¼ g 4 ¼ ð1  2j0 ÞU0 þ 2j0 U1 ;

g 6 ¼ 2j0 U2 :

ð45Þ ð46Þ

Here the functions Uk = Uk(jxj/h) and wk = wk(jxj/h) (k = 0, 1, 2) are defined in Eqs. (22)–(26). Similar to Eq. (28), the discretized Eq. (41) may be presented in the canonical form (30). In this case, the vector of unknowns X is composed from the values of the stress tensor rij(x) at the nodes, the vector of the right hand side F consists of the components of the external stress field r0ij ðxÞ at the nodes similar to Eqs. (31) and (32). The matrices bpq in Eq. (33) ðr;sÞ ðr;sÞ 1ðsÞ 1ðsÞ are defined in Eqs. (34)–(39), where Ppqij should be changed for Cpqij , and C ijkl for Bijkl . 3.3. Numerical solution of the system (30) It follows from Eqs. (21)–(25), (26), (35)–(39) that B in Eq. (30) is a non-sparse matrix which dimensions may be very large if high accuracy of the solution is required. For the solution of linear algebraic systems with such matrices, only iterative methods are efficient. For instance, if the minimal residue method (see, e.g., Press, Flannery, Teukolsky, & Vetterling, 1992) is used, the nth iteration X(n) of the solution of Eq. (30) is calculated as follows:

XðnÞ ¼ Xðn1Þ  aYðn1Þ ;

a¼

Yðn1Þ  ðI þ BÞY ðn1Þ ðI þ BÞYðn1Þ  ðI þ BÞYðn1Þ

YðnÞ ¼ Yðn1Þ  aðI þ BÞYðn1Þ

;

ð47Þ ð48Þ

with the initial values X(0), Y(0) of the vectors X and Y

Xð0Þ ¼ F;

Yð0Þ ¼ BF:

ð49Þ

Thus, the vector Y(n1) is to be multiplied by the matrix B at every step of the iteration process. For non-sparse matrices of large dimensions, calculation of such a product is an expensive computational operation. If, however, a regular grid of approximating nodes is used, the volume of calculations is reduced substantially. Let us consider the product BY in detail. For the matrix B that corresponds to Eqs. (29)–(35), this product is a combination of the following sums ðrÞ

Pij ¼

N X s¼1

Pijkl ðxðrÞ  xðsÞ ÞZ ðsÞ kl ;

ðsÞ

1ðsÞ

ðsþðj1ÞNÞ Z kl ¼ C klmn Y mn ;

ð50Þ

426

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Fig. 2. The correlation function W(r).

where the tensor Pijkl(x) is defined in Eq. (21). For a regular node grid with a step h, the coordinates of every node x(r) can be presented in the form:

 ðmÞ ðnÞ ðpÞ xðrÞ ¼ x1 ; x2 ; x3 ; ðmÞ

x1

ð51Þ ðnÞ

¼ L1 þ hðm  1Þ; x2 ¼ L2 þ hðn  1Þ;

ðpÞ

x3 ¼ L3 þ hðp  1Þ:

ð52Þ

Here m, n, p = 1, 2, . . . are integers, L1, L2, L3 are minimal values of the node coordinates in the cuboid W which sides are parallel to the axes x1, x2, x3. Thus, the position of every node may be defined by 3 integers (m, n, p). Connection between the one and three-index numerations of the nodes may be introduced by the equation

rðm; n; pÞ ¼ m þ N1 ðn  1Þ þ N1 N2 ðp  1Þ; 1 6 m 6 N1 ;

1 6 n 6 N2 ;

ð53Þ

1 6 p 6 N3 :

Here N1, N2, N3 are the number of the nodes along the corresponding sides of the cuboid W, N = N1N2N3. In the three-index numeration, the sum (50) is presented as follows: rðm;n;pÞ

Pij

¼

N1 X N2 X N3 X





ðqÞ ðnÞ ðsÞ ðpÞ ðtÞ ðq;s;tÞ Pijkl xðmÞ Z kl ; 1  x1 ; x 2  x 2 ; x3  x 3

ðq;s;tÞ

Z kl

rðq;s;tÞ

¼ Z kl

;

ð54Þ

q¼1 s¼1 t¼1

where





ðqÞ ðnÞ ðsÞ ðpÞ ðtÞ Pijkl xðmÞ ¼ Pijkl ðhðm  qÞ; hðn  sÞ; hðp  tÞÞ: 1  x 1 ; x 2  x2 ; x 3  x3

It is seen from the last equation that the object Pijkl



ð55Þ

ðmÞ ðqÞ ðnÞ ðsÞ ðpÞ ðtÞ x1  x1 ; x2  x2 ; x3  x3 has the Toeplitz structure: it depends on

the differences of the indices: m  q, n  s, and p  t. As a result, the fourier transform technique can be used for the calculation of the triple sums in Eq. (54), and therefore, of the matrix–vector products. Application of the FFT algorithms for calculation of these sums essentially accelerates the iterative process (47)–(49). In addition, one has to keep in the computer ðr;sÞ

memory not all the matrices Pijkl of the dimensions (N1N2N3)  (N1N2N3) but only one column and one row of every such a matrix: the object that has the dimensions (2N1  2N2  2N3) and is calculated once for all the iteration process. The details of the FFT algorithm for the calculation of the matrix–vector products and triple sums in (54) are described in Golub and Van Loan (1996), and Kanaun (2009). 3.4. An isolated layered inhomogeneity in a homogeneous elastic medium We now apply the described method to the calculation of elastic fields in a layered spherical inclusion of radius a embedded in a homogeneous and isotropic elastic medium with Young’s modulus E0 and Poisson ratio m0. In the example considered here, the inhomogeneity consists of a central kernel in the region jxj/a 6 0.5 with the Young modulus E = 0.2E0, and a spherical layer in the region 0.5 < jxj/a 6 1 with the Young modulus E = 0.5E0. The Poisson ratios of the medium and the inclusion are the same: m = m0 = 0.3. For the uniaxial constant external strain field e0ij ¼ e0 di1 dj1 , the corresponding distributions of the component e11 of the strain tensor along the x1 and x2-axes inside the inclusion are presented in Fig. 3. While the function e 11(x1, 0, 0) is on the right-hand side in this figure, the function e11(0, x2, 0) is on its left-hand side. Eq. (1) for

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

427

Fig. 3. The dependence of the component e11 of the elastic strain field inside a soft spherical inclusion with step-wise change of the elastic properties along the radius. The medium is subjected to a uniaxial strain field in the direction of the x1-axis.

strains and cubic node grids with the steps h/a = 0.1 (N = 9261), h/a = 0.05 (N = 68,921), and h/a = 0.01 (N = 8,120,601) were used for the construction of the graphs in this figure. The bold line in Fig. 3 is exact distribution of the strain tensor component e11 taken from Kanaun and Levin (2008). The distribution of the component r11 of the stress tensor in a layered spherical inclusion having Young’s moduli larger than the Young’s modulus of the medium are presented in Fig. 4. In this case, E = 10E0 when jxj/a 6 0.5 and E = 5E0 when 0.5 < jxj/a 6 1, m = m0 = 0.3. The figure presents the function r11(x1, 0, 0) and r11(0, x2, 0) for the applied uniaxial external stress field r011 . For the numerical solutions, Eq. (4) for stresses was used with the node grid steps h/a = 0.1, h/a = 0.02 (M = 1,030,301), and h/a = 0.01. The bold line in this figure is the exact distribution of the component r11 for this case (see Kanaun & Levin, 2008). Figs. 3 and 4 demonstrate convergence of the numerical solution to the exact distribution of the elastic fields inside the inclusions when the step h of the node grid decreases. The iterative scheme (47)–(49) converges for every finite values of the elastic constants of the inclusion and the matrix. The number of the iterations increases with increase of the elastic contrast as well as the number N of the approximating nodes. Note that if the contrast is not very large, the well-known conjugate gradient method (see, e.g., Press et al., 1992) applied to the solution of Eq. (30) converges faster than the minimal residue method (47)–(49). But for large contrasts in the properties, the conjugate gradient method diverges, meanwhile the minimal residue method keeps converging. Note that for the inclusions which Young’s modulus is smaller than Young’s modulus of the matrix, the iteration scheme based on Eq. (1) for strains converges faster than the same scheme based on the Eq. (4) for stresses. Meanwhile for the inclusions that are stiffer than the matrix, the numerical scheme based on Eq. (4) turns out to be more efficient.

Fig. 4. The dependence of the component r11 of the elastic stress field inside a stiff spherical inclusion with step-wise change of the elastic properties along the radius. The medium is subjected to a uniaxial stress field in the direction of the x1-axis.

428

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Fig. 5. The dependence of the relative bulk modulus K⁄/K0 of the composite with a simple cubic lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 27 inclusions; h is the step of the cubic node grid; + are the data for the composite with a cubic lattice of absolutely rigid spheres presented in Nunan and Keller (1984).

Fig. 6. The dependence of the relative shear modulus l⁄/l0 of the composite with a simple cubic lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 27 inclusions; h is the step of the cubic node grid; + are the data for the composite with a cubic lattice of absolutely rigid spheres presented in Nunan and Keller (1984).

4. The effective field method 4.1. The effective external field acting on a typical cell Let V(x) be a characteristic function of a homogeneous in space random set of inclusions. We consider a finite number of such inclusions in a spherical region W0 of the radius r0 (a typical cell of the composite, see Fig. 1). The characteristic function of the region occupied by the inclusions in W0 is denoted as V0(x). The stress field in the region W0 satisfies the equation that follows from (4):

rij ðxÞ 

Z W0

0

Sijkl ðx  x0 ÞB1klmn ðx0 Þrmn ðx0 ÞV 0 ðx0 Þdx ¼ rij ðxÞ;

ð56Þ

where

rij ðxÞ ¼ r0ij þ

Z W1

0

Sijkl ðx  x0 ÞB1klmn ðx0 Þrmn ðx0 ÞV 1 ðx0 Þdx :

ð57Þ

Here W1 is the complement of W0 to the entire space, V1(x) is the characteristic function of the region occupied by the inclusions in W1. The stress field rij ðxÞ in the right-hand side of Eq. (56) has sense of the external field that acts on the region W0. Let us fix the set of the inclusions inside W0 and average the field rij ðxÞ over the ensemble realizations of the random field of

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

inclusions in W1. The corresponding mean value of After averaging Eq. (57) we obtain

D

E

r ij ðxÞ ¼ rij ðxÞjW 0 ¼ r0ij þ

Z

429

rij ðxÞ is called the effective field acting on W0 and is denoted as rij ðxÞ.

D E 0 Sijkl ðx  x0 Þ B1klmn ðx0 Þrmn ðx0 ÞV 1 ðx0 ÞjW 0 dx :

ð58Þ

Here h  jW0i is the average over the ensemble realizations of the random set of inclusions in the region W1 subject to the condition that the set of the inclusions in W0 is fixed. The area of integration in the last integral may be spread over entire space since the function V1(x) is equal to zero in the region W0. We assume that the random function B1klmn ðx0 Þrmn ðx0 Þ is statistically independent on the spatial positions of inclusions. Then, the average under the integral in Eq. (58) takes the form

D

D E E B1klmn ðx0 Þrmn ðx0 ÞV 1 ðx0 ÞjW 0 ¼ p B1klmn ðx0 Þrmn ðx0 Þ Wðx0 ; x0 Þ;

ð59Þ

where

Wðx0 ; x0 Þ ¼

E 1D V 1 ðx0 ÞjW 0 : p

ð60Þ

Here W(x0, x) is a specific correlation function that is the normalized density of probability of finding a point x0 inside an inclusion in the region W1 if the point x0 is the center of the region W0, p = hV(x)i is the volume fraction of the inclusions. For a homogeneous and isotropic random set of inclusions, the function W(x0, x0 ) depends only on the distance between x0 and x0 (W(x0, x0 ) = W(jx0  x0 j)) and has the properties

Wðjx0  x0 jÞ ¼ 0;

jx0  x0 j 6 r 0 ;

Wðjx0  x0 jÞ ! 1;

when jx0  x0 j ! 1:

ð61Þ

The function W(r) is sketched in Fig. 2. For a homogeneous in space random set of inclusions and a constant external field r0, the average hB1klmn rmn i in Eq. (59) is constant, and the integral in Eq. (58) is calculated as follows:

Z

Z D E D E 0 0 Sijkl ðx  x0 Þ B1klmn ðx0 Þrmn ðx0 ÞV 1 ðx0 ÞjW 0 dx ¼ p Sijkl ðx  x0 ÞWðjx0  x0 jÞdx B1klmn rmn Z E   0D ¼ p Sijkl ðx  x0 Þ Wðjx0  x0 jÞ  1 dx B1klmn rmn D E ¼ pDijkl B1klmn rmn ; x 2 W 0 :

ð62Þ

Here we use Eq. (7) and the property of the operator S indicated in Kanaun and Levin (2008):

Z

0

Sðx  x0 Þdx ¼ 0;

ð63Þ

that holds if the external stress field r0 is fixed in the problem. The Cauchy principal value of the integral in (7) vanishes for  ij in Eq. (58) is constant in the region W0 and takes the function W possessing properties (61). As a result, the effective field r the form

D

E

r ij ¼ r0ij  pDijkl B1klmn rmn ; x 2 W 0 :

ð64Þ

The second term on the right is the average stress field induced in W0 by all the inclusions in the region W1 (outside W0). Note that this result depends neither on the radius r0 of the region W0 nor on details of behavior of the function W(jx0  x0 j) outside W0.  ij in (64) we obtain the following equation for the After changing the external field rij ðxÞ in Eq. (56) to the effective field r stress field in the region W0

rij ðxÞ 

Z W

0

D E 0 Sijkl ðx  x0 ÞB1klmn ðx0 Þrmn ðx0 Þdx þ pDijkl B1klmn rmn ¼ r0ij ;

x 2 W0:

ð65Þ

Note that if the cell W0 is sufficiently large, one can calculate the mean hB1klmn rmn i not over the region W1 but over the region W0 of a typical cell of the composite

D

E D E B1klmn rmn ¼ B1klmn rmn

V

0

¼

1

v0

Z W0

B1klmn ðxÞrmn ðxÞV 0 ðxÞdx:

ð66Þ

Here v0 is the volume of the region V0. This equation is the condition of self-consistency of the considered version of the effective field method. Finally, Eq. (65) for the stress field in the region W0 takes the form

rij ðxÞ 

Z W0

  p 0 Sijkl ðx  x0 Þ þ Dijkl B1klmn ðx0 Þrmn ðx0 ÞV 0 ðx0 Þdx ¼ r0ij ;

v0

x 2 W0:

ð67Þ

430

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Thus, the problem is reduced to the construction of the stress field r(x) inside the region W0 that contains a finite number of isolated inclusions of arbitrary shapes and properties. The solution of this equation may be found only numerically, and an appropriate numerical method is developed in Section 3. 4.2. The tensor of the effective elastic moduli Averaging Eqs. (1) and (2) over the ensemble realizations of the random field of inclusions yields





eij ðxÞ ¼ e0ij þ p

Z

rij ðxÞ ¼ r0ij þ p

Z

D E 0 K ijkl ðx  x0 Þdx C 0klmn B1mnrs rrs ;

ð68Þ

D E 0 Sijkl ðx  x0 Þdx B1klmn rmn ;

ð69Þ

where we take into account that the average hB1mnrs ðxÞrrs ðxÞi does not depend on x, and that the following equation holds

 C1 e ¼ C  C0 C1 r ¼ C0 B1 r:

ð70Þ

0

If the external stress field r is given, we find from Eqs. (69) and (63) that

hrðxÞi ¼ r0 :

ð71Þ

It follows from Eq. (63) and definition (5) of S(x) that the function K(x) in (68) has the property

Z

 1 0 ¼ B0 : Kðx  x0 Þdx ¼ C0

ð72Þ

As a result, from (68) we obtain the following equation for the average strain field he(x)i

D



E

eij ðxÞ ¼ e0ij þ p B1ijkl rkl ; e0ij ¼ B0ijkl r0kl :

ð73Þ

Because of linearity of the problem, the mean hB1klmn rmn i is a linear function of the constant external field

D E B1klmn rmn ¼ M ijkl r0kl ;

r0ij ð74Þ

where M is a fourth-rank tensor that depends on the microstructure of the composite. Then, from Eqs. 73, 74, and 71, we find

hei ¼ ðB0 þ pMÞr0 ¼ B hri;

ð75Þ

B ¼ B0 þ pM;

ð76Þ

where B⁄ is the effective elastic compliance tensor of the composite that connects the mean values of the strain and stress tensors. Thus, for the solution of the homogenization problem (calculation of the tensor B⁄), one has to solve the integral Eq. (67) for the stress field inside the region W0, find the mean value of the product B1klmn ðxÞrmn ðxÞ over the region V0 occupied by the inclusions in W0, and then construct the tensor M from Eq. (74). 4.3. The effective field method based on the Eq. (1) for strains The effective field method may be also developed on the basis of Eq. (1) for the strain tensor e(x). This equation may be rewritten in the form

eij ðxÞ þ

Z

0

W0

where the field

K ijkl ðx  x0 ÞC 1klmn ðx0 Þemn ðx0 ÞV 0 ðx0 Þdx ¼ eij ðxÞ;

ð77Þ

eij ðxÞ on the right

eij ðxÞ ¼ e0ij 

Z

W1

0

K ijkl ðx  x0 ÞC 1klmn ðx0 Þemn ðx0 ÞV 1 ðx0 Þdx ;

ð78Þ

may be interpreted as an external strain field acting on the region W0. Averaging Eq. (78) over the ensemble realization of the random set of inclusions by the condition that the inclusions in the region W0 are fixed leads to the following equation for the effective external field e⁄ = he ⁄(x)jW0iacting on W0

D E eij ¼ e0ij þ pAijkl C 1klmn emn ;

x 2 W 0:

ð79Þ

431

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

By deriving this equation we accept that the external field e0 is fixed in the problem and hence, the following equation holds (Kanaun & Levin, 2008):

Z

0

Kðx  x0 Þdx ¼ 0:

ð80Þ

The assumption that the mean hC 1ijkl ðxÞekl ðxÞi may be calculated over the region W0 instead of the region W1 leads to the equation for the field e(x) inside W0 in the form



Z

eij ðxÞ þ

W0

K ijkl ðx  x0 Þ 

p

v0

 0 Aijkl C 1klmn ðx0 Þemn ðx0 ÞV 0 ðx0 Þdx ¼ e0ij ;

x 2 W 0:

ð81Þ

This equation is reciprocal to Eq. (67). If the external strain field e0 is fixed in the problem, the following equations hold (Kanaun & Levin, 2008):

Z

0

Kðx  x0 Þdx ¼ 0;

Z

0

Sðx  x0 Þdx ¼ C0

ð82Þ

and averaging Eqs. (1) and (2) over the ensemble realizations of the random set of inclusions yields





eij ðxÞ ¼ e0ij  p 

Z

rij ðxÞ ¼ r0ij  p

D E 0 K ijkl ðx  x0 Þdx C 1klmn emn ¼ e0ij ;

Z

D E D E 0 1 1 0 Sijkl ðx  x0 Þdx ðC 0 Þ1 klmn C mnpq epq ¼ rij þ p C ijkl ekl :

ð83Þ

ð84Þ

Here we take into account Eq. (82) and the equation

  1   1 Ce ¼  C0 B1 r ¼ C1  C0 C1 e:

ð85Þ

For linearity of the problem, the mean hC 1ijkl ekl i is a linear function of the constant external field

D

E C 1ijkl ekl ¼ Q ijkl e0kl ;

e0kl , and we can write ð86Þ

where Q is a constant tensor that depends on the microstructure of the composite. Thus, from Eqs. (83) and (84) we have  r0ij ¼ C 0ijkl e0kl







rij ðxÞ ¼ C 0ijkl þ pQ ijkl e0kl ¼ C ijkl hekl ðxÞi;

ð87Þ

where C⁄ is the effective elastic stiffness tensor of the composite

C ¼ C0 þ pQ : ⁄

ð88Þ

⁄

Because C and B in Eqs. (88) and (76) are the effective stiffness and compliance tensors of the same composite material, the equation (B⁄)1 = C⁄ holds. As a result, the tensors Q and M in Eqs. (86) and (74) are connected by the equation

 1 Q ¼ C0 MC0 E1 þ pMC0 :

ð89Þ

4.4. Discretization of the integral equations of the effective field method Using the Gaussian approximating functions Eqs. (67) and (81) of the effective field method my be discretized similar to the case of a finite number of inclusions considered in Section 3. The discretized forms of these equations are

rðrÞ ij 

N  X



1ðsÞ 1ðsÞ ðsÞ 0 Cðr;sÞ ijkl Bklmn  psDijkl Bklmn rmn ¼ rij ;

r ¼ 1; 2; . . . ; N;

ð90Þ

s¼1

eijðrÞ þ

N  X



ðr;sÞ 1ðsÞ 1ðsÞ 0 Pijkl C klmn  psAijkl C klmn eðsÞ mn ¼ eij ;

r ¼ 1; 2; . . . ; N:

ð91Þ

s¼1 ðr;sÞ

ðr;sÞ

The objects Cijkl and Pijkl are defined in Eqs. (43) and (29), correspondingly, and s = h3/v0. Eqs. (90) and (91) can be rewritten in the canonical form similar to Eq. (30).

432

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

5. Effective elastic constants of periodic and random composites In this section, we consider an isotropic 3D-medium containing an infinite set of isotropic spherical inclusions. For such a medium, the tensors A and D in Eq. (8) take the forms

A¼

  1  j0 2 5  2j0 1 E þ E1  E2 ; 3 9l0 15l0

ð92Þ

D¼

  4l0 ð1  4j0 Þ 2 2l0 ð5 þ 4j0 Þ 1 1 2 E  E  E : 9 15 3

ð93Þ

5.1. Let the cell W0 contain only one spherical inclusion V0 In this case, Eqs. (1) and (4) have exact solutions that are constant inside the inclusion V0 and have the forms (Kanaun & Levin, 2008): 1 1 eij ¼  ijkl e0kl ;  1 ijkl ¼ Eijkl þ ð1  pÞAijmn C mnkl ;

ð94Þ

1 1 rij ¼ Kijkl r0kl ; K1 x 2 V 0: ijkl ¼ Eijkl  ð1  pÞDijmn Bmnkl ;

ð95Þ

Thus, for the averages hB1ri and hC1ei in Eqs. (74) and (86) we have

D E B1 r ¼ M ijkl r0kl ; ij

D E C 1 e ¼ Q ijkl e0kl ; ij

M ijkl ¼ B1ijmn Kmnkl ;

ð96Þ

Q ijkl ¼ C 1ijmn  mnkl :

ð97Þ

As the result, the tensors of the effective compliances B⁄ and effective moduli C⁄ of the composite in Eqs. (76) and (88) are

h i1 Bijkl ¼ B0ijkl þ pB1ijmn E1mnkl  ð1  pÞDmnpq B1pqkl ;

ð98Þ

h i1 C ijkl ¼ C 0ijkl þ pC 1ijmn E1mnkl þ ð1  pÞAmnpq C 1pqkl :

ð99Þ

It is easy to verify that (B⁄)1 = C⁄, and these equations lead to the formulas for the effective bulk and shear moduli K⁄,l⁄ of the composite in the forms

K ¼ K0 þ p



l0 l  l0 K 0 ðK  K 0 Þ

; ; l ¼ l þ p K þ ð1  pÞs1 ðK  K 0 Þ l0 þ ð1  pÞs2 l  l0

3K 0 s1 ¼ 3K 0 þ 4l0



6 K 0 þ 2l 0

: s2 ¼ 5 3K 0 þ 4l0

ð100Þ

ð101Þ

Here K0, l0 and K, l are the bulk and shear moduli of the matrix and inclusion materials. These equations coincide with the results of the well-known Mori–Tanaka method. 5.2. Simple cubic lattice of spherical inclusions If identical spherical inclusions compose a simple cubic lattice in the matrix material, the tensor of the effective elastic constants of the composite has cubic symmetry. A convenient basis for representation of such tensors consists of the following three linearly independent tensors

H1ijkl ¼ d1i d1j d1k d1l ;

H2ijkl ¼ E2ijkl  H1ijkl ;

 H3ijkl ¼ 2 E1ijkl  H1ijkl :

ð102Þ

In this basis, the tensors B0 and M in Eq. (76) for the effective compliance tensor B⁄ are presented in the form 0

0

0

B0 ¼ b1 H1 þ b2 H2 þ b3 H3 ; 0

M ¼ m1 H1 þ m2 H2 þ m3 H3 ; 0

ð103Þ

0

where the scalar coefficients b1 ; b2 ; b3 are: 0

b1 ¼

l0 þ 3K 0 2l0  3K 0 1 0 0 ; b2 ¼ ; b3 ¼ : 4l0 9K 0 l0 18K 0 l0

ð104Þ

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

433

As a result, the tensor B⁄ takes the form: 





B ¼ B0 þ pM ¼ b1 H1 þ b2 H2 þ b3 H3 ; 

0



0

b1 ¼ b1 þ pm1 ; b2 ¼ b2 þ pm2 ;

ð105Þ 0



b3 ¼ b3 þ pm3 :

ð106Þ

The coefficients m1, m2, m3 are to be calculated from Eq. (74)

D

B1 r

E

ij

¼ ðm1  2m3 Þr011 d1i d1j þ m2









r011 d2i d2j þ d3i d3j þ ðr022 þ r033 Þdij þ 2m3 r0ij :

ð107Þ

In order to find the value of m1, m2, m3, Eq. (67) for the stress field in the cell of the composite should be solved two times for the external stress fields r0ij in the forms

rij0ð1Þ ¼ d1i d1j ; r0ð2Þ ¼ d1i d2j jðijÞ : ij

ð108Þ

The products of these tensors with the tensor M are: 0ð1Þ

Mijkl rkl



¼ m1 d1i d1j þ m2 d2i d2j þ d3i d3j ;

0ð2Þ

M ijkl rkl 0(1)



¼ m3 d1i d2j þ d2i d1j : 0(2)

(1)

If the solution of Eq. (67) for the right-hand sides r and r are r through the components of the tensors hB1r (1)iand hB1r(2)ias follows

ð109Þ (2)

and r , the coefficients m1, m2, m3 are expressed

Fig. 7. The dependence of the relative shear modulus M⁄/l 0 of the composite with a simple cubic lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 27 inclusions; h is the step of the cubic node grid; + are the data for the composite with a cubic lattice of absolutely rigid spheres presented in Nunan and Keller (1984).

Fig. 8. The dependencies of the relative bulk modulus K⁄/K0 of the composite with a simple cubic lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 8 or 27 inclusions, h = 0.02; + are the data for the composite with a cubic lattice of absolutely rigid spheres presented in Nunan and Keller (1984).

434

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Fig. 9. The dependencies of the relative shear modulus l⁄/l0 of the composite with a simple cubic lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 8 or 27 inclusions, h = 0.02; + are the data for the composite with a cubic lattice of absolutely rigid spheres presented in Nunan and Keller (1984).

Fig. 10. The dependencies of the relative shear modulus M⁄/l0 of the composite with a simple cubic lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 8 or 27 inclusions, h = 0.02; + are the data for the composite with a cubic lattice of absolutely rigid spheres presented in Nunan and Keller (1984).

Fig. 11. The dependence of the relative bulk modulus K⁄/K0of the composite with a simple cubic lattice of soft spherical inclusions (E/E0 = 0.001, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 27 inclusions; h is the step of the cubic node grid;  are the data for the medium with a simple cubic lattice of spherical pores presented in Kushch (1997), + are the data presented in Iwakuma and Nemat-Nasser (1983).

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

m1 ¼ hB1 rð1Þ i11 ;

m2 ¼ hB1 rð1Þ i22 ;

m3 ¼ hB1 rð2Þ i12 : ⁄

⁄ 1

The tensor of effective elastic moduli C = (B ) 

C ¼

c1 H1

þ

c2 H2

þ

where scalar coefficients

c1 ¼

1   ðb þ b2 Þ; 4 1

ð110Þ

is presented in the form similar to (105)

c3 H3 ;

c1 ;

c2 ;

c2 ¼ 

435

ð111Þ c3  b2

4

are expressed via the coefficients

;

c3 ¼

 b1 ;

1      ; 4 ¼ ðb1  b2 Þðb1 þ 2b2 Þ: 4b3

 b2 ;

 b3

in Eq. (106)

ð112Þ

The tensor C⁄ may be also presented as follows

   1 C ¼ K  E2 þ 2l E1  H1 þ 2M H1  E2 ; 3 where the effective bulk modulus K⁄ and shear moduli l

1 K  ¼ ðc1 þ 2c2 Þ; 3

l ¼ c3 ;

1 M ¼ ðc1  c2 Þ: 2

ð113Þ ⁄

and M⁄ are

ð114Þ

Fig. 12. The dependence of the relative shear modulus l⁄/l0 of the composite with a simple cubic lattice of soft spherical inclusions (E/E0 = 0.001, m = m0 = 0.3) on their volume fractions p; the cuboid W (2  2  2) contains 27 inclusions; h is the step of the cubic node grid;  are the data for the medium with a simple cubic lattice of spherical pores presented in Kushch (1997), + are the data presented in Iwakuma and Nemat-Nasser (1983).

Fig. 13. The dependence of the relative shear modulus M⁄/l 0 of the composite with a simple cubic lattice of soft spherical inclusions (E/E0 = 0.001, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 27 inclusions; h is the step of the cubic node grid;  are the data for the medium with a simple cubic lattice of spherical pores presented in Kushch (1997), + are the data presented in Iwakuma and Nemat-Nasser (1983).

436

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Let the cube W of the sizes (2  2  2) contains 27 identical spherical inclusions which centers compose a simple cubic lattice. In the case of stiff inclusions (E/E0 = 1000, m = m0 = 0.3), the integral Eq. (67) was used for the calculation of the stress tensor r(x) inside the cube W. Here E and E0 are the Young moduli of the inclusion and the matrix materials. The results of the calculations of the relative effective bulk modulus K⁄/K0 and shear moduli l⁄/l0, M⁄/l0 are shown in Figs. 5–7 for the grid steps h = 0.1, 0.04, 0.02, and 0.01. The crosses in Figs. 5–7 are the values of the effective elastic constants of the composite with a simple cubic lattice of absolutely rigid spheres presented in Nunan and Keller (1984) (NK). The results of calculation of the effective elastic constants of the composites on the basis of a simple cubic cell that contains only 8 stiff inclusions are shown in Figs. 8–10. The effective elastic constants of the composite with a simple cubic lattice of soft inclusions (E/E0 = 0.001, m = m0 = 0.3) are shown in Figs. 11–13. Crosses in this figure are the values of the effective elastic moduli of the composites with a simple cubic lattice of spherical pores presented in Iwakuma and Nemat-Nasser (1983) (IN) and Kushch (1997) (Ku). In this case, the numerical scheme based on Eq. (91) for strains turns out to be more efficient than the scheme based on Eq. (90) for stresses. 5.3. Face-centered cubic lattices of spherical inclusions Let us consider the effective elastic moduli of the composite with a face-centered cubic lattice (FCC) of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3). In this case, the tensor of the effective elastic constants has cubic symmetry, and the dependencies of the relative effective bulk and shear moduli K⁄/K0, l ⁄/l0 and M⁄/l0 of the composite on the volume fraction p of the inclusions are presented in Figs. 14–16. Crosses in these figures are the values of the effective elastic constants of the FCC-composite with rigid spheres presented in Nunan and Keller (1984) (NK).

Fig. 14. The dependencies of the relative bulk modulus K⁄/K0 of the composite with the FCC-lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 14 inclusions; + are the data for the composite with the FCC-lattice of absolutely rigid spheres presented in Nunan and Keller (1984).

Fig. 15. The dependencies of the relative shear modulus l⁄/l0 of the composite with the FCC-lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 14 inclusions, h is the step of the cubic node grid; + are the data for the composite with the FCClattice of absolutely rigid spheres presented in Nunan and Keller (1984).

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

437

Fig. 16. The dependencies of the relative shear modulus M⁄/l0 of the composite with the FCC-lattice of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.3) on their volume fractions p; cuboid W (2  2  2) contains 14 inclusions, h is the step of the cubic node grid; + are the data for the composite with the FCClattice of absolutely rigid spheres presented in Nunan and Keller (1984).

5.4. Random set of spherical inclusions Let us consider a statistical model of a random set of spherical inclusions in the cubic region W of the  2  2).  sizes (2 ðiÞ ðiÞ ðiÞ First we take a regular lattice of identical spherical inclusions of the radii R0 centered at points Y ðiÞ ¼ y1 ; y2 ; y3 , where ðiÞ ðiÞ ðiÞ y1 ; y2 ; y3 are the Cartesian coordinates of the ith center. Then, the center of every ith inclusion is moved at the point ðiÞ ðiÞ ðiÞ y1 þ r1 ; y2 þ r2 ; y2 þ r2 , where r1, r2, r3 are independent random variables homogeneously distributed in the interval {(Rp  R0), (Rp  R0)}, and Rp is the radios of the inclusions corresponded to the percolation threshold of the initial regular lattice of inclusions embedded in the cube W. The inclusion overlapping is not allowed as well as their translations outside the region W. Spatial orientation of a cell with a fixed realization of the inclusions is also random, and the corresponding distribution over orientations is homogeneous. An example of the randomize FCC-cell is shown in Fig. 1 for the volume fraction of the inclusions p = 0.2. We consider firstly an example of stiff inclusions (E = 1000E0, m = 0.49, m0 = 0.25). In order to calculate the effective elastic constants of the composite, we take a fixed realization of a random field of inclusions in the cube W and solve Eq. (67) for stresses 6 times for the following external stress fields r0(k)

rij0ðkÞ ¼ dki dkj ; k ¼ 1; 2; 3; rij0ð4Þ ¼ d1i d2j jðijÞ ; rij0ð5Þ ¼ d1i d3j jðijÞ ; r0ð6Þ ¼ d2i d3j jðijÞ : ij

ð115Þ ðkÞ

We denote the corresponding solutions as rij ðxÞ. In order to calculate the tensor M in Eq. (76) for B⁄one has to multiply the ðkÞ solutions rij ðxÞ by the tensor B1ijkl ðxÞ and average the results over the region W. Finally, we obtain 6 constant tensors ðkÞ 1 ðkÞ T ij ¼ hB r iij (k = 1, . . . , 6) which components are connected with the components of the tensor M in Eq. (74) by the equations ðkÞ

M ijkk ¼ T ij ;

k ¼ 1; 2; 3; ð5Þ

M ij13 ¼ M ij31 ¼ T ij ;

ð4Þ

M ij12 ¼ Mij21 ¼ T ij ; ð6Þ

M ij23 ¼ M ij32 ¼ T ij :

ð116Þ

Now, we should average the tensor M over the spatial orientations of the cell with the fixed realization of the inclusions inside. As a result, for the average tensor hMi, we obtain the equation (see Appendix A)

  1 hMi ¼ m1 E2 þ m2 E1  E2 ; 3 m1 ¼ m2 ¼

1 ðM 1111 þ M 2222 þ M 3333 þ M 1122 þ M 2211 þ M 1133 þ M 3311 þ M 2233 þ M3322 Þ; 9

ð117Þ

ð118Þ

1 ½2ðM1111 þ M2222 þ M3333 Þ  ðM 1122 þ M 2211 þ M 1133 þ M 3311 þ M 2233 þ M 3322 Þ þ 3ðM 1212 þ M 1313 þ M 2323 Þ: 15 ð119Þ

In the next step, the values of the coefficients m1 ; m2 should be averaged over several independent realizations of the inclusions in W. The corresponding averages we denote as hm1 i and hm2 i. Finally, the tensor of effective elastic stiffness of the composite C⁄ is isotropic and takes the form

438

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Fig. 17. The dependencies of the relative bulk modulus K⁄/K0 of the composite with a random set of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.25) on their volume fractions p; cuboid W (2  2  2) contains 14 inclusions, h is the step of the cubic node grid; + are the data for the composite with a random set of absolutely rigid spheres presented in Segurado and Llorca (2002).

Fig. 18. The dependencies of the relative shear modulus l⁄/l0 of the composite with a random set of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.25) on their volume fractions p; cuboid W (2  2  2) contains 14 inclusions, h is the step of the cubic node grid; + are the data for the composite with a random set of absolutely rigid spheres presented in Segurado and Llorca (2002).

Fig. 19. The dependencies of the relative bulk modulus K⁄/K0 of the composite with a random set of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.25) on their volume fracture p, cuboid W (2  2  2) contains 1, 8, or 14 inclusions, h = 0.02, + are the data for the composite with a random set of absolutely rigid spheres presented in Segurado and Llorca (2002).

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

  1 C ¼ K  E2 þ 2l E1  E2 ; 3 K ¼

K0

; 1 þ 9pK 0 m1

l ¼

439

ð120Þ

l

0 : 1 þ 2pl0 m2

ð121Þ

The dependencies of the relative bulk K⁄/K0 and shear l⁄/l0 moduli of the composite with stiff spherical inclusions are presented in Figs. 17 and 18 for the randomized initial cubic cell with 27 spherical inclusions inside W (E/E0 = 1000, m = 0.49, m0 = 0.25). Crosses in these figures are the results of the calculations of the effective elastic constants based on the Finite Element Method for the RVE containing 30 identical absolutely rigid inclusions with periodical boundary conditions on the RVE sides (see Segurado & Llorca, 2002, (SL)). Figs. 19 and 20 present the dependencies of the effective elastic moduli of the composite with the random set of stiff inclusions on their volume fractions p. In these figures dashed lines correspond to the cell that contains 1 inclusion (Eq. (100)), lines with squares to the cell with 8 inclusions for the initial simple cubic cell, lines with circles to 14 inclusions for the initial FCC-cell. Crosses are the data presented in Segurado and Llorca (2002) (SL) for a random set of absolutely rigid spheres. The values of the bulk and shear moduli of the composite with a random set of soft inclusions (E = 0.001E0, m = m0 = 0.25) are compared with the numerical solution presented in Segurado and Llorca (2002) for a medium with spherical pores in Figs. 21 and 22. In this case, Eq. (81) of the EFM for strains was used, and it is seen that calculations based on a cell with

Fig. 20. The dependencies of the relative shear modulus l⁄/l0 of the composite with a random set of stiff spherical inclusions (E/E0 = 1000, m = m0 = 0.25) on their volume fractions p; cuboid W contains 1, 8, or 14 inclusions, h = 0.02, + are the data for the composite with a random set of absolutely rigid spheres presented in Segurado and Llorca (2002).

Fig. 21. The dependencies of the relative bulk modulus K⁄/K0 of the composite with a random set of soft spherical inclusions (E/E0 = 0.001, m = m0 = 0.25) on their volume fractions p; cuboid W (2  2  2) contains 1, 8, or 14 inclusions, h = 0.02, + are the data for the porous medium presented in Segurado and Llorca (2002).

440

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Fig. 22. The dependencies of the relative shear modulus l⁄/l0 of the composite with a random set of soft spherical inclusions (E/E0 = 0.001, m = m0 = 0.25) on their volume fractions p; cuboid W (2  2  2) contains 1, 8, or 14 inclusions, h = 0.02, + are the data for the porous medium presented in Segurado and Llorca (2002).

8 inclusions predict the values of the effective elastic constants of the composite that practically coincide with the results presented in Segurado and Llorca (2002). For the construction of the tensor Q in Eqs. (86) and (88), Eq. (81) for the strain tensor should be solved 6 times with the right hand side similar to Eq. (115) for every fixed realization of a random set of inclusions inside W. Then, the tensor Q is to be averaged over the orientations similar to tensor M (see Appendix A). In both cases of stiff and soft inclusions, the deviations of the effective elastic constants in different realizations were small and do not exceed 3%.

6. Conclusions An efficient numerical method for the solution of 3D-problems of elasticity for a medium with a finite number of isolated inclusions is proposed. The main difficulty in the numerical solution of 3D-integral equations of elasticity for a cell of a composite material is a large number of approximation functions (approximating nodes) that should be taken in order to achieve acceptable accuracy. In the present method, this difficulty is overcome by the following means. First, the usage of Gaussian approximating functions allow calculating the elements of the matrix of the discretized problems in closed analytical forms and thus, to calculate them fast. Second, the use of regular node grids provides Toeplitz properties to the matrix of the discretized problem. As a result, the number of independent elements of this matrix is essentially reduced, and in addition, the FFT algorithms may be applied to calculating matrix–vector products with such matrices. It accelerate essentially the process of iterative solution of the discretized problems for a cell. In the framework of the method, many parts of the problem may be calculated in parallel. For instance, the elements of the block-matrices in Eqs. (33)–(39) and products of these matrices with vectors may be carried out independently, and at the same time. This is another advantage of the method and a source of accelerating of the computation process. Combination of the proposed numerical technique and a (self-consistent) effective field method allows one to solve the homogenization problem for matrix composite materials with sufficient accuracy. It is shown that this solution may be obtained on the basis of a rather small cell of the composite that contains about a dozen of inclusions. In this work, only a selfconsistent effective field method is considered. The developed numerical technique may be also used in the framework of another self-consistent method, e.g., the effective medium method and its various versions (see, e.g., Kanaun & Levin, 2008 for details). In this case, a typical cell of the composite should be embedded in a homogeneous effective medium, and the tensor of the effective elastic stiffness of the latter should be chosen such as the average stress and strain fields over the cell are connected by the same fourth-rank tensor.

Appendix A Let us introduce a Cartesian basis of three mutually orthogonal vectors e1, e2, e3. In this basis, the tensor M in Eq. (74) is presented in the form

M¼

3 X i;j;k¼1

Mijkl ei  ej  ek  el :

ðA:1Þ

441

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

The rotated basis qi is defined by the equation

qi ¼

3 X

Q ij ðu; h; cÞej ;

i ¼ 1; 2; 3;

ðA:2Þ

j¼1

0 6 u 6 2p;

0 6 h 6 p;

0 6 c 6 2p;

ij

where Q (u, h, c) is an orthogonal matrix, and (u, h, c) are Euler angles that define orientation of the new basis qi. If a fixed configuration of inclusions is rotated by the Euler angles (u, h, c), the tensor M0 that corresponds to the new configuration takes the form

M0 ¼

3 X

M ijkl qi  qj  qk  ql ;

ðA:3Þ

i;j;k¼1

where Mijkl are the components of the tensor M in the original basis ei. The average value of the tensor M over the orientations (u, h, c) calculated as follows

hMi ¼

3 X

 M ijkl qi  qj  qk  ql ;

ðA:4Þ

i;j;k¼1

 1 qi  qj  qk  ql ¼ 8p 2

Z

qi  qj  qk  ql dl;

dl ¼ sinðhÞdudhdc:

ðA:5Þ

where the integration is spread over all the values of the Euler angles. It is assumed that the distribution of a fixed configuration of the inclusions over the orientations is homogeneous. By averaging the tensor M over the orientations the following equations should be taken into account:

 1  2 qi  qi  qi  qi ¼ E þ 2E1 ; i ¼ 1; 2; 3; 15

i  1  2 2E  E1 ; i – j; i; j ¼ 1; 2; 3; q  qi  qj  qj ¼ 15

i  1  1 j j i q q q q ¼ 3E  E2 ; i – j; i; j ¼ 1; 2; 3; 15

i  q  qj  qi  qj ¼ 0; i – j; i; j ¼ 1; 2; 3;

i   q  qi  qj  qk ¼ qi  qj  qk  qi

j  ¼ q  qk  qi  qi ¼ 0; i – j; i – k; j – k; i; j; k ¼ 1; 2; 3;

i  i  i i j q  q  q  q ¼ q  qi  qj  qi

i 

 ¼ q  qj  qi  qi ¼ qj  qi  qi  qi ¼ 0; i – j; i; j ¼ 1; 2; 3:

ðA:6Þ

As a result, the averaged tensor hMiturns out to be isotropic and takes the form

  1 hMi ¼ m1 E2 þ m2 E1  E2 ; 3

ðA:7Þ

m1 ¼

1 ðM 1111 þ M 2222 þ M 3333 þ M 1122 þ M 2211 þ M 1133 þ M 3311 þ M 2233 þ M3322 Þ; 9

m2 ¼

1 ½2ðM1111 þ M2222 þ M3333 Þ  ðM 1122 þ M 2211 þ M 1133 þ M 3311 þ M 2233 þ M 3322 Þ þ 3ðM 1212 þ M 1313 þ M 2323 Þ: 15

References Davis, I., Hatch, R., Yener, M., & Chompooming, K. (1993). Stress and strain fields in random pack of rigid bounded spheres. Physical Review B, 47(5), 2530–2545. Golub, G., & Van Loan, C. (1996). Matrix computations (3d ed.). The Johns Hopkins University Press. Gusev, A. (1999). Representative volume element size for elastic composites: A numerical study. Journal of the Mechanics and Physics of Solids, 45, 1449–1459. Iwakuma, T., & Nemat-Nasser, S. (1983). Composites with periodic microstructure. Computers and Structures, 16, 13–19. Kanaun, S. (2009). Fast calculation of elastic fields in a homogeneous medium with isolated heterogeneous inclusions. International Journal of Multiscale Computational Engineering, 7/4, 263–276. Kanaun, S., & Levin, V. (2008). Self-consistent methods for composites static problems vol. 1. Solids Mechanics and its Applications (vol. 148). Springer. Kunin, I. A. (1983). Elastic media with microstructure II. Springer Series in Solid-State Sciences (vol. 44). Berlin, Heidelberg, NY, Tokyo: Springer. Kushch, V. (1997). Microstresses and effective elastic moduli of a solid reinforced by periodically distributed spherical particles. International Journal of Solids and Structures, 34, 1336–1353. Maz’ya, V., & Schmidt, G. (2007). Approximate approximation. American mathematical society. Mathematical Surveys and Monographs, 141.

442

S. Kanaun, E. Pervago / International Journal of Engineering Science 49 (2011) 420–442

Mikhlin, S. (1965). Multidimensional singular integrals and integral equations. Oxford, NY: Pergamon Press. Moulinec, H., & Suquet, P. (1994). A fast numerical method for computing the linear and nonlinear properties of composites. Comptes Rendus de l’ Academie des Sciences Paris II, 318, 1417–1423. Moulinec, H., & Suquet, P. (1998). A numerical method for computing the overall response of nonlinear composites with complex microstructure. Computational Methods in Applied Mechanics and Engineering, 157, 69–94. Nunan, K., & Keller, J. (1984). Effective elasticity tensor of periodic composites. Journal of the Mechanics and Physics of Solids, 45, 1421–1448. Press, W., Flannery, B., Teukolsky, S., & Vetterling, W. (1992). Numerical recipes in FORTRAN: The art of scientific computing (2nd ed.). Cambridge University Press. Segurado, J., & Llorca, J. (2002). A numerical approximation to the elastic properties of sphere-reinforced composites. Journal of the Mechanics and Physics of Solids, 50, 2107–2121. Zohdi, T., & Wriggers, P. (2005). An introduction to computational micromechanics. Berlin, Heidelberg: Springer.