Determination of RVE size for random composites with local volume fraction variation

Accepted Manuscript Determination of RVE size for random composites with local volume fraction variation Dimitrios Savvas, George Stefanou, Manolis Papadrakakis PII: DOI: Reference:

S0045-7825(16)30082-2 http://dx.doi.org/10.1016/j.cma.2016.03.002 CMA 10873

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date: 28 October 2015 Revised date: 30 January 2016 Accepted date: 2 March 2016 Please cite this article as: D. Savvas, G. Stefanou, M. Papadrakakis, Determination of RVE size for random composites with local volume fraction variation, Comput. Methods Appl. Mech. Engrg. (2016), http://dx.doi.org/10.1016/j.cma.2016.03.002 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Determination of RVE size for random composites with local volume fraction variation Dimitrios Savvasa , George Stefanoub,∗, Manolis Papadrakakisa a

Institute of Structural Analysis & Antiseismic Research, National Technical University of Athens, 9 Iroon Polytechneiou, Zografou Campus, 15780 Athens, Greece b Institute of Structural Analysis & Dynamics of Structures, Department of Civil Engineering, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece

Abstract In this paper, a novel computational procedure is proposed for the determination of the representative volume element (RVE) size for random composites, which is the basis of homogenization methods. The proposed approach takes into account the local volume fraction variation by processing computersimulated images of composites with randomly scattered inclusions. A variable number of microstructure models are derived directly from the images using the moving window technique. Each microstructure model contains different amount of inclusions, measured by image analysis tools. The proposed computational procedure is based on the solution of multiple boundary value problems under kinematic and static uniform boundary conditions using the extended finite element method (XFEM) coupled with Monte Carlo simulation (MCS). In this way, the statistical characteristics of the upper and lower bounds of the apparent material properties are computed. The RVE is attained within a prescribed tolerance by examining the convergence of these two bounds with respect to the mesoscale size. The numerical results clearly demonstrate the significant effect of matrix/inclusion stiffness ratio and volume fraction on the convergence rate. The proposed computational procedure leads to RVEs which can give more realistic results when used in homogenization problems. Mesoscale random fields describing the spatial variation of the components of the apparent elasticity tensor are also determined using the proposed approach. ∗

Corresponding author Email address: [email protected] (George Stefanou)

Preprint submitted to Computer methods in applied mechanics and engineeringMarch 19, 2016

Keywords: Random Composites, Apparent Properties, Mesoscale random fields, RVE size, XFEM, Monte Carlo Simulation 1. Introduction The macroscopic mechanical and physical properties of heterogeneous materials such as composites, concrete, polycrystals, solid foams and soils can be efficiently determined using either analytical [e.g. 1, 2, 3, 4, 5, 6] or numerical homogenization techniques [e.g. 7, 8, 9, 10, 4, 11, 5, 12, 13]. The application of homogenization methods in the numerical analysis of multiphase materials is much more advantageous over direct simulation of the full microstructure in terms of required computational resources [14]. In addition, homogenization of continua can predict the full multiaxial properties and responses of heterogeneous materials, which are often anisotropic, and thus expensive testing methods for measuring such material properties are avoided. The key issue in homogenization methods is the linking of micromechanical characteristics with the random variation of material properties at the macroscale [15, 16, 17, 18]. This link can be established using the Hill-Mandel macro-homogeneity condition [19] among other asymptotic or heterogeneous multiscale methods [20, 21]. In all these methods, the identification of a representative volume element (RVE) of the heterogeneous material is required over which a fine-scale boundary value problem is solved. Then effective material properties can be calculated using the link between the fine and coarse scale equations. According to Hill [19], the RVE is a sample that is structurally entirely typical of the whole mixture on average, and contains a sufficient number of inclusions for the apparent overall moduli to be effectively independent of the surface values of traction and displacement, as long as these values are macroscopically uniform. Following Hill’s definition, it has been pointed out that the RVE is clearly defined in two situations only: (i) unit cell in a periodic microstructure, and (ii) volume containing a very large (mathematically infinite) set of microscale elements (e.g. grains), possessing statistically homogeneous and ergodic properties [22]. While the concept of unit cell is valid whenever periodicity of microstructure is present, in the case of spatial randomness (e.g. random position, shape, size of inclusions), the identification of RVE must be based on computational convergence schemes with respect to specific apparent properties. 2

The typical procedure consists in setting multiple realizations of the microstructure followed by finite element simulation and statistical analysis of the results [23, 24, 25, 26]. Ostoja-Starzewski [22] attained an approximate RVE using convergence of apparent properties calculated on mesoscale models with increasing size. This occurred in terms of two hierarchies of bounds stemming from Dirichlet and Neumann boundary value problems, respectively. Trias et al. [27] reviewed the criteria which are used to define the minimal required size for a statistical RVE (SRVE) for fibre-reinforced composites considering matrix cracking failure. In Gitman et al. [28], a procedure based on the Chi-square criterion has been presented in order to determine the RVE size. Xu and Chen [29] used generalized variational principles to quantify the RVE size for 2-D and 3-D random two-phase composites taking into account the effect of volume fraction and stiffness ratio. Ma et al. [30, 31] combined FEM and MCS in a numerical convergence scheme of a homogenized parameter in order to determine a suitable size for the RVE. They also investigated the influence of randomness and correlation of microstructural parameters on the effective properties. In the context of stochastic finite element analysis of composite structures, the RVE size determines the minimum element size that should be used for discretization at the macroscale [32, 33, 34]. Since the volume fraction is one of the primary microstructural features in homogenization problems, the constant volume fraction assumption can lead to inexact effective properties. In this work, a novel computational scheme is proposed for the determination of RVE size for particle-reinforced composites based on computer-simulated images of their microstructure. A variable number of microstructure models are directly constructed from the composite material image using the moving window technique. In this way, the scatter of inclusions for each sample model is not achieved in an arbitrary manner but is an exact geometrical representation of the corresponding window segment of the composite. Moreover, the volume fraction of inclusions is computed through digital image processing of the microstructure models. The proposed numerical procedure takes into account the local volume fraction variation and thus the obtained RVE can lead to more accurate homogenized properties. This paper couples the extended finite element method (XFEM) with Monte Carlo simulation (MCS) in order to analyze the microstructure models and obtain statistical information for the constitutive properties of the samples in each window size. The number of MCS depends on the moving 3

