A random field-based method to estimate convergence of apparent properties in computational homogenization

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270 www.elsevier.com/locate/cma

A random field-based method to estimate convergence of apparent properties in computational homogenization Kirubel Teferra a , ∗, Lori Graham-Brady b a Multifunctional Materials Branch Code 6350, U.S. Naval Research Laboratory, Washington, DC, USA b Department of Civil Engineering, Johns Hopkins University, MD, USA

Received 24 July 2017; received in revised form 29 September 2017; accepted 26 October 2017 Available online 8 November 2017

Highlights • • • • •

Convergence rate of apparent properties in computational homogenization is computed. Rate is determined from analysis of only one statistical volume element size. Higher order terms affecting variance of apparent properties are derived. Omission of higher order terms leads to extrapolation errors in predicting RVE size. Methodology is demonstrated for a 3D crystal plasticity finite element model.

Abstract This work addresses the challenge of efficiently determining the representative volume size in computational homogenization by exploiting the requirement that mechanical response quantities are statistically homogeneous and ergodic random fields. The proposed computational homogenization approach focuses on empirically determining the autocorrelation functions of output response quantities through post processing finite element analyses. Once the autocorrelation functions are known, only trivial computations are needed to determine the variance of apparent properties as a function of domain size without any additional finite element computation by utilizing a simple formula relating the autocorrelation function to the variance of spatially averaged quantities. This approach improves upon the current established approach by circumventing the need to analyze numerous successively larger domains in order to determine convergence of apparent properties, which can be computationally prohibitive. Furthermore, similar previous analytical expressions for the variance of apparent properties are asymptotic with respect to domain size while the expression in the proposed approach is exact. After presenting the formulation, the method is first demonstrated for a one-dimensional bar model where analytical expressions can be derived. Then the approach is demonstrated on two numerical examples: (1) a plane strain, linear elastic domain with a lognormal random field describing its compliance, and (2) a stochastic polycrystalline microstructure of a Nickel super alloy having a crystal plasticity model for its constitutive law. Published by Elsevier B.V.

Keywords: Computational homogenization; Crystal plasticity; Statistical volume element; Representative volume element; Stochastic simulation

∗ Corresponding author.

E-mail address: [email protected] (K. Teferra). https://doi.org/10.1016/j.cma.2017.10.027 0045-7825/Published by Elsevier B.V.

254

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

1. Introduction Computing and post processing the mesoscale mechanical response through the finite element solution of a boundary value problem (BVP) formulated on the scale of material heterogeneity is of great interest in the areas of computational homogenization, multiscale methods, and, more recently, materials by design. The primary goal of these methods is to relate the effects of material heterogeneity, e.g., polycrystalline morphology, crystallography, and constituents, to macroscopic material properties. For computational homogenization and multiscale methods, mesoscale mechanical response quantities are spatially averaged to obtain macroscopic properties. On the other hand, materials design is primarily focused on relating the effects of underlying deformation mechanisms on mesoscale mechanical response in order to understand how microstructural features give materials their idiosyncratic properties (i.e., microstructure–property correlations [1]). The mesoscale, as defined in this work, is the smallest length scale that encompasses material heterogeneity with statistical significance, which, for polycrystalline metals, is typically 100 µm to 1000 µm, containing on the order of thousands of grains. For this reason, the mesoscale will be denoted as the microscale for the remainder of this paper. This definition of the microscale, often described as the representative volume element (RVE), leads to an ambiguous definition of its size, which is strongly sensitive to the particular mechanical response quantity of interest. In response, the concept of the statistical volume element (SVE) was introduced by recognizing that spatially averaged response quantities, denoted as apparent properties, are random variables that converge to deterministic quantities asymptotically with respect to the size of the spatial domain of the microscale BVP [2,3]. The random apparent properties are a direct consequence of the stochastic nature of microstructure configurations, notably morphology, crystallography, and constitution. There is extensive literature in recent years on sampling based methods to compute apparent properties as a function of SVE size, e.g., [4–28]. These methods follow the general framework of simulating numerous statistically equivalent microstructures, solving a BVP under admissible set of BCs satisfying the Hill–Mandel condition, evaluating statistics of the apparent properties, and assessing the number of microstructural features required to obtain convergence of apparent properties to effective properties within some tolerance. For example, in [5] a thorough study was done computing variability of Von-Mises yield criteria through Monte Carlo (MC) simulation for ductile, porous metals using J2 -type plasticity for identically-sized spherical pores of varying number and volume fraction. One of the most commonly recommended procedures to determine RVE size is to observe the decrease in variation of apparent properties by performing MC simulation for increasing SVE sizes [9,12,29]. In [12], an empirical power law form for the variance of apparent properties as a function of SVE size was introduced whose parameters are fit via the series of numerical experiments described above for specifically chosen SVE sizes. This approach, and similar variants (e.g., [9]), are considered to be the only alternative to simply performing MC simulation of the BVP on increasing size SVEs until convergence is obtained. However, this approach is still a prohibitively large computational task for challenging problems and a more efficient means to determine the convergence of apparent to effective properties is desired. There are numerous sampling-based studies on computational homogenization and only a few have been highlighted above. For a more complete overview, the reader is referred to the review paper [19] for explicit treatment of material randomness, [30] for review of deterministic treatment of classical, second order, and discontinuous homogenization methods, and [6] for general overview of computational homogenization and multiscale methods. In this paper, the requirement in homogenization theory that micromechanical response quantities are statistically homogeneous random fields is exploited to introduce an efficient and systematic approach to determine convergence of apparent properties as a function of SVE size (e.g., see [19] for a discussion on the requirements of statistical homogeneity and ergodicity in computational homogenization). Micromechanical simulations are performed with the target of obtaining convergence of empirically determined autocorrelation functions of response quantities. Based on these results, the variance of the spatial average of a response quantity, (i.e., apparent property) is easily computed as a function of SVE size. The need to perform repeated, costly experiments of MC simulation of SVEs of increasing size is circumvented. Instead, only one judiciously selected SVE size is required to estimate the random field properties. This SVE size is significantly smaller than the RVE size. Therefore, the main contribution of this paper is the introduction of a more efficient approach to predict the convergence behavior of apparent properties than the established brute-force approach. This novel approach utilizes a simple analytical and exact expression relating the autocorrelation function of a random field with the variance of its spatial average. Additionally, this analytical form reveals that the known convergence rate, which states that the variance of apparent properties is proportional to L −d (e.g., [31]), is only

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

255

