International Journal of Engineering Science 58 (2012) 21–34
Contents lists available at SciVerse ScienceDirect
International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci
Effects of particle shape on the macroscopic and microscopic linear behaviors of particle reinforced composites Azra Rasool, Helmut J. Böhm ⇑ Institute of Lightweight Design and Structural Biomechanics, Vienna University of Technology, Vienna, Austria
a r t i c l e
i n f o
Article history: Received 14 November 2011 Accepted 16 January 2012 Available online 19 April 2012 Keywords: Composites Thermoelastic properties Conduction properties Particle shape effects
a b s t r a c t A systematic comparison of inhomogeneity shape effects on the linear elastic, thermoelastic and thermal conduction responses of particle reinforced composites is carried out. For this purpose, multi-particle unit cells that contain randomly positioned and, where applicable, oriented, identical particles having the shapes of spheres, regular octahedra, cubes or regular tetrahedra, respectively, and a volume fraction of 20% are employed. The macroscopic moduli and microscopic responses, such as phase averages, as well as phase-level standard deviations and distribution functions of the microfields are evaluated and compared to analytical estimates. The results indicate the presence of relatively small but consistent effects of the particle shape on the effective behavior of particulate composites. Effects on the microscopic stress and flux fields are predicted to be more pronounced. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction A well-known mean field expression for the effective elastic tensor of two-phase composites, E⁄, can be written in the form
E ¼ nðmÞ EðmÞ AðmÞ þ nðiÞ EðiÞ AðiÞ :
ð1Þ ⁄
(m)
(i)
From this equation it is evident that the leading parameters determining E are the elastic tensors, E and E , as well as the volume fractions, n(m) and n(i), of the matrix (m) and the inhomogeneity (i) phases, respectively. The phase averaged strain concentration tensors, AðmÞ and AðiÞ , contain information, on the one hand, on the macroscopic symmetry of the materials, i.e., on the two-point statistics of the phase arrangement. On the other hand, they also depend on details of the microgeometry, such as clustering, size distributions and shape(s) of the reinforcement phase. The latter descriptors can be captured by higher correlations and have more limited effects on the effective tensors. Analogous considerations hold for the thermoelastic and thermal conduction responses of inhomogeneous materials. The present work concentrates on one of these details, particle shape, in particle reinforced composites. Among the established methods for modeling the behavior of particle reinforced composites, the classical equivalent inclusion approach of Eshelby (1957) and the standard mean field methods based on it are capable of handling ellipsoidal inhomogeneities of different aspect ratios. Non-ellipsoidal inclusions are known to give rise to Eshelby tensors that are position dependent (Kang & Milton, 2008). Accordingly, extensions of mean field models to more general inhomogeneity shapes have, on the one hand, been based on inclusion averaged Eshelby tensors (Gao & Liu, 2012; Nozaki & Taya, 2001; Onaka, 2001; Rodin, 1996; Zheng, Zhao, & Du, 2006). On the other hand, numerically evaluated compliance contribution (Eroshkin & Tsukrov, 2005; Kachanov, Tsukrov, & Shafiro, 1994; Sevostianov, Kachanov, & Zohdi, 2008) or replacement concentration ⇑ Corresponding author. E-mail address:
[email protected] (H.J. Böhm). 0020-7225/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijengsci.2012.03.022
22
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
(Nogales & Böhm, 2008) tensors have been combined with mean field methods for studying composites reinforced by nonellipsoidal inhomogeneities. These two groups of approaches can be directly applied to materials reinforced by aligned particles and can be combined with orientational averaging procedures to obtain estimates on the effective behavior of composites reinforced by nonaligned reinforcements (Hashemi, Avazmohammadi, Shodja, & Weng, 2009). Bounds and estimates based on three-point statistics (Torquato, 2002) are in principle capable of handling particle shape effects. However, at present the availability of three-point microstructural parameters is largely restricted to spherical inhomogeneities. All of the above methods cannot provide information on the microfields that goes beyond phase averages. Micromechanical approaches that numerically evaluate the microfields of discrete phase arrangements, e.g., by the Finite Element (FE) method, can handle general particle shapes in a rather straightforward way and provide predictions for both the macroscopic and microscopic responses. Although work on particle reinforced composites making use of periodic homogenization and related methods has predominantly concentrated on spherical inhomogeneities, a number of studies involving non-spheroidal particles have been published. These include axisymmetric cell models, compare (Finot, Shen, Needleman, & Suresh, 1994), and unit cell models of cubic arrangements of cylindrical and cube-shaped inhomogeneities, see, e.g., (Meijer, Xia, & Ellyin, 1997; Weissenbek, Böhm, & Rammerstorfer, 1994). Recently, a number of contributions involving synthetic multi-particle volume elements containing randomly positioned and oriented polyhedral inhomogeneities (Asahina & Bolander, 2011; Fritzen & Böhlke, 2011; Galli, Botsis, & Janczak-Rusch, 2008; Han, Eckschlager, & Böhm, 2001; Hu & Chen, 2011; Nogales & Böhm, 2008; Shen & Lissenden, 2003; Zhu, Cai, & Tu, 2009) were published. In addition, for ‘‘real structure models’’, which are directly based on microgeometries obtained by tomography (Kenesei, Biermann, & Borbély, 2006; Mishnaevsky, Dong, Hönle, & Schmauder, 1999; Terada & Kikuchi, 1996) or serial sectioning (Chawla, Sidhu, & Ganesh, 2006; Lippmann, Steinkopff, Schmauder, & Gumbsch, 1996) non-spheroidal particles are the rule rather than the exception, and particles of non-convex shape may also be present. However, few of the above models, among them (Han et al., 2001; Nogales & Böhm, 2008; Zhu et al., 2009), were used for comparing the effects of different particle shapes, and, at this time, there appear to be no studies in the literature that report on systematic comparisons of this type. For this reason, the present contribution aims at elucidating some aspects of the influence of inhomogeneity shapes on the macroscopic and microscopic, thermoelastic and thermal conduction behavior of particle reinforced two-phase composites, i.e., materials of matrix-inclusion topology that contain (approximately) equiaxed inhomogeneities. 2. Methods 2.1. Modeling strategy The modeling strategy underlying the present study is based on employing relatively simple models that are nevertheless suitable for obtaining detailed information on particle shape effects on the behavior of particle reinforced composites. In keeping with this philosophy relatively small synthetic volume elements are studied by standard numerical methods. The discussion is restricted to macroscopically isotropic, inhomogeneous two-phase materials showing matrix-inclusion microtopology, for which influences of the particle shape cannot be masked by effects of macroscopic anisotropy. Of course, such an approach implies that inhomogeneities must be randomly positioned in the volume element and that particles of nonspherical shape must be randomly oriented. The effective elastic and conduction properties of such composites fulfill the classical Hashin–Shtrikman bounds (Hashin & Shtrikman, 1962, 1963), i.e., all considered phase arrangements follow essentially the same two-point statistics but differ in their three-point statistics. Furthermore, high quality reference solutions for the macroscopic behavior are available in the form of three-point bounds (Torquato, 1991) and three-point estimates (Torquato, 1998) for the special case of composites reinforced by randomly positioned identical spheres. For the phase contrasts considered here, the latter methods are sufficiently sensitive for assessing the quality of numerical results pertaining to random arrangements of identical spheres. The non-spherical particle shapes studied encompass regular octahedra, cubes and regular tetrahedra, i.e., a selection of equiaxed and convex, simple polyhedra. All of the above shapes pffiffiffi differ markedly from spheres, with the ratios of the volumes of circumscribed and inscribed spheres ranging between 3 3 (octahedra and cubes) and 27 (tetrahedra). For this reason bounds on the overall behavior of the resulting composites making direct use of Hill’s comparison theorem (Hill, 1963) as proposed by Huet (1990) are limited to low particle volume fractions and are too slack for practical use. All work reported in the following pertains to equally sized particles of a total volume fraction of n(i) = 0.2. Generic, isotropic material properties are prescribed for the composites’ constituents and the elastic, thermoelastic as well as thermal conduction responses are evaluated. The material parameters listed in Table 1 are used unless explicitly specified otherwise,
Table 1 Thermoelastic and conduction material parameters used in the simulations. E [GPa] Particles Matrix
10.0 1.0
m [] 0.10 0.33
a [K1] 6
10 105
k [Wm1K1] 10.0 1.0
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
23
with E standing for the Young’s modulus, m for the Poisson number, a for the coefficient of thermal expansion (CTE) and k for the thermal conductivity of a given constituent. These phase-level material properties correspond to an elasticity contrast of cE = E(i)/E(m) = 10 and a conductivity contrast of ck = k(i)/k(m) = 10. The interfaces between particles and matrix are assumed to be perfect. All volume elements discussed in the following consist of periodic arrangements of randomly positioned and, where applicable, randomly oriented particles that are embedded in the matrix. Different sets of particle positions were employed for each particle shape, the use of identical centers for particles of different shape as proposed in (Han et al., 2001; Nogales & Böhm, 2008) being impractical for the particle shapes and volume fraction considered. In order to limit the analysts’ effort and the computational costs, the responses of periodic volume elements containing some 21 particles each were evaluated. Although such phase arrangements can be expected to describe the targeted phase arrangements reasonably well (Drugan & Willis, 1996), they are too small to be proper Representative Volume Elements (RVEs) for composites reinforced by randomly positioned and oriented particles, giving rise to apparent moduli. Accordingly, ensemble averaging over a number of statistically equivalent volume elements was used to obtain the effective moduli as well as improved estimates of the microscopic behavior. In addition to allowing the use of relatively small model geometries, this approach has the advantage of providing measures such as standard deviations of the effective tensors that allow assessing the quality of the results extracted from a set of volume elements (Kanit, Forest, Galliet, Mounoury, & Jeulin, 2003). In addition to predicting the macroscopic moduli of the model composites the present work aims at obtaining information on the stress, strain, flux and gradient fields at the microscale, which are referred to as microfields in the following. 2.2. Generation and homogenization of volume elements The cube-shaped periodic volume elements, which are also referred to as multi-particle unit cells, were generated by a modified random sequential addition algorithm (Widom, 1966), in which the polyhedra were checked for overlap and for the smallest distance between neighboring particles by the separating axis method (Schneider & Eberly, 2003). A finite minimum distance between neighboring particles and conditions on intersections between particles and the boundaries of the unit cell were prescribed in order to ensure that at least two ‘‘layers’’ of finite elements can be provided in the matrix between neighboring particles and to avoid geometries that are awkward to mesh. The particle volume fractions that can be achieved by this approach depend fairly strongly on the particle shape; for the cases considered here the value of n(i) = 0.2 could be reached without major difficulty. Volume elements containing 20 to 22 particles were generated this way. Fig. 1 shows typical multi-particle unit cells containing spherical (SPH), octahedral (OCT), cube-shaped (CUB) and tetrahedral (TET) particles at this volume fraction. Five different multi-particle unit cells of this type were generated per particle shape; these five realizations are referred to as ‘‘equivalent volume elements’’ in the following. Meshing of the configurations was carried out with the preprocessor program HyperMesh (Altair Engineering Inc., Troy, MI). Ten-node quadratic tetrahedral elements were used for elastic, thermoelastic and conduction analysis, and phase boundaries always coincinded with element boundaries. The meshes at pairs of opposite faces of the volume elements were enforced to be compatible. Convergence tests were carried out in terms of the predictions for the macroscopic elastic responses of the unit cells, for details see (Rasool, 2011). On this basis a reference mesh size of 0.03 times the side length of the unit cells was selected, resulting in total element counts of the order of 150,000 per volume element. The macroscopic responses as well as the microfields within the volume elements were evaluated by Finite Element based periodic homogenization using the method of macroscopic degrees of freedom (Michel, Moulinec, & Suquet, 1999; Smit, Brekelmans, & Meijer, 1998). The geometrically and materially linear computations were carried out with the commercial FE code ABAQUS/Standard (3DS, Dassault Systèmes, Providence, RI). Periodicity boundary conditions were applied via multi-point constraint equations acting on homologous nodes on pairs of opposite faces of the volume elements. Six linearly independent elastic load cases (uniaxial tension in each coordinate direction and simple shear in three mutually normal planes), one thermoelastic load case (a uniform temperature excursion from the stress-free state) and three linearly independent conduction cases (prescribed flux in each coordinate direction) were evaluated. This way the linear elastic, thermoelastic and conduction behavior of each volume element was fully described. 2.3. Processing of results The above load cases allowed extracting the macroscopic elastic, thermal expansion and thermal conduction tensors for each multi-particle unit cell. Due to the relatively small size of the volume elements these effective tensors showed minor deviations from macroscopic isotropy. As can be seen from Table 2, which pertains to a material reinforced by randomly oriented tetrahedra, these deviations from isotropy of the ‘‘raw’’ effective elasticity tensor, E;raw TET , were not fully removed by ensemble averaging. For the sizes and numbers of equivalent volume elements used in the present study, the anisotropic contributions to these elasticity tensors were found to be of comparable magnitude to the standard deviations from ensemble averaging, i.e., of the order of 1% of the diagonal elements. Similar behavior was observed with respect to the effective thermal expansion and thermal conduction tensors. The knowledge that the effective tensors of the materials considered must be isotropic can be made use of for further refining the estimates. In the case of the elasticity tensors this was done by first reducing the raw tensors to orthotropic ones. For this purpose, the coupling terms were set to zero and the remaining tensor elements were averaged in accordance with
24
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
Fig. 1. Typical periodic multi-particle unit cells, each containing some 21 identical tetrahedral (TET), cube-shaped (CUB), octahedral (OCT) or spherical (SPH) particles of volume fraction n(i) = 0.2.
orthotropic symmetry. The resulting orthotropic elasticity tensors, in turn, were ‘‘isotropized’’ by extracting the closest isotropic tensors using the algorithm of He and Curnier, 1995. The effective elasticity tensor listed in Table 2, E;iso TET , is a result of this set of operations. Only the ‘‘zeroing out and averaging’’ process was applied to the thermal expansion and conduction tensors, which are of order 2. All effective moduli listed in Section 3 were evaluated from isotropized effective tensors.
Table 2 Predictions for the elasticity tensor of a composite reinforced by randomly oriented identical tetrahedral particles. Top: Raw tensor, ETET⁄,raw, obtained by ensemble averaging the results from five equivalent volume elements containing some 21 particles each. Bottom: isotropized elasticity tensor, ETET⁄,iso. 0
2:028 0:009 0:854 0:002 0:857 0:004 B 0:854 0:002 2:028 0:006 0:853 0:004 B B 0:857 0:004 0:853 0:004 2:026 0:005 ;raw ETET ¼B B 0:000 0:002 0:002 0:003 0:002 0:007 B @ 0:003 0:004 0:000 0:003 0:000 0:002 0:000 0:002 0:003 0:003 0:001 0:002 1 0 2:030 0:853 0:853 0:000 0:000 0:000 B 0:853 2:030 0:853 0:000 0:000 0:000 C C B B 0:853 0:853 2:030 0:000 0:000 0:000 C ;iso C ETET ¼B C B B 0:000 0:000 0:000 0:588 0:000 0:000 C @ 0:000 0:000 0:000 0:000 0:588 0:000 A 0:000 0:000 0:000 0:000 0:000 0:588
0:000 0:002 0:002 0:003 0:002 0:007 0:587 0:005 0:000 0:004 0:001 0:003
0:003 0:004 0:000 0:003 0:000 0:002 0:000 0:004 0:592 0:003 0:001 0:002
1 0:003 0:003 0:001 0:002 C C 0:000 0:002 C C 0:001 0:003 C C 0:001 0:002 A 0:590 0:003
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
25
In addition to extracting the effective tensors and moduli, the statistics of the microscopic fields were evaluated at the levels of the constituents and of the individual particles. Only a limited amount of work appears to be available in the literature that bears on the question whether the elastic fields at the vertices and edges of polyhedral inhomogeneities embedded in a matrix show singularities. Related configurations were, however, reported to display singular behavior (Ghahremani, Hutchinson, & Tvergaard, 1990; Glushkov, Glushkova, & Lapina, 1998), and Rodin (1996) showed that singularities are typically present at vertices formed by the intersection of three orthogonal faces. The FE meshes employed in the present contribution are obviously limited in their ability to fully resolve strong local gradients of the stress, strain, flux and gradient fields in the vicinity of edges and vertices. They definitely cannot provide detailed information on singularities. Whereas details of the shapes of edges and vertices, and thus details of the local fields at such features, are known to have little impact on the predicted macroscopic responses of composites of the type considered here (Kachanov & Sevostianov, 2005), some degradation of the data evaluated for the microfields must be expected to result from such issues. Phase averages, hfi, of the microfields at the phase and particle levels were evaluated by approximate numerical integration, in which volume integrals are replaced by weighted sums
hf i ¼
1
X
Z
f ðxÞdX
X
N 1P
f Xl :
X l¼1 l
ð2Þ
Here fl and Xl are the function value and the integration weight (in terms of the volume of the integration point), respectively, associated with the l-th integration point within a given integration volume X that contains N integration points. In analogy, the standard deviation of microfield f(x) over the volume X was approximated as
Sðf Þ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N 1P ðfl hf iÞ2 Xl :
X l¼1
ð3Þ
In evaluating the phase averages and standard deviations listed in Section 3 the volumes X were selected as the union of all appropriate phase volumes from the five unit cells pertaining to a given particle shape. In analogy to the evaluation of the macroscopic responses, use was made of the fact that a number of the load cases considered in generating the effective tensors are equivalent for isotropic materials, an example being the microscopic ‘‘axial’’ stress components rxx, ryy and rzz obtained by loading a given volume element by applied uniaxial stresses in x-, y- and z-directions, rax , ray and raz , respectively. Results obtained from such sets of ‘‘pertinent load cases’’ are included in the data on phase-level and particle-level averages and standard deviations presented in Section 3. Density distribution functions of stress, strain, flux and gradient components and measures in either phase (in the sense of ‘‘stress spectra’’ (Bornert, Hervé, Stolz, & Zaoui, 1994)) were obtained by binning the function values at the integration points, fl, and weighting them with the corresponding integration point volumes Xl. Again, contributions from all equivalent volume elements and all pertinent load cases were combined. The inter-particle and intra-particle fluctuations of the microfields can be visualized on the basis of the averages and standard deviations over the individual particles. In contrast to work on composites reinforced by nonaligned short fibers (Duschlbauer, Böhm, & Pettermann, 2006), in the present study the particles are not arranged with respect to their orientation, but rather in terms of the average values of the fields concerned. Each plot represents results obtained from all equivalent multi-particle unit cells and for all pertinent load cases. All data on microscopic stresses and fluxes given in Section 3 is normalized with respect to suitable macroscopic stresses and fluxes, respectively. 3. Results and discussion 3.1. Macroscopic responses The normalized effective elastic moduli evaluated by periodic homogenization for equally sized spherical (SPH), octahedral (OCT), cube-shaped (CUB) and tetrahedral (TET) particles are listed in Table 3. These results were extracted from ensemble averages of the isotropized elastic tensors obtained from the five equivalent volume elements studied for each particle shape. In addition, predictions generated with the Hashin–Shtrikman bounds (HSB), three-point bounds (3PB) and threepoint estimates (3PE) are given. The latter results were obtained with three-point microstructural parameters pertaining to non-interpenetrating spheres of equal size (Torquato, 2002) and thus can be directly compared to the unit cell results for spherical particles. G and K stand for the shear and bulk moduli, respectively. All unit-cell based predictions fulfill the Hashin–Shtrikman bounds and, as appropriate for materials consisting of stiff inhomogeneities in a more compliant matrix (Torquato, 1991), closely approach the lower bounds. The ensemble averaged moduli predicted for the spherical inhomogeneities lie within the pertinent three-point bounds and nearly coincide with the corresponding three-point estimates. Together with the small standard deviations, which remain below 1% of the corresponding moduli, this indicates a satisfactory quality of the results on the effective elastic responses. Fig. 2 shows a plot in the K⁄–G⁄ plane (Berryman & Milton, 1988) of the two-point and three-point bounds as well as the numerical estimates obtained for the four particle shapes. It is evident that, among the particle shapes considered in the present work, the effective elastic responses of the tetrahedra differ most from those of the spheres. Octahedral and cube-shaped inhomogeneities
26
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
Table 3 Normalized effective elastic moduli predicted from the Hashin–Shtrikman bounds (HSB) for macroscopically isotropic two-phase materials, the three-point bounds (3 PB) and three-point estimates (3PE) for composites reinforced by identical spheres, and from multi-particle unit cells containing spherical particles (SPH) or randomly oriented octahedral (OCT), cube-shaped (CUB) or tetrahedral (TET) particles of equal size. Standard deviations obtained from five equivalent volume elements each are listed for the unit cell results. All effective moduli are normalized by the matrix Young’s modulus.
m⁄
G⁄/G(m)
K⁄/K(m)
HSB
1.406–2.152
0.214–0.335
1.399–2.233
1.458–1.738
3 PB 3PE SPH
1.417–1.499 1.427 1.428 ± 0.009
0.295–0.309 0.305 0.305 ± 0.001
1.412–1.503 1.422 1.422 ± 0.008
1.460–1.483 1.462 1.462 ± 0.003
OCT CUB TET
1.468 ± 0.010 1.478 ± 0.010 1.525 ± 0.008
0.300 ± 0.001 0.301 ± 0.001 0.296 ± 0.001
1.479 ± 0.008 1.466 ± 0.008 1.529 ± 0.005
1.475 ± 0.003 1.472 ± 0.003 1.495 ± 0.003
NORMALIZED EFFECTIVE SHEAR MODULUS []
E⁄/E(m)
0.85
TET CUB OCT SPH
0.8
0.75
0.7
0.65
0.6
0.55 1.25
1.3
1.35
1.4
1.45
NORMALIZED EFFECTIVE BULK MODULUS [] Fig. 2. Normalized effective bulk moduli K⁄/K(m) and shear moduli G⁄/G(m) predicted for composites reinforced by identical spherical particles (SPH, circle) or identical randomly oriented octahedral (OCT, triangle), cube-shaped (CUB, pentagon) or tetrahedral (TET, square) particles of volume fraction of n(i) = 0.2 and elastic contrast cE = 10. The pertinent Hashin–Shtrikman bounds are shown as a dotted box and the three-point bounds for equally sized spherical particles as a dashed box.
give rise to fairly similar behavior between the extremes. The present periodic models predict that at the particle volume fraction of n(i) = 0.2 both the bulk and the shear moduli are higher for polyhedral than for spherical particles. Table 3 shows the effect to be considerably stronger for G⁄ than for K⁄, the maximum increase in an effective elastic modulus being found for GTET , which is about 7.5% higher than GSPH . These results differ from those of a Mori–Tanaka-based study (Hashemi et al., 2009) that predicted shape-independent bulk moduli for composites containing randomly oriented superspherical particles. Table 4 gives the corresponding results for the normalized effective coefficients of thermal expansion, a⁄/a(m), and thermal conductivities, k⁄/k(m). Again, good agreement with the bounds and analytical estimates can be observed, the effective CTEs being close to the upper bounds as expected for the material parameters listed in Table 1. The thermal conduction behavior is qualitatively similar to the elastic one, with the predictions pertaining to cube-shaped and octahedral particles
Table 4 Normalized effective coefficients of thermal expansion and effective conductivities predicted from the Hashin–Shtrikman bounds (HSB) for macroscopically isotropic two-phase materials, the three-point bounds (3 PB) and three-point estimates (3PE) for composites reinforced by identical spheres, and from multi-particle unit cells containing spherical particles (SPH) or randomly oriented octahedral (OCT), cube-shaped (CUB) or tetrahedral (TET) particles of equal size.
a⁄/a(m)
k⁄/k(m)
HSB
0.612–0.773
1.53–2.34
3 PB 3PE SPH
0.756–0.772 0.770 0.770
1.54–1.62 1.55 1.54
OCT CUB TET
0.764 0.762 0.749
1.60 1.61 1.71
27
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
lying between those obtained for spherical and tetrahedral inhomogeneities. Quantitatively, the particle shape effects in conduction are found to be rather more pronounced than in elastostatics, the conductivity for the composite reinforced by randomly oriented tetrahedra being predicted to exceed that of the sphere-reinforced material by more than 10%. 3.2. Phase-level microfields Predictions for the phase averages and, where applicable, the corresponding phase-level standard deviations of some typical descriptors for the microscopic stress fields in the composites under study are listed in Table 5 for the macroscopic load case of an applied uniaxial tensile stress, ran . The data includes results obtained from Benveniste’s (Benveniste, 1987) Mori– Tanaka method (MTM) and from Torquato’s three-point estimates, both pertaining to spherical particles. For the latter method the phase averages were evaluated by inserting the effective and phase compliance tensors, C⁄, C(i) and C(m), into mean field expressions closely related to Eq. (1),
C ¼ CðmÞ þ nðiÞ ½CðiÞ CðmÞ BðiÞ ¼ CðiÞ þ nðmÞ ½CðmÞ CðiÞ BðmÞ ;
ð4Þ ðiÞ
ðmÞ
and solving for the phase averaged stress concentration tensors, B and B , respectively. Results are provided for the axial ðpÞ ðpÞ ðpÞ stress component, rkk , the equivalent stress, reqv , and the volumetric stress, rvol , all of which are normalized by the maga nitude of the applied uniaxial stress, rn , where (p) refers to the phases, inhomogeneities and matrix. Excellent agreement in the phase averages of the axial stress component and the volumetric stress is evident between the three-point estimates and the unit cell results pertaining to identical spheres. This, of course, is to be expected from the nearly identical effective moduli by considering Eq. (4). The clear difference in the predictions for the phase averaged equivðpÞ alent stresses of the sphere-reinforced composites is due to the fact that the reqv quoted for the 3PE and MTM models were evaluated from the phase averages of the stress components and, accordingly, do not account for the effects of local stress fluctuations, compare (Buryachenko, 1996). The phase averages in the inhomogeneities (i) can be seen to be substantially higher for tetrahedral than for spherical particles, whereas the phase averages in the matrix ( m) show the opposite trend, albeit in a much less marked way. This behavior gives rise to the increase in the effective moduli predicted for reinforcing by polyhedral particles evident in Table 3. Octahedral and cube-shaped particles again lead to responses that lie between those predicted for spherical and tetrahedral inhomogeneities. These results clearly indicate that changes in the particle shape lead to changes in the stress partitioning between particles and matrix, with the former carrying increasingly high average stresses as their shape deviates more and more from that of a sphere. Mori–Tanaka expressions typically provide low estimates for the average stresses in the particles and high ones in the matrix. According to Table 5 the predicted phase-level standard deviations of the equivalent stress in the matrix show hardly any shape dependence for this load case, and for the volumetric stress and the axial stress component the standard deviations in the matrix are somewhat smaller for the polyhedral particles than for the spherical ones. In contrast, the standard deviations of these stress components and measures in the particles are predicted to be markedly higher for the polyhedra. This behavior is reflected in Fig. 3, which displays the predicted phase-level density distributions of the normalized axial stress comðpÞ ponent, rkk =ran , in particles and matrix for loading by uniaxial macroscopic stresses, ran . Whereas the distributions obtained for the matrix (dashed lines) are very similar for the four inhomogeneity shapes considered, the ones evaluated for the particles spread out dramatically as the particle shape is changed from spheres to tetrahedra. The density distributions evaluated for polyhedral particles generally not only tend to show marked high-stress tails (for all stress spectra shown, very small but nonzero contributions are typically present for stresses outside the plots’ range for the polyhedral particles), but also a higher probability of finding axial stresses smaller than the applied stress. The shape dependence of the phase averages (represented by vertical lines) in particles and matrix also comes out clearly in the plots. Table 5 ðpÞ ðpÞ ðpÞ Phase averages and phase-level standard deviations in phase (p) of the normalized axial stress rkk =ran , equivalent stress reqv =ran and volumetric stress rvol =ran predicted for composites reinforced by spherical (SPH) or randomly oriented octahedral (OCT), cube-shaped (CUB) or tetrahedral (TET) particles of volume fraction n(i) = 0.2 and elastic contrast cE = 10 that are subjected to an applied uniaxial stress ran . In addition, results obtained with the Mori–Tanaka estimates (MTM) and the three-point estimates (3PE) are given. Results are based on all equivalent volume elements and all pertinent load cases.
rðpÞ =ran kk
a rðpÞ eqv =rn
rðpÞ =ran vol
MTM
(i) (m)
1.514 0.871
1.641 0.840
0.420 0.312
3PE
(i) (m) (i) (m)
1.561 0.860 1.563 ± 0.216 0.859 ± 0.262
1.703 0.824 1.734 ± 0.228 0.858 ± 0.193
0.426 0.310 0.424 ± 0.093 0.311 ± 0.152
(i) (m) (i) (m) (i) (m)
1.638 ± 0.467 0.840 ± 0.255 1.655 ± 0.481 0.836 ± 0.255 1.767 ± 0.683 0.08 ± 0.240
1.881 ± 0.466 0.837 ± 0.190 1.921 ± 0.489 0.832 ± 0.191 2.127 ± 0.664 0.801 ± 0.192
0.437 ± 0.229 0.307 ± 0.146 0.440 ± 0.230 0.307 ± 0.145 0.464 ± 0.335 0.301 ± 0.131
SPH OCT CUB TET
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
TET
2.5
PROBABILITY DENSITY []
PROBABILITY DENSITY []
28
2 1.5 1 0.5 0
CUB
2.5 2 1.5 1 0.5 0
0
1
2
3
4
5
0
OCT
2.5
1
2
3
4
5
NORMALIZED AXIAL NORMAL STRESS []
PROBABILITY DENSITY []
PROBABILITY DENSITY []
NORMALIZED AXIAL NORMAL STRESS []
2 1.5 1 0.5 0
SPH
2.5 2 1.5 1 0.5 0
0
1
2
3
4
5
NORMALIZED AXIAL NORMAL STRESS []
0
1
2
3
4
5
NORMALIZED AXIAL NORMAL STRESS [] ðpÞ
Fig. 3. Density distributions of the normalized axial stress component rkk =ran in phase ( p) predicted for composites reinforced by spherical (SPH) or randomly oriented octahedral (OCT), cube-shaped (CUB) or tetrahedral (TET) particles of volume fraction n(i) = 0.2 and elastic contrast cE = 10 that are subjected to an applied uniaxial stress ran . Dashed lines denote matrix and solid lines particle behavior, thin vertical lines indicate corresponding phase averages. Results are based on all equivalent volume elements and all pertinent load cases.
ðpÞ
Fig. 4 shows the phase-level density distributions predicted for the normalized volumetric stress, rvol =ras , under loading by a macroscopic shear stresses, ras . For this macroscopically purely deviatoric load case the phase averages of the volumetric stress are zero in both matrix and particles for all configurations considered. At the microscale, however, finite volumetric contributions to the local stress tensors are obviously present, leading to fairly wide distributions that closely approach symmetry with respect to zero volumetric stress. In this case, too, the matrix stress spectra are nearly independent of the particle shape, whereas the spectra for the particles are much sharper for the spherical than for the tetrahedral particles. Composites subjected to unconstrained thermal expansion typically display non-zero self-equilibrated microstresses and, in the case of macroscopically isotropic materials, purely volumetric phase averaged stress tensors. It is well known that, nevertheless, the phase averages of the equivalent stress are nonzero due to local fluctuations of the stress tensor. The corðpÞ responding density distributions of the equivalent stress in matrix and particles, reqv , are presented in Fig. 5 for the particle reinforced composites considered here. The thermal stress of the matrix material under fully constrained conditions, E(m)a(m)DT, is used for normalizing the equivalent stress, with DT standing for a spatially homogeneous temperature excursion. A marked shape dependence of both the density distributions and the phase averages of the equivalent stress is evident for the particles. The phase averages of the equivalent stress in the matrix show only a limited dependence on the particle shape. Whereas the phase average of the particles is predicted to be slightly lower than that of the matrix for the spherical inhomogeneities, the former is much higher in the case of the tetrahedra. Interestingly, a clear secondary peak of the distribution is predicted for the matrix reinforced by spherical inhomogeneities, which degrades into a rather inconspicuous ‘‘shoulder’’ in the case of tetrahedral particles. This secondary peak is due to a concentric matrix region of elevated equivalent stress surrounding each spherical particle (Rasool, 2011), the spherical symmetry of the stress tensor around each inhomogeneity being perturbed only mildly by neighboring inhomogeneities. Polyhedral particles, in contrast, break the symmetry of the stress states in this matrix region, resulting in a much less pronounced signature in the stress spectra. ðpÞ Fig. 6 depicts the predicted density distributions of the normalized axial flux components, qk =qan , under loading by a maca a roscopic thermal flux, qn , of magnitude qn . The behavior of these microfields can be seen to be qualitatively similar to that of the axial stresses under uniaxial tensile loading in the elastic range, Fig. 3. Taken together, these results establish a clear pattern, which also applies to other stress and flux components (Rasool, 2011). The phase averages and density distributions of the microfields tend to be fairly independent of the particle shape for the matrix phase, but to show a marked dependence for the inhomogeneity phase. This indicates that particles with
29
TET
3 2.5 2 1.5 1
PROBABILITY DENSITY []
PROBABILITY DENSITY []
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
0.5
CUB
3 2.5 2 1.5 1 0.5
0
0 -2
-1
0
1
2
-2
OCT
3 2.5 2 1.5 1 0.5
-1
0
1
2
NORMALIZED VOLUMETRIC STRESS [] PROBABILITY DENSITY []
PROBABILITY DENSITY []
NORMALIZED VOLUMETRIC STRESS []
SPH
3 2.5 2 1.5 1 0.5
0
0 -2
-1
0
1
2
NORMALIZED VOLUMETRIC STRESS []
-2
-1
0
1
2
NORMALIZED VOLUMETRIC STRESS []
ðpÞ
Fig. 4. Density distributions of the normalized volumetric stress rvol =ras in phase (p) predicted for composites reinforced by spherical (SPH) or randomly oriented octahedral (OCT), cube-shaped (CUB) or tetrahedral (TET) particles of volume fraction n(i) = 0.2 and conductivity contrast of ck = 10 that are subjected to an applied shear stress ras . Dashed lines denote matrix and solid lines particle behavior, thin vertical lines indicate corresponding phase averages. Results are based on all equivalent volume elements and all pertinent load cases.
pronounced edges and vertices ‘‘attract’’ stresses and fluxes, giving rise to higher magnitudes of the components of the inhomogeneity stress and flux concentration tensors. If the particles are the stiffer and more highly conductive phase, these changes to the concentration tensors, in turn, lead to higher effective elastic moduli and conductivities, respectively. It is of interest to study how these particle shape effects evolve with growing elasticity and conduction contrasts. For this purpose computations on the same multi-particle unit cells were also carried out for inhomogeneities the elasticity and conductivity tensors of which were modified such that both the elasticity contrast, cE, and the conduction contrast, ck, took values of 4 and 25. Fig. 7 presents such results in a plot of the normalized axial stress in matrix (lower set of curves) and particles (upper set of curves) as functions of the elasticity contrast in the range 1 6 cE 6 30. The phase averages obtained with the Mori–Tanaka method and the three-point estimates are plotted as solid and dotted lines, respectively. A close inspection of these results shows that the curves describing matrix and inhomogeneities do not meet at cE = 1, which is a consequence of the different Poisson numbers used for particles and matrix, compare Table 1. The phase averages and phase-level standard deviations obtained by periodic homogenization for elasticity contrasts of 4, 10 and 25 are represented by symbols, which indicate the particle shape, and by error bars showing the phase-level standard deviations, respectively. As stated in Section 2.3 these phase averages and phase-level standard deviations were evaluated via Eqs. (2) and (3), respectively, from the results for all pertinent load cases applied to each of the five equivalent volume elements pertaining to a given particle shape. In the interest of the clarity of the plot the data for spherical particles are shifted to the left from their proper position, whereas those for cube-shaped and tetrahedral particles are displaced to the right. The predictions shown in Fig. 7 confirm that at the selected inhomogeneity volume fraction of n(i) = 0.2 the microstresses in the matrix have only a rather weak dependence on the particle shape for all elasticity contrasts considered, with the phase averages being close to the analytical results throughout the range covered in the plot and the standard deviations showing only minor growth with increasing elasticity contrast. The microstresses in the particle phase, in contrast, are characterized by increasing differences between the phase averages predicted for the different particle shapes and by a marked growth of the standard deviations for the polyhedral inhomogeneities, the behavior being especially pronounced for tetrahedral particles. A very similar picture emerges for the thermal conduction behavior, which is depicted in Fig. 8. For this case the threepoint estimates for the phase averages of the fluxes were evaluated via an analogue of Eq. (4) that follows from expressions given, e.g., in (Böhm, Pahr, & Daxner, 2009). The main quantitative differences to the elastic case are a wider gap between the
30
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
2
TET
PROBABILITY DENSITY []
PROBABILITY DENSITY []
2
1.5
1
0.5
0
CUB
1.5
1
0.5
0 0
1
2
3
4
5
6
7
8
0
NORMALIZED EQUIVALENT STRESS []
2
3
4
5
6
7
8
NORMALIZED EQUIVALENT STRESS [] 2
OCT
PROBABILITY DENSITY []
PROBABILITY DENSITY []
2
1
1.5
1
0.5
0
SPH
1.5
1
0.5
0 0
1
2
3
4
5
6
7
NORMALIZED EQUIVALENT STRESS []
8
0
1
2
3
4
5
6
7
8
NORMALIZED EQUIVALENT STRESS []
ðpÞ
Fig. 5. Density distributions of the normalized equivalent stress reqv =ðEðmÞ aðmÞ DTÞ in phase (p) predicted for composites reinforced by spherical (SPH) or randomly oriented octahedral (OCT), cube-shaped (CUB) or tetrahedral (TET) particles of volume fraction n( i) = 0.2, elastic contrast cE = 10 and expansion contrast a(i)/a(m) = 0.1 that are subjected to a uniform temperature excursion DT. Dashed lines denote matrix and solid lines particle behavior, and thin vertical lines indicate corresponding phase averages. Results are based on all equivalent volume elements and all pertinent load cases.
phase averages of the flux in particles and matrix and a tendency towards somewhat lower standard deviations of the fields in the reinforcement phase. 3.3. Particle-level microfields In Fig. 9 the particle-level averages (‘‘particle averages’’) and particle-level standard deviations of the normalized axial stress component predicted for inhomogeneities of the four shapes considered are plotted for loading by a uniaxial macroscopic stress ran . These results were evaluated for each relevant load case and for each particle contained in the five equivalent volume elements pertaining to a given particle shape. The particle averages, which are indicated by the dark, s-shaped, line-like sequence of symbols in each plot, serve to bring out the inter-particle fluctuations of the axial stress component. Data from individual particles is ordered by increasing particle averages, the resulting ‘‘particle numbers’’ on the abscissa reaching values slightly exceeding 300 (5 equivalent volume elements 21 particles 3 pertinent load cases). The results obtained for the different particle shapes are qualitatively similar, the particle averaged axial stresses being somewhat higher for the composites reinforced by convex polyhedral particles than for the sphere-reinforced ones. A pronounced shape dependence is evident for the standard deviations of the axial stress component within individual particles, showing that tetrahedral particles are subject to intra-particle stress fluctuations that are markedly higher than the ones present in spherical inhomogeneities. Fairly similar behavior is again predicted for cube-shaped and octahedral particles, with the latter tending towards somewhat higher inter-particle and intra-particle fluctuations. A qualitatively similar behavior was obtained for the particle-level averages and standard deviations of the normalized flux components under loading by macroscopic axial heat fluxes, which are shown in Fig. 10. However, some tendency can be detected for the shape dependence of the intra-particle fluctuations to be less pronounced than in the elastic case. The relatively limited variation of the intra-particle fluctuations of the axial flux in cube-shaped particles, for which the predicted particle-level standard deviations appear much less ‘‘frizzy’’ than for octahedral and, especially, tetrahedral particles, is not understood at present. A quantitative evaluation of the intra-particle fluctuations of the fields in terms of the ranges of the particle averages appears possible in principle. However, a closer inspection of Figs. 9 and 10 shows that these ranges are determined mainly by the smallest and largest phase averages within the given sample of particles. These extremal values tend to result from special local configurations of the particle positions, i.e., from ‘‘rare events’’. This implies that considerably larger numbers of
31
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
PROBABILITY DENSITY []
2 1.5 1 0.5
CUB
2.5
PROBABILITY DENSITY []
TET
2.5
2 1.5 1 0.5
0
0 0
1
2
3
4
5
6
0
1
NORMALIZED AXIAL FLUX []
PROBABILITY DENSITY []
3
4
5
6
2 1.5 1 0.5
SPH
2.5
PROBABILITY DENSITY []
OCT
2.5
2
NORMALIZED AXIAL FLUX []
2 1.5 1 0.5
0
0 0
1
2
3
4
5
6
0
1
NORMALIZED AXIAL FLUX []
2
3
4
5
6
NORMALIZED AXIAL FLUX [] ðpÞ
NORMALIZED PHASE AVG. AXIAL STRESS []
Fig. 6. Density distributions of the normalized axial flux component qk =qan in phase (p) predicted for composites reinforced by spherical (SPH) or randomly oriented octahedral (OCT), cube-shaped (CUB) or tetrahedral (TET) particles of volume fraction n(i) = 0.2 and conductivity contrast of ck = 10 that are subjected to an applied flux qan . Dashed lines denote matrix and solid lines particle behavior, thin vertical lines indicate corresponding phase averages. Results are based on all equivalent volume elements and all pertinent load cases.
3
2.5
TET CUB OCT SPH
2
1.5
1
0.5
0 5
10
15
20
25
ELASTIC CONTRAST [] ðpÞ
Fig. 7. Phase averages and phase-level standard deviations of the normalized axial stress components rkk =ran predicted for composites reinforced by spherical (SPH, circles) or randomly oriented octahedral (OCT, pentagons), cube-shaped (CUB, triangles) or tetrahedral (TET, squares) particles of volume fraction n(i) = 0.2. Loading is by an applied uniaxial stress ran and elasticity contrasts of cE = 4, cE = 10 and cE = 25 are considered. Results obtained with the Mori–Tanaka method and with three-point estimates for identical spheres are shown as solid and dotted lines, respectively. The upper branch of results pertains to the particles, the lower to the matrix. Numerical results are based on all equivalent volume elements and all pertinent load cases.
particles must be sampled for obtaining reliable results than was possible within the present work. As to the standard deviations of the fields in the individual particles, the considered cases show fairly large ranges of values and a limited correlation with the particle averages for a given inhomogeneity shape. This allows the interpretation that the intra-particle fluctuations are not only determined by shape effects, but also by interactions with neighboring particles. These
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
NORMALIZED PHASE AVG. AXIAL FLUX []
32
3.5
TET CUB OCT SPH
3 2.5 2 1.5 1 0.5 0
5
10
15
20
25
30
CONDUCTIVITY CONTRAST [] ðpÞ
Fig. 8. Phase averages and phase-level standard deviations of the normalized axial flux components qk =qan predicted for composites reinforced by spherical (SPH, circles) or randomly oriented octahedral (OCT, pentagons), cube-shaped (CUB, triangles) or tetrahedral (TET, squares) particles of volume fraction n(i) = 0.2. Loading is by an applied flux qan and conductivity contrasts of ck = 4, ck = 10 and ck = 25 are considered. Results obtained with the Mori–Tanaka method and with three-point estimates for identical spheres are shown as solid and dotted lines, respectively. The upper branch of results pertains to the particles, the lower to the matrix. Numerical results are based on all equivalent volume elements and all pertinent load cases.
3.5
TET
3
NORMALIZED NORMAL STRESS []
NORMALIZED NORMAL STRESS []
3.5
2.5 2 1.5 1 0.5 0
CUB
3 2.5 2 1.5 1 0.5 0
0
50
100
150
200
250
300
0
50
PARTICLE NUMBER []
100
150
200
250
300
PARTICLE NUMBER []
NORMALIZED NORMAL STRESS []
NORMALIZED NORMAL STRESS []
3.5
OCT
3 2.5 2 1.5 1 0.5 0
SPH
3 2.5 2 1.5 1 0.5 0
0
50
100
150
200
PARTICLE NUMBER []
250
300
0
50
100
150
200
250
300
PARTICLE NUMBER [] ðpÞ
Fig. 9. Normalized averages (dark symbols) and standard deviations (error bars) of stress component rkk within individual randomly oriented tetrahedral (TET), cube-shaped (CUB), octahedral (OCT) or spherical (SPH) particles predicted for an inhomogeneity volume fraction of n(i) = 0.2. An elastic contrast cE = 10 and loading by an applied unit uniaxial stress ran are considered. Results are based on all equivalent volume elements and all pertinent load cases. Particles are numbered by increasing particle averages of axial stress.
considerations also indicate that the number of particles used in the present contribution is insufficient for an in-depth study of intra-particle fluctuations. In addition, particle-level standard deviations must be expected to be liable to degradation due to mesh resolution issues, as mentioned in Section 2.3.
33
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34 3.5
TET
NORMALIZED AXIAL FLUX []
NORMALIZED AXIAL FLUX []
3.5 3 2.5 2 1.5 1 0.5 0
2.5 2 1.5 1 0.5 0
0
3.5
50
100 150 200 PARTICLE NUMBER []
250
300
0
50
100
150
200
250
300
PARTICLE NUMBER [] 3.5
OCT
NORMALIZED AXIAL FLUX []
NORMALIZED AXIAL FLUX []
CUB
3
3 2.5 2 1.5 1 0.5 0
SPH
3 2.5 2 1.5 1 0.5 0
0
50
100
150
200
250
PARTICLE NUMBER []
300
0
50
100
150
200
250
300
PARTICLE NUMBER [] ðpÞ
Fig. 10. Normalized averages (dark symbols) and standard deviations (error bars) of flux component qk within individual randomly oriented tetrahedral (TET), cube-shaped (CUB), octahedral (OCT) or spherical (SPH) particles predicted for an inhomogeneity volume fraction of n( i) = 0.2. A conductivity contrast of ck = 10 and loading by an applied unit flux qan are considered. Results are based on all equivalent volume elements and all pertinent load cases. Particles are numbered by increasing particle averages of axial flux.
The present results can be expected to hold for moderate particle volume fractions, qualitatively similar responses having been obtained for polydisperse particles of a volume fraction of n(p) = 0.3 (Rasool, 2011). For carrying out analogous studies at elevated particle volume fractions, however, advances in synthesizing suitable phase arrangements are required. 4. Conclusions A detailed study was carried out on particle shape effects on the macroscopic and microscopic elastic, thermoelastic and thermal conduction behavior of macroscopically isotropic, particle reinforced composites. For the cases considered, viz., a compliant matrix reinforced by stiffer inhomogeneities having the shapes of spheres or regular, randomly oriented octahedra, cubes and tetrahedra of a volume fraction of 20%, some clear trends could be identified. On the one hand, randomly oriented convex polyhedra lead to some increase in the effective moduli compared to spherical particles. On the other hand, polyhedral inhomogeneities give rise to much higher fluctuations of the microstresses and and microfluxes in the particles than do spherical ones, and thus to markedly higher probabilities of encountering high values of stress components and stress measures. The microfields in the matrix, in contrast, do not show a strong dependence on the particle shapes. The prospects for making use of particle shape effects in practical applications of particle reinforced composites appear to be limited. For the particle volume fraction of 20% considered here, gains of some 10% in terms of the effective thermal conductivity are offset by much higher local maximum stress levels – which imply a markedly higher vulnerability to microdamage – in tetrahedral particles when compared to spherical inhomogeneities. When mechanical strength is taken into account, the most desirable particle shapes appear to be spheres. Acknowledgement The financial support of A.R. by the Higher Education Commission of Pakistan is gratefully acknowledged. References Asahina, D., & Bolander, J. (2011]). Powder Technology, 213, 92–99. Benveniste, Y. (1987]). Mechanics of Materials, 6, 147–157. Berryman, J., & Milton, G. (1988]). Journal of Physics D, 21, 87–94.
34
A. Rasool, H.J. Böhm / International Journal of Engineering Science 58 (2012) 21–34
Böhm, H., Pahr, D., & Daxner, T. (2009]). In V. Silberschmidt (Ed.), Computational and experimental mechanics of advanced materials. CISM courses and lectures (Vol. 514, pp. 167–223). Vienna: Springer–Verlag. Bornert, M., Hervé, E., Stolz, C., & Zaoui, A. (1994]). Applied Mechanics Reviews, 47, 66–76. Buryachenko, V. (1996]). Acta Mechanica, 119, 93–117. Chawla, N., Sidhu, R., & Ganesh, V. (2006]). Acta Materialia, 54, 1541–1548. Drugan, W., & Willis, J. (1996]). Journal of the Mechanics and Physics of Solids, 44, 497–524. Duschlbauer, D., Böhm, H., & Pettermann, H. (2006]). Journal of Composite Materials, 40, 2217–2234. Eroshkin, O., & Tsukrov, I. (2005]). International Journal of Solid Structures, 42, 409–427. Eshelby, J. (1957]). Proceedings of the Royal Society of London, A241, 376–396. Finot, M., Shen, Y., Needleman, A., & Suresh, S. (1994]). Metallurgical Transactions, 25A, 2403–2420. Fritzen, F., & Böhlke, T. (2011]). International Journal of Solid Structures, 48, 706–718. Galli, M., Botsis, J., & Janczak-Rusch, J. (2008]). Computational Materials Science, 41, 312–321. Gao, X., & Liu, M. (2012]). Journal of the Mechanics and Physics of Solids, 60, 261–276. Ghahremani, F., Hutchinson, J., & Tvergaard, V. (1990]). Journal of the American Ceramic Society, 73, 1548–1554. Glushkov, E., Glushkova, N., & Lapina, O. (1998]). International Journal of Solid Structures, 36, 1105–1128. Han, W., Eckschlager, A., & Böhm, H. (2001]). In Y. Zhang (Ed.). Proceedings of the 13th International Conference on Composite Materials (ICCM-13). Beijing: Scientific and Technical Documents Publishing House, paper #1222. Hashemi, R., Avazmohammadi, R., Shodja, H., & Weng, G. (2009]). Philosophical Magazine Letters, 89, 439–451. Hashin, Z., & Shtrikman, S. (1962]). Journal of Applied Physics, 33, 3125–3131. Hashin, Z., & Shtrikman, S. (1963]). Journal of the Mechanics and Physics of Solids, 11, 127–140. He, Q., & Curnier, A. (1995]). International Journal of Solid Structures, 32, 1433–1457. Hill, R. (1963]). Journal of the Mechanics and Physics of Solids, 11, 357–372. Hu, L., & Chen, S. (2011]). Advances in Material Research, 194–196, 972–976. Huet, C. (1990]). Journal of the Mechanics and Physics of Solids, 38, 813–841. Kachanov, M., & Sevostianov, I. (2005]). International Journal of Solid Structures, 42, 309–336. Kachanov, M., Tsukrov, I., & Shafiro, B. (1994]). Applied Mechanics Reviews, 47, S151–S174. Kang, H., & Milton, G. (2008]). Archive for Rational Mechanics and Analysis, 188, 93–116. Kanit, T., Forest, S., Galliet, I., Mounoury, V., & Jeulin, D. (2003]). International Journal of Solid Structures, 40, 3647–3679. Kenesei, P., Biermann, H., & Borbély, A. (2006]). Advances in Engineering Materials, 8, 500–506. Lippmann, N., Steinkopff, T., Schmauder, S., Gumbsch, P. (1996). In M. Rappaz, M. Kedro (Eds.), Proceedings of the general COST 512 workshop on modelling in materials science and processing, European Commission, D.G. XII (pp. 404–410). Brussels, Belgium. Meijer, G., Xia, Z., & Ellyin, F. (1997]). Acta Materialia, 45, 3237–3249. Michel, J., Moulinec, H., & Suquet, P. (1999]). Computational Methods in Applied Mechanics and Engineering, 172, 109–143. Mishnaevsky, L., Dong, M., Hönle, S., & Schmauder, S. (1999]). Computational Materials Science, 16, 133–143. Nogales, S., & Böhm, H. (2008]). International Journal of Engineering Science, 46, 606–619. Nozaki, H., & Taya, M. (2001]). Journal of Applied Mechanics, 68, 441–452. Onaka, S. (2001]). Philosophical Magazine Letters, 81, 265–272. Rasool, A. (2011). Improved multi-particle unit cell models for studying particle reinforced composites, Ph.D. thesis, Vienna University of Technology, Vienna, Austria. Rodin, G. (1996]). Journal of the Mechanics and Physics of Solids, 44, 1977–1995. Schneider, P., & Eberly, D. (2003]). Geometric tools for computer graphics. San Francisco: Morgan-Kaufmann. Sevostianov, I., Kachanov, M., & Zohdi, T. (2008]). International Journal of Solid Structures, 45, 4375–4383. Shen, H., & Lissenden, C. (2003]). Metallurgical and Materials Transactions, 36A, 1653–1660. Smit, R., Brekelmans, W., & Meijer, H. (1998]). Computational Methods in Applied Mechanics and Engineering, 155, 181–192. Terada, K., & Kikuchi, N. (1996]). Material Science Research Institute, 2, 65–72. Torquato, S. (1991]). Applied Mechanics Reviews, 44, 37–75. Torquato, S. (1998]). Journal of the Mechanics and Physics of Solids, 46, 1411–1440. Torquato, S. (2002]). Random heterogeneous media. New York, NY: Springer-Verlag. Weissenbek, E., Böhm, H., & Rammerstorfer, F. (1994]). Computational Materials Science, 3, 263–278. Widom, B. (1966]). Journal of Chemical Physics, 44, 3388–3394. Zheng, Q., Zhao, Z., & Du, D. (2006]). Journal of the Mechanics and Physics of Solids, 54, 368–383. Zhu, L., Cai, W., & Tu, S. (2009]). Journal of Computation, 4, 1237–1242.