Three-dimensional reconstruction and homogenization of heterogeneous materials using statistical correlation functions and FEM

Three-dimensional reconstruction and homogenization of heterogeneous materials using statistical correlation functions and FEM

Computational Materials Science 51 (2012) 372–379 Contents lists available at SciVerse ScienceDirect Computational Materials Science journal homepag...

3MB Sizes 0 Downloads 20 Views

Computational Materials Science 51 (2012) 372–379

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Three-dimensional reconstruction and homogenization of heterogeneous materials using statistical correlation functions and FEM M. Baniassadi a,b,⇑, B. Mortazavi a,b, H. Amani Hamedani c, H. Garmestani c, S. Ahzi a,c, M. Fathi-Torbaghan d, D. Ruch b, M. Khaleel e a

University of Strasbourg, IMFS, 2 Rue Boussingault, 67000 Strasbourg, France Centre de Recherche Public Henri Tudor, AMS, 70 Rue de Luxembourg, L-4221 Esch-sur-Alzette, Luxembourg School of Materials Science and Engineering, Georgia Institute of Technology, 771 Ferst Dr. N.W. Atlanta, GA 30332-0245, USA d Faculty of Electrical and Computer Science, University of Siegen, Holderlinstr. 3, 57068 Siegen, Germany e Pacific Northwest National Laboratory, 790 6th Street, Richland, WA, USA b c

a r t i c l e

i n f o

Article history: Received 2 July 2011 Received in revised form 28 July 2011 Accepted 1 August 2011 Available online 31 August 2011 Keywords: Heterogeneous media Three-dimensional microstructure reconstruction Two-point correlation functions Two-point cluster functions FEM analysis

a b s t r a c t In this study, a previously developed reconstruction methodology is extended to three-dimensional reconstruction of a three-phase microstructure, based on two-point correlation functions and two-point cluster functions. The reconstruction process has been implemented based on hybrid stochastic methodology for simulating the virtual microstructure. While different phases of the heterogeneous medium are represented by different cells, growth of these cells is controlled by optimizing parameters such as rotation, shrinkage, translation, distribution and growth rates of the cells. Based on the reconstructed microstructure, finite element method (FEM) was used to compute the effective elastic modulus and effective thermal conductivity. A statistical approach, based on two-point correlation functions, was also used to directly estimate the effective properties of the developed microstructures. Good agreement between the predicted results from FEM analysis and statistical methods was found confirming the efficiency of the statistical methods for prediction of thermo-mechanical properties of three-phase composites. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction Development of advanced microstructure reconstruction methodologies is essential to access variety of analytical information associated with complexities in the microstructure of multi-phase materials. Several experimental and theoretical techniques such as X-ray computed tomography (CT) scanning and computer generated micrographs have been used to obtain a sequence of twodimensional (2D) images that can be further reconstructed in a 3D space. However, due to cost of sample preparation processes, simulation methods are often more applicable in reconstruction of heterogeneous microstructures in different areas [1–7]. Using lower-order statistical correlation functions, Torquato [8] established the reconstruction of one- and two-dimensional microstructures with short-range order using stochastic optimization. However, he later showed that the lower-order correlation functions cannot solely represent a two-phase heterogeneous material and therefore more than one solution may exist for a specific loworder correlation function [8]. Sheehan and Torquato [9] later introduced more orientations in the correlation function to effectively eliminate the effect of artificial anisotropy. In the case of ⇑ Corresponding author at: Université de Strasbourg, IMFS, 2 Rue Boussingault, 67000 Strasbourg, France. Tel.: +33 2545580628; fax: +33 2425991555. E-mail address: [email protected] (M. Baniassadi). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.08.001