accurate in the limit as L is large. Here, L and d are the edge length and dimension of the domain, respectively. This asymptotic form omits additional higher order terms, which can lead to incorrect prediction of the convergence rate of apparent properties for SVE sizes typically analyzed in practice. The outline to this paper is as follows. A brief review of stochastic homogenization is given in Section 2. Section 3 derives a simple expression for the variance of spatially averaged quantities of statistically homogeneous random fields and introduces the proposed approach to determine the convergence of apparent properties as a function of SVE size. An analytical example is given for a hyperelastic bar in Section 4. Numerical examples in Sections 5 and 6 demonstrate the extent of applicability to more challenging problems, namely a plane strain linear elastic model and a three dimensional crystal plasticity finite element model of a polycrystalline material. The concluding statements summarize the findings in the paper. 2. Review of homogenization theory The main results of stochastic homogenization for hyperelastic materials are reviewed below. For more details, see [32,33] for periodic, deterministic microstructures and [19,34,35] for stochastic microstructures. Under quasistatic and monotonic loading, the hyperelastic assumption extends to a broad class of constitutive laws by equating the macroscopic strain energy density to the spatial average of the microscopic strain energy density. In the following, bold text is used to denote vector and tensor quantities without distinction, while scalar quantities are not in bold. The equilibrium condition is a stochastic PDE given by ( ) ∇ · σ δ (x, ω) + fδ (x) = 0 in D (1) uδ (x, ω) = u0 (x) on ∂Du σ δ (x, ω) · n(x) = t0 (x) on ∂Dt , where x ∈ D is bounded in R3 , the displacement and traction boundary conditions are prescribed on ∂Du and ∂Dt , respectively, and the boundary ∂D = ∂Du ∪ Dt . The parameter ω represents a sample of an underlying random variable within the probability space (Ω , F, P), which is comprised of sample space Ω , the collection of events F being a σ -algebra on Ω , P ∈ [0, 1] denotes probability, and the space is endowed with a probability measure d P(ω). The parameter δ > 0 is a length-scale parameter such that the limit as δ → 0 corresponds to the ratio of the size of the domain divided by the length scale of material property fluctuation going to ∞, with the notation uδ (x, ·) = u(δ −1 x, ·). The displacement uδ (x, ω) is a vector valued random field and the stress σ δ (x, ω) is a symmetric second order tensor valued random field. The deterministic body force vector fδ (x) is restricted to vary slowly for large δ. The strain energy density W δ (x, ω) defines the constitutive behavior as ∂ W δ (x, ω) = σ δ (x, ω) ∂ϵ δ

(2a)

∂ 2 W δ (x, ω) δ ( )2 = Ct (x, ω), ∂ ϵδ

(2b)