window technique and is increasing for small window sizes. The implementation of XFEM is particularly suitable for this type of problems since there is no need to generate conforming finite element (FE) meshes at each MCS [35, 36, 37]. The weak discontinuity in the solution due to material interfaces is captured using nodal enrichment functions within the framework of the partition of unity method to augment the FE approximations over a structured mesh [38]. The XFEM analysis of the microstructure models is then used to compute upper and lower bounds for their constitutive behavior by solving multiple boundary value problems under kinematic and static uniform boundary conditions [10]. The RVE is attained within a prescribed tolerance by examining the convergence of these two bounds with respect to the mesoscale size. In the numerical examples, the effect of matrix/inclusion stiffness ratio as well as of the volume fraction of inclusions on the RVE size is assessed. The proposed computational scheme is applied to composites containing either stiff or compliant inclusions leading to much slower convergence in the latter case. Mesoscale random fields describing the spatial variation of the components of the apparent elasticity tensor are also determined. Since the numerical procedure is based on image analysis of computer-simulated microstructures, it can provide with a more realistic RVE and can be extended to any filler-reinforced composite. The paper is organized as follows. The computation of local volume fraction variation in composite materials by using the moving window technique is described in section 2. The computational procedure for the determination of the RVE size using XFEM is detailed in section 3. Numerical results along with a detailed discussion on the effect of volume fraction and stiffness ratio on the RVE size are provided in section 4. The major conclusions of the paper are finally presented in section 5. 2. Computation of local volume fraction variation The variability of local volume fraction (vf) in computer-simulated composites is studied in this section exploiting image analysis techniques available in MATLAB software. Specifically, images of two composites containing an initial amount (vf) of 40% and 20% of inclusions, respectively, are processed (see Fig. 1). Both images have dimensions Limage × Limage (Limage = 640µm) and an approximate analysis of 5800 × 5800 pixels. The represented composites contain 4096 and 2048 non-overlapping randomly scattered inclusions, 4

Figure 1: Computer-simulated images of composite materials containing circular inclusions with a) vf=40% (image 1) and b) vf=20% (image 2).

respectively. Each inclusion has circular shape with constant area 40 µm2 which is equivalent to a constant diameter d ' 7.14 µm. The vf computed from images 1 and 2 are 37.61% and 18.94% respectively, which is less than the initial values of vf for two reasons: (i) some of the randomly positioned inclusions are cut by the square border of the images and (ii) computational errors due to image resolution. Note that, instead of these computer-simulated images used here for convenience reasons, any SEM image of the composite could be used but in that case the error in the computation of vf would be amplified due to image noise. The computation of the local volume fraction variation is based on the moving window technique. In this method, an initial window of area L × L is set at a starting point O of the image and then, by moving this window over the image by a vector ξ~p = ξxp e~x + ξyp e~y a set of windows of the same size can be extracted (see Fig. 2). In this paper, the moving distance step ∆ξ is assumed to be the same along both directions. This means that ξxp = i∆ξ and ξyp = j∆ξ, with i, j the number of steps along the x and y directions, respectively. The total number of windows nw belonging to each set depends on the image size Limage , the investigated window size L and the selected

5

moving distance step ∆ξ as follows: 2  Limage − L +1 nw = ∆ξ

(1)

By choosing ∆ξ = L with L an integer divisor of Limage the whole image is segmented into nw = (Limage /L)2 non-overlapping windows (mesoscale models), as shown in Fig. 3. Specifically 10 sets containing nw =4, 16, 25, 64, 100, 256, 400, 1024, 1600, 4096 windows have been generated by choosing ∆ξ = L= 320, 160, 128, 80, 64, 40, 32, 20, 16, 10 µm, respectively. The vf of inclusions is computed in each window by using image analysis tools. The histograms of the vf corresponding to each set of windows are shown in Fig. 4 for the two images of Fig. 1. As the dispersion of the inclusions within the composite materials is not uniform, it is obvious that there are regions rich or poor in inclusions. This is clearly illustrated in the histograms where a very large variability of vf is observed for the sets containing windows of small size. In Fig. 5, the minimum and maximum vf computed from each set of windows are plotted and compared with the average vf in the set. The variation of vf is obvious in this figure because of the growing difference between the minimum and maximum vf as L is decreasing. In particular, in Fig. 5(b) corresponding to image 2 of Fig. 1, the minimum vf is zero for sets of L=10, 16 and 20 µm, which means that there is at least one window belonging to these sets which does not contain any inclusion. On the other hand, there are windows very rich in inclusions (vf>60%). As the size L increases, which means that the sets contain larger windows, the minimum and maximum vf are approaching the average value. For these cases, the empirical PDFs of vf become very narrow and tend to a spike (see Fig. 4). Note also in Fig. 5 that the average vf is approximately equal to the vf computed on the whole image, irrespectively of the size L. As shown above, for the composite materials studied in this paper, the vf varies significantly with the size and location of the windows. Therefore, in the context of homogenization methods, choosing a RVE size for composite materials with spatial randomness without considering local vf variability, can lead to unrealistic estimations of their mechanical behavior [39, 40].

6

Figure 2: Illustration of the moving window technique.

7

Figure 3: Generation of a set of nw = 256 non-overlapping windows of area L × L (L = 40 µm) and detail of the local microstructure in a window (mesoscale model).

Figure 4: Histograms of vf for various window sizes L for a) image 1, b) image 2.

8

Figure 5: Variation of vf with respect to window size L for a) image 1, b) image 2

3. Computational procedure for RVE determination According to Hashin [41] the homogenization process is based on the fundamental assumption of statistical homogeneity of the heterogeneous medium. This means that all statistical properties of the state variables are the same at any material point and thus a RVE can be identified. While in case of composites with periodic or nearly periodic geometry the RVE is explicitly defined, in case of spatial randomness, the existence of RVE has to be sought using computational methods. In this paper, identification of RVE is performed using a computational procedure based on Hill’s definition [19] which postulates separation of scales in the form: d  L  Lmacro

(2)

In the above inequality, the microscale parameter d denotes a characteristic size of the fillers, e.g. their diameter in case of circular inclusions, the mesoscale parameter L denotes the size of the RVE and the macroscale parameter Lmacro denotes the characteristic length over which the macroscopic loading varies in space, or in the case of complete scale separation, the size of the macroscopic homogeneous medium (see Fig. 6). In the proposed computational procedure, the RVE will be defined with respect to the non-dimensional parameter δ = L/d with δ ∈ [1, ∞]. Note that in the absence of spatial periodicity in composite materials the RVE is exactly obtained only in the limit δ → ∞. However, the RVE can be 9

Figure 6: Illustration of scale separation.

attained at a finite mesoscale size δ when a particular apparent property derived from statistical volume elements (SVEs) is almost approaching to a constant effective property. At this mesoscale size, the estimation of the particular effective property is not changing with increasing the number of realizations. In this paper, the RVE is identified in terms of the convergence of the components of the spatially averaged apparent elasticity tensor, which are calculated on mesoscale models of increasing size by applying both kinematic and static uniform boundary conditions. In each investigated mesoscale size δ, variability of vf is taken into account in the computational procedure by processing microstructure models directly extracted from images of the particular composite material (section 2). The RVE obtained from this computational approach is expected to lead to more realistic homogenization results. In the numerical examples of section 4, it will be shown that the convergence rate of the mesoscale models to the RVE is mainly depending on the matrix/inclusion stiffness ratio and volume fraction of the composite material.

10