multi-phase materials, Kröner [10,11] and Beran [12] have developed a model in which n-point statistical functions are used to account for the volume fraction and the number of phases in the microstructure. Using higher-order correlation functions, one can account for the contribution of shape and geometry effects [8]. Torquato and Yeong [13] also developed a new hybrid stochastic reconstruction technique for reconstruction of three-dimensional (3D) random media by using the information from the lineal path function and the two-point correlation functions during the nucleation annealing process. Different optimization techniques such as simulated annealing and maximum entropy have been applied in order to improve the simulation process [14]. Torquato and coworkers have used two-point cluster function to capture the salient structural features of the media for verification of reconstruction procedure [15]. For example, other researchers have shown the necessity of employing the two-point probability, two-point cluster, and lineal path functions to produce a generally accurate statistical reconstruction [16]. In addition to 3D reconstruction processes based on probability functions, these functions can be used to account for more details of microstructure heterogeneities and for the relationships between microstructure, local and effective properties of a multi-phase material. The effective properties can be obtained via perturbation expansions [17,18]. One general approach for prediction of the effective properties of a two-phase material with properties of each

M. Baniassadi et al. / Computational Materials Science 51 (2012) 372–379

373

Fig. 1. Basic steps in the realization algorithm (OF = objective function; MC = Monte Carlo).

phase near the average ones is called ‘‘weak-contrast’’ expansion. However, in materials with a high degree of contrast between the properties of their phases, ‘‘strong-contrast’’ theory is applied. Brown [19] suggested an expansion for effective dielectric property of two-phase heterogeneous materials. This expansion for perturbation homogenization was modified and extended for elasticity by Torquato [20] for two-phase materials and later the solution was extended to multi-phase materials by others [18,21]. Several numerical methods can be used to obtain the effective thermal/ electrical conductivity as well as effective elastic properties of multi-phase composites of complex geometries containing arbitrary oriented inhomogeneities [22–24]. The herein proposed work extended a previously developed reconstruction methodology to 3D microstructure reconstruction based on two-point correlation function and two-point cluster function [25,26]. Using a hybrid stochastic methodology for simulating the virtual microstructure, growth of the phases represented by different cells is controlled by optimizing parameters such as rotation, shrinkage, translation, distribution and growth rates of the cells [26]. We used the finite element method (FEM) to predict the effective thermo-mechanical properties such as the elastic modulus and thermal conductivity of the reconstructed microstructure. We also used the strong contrast statistical method,

based on two-point correlation functions calculated for the initial computer generated microstructure, to estimate the effective elastic and thermal properties. Comparison of the results from both approaches show a good agreement. 2. Statistical functions used in microstructure reconstruction 2.1. Two-point correlation function The N-point correlation functions for phase (a) can be interpreted as the probability that N points at positions x1, x2,. . .xN are found in phase (a). For statistically homogeneous media, N-point correlation functions are a translational invariant and therefore depend only on relative positions of the N points. In this study, we used two-point correlation (N = 2) as statistical microstructure descriptor of heterogeneous materials. For a three-phase composites, the indices (i, j) in the probability take the values of 1, 2, 3. Therefore, we have nine probability functions Pij (i, j = 1.3). Due to normality conditions, the following equations are satisfied: 3 X 3 X i¼1

j¼1

Pij ðr Þ ¼ 1

ð1Þ

374

M. Baniassadi et al. / Computational Materials Science 51 (2012) 372–379

correlation functions [26]. The realization process includes three steps: (1) generation, (2) distribution, and (3) growth of cells. Here, cells (or alternately grains or particles) refer to initial geometries assigned to each phase before the growth step. During the initial microstructure generation, basic cells are created from the random nucleation points and then the growth occurs as in crystalline grain growth in real materials [26]. After distribution of nucleation points and assignment of basic cell geometries, the growth of cells starts according to the cellular automaton approach. The three steps of realization algorithm are repeated continuously to satisfy the optimization parameters until an adequately realistic microstructure is developed as compared statistically to the true microstructure. It is worth noting that in various steps of algorithm execution; several controlling parameters are developed that facilitate minimization of the objective function (OF) which is an index of successful realization. This objective function is defined based on the three independent two-point correlation functions (Pii) and two-point cluster functions (P cii ) as follows: Fig. 2. Representation of vectors in spherical coordinate. 3 X

