INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2010; 84:972–988 Published online 17 May 2010 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/nme.2928 The statistical second-order two-scale method for mechanical properties of statistically inhomogeneous materials Fei Han1, ∗, † , Junzhi Cui1,2 and Yan Yu1 1 Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710129, People’s Republic of China 2 LSEC, ICMSEC, Academy of Mathematics and System Sciences, Chinese Academy of Sciences, Beijing 100190, People’s Republic of China SUMMARY A class of random composite materials with statistically inhomogeneous microstructure, for example, functionally graded materials is considered in this paper. The microstructures inside a component are gradually varying in the statistical sense. In view of this particularity, a novel statistical second-order two-scale (SSOTS) method is presented to predict the mechanical properties, including stiffness, and elastic limit. To develop this method, the microstructures of statistically homogeneous, and inhomogeneous materials are represented. In addition the SSOTS formulas are derived based on normalized cell depending on the position variables by a constructing way, and the algorithm procedure is described. The mechanical properties of the different inhomogeneous materials are evaluated. The numerical results are compared with the experimental findings. It shows that the SSTOS method is effective. Copyright 䉷 2010 John Wiley & Sons, Ltd. Received 31 March 2009; Revised 26 February 2010; Accepted 24 March 2010 KEY WORDS: multi-scale; statistically inhomogeneous materials; mechanical properties; homogenization 1. INTRODUCTION The statistically inhomogeneous materials are a kind of random distribution composites, in which the model of random inclusion distribution is gradually changing in statistics. As shown in Figure 1, the uniform random microstructure decides that material ‘A’ is the statistically homogeneous medium in macroscopic view. Material ‘B’, however, has apparently non-uniform random dispersions at the micro level, which leads to the macroscopic inhomogeneous feature in the statistical sense. Thus we call it as the statistically inhomogeneous material. For example, functionally ∗ Correspondence to: Fei Han, Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710129, People’s Republic of China. † E-mail: cjzgroup@mail.nwpu.edu.cn Contract/grant sponsor: National Basic Research Program of China; contract/grant number: 2010CB832702 Contract/grant sponsor: National Natural Science Foundation of China; contract/grant number: 90916027 Copyright 䉷 2010 John Wiley & Sons, Ltd. THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD 973 Figure 1. The schematic representation for the microstructure of random composites. A: statistically homogeneous material; and B: statistically inhomogeneous material. graded materials (FGM) are special cases of statistically inhomogeneous media. The microstructures that the inclusion dispersions gradually vary result in changes of the macroscopic physical, and mechanical properties of the materials. Therefore, the statistically inhomogeneous materials are of multifunction in physics and mechanics. So far, some methods to predict the physical and mechanical properties of random composite materials have been developed, such as the law of mixtures [1, 2], the Hashin–Shtrikman bounds [3], the self-consistent method [4, 5], the Eshelby’s equivalent inclusion method [6], and the Mori– Tanaka homogenization method [7]. Although these methods effectively promoted the development of computational material science, the microstructure of real materials was greatly simplified to reduce the theoretical complexity. Later, with the appearance of FGM, the early methods were improved to predict the effective properties of FGM. For instance, the traditional Mori–Tanaka method was extended to linearly variable fields [8]. According to the representative volume element, Aboudi et al. [9] developed a higher-order theory by coupling the local and global responses. The other micro-mechanics-based elastic models were constructed by considering the locally pair-wise interactions between particles and the gradient effects of particle volume fraction [10]. Recently, Ferrante et al. [11] and Rahman et al. [12] introduced the random variations into micro-mechanical models in view of the statistical uncertainties for FGM. On the other hand, some finite element methods based on the hexagonal cells [13], voronoi cells [14], and rectangular cells [15] were employed to predict the mechanical properties of functionally graded plates. Guild et al. [16–18] developed the representative cell model based on the finite element method for the particle-filled composites. On this basis, Poon et al. [19] presented the statistical spherical cell model for the elastic properties of composites. In addition, some methods involved three-dimensional cell models [18, 20, 21], but there was only one sphere in the cell. Later, the representative cell model with regular arrangements of some spheres was used in some numerical methods based on finite element analysis [22–24]. The regular arrangements of spheres in a representative cell are simpler than the microstructure of real materials. Recently, the cell model with random sphere dispersions has appeared [25, 26]. A few methods, however, take into account all the variation characteristics, including inclusion’s shape, size, orientation, spatial distribution, and volume fraction and so on. Besides, these micromechanics-based methods can be used to predict the macroscopic constitutive parameters, but cannot be employed to calculate the local strain and stress fields. Thus, the strength performances Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 974 F. HAN, J. CUI AND Y. YU Figure 2. The representative cells of statistically inhomogeneous materials. cannot yet be predicted by these methods. Recently, the authors proposed an effective computer algorithm to produce complex and realistic representative cells for random composite materials [27]. This algorithm is based on the probability distribution model, which takes into account the inclusion’s shape, size, orientation, spatial distribution, and volume fraction. Figure 2 shows the representative cells for the microstructures of statistically inhomogeneous materials. On the basis of the asymptotic expansion homogenization approach [28, 29], scientists have made some worthwhile contributions in mechanical calculation. Fish et al. generalized the mathematical homogenization theory for damage mechanical problems [30–33], molecular dynamics equations [34–36], and non-periodic heterogeneous media [37]. Chung et al. [38] and Zhang et al. [39], respectively, extended the asymptotic expansion approach for transient elasto-plastic response and thermodynamic wave propagation in periodic composite materials. Meanwhile, Cui et al. proposed a multi-scale method to predict the physical and mechanical properties of periodic composites [40]. For physics field problems of random composite materials, Jikov et al. [41] proved the existences of the homogenization coefficients and the homogenization solution. On the basis of the Monte Carlo method, Cui et al. presented a statistical second-order two-scale (SSOTS) analysis method to predict the physical or mechanical properties of statistically homogeneous materials [42, 43]. However, the previous second-order two-scale asymptotic expansion cannot be employed to the statistically inhomogeneous materials because the inclusion’s probability distributions are changing in space. In this paper, we introduce a correction term into the asymptotic expansion of the displacement field, and then derive a family of cell problems. A new SSOTS method is presented to predict the mechanical performances of the statistically inhomogeneous materials. This paper is organized as follows. In Section 2 the microstructure of the statistically inhomogeneous material is represented by setting up a statistic screen pointwise. Section 3 is devoted to the formulation of the SSTOS method and the algorithm procedure. Finally, the numerical results for the mechanical performances of some statistically inhomogeneous materials are shown. 2. MICROSTRUCTURE REPRESENTATION For composite materials with random inclusion dispersions, all the inclusions in the geometry are considered as ellipsoids or the polyhedrons inscribed inside ellipsoids. The size of each inclusion Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD 975 is denoted by the long axis a of ellipsoid. From the engineering survey and the statistic fitting method of data, the microstructure of random composites is represented as follows: (1) There exists a constant ε, a ≪ ε ≪ L, L is the size of structure . If at any point inside  there is a cell with the same size ε and the same probability distribution of inclusions, this composite material is called as the statistically homogeneous material. If at an arbitrary point x0 of  there exists a cell with size εx0 and a specific probability distribution model P(x0 ) depended on x0 , then the expected stiffness tensor is a function of x0 in the statistical sense. We call it the statistically inhomogeneous material. In addition, the size εx0 and the probability distribution are assumed to be continuously changing in space. (2) Each ellipsoid in the three-dimensional space is denoted by nine random parameters, including the coordinates (xc , yc , z c ) of the center, three Euler angles (, , ) of the rotations, and the sizes (a, b, c) of three axes. Let the random vector  = (x c , yc , z c , , , , a, b, c). Suppose that a cell ε Q sx0 with center at x0 contains N ellipsoids, where Q sx0 represents a normalized cell and s = 1, 2, 3, . . . denotes the index of samples, then its random sample is defined as sx0 = (s1 , s2 , s3 , . . . , sN −1 , sN ) ∈ P(x0 ). The stiffness tensor in cell ε Q sx0 is expressed as {aiεj hk (x, sx0 )}. If the composites are made from particles and matrix, then the material parameters of a sample can be defined as follows: ⎧ ⎪ 1 ⎪ ⎪ ai j hk ⎨ aiεj hk (x, sx0 ) = ⎪ ⎪ ⎪ ⎩ ai2j hk if x ∈ N  ei i=1 if x ∈ ε Q sx0 − N  ei i=1 where i, j, h, k = 1, . . . , n, ei denotes the ith ellipsoid inside ε Q sx0 , and ai1j hk and ai2j hk are the material coefficients of particles and matrix, respectively, and satisfy that max{|ai1j hk |, |ai2j hk |} < M, M > 0. Thus, the {aiεj hk (x, sx0 )} is a bounded and measurable tensor of random variables sx0 . 3. STATISTICAL SECOND-ORDER TWO-SCALE METHOD 3.1. Statistical second-order two-scale formulation In this section, a novel SSOTS formulation is derived formally at an arbitrary point of the structure for calculating the mechanical properties, including stiffness tensor, strain and stress tensors. For the structure  made from the statistically inhomogeneous material, the boundary value problem of its mechanical behavior can be expressed as follows:   * 1 *u εh (x, ) *u εk (x, ) ε − a (x, ) + 2 *x j i j hk *xk *x h = f i (x), uε (x, ) = ū(x),  j a εji hk (x, ) 1 2  Copyright 䉷 2010 John Wiley & Sons, Ltd. *u εh (x, ) *u εk (x, ) + *xk *x h = pi (x), x∈ x ∈ 1 (1) x ∈ 2 Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 976 F. HAN, J. CUI AND Y. YU where {aiεj hk (x, )} satisfies the conditions of symmetry ε aiεj hk (x, ) = a εji hk (x, ) = aiεjkh (x, ) = ahki j (x, ) and ellipticity ∃1 , 2 > 0, 1 ε i h i h ai j hk (x, ) i h jk 2 i h i h ∀ i h = hi And  = {sx0 , x0 ∈ ε Q sx0 ⊂ }; f i (x) is a body force assumed to be independent of ε; uε (x, ) is the displacement solution vector; 1 and 2 denote boundary portions where displacement ū(x) and loads pi (x) are prescribed, respectively, such that 1 ∩2 = , 1 ∪2 = *; v j ( j = 1, . . . , n) are the normal direction cosines of 2 . We introduce a variable = (x−x0 )/ε ∈ Q sx0 which denotes local coordinates on the normalized cell Q sx0 , and then the stiffness tensor can be expressed as aiεj hk (x, ) = ai j hk ( , sx0 ). As the displacement solution uε (x, ) of the boundary value problem (1) depends on both the global behaviors of the structure  and the local configuration in Q sx0 , it can be expressed as uε (x, ) = u(x, , sx0 ). Further, enlightened by [42] uε (x, ) is assumed to be the following expansion in two-scale variables x, 2 uε (x, ) = u0 (x)+εN1 ( , sx0 ) +ε 3 g(x, , sx0 ), * u0 (x) *u0 (x) *u0 (x) +ε 2 M1 (x, , sx0 ) +N1 2 (x, , sx0 ) *x1 *x1 *x1 *x2 x ∈ ε Q sx0 ⊂ , ∈ Q sx0 , sx0 ∈ P(x0 ) (2) where u0 (x) only reflects the macroscopic behaviors of the structure, and is called as the homogenization solution; N1 ( , sx0 ), M1 (x, , sx0 ), and N1 2 (x, , sx0 ) (1 , 2 = 1, . . . , n) involve the microstructural information, they are matrix-valued functions, and are expressed as follows: ⎞ ⎛ N1 11 ( , sx0 ) · · · N1 1n ( , sx0 ) ⎟ ⎜ ⎟ ⎜ .. .. N1 ( , sx0 ) = ⎜ ⎟ = (N1 1 , . . . , N1 n ) . · · · . ⎠ ⎝ s s N1 n1 ( , x0 ) · · · N1 nn ( , x0 ) ⎞ ⎛ M1 11 (x, , sx0 ) · · · M1 1n (x, , sx0 ) ⎟ ⎜ ⎟ ⎜ .. .. M1 (x, , sx0 ) = ⎜ ⎟ = (M1 1 , . . . , M1 n ) . · · · . ⎠ ⎝ M1 n1 (x, , sx0 ) · · · ⎛ ⎜ ⎜ N1 2 (x, , sx0 ) = ⎜ ⎝ N1 2 11 (x, , sx0 ) · · · .. . N1 2 n1 (x, , sx0 ) M1 nn (x, , sx0 ) N1 2 1n (x, , sx0 ) ··· .. . ··· N1 2 nn (x, , sx0 ) ⎞ ⎟ ⎟ ⎟ = (N1 2 1 , . . . , N1 2 n ) ⎠ N1 m ( , sx0 ), M1 m (x, , sx0 ), N1 2 m (x, , sx0 ) (1 , 2 , m = 1, . . . , n), and u0 (x) will be determined by a constructing way below. Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD 977 It is worth noting that the expansion (2) is different from the form given by [42], the differences are that a correction term ε 2 M1 (x, , sx0 )(*u0 (x))/(*x1 ) is constructed into the asymptotic expansion and M1 (x, , sx0 ) and N1 2 (x, , sx0 ) depend on the macroscopic variable x. Owing to = (x−x0 )/ε, the indirect macroscopic spatial derivatives can be calculated by the chain rule as */*xi = (*/*xi )+(1/ε)(*/* i ). Substituting (2) into (1) yields the equality   *u εh (x, sx0 ) *u εk (x, sx0 ) * s 1 − ai j hk ( , x0 ) + 2 *x j *xk *x h  *ai j hk ( , sx0 ) 1 *u 0h (x) *u 0k (x) −1 + =ε − 2 * j *xk *x h    *N1 hm ( , sx0 ) *N1 km ( , sx0 ) *u 0m (x) * s 1 − + ai j hk ( , x0 ) 2 * j * k * h *x1    *M1 hm (x, , sx0 ) *M1 km (x, , sx0 ) * *u 0m (x) 0 s 1 +ε − + ai j hk ( , x0 ) 2 * j * k * h *x1 2 2 −ai j hk ( , sx0 ) 1 2 * u 0h (x) * u 0k (x) + *xk *x j *x h *x j −ai j hk ( , sx0 ) 1 2 *N1 hm ( , sx0 ) * + *N1 km ( , sx0 ) k * h 2 * u 0m (x) *x1 *x j   2 0 2 0 * s * u m (x) s 1 s * u m (x) N1 hm ( , x0 ) − ai j hk ( , x0 ) + N1 km ( , x0 ) 2 * j *x1 *xk *x1 *x h   2  *N1 2 hm (x, , sx0 ) *N1 2 km (x, , sx0 ) * * u 0m (x) s 1 − ai j hk ( , x0 ) + 2 * j * k * h *x1 *x2 +O(ε) = f i (x) (3) Suppose that the equality (3) holds for any ε > 0, then from the coefficients of ε −1 in both sides of equality (3), the following equality is obtained: − *ai j hk ( , sx0 ) 1 * j 2 *u 0h (x) *u 0k (x) + *xk *x h  * 1 − ai j hk ( , sx0 ) 2 * j Copyright 䉷 2010 John Wiley & Sons, Ltd. *N1 hm ( , sx0 ) * k + *N1 km ( , sx0 ) * h  *u 0m (x) =0 *x1 (4) Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 978 F. HAN, J. CUI AND Y. YU By virtue of the symmetry of stiffness tensor, one obtains that  − *ai j 1 m ( , sx0 ) * j  * 1 − ai j hk ( , sx0 ) 2 * j *N1 hm ( , sx0 ) * + *N1 km ( , sx0 ) * k h  *u 0m (x) =0 *x1 (5) As *u 0m (x)/*x1 for m, 1 = 1, . . . , n are not identically zero, in order that equality (5) always holds one can construct local functions N1 m ( , sx0 ) for 1 , m = 1, . . . , n so that they are the solutions to the following problems:  1 * ai j hk ( , sx0 ) − 2 * j *N1 hm ( , sx0 ) * + *N1 km ( , sx0 ) k * h N1 m (  , sx0 ) = *ai j 1 m ( , sx0 ) * , j ∈ Q sx0 (6) ∈ *Q sx0 = 0, It is proved that problems (6) has a unique solution for any specified sample. Inspired by [42], for any sample sx0 s = 1, 2, 3, . . ., the homogenized stiffness tensor is defined as âi j hk (sx0 ) =  Q sx0 ai j hk ( , sx0 )+ai j pq ( , sx0 ) 1 2 *Nhpk ( , sx0 ) * + *Nhqk ( , sx0 ) q * d (7) p From the boundedness of ai j hk ( , sx0 ) and Lemma 3.2 in [42], it follows that âi j hk (sx0 ) is bounded with respect to random variable. Then there exists an expected value of âi j hk (sx0 ). Then from Kolmogorov’s strong law of large numbers, the expected homogenized stiffness tensor at x0 can be calculated by ⌢ a i j hk (x)|x=x0 = lim M→+∞ M s s=1 âi j hk (x0 ) (8) M ⌢ After the expected homogenized stiffness function a i j hk (x) is obtained in the structure , the homogenized problem is defined as follows:  * ⌢ 1 − a i j hk (x) 2 *x j *u 0h (x) *u 0k (x) + *xk *x h  = f i (x), u0 (x) = ū(x), ⌢  j a i j hk (x) 1 2 *u 0h (x) *u 0k (x) + *xk *x h = pi (x), x∈ x ∈ 1 (9) x ∈ 2 ⌢ Similar to [42], it is proved that a i j hk (x) is symmetric and positive definite. Thus the homogenized problem (9) has the unique homogenization solution u0 (x). Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 979 THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD As the right side of equality (3) is independent of ε, the coefficient of ε 0 on the left side of (3) equals f i (x), that is  * 1 − ai j hk ( , sx0 ) 2 * j *M1 hm (x, , sx0 ) * 2 + *M1 km (x, , sx0 ) * k h  *u 0m (x) *x1 2 −ai j hk ( , sx0 ) 1 2 * u 0h (x) * u 0k (x) + *xk *x j *x h *x j −ai j hk ( , sx0 ) 1 2 *N1 hm ( , sx0 ) * + *N1 km ( , sx0 ) * k h 2 * u 0m (x) *x1 *x j   2 0 2 0 * s * u m (x) s 1 s * u m (x) + N1 km ( , x0 ) ai j hk ( , x0 ) N1 hm ( , x0 ) − 2 * j *x1 *xk *x1 *x h   2 *N1 2 hm (x, , sx0 ) *N1 2 km (x, , sx0 ) * * u 0m (x) s 1 − ai j hk ( , x0 ) + 2 * j * k * h *x1 *x2 = f i (x) (10) Making use of the symmetry of {ai j hk ( , sx0 )} and the homogenized equation (9), equality (10) can be rewritten as:     ⌢ *M1 hm (x, , sx0 ) *M1 km (x, , sx0 ) *a i j 1 m (x) *u 0m (x) * s 1 − ai j hk ( , x0 ) + + 2 * j * k * h *x j *x1    *N1 2 hm (x, , sx0 ) *N1 2 km (x, , sx0 ) * s 1 ai j hk ( , x0 ) + + − 2 * j * k * h ⌢ +a i 2 m 1 (x)−ai 2 m 1 ( , sx0 )−ai 2 hk ( , sx0 ) − *N1 hm ( , sx0 ) * k  *2 u 0m (x) *  =0 ai j h 2 ( , sx0 )N1 hm ( , sx0 ) * j *x1 *x2  (11) By the constructing way analogous to determining N1 m ( , sx0 ) for 1 , m = 1, . . . , n in order that equality (11) always holds, one can let M1 m (x, , sx0 ), and N1 2 m (x, , sx0 ) (1 , 2 , m = 1, . . . , n) are the solutions of the following problems, respectively.   ⌢ *M1 hm (x, , sx0 ) *M1 km (x, , sx0 ) *a i j 1 m (x) * s 1 ai j hk ( , x0 ) =− + , ∈Q sx0 − 2 * j * k * h *x j (12) M1 m (x, , sx ) = 0, Copyright 䉷 2010 John Wiley & Sons, Ltd. ∈ *Q sx0 Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 980 F. HAN, J. CUI AND Y. YU  * 1 − ai j hk ( , sx0 ) 2 * j *N1 2 hm (x, , sx0 ) * * h ⌢ = −a i 2 m 1 (x) , sx0 ) *N1 hm ( (13) k  *  ai j h 2 ( , sx0 )N1 hm ( , sx0 ) , * j N1 2 m (x, , sx0 ) = 0, * k +ai 2 m 1 ( , sx0 )+ai 2 hk ( , sx0 ) + +  *N1 2 km (x, , sx0 ) ∈ Q sx0 ∈ *Q sx0 It is proved that problems (12) and (13) have the unique solutions M1 m (x, , sx0 ), and N1 2 m (x, , sx0 ), respectively. Summing up, one acquires the following theorem: Theorem The elasticity boundary value problem (1) for the structure  made from the statistically inhomogeneous material formally has the second-order two-scale approximate solution as follows: uε (x, ) ∼ = u0 (x)+εN1 ( , sx0 ) *u0 (x) *x1 2 +ε 2 M1 (x, , sx0 ) x ∈ ε Q sx0 ⊂ , * u0 (x) *u0 (x) +N1 2 (x, , sx0 ) *x1 *x1 *x2 ∈ Q sx0 , sx0 ∈ P(x0 ) (14) where u0 (x) is the solution of the homogenized problem (9), called as the homogenization solution, N1 m ( , sx0 ), M1 m (x, , sx0 ) and N1 2 m (x, , sx0 ) are local solutions of the localized problems (6), (12), and (13) defined on the normalized cell Q sx0 , respectively. Furthermore, the strains inside ε Q sx0 can be expressed in terms of the displacement expansion (14) as follows: ε εhk (x, sx0 ) = 1 2 *u εh (x, sx0 ) = 1 2 *u 0h (x) *u 0k (x) 1 + + 2 *xk *x h +ε 1 * u 0m (x) * u 0m (x) + N1 km ( , sx0 ) N1 hm ( , sx0 ) 2 *x1 *xk *x1 *x h +ε 1 2 *xk + *u εk (x, sx0 ) *x h *N1 hm ( , sx0 ) * * k Copyright 䉷 2010 John Wiley & Sons, Ltd. *N1 km ( , sx0 ) * k h *u 0m (x) *x1 2 2 *N1 2 hm (x, , sx0 ) + + *N1 2 km (x, , sx0 ) * h 2 * u 0m (x) *x1 *x2 Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 981 THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD *N1 2 hm (x, , sx0 ) *N1 2 km (x, , sx0 ) 2 * u 0m (x) *x1 *x2 +ε 2 1 2 +ε 2 * u 0m (x) 1 * u 0m (x) + N1 2 km (x, , sx0 ) N1 2 hm (x, , sx0 ) 2 *x1 *x2 *xk *x1 *x2 *x h *xk + *x h 3 3 +ε 1 2 *M1 hm (x, , sx0 ) * + *M1 km (x, , sx0 ) * k *M1 hm (x, , sx0 ) h *M1 km (x, , sx0 ) *u 0m (x) *x1 *u 0m (x) *x1 +ε 2 1 2 +ε 2 * u 0m (x) * u 0m (x) 1 + M1 km (x, , sx0 ) M1 hm (x, , sx0 ) 2 *x1 *xk *x1 *x h *xk + *x h 2 2 (15) From Hooke’s law, the stresses inside ε Q sx0 are evaluated as ε s ε s ε s i j (x, x0 ) = ai j hk (x, x0 )εhk (x, x0 ) (16) Then the strains and stresses at any point x are naturally adopted to evaluate the elastic limit value S(sx0 ) of statistically inhomogeneous materials. 3.2. Elastic limit of statistically inhomogeneous materials The elastic limit is obtained by evaluating the critical point of strains/or stresses inside numerically experimental component , for example, tensile column, bending beam or twisting bar, according to the proper strength criterion of materials. The strength criteria are usually different for particle’s material and matrix. On the basis of the microstructural characteristic at an arbitrary position x of component , we generate the sample distribution Q x (sx ) (s = 1, 2, 3, . . .) and calculate strains and stresses inside ε Q sx ⊂ , and then evaluate the elastic limit value S(sx ). Because plenty of samples Q x (sx ) (s = 1, 2, 3 . . .) are needed to be generated in statistics, a set of elastic limit values S(sx ) can be obtained. Naturally, through computing the set of S(sx ) by using Kolmogorov’s strong law of large numbers, the expected elastic limit Ŝ of the statistically inhomogeneous material is obtained as follows: M S(sx ) (17) Ŝ = s=1 M where M is the number of samples. However, the expected elastic limit Ŝ is only an average value; it does not sufficiently manifest the elastic limit performance of materials. Therefore, the minimum value of all samples’ elastic limits S(sx ) should be defined in the engineering practice. It is expressed as follows: Smin = Copyright 䉷 2010 John Wiley & Sons, Ltd. min s=1, ...,M S(sx ) (18) Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 982 F. HAN, J. CUI AND Y. YU 3.3. Algorithm procedure The algorithm procedure of the SSTOS method for predicting the mechanical properties of statistically inhomogeneous materials is stated as follows: 1. Select N points xi (i = 1, 2, . . . , N ) inside the structure  according to the structural topology and the random distribution model P(x) (x ∈ ) of inclusions. 2. Generate a sample distribution Q xi (sxi ) at point xi (i = 1, 2, . . . , N ) according to P(x), and then partition Q xi (sxi ) into its finite element set. 3. Solve problem (6) in Q xi (sxi ) to obtain N1 m ( , sxi ) (1 , m = 1, . . . , n), and then evaluate the homogenized coefficient âi j hk (sxi ) by the formula (7). 4. Repeat steps 2 and 3 for different samples Q xi (sxi ) (s = 1, 2 . . . , M), M is the number of samples, and M homogenized coefficients âi j hk (sxi ) are obtained. Then the expected ⌢ homogenized coefficient a i j hk (xi ) is calculated by (8). 5. Repeat steps 2 to 4 for every point xi (i = 1, 2, . . . , N ). Then N expected homogenized ⌢ coefficients a i j hk (xi ) are obtained. And then the expected homogenized coefficient function ⌢ a i j hk (x)(x ∈ ) is constructed by N -points interpolation. ⌢ 6. Define problem (9) in the whole structure  using a i j hk (x)(x ∈ ) obtained in step 5, and solve it to obtain the homogenization displacement field u0 (x) by the finite element method. 7. For sample Q xi (sxi ), evaluate M1 m (x, , sxi ), and N1 2 m (x, , sxi ) by solving problems (12) and (13), using the same finite element meshes and stiffness matrix as in step 3. 8. For sample Q xi (sxi ), using N1 m ( , sxi ), M1 m (x, , sxi ), N1 2 m (x, , sxi ), and u0 (x) to calculate the stain field and stress field by (15) and (16). 9. Evaluate the elastic limit value S(sxi ) according to the strength criteria of matrix and particles. 10. Repeat the steps from 7 to 9 for different samples Q xi (sxi ) (s = 1, 2, . . . , M), and M elastic limits S(sxi ) are obtained. Then the expected elastic limit Ŝ and minimum elastic limit Smin are evaluated by (17) and (18). 4. NUMERICAL EXPERIMENTS AND RESULTS In order to verify that the presented algorithm is feasible and effective to predict the mechanical properties of the statistically inhomogeneous materials, we have developed computer programs and performed some numerical experiments. Here some numerical results are shown and compared with experimental findings. Example 1: This example is on the stiffness calculation for the tensile column made of cenospheres reinforced polyester resin [44]. The volume fraction of cenospheres decreases from about 50 to 0% along the tensile direction, shown as in Figure 3. The material properties of cenospheres and resin matrix are: cenospheres E C = 6.0 GPa; C = 0.35; polyester E P = 3.6 GPa; P = 0.41; polyester-plasticizer E PP = 2.5 GPa; PP = 0.33. According to the variation curve of the cenosphere volume fraction in paper [44], the changing curve of stiffness of this FGM is obtained by using the SSTOS method. The following results are based on the statistics of 50 samples. The tensile moduli at some positions of FGM are calculated first, and then the spatial variations of modulus are drawn by interpolation as shown in Figure 4. The figure shows that the predicted Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD 983 Figure 3. The tensile column of FGM. Figure 4. Comparisons of Young’s moduli between the numerical results and experimental data. modulus increases as the position extends, also as the volume faction increases. The numerical curves are close to the simulative curves by the micro-mechanics-based elastic model (MEM) in [10], and agree with the experimental data [44]. Example 2: This example is on the polymer blends in injection molded part [45, 46]. The shape and size of PET particles are gradually changing in PE matrix because of melt flow, shown as in Figure 5. The left of the figure is close to the injection gate and the right is at the injection end. The schematic figure is divided into three main zones: sub-skin, intermediate, and core zone [45, 46]. Because of limited observed data, we choose only six representative points, which are located in three different zones near the injection gate and end, as shown in Figure 5. According to the observed results [45, 46] of particle distributions in low injection speed, the three-dimensional cells are generated, and then the numerical results of six positions in a tensile part are calculated, as shown in Table I. The mechanical properties of PET and PE are: E PET = 2900 MPa; vPET = 0.40; PET strength= 78 MPa; E PE = 985 MPa; vPE = 0.42; PE strength= 20.5 MPa. Because the PET is Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 984 F. HAN, J. CUI AND Y. YU Figure 5. The schematic microconfiguration of PET/PE polymer blends in injection molded part. Table I. SSOTS results of mechanical properties. G–G section N–N section Cell’s location Sub-skin (G1) Intermediate (G2) Core (G3) Sub-skin (N1) Intermediate (N2) Core (N3) Average diameter (m) 0.8 1121.9 1121.1 1138.8 393.7 395.9 396.5 0.419 0.410 0.410 0.9 1128.1 1127.3 1136.0 396.4 397.8 398.4 0.417 0.413 0.413 1.2 1132.3 1131.2 1134.9 398.4 398.9 399.4 0.416 0.414 0.414 1.2 1127.9 1127.4 1139.7 396.3 398.3 398.5 0.417 0.412 0.412 1.5 1131.6 1131.3 1140.8 397.9 399.6 399.9 0.417 0.412 0.412 1.9 1140.2 1140.1 1141.2 401.8 401.9 402.1 0.415 0.414 0.414 31.70 30.94 30.36 31.19 30.61 28.98 31.18 30.39 28.83 30.36 29.24 27.64 E1 (MPa) E2 (MPa) E3 (MPa) G12 (MPa) G23 (MPa) G13 (MPa) v12 v23 v13 Expected elastic limit (MPa) Minimum elastic limit (MPa) mixed with PE in a fixed weight ratio of 15/85 in the experiments [45], the volume fraction of PET is about 11.4% evaluated by the densities of PET and PE, which are 1.32 and 0.96 g/cm3 , respectively. The results in the columns of Table I show that the homogenized material properties are orthogonal-anisotropic. In fact, subscript ‘3’ denotes the melt flow direction. The PET particles are elongated along the ‘3’ direction during the melt flow. The modulus of particle is larger than that of the matrix, and the Poisson ratio of particle is smaller than that of matrix. Therefore, the Young modulus and shear modulus in the ‘3’ direction are all larger and the Poisson ratio is smaller. On the other hand, the data in the rows show that the moduli at the injection end are all larger than that at the injection gate with the size of particles increasing. The subscript ‘3’ also denotes the main loading direction of the injection modeled part. In this direction, the expected elastic limit and Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 985 THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD Figure 6. The flexural beam of hollow microballoons-filled epoxy resin. Table II. Mechanical properties with different particle volume fraction and wall thickness ratio. Particle volume fraction 30% Wall thickness E ratio (%) (MPa) (MPa) v Expected flexural elastic limit (MPa) 3 4.5 5.5 6.5 1196 1364 1471 1567 0.3202 0.3187 0.3170 0.3151 47.16 55.12 58.95 62.82 3115 3526 3791 4030 G Particle volume fraction 40% Minimum flexural elastic limit (MPa) E G (MPa) (MPa) 36.65 43.52 53.08 51.35 3136 3672 4018 4342 1220 1444 1587 1721 v Expected flexural elastic limit (MPa) Minimum flexural elastic limit (MPa) 0.3114 0.3099 0.3082 0.3063 43.52 49.31 55.71 59.72 35.35 39.63 48.13 50.54 minimum elastic limit increase as the particle is smaller and more elongated in the three different morphology zones. This example simulated the real experiments [45, 46] and supported the experimental investigation effectively. Table I shows that the moduli are indeed orthotropic, but only very weakly so. The difference between moduli in different directions is only about 1%. This is due to small differences between input data, where the ratio of moduli of particle and matrix is about 3:1, the average length-to-diameter ratio of particles at ‘G1’ position is less than 2:1. In fact, we changed the modulus ratio of particle and matrix into 10:1, and stretched the average length-to-diameter ratio of particle to 10:1. The new results show that the difference between moduli in different directions can reach 33%, which shows obviously orthotropic. These findings suggest that the morphology of mixed particles is important for the mechanical properties of blends. Example 3: The flexural properties of functionally graded beams are calculated in this example [47]. These FGM beams made of hollow glass microballoon-filled epoxy resins, and the wall thickness of the microballoon has a gradient change, as shown in Figure 6. The flexural moduli and strengths are calculated at four positions, where the ratios of wall thickness to radius are 3, 4.5, 5.5, and 6.5%, respectively. We consider two beams with volume fractions 30 and 40% of microballoons. The Young modulus, the Poisson ratio and flexural strength of epoxy matrix are 3.0 GPa, 0.35, and 68 MPa, respectively. In addition, the Young modulus, Poisson ratio, and the strength of glass are 76.0 GPa, 0.23, and 20 MPa. In this paper, the finite element meshes of hollow microballoons are generated similar to that of core-shell particles in [43]. Table II shows that the mechanical properties of the beam are affected by the microstructure and volume fraction of microballoons. On the one hand, the data in the rows show that the Young modulus and the shear modulus increase, but the Poisson ratio, the expected flexural elastic limit, Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 986 F. HAN, J. CUI AND Y. YU and the minimum limit decrease as the volume fraction of microballoon increases. On the other hand, as the volume fraction is constant in the columns, the moduli increase and the Poisson ratio decreases, but the expected flexural elastic limit and the minimum limit increase with the wall thickness rising. The changing trends of the numerical results agree with that of the experimental data [47]. 5. CONCLUSION A SSOTS method has been developed for predicting the mechanical properties of statistical inhomogeneous materials. The solution of the displacement obtained by the SSOTS method is equivalent to the true solution of the original equations in the nearly pointwise sense with ε order in a cell. The expected homogenized stiffness tensor has been found to be a function of the space coordinates. It means that the macroscopic properties of the investigated materials are not constant in the statistical sense. In other words, it is the statistically inhomogeneous material. In addition, a workable SSOTS algorithm has been constructed. The SSOTS numerical results have shown that the microstructure has a marked effect on the macroscopic mechanical properties. Specifically, these properties vary with the probability model of random inclusion dispersions. As a result, the information about the microstructure can be captured exactly by the SSTOS method. The SSTOS method can be extended to deal with other statistically inhomogeneous materials, such as cancellous bone and short fiber reinforced composites. ACKNOWLEDGEMENTS This paper is supported by the special funds for the National Basic Research Program of China (973 Program 2010CB832702), the National Natural Science Foundation of China (90916027), and also by the State Key Laboratory of Science and Engineering Computing, and the Center for High Performance Computing, Northwestern Polytechnical University in China. REFERENCES 1. Voigt W. Lehrbuch der Kristallphysik. Teubner: Leipzig, 1928. 2. Reuss A. Berechnung der Fliessgrenze von Mischkristallen auf Grunder Plastizitätsbedingung fRr Einkristalle. Zeitschrift für Angewandte Mathematik und Mechanik 1929; 9:49–58. 3. Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behavior of multiphase materials. Journal of the Mechanics and Physics of Solids 1963; 11:127–140. 4. Hill R. A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids 1965; 13:213–222. 5. Budiansky B. On the elastic moduli of some heterogeneous materials. Journal of the Mechanics and Physics of Solids 1965; 13:223–227. 6. Eshelby JD. The elastic field outside an ellipsoidal inclusion. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 1959; 252:561–569. 7. Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusion. Acta Metallurgica 1973; 21:571–574. 8. Zuiker J, Dvorak G. The effective properties of functionally graded composites—I. Extension of the Mori–Tanaka method to linearly varying fields. Composites Engineering 1994; 4:19–35. 9. Aboudi J, Pindera MJ, Arnold SM. Higher-order theory for functionally graded materials. Composites Part B: Engineering 1999; 30:777–832. Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme THE STATISTICAL SECOND-ORDER TWO-SCALE METHOD 987 10. Yin HM, Sun LZ, Paulino GH. Micromechanics-based elastic model for functionally graded materials with particle interactions. Acta Materialia 2004; 52:3535–3543. 11. Ferrante FJ, Graham-Brady LL. Stochastic simulation of non-Gaussian/non-stationary properties in a functionally graded plate. Computer Methods in Applied Mechanics and Engineering 2005; 194:1675–1692. 12. Rahman S, Chakraborty A. A stochastic micromechanical model for elastic properties of functionally graded materials. Mechanics of Materials 2007; 39:548–563. 13. Reiter T, Dvorak GJ, Tvergaard V. Micromechanical models for graded composite materials. Journal of the Mechanics and Physics of Solids 1997; 45:1281–1302. 14. Grujicic M, Zhang Y. Determination of effective elastic properties of functionally graded materials using Voronoi cell finite element method. Materials Science and Engineering A 1998; 251:64–76. 15. Cho JR, Ha DY. Averaging and finite-element discretization approaches in the numerical analysis of functionally graded materials. Materials Science and Engineering A 2001; 302:187–196. 16. Guild FJ, Young RJ, Lovell PA. The influence of material properties on the predicted behaviour of rubbertoughened polymers. Journal of Materials Science Letters 1994; 13:10–14. 17. Guild FJ, Kinloch AJ. Modelling the properties of rubber-modified epoxy polymers. Journal of Materials Science 1995; 30:1689–1697. 18. Guild FJ, Bonfield W. Predictive modelling of the mechanical properties and failure processes in hydroxyapatitepolyethylene (HapexTM ) composite. Journal of Materials Science: Materials in Medicine 1998; 9:497–502. 19. Poon YM, Luk WL, Shin FG. Statistical spherical cell model for the elastic properties of particulate-filled composite materials. Journal of Materials Science 2002; 37:5095–5099. 20. Tsui CP, Tang CY, Lee TC. Finite element analysis of polymer composites filled by interphase coated particles. Journal of Materials Processing Technology 2001; 117:105–110. 21. Tsui CP, Chen DZ, Tang CY, Uskokovi¢ PS, Fan JP, Xie XL. Prediction for debonding damage process and effective elastic properties of glass-bead-filled modified polyphenylene oxide. Composites Science and Technology 2006; 66:1521–1531. 22. Marur PR. Estimation of effective elastic properties and interface stress concentrations in particulate composites by unit cell methods. Acta Meterialia 2004; 52:1263–1270. 23. Balać I, Milovančević M, Tang CY, Uskoković PS, Uskoković DP. Estimation of elastic properties of a particulate polymer composite using a face-centered cubic FE model. Materials Letters 2004; 58:2437–2441. 24. Kang GZ, Shao XJ, Guo SJ. Effect of interfacial bonding on uniaxial ratchetting of SiC p /6061Al composites: finite element analysis with 2-D and 3-D unit cells. Materials Science and Engineering A 2008; 487:431–444. 25. Sun CJ, Saffari P, Sadeghipour K, Baran G. Effects of particle arrangement on stress concentrations in composites. Materials Science and Engineering A 2005; 405:287–295. 26. Sun CJ, Saffari P, Ranade R, Sadeghipour K, Baran G. Finite element analysis of elastic property bounds of a composite with randomly distributed particles. Composites Part A: Applied Science and Manufacturing 2007; 38:80–86. 27. Yu Y, Cui JZ, Han F. An effective computer generation method for the composites with random distribution of large numbers of heterogeneous grains. Composites Science and Technology 2008; 68:2543–2550. 28. Bensousson A, Lions JL, Papanicolaou G. Asymptotic Analysis for Periodic Structure. North-Holland: Amsterdam, 1978. 29. Sanchez-Palencia E. Non-homogeneous Media and Vibration Theory. Lecture Notes in Physics, vol. 127. Springer: Berlin, 1980. 30. Fish J, Yu Q, Shek K. Computational damage mechanics for composite materials based on mathematical homogenization. International Journal for Numerical Methods in Engineering 1999; 45:1657–1679. 31. Fish J, Yu Q. Two-scale damage modeling of brittle composites. Composites Science and Technology 2000; 61:2215–2222. 32. Fish J, Yu Q. Multiscale damage modelling for composite materials: theory and computational framework. International Journal for Numerical Methods in Engineering 2001; 52:161–191. 33. Fish J, Yu Q. Computational mechanics of fatigue and life predictions for composite materials and structures. Computer Methods in Applied Mechanics and Engineering 2002; 191:4827–4849. 34. Chen W, Fish J. A mathematical homogenization perspective of virial stress. International Journal for Numerical Methods in Engineering 2006; 67:189–207. 35. Fish J, Chen W, Li R. Generalized mathematical homogenization of atomistic media at finite temperatures in three dimensions. Computer Methods in Applied Mechanics and Engineering 2007; 196:908–922. 36. Li A, Li R, Fish J. Generalized mathematical homogenization: from theory to practice. Computer Methods in Applied Mechanics and Engineering 2008; 197:3225–3248. Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme 988 F. HAN, J. CUI AND Y. YU 37. Fish J, Fan R. Mathematical homogenization of nonperiodic heterogeneous media subjected to large deformation transient loading. International Journal for Numerical Methods in Engineering 2008; 76:1044–1064. 38. Chung PW, Tamma KK, Namburu RR. A computational approach for multi-scale analysis of heterogeneous elasto-plastic media subjected to short duration loads. International Journal for Numerical Methods in Engineering 2004; 59:825–848. 39. Zhang HW, Zhang S, Bi JY, Schrefler BA. Thermo-mechanical analysis of periodic multiphase materials by a multiscale asymptotic homogenization approach. International Journal for Numerical Methods in Engineering 2007; 69:87–113. 40. Yu XG, Cui JZ. The prediction on mechanical properties of 4-step braided composites via two-scale method. Composites Science and Technology 2007; 67:471–480. 41. Jikov VV, Kozlov SM, Oleinik OA. Homogenization of Differential Operators and Integral Functions. Springer: Berlin, 1994. 42. Li YY, Cui JZ. The multi-scale computational method for mechanics parameters of the materials with random distribution of multi-scale grains. Composites Science and Technology 2005; 65:1447–1458. 43. Han F, Cui JZ, Yu Y. The statistical two-order and two-scale method for predicting the mechanics parameters of core-shell particle-filled polymer composites. Interaction and Multiscale Mechanics 2008; 1(2):231–250. 44. Parameswaran V, Shukla A. Processing and characterization of a model functionally gradient material. Journal of Materials Science 2000; 35:21–29. 45. Li ZM, Yang W, Yang SY, Huang R, Yang MB. Morphology-tensile behavior relationship in injection molded poly(ethylene terephthalate)/polyethylene and polycarbonate/polyethylene blends (I) Part I Skin-core structure. Journal of Materials Science 2004; 39:413–431. 46. Li ZM, Xie BH, Yang SY, Huang R, Yang MB. Morphology-tensile behavior relationship in injection molded poly(ethylene terephthalate)/polyethylene and polycarbonate/polyethylene blends (II) Part II Tensile behavior. Journal of Materials Science 2004; 39:433–443. 47. Gupta N, Gupta SK, Mueller BJ. Analysis of a functionally graded particulate composite under flexural loading conditions. Materials Science and Engineering A 2008; 485:439–447. Copyright 䉷 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2010; 84:972–988 DOI: 10.1002/nme