3.1. Problem formulation Let us denote a mesoscale realization of a composite material as Bδ (θ) and its boundary surface as ∂Bδ (θ) , θ ∈ Θ where Θ is the sample space (see Fig. 7). The governing equilibrium equation for the elastostatic problem of the medium is divσ (θ, Y ) + b = 0 in Bδ (θ)

(3)

where b are the body forces acting on the medium and Y denotes the coordinate system on the mesoscale model Bδ (θ) with the stress and strain fields connected by Hooke’s elasticity tensor C (θ, Y ) σ (θ, Y ) = C (θ, Y ) : ε (θ, Y )

(4)

These fields can be expressed as a superposition of means (¯ σ and ε¯) and of 0 0 zero-mean fluctuations (σ and ε ) as follows: σ (θ, Y ) = σ ¯ + σ 0 (θ, Y ) , ε (θ, Y ) = ε¯ + ε0 (θ, Y )

(5)

The means of stress and strain tensors at some point X of the macrocontinuum (see Fig. 6) are computed as volume averages over Bδ (θ) in the form [19]: Z Z 1 1 σ (θ, Y ) dVδ , ε¯(X) = ε (θ, Y ) dVδ (6) σ ¯ (X) = Vδ Vδ with Vδ =

Z

Bδ (θ)

Bδ (θ)

dVδ .

Bδ (θ)

Also the volume average of the strain energy can be calculated as: 1 U¯ = 2Vδ

Z

1 1 1 ¯ :ε ¯ + σ 0 : ε0 σ (θ, Y ) : ε (θ, Y ) dVδ = σ : ε = σ 2 2 2

(7)

Bδ (θ)

Note that for an unbounded spacedomain (δ → ∞) the fluctuation terms in  0 (7) become negligible σ : ε0 = 0 and thus the following equation holds: ¯ :ε ¯ σ:ε=σ 11

(8)

which is known as Hill’s condition. However, at a finite mesoscale, Hill’s condition is valid provided that the following constraint is satisfied [42]: Z (t − σ ¯ · n) · (u − ε¯ · Y ) dS = 0 (9) ∂Bδ

The constraint equation (9) is a priori satisfied by the following types of boundary conditions: a) uniform strains (kinematic, essential or Dirichlet): u (Y ) = ε¯ · Y , ∀ Y ∈ ∂Bδ

(10)

b) uniform stresses (static, natural or Neumann): t (Y ) = σ ¯ · n , ∀ Y ∈ ∂Bδ

(11)

c) uniform orthogonal-mixed: (t (Y ) − σ ¯ · n) · (u (Y ) − ε¯ · Y ) , ∀ Y ∈ ∂Bδ

(12)

d) periodic:




u Y + − u Y − = ε¯ · Y + − Y − and t Y + + t Y − = 0 , ∀ Y + ∈ ∂Bδ+ and ∀ Y − ∈ ∂Bδ− with ∂Bδ = ∂Bδ+ ∪ ∂Bδ−

(13)

3.2. Weak form and XFEM discrete system Using the principle of total potential energy, the weak form of Eq. (3) is given by Z

Bδ

σ (u) : ε (δu) dVδ =

Z

∂Bδ t

t · δu dS +

Z

Bδ

b · δu dVδ , ∀δu ∈ V

(14)

The virtual displacements δu belong to the functional space V which contains any set of kinematically admissible test functions and is defined as V := {δu ∈ H1 (Bδ ) : δu = 0 on ∂Bδ u } 12

(15)

Figure 7: Discretization of a mesoscale model of a spatially random composite material.

13

where H1 (Bδ ) is the Sobolev space of square integrable functions with square integrable first derivatives on Bδ . The static problem can then be stated as: find u ∈ U := {u ∈ H1 (Bδ ) : u = ¯u on ∂Bδ u } such that ∀δu ∈ V

(16)

a (u, δu) = l (δu) where the bilinear form a (·, ·) and the linear form l (·) are defined as: Z a (u, δu) := σ (u) : ε (δu) dVδ l (δu) :=

Z

Bδ

∂Bδ t

t · δu dS +

Z

Bδ

(17)

b · δu dVδ

In the context of XFEM, weak discontinuities (material interfaces) can be captured by a discontinuous approximation of the displacement function uh (Y ) [43]. This function can be decomposed into a continuous (fem) and a discontinuous (enriched) part as follows: uh (Y ) = uhfem (Y ) + uhenr (Y ) = X

Ni (Y ) ui +

X

Nj (Y )

j∈J

i∈I

n0 X

ψk (Y ) αjk

k=1

!

(18)

where I is the set of all nodes in the mesh, J is the set of nodes (denoted as red circles in Fig. 7) that are enriched with the enrichment functions ψk and αjk are the enriched nodal variables corresponding to node j whose support is cut by the k th inclusion. To improve the accuracy and convergence of XFEM solution, Mo¨es et al. [35] have proposed the following enrichment function: X k X ψk (Y ) = Ni (Y ) φi − Ni (Y ) φki (19) i∈I

i∈I

This ridge function is centered on the interface, has discontinuous first derivative and zero value on the elements which are not crossed by the interface. φki are the nodal values of a level set function corresponding to the k th inclusion. The level set function is used here to represent the circular shape of the inclusions and is expressed in the form: φ (Y ) = kY − ck − r 14

(20)

where Y = (Y1 , Y2 ) is the spatial location of a point in the meshed domain (see Fig. 7) and c is the center of the circular inclusion of radius r. The isozero of Eq. (20) defines the location of the inclusion interface. Substituting Eq. (18) into the weak form of Eq. (14), we get a discrete system of algebraic equations:      Kuu Kuα u F = u (21) Kαu Kαα α Fα Pn0 Pn0 where [Kuu ]ij = a (NP i , Nj ), [Kαα ]ij = a (Ni k=1 ψk , Nj k=1 ψk ) and [Kuα ]ij = 0 [Kαu ]ji = a (Ni , Nj nk=1 ψk ) are the stiffness matrices associated with the standard FE approximation, the enriched approximation and the coupling between them, The forces are expressed as [Fu ]i = l (Ni ) and Prespectively. n0 [Fα ]j = l (Nj k=1 ψk ). From the solution of the system, we finally obtain the nodal displacements u and the enriched variables α. 3.3. Computation of apparent properties on mesoscale Miehe and Koch [10] proposed a computational procedure to define apparent properties (homogenized stresses and overall tangent moduli) of microstructures undergoing small strains. They have shown that apparent properties can be defined in terms of discrete forces and stiffness properties on the boundary of discretized microstructures. Using these deformation-driven algorithms, the apparent stiffness or compliance tensor of a mesoscale model of size δ can be calculated by solving a uniform strain or a uniform stress boundary value problem, respectively. The adopted computational procedure is outlined below for the two cases of uniform boundary conditions. 3.3.1. Uniform strains A prescribed uniform strain tensor ε¯ = [¯ ε11 ε¯22 2¯ ε12 ]T is applied on the boundary ∂Bδ of a discretized mesoscale model, as that of Fig. 7 through displacement boundary conditions (Dirichlet) in the form: ub = DbT ε¯