Pij ðr Þ ¼ v i

ðFor i ¼ 1; 2 and 3Þ

OF ¼ ððP12 Þreal  ðP12 Þsim Þ2 þ ð2Þ

j¼1

Knowing that the probability functions are symmetric (Pij = Pji) reduces the above conditions to only three independent relations. This results in the important conclusion that only three of the nine probabilities should be considered as independent variables. For instance, we can choose P11, P12, and P22 as the three independent probability parameters. 2.2. Two-point cluster function Torquato [8] have introduced the two-point cluster function (TPCCF(x, x0 )) as the probability of finding both points x and x’ in the same cluster of one of the phases. This quantity is a useful signature of the microstructure insofar as it reflects clustering information. Incorporation of such information rather than the lowerorder two-point cluster functions have led to the formulation of rigorous bounds on transport and mechanical properties of twophase media [8]. 3. Development of a Monte Carlo reconstruction methodology A new algorithm is presented based on Monte Carlo methodology for the reconstruction of microstructures using two-point

2 X ððPii Þreal  ðPii Þsim Þ2 i¼1

3 X þ ððPcii Þreal  ðP cii Þsim Þ2

ð3Þ

i¼1

where the subscripts real and sim indicate, respectively, the values from the real and reconstructed microstructures. The procedure of reconstruction and optimization is repeated until the objective function takes a value that is of the same order as the Monte Carlo (M-C) repeat error. In this paper, material’s heterogeneity is represented by statistical two-point correlation functions and two-point cluster functions. Hypothetical statistical functions are optimized and compared to the initial statistical distribution functions of the sample microstructure. Stochastic optimization methodologies incorporate probabilistic (random) elements, either in the input data (the objective function, the constraints, etc.), or in the algorithm itself (through random parameters, etc.) or both [26]. By applying different optimization parameters to the simulations, a minimum error is achieved through minimization of the objective function that is constructed from the comparison of the two-point correlation function and two-point cluster function of the experimental (sample) and simulated (realization) microstructures. A direct simple search optimization technique was used for finding the minimum objective function [26]. Fig. 1 depicts a schematic of the reconstruction algorithm.

Fig. 3. Finite element illustration and boundary condition of computer generated and reconstructed microstructure (left: computer generated and right: reconstructed microstructures), for thermal and mechanical loading.

M. Baniassadi et al. / Computational Materials Science 51 (2012) 372–379

375

Fig. 4. 2D arbitrary sections in the z-direction of the 3D reconstructed microstructure. (a) layer close to the bottom surface, (b) layer in the middle area, and (c) layer close to the top surface.

Fig. 5. (a) Two-point correlation function (P11) for the red-phase, (b) two-point correlation function (P22) for the black-phase and (c) two-point correlation function (P12) for the black–red phases for the computer generated and reconstructed microstructures. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

376

M. Baniassadi et al. / Computational Materials Science 51 (2012) 372–379

Fig. 6. Two-point cluster function P c11 for the red-phase. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. (a) An arbitrary 2D section of the 3D reconstructed microstructure (black = porosity); (b–d) the corresponding percolation of voids (porosity) showing the percolation clusters by similar colors other than white. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4. Statistical characterization of microstructures

u2a

4.1. Mechanical model Exact perturbation expansions were used to predict the effective stiffness tensor of a macroscopically isotropic two phase composite media [20]. In this method, an integral equation for the local ‘‘cavity’’ strain field, relate to the local ‘‘polarization’’ stress field. Manipulating integral equations for the local ‘‘cavity’’ strain field and polarization homogenization leads to find series’ expansions for the effective stiffness tensor in powers of the elastic polarizabilities. The general term of the expansion for macroscopically isotropic media can be simplified as follows [20]:

" # 1 X l kaq Kh þ aq Ks ¼ ua I  BðaÞ n keq leq n¼2

ð4Þ

For simplification, the calculations have been performed using only the first and second order terms of the n-point tensors BðaÞ n for phase a [20].

