Micromechanically based stochastic finite elements: length scales and anisotropy

Probabilistic Engineering Mechanics 11 (1996) 205-214 PII:

ELSEVIER

SO266-S920(96)o0015-X

Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0266-8920/96 $15.00

Micromechanically based stochastic finite elements: length scales and anisotropy K. AizeMeh Department of Materials Science and Mechanics, Michigan State University, East Lansing, MI 48824-1226, USA

&

M. Ostoja-Starzewski Institute of Paper Science and Technology, and Georgia Institute of Technology, 500 lOth Street, Atlanta, GA 30318-5794, USA The present stochastic finite element (SFE) study amplifies a recently developed micromechanically based approach in which two estimates (upper and lower) of the finite element stiffness matrix and of the global response need first to be calculated. These two estimates correspond, respectively, to the principles of stationary potential and complementary energy on which the SFE is based. Both estimates of the stiffness matrix are anisotropic and tend to converge towards one another only in the infinite scale limit; this points to the fact that an approximating meso-scale continuum random field is neither unique nor isotropic. The SFE methodology based on this approach is implemented in a Monte Carlo sense for a conductivity (equivalently, out-of-plane elasticity) problem of a matrix-inclusion composite under mixed boundary conditions. Two versions are developed: in one an exact calculation of all the elements' stiffness matrices from the microstructure over the entire finite element mesh is carried out, while in the second one a second-order statistical characterization of the mesoscale continuum random field is used to generate these matrices. Copyright © 1996 Elsevier Science Ltd.

1 INTRODUCTION

this challenge rather than as a criticism of the conventional SFE methods. It has been shown in our recent works4-7that the development of SFE accounting for spatial random material fluctuations hinges on two concepts: determination of scale-dependent stochastic constitutive laws from microstructure, and setting up of a variational formulation of the problem consistent with such laws. It follows from the first concept that the stochastic finite element calculations must involve three length scales: micro (level of heterogeneity); meso (finite element level); and macro (global structure). Our procedure is thus divided into two steps: a micro-meso step and a meso-macro step. The first step involves a calculation of the stiffness matrices at the meso-scale for a specific type of material microstructure; the second step involves a solution of a global SFE problem. Here, every finite element is 'cut out' of a heterogeneous medium and characterized by two different stiffness matrices, according to whether essential or natural