(22)

where Db is a geometric matrix which depends on the coordinates of the boundary node b and is defined as:   2Y1 0 1 (23) Db =  0 2Y2  with (Y1 , Y2 ) ∈ Y 2 Y2 Y1 15

Note that the stiffness matrix K of the model (see eq. (21)) can be rearranged into sub-matrices associated with interior nodes i and boundary nodes b. Thus the static problem is denoted by:      ui Fui Kui ui Kui αi Kui ub Kui αb Kαi ui Kαi αi Kαi u Kαi α   αi  Fαi  b b      (24) Kub ui Kub αi Kub ub Kub αb   ub  = Fub  Fαb Kαb ui Kαb αi Kαb ub Kαb αb αb where by setting

     Kub ub Kub αb Kui ub Kui αb Kui ui Kui αi , Kbb = , Kib = (25) Kii = Kαi ui Kαi αi Kαb ub Kαb αb Kαi ub Kαi αb 

        ui ub Fui F Ui = , Ub = , Fi = , Fb = ub αi αb F αi Fαb

it can be simplified as follows:      Kii Kib Ui F = i Kbi Kbb Ub Fb

(26)

(27)

Then the apparent stiffness tensor CδD (θ) of the mesoscale model of size δ under Dirichlet boundary conditions can be calculated in terms of the ˜ bb in the form: condensed stiffness matrix K CδD (θ) =

1 ˜ T DKbb D Vδ

(28)

˜ bb = Kbb − Kbi K −1 Kib and D = [D1 D2 ... DM ] with M the total where K ii number of boundary nodes b. 3.3.2. Uniform stresses A prescribed uniform stress tensor σ ¯ = [¯ σ11 σ ¯22 2¯ σ12 ]T is applied on the boundary surface ∂Bδ of a discretized mesoscale model, as that of Fig. 7 through traction boundary conditions (Neumann) as follows: Fb = SbT σ ¯

16

(29)

where Sb is a matrix depending on the components of the discrete area vector ab which is given in terms of the nodal coordinates of the neighboring boundary nodes b − 1, b and b + 1 in the form: 1 [Yb+1 − Yb−1 ] × en (30) 2 These nodes are oriented so that the cross product with the Cartesian, out-ofplane, base vector en yields ab as an outward normal vector at the boundary node b. Thus the matrix Sb is defined as:   ab1 0 Sb =  0 a b2  (31) ab2 ab1 ab =

Then the apparent compliance tensor SδN (θ) of the mesoscale model of size δ under Neumann boundary conditions can be calculated in terms of the ˜ bb in the form: condensed stiffness matrix K SδN (θ) =

1 ˜ −1 T SKbb S Vδ

(32)

where S = [S1 S2 ... SM ] with M the total number of boundary nodes b. For this type of boundary conditions the apparent stiffness tensor is obtained by inverting the compliance tensor of Eq. (32):  −1 CδN (θ) = SδN (θ)

(33)

Dirichlet and Neumann boundary conditions provide upper and lower bounds of the strain energy which converge to each other as the mesoscale size δ is increasing. Thus the following relation holds [44]:  1   1   ε¯ : CδD (θ) : ε¯ for δ finite   ε¯ : CδN (θ) : ε¯ < 2 2 (34)    1 ε¯ : C N (θ) : ε¯ = 1 ε¯ : C D (θ) : ε¯ for δ → ∞ ∞ ∞ 2 2   In case of uniaxial strains ε¯ = [1 0 0]T or [0 1 0]T or simple shear ε¯ = [0 0 1]T the following notation can be used:

CδN (θ) ≤ CδD (θ) 17

(35)

4. Results and discussion In this section, the spatial average (mean) of the apparent moduli is calculated at each size δ from nw mesoscale models extracted by using the window technique on the composite images of section 2. Thus the following formula for the mean is used: nw   1 X Cδ (θ) = Cδ ξ~p , θ nw p=1

(36)

with ξ~p the position vector of the p-th mesoscale window model on the image (see Fig. 2). Following the notation of Eq. (35) the bounds CδN (θ) and CδD (θ) are calculated for the axial and shear components of the apparent elasticity tensor. Note that, while the constituent materials of the composite are considered isotropic, the computational method of section 3.3 results in an anisotropic apparent elasticity tensor due to spatial randomness. Composites containing either stiff or compliant inclusions are investigated which are obtained by keeping the Young’s modulus of the matrix constant at Em = 1 GPa and varying the Young’s modulus of the inclusions Ein . The Poisson ratio is assumed to be the same for the two phases, νm = νin = 0.3. All components of the apparent elasticity tensor computed in this section are in GPa. 4.1. Comparison with analytical bounds The convergence of CδN (θ) and CδD (θ) with respect to the non-dimensional parameter δ is illustrated in Fig. 8 for stiff (Ein /Em =10) and compliant (Ein /Em =0.1) inclusions. The results refer to the average axial stiffness Cii /2 (i = 1, 2) and shear stiffness C33 ≡ G components of the mean elasticity tensor. Moreover, analytical bounds such as those of Voigt and Reuss CδV (θ) = fm (θ) · Cm + fin (θ) · Cin CδR



1 1 + fin (θ) · (θ) = fm (θ) · Cm Cin

−1

as well as those of Hashin-Shtrikman  −1 1 1 Hu Cδ (θ) = Cm + fin (θ) · + fm (θ) · Cin − Cm 2Cm 18

(37) (38)

(39)

CδH`



1 1 + fin (θ) · (θ) = Cin + fm (θ) · Cm − Cin 2Cin

−1

(40)

have been calculated on each mesoscale model of size δ and the corresponding spatially averaged results are plotted in the same figure for comparison purposes. In Eqs. (37)-(40) the symbol f (θ) stands for the random volume fraction, with fm (θ) + fin (θ) = 1. Note that Eqs. (37)-(40) refer to bulk and shear modulus as the corresponding elasticity tensors are isotropic. On the other hand the computational method of section 3.3 can calculate anisotropic elasticity tensors of the form:   C11 C12 C13 C22 C23  Cδ (θ) =  (41) sym C33

It can be observed from Fig. 8 that for all mesoscale sizes δ in case of stiff inclusions and for δ > 5 in case of compliant inclusions the following inequalities hold: CδR (θ) ≤ CδH` (θ) ≤ CδN (θ) ≤ CδD (θ) ≤ CδHu (θ) ≤ CδV (θ)

(42)