u2a

" # l kaq ðaÞ Kh þ aq Ks ¼ ua I  BðaÞ 2  B3 keq leq

ð5Þ

The subscript/superscript q stands for the reference phase and the subscript/superscript a stands for the other phase and the subscript/superscript e stands for effective properties. where /a is the volume fraction of phase a and I is the fourthorder identity tensor, Kh and Ks are the fourth-order hydrostatic

377

M. Baniassadi et al. / Computational Materials Science 51 (2012) 372–379

Fig. 8. Temperatures and Mises stress contours (left: computer generated and right: reconstructed microstructure).

projection tensor and shear projection tensor and kaq and laq are the bulk and shear modulus polarizability, respectively. In Eqs. (6) and (7) the n-points tensor coefficients (Bn) are given by the following integrals over products of the Uq tensors in Eq. (8), for the reference phase q and the N-point probability functions P an for phase a. The tensor Uq is calculated based on the position-dependent fourth-order tensor H(r) in Eq. (9) and related tensor Lq [20] for phase q: ðaÞ

B2 ¼ ðaÞ

B3 ¼ ðqÞ U ijkl ðrÞ

Z

0

e



ðaÞ

dx U ðqÞ ðx; x0 Þ½P2 ðx; x0 Þ  u2a  1

ua

Z

0

dx

Z

ðqÞ ðqÞ Lijmn Hmnkl ðrÞ

00

ðaÞ

dx U ðqÞ ðx; x0 Þ : U ðqÞ ðx0 ; x00 ÞD3 ðx; x0 ; x00 Þ 

¼ 3K q þ 4Gq 5Gq þ laq HðqÞ ðrÞ ijkl 3ðK q þ 2Gq Þ

¼

ð6Þ

 

Sample

Sample 1 (Thermal conductivity and elastic modules)

Sample 2 (Thermal conductivity and elastic modules)

Sample 3 (Thermal conductivity and elastic modules)

Phase 1 Phase 2 Phase 3

1 10 10

1 1 10

1 10 1

ð7Þ

 5Gq d K aq  laq ij HðqÞ ðrÞ 3ðK q þ 2Gq Þ 3 mmkl ð8Þ

Here, H(r) is the symmetric double gradient tensor which is given as follows [8]:

1 1 ðqÞ Hijkl ðrÞ ¼ 2X½3K q þ 4Gq  r 3     3ð3  aq Þ aq dij dkl  d½dik dil þ dil djk  daq dij tk tl þ dkl ni nj þ 2    dik t j t l þ dil t j t k þ dik t i nl þ dij t i t k þ 15aq ti t j t k t l

Table 1 The thermal conductivity and the elastic modulus associated with each phase for the cases considered in this study.

The Kq and Gq are bulk and shear modulus of the reference phase (phase q) and The Ka and Ga are bulk and shear modulus of the other phase (phase a) [8]. kaq and laq are the bulk and shear, modulus polarizability respectively which are given as follows:

kaq ¼

laq ¼ ð9Þ

Ka  Kq K a þ 43 Gq

Ga þ

ð10Þ

Ga  Gq Gq ½3K q =2þ4Gq =d

ð11Þ

K q þ2Gq

And:

aq ¼ 3kq =Gq þ 1

ð12Þ

378

M. Baniassadi et al. / Computational Materials Science 51 (2012) 372–379

Fig. 9. Elastic module of reconstructed microstructure using FEM and strong contrast technique (left), thermal conductivity of reconstructed microstructure using FEM and strong contrast technique (right).

ðaÞ

P ðx; x0 Þ

2 0 00 DðaÞ ðx; x ; x Þ ¼

ðaÞ 3

P ðx; x0 ; x00 Þ

ðaÞ

P1 ðx0 Þ ðaÞ P2 ðx0 ; x00 Þ

3







ðaÞ

ð13Þ

The following approximation is used for estimation of the threepoint correlation function from the two-point function [27–30]:



jxx0 j Pa ðx; x00 Þ 0 jxx j þ jxx00 j 2  a 0 00 jxx00 j P2 ðx ; x Þ a 0 þ 0 ðx; x Þ P jxx j þ jxx00 j 2 Pa1 ðxÞ

Pa3 ðx; x0 ; x00 Þ ffi

ð14Þ

Fig. 2 defines the variables used in this approximation in local coordinates. 4.2. Thermal conductivity The effective conductivity tensor ke of the composite materials can be determined using the strong-contrast formulation of the statistical continuum theory [8,31] considering the isotropic properties of the phases:

1 I  3kq baq Pa1 ðxÞ  Z  a P2 ðx; x0 Þ  Pa1 ðxÞPa1 ðx0 Þ q 0 2  M ðx; x0 Þdx  k2q d baq a a 0 P 1 ðxÞP1 ðx Þ  Z Z  a P3 ðx; x0 ; x00 Þ P a2 ðx; x0 ÞPa2 ðx; x00 Þ q 0  a a 0  a a 0 a 00 M ðx; x Þ P1 ðxÞP1 ðx Þ P1 ðxÞP1 ðx ÞP1 ðx Þ

fke  kq Ig1  fke þ 2kq Ig ¼

0

00

 Mq ðx0 ; x00 Þdx dx    

ð15Þ

Here, I is the second-order identity tensor and baq is the polarizability given by:

baq ¼

ka  kq ka þ 2kq

ð16Þ

with kq representing the reference conductivity. The subscript q stands for the reference phase and the subscript a stands for the other phase. The second-order tensor Mq(x, x0 ) is defined by:

M q ðx; x0 Þ ¼

1

3tt  I

ð17Þ

Xkq ðjx  x0 jÞ3

where X is the total solid angle contained in a 3-dimensional sphere 0Þ and t ¼ ðxx . jxx0j ðaÞ

ðaÞ

ðaÞ

P1 ðxÞ, P 2 ðx; x0 Þ and P 3 ðx; x0 ; x00 Þ are the probability functions that contain microstructure information. The one-point probability ðaÞ function, P 1 ðxÞ is the volume fraction of phase a. The two-point ðaÞ correlation function, P 2 ðx; x0 Þ, is calculated from the Monte Carlo

simulation. The three-point probability function, P3 ðx; x0 ; x00 Þ, is calculated from the analytical approximation given by Eq. (14) [27]. 5. FEM characterization of multi-phase heterogeneous materials The computer generated sample and the 3D reconstructed microstructure based on two-point correlation functions and two-point cluster functions are used for our FEM characterization. Finite Element simulations were carried out using ABAQUS/Standard (Version 6.10). Due to the extensive computationally time, only ten layers of the real specimens were included in the modeling. For the purpose of thermal modeling, the specimen was meshed using eight-node linear heat transfer brick (DC3D8-type) elements. For the mechanical modeling, the eight-node linear brick, 3D stress with reduced integration (C3D8R-type) elements, were used. Each mesh element was assigned to the corresponding phase. In order to obtain the thermal conductivity of the specimen, constant heating surface heat flux was applied to a plane in the X direction while cooling surface heat flux equal to the cooling heating flux was applied in the opposite surface. In this way, steady state heat transfer criteria will be fully observed and by averaging the temperatures in each surface, the created temperature gradient as a function of distance in the specimen can be evaluated. The loading condition for the thermal and mechanical models are illustrated in Fig. 3. Using a one-dimensional form of the Fourier law, the thermal conductivity of the specimens was obtained. In order to obtain the elastic modulus of the specimen, a small strain was applied to the loading surface in its normal direction while the opposite surface was fixed only in its normal direction. By summation of the reaction forces in the fixed surface, the applied stress was calculated. Then, the elastic modulus of the specimen was obtained using Hook’s law. 6. Result and discussion The computer-generated three-phase sample is assumed to contain 10% volume fraction of red phase, 30% of the green phase and 60% of the black phase. This computer-generated three-phase sample is reconstructed and imported to the ABAQUS package for the FEM characterization. Fig. 4 shows 2D sections of three arbitrary layers taken from through-the-depth of the corresponding reconstructed microstructure in the ABAQUS package. These