where the strain tensor ϵ δ (x, ω) is a symmetric second order tensor valued random field and the tangent modulus Cδt (x, ω) is a 4th order tensor valued random field having major and minor symmetries. Cδt (x, ω) is necessarily a statistically homogeneous (Eqs. (3a) and (3b)) and ergodic (Eq. (3c)) random field, that is ∫ [ δ ] E Ct (x, ω) = Cδt (·, ω)d P(ω) = constant (3a) Ω

] E Cδt (x1 , ω)Cδt (x2 , ω) = RCδt (x2 − x1 ), [

(3b)

∫ ∫ 1 Cδt (x, ·)dx = Cδt (·, ω)d P(ω) (3c) δ→0 |D| D Ω holds component-wise, where E[·] denotes the expectation operator, R denotes the autocorrelation function, and |·| denotes volume. In other words, the ergodicity property states that the spatial average of any sample of a random field converges to its expected value as the ratio of the size of the sample with respect to its correlation length increases to infinity (i.e., the length-scale of fluctuation of Cδt (x, ω) is ≪ |D|1/3 ). As can be seen in Eqs. (3a) and (3b), the paper utilizes the common presumption that random field quantities can be approximated as weakly statistically lim

256

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

homogeneous second order random fields. This type of random field has the property that its mean is constant, its variance is bounded, and its autocorrelation function only depends on the lag between points. For further details on random field theory, the reader is referred to [36,37]. Following from above, there exists a homogenized strain energy density Whom that does not explicitly depend on x, is deterministic, and defined as ∫ ∫ 1 Whom = lim W δ (x, ·)dx = W δ (·, ω)d P(ω). (4) δ→0 |D| D Ω The second equality in Eq. (4) follows from the ergodicity property. Whom is the strain energy density for the deterministic analog of Eq. (1). The Hill–Mandel condition [38–40] provides the framework to be able to solve for a work conjugate pair ϵ δ (x, ω), σ δ (x, ω) to estimate Whom . There is a large body of work dating back to [38] on this topic, and the reader is referred to the review paper [19] and references therein for further reading. Briefly, a separation of scales is proposed as δ −1 x = y ∈ B as δ → 0 for B being a microstructure of D. To clarify, B is a subdomain of D whose resolution is on the order of the length scale of material heterogeneity such that the spatially varying constitutive properties are smoothly captured. For pragmatic reasons, the size of B is typically < 100× the correlation length of the constitutive properties, which is defined below in Eq. (17). The subsequent analysis in the remainder of the paper considers the restriction of random field quantities to subdomain B, while it is understood that random field quantities are defined over D [41, Chpt. 1.8]. This is a subtle, yet necessary, detail such that the properties of statistical homogeneity and ergodicity remain valid. The following notation is used for spatial and ensemble averages ∫ ∫ 1 Θ(y, ω)d P(ω). (5) Θ(ω) = Θ(y, ω)dy , ⟨Θ(y)⟩ = |B| B Ω The stress and strain tensors are parsed into their macroscale spatial average and microscale spatial fluctuations as σ (y, ω) = σ (ω) + σ ′ (y, ω), ϵ(y, ω) = ϵ(ω) + ϵ ′ (y, ω),

σ ′ (y, ω) = 0

(6a)

ϵ ′ (y, ω) = 0.

(6b)

The microscale spatial fluctuations σ (y, ω) and ϵ (y, ω) are defined as random perturbations that do not affect the average strain energy (i.e., their spatial average over B is zero). Equating the strain energy between the macroscale and microscale responses gives σ (y, ω) : ϵ(y, ω) = σ (ω) : ϵ(ω), which implies that ∫ ( ) ( ) t(y) − σ (ω) · n(y) · u(y) − ϵ(ω) · y dy = 0, σ ′ (y, ω) : ϵ ′ (y, ω) = (7) ′

′

∂B

where the first equality follows from the Divergence theorem. A BVP on B is formulated to solve the microscale stress σ (y, ω) and strain ϵ(y, ω) with boundary conditions chosen to satisfy Eq. (7). These boundary conditions are the stress uniform boundary conditions (SUBC), kinematically uniform boundary conditions (KUBC), or a combination of the two denoted as orthogonal-mixed. The SUBC prescribes uniform traction t(y) = σ 0 · n(y) on ∂B, while the KUBC prescribes uniform displacement u(y) = ϵ 0 · y on ∂B. Solving the KUBC-BVP on B to get σ (ω) for a number of judiciously selected∫values of ϵ 0 leads to a linear system of equations begetting the so-called apparent tangent modulus, ϵ Cδt (ω), as σ (ω) = 0 0 Cδt (ω) : dϵ. Likewise, solving the SUBC-BVP on B to get ϵ(ω) for a number of judiciously selected values the so-called apparent tangent compliance, Sδt (ω), ∫ σ 0 ofδ σ 0 leads to a linear system of equations begetting δ as ϵ(ω) = 0 St (ω) : dσ . The homogenized properties Ct (ω) and Sδt (ω) are denoted as apparent properties as they ef f ef f are random variables and converge to deterministic constants, denoted as effective properties (i.e., Ct and St , respectively), by virtue of ergodicity, that is ∫ ϵ0 ef f ef f lim σ (ω) = ⟨σ ⟩ = Ct : dϵ ⇒ lim Cδt (ω) = Ct (8a) δ→0

lim ϵ(ω) = ⟨ϵ⟩ =

δ→0

δ→0

0 σ0

∫

ef f

St

: dσ ⇒ lim Sδt (ω) = St . ef f

δ→0

0

(8b)

The rate of convergence of Cδt (ω) → Ct (Sδt (ω) → St ) as a function of the size of |B|1/3 is discussed in Section 3. It is noted that the applied boundary conditions induce statistical inhomogeneity in the response quantities, violating the requirement of statistical homogeneity. However, this effect of the boundary conditions is mostly ef f

ef f

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

257

Fig. 1. Depiction of the state-of-the-art strategy to determine the convergence rate of apparent properties as a function of SVE size. Brute-force MC simulation of numerous SVEs of increasing size is performed. The computed variance from the MC results are the black diamonds. A curve is fit to these points to establish the convergence rate, while the RVE size is the size where the variance of apparent properties is less than some pre-determined threshold.

confined to the region of the domain near the boundaries, while the response in the interior of the domain for reasonably sized SVEs is statistically homogeneous. This conclusion can easily be drawn from the numerous studies that show that regardless of the boundary conditions chosen (i.e., either KUBC, SUBC, or orthogonal-mixed), (1) the apparent properties converge to the same effective properties and (2) the response over this domain is approximately ergodic. This can only occur if the influence of the boundary conditions diminishes and the vast majority of the response is statistically homogeneous. 3. Convergence behavior of homogenized properties The current established approach to determine the convergence of apparent properties as a function of SVE size is given pictorially in Fig. 1. The aim is to determine the size where the variance of apparent properties is considered negligible, denoted as the RVE. A series of SVE sizes are selected and the variance of apparent properties for each size is computed. Thus, for each SVE size, a MC simulation is performed by generating numerous statistically equivalent microstructures and solving a set of BVPs as described above in Section 2 for each statistically equivalent microstructure. These results correspond to the black diamonds in Fig. 1. A curve is then fit to the computed variances to extrapolate to the size of the RVE. This brute-force approach can be computationally prohibitive for many problems of interest. There are a few studies and proofs to determine the rate at which [the ]variance of an apparent property converges to zero. Ref. [31] proposes a theorem bounding the variance as Var Cδt ≤ θ/L d : d > 2 component-wise for L being the size of one dimension in B, d denoting the dimension of the domain, and θ being a positive scalar. This is for a discrete linear elliptic operator on a scalar valued solution field whose constitutive matrix is constructed from independent identically distributed coefficients, while in [35] the theorem is extended for continuum linear elliptic operators. This is considered to be an improvement over the algebraic convergence rates derived in [42]. For fully nonlinear elliptic problems a generic logarithmic convergence rate is derived in[ [43]. ] One of the most well known formulas is from [12] which states that Var Cδt = A/|B| for any component of the elastic tensor and parameter A is the so-called integral range of the random field that is[ spatially averaged to ] obtain Cδt (ω). Rather than directly computing A, parameters θ and α are fit empirically to Var Cδt ≈ θ/|B|α through MC simulation of increasing size SVEs as described in Fig. 1. Based on extensive examples, that work recommends performing MC simulation on 4–5 different size SVEs to determine θ and α. Therefore, the computational burden associated with their approach can be significant.

258

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

In this work, it is proposed to directly compute the variance of an apparent property from the autocorrelation function of its corresponding response quantity of interest, which is calculated from MC simulation at a single SVE size. It is emphasized that this approach expresses the variance in terms of the autocorrelation function of a response quantity and not that of the input uncertainty. The convergence of spatially averaged quantities of statistically homogeneous random fields is given exactly as follows. Let ζ (y, ω) denote the restriction to B of a zero-mean statistically ∫homogeneous random field (e.g., any component of ϵ(y, ω), σ (y, ω), or W (y, ω)) with spatial average 1 ζ (ω) = |B| B ζ (y, ω)dy, then the following results hold (derived for d = 3, while it is straightforward to modify for d = 1 or d = 2) ∫ ∫ [ ] 1 Var ζ = Rζ (y1 , y2 )dy1 dy2 |B|2 B B (9) ∫ L∫ L∫ L 8 = 6 (L − ξ1 )(L − ξ2 )(L − ξ3 )Rζ (ξ )dξ1 dξ2 dξ3 , L 0 0 0 where Rζ is the autocorrelation function of random field ζ (y, ω). For simplicity, but without loss of generality, it is assumed that B is a cube with side length L. The first equality follows from the definition of the autocorrelation function, while the second equality utilizes the statistical homogeneity property with lag ξ = y2 − y1 . The second equality is obtained by substituting ξ for one of y1 or y2 , rearranging the order of integration, and modifying the integration limits [44, Chpt. 3]. Assuming the common property of autocorrelation functions that limh→∞ Rζ (hn) ≤ c1 h −1 for arbitrary unit normal n and positive, finite constant c1 , then the variance can be expressed as [ ] θ1 (L) θ2 (L) θ3 (L) θ4 (L) Var ζ = + + + . (10) L3 L4 L5 L6 The θi : i = 1 : 4 are finite scalar functions of L and are simply the integral expressions obtained from expanding the integrand in Eq. (9), given as ∫ L∫ L∫ L θ1 (L) = 8 Rζ (ξ )dξ1 dξ2 dξ3 (11a) 0

θ2 (L) = −8

0

∫ 0

θ3 (L) = 8

L

∫

L

L

∫

0

L

(ξ1 + ξ2 + ξ3 )Rζ (ξ )dξ1 dξ2 dξ3

(11b)

0 L

∫

0

∫

L

∫

0

0

θ4 (L) = −8

0

L∫

(ξ1 ξ2 + ξ1 ξ3 + ξ2 ξ3 )Rζ (ξ )dξ1 dξ2 dξ3

(11c)

0

∫ 0

L

L

∫

ξ1 ξ2 ξ3 Rζ (ξ )dξ1 dξ2 dξ3 .

(11d)

0

Eq. (10) reveals that the variance of the spatial average of a statistically homogeneous random field converges to 0 at a rate that has contributions from terms proportional to L14 , L15 , L16 , in addition to L13 . As L increases beyond the correlation length, the[ θi]’s become constant, and the 1/L 3 term dominates the other terms. Therefore, the known convergence rate Var ζ ∝ 1/L d is an asymptotic solution while Eqs. (9) and (10) are exact for all L. By knowing the autocorrelation of the response quantity of [interest, Eq. (11) provides a simple way to assess if the SVE size ] analyzed is within the asymptotic range for Var ζ = c1 /L d to be a valid approximation, with c1 being a positive, finite constant. An analytically derived solution for a bar in Section 4 supplements these findings. In general, numerical simulation is required to determine apparent properties, and the accuracy of the approach is demonstrated for finite element analysis for a two-dimensional problem in Section 5 and a three-dimensional problem in Section 6. MC simulation is conducted with the aim of determining the autocorrelation function of the response quantities from which apparent properties are derived in order to evaluate Eq. (9). The strategy is as follows: • Judiciously determine an SVE size that is large enough to encompass enough fluctuations of the response quantities such that spatial coordinates far apart within the SVE are uncorrelated (i.e., the domain over which the empirically computed autocorrelation function of a response quantity should be large enough to observe the autocorrelation function’s decay to zero).

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

259

Fig. 2. Bar subjected to isostress condition.

• Perform MC simulation by generating random field samples of the input uncertainty and repeatedly solving the BVP. • Estimate the autocorrelation function of appropriate response quantities empirically by computing the covariance matrix of the response quantity samples and averaging all coordinate pairs having the same lag ξ . • If necessary, fit a functional form to the empirical autocorrelation function in order to extrapolate beyond the analysis domain.

4. Analytical example: homogenization of bars Analytical expressions for the apparent tangent compliance of hyperelastic bars corresponding to the isostress scenario are derived below in order to demonstrate the role of correlation length on the convergence of spatially averaged quantities. Analytical expressions also exist for the apparent tangent modulus derived from the isostrain condition but are omitted for brevity. Since these are analytical solutions, the variance of the spatially averaged quantity can be directly related to the autocorrelation of the input uncertainty (i.e., elastic modulus or compliance). Consider the bar in Fig. 2 with length L = δ −1 L 0 and external load σ0 . The length L 0 is an arbitrary reference length and δ varies to provide varying beam lengths with respect to L 0 . The beam is modeled by the small strain, monomial constitutive law σ = Eϵ α for α > 0 with its compliance modeled as a random field defined over spatial coordinate y given as 1 = D0 (1 + f (y)). E(y)

(12)

Here, f (y) is a zero-mean statistically homogeneous random field with variance σ 2f , D0 is the mean value of the compliance, and E(y) is the elastic modulus. For the remainder of the paper the explicit dependence of random quantities on the underlying random sample ω is dropped for notational convenience. The apparent tangent compliance Dtδ is obtained by Dtδ dσ δ = dϵ δ . Since the bar is statically determinate, ∫ ∫ 1/α L 1 1/α δ δ δ δ = 1 σ 1/α−1 dσ D 1/α L (1 + f (y))1/α dy. σ = σ0 ⇒ dσ = dσ0 , and ϵ = L (σ0 D0 ) (1 + f (y)) dy ⇒ dϵ 0 0 0 0 αL 0 Solving for Dtδ gives ∫ L 1 1/α−1 1/α δ Dt = D0 (1 + f (y))1/α dy. (13) (σ0 ) αL 0 The variance of Dtδ is a function of the parameter α. For demonstration purposes the values of α = 1/2 and α = 1 are selected. For α = 1, corresponding to the linear elastic case, the variance is ( )2 ∫ L [ ] D0 Var Dtδ = 2 (L − ξ )R f (ξ )dξ, (14) L 0 for [R f (ξ]) being the autocorrelation function of f (y). Following the form in Eq. (10) this can be re-written as Var Dtδ = θ1L(L) + θ2L(L) 2 , where ∫ L θ1 (L) = 2D02 R f (ξ )dξ 0 (15) ∫ L

θ2 (L) = −2D02

ξ R f (ξ )dξ.

0

260

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

Likewise, for α = 1/2, the results are [∫ L ] [ ] σ 2 D4 L2 4 (L − ξ )(4R11 (ξ ) + 4R12 (ξ ) + R22 (ξ ))dξ − Var Dtδ = 8 0 2 0 σf L 2 0 θ1 (L) =

8σ02 D04

L

∫

4R11 (ξ ) + 4R12 (ξ ) + R22 (ξ ) − σ 4f dξ

0

θ2 (L) =

−8σ02 D04

(16a)

(16b)

L

∫

ξ 4R11 (ξ ) + 4R12 (ξ ) + R22 (ξ ) − (

σ 4f

)

dξ,

0

∫L∫L ∫L [ ] where Ri j (ξ ) = E f (y1 )i f (y2 ) j and 0 0 Ri j (ξ )dy1 dy2 = 2 0 (L − ξ )Ri j (ξ )dξ . The convergence behavior of the apparent tangent compliance as a function of bar length is shown in Figs. 3(a) and 3(b) for α = 1 and α = 1/2, respectively. Random field f (y) in Eq. (12) is modeled by a zero-mean, scaled and shifted lognormal distribution with autocorrelation function R11 (ξ ) = exp(Rg (ξ ))−1. This lognormal random field is given as f (y) = exp(σg g(y)−σg2 /2)−1 where g(y) is an underlying zero-mean Gaussian random field with variance σg2 = .01 √ [ ] [ ] [ ] δ 2 2 and autocorrelation function Rg = σg exp(−ξ /100). The coefficient of variation Cov Dt = Var Dtδ /E Dtδ is plotted, rather than the variance, since it is independent of load level σ0 . It is plotted as a function of bar length normalized by correlation length η. The correlation length for random field f (y) : y ∈ Rd with autocorrelation function R f (ξ ) is defined in this work to be )1/d ( ∫ ∫ ∫ ∞ ∞ ∞ 1 R f (ξ )dξ1 dξ2 ...dξd . (17) ... η= σ 2f 0 0 0 The correlation length for the random field in this example equals 9.04. Note that units are not provided in this example since the results presented are normalized to be unitless. Figs. 3(a) and 3(b) plot the exact coefficient of variation of the apparent compliance in solid lines and the coefficient of variation approximated by only considering the term involving θ1 (i.e., omitting the contribution of the higher order term involving θ2 ) in dashed lines. The following observations can be made: (1) the decay rate is slow for d = 1, (2) the decay rate depends on both the constitutive law parameter α and the correlation length of the input random field, and (3) for bar lengths less than 15× the correlation length of the compliance, there are noticeable errors in the approximated coefficient of variation when only the contribution from the term involving θ1 is considered. The α = 1/2 case has larger coefficient of variation for a given size than for α = 1, and this is common for nonlinear materials (i.e., the RVE size increases due to material nonlinearity). For this special case of a monotonically loaded statically determinate bar, the convergence behavior of the coefficient of variation is independent of load level σ0 , but in general this is not the case. It is worth noting that the variance converges proportional to the square of the coefficient of variation and thus will converge faster; but, nonetheless, Fig. 3 illustrates that the SVE must contain many fluctuations of the random field to have a small variance in the apparent properties. Although this is well known, extremely large SVEs are rarely analyzed in practice. Thus, it is important to understand the impact of the higher order term on the convergence rate of the variance. 5. Plane strain elasticity The analytical derivations in Section 4 are only attainable for statically determinate structures. The examples in this section and in Section 6 demonstrate how to establish the convergence of apparent properties for general problems. Fig. 4 shows an L × L square domain plane strain model subject to a uniform traction σ22 = 51 kPa at the n = (0, 1) face with displacement boundary conditions that restrict motion in the y2 direction at the n = (0, −1) face and traction-free boundary conditions at the n = (−1, 0) and n = (1, 0) faces. The contours in the figure correspond to one sample of the stochastic elastic modulus, where the compliance is modeled by a scaled and shifted, zero-mean lognormally distributed random field f (y1 , y2 ), such that the elastic modulus is E(y1 , y2 ) = E 0 /(1 + f (y1 , y2 )). The lognormal distribution is defined as a mapping from a standard normal Gaussian distribution g(y1 , y2 ) as f (y1 , y2 ) = exp(m + g(y1 , y2 )s) + a,

(18)

where s = .37, m = −1.23, and a = −.313. The autocorrelation function of the underlying √Gaussian field is given as Rg (ξ1 , ξ2 ) = exp(−c(ξ12 + ξ22 )) for c =2/5 mm−2 . The correlation length for this field is π/(4c). The remaining

261

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

(a) Linear case α = 1.

(b) Nonlinear case α = 1/2.

[ ] Fig. 3. Coefficient of variation Cov Dtδ as a function of bar length L normalized by correlation length η, defined in Eq. (17), for (a) α = 1 and (b) α = 1/2. The solid line is the exact coefficient of variation while the dashed line is the approximate coefficient of variation by only considering term θ1 and omitting the contribution from higher order term θ2 . Noticeable errors of the asymptotic form exist up to L/η ≈ 15.

Fig. 4. Plane strain domain analyzed. The contour is one sample of the elastic modulus providing a qualitative depiction of the spatial heterogeneity.

model parameters are Poisson’s ratio ν = .3, and E 0 = 1 GPa. Three MC simulations are conducted for domain sizes L = 14.25 mm, L = 17.80 mm, and L = 21.42 mm using 750 samples each. The FE meshes are discretized into 74 × 74, 93 × 93 and 112 × 112 square 4 node linear elements for the small, medium, and large domains, respectively. The ensuing element size of 0.19 mm × 0.19 mm, which is the same for all 3 domain sizes studied, is sufficiently refined to smoothly represent the input random field and mechanical response quantities. Samples of the underlying Gaussian field g(y1 , y2 ) are generated via the spectral representation method [45], and then samples of the lognormal random field f (y1 , y2 ) are obtained by Eq. (18). The finite element analysis is done using FEAP [46]. Following the MC simulations of the smaller domain, samples of the strain energy density W are calculated for ∑ Nint (σ j : ϵ j ) for the number of each finite element as the average over the element integration points as W = 2N1int j=1 integration points per element Nint = 4. The autocorrelation function of W is obtained empirically by computing the covariance matrix (i.e., covariance of every pairwise product in the FE domain) for each point over the total FE domain and averaging among pairwise products having the same lag. In order to extrapolate beyond the analysis domain, this is then approximated using a function of the form RW (ξ1 , ξ2 ) = σW2 exp(−|a1 ξ1 + a2 ξ2 + a3 ξ1 ξ2 + a4 ξ12 + a5 ξ22 |), σW2

(19)

where is the computed variance of W and parameters ai : i = 1, 5 are optimized by minimizing the ℓ -norm of RW and the empirical autocorrelation function. Fig. 5(a) plots the empirical autocorrelation function in solid lines and the fitted autocorrelation function in dashed lines along the y1 direction for y2 = 0 (blue lines), the y2 direction 2

262

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

(a) Empirical (solid) and fitted (dashed) autocorrelation function using Eq. (19).

(b) Coefficient of variation of W versus domain size.

Fig. 5. Fig. 5(a): the empirical autocorrelation function of W (solid lines) and the fitted autocorrelation using Eq. (19) (dashed lines) along the y1 direction for y2 = 0 (blue lines), the y2 direction for y1 = 0 (red lines), and the direction y1 = y2 (black lines). The empirical autocorrelation function is from post processing the MC simulation of the small domain (i.e., the L = 14.25 mm domain). Fig. 5(b): the coefficient of variation of W as a function of edge length L normalized by correlation length η of the compliance (i.e., input uncertainty). The solid line is using the formula of Eq. (9) while the dashed line is from considering the θ1 term’s contribution only. The computed variance for L = 17.80 mm and L = 21.42 mm is very accurately predicted when using Eq. (9), while noticeable errors exist for the asymptotic form up to L/η ≈ 20. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

for y1 = 0 (red lines), and the direction y1 = y2 (black lines). Following these results, the variance of the spatially averaged strain energy density W is computed as a function of domain size substituting the autocorrelation function in Eq. (19) in Eq. (9). Additionally, the convergence rate is when omitting the contribution of the higher ∫ Lpredicted ∫L [ ] order terms, that is Var W ≈ θ1 (L)/L 2 with θ1 (L) = 4 0 0 RW (ξ1 , ξ2 )dξ1 dξ2 . Fig. 5(b) plots the coefficient of variation of these two curves with the solid line for the exact solution and the dashed line for the single-term solution. The mean value taken for the denominator for the coefficient of variation is the mean of W computed from the MC simulation of the largest domain. The coefficient of variation is plotted rather than variance because it is easier to interpret as a measure of spread since it is non-dimensionalized. Noticeable discrepancies exist between the exact and single-term prediction for edge lengths up to 20× the correlation length of the input uncertainty. Observe that the convergence rate for the two-dimensional example is slightly faster than the one-dimensional, linear example. For example at L/η = 10, the coefficient of variation of W is 2.5% for the plane strain example while for Dtδ it is 4.5% in the linear, bar example. Note that the coefficient of variation for Dtδ is the same as the strain energy density for the bar example because the stress is a constant that factors out of the calculation. 6. Application to crystal plasticity micromechanics In this example, the convergence behavior of apparent properties as a function of SVE size is evaluated for a morphology-resolved, crystal plasticity model of metallic microstructures. The IN100 Nickel alloy three-dimensional data set available through the DREAM.3D project [47] is used to fit the parameters of a tessellation through which statistically equivalent microstructures are generated using stochastic simulation. A classical crystal plasticity constitutive law (e.g., [48]) is utilized in a large strain finite element setting to analyze each simulated microstructure sample. Post processing the results reveals the convergence behavior of apparent properties. 6.1. Stochastic simulation of microstructure morphology In a separate paper [49], the authors developed a type of ellipsoidal growth tessellation (EGT) model to represent polycrystalline morphology and crystallography. That paper introduced an efficient method to obtain the probability distributions of the EGT model parameters given an Electron Backscatter Diffraction (EBSD) data set as well as a method to simulate the parameters consistent with their distributions. The EGT consists of ellipsoids growing in

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

263

Fig. 6. A sample realization of the EGT model whose parameters are fit to an IN100 Nickel data set. Each color represents a grain, which is a contiguous region of material having the same crystallographic orientation.

space from points, denoted as nucleation sites, until they abut each other. It is modeled as a vector-valued marked point process [50, Chpt. 4], which is a point process with a vector of parameters associated with each point. The nucleation sites of the grains are modeled as a point process. The associated mark vectors represent the velocity, spatial orientation, and crystallographic orientation of the growing ellipsoids defining the polycrystalline structure. A sample of a simulated microstructure using this model is shown in Fig. 6. Details of the approach are provided in [49]. This approach is used in the current work to generate sample SVE microstructures. 6.2. Constitutive formulation A widely used crystal plasticity constitutive and kinematic model is implemented for this study and briefly described as follows. The deformation gradient F ≡ ∇Y y is multiplicatively decomposed into its elastic and plastic parts, F∗ and F p , respectively, as F = F∗ F p ,

(20)

with y and Y denoting the spatial and reference configuration frames, respectively. This gives the velocity gradient to be ∇y v ≡ L = L∗ + L p , with its plastic component defined as Nsli p

L p = F˙ p F− p =

∑ α=1

Nsli p

γ˙ α s∗α ⊗ m∗α =

∑

γ˙ α S∗α 0 ,

(21)

α=1

where γ α , sα and mα denote the shear strain, slip direction, and slip plane normal, respectively, for slip system α, the asterisk superscript denotes the crystallographic lattice configuration frame, and (˙) denotes the time rate of change. ∗α The Schmid tensor S∗α ⊗ m∗α , defined as the tensor product of the slip system direction with the slip system 0 = s normal, maps the Cauchy stress σ to the slip system resolved shear stress τ α . The constitutive relation of the slip system shear strain rate and resolved shear stress follows the classical power law formula, given as ⏐ α ⏐n ⏐τ ⏐ (22a) γ˙ α = γ˙0 ⏐⏐ α ⏐⏐ sign(τ α ) g Nsli p α

g˙ =

∑

h αβ |γ˙ β |

(22b)

β=1

h αβ

⏐r ⏐ ( ) ⏐ β ⏐ gβ g ⏐ ⏐ αβ β αβ β = q h = q h 0 ⏐1 − β ⏐ sign 1 − β , ⏐ gsat ⏐ gsat

(22c)

where g α is the critically resolved shear stress on α, h αβ is the hardening law of g α . Model parameter q αβ provides α tuning of latent and self hardening of slip system α while gsat represents the saturation stress, which is the upper

264

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270 Table 1 Slip systems for FCC materials in local coordinates of the crystal. The ¯ denotes a negative sign. overbar (·) Slip system index α

Slip plane normal mα

Slip direction sα

1

√1 (1 3 √1 (1 3 √1 (1 3 ¯ √1 (1 3 ¯ √1 (1 3 1 ¯ √ (1 3 1 √ (1 3 √1 (1 3 √1 (1 3 √1 (1 3 √1 (1 3 √1 (1 3

¯ √1 [1 2 1 ¯ √ [1 2 √1 [0 2 √1 [1 2 √1 [1 2 √1 [0 2 ¯ √1 [1 2 √1 [1 2 √1 [0 2 √1 [1 2 √1 [0 2 ¯ √1 [1 2

2 3 4 5 6 7 8 9 10 11 12

Table 2 Model constitutive parameters. Quantities with units are in MPa. ⏐ (α) (α) (α) γ0 q (α,α) q (α,β) ⏐α̸=β h0 gsat 1500

.001

1.0

1.4

550

1 1) 1 1) 1 1) 1 1) 1 1) 1 1) ¯ 1 1) ¯ 1 1) ¯ 1 1) 1¯ 1) 1¯ 1) 1¯ 1)

n 20.0

1 0] 0 1] 1¯ 1] 1 0] 0 1] 1¯ 1] 1 0] 0 1] 1 1] 1 0] 1 1] 0 1]

r

C11

C12

C44

1.3

201.7e3

134.4e3

104.5e3

β

limit of shear stress on α. Model parameters r , n, h 0 , and γ˙0 control the shape of the stress–strain curve. The time integration method developed in [51] is directly followed and implemented in a parallel version of FEAP [46] as a UMAT. The IN100 Nickel alloy has a face centered cubic (FCC) crystal structure, and its slip systems are given in Table 1. The elastic constants C11 = C22 = C33 ; C12 = C23 = C13 ; C44 = C55 = C66 with all other entries in C being zero. The UMAT is implemented with the commonly used constitutive model parameter values for Nickel super alloys and are given in Table 2. 6.3. Numerical results Crystal plasticity finite element simulation results are presented for 231 samples of statistically equivalent IN100 Nickel alloy microstructures. The domain is a cube with edge length 12 µm and the average number of grains in this domain is 52. The boundary conditions are applied tensile forces at all nodes on the n = (0, 0, 1) face corresponding to a uniform normal traction, and normal direction displacements are enforced to be zero on all nodes on the n = (−1, 0, 0), n = (0, −1, 0), and n = (0, 0, −1) faces. The domain is discretized using a structured mesh that has 32 × 32 × 32 cube elements with edge length .375 µm. The finite elements are standard 8 node brick elements with tri-linear shape functions and 2 × 2 × 2 integration points. This gives on average 64 elements for a (1.5 × 1.5 × 1.5) µm3 grain, which is approximately the smallest grain size modeled. The stress controlled loading is incremented with the target of inducing an average strain rate of 1.6e−4 /s for 300 s, leading to 5% average strain. The conjugate gradient iterative solver with a block Jacobi pre-conditioner is used to solve the system of equations for each time increment using PETSC [52], and each time increment is allowed to vary according to FEAP’s built in time stepping algorithm. Each simulation takes approximately 8 h on 6 processors. Fig. 7 plots the spatially averaged stress–strain for all 231 simulations. In addition to these simulations, another MC simulation study is performed on a larger cubic domain microstructure of edge length 18 µm (i.e., 1.5× the edge length of the previous study) having an average of 175 grains. The element size is unchanged, resulting in an FE mesh size of 48 × 48 × 48. Each simulation takes approximately 12 h on 8 processors. A total of 150 simulations are performed and the spatially averaged stress–strain curves for all simulations

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

265

Fig. 7. Spatially averaged stress–strain curves for all 231 simulations for the (12 µm)3 sized SVE. The dashed horizontal lines are selected load levels where the variance of the homogenized strain energy density is analyzed. ‘L.L.’ denotes load level snapshot. The images are the total 33 strain for an example simulation at L.L. 1 and L.L.6.

Fig. 8. Spatially averaged stress–strain curves for all 150 simulations for the (18 µm)3 sized SVE. The dashed horizontal lines are selected load levels where the variance of the homogenized strain energy density is analyzed. ‘L.L.’ denotes load level snapshot. The images are the total 33 strain for an example simulation at L.L. 1 and L.L.6.

are plotted in Fig. 8. These simulations serve as reference data to evaluate the extrapolated results obtained from the 32 × 32 × 32 mesh. The variance of the spatially averaged strain energy density W is studied as a function of size based on its autocorrelation function. The empirical autocorrelation is computed by approximating W as a statistically homogeneous random field. First, a single W is computed∑ for each element by averaging over the element’s integration Nint j j j j points and computing for each load step i: Wi = 2N1int j=1 (σ i + σ i−1 ) : (ϵ i − ϵ i−1 ) + Wi−1 , for Nint being the number of integration points per element, which is 8 for this example. The covariance matrix of Wi (i.e., covariance of every pairwise product in the FE domain) is computed, and then an average is taken among all products with the √ 2 same separation. Fig. 9 plots the empirical autocorrelation for (ξ1 , 0, 0) , (0, ξ2 , 0) , (0, 0, ξ3 ), and r = ξ1 + ξ22 + ξ32 for load level snapshots 1 (Fig. 9(a)), 3 (Fig. 9(b)), 4 (Fig. 9(c)) and 6 (Fig. 9(d)). The empirical autocorrelation function suffers from numerical error at the larger lags due to the fewer number of pairwise data points at larger lags. Nonetheless, both the SVE and number of MC simulations are large enough to ascertain the pattern of the

266

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

(a) Load level 1.

(c) Load level 4.

(b) Load level 3.

(d) Load level 6.

Fig. 9. Empirical autocorrelation function for W (dashed lines) along lines (ξ1 , 0, 0), (0, ξ2 , 0), (0, 0, ξ3 ), and ξ1 = ξ2 = ξ3 , along with its isotropic approximation RW in Eq. (23) (solid lines) for load level snapshots 1, 3, 4, and 6. RW is used to extrapolate the variance of W as a function of SVE size.

autocorrelation function. Observe that the autocorrelation function is very close to being isotropic. In order to extrapolate beyond the analysis domain, it is approximated by ( ) √ RW (ξ1 , ξ2 , ξ3 ) = σW2 exp −c0 ξ12 + ξ22 + ξ32 , (23) where σW2 is the computed variance of W . This approximation is plotted along with the empirical autocorrelation function in Fig. 9 as the thick black line. The parameter c0 is determined through optimization for each load level snapshot by minimizing the square of the ℓ2 -norm between the empirical autocorrelation and Eq. (23). Due to the numerical error in computing the empirical autocorrelation for larger lags, the minimization is performed on a 16 × 16 × 16 subdomain (i.e., up to one-half the total lag for each direction). By having a functional form for the autocorrelation function, the variance of W as a function of domain size can be computed using Eq. (9). The predicted variance for the spatially averaged strain energy density for an SVE size of (18 µm)3 is compared to the variance computed from the MC simulations of the larger SVE. Fig. 10 plots the computed variance of the larger SVE in red and the predicted variance using the autocorrelation function in Eq. (23) in black versus the load level snapshots. As a comparison, Fig. 10 also shows these results for the (12 µm)3 SVE. The error bars are the 95% confidence intervals of the prediction due to sampling error (i.e., not enough MC simulations performed) computed through bootstrapping [53]. The predictions from the autocorrelation function are very close to the estimates computed from MC simulation in comparison to the amount of spread in the estimator. These figures show increased spread in the estimated variance of the spatially averaged strain energy density function as a function of load level (this is also reflected in the stress–strain curves in Figs. 7 and 8). Fig. 11 plots the coefficient of variation as a function of edge length L normalized by average grain size ⟨D⟩ = 2 µm for load level snapshots 1, 3, 4, and 6. The solid lines are derived from the variances predicted using Eq. (9) while

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

(a) Smaller SVE.

267

(b) Larger SVE.

Fig. 10. Comparison of the predicted variance of the strain energy density (black curve) of an SVE size of (a) (12 µm)3 and (b) (18 µm)3 with that computed through MC simulation (red curve). The predictions are computed via Eq. (9) using the autocorrelation function in Eq. (23), which is fitted to the empirical autocorrelation function from MC results of an SVE size of (12 µm)3 . The error bars are the 95% confidence interval of the variance estimator using the statistical bootstrapping technique. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the [dashed ∫ L from ∫ L ∫ Lthe variances predicted by omitting contributions from the higher order terms, that is, ] lines are Var W ≈ 8/L 3 0 0 0 RW (ξ1 , ξ2 , ξ3 )dξ1 dξ2 dξ3 . The denominator for the coefficient of variation for each load level is taken to be the mean of W computed from the MC simulation of the larger domain. Fig. 11 shows that the asymptotic form for the variance of apparent properties is valid only after a very large SVE size is reached, typically much larger than what is analyzed in practice. The error of the asymptotic form increases as the plasticity increases, since the correlation length of response quantities tends to increase with increasing plasticity. The slight drop in the coefficient of variation between load level 4 to 6, is predominantly due to the increasing mean value, while the variance of load level 6 is much larger than that of load level 4. These results suggest that using the asymptotic form in practice leads to an incorrect prediction of the RVE size. It is difficult to make comparisons between these results and the oneand two-dimensional examples because ⟨D⟩ is not identical to the correlation length. Nonetheless, the trend and scale of the error of the asymptotic approximation for the variance as function of domain size appear to be very similar for linear materials. To recall, 231 and 150 simulations were performed for the small and large SVE, respectively. As these simulations are computationally intensive, it is currently impractical to perform MC of simulations of much larger SVEs in order to observe the convergence behavior with high confidence, especially for larger load levels. Therefore, this motivates the need for alternatives to brute force techniques, such as that proposed here, to evaluate convergence of properties to deterministic values as a function of SVE size. 7. Conclusions In this work, the evaluation of the RVE size in computational homogenization is approached by treating mechanical response quantities as statistically homogeneous random fields and exploiting their properties. A novel approach to determine the convergence of apparent properties to effective properties as a function of domain size is introduced, and its accuracy is demonstrated through 3 examples covering problems in one-, two-, and three-dimensions. Rather than performing MC simulation on a number of microstructures of increasing size SVEs to determine the convergence of spatially averaged quantities, which can be computationally intractable for many problems, efforts are targeted towards estimating the autocorrelation function of the response quantities through MC simulations of one size SVE. The size of the SVE and number of MC simulations must be determined through scientific judgment such that (1) the empirical autocorrelation function is sufficiently described such that a function form can be fit to it, and (2) the statistical inhomogeneity caused by the boundary conditions has sufficiently diminished. Once the autocorrelation function is established, then the variance of apparent properties as a function of size can be predicted through a trivial computation of the integration of the autocorrelation function over the domain of the SVE.

268

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

(a) Load level 1.

(c) Load level 4.

(b) Load level 3.

(d) Load level 6.

[ ] Fig. 11. Coefficient of variation Cov W for load level snapshots 1, 3, 4, and 6 as a function of edge length L normalized by average grain diameter ⟨D⟩ = 2 µm. The solid lines are the predictions using the formula of Eq. (9), while the dashed lines are the predictions by omitting the contributions of higher order terms θ2 (L), θ3 (L), θ4 (L) in (11). The computed variance for the large domain is very accurately predicted when using Eq. (9) based on the autocorrelation function determined from the smaller domain. Very large errors exist when predicting the variance by using the asymptotic form, and the error increases with plasticity.

Additionally, through a simple derivation, the variance of apparent properties can be expressed as a sum of separate terms of order 1/L 3 , 1/L 4 , 1/L 5 , and 1/L 6 , which the 1/L 3 term dominates as L increases for three-dimensional problems. It is straightforward to derive the corresponding formulas for one- and two-dimensional problems. Once the autocorrelation function of the response quantity of interest is known, it is trivial to determine the contributions of the higher order terms in order to assess if the SVE size analyzed is large enough such that the asymptotic formula for convergence to effective properties is valid. Based on the examples in this paper, the higher order terms become negligible for SVE sizes approximately > 15× the correlation length of the input uncertainty for both the one- and two-dimensional problems. For the three-dimensional problem, the length scale where the higher order terms are negligible depends on the amount on nonlinearity, ranging from > 20× to > 30× the average grain diameter from linear to 4% average strain. The coefficient of variation drops to 1% for the polycrystalline sample at approximately 50× the average grain diameter, which results in an RVE size on the order of [3500–4000] grains. Acknowledgments This work has been partially supported through a grant (No. FA9550-12-1-0445) to the Center of Excellence on Integrated Materials Modeling (CEIMM) at Johns Hopkins University (partners JHU, UIUC, UCSB), awarded by the AFOSR/RSL (Computational Mathematics Program, Manager Dr. A. Sayir) and AFRL/RX (Monitors Dr. C. Woodward and C. Przybyla). Partial funding for this project was provided by the Office of Naval Research (ONR)

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

269

through the Naval Research Laboratory’s Basic Research Program (No. N0001418WX00093). The finite element analysis was supported by the Department of Defense High Performance Computing Modernization Program using the Army Research Laboratory and Air Force Research Laboratory Major Shared Resource Centers under project 416, subproject 231. References [1] J.H. Panchal, S.R. Kalidindi, D.L. McDowell, Key computational modeling issues in integrated computational materials engineering, Comput. Aided Des. 45 (1) (2013) 4–25. [2] S. Hazanov, C. Huet, Order relationships for boundary conditions effect in heterogeneous bodies smaller than the representative volume, J. Mech. Phys. Solids 42 (12) (1994) 1995–2011. [3] C. Huet, Application of variational concepts to size effects in elastic heterogeneous bodies, J. Mech. Phys. Solids 38 (6) (1990) 813–841. [4] S.R. Arwade, G. Deodatis, Variability response functions for effective material properties, Probab. Eng. Mech. 26 (2) (2011) 174–181. [5] F. Fritzen, S. Forest, T. Böhlke, D. Kondo, T. Kanit, Computational homogenization of elasto-plastic porous metals, Int. J. Plast. 29 (2012) 102–119. [6] S. Ghosh, Micromechanical Analysis and Multi-Scale Modeling Using the Voronoi Cell Finite Element Method, CRC Press, 2011. [7] S. Ghosh, K. Lee, S. Moorthy, Multiple scale analysis of heterogeneous elastic structures using homogenization theory and Voronoi cell finite element method, Int. J. Solids Struct. 32 (1) (1995) 27–62. [8] E. Ghossein, M. Lévesque, Homogenization models for predicting local field statistics in ellipsoidal particles reinforced composites: Comparisons and validations, Int. J. Solids Struct. 58 (2015) 91–105. [9] T. Hoang, M. Guerich, J. Yvonnet, Determining the size of RVE for nonlinear random composites in an incremental computational homogenization framework, J. Eng. Mech. 142 (5) (2016) 04016018. [10] C. Huet, Coupled size and boundary-condition effects in viscoelastic heterogeneous and composite bodies, Mech. Mater. 31 (12) (1999) 787–829. [11] M. Jiang, M. Ostoja-Starzewski, I. Jasiuk, Scale-dependent bounds on effective elastoplastic response of random composites, J. Mech. Phys. Solids 49 (3) (2001) 655–673. [12] 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, Int. J. Solids Struct. 40 (13) (2003) 3647–3679. [13] Z. Khisaeva, M. Ostoja-Starzewski, On the size of RVE in finite elasticity of random composites, J. Elasticity 85 (2) (2006) 153–173. [14] V. Kouznetsova, M.G. Geers, W.M. Brekelmans, Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme, Internat. J. Numer. Methods Engrg. 54 (8) (2002) 1235–1260. [15] V. Kouznetsova, M. Geers, W. Brekelmans, Size of a representative volume element in a second-order computational homogenization framework, Int. J. Multiscale Comput. Eng. 2 (4) (2004). [16] 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, Comput. Mech. 55 (5) (2015) 819–835. [17] H. Moussaddy, D. Therriault, M. Lévesque, Assessment of existing and introduction of a new and robust efficient definition of the representative volume element, Int. J. Solids Struct. 50 (24) (2013) 3817–3828. [18] V.P. Nguyen, O. Lloberas-Valls, M. Stroeven, L.J. Sluys, On the existence of representative volumes for softening quasi-brittle materials–a failure zone averaging scheme, Comput. Methods Appl. Mech. Engrg. 199 (45) (2010) 3028–3038. [19] M. Ostoja-Starzewski, Material spatial randomness: from statistical to representative volume element, Probab. Eng. Mech. 21 (2) (2006) 112–132. [20] M. Ostoja-Starzewski, J. Schulte, Bounding of effective thermal conductivities of multiscale materials by essential and natural boundary conditions, Phys. Rev. B 54 (1) (1996) 278. [21] S.I. Ranganathan, M. Ostoja-Starzewski, Scaling function, anisotropy and the size of RVE in elastic random polycrystals, J. Mech. Phys. Solids 56 (9) (2008) 2773–2791. [22] Z.-Y. Ren, Q.-S. Zheng, A quantitative study of minimum sizes of representative volume elements of cubic polycrystals–numerical experiments, J. Mech. Phys. Solids 50 (4) (2002) 881–893. [23] D. Savvas, G. Stefanou, M. Papadrakakis, Determination of RVE size for random composites with local volume fraction variation, Comput. Methods Appl. Mech. Engrg. 305 (2016) 340–358. [24] G. Stefanou, D. Savvas, M. Papadrakakis, Stochastic finite element analysis of composite structures based on material microstructure, Compos. Struct. 132 (2015) 384–392. [25] S. Swaminathan, S. Ghosh, Statistically equivalent representative volume elements for unidirectional composite microstructures: Part II-With interfacial debonding, J. Compos. Mater. 40 (7) (2006) 605–621. [26] K. Teferra, S.R. Arwade, G. Deodatis, Stochastic variability of effective properties via the generalized variability response function, Comput. Struct. 110 (2012) 107–115. [27] K. Teferra, S.R. Arwade, G. Deodatis, Generalized variability response functions for two-dimensional elasticity problems, Comput. Methods Appl. Mech. Engrg. 272 (2014) 121–137. [28] X. Yin, W. Chen, A. To, C. McVeigh, W.K. Liu, Statistical volume element method for predicting microstructure–constitutive property relations, Comput. Methods Appl. Mech. Engrg. 197 (43) (2008) 3516–3529. [29] J. Dirrenberger, S. Forest, D. Jeulin, Towards gigantic RVE sizes for 3D stochastic fibrous networks, Int. J. Solids Struct. 51 (2) (2014) 359–376.

270

K. Teferra, L. Graham-Brady / Comput. Methods Appl. Mech. Engrg. 330 (2018) 253–270

[30] M.G. Geers, V.G. Kouznetsova, W. Brekelmans, Multi-scale computational homogenization: Trends and challenges, J. Comput. Appl. Math. 234 (7) (2010) 2175–2182. [31] A. Gloria, F. Otto, et al., An optimal variance estimate in stochastic homogenization of discrete elliptic equations, Ann. Probab. 39 (3) (2011) 779–856. [32] A. Bensoussan, J.-L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structures, Vol. 374, American Mathematical Soc., 2011. [33] G. Geymonat, S. Müller, N. Triantafyllidis, Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity, Arch. Ration. Mech. Anal. 122 (3) (1993) 231–290. [34] A. Gloria, Stochastic diffeomorphisms and homogenization of multiple integrals, Appl. Math. Res. Express 2008 (abn001) (2008) 1–23. [35] A. Gloria, Qualitative and Quantitative Results in Stochastic Homogenization, Université des Sciences et Technologie de Lille-Lille I, 2012. [36] M. Grigoriu, Stochastic Systems: Uncertainty Quantification and Propagation, Springer Science & Business Media, 2012. [37] M. Grigoriu, Stochastic Calculus: Applications in Science and Engineering, Springer Science & Business Media, 2013. [38] R. Hill, Elastic properties of reinforced solids: some theoretical principles, J. Mech. Phys. Solids 11 (5) (1963) 357–372. [39] R. Hill, The essential structure of constitutive laws for metal composites and polycrystals, J. Mech. Phys. Solids 15 (2) (1967) 79–95. [40] P. Suquet, Elements of homogenization for inelastic solid mechanics, Homogenization Tech. Compos. Media 272 (1987) 193–278. [41] R.R. Stoll, Sets, Logic, and Axiomatic Theories, 1961. [42] V. Yurinskii, Averaging of symmetric diffusion in random medium, Sib. Math. J. 27 (4) (1986) 603–613. [43] L.A. Caffarelli, P.E. Souganidis, Rates of convergence for the homogenization of fully nonlinear uniformly elliptic pde in random media, Invent. Math. 180 (2) (2010) 301–360. [44] E. Vanmarcke, Random Fields: Analysis and Synthesis, World Scientific, 2010. [45] M. Shinozuka, G. Deodatis, Simulation of multi-dimensional Gaussian stochastic fields by spectral representation, Appl. Mech. Rev. 49 (1) (1996) 29–53. [46] R.L. Taylor, FEAP Finite Element Analysis Program, Ing.-Gemeinschaft Klee & Wrigges, 1987. [47] M.A. Groeber, M.A. Jackson, DREAM. 3D: a digital representation environment for the analysis of microstructure in 3D, Integr. Mater. Manuf. Innov. 3 (1) (2014) 1–17. [48] D. Peirce, R.J. Asaro, A. Needleman, Material rate dependence and localized deformation in crystalline solids, Acta Metall. 31 (12) (1983) 1951–1976. [49] K. Teferra, L. Graham-Brady, Tessellation growth models for polycrystalline microstructures, Comput. Mater. Sci. 102 (2015) 57–67. [50] S.N. Chiu, D. Stoyan, W.S. Kendall, J. Mecke, Stochastic Geometry and its Applications, John Wiley & Sons, 2013. [51] Y. Huang, A User-Material Subroutine Incroporating Single Crystal Plasticity in the ABAQUS Finite Element Program, Harvard Univ., 1991. [52] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, V. Eijkhout, W. Gropp, D. Kaushik, M. Knepley, (2014) et al., PETSc users manual revision 3.5, Argonne National Laboratory (ANL). [53] B. Efron, R.J. Tibshirani, An Introduction to the Bootstrap, CRC press, 1994.