This means that much more strict bounds are calculated through numerical analysis. The lower Hashin-Shtrikman bound is larger than the lower numerical bound only in the case of compliant inclusions for small windows (δ ≤ 5). However, even in this case, the numerical bounds are inside the analytical bounds of Voigt and Reuss, which means that the proposed computational scheme gives correct results. 4.2. Effect of volume fraction and stiffness ratio on RVE size In Fig. 9, results related to the determination of RVE size are presented for the composite materials of images 1 and 2 (see Fig. 1). As mentioned before in section 2, image 1 corresponds to a composite with initial vf=40%, whereas image 2 corresponds to a composite with initial vf=20%. Note that both composites contain circular inclusions of constant size (d ' 7.14 µm). The effect of volume fraction on the mechanical properties of the composite is also illustrated in Fig. 9. It can be observed that the convergence rate of the two numerical bounds with respect to δ is slower as vf increases. This means that, in order to perform homogenization in composites with high volume fraction, a much larger RVE has to be used. 19

Figure 8: Analytical and numerical bounds of apparent properties with respect to mesoscale size δ for the composite of image 1 with (a-b) stiff (Ein /Em = 10) and (cd) compliant (Ein /Em = 0.1) inclusions.

Figure 9: Comparison of numerical bounds of apparent properties for the composites of images 1 and 2 containing stiff inclusions (Ein /Em = 10).

20

The results of a parametric study with respect to matrix/inclusion stiffness ratio are shown in Figs. 10 and 11 for the composite of image 1. The cases of stiff (Ein /Em =10, 100, 1000) and compliant (Ein /Em =0.1, 0.01, 0.001) inclusions have been tested. As the Young’s modulus of the inclusions increases, the computed apparent properties of the composite increase as well. The effect of stiffness ratio on the convergence rate of the apparent properties is also illustrated in these figures. This effect is better illustrated in Fig. 12 where the discrepancy (tolerance) eδ between the numerical bounds is plotted against the non-dimensional parameter δ for both composites examined (images 1 and 2). The discrepancy is calculated as: CD − CN ij ij (43) eδ = CijD δ

CijD

CijN

where and are the components of the spatially averaged apparent elasticity tensor of Eq. (41) obtained by applying Dirichlet and Neumann boundary conditions, respectively. As it can be observed from Fig. 12, in order to attain the RVE within a tolerance of 20%, a mesoscale window of size δ > 2.8 is required for composites containing stiff inclusions with stiffness ratio equal to 10, while for stiffness ratio equal to 100, a size δ > 18 is needed. For stiffness ratio equal to 1000, it is not possible to attain the RVE within a tolerance of 20%. In this case, the convergence rate is very slow and thus a very large value of δ is needed, which practically corresponds to the full microstructural analysis of the composite. From Fig. 12, it can also be observed that the compliant inclusions lead to slower convergence compared to the corresponding stiff ones. Another important observation is that the RVE size has to be determined with respect to the convergence of an individual apparent property. Note that the convergence of the apparent shear modulus C33 is much faster than the convergence of the average apparent axial stiffness Cii /2 (i=1, 2). 4.3. Determination of mesoscale random fields In section 4.2 it was shown that the RVE can not always be attained within an acceptable tolerance. It was also shown that the RVE size mainly depends on the volume fraction, the stiffness ratio between the constituent materials and the specific apparent property, the convergence of which is studied with respect to the boundary conditions. In the context of stochastic multi-scale analysis of composites with the stochastic finite element method, 21

Figure 10: Convergence of numerical bounds on the average axial stiffness (Cii /2) with respect to mesoscale size δ for stiffness ratio a) Ein /Em = 10, b) Ein /Em = 100, c) Ein /Em = 1000, d) Ein /Em = 0.1, e) Ein /Em = 0.01, f) Ein /Em = 0.001.

Figure 11: Convergence of numerical bounds on the shear stiffness (C33 ) with respect to mesoscale size δ for stiffness ratio a) Ein /Em = 10, b) Ein /Em = 100, c) Ein /Em = 1000, d) Ein /Em = 0.1, e) Ein /Em = 0.01, f) Ein /Em = 0.001.

22

Figure 12: Discrepancy between upper and lower bounds with respect to mesoscale size δ for the composites of image 1 (a-b) and image 2 (c-d) and various cases of stiffness ratio.

23

the material properties provided by the up-scaling procedure appear either as random fields (apparent properties in case of SVE) or as random variables (effective properties in case of RVE) [45]. In this section, random fields for the volume fraction and the apparent elasticity tensor of the composite of image 1 (see section 2) are obtained in mesoscale sizes δ from the simulation of a large number of SVE models extracted using the moving window technique. In order to obtain accurate statistical properties, the moving distance step of the method has been chosen as ∆ξ = 0.25L, with L the moving window size. Specifically, random fields for δ =5.605, 11.210 and 22.420 have been derived by simulating nw =3721, 841 and 169 SVEs, respectively (see Eq. (1) in section 2 for nw ). All the presented results refer to kinematic uniform boundary conditions. Figs. 13, 14 and 15 depict the computed random fields along with the respective empirical distributions and 2-D spatial correlations of the volume fraction and the average axial and shear components of the apparent elasticity tensor, for the three mesoscale sizes δ examined. Note that the spatial correlations ρB A have been calculated for every lag (ξx , ξy ) according to the following formula: √

√

 nw nw  ¯ B (xi + ξx , yj + ξy ) − B 1 X X A (xi , yj ) − A¯ B ρA (ξx , ξy ) = nw − 1 i=0 j=0 σA σB √ √ √ √ − nw ∆ξ ≤ ξx ≤ nw ∆ξ , − nw ∆ξ ≤ ξy ≤ nw ∆ξ

(44) with denoting auto-correlations when quantity A ≡ B, otherwise cross¯ B ¯ are the spatial average values while σA , σB are correlations are defined. A, the standard deviations of quantities A, B, respectively. The mesoscale random fields of the components of the apparent elasticity tensor in Figs. 14 and 15 have been derived for stiffness ratio Ein /Em = 10. A general observation for all mesoscale random fields is that their empirical PDFs become narrower as the mesoscale size δ increases. In other words, the random field tends to a random variable and thus the SVE tends to the RVE as δ increases. The auto-correlations for lag (ξx = 0, ξy = 0) are 1 and tend to zero for lag values |ξx | > L and |ξy | > L. 0.5(C +C ) 33 The 1-D cross correlations ρνf 11 22 and ρC νf are depicted in Figs. 16 and 17, respectively, for lag values along x (ξx , ξy = 0) and y (ξx = 0, ξy ) for stiffness ratio Ein /Em =10 and 1000. All cross-correlations are between ρB A

24

Figure 13: Mesoscale random fields for volume fraction.

Figure 14: Mesoscale random fields for average axial stiffness (Ein /Em = 10).

25

Figure 15: Mesoscale random fields for shear stiffness (Ein /Em = 10).

