- Email: [email protected]
Scale effects in materials with random distributions of needles and cracks
Mechanics of Materials 31 (1999) 883±893
www.elsevier.com/locate/mechmat
Scale eects in materials with random distributions of needles and cracks Martin Ostoja-Starzewski * Institute of Paper Science and Technology, Georgia Institute of Technology, 500 10th St., N.W., Atlanta, GA 30318-5794, USA Received 14 September 1998; received in revised form 6 May 1999
Abstract According to a classical prescription of micromechanics (R. Hill, J. Mech. Phys. Solids 11, 1963), a representative volume element (RVE) is well de®ned when the response under uniform displacement (Dirichlet) boundary condition becomes the same as that under uniform stress (Neumann) boundary condition. We study the convergence of both responses in anti-plane elasticity of sheets with non-periodic, random distributions of thin needle-shaped inclusions. By lowering the stiness of inclusions and increasing their aspect ratio (up to 100), we approach the situation of cracks embedded in a matrix. We show that, with the needles' stiness decreasing and their slenderness growing, the RVE tends to be very large. The statistics of the ®rst and second invariants of both response tensors are very well modeled by a Beta probability distribution. For moderate aspect ratio needles, the coecient of variation of the second invariant is found to stay at about 0.5 irrespective of the window size, the mismatch in stiness between the inclusions and the matrix, and the needle aspect ratio. Ó 1999 Elsevier Science Ltd. All rights reserved. Keywords: Scale eects; Mesoscale; Micromechanics; Random composites; Cracked media
1. Introduction Dating back to the seminal work of (Budiansky and O'Connell, 1976), the problem of eective moduli of materials with randomly arranged ellipsoidal inclusions or cracks has received considerable attention over the past two decades; for extensive reviews see Kachanov (1993), NematNasser and Hori (1993), Taya and Arsenault (1989). The works on this subject can be classi®ed into two categories: eective medium theories and bounds. Those in the ®rst category include self-
*
Tel.: +1-404-894-5700; fax: +1-404-894-4778. E-mail address: [email protected] (M. Ostoja-Starzewski)
consistent models, Mori±Tanaka models, mean®eld models, and computer simulations. The works in the second category are based on the Hashin±Shtrikman bounds. All these studies strive to establish a macroscopically eective stiness tensor on scales corresponding to a representative volume element (RVE) of an equivalent continuum. However, no information has been obtained on the size of this RVE as a function of the density of the inclusions or cracks, the mismatch (contrast) between the inclusion material and the matrix, and various types of geometric disorders. This paper attempts to answer some of these issues. The study is conducted in the setting of antiplane classical elasticity, that leads directly to twodimensional ®eld problems of two-phase media where one phase is the matrix, while the other is
0167-6636/99/$ - see front matter Ó 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 1 6 7 - 6 6 3 6 ( 9 9 ) 0 0 0 3 9 - 3
884
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
the inclusion. Both phases are locally isotropic with C m and C i being the stinesses of the matrix and the inclusion, respectively. The inclusions are taken as thin needles, generated as germ-grain Boolean processes (Stoyan et al., 1987). We focus on the soft inclusion case (C i =C m < 1), and also try to assess the limit of zero stiness C i =C m ! 0. We investigate the RVE problem in the spirit of Hill's de®nition (Hill, 1963), which says that the RVE is well de®ned when the response under uniform displacement (Dirichlet) boundary condition becomes the same as that under uniform stress (Neumann) boundary condition. While the perfect coincidence of both responses for non-periodic, random microstructures occurs on scales in®nitely larger than the size of an inhomogeneity, their convergence may be anywhere from rapid through slow depending on the parameters. This issue had been studied earlier for 2D and 3D elasticity of several composite materials; see e.g. Amieur et al., 1995; Huet, 1990, 1991, 1997; Ostoja-Starzewski and Wang, 1989; OstojaStarzewski, 1994, 1998a; Ostoja-Starzewski and Schulte, 1996, and references therein.
2. Needle systems 2.1. Random geometry In our two-dimensional problem, needles are taken as narrow d w rectangles, where d and w are the length and width, respectively. In order to generate a ®eld of needles, such as the one shown
in Fig. 1(a), we proceed in the following steps. We employ a Poisson point process to generate the needles' center points. Strictly speaking, the Poisson point process of intensity l is speci®ed by (|A| is a Lebesgue measure of the set A) P fN A kg
kjAj ÿk e jAj; k!
A R2 :
1
It occurs on the entire R2 plane, and so, its numerical simulation must necessarily be done on a compact subset, which actually involves a binomial point process. To specify a needle's orientation we generate an angle h from 0; p according to a probability density f h 1=p, where h is the angle of inclination to the x-axis. This results in a statistically isotropic con®guration of needles, and hence, in a macroscopically isotropic (in the ensemble average sense) response. The subset of R2 on which the binomial process is conducted must be larger than the actual L L window domain which will be tested for apparent moduli ± lines that originated from Poisson points outside the window, located up to a distance d=2, should be accounted for. With L and d at hand, we are ready to introduce a non-dimensional parameter d
L d
2
specifying the scale of observation relative to the size of a heterogeneity, i.e. a mesoscale. Another basic parameter specifying the composite is the aspect ratio of needles d=w. By varying the aspect ratio from 1 up through higher values we can
Fig. 1. (a) A realization of random microstructure of a two-phase, needle-matrix composite material, and the window-scale concept; d is the typical needle length and L the window size. (b) A spring network as a basis for resolution of needles and ellipses.
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
model systems having rectangle-type through needle-type inclusions. A realization of an isotropic ®eld of one hundred needles of length d 100 units and width w 1 unit in a 1000 1000 window is shown in Fig. 2; k 10ÿ4 . Already at this rather low density, random clustering of needles is easily seen. In this paper we focus on two cases: long (d 100) and short (d 10) needles, both at w 1. Another important parameter is the volume fraction of inclusions de®ned by v i V i =Vtotal , where the volume V is really just the area in our planar system. The range through which v i may be varied will depend on the chosen model of placement of inclusions in the matrix. For example, in the case of no-overlap condition of inclusions we would be looking at a random packing problem, while in the case of no such condition ± as implied by the Poisson point process of needle centers ± we can go theoretically up to 100%. However, the latter situation is of no interest per se when studying eective elastic properties. Rather, of interest is the dependence of eective and apparent moduli on the mismatch of individual phases' moduli and on the window
885
scale d, possibly up to the point of percolation of needles. 2.2. Basic mechanics We focus on matrix-inclusion, locally isotropic, elastic materials with inclusions aligned with the x3 , subjected to anti-plane shear. Inclusions are needle shaped. That is, locally (for either phase) the Hooke's law of the matrix (m) and the inclusion (i) material is given in the x1 ; x2 -plane of loading by ri Cij ej ; Cij C
m
i; j 1; 2;
dij
or
C i dij ;
3
where, for simplicity of notation, we denote ri ri3 ;
ei ei3 ;
i; j 1; 2:
4
Parallel to the subscript notation we shall occasionally use the symbolic notation, so that, for example: Cij C. The governing equation of the piecewise-constant material is 2 o u o2 u 0; C C i or C m ; C ox21 ox22 5 u u3 : The isotropy of both stiness tensors C m and in (3) leads to a so-called contrast C a C i =C m , also called a mismatch in other studies. It is clear that by increasing the contrast we can go to very sti inclusions, and approximate a rigid case. Similarly, by decreasing the contrast, we can go to very soft inclusions and approach a system with holes, and, for high d=w, a system with cracks. C i =C m < 1 is, in fact, the case of interest in the present study. Throughout the paper we shall adopt C m 1. Once a speci®c geometry (e.g. Fig. 2) of needles is generated and the assignment of local phase properties to the matrix and needles is made, we are dealing with one realization B x of a random medium B. The latter is taken as a set fB x; x 2 Xg; each B x is treated as a deterministic problem. Fig. 3 depicts the parameter plane ± window scale d L=d and the contrast a C i =C m ± that is investigated in this paper. See Section 4 for further details. i
Fig. 2. A ®eld of long (1 100) needles in a 1000 1000 square, generated from a Poisson point process of density k 10ÿ4 , i.e. there are 100 needles.
886
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
Fig. 3. A parameter plane: window scale d L=d and the contrast a C i =C m . Bold points show the particular cases that were investigated in support of the conclusions in this paper. For long needles (d 100) this translates to d L=100, while for short needles (d 10), to d L=10. Numbers next to the bold points indicate the sample sizes analyzed.
2.3. Numerical solution method In order to solve the ®eld equation (5) of our two-phase composite, with all the interactions, we employ a spring network method, Fig. 1(b). The idea is to approximate the planar, piecewise-constant continuum by a very ®ne square mesh of springs, e.g. (Ostoja-Starzewski and Schulte, 1996). We employ a square mesh for discretization of the displacement ®eld u, so that the governing equations become u i; jkr kl ku kd ÿ u i 1; jkr ÿ u i ÿ 1; jkl ÿ u i; j 1ku
6
ÿ u i; j ÿ 1kd 0: Here i and j are the coordinates of mesh points, and kr ; kl ; ku and kd are de®ned from the series spring model ÿ1
kr 1=C i; j 1=C i 1; j ; ÿ1
kl 1=C i; j 1=C i ÿ 1; j ; ku 1=C i; j 1=C i; j 1ÿ1 ;
7
ÿ1
kd 1=C i; j 1=C i; j ÿ 1 ; where C i; j is the stiness ± i.e., C i or C m ± at the mesh point i; j. The spring network of Fig. 1(b) can be used to model needles of width equal to (or greater than) the mesh spacing, or other inclusions such as ellipses. Just like every numerical method, the spring
network has its pros and cons. Here they are our main considerations. The generation of each and every con®guration of needles is extremely rapid and it does not require any remeshing as the same underlying network 1000 1000, say) is always being used. In particular, no remeshing is needed to handle intersecting needles, as these just happen to appear as certain pixel patterns on the network, a sample of which is shown in Fig. 2. This is an important consideration when dealing with hundreds of realizations (con®gurations) from the sample space X. Further, to better resolve the local stress/strain ®elds within and around the inclusions one may decrease the mesh spacing, just as one would proceed with a ®nite element method (FEM). However, the FEM is limited by a higher price of costly and cumbersome remeshing for each and every new con®guration required in statistical studies. Additionally, we are not concerned with local stress/strain ®eld details as they are known to have little eect on the resulting, eective/apparent moduli. Finally, both methods ± spring networks and ®nite elements ± suer from a need to work with a very large number of nodes to cover the free spaces between the needles. This wasting of memory, and the associated increase in size of the resulting algebraic Ax b problems, can be avoided by employing a boundary element method, whereby the nodes are placed at the needles' and external window boundaries only. We have just such a computer program under development, and we have been able to establish, in speci®c cases, that the dierences between the spring network and the boundary element methods are within 5% for contrasts as low as 10ÿ4 . This accuracy of spring networks is perfectly sucient to rapidly establish the elastic moduli of a number of dierent B x realizations from the random medium B, and their corresponding statistics. 3. Scale-dependent hierarchies of bounds on apparent moduli As stated in the introduction, we assess the apparent moduli for any ®nite window domain of
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
887
mesoscale d depending on whether a uniform strain e0j or a uniform stress r0j is prescribed. In the ®rst case, we should choose essential (or, displacement, Dirichlet) boundary conditions while in the second case, we should choose natural (or, traction, Neumann) boundary conditions. The ®rst setup is ri Cije e0j
under
u x e0j Lj
8x 2 oB;
~
~
8
where rij is the resulting mean (volume average) stress, which leads to a mesoscale stiness C ed (e stands for essential); oB is the boundary of B. The second setup is ei Sijn r0j
under t x r0j nj x 2 oB; ~
~
9
where eij is the resulting mean (volume average) strain, and leads to a mesoscale compliance S nd (n stands for natural). Determination of either second-rank tensor, C ed or S nd , requires three tests. Henceforth, we follow Huet (1990, 1991, 1997), by adopting the terminology ``apparent moduli'' for any mesoscale situation (d ®nite), and reserve ``eective moduli'' for the macroscopic limit (d ! 1). For any realization B x, a window's response on the mesoscale (d ®nite) is, under these de®niÿ1 tions, non-unique ± because Cije 6 Sijn almost surely ± and anisotropic. At this stage, we should note the mixing (and hence, ergodic) as well as strict-sense stationarity property of the Poisson point process underlying our two-phase microstructure (Stoyan et al., 1987). It can then be shown from the variational principles that the ensemble averages of these two tensors provide, with the increasing scale d, an ever tighter pair of bounds on the macroscopic apparent stiness tensor C eff , C R S R
ÿ1
hS n1 i
ÿ1
6 hS nd0 i
ÿ1
6 hS nd i
6 C eff 6 hC ed i 6 hC ed0 i 6 hC e1 i C V
ÿ1
1990; Sab, 1992; Ostoja-Starzewski and Schulte, 1996) to correspond to commensurate partitions of a body domain in the sense that d0 is a multiple of d. However, all the numerical results presented here and those given in our previous studies (as well as those of Huet and his coworkers for 3D problems) speak in favor of (10) being true for arbitrary d0 < d. This is now shown in Appendix A for any microstructures described by wide-sense stationary random ®elds. Hierarchy (10) is depicted in Figs. 4 and 7 for two random needle systems; see Section 4 below for further details. At this point we note that the scale (i.e., d) and contrast (i.e., a) dependencies of both tensors in (10) follow, with very good accuracy, the laws ®rst found for 2-D Bernoulli lattices (Ostoja-Starzewski and Schulte, 1996), namely htr C ed i a0 a1 exp a2 dÿa3 a ; ÿ ÿ1 htr S nd i b0 b1 exp b2 dÿb3 a :
11
8d0 < d: 10
eff
Fig. 4. Eect of decreasing contrast a C i =C m on the divergence of bounds hC e i and hS n iÿ1 , at the window size d 5, for the long needle system one realization of which is shown in Fig. 2; C m 1. The eective stinesses C SC 0:61 and C holes 0:65 computed by Eqs. (B.1) and (B.2) of the selfconsistent and the mean ®eld methods are also shown.
Clearly, C is understood as C 1 C d jd1 . The order relation employed in (10) means that t B t 6 t A t for any vector t 6 0 and two second rank tensors A and B. Furthermore, d and d0 in (10) were shown in the previous works (Huet,
4. Scale-dependent hierarchies of bounds: numerical results We investigate two types of systems: (i) those made of long needles (1 100), and (ii) those made of short needles (1 10). The probability distri-
888
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
butions of mesoscale moduli that are involved in the hierarchy (10) are assessed in terms of the statistics of their two invariants ± half-traces and radii of the corresponding Mohr's circles (C C ed or S nd ) tr C=2
C11 C22 ; 2q
R C12;max
2 C11 ÿ C22 2 =4 C12 :
12
We ®rst focus on the long needles. Fig. 2 displays a typical ®eld of such needles, with a volume fraction 0.01, which is a result of placing 100 needles of size 1 100 in a 1000 1000 window. Fig. 4 displays the eect of decreasing contrast ÿ1 C i =C m on the divergence of bounds hC e ihS n i and at the window size d 5. For the sake of comparison we also indicate the eective stinesses C SC and C holes computed by Eqs. (B.1) and (B.2) of the self-consistent and eective medium models reported in Appendix B. Figs 5 and 6 display, for three scales d 1; 4; and 5, the statistics of tr C=2 and R for tensors C ed and S nd at contrast a 10ÿ4 , respectively. They have been obtained by solving, through a spring
network method of Section 2.3, Dirichlet and Neumann boundary value problems for a number of composites B x of the set B, all being generated in a Monte Carlo sense. Also shown in Figs. 5 and 6 are the Beta function probability density ®ts to the traces and radii; note that the distributions of tr C=2 for C ed are quite dierent from those of S nd . A number of other distributions were also tried ± in particular, Chi, Chi-Square, Exponential, Frechet, Gumbel Min, Log-Normal, Rayleigh and Weibull. However, none of them oered the same overall goodness-of-®t as Beta. It is indeed of interest that Beta is also very satisfactory for other types of composites, both in the elastic (Ostoja-Starzewski, 1998a,b) as well as inelastic (Alzebdeh et al., 1998) response ranges (see Fig. 7). Note that while the traces have asymmetric distributions, their character is very similar for both tensors. The Mohr's circles' radii have even stronger skewness, and, most interestingly, their coecients of variation (CV), for a given type of boundary conditions (either Dirichlet or Neumann), are practically constant (!) with the changing mesoscale d.
e Fig. 5. Probability densities of Ciie =2 (left column) and C12;max (right column) of the needle-matrix composite (at a 10ÿ4 ) whose one realization is shown in Fig. 2 at d 1; 4; and 5.
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
889
n Fig. 6. Probability densities of Siin =2 (left column) and S12;max (right column) of the needle-matrix composite (at a 10ÿ4 ) whose one realization is shown in Fig. 2 at d 1; 4 and 5.
Fig. 7. Eect of increasing widow scale on the convergence of the hierarchy of scale dependent bounds for the composite whose one realization is shown in Fig. 2; C m 1; C i 10ÿ4 . The eective stiness C holes 0:65 computed by Eq. (B.2) of the mean ®eld method is also shown.
e CV C12;max 1:3;
n CV S12;max 0:55:
13
Turning to the short needles, in Fig. 8 we present a typical realization B x. Next, Fig. 9 displays the eect of increasing contrast C i =C m ÿ1 on the divergence of bounds hC e i and hS n i , at the window size d 50. The Poisson process density has been chosen so as to result in roughly the same
Fig. 8. A ®eld of short (1 10) needles in a 1000 1000 square, generated from a Poisson point process of density k 10ÿ4 , i.e. there are 1000 needles.
volume fraction of the soft (i.e., needle) phase as was the case with long needles; here it equals 0.015. The eective stiness C holes 0:43 computed by Eq. (B.2) of the mean ®eld method is also
890
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
C 3 should decrease with delta due to a decrease of C 1 with delta. In eect, (15) allows a determination of any statistic of C 3 from those of C 1 and C 2 . We do not pursue this here. Figs. 4 and 9 lead to two more conclusions. Firstly, the Dirichlet bounds hC ed i are seen to asymptote very quickly with increasing contrast, ÿ1 while the Neumann bounds hS nd i keep falling o.
Secondly, the behavior of both bounds is indicative of the goodness of approximation of a twophase system, with needles having a low albeit ®nite property, relative to a material with holes: we see that the higher the aspect ratio, the closer is the system to the material with holes. In the studies of eective moduli of heterogeneous materials, the resulting C eff is presented graphically versus the volume fraction of one of the phases. For our system of short needles, this is shown in terms of the C holes of Eq. (B.2) of Appendix B against x nL2eff almost all the way to the percolation point at 5.9 (see Fig. 10). In addition, we also plot here the Dirichlet and Nuemann moduli hC ed i and hS nd iÿ1 at two mesoscales: d 10 and 50. Clearly, the three data points of Fig. 9 at C i =C m 10ÿ5 are at x 12:1 (since n 0:01; a 5; b 0:5), while all the data at (two, three and four times) higher volume fractions ± x 2:42; 3:63; and 4.84 ± are new. Thus, this ®gure displays (i) a very slow approach of hC ed i and hS nd iÿ1 to the RVE (i.e., curve C holes ), and (ii) a discrepancy between C holes and the Dirichlet, as well as C holes and the Neumann bounds, respectively. Note that C holes corresponds to an eective (macroscopic) response of a very large random system under periodic boundary conditions (Garboczi et al., 1991).
Fig. 9. Eect of increasing contrast a C i =C m on the divergence of bounds hC e i and hS n iÿ1 , at the window size d 50, for the short needle system one realization of which is shown in Fig. 8; C m 1. The eective stiness C holes 0:43 computed by Eq. (B.2) of the mean ®eld method is also shown.
Fig. 10. Normalized overall moduli hC ed i and hS nd iÿ1 , at d 10 and 50, and the eective stiness C holes from Eq. (B.2), for a random ®eld of short 1 10 needles (such as that of Fig. 8), as functions of the volume fraction x nL2eff . Data were computed only at discrete intervals x 1:21; 2:42; 3:63 and 4.84.
shown. The CVs of Mohr's circles' radii are now not only constant for a given type of boundary condition (either Dirichlet or Neumann) and independent of the changing mesoscale d, but also much closer to each other, e 0:45; CV C12;max
n CV S12;max 0:55:
14
Most interestingly, these numbers are very close to those reported recently for two-phase planar Voronoi mosaics (Ostoja-Starzewski, 1998b), or those found for matrix-inclusion composites, where inclusions are circular. Thus, we conjecture e n is close to CV S12;max for matethat CV C12;max rials with heterogeneities of moderate rather than extreme aspect ratios. The third invariant Cij Cjk Cki also appears to have similar universal properties. In fact, observing that, in 2D problems i; j 1; 2, h i 2 =2 15 C 3 Cij Cjk Cki C 1 3C 2 ÿ C 1
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
5. Closure The problem of scale dependence and of approach to the RVE in heterogeneous materials is a dicult one ± even in the simple setting of twophase composites in anti-plane shear ± and we ®nd the computational mechanics to be the most effective method of attack. In this paper we have investigated scale dependence of apparent moduli in systems having up to one hundred needles of aspect ratio 100, or up to 1000 needles of aspect ratio 10. While arbitrary geometries, from dilute right down to the percolation point, could be handled, with all the interactions being taken into account, we have focused on the volume fraction of about 1% to 1.5%. We believe this study is useful because: (1) The results provide guidance as to the size of the RVE, up to whatever desired discrepancy between the Dirichlet and Neumann type bounds (C ed and S nd ) one might want to choose. This, in turn, provides the key criterion in, say, the choice of ®nite element sizes in solving large (macroscopic) boundary value problems. (2) In the case when the above mentioned discrepancy is substantial and when one is interested in a macroscopic boundary value problem, one may have to choose to work with a stochastic, rather than a deterministic, ®nite element method. The window is then recognized as a mesoscale ®nite element, whose properties are statistical. The Dirichlet and Neumann type stinesses enter, respectively, the stiness and ¯exibility matrices, which themselves are derived from the principles of minimum potential and complementary energy principles, which then jointly result in bounds on macroscopic response, e.g., Ostoja-Starzewski and Wang (1999) and Ostoja-Starzewski (1999). (3) Probability distributions of C ed or S nd have been found to follow, in terms of their ®rst and second invariants, the Beta distribution to a high accuracy. The coecient of variation of the second invariant is seen to behave like a universal material constant: for short needles, it is constant irrespective of the contrast and the mesoscale, and equals 0.45 for the Dirichlet and 0.55 for the Neumann bound. The closeness of these values to those found for other composites (circular disk
891
inclusions in a matrix, or two-phase Voronoi mosaics), also for a wide range of contrasts and mesoscales, is striking. Acknowledgements Extensive comments of Prof. C. Huet (Ecole Polytechnique Federale de Lausanne) on the research reported here are gratefully acknowledged. We also bene®ted from the discussions with Prof. I. Jasiuk (Georgia Tech) and Mr. R. Muzzolini (Ecole des Mines de Paris). This research was made possible by the grant CMS-9713764 from the National Science Foundation. Appendix A. Validity of scale dependent bounds for non-commensurate partitions We want to prove the inequalities (10) for an arbitrary pair of (meso)scales d0 < d, and not just for commensurate ones, i.e., for partitions in which d nd0 ; n being a natural number. It will suce to focus on the Dirichlet bounds, because then the Neumann bounds follow by an analogous argument. We know that the hierarchy holds rigorously for the commensurate sequences of scales. Let us therefore consider two separate cases of the hierarchy (10) for commensurate partitions: one at an arbitrary d1 and another at d2 , whereby d1 < d2 < 2d1 . Thus, in the ®rst case, we have a sequence of inequalities 6 hC e4d1 i 6 hC e2d1 i 6 hC ed1 i 6
A:1
while in the second case, we have another sequence 6 hC e4d2 i 6 hC e2d2 i 6 hC ed2 i 6
A:2
Let us proceed by contradiction: we assume that the inequalities (A.1) and (A.2) are not consistent with each other in the following way, > hC e4d2 i < hC e4d1 i > hC e2d2 i < hC e2d1 i > hC ed2 i < hC ed1 i >
A:3
This, however, would imply that the microstructure is characterized by a Dirichlet bound which,
892
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
while displaying an overall decreasing behavior with the scale d, also ¯uctuates with a period increasing like nd. One can certainly construct random, ergodic microstructures with such scale dependent ¯uctuation. This would, however, contradict the assumption of statistical homogeneity of the material properties that is always taken in the de®nition of an RVE. From the standpoint of stochastic processes, there are two types of the statistical homogeneity: strict-sense stationarity, and wide-sense stationarity. In the ®rst case we are assured of the invariance of any probability distribution of material properties with respect to arbitrary shifts, while in the second case, of the invariance of the mean only with respect to such shifts along with the dependence of two-point correlation functions on the interpoint separations only. The classical RVE de®nition does not say exactly which of these is involved. Nevertheless, considering that the hierarchy (10) is stated in terms of the averages, it suces to choose the setting of material properties speci®ed via wide-sense stationary random ®elds. The ¯uctuation of Dirichlet bound in (A.3) would not even be consistent with the wide-sense stationarity ± as a result, strict-sense stationarity would also be violated. Summarizing, we conclude that the hierarchy (10) holds for non-commensurate mesoscale sequences. Other than (A.3) types of inconsistencies between (A.1) and (A.2) may also be considered, but then one is led to similar contradictory conclusions as above. Appendix B. Eective medium models and bounds As mentioned in the introduction, there is a host of eective medium models as well as formulas for rigorous bounds. Of these the case which is applicable to our system is that of elongated holes in a matrix ± apparently, no results are available for elongated, soft inclusions. First, Nemat-Nasser and Hori (1993) give a self-consistent estimate CSC for a random dilute distribution of frictionless cracks of density f C SC p 1ÿf : C m 2
B:1
As con®rmed by Fig. 4, this formula works well for our long 1 100 needles of negligible stiness, even at contrast 10ÿ3 , but, understandably, cannot be applied to the short needles. However, there exists in the physics literature a result valid for systems with needles of any aspect ratio, and having arbitrarily strong interactions (due to the presence of needle±needle intersections). According to Garboczi et al. (1991), the eective conductivity of a plane with random holes (equivalent, by a mathematical analogy, to antiplane shear) is given by a formula C holes x x x2 1ÿ 1 ÿ C m 5:90 5:90 24:97 1:3 x ÿ1 1 ; x nL2eff ; B:2 3:31 where C1 is the conductivity of the matrix, n is the number of holes per unit area, and Leff is a parameter dependent on the hole geometry. In the case of ellipses, Leff 2 a b, where a and b are the semimajor and semiminor axes. This is a universal formula in the sense that it also holds for other than elliptical holes, and right down to the percolation threshold xc 5:9 0:4. Clearly, it grasps the strong interaction eects very well, its only drawback being that it does not work for holes ®lled with a second phase. Therefore, we employ it for our systems at very low contrast, and, in fact, use it as an indication whether a system at low contrast (soft needles), say 10ÿ4 , is equivalent to the one of zero contrast (voids, cracks). For long needles (Fig. 4) the answer is in the armative, while for the short ones (Fig. 9), one needs to go to contrast 10ÿ5 or lower. Finally, we can consider variational bounds of one kind or another. Since we are dealing here with an isotropic system where the connected phase is stier than the inclusions, we could employ the upper bound for planar two-phase materials (Hashin, 1983). Unfortunately, a quick inspection reveals that this bound is so wide as to be of no practical use compared to the numbers we get via computational mechanics or from the effective medium formula (B.2). Therefore, we do not plot it.
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
References Alzebdeh, K., Al Ostaz, A., Jasiuk, I., Ostoja-Starzewski, M., 1998. Fracture of random matrix-inclusion composites: scale eects and statistics. Int. J. Solids Struct. 35, 2537. Amieur, M., Hazanov, S., Huet, C., 1995. Numerical and experimental assessment of the size and boundary conditions eects for the overall properties of granular composite bodies smaller than the representative volume. In: Parker, D.F., England, A.H. (Eds.), IUTAM Symposium on Anisotropy, Inhomogeneity and Nonlinearity in Solid Mechanics. Kluwer Academic Publisher, Dordrecht, pp. 149±154. Budiansky, B., O'Connell, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct. 12, 81±97. Garboczi, E.J., Thorpe, M.F., DeVries, M.S., Day, A.R., 1991. Universal conductance curve for a plane containing random holes. Phys. Rev. A 43, 6473. Hashin, Z., 1983. Analysis of composite materials ± a survey. J. Appl. Mech. 50, 481. Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. J. Mech. Phys. Solids 11, 357. Huet, C., 1990. Application of variational concepts to size eects in elastic heterogeneous bodies. J. Mech. Phys. Solids 38 (6), 813. Huet, C., 1991. Hierarchies and bounds for size eects in heterogeneous bodies. In: Maugin, G.A. (Ed.), Continuum Models and Discrete Systems, vol. 2, Longman Scienti®c & Technical, pp. 127±134. Huet, C., 1997. An integrated micromechanics and statistical continuum thermodynamics approach for studying the fracture behaviour of microcracked heterogeneous materials with delayed response. Eng. Fract. Mech. 58 (5±6), 459.
893
Kachanov, M., 1993. Elastic solids with many cracks and related problems. Advance Applied Mechanics, vol. 30. Academic Press, New York, p. 259. Nemat-Nasser, S., Hori, M., 1993. Micromechanics: Overall Properties of Heterogeneous Solids. North-Holland, Amsterdam. Ostoja-Starzewski, M., Wang, C., 1989. Linear elasticity of planar Delaunay networks: random ®eld characterization of eective moduli. Acta Mech. 80, 61. Ostoja-Starzewski, M., 1994. Micromechanics as a basis of continuum random ®elds. Appl. Mech. Rev. 47 (1, Pt. 2), 221. Ostojas-Starzewski, M., Schulte, J., 1996. Bounding of eective thermal conductivities of multiscale materials by essential and natural boundary conditions. Phys. Rev. B 54 (1), 278. Ostoja-Starzewski, M., 1998a. Random ®eld models of heterogeneous materials. Int. J. Solids Struct. 35, 2429. Ostoja-Starzewski, M., 1998b. Random ®eld models and scaling laws of heterogeneous media. Arch. Mech. 50 (3), 549. Ostoja-Starzewski, M., 1999. Microstructural Disorder, Mesoscale Finite Elements, and Macroscopic Response. Proc. Roy. Soc. Lond. A 455, 3189±3199. Ostoja-Starzewski, M., Wang, X., 1999. Stochastic ®nite elements as a bridge between random material microstructure and global response. Comp. Meth. Appl. Mech. Eng. 168 (1±4), 35. Sab, K., 1992. On the homogenization and the simulation of random materials. European J. Mech. A Solids 11, 585±607. Stoyan, D., Kendall, W.S., Mecke, J., 1987. Stochastic Geometry and its Applications. Wiley, New York. Taya, M., Arsenault, R.J., 1989. Metal Matrix Composites: Thermomechanical Behavior. Pergamon Press, Oxford.
www.elsevier.com/locate/mechmat
Scale eects in materials with random distributions of needles and cracks Martin Ostoja-Starzewski * Institute of Paper Science and Technology, Georgia Institute of Technology, 500 10th St., N.W., Atlanta, GA 30318-5794, USA Received 14 September 1998; received in revised form 6 May 1999
Abstract According to a classical prescription of micromechanics (R. Hill, J. Mech. Phys. Solids 11, 1963), a representative volume element (RVE) is well de®ned when the response under uniform displacement (Dirichlet) boundary condition becomes the same as that under uniform stress (Neumann) boundary condition. We study the convergence of both responses in anti-plane elasticity of sheets with non-periodic, random distributions of thin needle-shaped inclusions. By lowering the stiness of inclusions and increasing their aspect ratio (up to 100), we approach the situation of cracks embedded in a matrix. We show that, with the needles' stiness decreasing and their slenderness growing, the RVE tends to be very large. The statistics of the ®rst and second invariants of both response tensors are very well modeled by a Beta probability distribution. For moderate aspect ratio needles, the coecient of variation of the second invariant is found to stay at about 0.5 irrespective of the window size, the mismatch in stiness between the inclusions and the matrix, and the needle aspect ratio. Ó 1999 Elsevier Science Ltd. All rights reserved. Keywords: Scale eects; Mesoscale; Micromechanics; Random composites; Cracked media
1. Introduction Dating back to the seminal work of (Budiansky and O'Connell, 1976), the problem of eective moduli of materials with randomly arranged ellipsoidal inclusions or cracks has received considerable attention over the past two decades; for extensive reviews see Kachanov (1993), NematNasser and Hori (1993), Taya and Arsenault (1989). The works on this subject can be classi®ed into two categories: eective medium theories and bounds. Those in the ®rst category include self-
*
Tel.: +1-404-894-5700; fax: +1-404-894-4778. E-mail address: [email protected] (M. Ostoja-Starzewski)
consistent models, Mori±Tanaka models, mean®eld models, and computer simulations. The works in the second category are based on the Hashin±Shtrikman bounds. All these studies strive to establish a macroscopically eective stiness tensor on scales corresponding to a representative volume element (RVE) of an equivalent continuum. However, no information has been obtained on the size of this RVE as a function of the density of the inclusions or cracks, the mismatch (contrast) between the inclusion material and the matrix, and various types of geometric disorders. This paper attempts to answer some of these issues. The study is conducted in the setting of antiplane classical elasticity, that leads directly to twodimensional ®eld problems of two-phase media where one phase is the matrix, while the other is
0167-6636/99/$ - see front matter Ó 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 1 6 7 - 6 6 3 6 ( 9 9 ) 0 0 0 3 9 - 3
884
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
the inclusion. Both phases are locally isotropic with C m and C i being the stinesses of the matrix and the inclusion, respectively. The inclusions are taken as thin needles, generated as germ-grain Boolean processes (Stoyan et al., 1987). We focus on the soft inclusion case (C i =C m < 1), and also try to assess the limit of zero stiness C i =C m ! 0. We investigate the RVE problem in the spirit of Hill's de®nition (Hill, 1963), which says that the RVE is well de®ned when the response under uniform displacement (Dirichlet) boundary condition becomes the same as that under uniform stress (Neumann) boundary condition. While the perfect coincidence of both responses for non-periodic, random microstructures occurs on scales in®nitely larger than the size of an inhomogeneity, their convergence may be anywhere from rapid through slow depending on the parameters. This issue had been studied earlier for 2D and 3D elasticity of several composite materials; see e.g. Amieur et al., 1995; Huet, 1990, 1991, 1997; Ostoja-Starzewski and Wang, 1989; OstojaStarzewski, 1994, 1998a; Ostoja-Starzewski and Schulte, 1996, and references therein.
2. Needle systems 2.1. Random geometry In our two-dimensional problem, needles are taken as narrow d w rectangles, where d and w are the length and width, respectively. In order to generate a ®eld of needles, such as the one shown
in Fig. 1(a), we proceed in the following steps. We employ a Poisson point process to generate the needles' center points. Strictly speaking, the Poisson point process of intensity l is speci®ed by (|A| is a Lebesgue measure of the set A) P fN A kg
kjAj ÿk e jAj; k!
A R2 :
1
It occurs on the entire R2 plane, and so, its numerical simulation must necessarily be done on a compact subset, which actually involves a binomial point process. To specify a needle's orientation we generate an angle h from 0; p according to a probability density f h 1=p, where h is the angle of inclination to the x-axis. This results in a statistically isotropic con®guration of needles, and hence, in a macroscopically isotropic (in the ensemble average sense) response. The subset of R2 on which the binomial process is conducted must be larger than the actual L L window domain which will be tested for apparent moduli ± lines that originated from Poisson points outside the window, located up to a distance d=2, should be accounted for. With L and d at hand, we are ready to introduce a non-dimensional parameter d
L d
2
specifying the scale of observation relative to the size of a heterogeneity, i.e. a mesoscale. Another basic parameter specifying the composite is the aspect ratio of needles d=w. By varying the aspect ratio from 1 up through higher values we can
Fig. 1. (a) A realization of random microstructure of a two-phase, needle-matrix composite material, and the window-scale concept; d is the typical needle length and L the window size. (b) A spring network as a basis for resolution of needles and ellipses.
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
model systems having rectangle-type through needle-type inclusions. A realization of an isotropic ®eld of one hundred needles of length d 100 units and width w 1 unit in a 1000 1000 window is shown in Fig. 2; k 10ÿ4 . Already at this rather low density, random clustering of needles is easily seen. In this paper we focus on two cases: long (d 100) and short (d 10) needles, both at w 1. Another important parameter is the volume fraction of inclusions de®ned by v i V i =Vtotal , where the volume V is really just the area in our planar system. The range through which v i may be varied will depend on the chosen model of placement of inclusions in the matrix. For example, in the case of no-overlap condition of inclusions we would be looking at a random packing problem, while in the case of no such condition ± as implied by the Poisson point process of needle centers ± we can go theoretically up to 100%. However, the latter situation is of no interest per se when studying eective elastic properties. Rather, of interest is the dependence of eective and apparent moduli on the mismatch of individual phases' moduli and on the window
885
scale d, possibly up to the point of percolation of needles. 2.2. Basic mechanics We focus on matrix-inclusion, locally isotropic, elastic materials with inclusions aligned with the x3 , subjected to anti-plane shear. Inclusions are needle shaped. That is, locally (for either phase) the Hooke's law of the matrix (m) and the inclusion (i) material is given in the x1 ; x2 -plane of loading by ri Cij ej ; Cij C
m
i; j 1; 2;
dij
or
C i dij ;
3
where, for simplicity of notation, we denote ri ri3 ;
ei ei3 ;
i; j 1; 2:
4
Parallel to the subscript notation we shall occasionally use the symbolic notation, so that, for example: Cij C. The governing equation of the piecewise-constant material is 2 o u o2 u 0; C C i or C m ; C ox21 ox22 5 u u3 : The isotropy of both stiness tensors C m and in (3) leads to a so-called contrast C a C i =C m , also called a mismatch in other studies. It is clear that by increasing the contrast we can go to very sti inclusions, and approximate a rigid case. Similarly, by decreasing the contrast, we can go to very soft inclusions and approach a system with holes, and, for high d=w, a system with cracks. C i =C m < 1 is, in fact, the case of interest in the present study. Throughout the paper we shall adopt C m 1. Once a speci®c geometry (e.g. Fig. 2) of needles is generated and the assignment of local phase properties to the matrix and needles is made, we are dealing with one realization B x of a random medium B. The latter is taken as a set fB x; x 2 Xg; each B x is treated as a deterministic problem. Fig. 3 depicts the parameter plane ± window scale d L=d and the contrast a C i =C m ± that is investigated in this paper. See Section 4 for further details. i
Fig. 2. A ®eld of long (1 100) needles in a 1000 1000 square, generated from a Poisson point process of density k 10ÿ4 , i.e. there are 100 needles.
886
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
Fig. 3. A parameter plane: window scale d L=d and the contrast a C i =C m . Bold points show the particular cases that were investigated in support of the conclusions in this paper. For long needles (d 100) this translates to d L=100, while for short needles (d 10), to d L=10. Numbers next to the bold points indicate the sample sizes analyzed.
2.3. Numerical solution method In order to solve the ®eld equation (5) of our two-phase composite, with all the interactions, we employ a spring network method, Fig. 1(b). The idea is to approximate the planar, piecewise-constant continuum by a very ®ne square mesh of springs, e.g. (Ostoja-Starzewski and Schulte, 1996). We employ a square mesh for discretization of the displacement ®eld u, so that the governing equations become u i; jkr kl ku kd ÿ u i 1; jkr ÿ u i ÿ 1; jkl ÿ u i; j 1ku
6
ÿ u i; j ÿ 1kd 0: Here i and j are the coordinates of mesh points, and kr ; kl ; ku and kd are de®ned from the series spring model ÿ1
kr 1=C i; j 1=C i 1; j ; ÿ1
kl 1=C i; j 1=C i ÿ 1; j ; ku 1=C i; j 1=C i; j 1ÿ1 ;
7
ÿ1
kd 1=C i; j 1=C i; j ÿ 1 ; where C i; j is the stiness ± i.e., C i or C m ± at the mesh point i; j. The spring network of Fig. 1(b) can be used to model needles of width equal to (or greater than) the mesh spacing, or other inclusions such as ellipses. Just like every numerical method, the spring
network has its pros and cons. Here they are our main considerations. The generation of each and every con®guration of needles is extremely rapid and it does not require any remeshing as the same underlying network 1000 1000, say) is always being used. In particular, no remeshing is needed to handle intersecting needles, as these just happen to appear as certain pixel patterns on the network, a sample of which is shown in Fig. 2. This is an important consideration when dealing with hundreds of realizations (con®gurations) from the sample space X. Further, to better resolve the local stress/strain ®elds within and around the inclusions one may decrease the mesh spacing, just as one would proceed with a ®nite element method (FEM). However, the FEM is limited by a higher price of costly and cumbersome remeshing for each and every new con®guration required in statistical studies. Additionally, we are not concerned with local stress/strain ®eld details as they are known to have little eect on the resulting, eective/apparent moduli. Finally, both methods ± spring networks and ®nite elements ± suer from a need to work with a very large number of nodes to cover the free spaces between the needles. This wasting of memory, and the associated increase in size of the resulting algebraic Ax b problems, can be avoided by employing a boundary element method, whereby the nodes are placed at the needles' and external window boundaries only. We have just such a computer program under development, and we have been able to establish, in speci®c cases, that the dierences between the spring network and the boundary element methods are within 5% for contrasts as low as 10ÿ4 . This accuracy of spring networks is perfectly sucient to rapidly establish the elastic moduli of a number of dierent B x realizations from the random medium B, and their corresponding statistics. 3. Scale-dependent hierarchies of bounds on apparent moduli As stated in the introduction, we assess the apparent moduli for any ®nite window domain of
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
887
mesoscale d depending on whether a uniform strain e0j or a uniform stress r0j is prescribed. In the ®rst case, we should choose essential (or, displacement, Dirichlet) boundary conditions while in the second case, we should choose natural (or, traction, Neumann) boundary conditions. The ®rst setup is ri Cije e0j
under
u x e0j Lj
8x 2 oB;
~
~
8
where rij is the resulting mean (volume average) stress, which leads to a mesoscale stiness C ed (e stands for essential); oB is the boundary of B. The second setup is ei Sijn r0j
under t x r0j nj x 2 oB; ~
~
9
where eij is the resulting mean (volume average) strain, and leads to a mesoscale compliance S nd (n stands for natural). Determination of either second-rank tensor, C ed or S nd , requires three tests. Henceforth, we follow Huet (1990, 1991, 1997), by adopting the terminology ``apparent moduli'' for any mesoscale situation (d ®nite), and reserve ``eective moduli'' for the macroscopic limit (d ! 1). For any realization B x, a window's response on the mesoscale (d ®nite) is, under these de®niÿ1 tions, non-unique ± because Cije 6 Sijn almost surely ± and anisotropic. At this stage, we should note the mixing (and hence, ergodic) as well as strict-sense stationarity property of the Poisson point process underlying our two-phase microstructure (Stoyan et al., 1987). It can then be shown from the variational principles that the ensemble averages of these two tensors provide, with the increasing scale d, an ever tighter pair of bounds on the macroscopic apparent stiness tensor C eff , C R S R
ÿ1
hS n1 i
ÿ1
6 hS nd0 i
ÿ1
6 hS nd i
6 C eff 6 hC ed i 6 hC ed0 i 6 hC e1 i C V
ÿ1
1990; Sab, 1992; Ostoja-Starzewski and Schulte, 1996) to correspond to commensurate partitions of a body domain in the sense that d0 is a multiple of d. However, all the numerical results presented here and those given in our previous studies (as well as those of Huet and his coworkers for 3D problems) speak in favor of (10) being true for arbitrary d0 < d. This is now shown in Appendix A for any microstructures described by wide-sense stationary random ®elds. Hierarchy (10) is depicted in Figs. 4 and 7 for two random needle systems; see Section 4 below for further details. At this point we note that the scale (i.e., d) and contrast (i.e., a) dependencies of both tensors in (10) follow, with very good accuracy, the laws ®rst found for 2-D Bernoulli lattices (Ostoja-Starzewski and Schulte, 1996), namely htr C ed i a0 a1 exp a2 dÿa3 a ; ÿ ÿ1 htr S nd i b0 b1 exp b2 dÿb3 a :
11
8d0 < d: 10
eff
Fig. 4. Eect of decreasing contrast a C i =C m on the divergence of bounds hC e i and hS n iÿ1 , at the window size d 5, for the long needle system one realization of which is shown in Fig. 2; C m 1. The eective stinesses C SC 0:61 and C holes 0:65 computed by Eqs. (B.1) and (B.2) of the selfconsistent and the mean ®eld methods are also shown.
Clearly, C is understood as C 1 C d jd1 . The order relation employed in (10) means that t B t 6 t A t for any vector t 6 0 and two second rank tensors A and B. Furthermore, d and d0 in (10) were shown in the previous works (Huet,
4. Scale-dependent hierarchies of bounds: numerical results We investigate two types of systems: (i) those made of long needles (1 100), and (ii) those made of short needles (1 10). The probability distri-
888
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
butions of mesoscale moduli that are involved in the hierarchy (10) are assessed in terms of the statistics of their two invariants ± half-traces and radii of the corresponding Mohr's circles (C C ed or S nd ) tr C=2
C11 C22 ; 2q
R C12;max
2 C11 ÿ C22 2 =4 C12 :
12
We ®rst focus on the long needles. Fig. 2 displays a typical ®eld of such needles, with a volume fraction 0.01, which is a result of placing 100 needles of size 1 100 in a 1000 1000 window. Fig. 4 displays the eect of decreasing contrast ÿ1 C i =C m on the divergence of bounds hC e ihS n i and at the window size d 5. For the sake of comparison we also indicate the eective stinesses C SC and C holes computed by Eqs. (B.1) and (B.2) of the self-consistent and eective medium models reported in Appendix B. Figs 5 and 6 display, for three scales d 1; 4; and 5, the statistics of tr C=2 and R for tensors C ed and S nd at contrast a 10ÿ4 , respectively. They have been obtained by solving, through a spring
network method of Section 2.3, Dirichlet and Neumann boundary value problems for a number of composites B x of the set B, all being generated in a Monte Carlo sense. Also shown in Figs. 5 and 6 are the Beta function probability density ®ts to the traces and radii; note that the distributions of tr C=2 for C ed are quite dierent from those of S nd . A number of other distributions were also tried ± in particular, Chi, Chi-Square, Exponential, Frechet, Gumbel Min, Log-Normal, Rayleigh and Weibull. However, none of them oered the same overall goodness-of-®t as Beta. It is indeed of interest that Beta is also very satisfactory for other types of composites, both in the elastic (Ostoja-Starzewski, 1998a,b) as well as inelastic (Alzebdeh et al., 1998) response ranges (see Fig. 7). Note that while the traces have asymmetric distributions, their character is very similar for both tensors. The Mohr's circles' radii have even stronger skewness, and, most interestingly, their coecients of variation (CV), for a given type of boundary conditions (either Dirichlet or Neumann), are practically constant (!) with the changing mesoscale d.
e Fig. 5. Probability densities of Ciie =2 (left column) and C12;max (right column) of the needle-matrix composite (at a 10ÿ4 ) whose one realization is shown in Fig. 2 at d 1; 4; and 5.
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
889
n Fig. 6. Probability densities of Siin =2 (left column) and S12;max (right column) of the needle-matrix composite (at a 10ÿ4 ) whose one realization is shown in Fig. 2 at d 1; 4 and 5.
Fig. 7. Eect of increasing widow scale on the convergence of the hierarchy of scale dependent bounds for the composite whose one realization is shown in Fig. 2; C m 1; C i 10ÿ4 . The eective stiness C holes 0:65 computed by Eq. (B.2) of the mean ®eld method is also shown.
e CV C12;max 1:3;
n CV S12;max 0:55:
13
Turning to the short needles, in Fig. 8 we present a typical realization B x. Next, Fig. 9 displays the eect of increasing contrast C i =C m ÿ1 on the divergence of bounds hC e i and hS n i , at the window size d 50. The Poisson process density has been chosen so as to result in roughly the same
Fig. 8. A ®eld of short (1 10) needles in a 1000 1000 square, generated from a Poisson point process of density k 10ÿ4 , i.e. there are 1000 needles.
volume fraction of the soft (i.e., needle) phase as was the case with long needles; here it equals 0.015. The eective stiness C holes 0:43 computed by Eq. (B.2) of the mean ®eld method is also
890
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
C 3 should decrease with delta due to a decrease of C 1 with delta. In eect, (15) allows a determination of any statistic of C 3 from those of C 1 and C 2 . We do not pursue this here. Figs. 4 and 9 lead to two more conclusions. Firstly, the Dirichlet bounds hC ed i are seen to asymptote very quickly with increasing contrast, ÿ1 while the Neumann bounds hS nd i keep falling o.
Secondly, the behavior of both bounds is indicative of the goodness of approximation of a twophase system, with needles having a low albeit ®nite property, relative to a material with holes: we see that the higher the aspect ratio, the closer is the system to the material with holes. In the studies of eective moduli of heterogeneous materials, the resulting C eff is presented graphically versus the volume fraction of one of the phases. For our system of short needles, this is shown in terms of the C holes of Eq. (B.2) of Appendix B against x nL2eff almost all the way to the percolation point at 5.9 (see Fig. 10). In addition, we also plot here the Dirichlet and Nuemann moduli hC ed i and hS nd iÿ1 at two mesoscales: d 10 and 50. Clearly, the three data points of Fig. 9 at C i =C m 10ÿ5 are at x 12:1 (since n 0:01; a 5; b 0:5), while all the data at (two, three and four times) higher volume fractions ± x 2:42; 3:63; and 4.84 ± are new. Thus, this ®gure displays (i) a very slow approach of hC ed i and hS nd iÿ1 to the RVE (i.e., curve C holes ), and (ii) a discrepancy between C holes and the Dirichlet, as well as C holes and the Neumann bounds, respectively. Note that C holes corresponds to an eective (macroscopic) response of a very large random system under periodic boundary conditions (Garboczi et al., 1991).
Fig. 9. Eect of increasing contrast a C i =C m on the divergence of bounds hC e i and hS n iÿ1 , at the window size d 50, for the short needle system one realization of which is shown in Fig. 8; C m 1. The eective stiness C holes 0:43 computed by Eq. (B.2) of the mean ®eld method is also shown.
Fig. 10. Normalized overall moduli hC ed i and hS nd iÿ1 , at d 10 and 50, and the eective stiness C holes from Eq. (B.2), for a random ®eld of short 1 10 needles (such as that of Fig. 8), as functions of the volume fraction x nL2eff . Data were computed only at discrete intervals x 1:21; 2:42; 3:63 and 4.84.
shown. The CVs of Mohr's circles' radii are now not only constant for a given type of boundary condition (either Dirichlet or Neumann) and independent of the changing mesoscale d, but also much closer to each other, e 0:45; CV C12;max
n CV S12;max 0:55:
14
Most interestingly, these numbers are very close to those reported recently for two-phase planar Voronoi mosaics (Ostoja-Starzewski, 1998b), or those found for matrix-inclusion composites, where inclusions are circular. Thus, we conjecture e n is close to CV S12;max for matethat CV C12;max rials with heterogeneities of moderate rather than extreme aspect ratios. The third invariant Cij Cjk Cki also appears to have similar universal properties. In fact, observing that, in 2D problems i; j 1; 2, h i 2 =2 15 C 3 Cij Cjk Cki C 1 3C 2 ÿ C 1
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
5. Closure The problem of scale dependence and of approach to the RVE in heterogeneous materials is a dicult one ± even in the simple setting of twophase composites in anti-plane shear ± and we ®nd the computational mechanics to be the most effective method of attack. In this paper we have investigated scale dependence of apparent moduli in systems having up to one hundred needles of aspect ratio 100, or up to 1000 needles of aspect ratio 10. While arbitrary geometries, from dilute right down to the percolation point, could be handled, with all the interactions being taken into account, we have focused on the volume fraction of about 1% to 1.5%. We believe this study is useful because: (1) The results provide guidance as to the size of the RVE, up to whatever desired discrepancy between the Dirichlet and Neumann type bounds (C ed and S nd ) one might want to choose. This, in turn, provides the key criterion in, say, the choice of ®nite element sizes in solving large (macroscopic) boundary value problems. (2) In the case when the above mentioned discrepancy is substantial and when one is interested in a macroscopic boundary value problem, one may have to choose to work with a stochastic, rather than a deterministic, ®nite element method. The window is then recognized as a mesoscale ®nite element, whose properties are statistical. The Dirichlet and Neumann type stinesses enter, respectively, the stiness and ¯exibility matrices, which themselves are derived from the principles of minimum potential and complementary energy principles, which then jointly result in bounds on macroscopic response, e.g., Ostoja-Starzewski and Wang (1999) and Ostoja-Starzewski (1999). (3) Probability distributions of C ed or S nd have been found to follow, in terms of their ®rst and second invariants, the Beta distribution to a high accuracy. The coecient of variation of the second invariant is seen to behave like a universal material constant: for short needles, it is constant irrespective of the contrast and the mesoscale, and equals 0.45 for the Dirichlet and 0.55 for the Neumann bound. The closeness of these values to those found for other composites (circular disk
891
inclusions in a matrix, or two-phase Voronoi mosaics), also for a wide range of contrasts and mesoscales, is striking. Acknowledgements Extensive comments of Prof. C. Huet (Ecole Polytechnique Federale de Lausanne) on the research reported here are gratefully acknowledged. We also bene®ted from the discussions with Prof. I. Jasiuk (Georgia Tech) and Mr. R. Muzzolini (Ecole des Mines de Paris). This research was made possible by the grant CMS-9713764 from the National Science Foundation. Appendix A. Validity of scale dependent bounds for non-commensurate partitions We want to prove the inequalities (10) for an arbitrary pair of (meso)scales d0 < d, and not just for commensurate ones, i.e., for partitions in which d nd0 ; n being a natural number. It will suce to focus on the Dirichlet bounds, because then the Neumann bounds follow by an analogous argument. We know that the hierarchy holds rigorously for the commensurate sequences of scales. Let us therefore consider two separate cases of the hierarchy (10) for commensurate partitions: one at an arbitrary d1 and another at d2 , whereby d1 < d2 < 2d1 . Thus, in the ®rst case, we have a sequence of inequalities 6 hC e4d1 i 6 hC e2d1 i 6 hC ed1 i 6
A:1
while in the second case, we have another sequence 6 hC e4d2 i 6 hC e2d2 i 6 hC ed2 i 6
A:2
Let us proceed by contradiction: we assume that the inequalities (A.1) and (A.2) are not consistent with each other in the following way, > hC e4d2 i < hC e4d1 i > hC e2d2 i < hC e2d1 i > hC ed2 i < hC ed1 i >
A:3
This, however, would imply that the microstructure is characterized by a Dirichlet bound which,
892
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
while displaying an overall decreasing behavior with the scale d, also ¯uctuates with a period increasing like nd. One can certainly construct random, ergodic microstructures with such scale dependent ¯uctuation. This would, however, contradict the assumption of statistical homogeneity of the material properties that is always taken in the de®nition of an RVE. From the standpoint of stochastic processes, there are two types of the statistical homogeneity: strict-sense stationarity, and wide-sense stationarity. In the ®rst case we are assured of the invariance of any probability distribution of material properties with respect to arbitrary shifts, while in the second case, of the invariance of the mean only with respect to such shifts along with the dependence of two-point correlation functions on the interpoint separations only. The classical RVE de®nition does not say exactly which of these is involved. Nevertheless, considering that the hierarchy (10) is stated in terms of the averages, it suces to choose the setting of material properties speci®ed via wide-sense stationary random ®elds. The ¯uctuation of Dirichlet bound in (A.3) would not even be consistent with the wide-sense stationarity ± as a result, strict-sense stationarity would also be violated. Summarizing, we conclude that the hierarchy (10) holds for non-commensurate mesoscale sequences. Other than (A.3) types of inconsistencies between (A.1) and (A.2) may also be considered, but then one is led to similar contradictory conclusions as above. Appendix B. Eective medium models and bounds As mentioned in the introduction, there is a host of eective medium models as well as formulas for rigorous bounds. Of these the case which is applicable to our system is that of elongated holes in a matrix ± apparently, no results are available for elongated, soft inclusions. First, Nemat-Nasser and Hori (1993) give a self-consistent estimate CSC for a random dilute distribution of frictionless cracks of density f C SC p 1ÿf : C m 2
B:1
As con®rmed by Fig. 4, this formula works well for our long 1 100 needles of negligible stiness, even at contrast 10ÿ3 , but, understandably, cannot be applied to the short needles. However, there exists in the physics literature a result valid for systems with needles of any aspect ratio, and having arbitrarily strong interactions (due to the presence of needle±needle intersections). According to Garboczi et al. (1991), the eective conductivity of a plane with random holes (equivalent, by a mathematical analogy, to antiplane shear) is given by a formula C holes x x x2 1ÿ 1 ÿ C m 5:90 5:90 24:97 1:3 x ÿ1 1 ; x nL2eff ; B:2 3:31 where C1 is the conductivity of the matrix, n is the number of holes per unit area, and Leff is a parameter dependent on the hole geometry. In the case of ellipses, Leff 2 a b, where a and b are the semimajor and semiminor axes. This is a universal formula in the sense that it also holds for other than elliptical holes, and right down to the percolation threshold xc 5:9 0:4. Clearly, it grasps the strong interaction eects very well, its only drawback being that it does not work for holes ®lled with a second phase. Therefore, we employ it for our systems at very low contrast, and, in fact, use it as an indication whether a system at low contrast (soft needles), say 10ÿ4 , is equivalent to the one of zero contrast (voids, cracks). For long needles (Fig. 4) the answer is in the armative, while for the short ones (Fig. 9), one needs to go to contrast 10ÿ5 or lower. Finally, we can consider variational bounds of one kind or another. Since we are dealing here with an isotropic system where the connected phase is stier than the inclusions, we could employ the upper bound for planar two-phase materials (Hashin, 1983). Unfortunately, a quick inspection reveals that this bound is so wide as to be of no practical use compared to the numbers we get via computational mechanics or from the effective medium formula (B.2). Therefore, we do not plot it.
M. Ostoja-Starzewski / Mechanics of Materials 31 (1999) 883±893
References Alzebdeh, K., Al Ostaz, A., Jasiuk, I., Ostoja-Starzewski, M., 1998. Fracture of random matrix-inclusion composites: scale eects and statistics. Int. J. Solids Struct. 35, 2537. Amieur, M., Hazanov, S., Huet, C., 1995. Numerical and experimental assessment of the size and boundary conditions eects for the overall properties of granular composite bodies smaller than the representative volume. In: Parker, D.F., England, A.H. (Eds.), IUTAM Symposium on Anisotropy, Inhomogeneity and Nonlinearity in Solid Mechanics. Kluwer Academic Publisher, Dordrecht, pp. 149±154. Budiansky, B., O'Connell, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct. 12, 81±97. Garboczi, E.J., Thorpe, M.F., DeVries, M.S., Day, A.R., 1991. Universal conductance curve for a plane containing random holes. Phys. Rev. A 43, 6473. Hashin, Z., 1983. Analysis of composite materials ± a survey. J. Appl. Mech. 50, 481. Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. J. Mech. Phys. Solids 11, 357. Huet, C., 1990. Application of variational concepts to size eects in elastic heterogeneous bodies. J. Mech. Phys. Solids 38 (6), 813. Huet, C., 1991. Hierarchies and bounds for size eects in heterogeneous bodies. In: Maugin, G.A. (Ed.), Continuum Models and Discrete Systems, vol. 2, Longman Scienti®c & Technical, pp. 127±134. Huet, C., 1997. An integrated micromechanics and statistical continuum thermodynamics approach for studying the fracture behaviour of microcracked heterogeneous materials with delayed response. Eng. Fract. Mech. 58 (5±6), 459.
893
Kachanov, M., 1993. Elastic solids with many cracks and related problems. Advance Applied Mechanics, vol. 30. Academic Press, New York, p. 259. Nemat-Nasser, S., Hori, M., 1993. Micromechanics: Overall Properties of Heterogeneous Solids. North-Holland, Amsterdam. Ostoja-Starzewski, M., Wang, C., 1989. Linear elasticity of planar Delaunay networks: random ®eld characterization of eective moduli. Acta Mech. 80, 61. Ostoja-Starzewski, M., 1994. Micromechanics as a basis of continuum random ®elds. Appl. Mech. Rev. 47 (1, Pt. 2), 221. Ostojas-Starzewski, M., Schulte, J., 1996. Bounding of eective thermal conductivities of multiscale materials by essential and natural boundary conditions. Phys. Rev. B 54 (1), 278. Ostoja-Starzewski, M., 1998a. Random ®eld models of heterogeneous materials. Int. J. Solids Struct. 35, 2429. Ostoja-Starzewski, M., 1998b. Random ®eld models and scaling laws of heterogeneous media. Arch. Mech. 50 (3), 549. Ostoja-Starzewski, M., 1999. Microstructural Disorder, Mesoscale Finite Elements, and Macroscopic Response. Proc. Roy. Soc. Lond. A 455, 3189±3199. Ostoja-Starzewski, M., Wang, X., 1999. Stochastic ®nite elements as a bridge between random material microstructure and global response. Comp. Meth. Appl. Mech. Eng. 168 (1±4), 35. Sab, K., 1992. On the homogenization and the simulation of random materials. European J. Mech. A Solids 11, 585±607. Stoyan, D., Kendall, W.S., Mecke, J., 1987. Stochastic Geometry and its Applications. Wiley, New York. Taya, M., Arsenault, R.J., 1989. Metal Matrix Composites: Thermomechanical Behavior. Pergamon Press, Oxford.