M. Baniassadi et al. / Computational Materials Science 51 (2012) 372–379

sections are arbitrarily chosen from the top, middle and bottom parts of the 3D reconstructed domain. The corresponding two-point correlation functions (P11 for red– red and P22 for black-black and P12 for red–black) are calculated for both computer-generated and reconstructed microstructures shown in Fig. 5. As shown in this figure, there is a good agreement between the two-point correlation functions of the reconstructed and computer-generated microstructures. The reconstruction process is performed based on the two-point correlation and the twopoint cluster functions which had been extracted from computer generated microstructure. To check the validity of the reconstruction process, two-point cluster function for non-percolated phase (red–red) is calculated and shown in Fig. 6. Good agreement between the calculated two-point cluster functions for the two microstructures is obtained which strongly confirms the validity of the reconstruction process. The boundaries of the percolated regions of different phases are identified for one of the phases (red-phase) in 2D section shown in (Fig. 7a) in which the phase percolation is less than the percolation threshold. The percolated aggregates have been recognized using different colors in Fig. 7b. In Fig. 7b and c, wide percolated clusters have been observed in the cut section images. As the other two phases are intrinsically percolated and their corresponding twopoint correlation functions and two-point cluster functions are identical, there was no need to analyze the percolation in these phases. In Fig. 8, the temperatures and Von Mises stress contours are represented for different cases for the local properties. The normalized properties (both thermal and elastic) for the three phases are taken (1, 10, 10) for case 1, (1, 1, 10) for case 2 and (1, 10, 1) for case 3. Table 1 summarizes values of the phase properties assumed for different target samples. As it can be clearly seen, there is fine agreement between these three cases with respect of the obtained fields of temperature and stress values. As a result, the differences in the obtained thermal conductivity and elastic modulus of the two microstructures are less than 1% error. Elastic properties and thermal conductivity of these microstructures have been compared using strong contrast and FEM analysis of 3D reconstructed microstructures (cases). The elastic properties for the three target samples have been shown in Fig. 9 (left). The FEM results show very good agreement with strong contrast results which were obtained using statistical correlation functions of the microstructure. Similarly, thermal conductivity for the samples has been calculated using FEM analysis of 3D reconstructed microstructure and strong contrast technique. A small gap has been observed between the results obtained from the FEM and strong contrast method (Fig. 9 (right)). This shows the validity of the proposed statistical homogenization technique for three-phase heterogeneous materials. 7. Conclusion A Monte Carlo methodology is developed to reconstruct 3D microstructures of a three-phase microstructure. Two-point correlation functions and two-point cluster functions are used as microstructure descriptors in the reconstruction procedure. Using a hybrid stochastic reconstruction technique, optimization of the two-point correlation functions during different 3D realizations is performed repeatedly. The main challenge in the 3D reconstruction is incorporating two-point cluster function as complimentary

379