-1 and +1 and tend to zero for |ξx | > L or |ξy | > L. This can be attributed to the fact that the probability of SVE models sharing common inclusions decreases as the length of the vector ξ~p (see Fig. 2) increases. A decrease on the values of cross-correlations with the increase of stiffness ratio is also observed, as in Sena et al. [46], especially for |ξx | ≤ L and |ξy | ≤ L for all sizes δ. 0.5(C +C ) This is more obvious in Fig. 18 where the cross-correlations ρνf 11 22 33 and ρC νf at zero lag (ξx = 0, ξy = 0) are plotted with respect to δ for Ein /Em =10 and 1000. This could be an explanation of the fact that the RVE size can not be defined within a reasonable tolerance for the case of high stiffness ratio (see results of section 4.2). In addition, the effect of mesoscale size δ on the examined cross-correlations is visible for large stiffness ratio. In other words, this means that the effect of volume fraction variability on the apparent properties as δ decreases is more important in case of large stiffness ratio. The cross-correlations between the axial stiffness components of the ap22 parent elasticity tensor ρC C11 are illustrated in Fig. 19 for lag values along x (ξx , ξy = 0) and y (ξx = 0, ξy ) for stiffness ratio Ein /Em =10 and 1000. Conclusions similar to those derived previously from Figs. 16 and 17 can be deduced. Fig. 20 illustrates the effect of stiffness ratio and mesoscale size δ

26

Figure 16: Cross correlations at (ξx , ξy = 0) and (ξx = 0, ξy ) for Ein /Em = 10 and Ein /Em = 1000.

27

Figure 17: Cross correlations at (ξx , ξy = 0) and (ξx = 0, ξy ) for Ein /Em = 10 and Ein /Em = 1000.

Figure 18: Cross correlations at (ξx = 0, ξy = 0) for Ein /Em = 10 and Ein /Em = 1000.

28

Figure 19: Cross correlations at (ξx , ξy = 0) and (ξx = 0, ξy ) for Ein /Em = 10 and Ein /Em = 1000.

on the level of anisotropy in the apparent elasticity tensor. It can be observed C22 that the cross-correlation ρC at (ξx = 0, ξy = 0) decreases with increasing 11 22 stiffness ratio and the anisotropy indicator eC C11 = mean(|C11 − C22 |/2C22 ) decreases with increasing δ. This last observation is in accordance with the assumption that, at the RVE level, the material can be considered isotropic in the mean sense.

29

Figure 20: Effect of stiffness ratio Ein /Em and macroscale size δ on the degree of anisotropy of the apparent elasticity tensor.

5. Conclusions The determination of RVE size for composite materials constitutes the basis of homogenization methods. In this paper, a novel computational procedure based on XFEM and MCS has been proposed for this purpose. The proposed approach takes into account the local volume fraction variation by processing computer-simulated images of composites with randomly dispersed inclusions. A variable number of microstructure models were derived directly from the images using the moving window technique. The XFEM analysis of the microstructure models was used to compute upper and lower bounds on the apparent material properties by solving multiple boundary value problems under kinematic and static uniform boundary conditions. The RVE was attained within a prescribed tolerance by examining the convergence of these two bounds with respect to the mesoscale size. The effect of matrix/inclusion stiffness ratio as well as of the volume fraction of inclusions on the RVE size was highlighted. It has been shown that the convergence rate is slower in the case of compliant inclusions and high volume fraction, and that the variability of the apparent properties is affected by the stiffness ratio and mesoscale size. Mesoscale random fields describing the spatial variation of the components of the apparent elasticity tensor were also determined. The proposed computational procedure can provide with a more realistic 30

RVE because it incorporates the local volume fraction variation present in real composites and can be extended to any filler-reinforced composite as it is based on image analysis of computer-simulated microstructures. Finally, the obtained results can be used in the framework of homogenization methods and stochastic finite element analysis of composite structures. Acknowledgements The financial support provided by the European Research Council Advanced Grant ”MASTER - Mastering the computational challenges in numerical modeling and optimum design of CNT reinforced composites” (ERC2011-ADG 20110209) is gratefully acknowledged by the authors. References [1] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behaviour of multiphase materials, Journal of the Mechanics and Physics of Solids 11 (2) (1963) 127–140. [2] R. Hill, A self-consistent mechanics of composite materials, Journal of the Mechanics and Physics of Solids 13 (4) (1965) 213–222. [3] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta metallurgica 21 (5) (1973) 571–574. [4] S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, New York, 2002. [5] T. I. Zohdi, P. Wriggers, An Introduction to Computational Micromechanics, 2nd Edition, Lecture Notes in Applied and Computational Mechanics, vol. 20, Springer, Heidelberg, 2008. [6] J. Ma, I. Temizer, P. Wriggers, Random homogenization analysis in linear elasticity based on analytical bounds and estimates, International Journal of Solids and Structures 48 (2) (2011) 280–291. [7] F. Feyel, J.-L. Chaboche, FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials, Computer methods in applied mechanics and engineering 183 (3) (2000) 309–330. 31

[8] K. Terada, N. Kikuchi, A class of general algorithms for multi-scale analyses of heterogeneous media, Computer methods in applied mechanics and engineering 190 (40) (2001) 5427–5464. [9] S. Baxter, M. I. Hossain, L. Graham, Micromechanics based random material property fields for particulate reinforced composites, International journal of solids and structures 38 (50) (2001) 9209–9220. [10] C. Miehe, A. Koch, Computational micro-to-macro transitions of discretized microstructures undergoing small strains, Archive of Applied Mechanics 72 (4-5) (2002) 300–317. [11] Z. Yuan, J. Fish, Toward realization of computational homogenization in practice, International Journal for Numerical Methods in Engineering 73 (3) (2008) 361–380. [12] N. Charalambakis, Homogenization techniques and micromechanics. A survey and perspectives, Applied Mechanics Reviews 63 (3) (2010) 803. [13] M. Geers, V. Kouznetsova, W. Brekelmans, Multi-scale computational homogenization: Trends and challenges, Journal of Computational and Applied Mathematics 234 (7) (2010) 2175–2182. [14] J. E. Bishop, J. M. Emery, R. V. Field, C. R. Weinberger, D. J. Littlewood, Direct numerical simulations in solid mechanics for understanding the macroscale effects of microscale material variability, Computer Methods in Applied Mechanics and Engineering 287 (2015) 262–289. [15] S. Baxter, L. Graham, Characterization of Random Composites Using Moving-Window Technique, Journal of Engineering Mechanics 126 (4) (2000) 389–397. [16] X. F. Xu, L. Graham-Brady, A stochastic computational method for evaluation of global and local behavior of random elastic media, Computer Methods in Applied Mechanics and Engineering 194 (42) (2005) 4362–4385. [17] D. Charmpis, G. Schu¨eller, M. Pellissetti, The need for linking micromechanics of materials with stochastic finite elements: A challenge for materials science, Computational Materials Science 41 (1) (2007) 27–37. 32