Conventional studies in the area of stochastic finite element methods (SFEM or SFE) rely on somewhat hypothetical assumptions concerning the variability in the material properties of the system. In particular, the linear elastic SFE analyses o f structures (e.g. Refs 1-3) assume a locally isotropic Hooke's law: O'ij : A(X, ~)(~ij ~kk -P 2#(x, w)eij

(1)

where the Lam6 constants A, # are taken to be random fields, x E B body domain and w E ~ the sample space of all realizations. Various assumptions concerning the constitutive coefficients are made in these studies such as taking Young's modulus as a random field o f Gaussian type with Poisson's ratio as a constant. This state of affairs is reflected on page 99 o f Ref. 3 as follows: "the issue o f random field meshes remains largely unexplored and much work has to be done on this particular aspect of S F E M . " The present paper is written as a response to 205

206

K. Alzebdeh, M. Ostoja-Starzewski

boundary conditions are employed. Correspondingly, two different variational principles need to be employed in order to bound the global structure's response, namely the potential energy principle and the complementary energy principle. To better display the connection to micromechanics, we present a quantitative comparison of the upper and lower estimates of the effective stiffnesses with the classical Voigt, Reuss, as well as the upper and lower Hashin bounds. This is done in the context of a simply formulated, yet very relevant material problem: a matrix-inclusion composite with round disks that are stiffer than the matrix. It is shown that with the increasing size of a finite element, the ensemble averages of the upper and lower stiffnesses of the composite tend to converge to the lower Hashin bound, which is the correct macroscopic solution. The proposed micromechanics based SFE method is finally illustrated on a classical boundary value problem for which the homogeneous medium solution is given in Ref. 8. We randomize this problem by assuming the microstructure of a matrix-inclusion composite. Next we develop two methodologies, both of them in a Monte Carlo sense, in that they involve a solution of a (finite) number of heterogeneous media problems. The first methodology - - termed an exact method - - depends on an exact micro-meso calculation of the stiffness matrix over each finite element of the global domain and for each realization B(w). The second method - - termed a covariance-based method - - relies on a statistical second-order specification of the continuum random field approximating a given material microstructure on the same scale as that given by the finite elements. This is then used to generate the stiffness matrix for each B(w) over the entire finite element mesh.

2 MICROMECHANICALLY BASED VARIATIONAL F O R M U L A T I O N 2.1 Principle of minimum potential energy Development of a stochastic finite element (SFE) methodology for structural mechanics applications accounting for the randomness in material constitutive response defines the goal of this communication. The basic tenet of this methodology is based on the following observation: randomness in material response is the result of its discrete nature. As expounded below, such a starting point allows a formulation of random continuum fields approximating the effective response at some larger scales as well as an SFE analysis related to a particular material. Let B(w) denote a heterogenous body, whose characteristic microscale, d, is much smaller than the body's macroscopic dimension L, that is, d << L. A paradigm is given schematically in Fig. 1 wherein w is an

Fig. 1. A macroscopic body B(w) with a boundary S = St + Su, the finite element mesh (mesoscale), and the microstructure (microscale d); a typical finite element with an inclusionmatrix microstructure is shown. index from the probability (sample) space fL A random medium B is a set of all B(w), all of them occupying the same spatial domain of volume V. The microstructure is chosen here to be of a (circular) inclusion-matrix type, and is shown in a subregion of the body; both phases are linear elastic. However, it should be pointed out that any other elastic microstructure with a single characteristic microscale - - such as a polycrystal or concrete - can be treated in the same manner as outlined in this paper. Let us now consider a boundary value problem for B(w). From the principle of minimum potential energy for each B(w) of B we have 6II(w) =

lvO'(W).6F.dx - Ivf.6udx f t-6udS=O J st

(2)

where the three terms stand for the (random) strain energy, the work done by the body forces, and the surface tractions' work, respectively. Also, S stands for the boundary of B(w), while St and Su for those parts of it where natural and essential boundary conditions, respectively, are being applied. For a typical B(w), we introduce a triangulation of its body domain (Fig. 1), a closed set, into a finite number of (closed) triangles V i (i = 1 , 2 , . . . , N ) such that N

v=Uv

i

(3)

i=1

where eqn (2) is replaced by

N

Nj I v,

e dV-

v, f. endV

N

Here iCe(w) is the stiffness tensor of the ith particular

Micromechanically based stochastic finite elements

207

finite element which should be determined from the micromechanics. The superscript 'e' indicates that an essential boundary condition is employed in its calculation as will be explained below. It is well-known that the principle [eqn (4)] will result in a set of linear algebraic random equations: [K(w)]{ U} = {F}

(5)

where [K(w)] is the global stiffness matrix, {U} is the displacement vector and {F} is the force vector. Equation (5) will then provide the basis for a Monte Carlo solution of the SFE problem, as well as for any approximate solution in terms of the moments (such as the second-order solution). However, the correctness of this solution will fundamentally depend on the input in terms of the constitutive response, i.e. the specification of iCe(w) for all the elements, and on the choice of interpolation functions. Let us observe that the fact that we work with a minimum potential energy formulation is consistent with an imposition of displacement-controlled boundary conditions, u = ~. x

(6)

on the boundary iS of each element, where ~ is a volume-averaged strain. This is one of two classical problem statements of micromechanics which are used to specify the effective response (e.g. Refs 9 and 10); the second one is based on the stress-controlled boundary conditions and will be used later. It is important to note that the interpolation functions over the finite element, as dictated by eqn (6), are linear:

u(xl, x2) = a + bXl + C.X2

(7)

In fact, taking either higher-order interpolation functions or other than triangular finite elements would result in effective stiffnesses dependent on the mean strain ~. This dependence would only vanish in the case of finite elements being smaller than the size of either phase, disks or matrix between the disks, and falling entirely into either phase. Problem (6) leads to an effective tensor C~(w), whereby we use the subscript '6' as a dimensionless parameter specifying the ratio 6 = Lid between the meso-scale (finite element) and the micro-scale (inclusion size) lengths. Clearly we have C~(w)= ice(w), which shows a link between notations used in micromechanics and finite elements, respectively. Determination of C~(w) for any given finite element in a closed form is practically impossible and a resort has to be made to any of a number of numerical or hybrid analytical-numerical methods. One may, for example, employ a fine finite element mesh to discretize the domain of a given finite dement, such as the one shown in Fig. 2. Our own preference is for a very fine finite difference mesh (also called spring-network), that is, a mesh whose spacing is several times smaller than the size

.d-H++¢++ Fig. 2. A finite element discretized by a fine finite difference (spring-network) mesh. of a heterogeneity, see e.g. Ref. 11. Several different spring-network schemes may be set up; they are typically based on periodic lattices and very suited for dasticity-type problems, but the particular choice depends on the case at hand. Recently they have become increasingly popular in an even harder field: fracture in heterogeneous media, see e.g. Refs 12, 13, and 20 and the references therein. Our spring-network is also illustrated in Fig. 2 where a square mesh is used. Thus, for instance, in the case of an out-of-plane elasticity problem, the displacement field u is governed locally by

u(i,j)[kr + k l + k u +kd] - u(i+ 1, j)kr - u ( i - 1,j)k 1 - u(i,j+ 1)ku - u ( i , j - 1)kd --- 0

(8) whereby kr, kl, ku and kd are the spring constants in all four directions (right, left, up, and down) from an (i, j) node. These constants are calculated from a series spring model, e.g. we have k r = [1/C(i,j) + 1/C(i+ l,j)] -1, where C(i,j) and C(i+ l, j) are the local elastic stiffnesses at the respective nodes. The advantage of the spring-network method lies in a very rapid assignment of all the spring constants, which obviates the needs for remeshing that would be required in a finite element method. Once the displacement field has been found, C~(w) is determined from the effective medium postulate: V Uspring-networ k = Ucontinuum ~ ~-~" C~(w). ~ (9) The important thing to note here is that this effective random tensor is, in general, anisotropic. This is obvious from the following observation: disks would need to be arranged in very special ways for the effective response to be isotropic. It follows, by rather conventional arguments based on Lebesgue measures of line patterns of plane figures, 14 that C~(w) is almost surely anisotropic. The global stiffness matrix [K(w)] in eqn (5) can now be synthesized from the stiffness matrices of all the elements: N [K(w)] = Z [ i K ( w ) ] (10) i=1

K. Alzebdeh, M. Ostoja-Starzewski

208

whereby each [iK(w)] is obtained as follows: [iK(w)] =

Jiv[io]T[iC(w)][iG] d V

= [iG]T[iC(w)][iG]iV

(11)

In the above [ic(w)] is the stiffness tensor of a homogenized ('smeared-out') finite element, given by ice(w); [iG] is the gradient matrix. 2.2 Principle of minimum complementary energy

An alternative finite element approach is based on a force formulation, which is obtained by starting from the principle of minimum complementary energy for each B(w) of B:

postulate of presence of a very large number of microscale constituents (molecules in a fluid, crystals in a metal, disks in a composite, grains in concrete, etc.) in a given volume. Equations (15) and (16) also demonstrate the basis for two approximating random continuum fields that will be discussed in Section 2.3. The variational principle [eqn (13)] results in a set of linear algebraic random equations:

[f(w)]{F}

= {V}

(17)

where [f(w)] is the global flexibility matrix, which is to be assembled from the iSn(w) tensors. As eqn (15) indicates, the response based on the natural boundary conditions is softer than that based on the essential boundary conditions. In practical terms, this means that the global response of a given specimen B(w) can be bounded by solutions of eqns (5) and (17).

6II'(w) = I v 6 o ' . ~ ( w ) d V - J v 6 f . u d V 2.3 Ensemble average formulation - [

J su

6t.udS= 0

(12)

With the discretization [eqn (3)] we obtain m*(w) =

dV i=1

Vi

dV _

Vi

N

where iSn (w) should be an effective compliance tensor of every ith particular finite element. The consistent way to determine iSn(w) in this force formulation is to impose a natural boundary condition (hence the superscript 'n'): t = a.n

(14)

on each finite element's boundary and find S~(w); O is the volume-averaged stress. Again, 6 indicates the inherent scale dependence and we see that S ~ ( w ) isn(w). Also, by the same arguments as in Section 2.1, S~(w) is almost surely anisotropic. The very important thing to note here is that the response tensor is not unique, but rather [S~(w)]-t < Cg(w)

(15)

This discrepancy vanishes in the limit 6 ~ ~ which defines the classical representative volume element (RVE) of deterministic continuum mechanics; see Refs 15 and 16 for extensive discussions. The inequality [eqn (15)] is understood in the following sense: for two fourth-rank tensors A and B, an order relation B < A means tijBijkltkl <_ tijAijkltkl for every symmetric tensor t. Furthermore, it may be shown that a following scaledependent ordering of bounds holds:

<

<__c (w) <__cb(w) v6' < (16)

In physical terms, the existence of the RVE is the

The foregoing developments provide a stepping stone for an SFE analysis of the response of an ensemble B = {B(w); w e f~}. First, we observe that the window (or finite elemen0 may be placed anywhere in the body domain, and thus two estimates of constitutive response related by eqn (15) may be derived. Consequently, passing to an ensemble, we arrive at a concept of two random continuum approximations: B~ = {B$(w); w E [~} and B~ = {B~(w), w E fl}; see also Refs 4 and 7. That is, we have two random continuum fields:

c (w, x),

s (w, x)

(18)

both of which have almost surely anisotropic realizations. Given the statistical characterization of these fields one may solve the global problem w by w. Furthermore, providing the random medium has the property of statistical homogeneity (invariance of probability distributions of microstructural properties with respect to arbitrary translations) and ergodieity, one can set up a hierarchy of scale (i.e. 6) dependent bounds 15'16 on the effective tensor C eH of an infinitely large window ((.) stands for the ensemble average): (sg)

<

_<

< (C~,) < C v

__ c

<_

VS' < 6

(19)

where SR-= (S~) is a so-called Reuss bound, and C v -= (C~) is a so-called Voigt bound, e.g. see Ref. 9. The main point is as follows: the larger the window scale 8 the closer are both bounds, and in the limit 6 ~ oo they coincide, i.e. the difference betw~n both types of boundary conditions disappears. One has thus two options: either (i) use both variational principles, eqns (4) and (13), for each specimen B(w) and repeat it in a Monte Carlo sense to get an assessment of the entire ensemble; or (ii) derive variational principles valid for the average solution.

209

Micromechanically based stochastic finite elements Regarding the latter option, two ensemble average formulations may be obtained from the above stochastic ones. 6 These are useful when, for example, one intends to solve the global problem (i.e. the meso-macro step) by an FE program based on the displacement method only. In that case one has to do the following: (i) Employ the stiffness matrix [K(w)] based on the ice(w) tensors and their distributions in order to calculate an upper bound on response. (ii) Carry out the averaging on S~(w) in order to find C~ = (S~) -I so that the inhomogeneity in surface traction distribution on any given finite element's boundary as well as the randomness in the compliance tensor are removed. The resulting stiffness matrix, based on (S~) -1, may now be employed to set up a deterministic stiffness matrix, which allows one to find the lower bound on the global response via the displacement method. This is, in fact, the method we shall employ in the sequel. Dually, an analogous approach based on the flexibility method could be employed. Here one would do the following: (i) Employ the flexibility matrix If(w)] based on the iSe(0~) tensors and their distributions in order to calculate a lower bound on response. (ii) Carry out the averaging on C~ (w) in order to find S~ = (C~) -1 so that the inhomogeneity in displacement distribution on any given finite element's boundary as well as the randomness in the stiffness tensor are removed. The resulting compliance matrix, based on (C~) -l, may now be employed to set up a deterministic flexibility matrix, which allows one to find the lower bound on the global response via the displacement method. Additionally, one has to keep in mind the pros and cons of the displacement and the flexibility methods known from the deterministic finite element analyses, e.g. Refs 8 and 17.

3 M E T H O D O L O G Y OF SOLUTION 3.1 Exact procedure This procedure is termed exact because the stiffness matrix of each finite element is determined for a specific heterogeneous microstructure B(w) within its domain; in practice, this is done by using a very fine finite difference (or element) mesh in order to capture the randomness at the micro level. The method can be summarized as follows: (i) Generation of a (finite) number of realizations B(w) of the random matrix-inclusion composite under study. A particular geometric model of randomly placed disk inclusions following a planar Poisson process with sequential inhibition is employed to represent the actual microstructure, see e.g. Refs 6 and 7.

(ii) Finite element discretization of the domain of the body; this specifies the mesoscale 6. (iii) Micromechanics calculation: considering one realization, for each finite element, a very fine finite difference mesh at micro level (Fig. 2) is adapted to calculate the stiffness tensors C~(oJ) and C~(w). We note here that these two tensors depend on the 6 parameter. (iv) Finite element solution: feeding these two stiffness tensors, separately, into a finite element code, two bounds (upper and lower) on the actual response of B(w) will be obtained (recall eqns. (5) and (17)). (v) By repeating this process for a sufficient number of realizations, any required moment of the bound can be calculated. 3.2 Covariance-based method The covariance-based method is an approximate method that decreases the CPU time of calculations very dramatically in that it takes advantage of the mesocontinuum random field model of a given microstructure. The method can be summarized as follows: (i) First, there must be developed two continuum random field approximations of the medium on a given meso-scale 6, which is defined by the chosen finite element mesh. This description includes, at least, the first and second-order moments of C~(oJ) and C~(,~). As an example presented in Ref. 7 indicates, for diskinclusion composites at moderate volume fractions, correlations of stiffnesses (and compliances) between contiguous elements are practically zero. A connection may be established here with the classical, deterministic micromechanics, by comparing the averages of our 6-dependent bounds, described by the hierarchy [eqn (9)], with the classical Voigt C v, Reuss C R, as well as the upper and lower Hashin bounds C~, C~. These are given, for a two-phase composite, by the following well-known equations (see e.g. Ref. 18): c v = c R

+ c(2)j 2)

=

CuH =

= c(,) Jr

-1

1/[c(1) - c(2)] + , % / 2 c ( 2 )

J(2)

- c(,)]

(20)

Herein we assume both phases (k = 1,2) to be locally isotropic - - Cij(k)= 6ij C(k) - - where C(k) stands for the property (transverse conductivity or anti-plane shear modulus) of either phase; their respective volume fractions are denoted by f(k). (ii) Knowledge of all the averages plus auto(and cross-) covariances between all the constitutive coefficients allows a generation (e.g. Ref. 19) of the

K. Alzebdeh, M. Ostoja-Starzewski

210 X2

T=0

0000000

~... •

>:.': •

•oo_•.

dT/dx=0

O0 01

.:.

-:.--

00 • • • O

•

•

T=0

dT/dy = 0

x1

Fig. 3. The domain, boundary conditions, and finite element discretization of the domain with 6 = 20; a 400 x 400 pixel mesh of the composite is employed.

stiffness tensor of an ith finite element via a linear transformation:

.fLit 0

C22(03)

= /L21

C12(~)

[L31 +

L22

0

L32 L33 (C22)

Z2(03) Z3(o)) (21)

(c 2) In eqn (21) [L] is a lower triangular matrix of the corresponding Cholesky-Banachiewicz decomposition. (iii) An approximating assumption with respect to the exact method is made here, namely {Z} in eqn (21) is a vector of independent random numbers taken to follow a standard normal distribution N(0, 1). However, it is noted here that according to micromechanical analyses7 the stiffnesses and compliances at the meso-scale 6 are not Gaussian, so that a higher-order, rather than just second-order, information on approximating continuum random fields is needed to obtain truly good results. (iv) With the global stiffness matrix specified, the same tasks as (iv) and (v) in Section 3.1 are carried out here.

4 NUMERICAL RESULTS The stochastic finite element methodology described above is used here to solve a conductivity problem of

a square plate (global structure) of 400Al × 400Al dimensions (AI is a ,"tit length) under mixed boundary conditions and a , ..orm unit heat source distribution. Note that, by virtue of mathematical analogies, this is equivalent to solving several other problems: out-ofplane elasticity, membrane displacements, etc. This problem has been posed and solved, in the case of a homogeneous medium, as in Example 7-14 in Ref. 8 on pages 421-426, and so it serves as a reference case here. The problem statement is as follows: -Au = 1

(0,y) =

(0,x) = 0

(x,y) ~ (0,1) × (0,1) u(1,x)=u(x, 1)=O (22)

In the sequel, the microstructure is taken as a matrixinclusion composite with inclusions' diameters being 5Al, see Fig. 3. Both matrix and inclusions are locally homogeneous isotropic continua, where C(]) = 0.5 for k = 1 (inclusion), and C(l) = 5.0 for k = 2 (matrix); this corresponds to a contrast C(2)/C(]) = 10. By taking the contrast as unity we recover the same numbers as those reported in Ref. 8. All the inclusions are modeled as circular disks embedded in the matrix according to a planar Poisson process with sequential inhibition; that is, their center points are placed sequentially, with a minimal distance between any two adjacent disks centers being 140% of the disk's diameter. Volume fraction equals 20%.

Micromechanically based stochastic finite elements

211

2.0 8 = 4 (upper bound) 8 = 4 (lower bound)

O

8 = 10 (upper bound) 1.5-

[]

8 = 10 (lower bound)

--4k--

8 = 20 (upper bound) 8 = 20 (lower bound)

/"/

................ Hashin-upper

.......

~.) 1.0-

Hashin-lower

V

Voigt

,~

Reuss

.j-tl / 1i ,f

0.5

0.0

'

0.0

I

I

I

0.1

0.2

0.3

'

I

0.4

'

0.5

f'2 Fig. 4. Calculation of scale-dependent upper and lower bounds, vs volume fraction f(2) of inclusions, on the effective stiffness tensor Celf of a triangular element for 6 = 4, 10, 20; the Voigt, Reuss and Hashin bounds are also shown.

In accordance with Section 3, the exact method relies on the 'cutting' of finite elements out of the microstructure, and, next, on calculating their effective stiffnesses Cg(w) and C~(a;) = [S~(t.,d)] - 1 using the fine finite difference mesh discussed in Section 2.1. On the other hand, the covariance-based method relies on the calculation of the first and second moments of these tensors, which provides the second-order characterization of the random fields C~(~, x) and S~(w,x). This information is then used, through the matrix [L], as input in eqn (21) to generate the global stiffness matrix. Thus Fig. 4 depicts the scale-dependent (i.e. 6) upper and lower bounds on the effective stiffness tensor of a triangular element for 6 = 4, 10,20 for a range of volume fractions of inclusions. The classical Voigt, Reuss, as well as the upper and lower Hashin bounds calculated according to eqn (20) are also shown, and in accordance with the hierarchy [eqn (19)], our upper and lower bounds get ever closer onto the C H bound as 6 increases. O f course, one would need a window of infinite extent to actually obtain this homogenization limit. Let us note here that the minimum spacing of inclusion's centers - - 1.4 o f their diameter - - permits

volume fractions no greater than just over 30%, and this is where we terminate our graph. Our global (i.e. meso-macro) solution to this problem is based on the displacement formulation outlined in Section 2.3. Thus it allows the determination of the distribution of the upper bound and the mean of the lower bound on the overall system response. As a measure of this response we choose the temperature at the lower left comer of the plate. Thus two sets of bounds on temperature, obtained using exact and covariance-based methods, respectively, are plotted in Fig. 5. In each case these results are based on 50 realizations B(w) of B. As expected, with the mesh size increasing, each set of the bounds converges. Two aspects merit discussion here: the limit t5 ~ e~ and the difference between both sets of bounds. The infinite 6 limit is, strictly speaking, unattainable: the finite element has to remain finite. We can only hope to approach it in the case when the material randomness (as measured by the contrast in the microstructure) becomes very small compared to the size t5 of the finite elements. Indeed this is the situation encountered in numerous finite element applications, whereby the error

K. Alzebdeh, M. Ostoja-Starzewski

212 0.4,

0~11) e" .....

~

0" ............

--e-- ....

-0 ......................

tl-

. . . . . . . . . . . . . .,

0.2-

I!.1-

--l---

Exact, eslmntlal Exact, natural

0.0 0

i 5

•- - B - - .

Covartance, essential

•- O - - -

Covaaance, natural

I 10

t 15

20

8 Fig. 5. Bounds on response (temperature of comer point) as a function of 6 for both methods.

due to the neglect of randomness in Hooke's law is acceptable. This also provides an explanation of a rather low interest in SFE problems in the applied mechanics community. Regarding the difference between both sets of bounds we first note that the exact method yields a somewhat lower result compared to that of the covariance-based method. This discrepancy is due to two simplifications: (i) assumption of zero correlation between the contiguous elements; (ii) neglect of the higher-order statistical information in the covariance-based method and the assumption of Gaussianity of the stiffnesses. Thus, to make a final assessment of this method, one should remove the first simplification and generate the global stiffness matrix in accordance with all the spatial correlations. Fortunately, the difference between both methods is not that significant, and all the merits of this approximate second-order method still remain to be explored. Let us note once again that all the information, per each finite element per each w, is being retained in the exact method. Furthermore, since we used a displacement FE formulation, the randomness in response is incorporated in the upper bound only, i.e. the standard deviation, and is plotted in Fig. 6. One can observe its convergence to zero with 6 going to infinity, although this limit, just as the limit in Fig. 5, is unattainable: the finite element has to remain finite.

Finally, the CPU time is an important criterion in determining the practicality of any numerical technique. With that objective in mind, we measured and compared the CPU time for both the micromechanics (MM step) and the finite element (FE step) calculations in the exact method, see Fig. 7. As 6 increases, the FE computational time decreases very rapidly due to a decreasing number of finite elements, while that of the MM increases dramatically due to the complexity (number of inclusions) of boundary value problems that have to be solved for ever larger windows. Thus an optimal 6 = 3 could be deduced here. We end by noting that, on the other hand, the CPU time of the covariance-based method is very small since the generation of the stiffness matrix according to eqn (21) is rapid - - this replaces the MM cost, and only the FE cost applies.

5 CONCLUSIONS (i) A micromechanics analysis of the effective constitutive response indicates that there is no unique locally isotropic random field of elastic properties; rather, two random continuum fields C~(w,x) and S~(w,x) may be introduced corresponding, respectively, to the displacement-controlled or tractioncontrolled boundary value problems at the mesoscale.

Micromechanically based stochastic finite elements

213

co

b "r'-

0 .

o

B

t~ D

c~ "O tt~ 0 0

I

I

1

5

10

15

20

8 Fig. 6. Standard deviation of the upper bound (for the reference temperature) as a function of 6.

jo oo ~ po

5-

f I.......

FE t MM

na

J

s 4

3-

J

J

ip JI

os Sp

.

sS os oS

os S wo ~

1-

s~ S°

0

I

I

I

I

I

1

I

;

i

2

4

6

8

10

12

14

16

18

20

6 Fig. 7. The CPU computation time of the micromechanics (MM) and the finite element (FE) steps plotted as a function of 6.

214

K. Alzebdeh, M. Ostoja-Starzewski

(ii) The displacement formulation of the finite element scheme for a given B(~o) realization of the random medium should be based on the C~(~o, x) field, with the force formulation based on the S~(~o, x) field. (iii) For an ensemble average solution, averages of these two random fields may be used as input to a finite element program in either a displacement or a force formulation. (iv) In either case of point (iii) above, there are two options: an exact method retaining all the information per finite element for each B(0J), or a covariancebased method retaining the second-order information only. (v) These methods have been illustrated on a specific boundary value problem.

ACKNOWLEDGMENT We benefitted from the suggestions of two anonymous referees. Support by the National Science Foundation under grant no. MSS 9202772 and by the Research Excellence Fund from the State of Michigan is gratefully acknowledged.

REFERENCES 1. Brenner, C. E., Stochastic finite elements (literature review). Internal Working Report no. 35-91, Institute of Engineering Mechanics, University of Innsbruck, Innsbruck, 1991. 2. Ghanem, R. G. and Spanos, P. D., Stochastic Finite Elements: A Spectral Approach. Springer, Berlin, 1991. 3. Kleiber, M. and Hien, T. D., The Stochastic Finite Element Method. Wiley, New York, 1993. 4. Ostoja-Starzewski, M., Micromechanics as a basis of random elastic continuum approximations. Probabilistic Engineering Mechanics, 1993, 8, 107-114. 5. Alzebdeh, K. and Ostoja-Starzewski, M., Micromechanically based stochastic finite elements. Finite Element Analysis and Design, 1993, 15, 35-41.

6. Ostoja-Starzewski, M., Micromechanics as a basis of stochastic finite elements and differences - - an overview. Applied Mechanics Reviews, 1993, 46, S136-S147. 7. Ostoja-Starzewski, M., Micromechanics as a basis of continuum random fields. Applied Mechanics Reviews, 1994, 47, 221-230. 8. Reddy, J. N., Applied Functional Analysis and Variational Methods in Engineering. McGraw-Hill, New York, 1986. 9. Willis, J. R., Variational and related methods for the overall properties of composites, Advances in Applied Mechanics, 1981, 21, 1-78. 10. Torquato, S., Random heterogeneous media: microstructure and improved bounds on effective properties. Applied Mechanics Reviews, 1991, 44, 37-76. 11. Day, A. R., Snyder, K. A., Garboczi, E. J. and Thorpe, M. F., The elastic moduli of a sheet containing circular holes. Journal of the Mechanics and Physics Solids, 1992, 40, 1031-1051. 12. Ostoja-Starzewski, M., Sheng, P. Y. and Jasiuk, I., Influence of random geometry on effective properties and damage formation in composite materials. ASME Journal of Engineering Material Technology, 1994, 116, 384-391. 13. Grah, M., Alzebdeh, K., Sheng, P. Y., Vaudin, M. D., Bowman, K. J. and Ostoja-Starzewski, M., Brittle intergranular failure in two-dimensional microstructures: experiments and computer simulations. Acta Materialia, 1996 (in press). 14. Stoyan, D., Kendall, W. S. and Mecke, J., Stochastic Geometry and its Applications. Wiley, New York, 1987. 15. Huet, C., Application of variational concepts to size effects in elastic heterogeneous bodies. Journal of the Mechanics and Physics Solids, 1990, 38, 813-841. 16. Sab, K., On the homogenization and the simulation of random materials. European Journal of Mechanics, A/Solids, 1992, 11, 585-607. 17. Ne~as, J. and Hlavfi6ek, I., Mathematical Theory of Elastic and Elasto-Plastic Bodies: An Introduction. Elsevier, New York, 1981. 18. Hashin, Z., Analysis of composite materials; a survey. Journal of Applied Mechanics, 1983, 50, 481-505. 19. Elishakoff, I., Probabilistic Methods in the Theory of Structure. Wiley, New York, 1983. 20. Ostoja-Starzewski, M., Sheng, P. Y. and Alzebdeh, K., Spring network models in elasticity and fracture of composites and polycrystals. Computational Materials Science, 1996 (in press).