statistical descriptor to perform reconstruction technique. Comparison of the two-point correlation functions from different sections of the final 3D reconstructed microstructure with the initial computer generated microstructure shows good agreement. In addition, we have shown that the thermo-mechanical properties of the generated and reconstructed microstructures are close by means of FEM simulations. This supports the capability of our proposed methodology to reconstruct 3D microstructure. We have also used the statistical homogenization technique to compute the effective elastic and thermal properties. The comparison of the results with those of the FEM simulations shows a good agreement. This agreement between the two approaches suggest that the statistical approach is a reliable approach. Acknowledgment The authors would like to acknowledge the financial support from the Fond National de la Recherche (FNR) of the Grand Duché of Luxembourg (AFR Research Assistantship for Majid Baniassadi) and (AFR COTCH Project Grants for Bohayra Mortazavi). References [1] B. Bochenek, R. Pyrz, Computational Materials Science 31 (2004) 93. [2] S.-Y. Chung, T.-S. Han, Computational Materials Science 49 (2010) 705. [3] Z.R. Liang, C.P. Fernandes, F.S. Magnani, P.C. Philippi, Journal of Petroleum Science and Engineering 21 (1998) 273. [4] A. Pierret, Y. Capowiez, L. Belzunces, C.J. Moran, Geoderma 106 (2002) 247. [5] V. Sundararaghavan, N. Zabaras, Computational Materials Science 32 (2005) 223. [6] M.S. Talukdar, O. Torsaeter, Journal of Petroleum Science and Engineering 33 (2002) 265. [7] D.S. Li, M. Baniassadi, H. Garmestani, S. Ahzi, M. Taha, D. Ruch, Journal of Computational and Theoretical Nanoscience 7 (2010) 1462. [8] S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, New York, 2002. [9] N. Sheehan, S. Torquato, Journal of Applied Physics 89 (2001) 53. [10] E. Kröner, Statistical Continuum Mechanics, Springer-Verlag, Wien, 1977. [11] E. Kröner, Journal of the Mechanics and Physics of Solids 25 (1977) 137. [12] M.J. Beran, Statistical Continuum Theories, Interscience, Publishers, New York, 1968. [13] C.L.Y. Yeong, S. Torquato, Physical Review E 57 (1998) 495. [14] C. Manwart, R. Hilfer, Physical Review E 59 (1999) 5596. [15] Y. Jiao, F.H. Stillinger, S. Torquato, Physical Review E 77 (2008) 031135. [16] M.A. Davis, S.D.C. Walsh, M.O. Saar, Physical Review E 83 (2011) 026706. [17] D.T. Fullwood, B.L. Adams, S.R. Kalidindi, Journal of the Mechanics and Physics of Solids 56 (2008) 2287. [18] A. Tewari, A.M. Gokhale, J.E. Spowart, D.B. Miracle, Acta Materialia 52 (2004) 307. [19] J.W.F. Brown, The Journal of Chemical Physics 23 (1955) 1514. [20] S. Torquato, Journal of the Mechanics and Physics of Solids 45 (1997) 1421. [21] A. Mikdam, A. Makradi, S. Ahzi, H. Garmestani, D.S. Li, Y. Remond, Journal of the Mechanics and Physics of Solids 57 (2009) 76. [22] A. Giraud, C. Gruescu, D.P. Do, F. Homand, D. Kondo, International Journal of Solids and Structures 44 (2007) 2627. [23] P.D. Spanos, A. Kontsos, Probabilistic Engineering Mechanics 23 (2008) 456. [24] M. Wang, N. Pan, Journal of Computational Physics 228 (2009) 5978. [25] H. Garmestani, M. Baniassadi, D.S. Li, M. Fathi, S. Ahzi, International Journal of Theoretical and Applied Multiscale Mechanics 1 (2009) 134. [26] M. Baniassadi, H. Garmestani, D.S. Li, S. Ahzi, M. Khaleel, X. Sun, Acta Materialia 59 (2011) 30. [27] A. Mikdam, A. Makradi, S. Ahzi, H. Garmestani, D.S. Li, Y. Remond, International Journal of Solids and Structures 46 (2009) 3782. [28] M. Baniassadi, A. Laachachi, A. Makradi, S. Belouettar, D. Ruch, R. Muller, H. Garmestani, V. Toniazzo, S. Ahzi, Thermochimica Acta 520 (2011) 33. [29] H.A. Hamedani, M. Baniassadi, M. Khaleel, X. Sun, S. Ahzi, D. Ruch, H. Garmestani, Journal of Power Sources 196 (2011) 6325. [30] M. Baniassadi, F. Addiego, A. Laachachi, S. Ahzi, H. Garmestani, F. Hassouna, A. Makradi, V. Toniazzo, D. Ruch, Acta Materialia 59 (2011) 2748. [31] D.C. Pham, S. Torquato, Journal of Applied Physics 94 (2003) 6591.