[18] A. Cl´ement, C. Soize, J. Yvonnet, Uncertainty quantification in computational stochastic multiscale analysis of nonlinear elastic materials, Computer Methods in Applied Mechanics and Engineering 254 (2013) 61–82, ISSN 0045-7825. [19] R. Hill, Elastic properties of reinforced solids: some theoretical principles, Journal of the Mechanics and Physics of Solids 11 (5) (1963) 357–372. [20] A. Bensoussan, J.-L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, vol. 374, American Mathematical Soc., 2011. [21] W. E, B. Engquist, X. Li, W. Ren, E. Vanden-Eijnden, Heterogeneous multiscale methods: a review, Commun. Comput. Phys 2 (3) (2007) 367–450. [22] M. Ostoja-Starzewski, Material spatial randomness: From statistical to representative volume element, Probabilistic engineering mechanics 21 (2) (2006) 112–132. [23] T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: statistical and numerical approach, International Journal of solids and structures 40 (13) (2003) 3647–3679. [24] J. Zeman, M. Sejnoha, From random microstructures to representative volume elements, Modelling and Simulation in Materials Science and Engineering 15 (4) (2007) 325–335. [25] M. Silani, H. Talebi, S. Ziaei-Rad, P. Kerfriden, S. P. Bordas, T. Rabczuk, Stochastic modelling of clay/epoxy nanocomposites, Composite Structures 118 (2014) 241–249. [26] J. Wimmer, B. Stier, J.-W. Simon, S. Reese, Computational homogenisation from a 3D finite element model of asphalt concrete–linear elastic computations, Finite Elements in Analysis and Design 110 (2016) 43–57. [27] D. Trias, J. Costa, A. Turon, J. Hurtado, Determination of the critical size of a statistical representative volume element (SRVE) for carbon reinforced polymers, Acta Materialia 54 (13) (2006) 3471–3484. 33

[28] I. Gitman, H. Askes, L. Sluys, Representative volume: existence and size determination, Engineering fracture mechanics 74 (16) (2007) 2518– 2534. [29] X. F. Xu, X. Chen, Stochastic homogenization of random elastic multiphase composites and size quantification of representative volume element, Mechanics of Materials 41 (2) (2009) 174–186. [30] J. Ma, J. Zhang, L. Li, P. Wriggers, S. Sahraee, Random homogenization analysis for heterogeneous materials with full randomness and correlation in microstructure based on finite element method and Monte-carlo method, Computational Mechanics 54 (6) (2014) 1395–1414. [31] J. Ma, S. Sahraee, P. Wriggers, L. De Lorenzis, Stochastic multiscale homogenization analysis of heterogeneous materials under finite deformations with full uncertainty in the microstructure, Computational Mechanics 55 (2015) 819–835. [32] M. Ostoja-Starzewski, X. Wang, Stochastic finite elements as a bridge between random material microstructure and global response, Computer Methods in Applied Mechanics and Engineering 168 (1) (1999) 35–49. [33] V. Lucas, J.-C. Golinval, S. Paquay, V.-D. Nguyen, L. Noels, L. Wu, A stochastic computational multiscale approach; Application to MEMS resonators, Computer Methods in Applied Mechanics and Engineering 294 (2015) 141–167. [34] G. Stefanou, D. Savvas, M. Papadrakakis, Stochastic finite element analysis of composite structures based on material microstructure, Composite Structures 132 (2015) 384–392. [35] N. Mo¨es, M. Cloirec, P. Cartraud, J.-F. Remacle, A computational approach to handle complex microstructure geometries, Computer Methods in Applied Mechanics and Engineering 192 (28) (2003) 3163–3177. [36] B. Hiriyur, H. Waisman, G. Deodatis, Uncertainty quantification in homogenization of heterogeneous microstructures modeled by XFEM, International Journal for Numerical Methods in Engineering 88 (3) (2011) 257–278.

34

[37] D. Savvas, G. Stefanou, M. Papadrakakis, G. Deodatis, Homogenization of random heterogeneous media with inclusions of arbitrary shape modeled by XFEM, Computational Mechanics 54 (5) (2014) 1221–1235. [38] J. M. Melenk, I. Babuˇska, The partition of unity finite element method: basic theory and applications, Computer Methods in Applied Mechanics and Engineering 139 (1) (1996) 289–314. [39] Z. Shan, A. M. Gokhale, Representative volume element for non-uniform micro-structure, Computational Materials Science 24 (3) (2002) 361– 379. [40] S. H. R. Sanei, R. S. Fertig, Uncorrelated volume element for stochastic modeling of microstructures based on local fiber volume fraction variation, Composites Science and Technology 117 (2015) 191–198. [41] Z. Hashin, Analysis of composite materials: A survey, Journal of Applied Mechanics 50 (2) (1983) 481–505. [42] S. Hazanov, C. Huet, Order relationships for boundary conditions effect in heterogeneous bodies smaller than the representative volume, Journal of the Mechanics and Physics of Solids 42 (12) (1994) 1995–2011. [43] T. Belytschko, R. Gracie, G. Ventura, A review of extended/generalized finite element methods for material modeling, Modelling and Simulation in Materials Science and Engineering 17 (4) (2009) 1–24. [44] P. M. Suquet, Elements of homogenization for inelastic solid mechanics, in: E. Sanchez-Palencia and A. Zaoui, Eds., Homogenization Techniques for Composite Media, Springer, 193–278, 1987. [45] A. Kucerov´a, H. Matthies, Uncertainty updating in the description of heterogeneous materials, Technische Mechanik 30 (1-3) (2010) 211–226. [46] M. P. Sena, M. Ostoja-Starzewski, L. Costa, Stiffness tensor random fields through upscaling of planar random materials, Probabilistic Engineering Mechanics 34 (2013) 131–156. [1] Z. Hashin, S. Shtrikman, A variational approach to the theory of the elastic behaviour of multiphase materials, Journal of the Mechanics and Physics of Solids 11 (2) (1963) 127–140. 35

[2] R. Hill, A self-consistent mechanics of composite materials, Journal of the Mechanics and Physics of Solids 13 (4) (1965) 213–222. [3] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta metallurgica 21 (5) (1973) 571–574. [4] S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, New York, 2002. [5] T. I. Zohdi, P. Wriggers, An Introduction to Computational Micromechanics, 2nd Edition, Lecture Notes in Applied and Computational Mechanics, vol. 20, Springer, Heidelberg, 2008. [6] J. Ma, I. Temizer, P. Wriggers, Random homogenization analysis in linear elasticity based on analytical bounds and estimates, International Journal of Solids and Structures 48 (2) (2011) 280–291. [7] F. Feyel, J.-L. Chaboche, FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials, Computer methods in applied mechanics and engineering 183 (3) (2000) 309–330. [8] K. Terada, N. Kikuchi, A class of general algorithms for multi-scale analyses of heterogeneous media, Computer methods in applied mechanics and engineering 190 (40) (2001) 5427–5464. [9] S. Baxter, M. I. Hossain, L. Graham, Micromechanics based random material property fields for particulate reinforced composites, International journal of solids and structures 38 (50) (2001) 9209–9220. [10] C. Miehe, A. Koch, Computational micro-to-macro transitions of discretized microstructures undergoing small strains, Archive of Applied Mechanics 72 (4-5) (2002) 300–317. [11] Z. Yuan, J. Fish, Toward realization of computational homogenization in practice, International Journal for Numerical Methods in Engineering 73 (3) (2008) 361–380. [12] N. Charalambakis, Homogenization techniques and micromechanics. A survey and perspectives, Applied Mechanics Reviews 63 (3) (2010) 803. 36

[13] M. Geers, V. Kouznetsova, W. Brekelmans, Multi-scale computational homogenization: Trends and challenges, Journal of Computational and Applied Mathematics 234 (7) (2010) 2175–2182. [14] J. E. Bishop, J. M. Emery, R. V. Field, C. R. Weinberger, D. J. Littlewood, Direct numerical simulations in solid mechanics for understanding the macroscale effects of microscale material variability, Computer Methods in Applied Mechanics and Engineering 287 (2015) 262–289. [15] S. Baxter, L. Graham, Characterization of Random Composites Using Moving-Window Technique, Journal of Engineering Mechanics 126 (4) (2000) 389–397. [16] X. F. Xu, L. Graham-Brady, A stochastic computational method for evaluation of global and local behavior of random elastic media, Computer Methods in Applied Mechanics and Engineering 194 (42) (2005) 4362–4385. [17] D. Charmpis, G. Schu¨eller, M. Pellissetti, The need for linking micromechanics of materials with stochastic finite elements: A challenge for materials science, Computational Materials Science 41 (1) (2007) 27–37. [18] A. Cl´ement, C. Soize, J. Yvonnet, Uncertainty quantification in computational stochastic multiscale analysis of nonlinear elastic materials, Computer Methods in Applied Mechanics and Engineering 254 (2013) 61–82, ISSN 0045-7825. [19] R. Hill, Elastic properties of reinforced solids: some theoretical principles, Journal of the Mechanics and Physics of Solids 11 (5) (1963) 357–372. [20] A. Bensoussan, J.-L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, vol. 374, American Mathematical Soc., 2011. [21] W. E, B. Engquist, X. Li, W. Ren, E. Vanden-Eijnden, Heterogeneous multiscale methods: a review, Commun. Comput. Phys 2 (3) (2007) 367–450. [22] M. Ostoja-Starzewski, Material spatial randomness: From statistical to representative volume element, Probabilistic engineering mechanics 21 (2) (2006) 112–132. 37

[23] T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: statistical and numerical approach, International Journal of solids and structures 40 (13) (2003) 3647–3679. [24] J. Zeman, M. Sejnoha, From random microstructures to representative volume elements, Modelling and Simulation in Materials Science and Engineering 15 (4) (2007) 325–335. [25] M. Silani, H. Talebi, S. Ziaei-Rad, P. Kerfriden, S. P. Bordas, T. Rabczuk, Stochastic modelling of clay/epoxy nanocomposites, Composite Structures 118 (2014) 241–249. [26] J. Wimmer, B. Stier, J.-W. Simon, S. Reese, Computational homogenisation from a 3D finite element model of asphalt concrete–linear elastic computations, Finite Elements in Analysis and Design 110 (2016) 43–57. [27] D. Trias, J. Costa, A. Turon, J. Hurtado, Determination of the critical size of a statistical representative volume element (SRVE) for carbon reinforced polymers, Acta Materialia 54 (13) (2006) 3471–3484. [28] I. Gitman, H. Askes, L. Sluys, Representative volume: existence and size determination, Engineering fracture mechanics 74 (16) (2007) 2518– 2534. [29] X. F. Xu, X. Chen, Stochastic homogenization of random elastic multiphase composites and size quantification of representative volume element, Mechanics of Materials 41 (2) (2009) 174–186. [30] J. Ma, J. Zhang, L. Li, P. Wriggers, S. Sahraee, Random homogenization analysis for heterogeneous materials with full randomness and correlation in microstructure based on finite element method and Monte-carlo method, Computational Mechanics 54 (6) (2014) 1395–1414. [31] J. Ma, S. Sahraee, P. Wriggers, L. De Lorenzis, Stochastic multiscale homogenization analysis of heterogeneous materials under finite deformations with full uncertainty in the microstructure, Computational Mechanics 55 (2015) 819–835.

38

[32] M. Ostoja-Starzewski, X. Wang, Stochastic finite elements as a bridge between random material microstructure and global response, Computer Methods in Applied Mechanics and Engineering 168 (1) (1999) 35–49. [33] V. Lucas, J.-C. Golinval, S. Paquay, V.-D. Nguyen, L. Noels, L. Wu, A stochastic computational multiscale approach; Application to MEMS resonators, Computer Methods in Applied Mechanics and Engineering 294 (2015) 141–167. [34] G. Stefanou, D. Savvas, M. Papadrakakis, Stochastic finite element analysis of composite structures based on material microstructure, Composite Structures 132 (2015) 384–392. [35] N. Mo¨es, M. Cloirec, P. Cartraud, J.-F. Remacle, A computational approach to handle complex microstructure geometries, Computer Methods in Applied Mechanics and Engineering 192 (28) (2003) 3163–3177. [36] B. Hiriyur, H. Waisman, G. Deodatis, Uncertainty quantification in homogenization of heterogeneous microstructures modeled by XFEM, International Journal for Numerical Methods in Engineering 88 (3) (2011) 257–278. [37] D. Savvas, G. Stefanou, M. Papadrakakis, G. Deodatis, Homogenization of random heterogeneous media with inclusions of arbitrary shape modeled by XFEM, Computational Mechanics 54 (5) (2014) 1221–1235. [38] J. M. Melenk, I. Babuˇska, The partition of unity finite element method: basic theory and applications, Computer Methods in Applied Mechanics and Engineering 139 (1) (1996) 289–314. [39] Z. Shan, A. M. Gokhale, Representative volume element for non-uniform micro-structure, Computational Materials Science 24 (3) (2002) 361– 379. [40] S. H. R. Sanei, R. S. Fertig, Uncorrelated volume element for stochastic modeling of microstructures based on local fiber volume fraction variation, Composites Science and Technology 117 (2015) 191–198. [41] Z. Hashin, Analysis of composite materials: A survey, Journal of Applied Mechanics 50 (2) (1983) 481–505. 39

[42] S. Hazanov, C. Huet, Order relationships for boundary conditions effect in heterogeneous bodies smaller than the representative volume, Journal of the Mechanics and Physics of Solids 42 (12) (1994) 1995–2011. [43] T. Belytschko, R. Gracie, G. Ventura, A review of extended/generalized finite element methods for material modeling, Modelling and Simulation in Materials Science and Engineering 17 (4) (2009) 1–24. [44] P. M. Suquet, Elements of homogenization for inelastic solid mechanics, in: E. Sanchez-Palencia and A. Zaoui, Eds., Homogenization Techniques for Composite Media, Springer, 193–278, 1987. [45] A. Kucerov´a, H. Matthies, Uncertainty updating in the description of heterogeneous materials, Technische Mechanik 30 (1-3) (2010) 211–226. [46] M. P. Sena, M. Ostoja-Starzewski, L. Costa, Stiffness tensor random fields through upscaling of planar random materials, Probabilistic Engineering Mechanics 34 (2013) 131–156.

40