Role of exterior statistics-based boundary conditions for property-based statistically equivalent representative volume elements of polydispersed elastic composites

Abstract

The property-based statistically equivalent RVE or P-SERVE has been introduced in the literature as the smallest microstructural volume element in non-uniform microstructures that has effective material properties equivalent to those of the entire microstructure. An important consideration is the application of appropriate boundary conditions for optimal property-based statistically equivalent representative volume element domains. The exterior statistics-based boundary conditions have been developed, accounting for the statistics of fiber distributions and interactions in the domain exterior to the property-based statistically equivalent representative volume element. This paper is intended to validate the efficacy of the exterior statistics-based boundary condition-based property-based statistically equivalent representative volume elements for evaluating homogenized stiffnesses of a unidirectional polymer matrix composite with a polydispersed microstructure characterized by nonuniform dispersion of carbon fibers of varying sizes in an epoxy matrix. Experimental tests and microstructural characterization of the polymer matrix composite are conducted for calibration and validation of the model. Statistically equivalent microstructural volume elements are constructed from experimental micrographs for direct numerical simulations. The performance of the property-based statistically equivalent representative volume element with exterior statistics-based boundary conditions is compared with other boundary conditions, as well as with the statistical volume elements. The tests clearly show the significant advantages of the exterior statistics-based boundary conditions in terms of accuracy of the homogenized stiffness and efficiency.

Keywords

Polymer matrix composites statistically equivalent representative volume element exterior statistic-based BC two-point correlation functions polydispersed microstructures

Introduction

Effective mechanical properties, like stiffness and strength, of composite materials are strongly affected by morphological features of the microstructure, such as fiber volume fraction, fiber size, shape and orientation, spatial dispersion of fibers, etc. These properties are typically evaluated by homogenization methods that involve averaging microscopic variables viz. stresses and strains from direct numerical simulations (DNS) of the microstructural representative volume element (RVE).^1–6 Typically, an RVE is a microstructural sub-domain that is representative of the morphology and/or properties of the entire microstructure in an averaged sense. For non-uniformly dispersed microstructures, as shown in Figure 1(a), RVE definitions that pertain to perfectly uniform or periodic microstructures may not apply. The concept of statistically equivalent RVE or SERVE has been introduced in Swaminathan et al.,¹ McDowell et al.⁷ and Ghosh² for non-uniform microstructures. Combining statistical analysis and DNS of representative microstructural domains, an SERVE can be identified as the smallest microstructural volume element exhibiting the following characteristics.

Effective material properties of the SERVE should be equivalent to those of the entire microstructure. This is classified as a property-based SERVE or property-based statistically equivalent (P-SERVE).

Distribution functions of parameters reflecting the statistics of local morphology in the SERVE should be equivalent to those for the overall microstructure. This can be classified as a microstructure-based SERVE or M-SERVE.

The SERVE should be independent of location in the local microstructure or loading.

Figure 1.

(a) Experimental micrograph of carbon fiber epoxy matrix unidirectional composite, (b) Voronoi tessellation of the polydispersed microstructural volume element (MVE) with gray-scale shading corresponding to the local volume fraction Φ. Coordinates x (or 1) and y (or 2) correspond to transverse directions in the section, while z (or 3) corresponds to the longitudinal direction of the fibers.

Figure 1(a) shows an experimental micrograph of an unidirectional carbon fiber, epoxy polymer matrix composite. A microstructural volume element (MVE) is a large microstructural domain representing the micrograph, for which statistical distributions of characteristic functions can be generated, but DNS is computationally prohibitive. Figure 1(b) shows an MVE of the polydispersed microstructure that is tessellated into Voronoi cells for determining the distribution of volume fraction Φ as shown with gray-scale shading. The goal of this paper is to devise and test an efficient and accurate method of determining effective properties of non-uniform composite microstructures. A variety of methods with these goals exist in the literature, where DNS of estimated SERVEs are performed subjected to certain boundary conditions. Special RVEs and homogenization methods for random media have been developed in Kanit et al.,⁸ Trovalusci et al.⁹ and Recciaa et al.¹⁰ Conventional boundary conditions include the constant strain or affine-transformation-based displacement (ATDBCs), constant stress or uniform traction (UTBCs) and periodic boundary conditions (PBCs).^11–13 Assumptions in these boundary conditions either imply that the SERVE is immersed in a continuum with spatially invariant strain-energy-density, ignoring the interaction of exterior fibers, or that the deformation patterns in the domain exterior to the SERVE are homologous. These assumptions are however not valid for microstructures with non-uniform dispersions and can result in over-estimation of the RVE region from convergence requirements. To overcome the deficiencies in boundary conditions, the exterior statistics-based boundary condition or ESBC has been developed for P-SERVEs in Ghosh and Kubair³ and Kubair and Ghosh,⁴ accounting for statistics of fiber distributions and interactions in the MVE domain exterior to the P-SERVE. Using statistically informed Green’s functions and the Eshelby equivalent inclusion method and the two-point correlation function $S_{2} (r, θ)$ of the MVE, ESBCs are derived to describe the interactions of exterior domain with the P-SERVE. The resulting improved boundary conditions lead to the reduced P-SERVE size and efficient computations.

This paper is intended to validate the efficacy of the ESBC-based P-SERVEs for a unidirectional polymer matrix composite (IM7/977-3 PMC) with a polydispersed microstructure containing a nonuniform distribution of IM7 carbon fibers of varying sizes in a 977-3 epoxy matrix. Only elastic materials are considered in this paper. The second section discusses experimental tests and microstructural characterization of the PMC. From experimental micrographs, statistically equivalent MVEs of the material are constructed to provide microstructural volumes for DNS in the section on microstructure reconstruction. The next section summarizes the formulation and algorithms for micromechanical analysis of the P-SERVE with exterior statistics-based boundary conditions (ESBCs). The selection of the P-SERVE with ESBCs from the convergence of effective properties is discussed in the following section. Next, the performance of the P-SERVE with ESBC is compared with the statistical volume elements (SVEs). The predicted effective stiffness by different approaches is compared with the results of mechanical tests. The paper concludes with a summary of the advantages achieved by the ESBC-based P-SERVE.

Microstructure imaging and characterization, and mechanical testing of the carbon fiber PMC

The unidirectional polymer matrix composite (IM7/977-3 PMC) studied here, has a composition of IM7 carbon fibers of varying sizes (mean radius of $2.43 μ m$ and standard deviation of $0.11 μ m$ ) nonuniformly distributed in a 977-3 epoxy matrix. Due to unidirectionality, microstructural sections perpendicular to the fiber direction are cut with a lubricated diamond saw blade and mounted in SamplKwick acrylic. After sufficient cure time, the sample is polished in 9 µm, 3 µm, and 1 µm diamond slurries for 20 min each, followed by a final polishing in 0.3 µm alumina slurry for 30 min. The sample is subsequently cleaned for imaging using acetone and isopropanol. Microstructural images are taken using a Zeiss Axio Observer.Z1 inverted microscope under bright-field light.

A 200 µm × 250 µm montage, acquired by the software autofocus with a

500 \times

objective lens, is depicted in Figure 1(a). The overall fiber volume fraction of this microstructure is

\sim 62 %

. Mechanical tests on the PMC are performed following the ASTM-D3039 standardized test protocols and have been reported in Clay and Knoth.¹⁴ The elastic properties of the IM7/977-3 polymer matrix composite are selectively documented in Table 1. The effective stiffness coefficient can be calculated from the values given in Table 1. For example, the transverse stiffness component

{\bar{C}}_{1111}^{exp} = 11.4 \pm 0.3597 GPa

, which is of the same order as

E_{2 T}

Table 1.

Experimentally obtained elastic properties of the IM7/977-3 polymer matrix composite Clay and Knoth.¹⁴

Property	Mean	St. Dev.
Tens. longitudinal modulus ( $E_{1 T}$ )	164 GPa	4.12
Comp. longitudinal modulus ( $E_{1 C}$ )	137 GPa	0.608
Tens. transverse modulus ( $E_{2 T}$ )	8.98 GPa	0.284
Compressive transverse modulus ( $E_{2 C}$ )	8.69 GPa	0.367
In-plane shear modulus (G₁₂)	5.01 GPa	0.249
Major Poisson’s ratio (ν₁₂)	0.320	0.0266

For micromechanical analysis in the subsequent sections, material properties of the microstructural constituents have been obtained from various sources. The Young’s modulus and Poisson’s ratio of the isotropic 977-3 epoxy matrix are E^M= 2.5 GPa and ν^M= 0.43, respectively. The IM7 carbon fibers are assumed to be transversely isotropic, for which the modulus in the longitudinal direction is recorded in the manufacturer database HEXCEL¹⁵ as $E_{long}^{F}$ = 276 GPa. The transverse direction modulus of the IM7 carbon fibers has been experimentally obtained in the Sottos group Montgomery and Sottos¹⁶ as $E_{trans}^{F}$ = 19.5 GPa, while the Poisson’s ratio is ν^F= 0.23.

Statistical characterization of the polydispersed microstructure

An image-processed micrograph of the cross-section of the unidirectional PMC is shown in Figure 1(a). A large region in the micrograph is designated as the MVE and tessellated into a network of Voronoi cells,² based on the fiber centroids, as shown in Figure 1(b). The selected MVE consists of 1239 fibers with a median fiber radius of 2.4624 µm. Voronoi cells provide a basis for microstructural characterization. The local fiber volume fraction Φ is defined as the ratio of the fiber cross-sectional area to the area of the associated Voronoi cell, and the shading represents the level of Φ. Brighter cells with lower values of Φ indicate regions that are matrix rich, while darker cells with large Φ indicate regions of fiber clustering.

The probability density functions (PDFs) of the local volume fraction Φ and the fiber size of the MVE are, respectively, plotted in Figures 2(a) and (b). The median volume fraction is evaluated to be $Φ = 0.63$ . The distribution of the normalized two-point correlation function $\frac{S_{2} (r, θ)}{S_{1}^{2}}$ is shown in Figure 2(c), where S₁ is the one-point correlation function corresponding to the overall volume fraction, and $S_{2} (r, θ)$ is the radial-distance and orientation dependent, two-point correlation function.¹⁷ The function $\frac{S_{2} (r, θ)}{S_{1}^{2}}$ has a peak near the center and reaches a far-field value of $S_{1}^{2}$ with some oscillations.

Figure 2.

Probability distribution functions of: (a) volume fraction Φ and (b) fiber size, and (c) the radial-distance and orientation dependent two-point correlation function $S_{2} (r, θ)$ of the polydispersed MVE. MVE: microstructural volume element.

Creating statistically equivalent MVEs from experimental micrographs

Prior to conducting DNS of the SERVE for effective properties, it is important to create statistically equivalent microstructural volume elements (SEMVEs) that are representative of the entire experimental micrograph as in Figure 1(a) with identical statistical distributions. The SEMVE forms the parent domain from which P-SERVEs may be extracted.

In this work, SEMVEs are obtained by a Monte-Carlo type reconstruction method using the one-point correlation function S₁, two-point correlation function S₂ and the probability distribution of fiber radius in the PMC micrograph, depicted in Figure 2. The process initializes the positions of N_f chosen fibers in hexagonal close packing arrangements within the MVE, with the initial spacing between fibers determined from the global volume fraction S₁. The radii of these fibers are assigned by randomly sampling from the fiber size distribution in Figure 2(b). Values of parameters used to start the process are: $N_{f} \approx 2000, S_{1} = 0.65$ , mean fiber radius $= 2.43 μ m$ and standard deviation of $= 0.11 μ m$ . The reconstruction process executes a large number of iterations ( $\sim 1000 N_{f}$ ). The centroid of a randomly selected fiber is perturbed in each iteration to transform the hexagonal arrangement to an amorphous microstructure. Non-overlapping perturbations are accepted for the first $\sim 60 %$ of the iterations towards this amorphization. The subsequent 40% iterations optimize the configuration of the SEMVE by minimizing the L₂-norm of the difference between the SEMVE instantiation and the reference MVE.

\begin{matrix} | | S_{2}^{exp} (r, θ) - S_{2}^{SEMVE} (r, θ) | |_{L_{2}} \\ = \sqrt{\int_{Ω^{MVE}} (S_{2}^{exp} (r, θ) - S_{2}^{SEMVE} (r, θ))^{2} d Ω} \end{matrix}

(1)

After each proposed non-overlapping fiber shuffle, $S_{2} (r, θ)$ of the SEMVE is evaluated and the resulting L₂-norm is compared with that for the unperturbed state. The move is accepted only if there is an improvement with respect to the L₂-norm. Figure 3(a) shows the convergence in $S_{2} (r, θ)$ as a function of iterations. The rate of convergence is rapid in the beginning of the amorphization. At $\sim 60 %$ of iterations, the optimization algorithm begins with a sharp increase in convergence rate and stabilizes at around 70% of the iterations. The map of error in the $S_{2} (r, θ)$ field between the experimental MVE and the reconstructed final SEMVE of Figure 3(c) is given in Figure 3(b). Figure 3(c) shows the square window containing the SEMVE represented by the fibers.

Figure 3.

(a) Convergence with respect to the L₂-norm of error with progressive iterations, and (b) map of error in the $S_{2} (r, θ)$ field between the experimental MVE and reconstructed SEMVEs, and (c) square window containing the final SEMVE microstructure with $N_{f} \sim 2000$ fibers. MVE: microstructural volume element; SEMVE: statistically equivalent microstructural volume element.

Micromechanical analysis of the P-SERVE with ESBCs

An essential step in the property evaluation process is the identification of the property-based SERVE and the associated boundary conditions. The ESBCs have been developed in Ghosh and Kubair³ and Kubair and Ghosh⁴ as a special conjugate to the optimal P-SERVE size. In Ghosh and Kubair,³ the ESBC developed was for statistically homogeneous monodispersed MVEs that allowed the use of the distance-based two-point correlation functions $S_{2} (r)$ in its formulation. The ESBC model in Kubair and Ghosh⁴ extended this to statistically nonhomogeneous (clustered) monodispersed MVEs and introduced the joint distance and orientation-based two-point correlation function $S_{2} (r, θ)$ in its development. The present paper further extends the scope of ESBCs to statistically nonhomogeneous (clustered), polydispersed microstructures. The detailed formulation in Ghosh and Kubair³ and Kubair and Ghosh⁴ is summarized here.

The polydispersed MVE with nonuniformly distributed fibers, for which homogenized properties are sought, occupies a large microscopic domain $Ω^{mve} \to Ω^{\infty}$ as shown in Figure 4(a). The MVE boundary value problem in $Ω^{mve}$ is partitioned into the complementary SERVE domain $Ω^{serve}$ and a domain $Ω^{ext}$ exterior to it, i.e. $Ω^{mve} = Ω^{ext} \cup Ω^{serve}$ . The effect of the exterior domain $Ω^{ext}$ is manifested through equivalent conditions on the SERVE boundary $Γ^{serve}$ that adequately reflects the interaction of $Ω^{ext}$ with $Ω^{serve}$ . It should result in the same invariant strain energy for the SERVE as for the entire MVE with the applied affine displacement conditions $u_{i}^{A} (x^{\infty}) = ε_{ij}^{0} x_{j}^{\infty}$ on $Γ^{\infty}$ . The principle of virtual work is written as the sum of the respective virtual work terms in the complementary domains in Figure 4(a) as

\begin{matrix} \int_{Ω^{ext}} σ_{ij}^{ext} (x) δ ε_{ij}^{ext} (x) d Ω + \int_{Ω^{serve}} σ_{ij}^{serve} (x) δ ε_{ij}^{serve} (x) d Ω = 0 \end{matrix}

(2)

Figure 4.

(a) Schematic view of the MVE containing the P-SERVE and its complementary exterior domain, i.e. $Ω^{mve} = Ω^{serve} \cup Ω^{ext}$ , and (b) effect of an interacting fiber-pair I-J on a field point O at the P-SERVE boundary $Γ^{serve}$ . MVE: microstructural volume element; P-SERVE: property-based statistically equivalent RVE.

Applying the divergence theorem to the first term containing the integral over $Ω^{ext}$ , incorporating equilibrium conditions in the absence of body forces and with ${\bar{ε}}_{ij}^{mve} = ε_{ik}^{0}$ on $Γ^{\infty}$ (where ${\bar{ε}}_{ij}^{mve} = á ε_{ij}^{mve} (x) ñ$ is the homogenized strain in the MVE), the principle of virtual work (2) reduces to that of the SERVE as

\begin{matrix} \int_{Ω^{serve}} σ_{ij}^{serve} (x) δ ε_{ij}^{serve} (x) d Ω - \int_{Γ^{serve}} T_{i}^{ext} (x) δ u_{i}^{ext} (x) d Γ = 0 \end{matrix}

(3)

$T_{i}^{ext} (x)$ is the traction on $Γ^{serve}$ resulting from the stresses in the domain $Ω^{ext}$ exterior to the SERVE. The second term in equation (3) will drop out if an effective displacement field can be prescribed on $Γ^{serve}$ , since $δ u_{i}^{ext} = 0$ on $Γ^{serve}$ . This can be incorporated by augmenting the affine transformation-based boundary displacement field $u_{i}^{A} (x^{serve}) = ε_{il}^{0} x_{l}^{serve}$ with an additional perturbation term, which represents the effects of the heterogeneous $Ω^{ext}$ on $Γ^{serve}$ . The solution process will not involve an explicit numerical solution of the exterior domain problem in $Ω^{ext}$ . Hence, a special analytical solution that involves the statistics of the exterior domain is developed in the form of ESBCs or $u_{i}^{ESBC}$ , that are expressed as

u_{i}^{ESBC} (x^{serve}) = u_{i}^{A} (x^{serve}) + u_{i}^{*} (x^{serve}) on Γ^{serve}

(4)

Here $u_{i}^{*}$ is an enhancement due to the interaction of the non-uniform exterior domain $Ω^{ext}$ with the SERVE.

The method of obtaining the displacement enhancement $u_{i}^{*} (x)$ (equation (4)) due to the interaction with the exterior domain is briefly described. The presence of heterogeneities in the form of inclusions or fibers alters the spatially invariant, homogeneous fields of matrix stress $σ_{ij}^{M}$ , strain $ε_{ij}^{M}$ and displacement $u_{i}^{M}$ in $Ω^{mve}$ by the respective perturbed fields of stress $σ_{ij}^{*}$ , strain $ε_{ij}^{*}$ and displacement $u_{i}^{*}$ that depend on the morphological characteristics of the microstructure. The total stress, strain and displacement fields in the MVE are defined as the sum of the homogeneous and perturbed parts as

\begin{matrix} σ_{ij} (x) = σ_{ij}^{M} + σ_{ij}^{*} (x), ε_{ij} (x) = ε_{ij}^{M} + ε_{ij}^{*} (x) \\ u_{i} (x) = u_{i}^{M} + u_{i}^{*} (x) \in Ω^{mve} \end{matrix}

(5)

Since the homogeneous stress $σ_{ij}^{M}$ is divergence-free, the equilibrium condition for the perturbed stress fields in the absence of body forces is $σ_{ij, j}^{*} (x) = 0$ . The solution to the problem of a heterogeneous medium is simplified to that of a homogeneous medium by introducing an equivalent inclusion approach, where an eigenstrain $ε_{ij}^{Λ} (x)$ is introduced in the inclusion domain to account for the constraint that the matrix imposes on the inclusion due to autonomous deformation. The eigenstrains with n_p interacting inclusions are obtained by statistically informed Green’s functions as detailed in Ghosh and Kubair³ and Kubair and Ghosh,⁴ as well as by applying the Eshelby’s stress consistency condition, which requires the total stress inside the inclusion $Ω^{F}$ to be equal to the total stress in the equivalent matrix domain.

For the domain $Ω^{mve}$ consisting of polydispersed fibers, the distribution of fiber radius and spatial location is represented by the PDF of the fiber-size $PDF (a)$ and the two-point correlation function $S_{2} (r)$ , where $r = x - x^{I}$ is the position vector separating the source and field points in the domain. This can be represented in a parametric form as $(r, θ)$ , where the parameter $r = | r |$ is the separation distance and $θ = ∠ r$ is the orientation of the line joining these points with a reference direction. The eigenstrain $ε_{ij}^{Λ}$ in a reference fiber occupying a domain $Ω^{F}$ is written as

\begin{matrix} ε_{ij}^{Λ} (x) = [ι^{F} (x) (S_{ijab} + M_{ijab}) - \\ \int_{Ω^{mve} \ Ω^{F}} PDF (a) S_{2} (r) {\hat{G}}_{ijmn} (r, a) (S_{mnpq} + M_{mnpq})^{- 1} \\ {\hat{G}}_{pqab} (r, a) d Ω]^{- 1} [((S_{abmn} + M_{abmn})^{- 1} \end{matrix}

(6)

\begin{matrix} \int_{Ω^{mve} \ Ω^{F}} PDF (a) S_{2} (r) {\hat{G}}_{mnkl} (r, a) d Ω) - \\ \frac{1}{2} (δ_{ak} δ_{bl} + δ_{al} δ_{bk})] ε_{kl}^{M} = A_{ijkl} (x) ε_{kl}^{M} \forall x \in Ω^{mve} \end{matrix}

S_ijkl are the spatially invariant interior, ${\hat{G}}_{ijkl} (r)$ are the position dependent exterior Eshelby tensors, respectively, described in Ghosh and Kubair³ and Kubair and Ghosh⁴ and ι is a phase-based unit function. Furthermore, $M_{ijkl} = (C_{ijpq}^{F} - C_{ijpq}^{M})^{- 1} C_{pqkl}^{M}$ , with $C_{ijkl}^{M}$ and $C_{ijkl}^{F}$ being the elastic stiffness of the matrix and inclusion materials, respectively. The perturbed displacements at an observation point O in Figure 4(b) can be written in terms of the matrix strain as

\begin{matrix} u_{i}^{*} (x) = [\int_{Ω^{mve} \ Ω^{F}} PDF (a) S_{2} (r') \\ \times L_{imn} (r', a) A_{mnkl} (r') d Ω] ε_{kl}^{M} \end{matrix}

(7)

where

L_{ikl} (r)

is a unified displacement-transfer tensor, described in Ghosh and Kubair³ and Kubair and Ghosh.⁴ Finally, using equation (4), the affine transformation-based displacement fields are superposed on the perturbed displacements to prescribe the ESBCs. The ESBCs are applied on the boundary

Γ^{serve}

of a P-SERVE domain

Ω^{serve}

of size L using the following steps.

Discretize the P-SERVE domain $Ω^{serve}$ into a finite element mesh, e.g. with 4-noded tetrahedral elements in 3D.

Extract the positions and coordinates x_i of candidate boundary-nodes on $Γ^{serve}$ .

Compute the affine transformation-based displacements $u_{i}^{A} (x)$ on all boundary nodes with the applied far-field strain $ε_{ij}^{0}$ as $u_{i}^{A} (x) = ε_{ij}^{0} x_{j}$ where x_j is measured relative to the centroid of the P-SERVE.

Compute the probability distribution function PDF(a) and the two-point correlation function $S_{2} (r, θ)$ for the MVE domain $Ω^{mve}$ .

Compute the perturbed displacements $u_{i}^{*}$ using equation (7) incorporating PDF(a) and $S_{2} (r, θ)$ for all the boundary nodes, using the following steps.

Each radial orientation is discretized into N_r number of equally spaced segments with increment $Δ r = \frac{R - a}{N_{r}}$ , where a is the fiber radius, R is the radius of horizon (extent of the MVE) and the lower limit of the integration is r = a. The αth radial point is given as $r_{α} = α \frac{R - a}{N_{r}}$ .

The angular orientation is discretized into $N_{θ}$ equally spaced points of $Δ θ = \frac{2 π}{N_{θ}}$ . The βth angular point is $θ_{β} = β \frac{2 π}{N_{θ}}$ .

The fiber size distribution is discretized into N_a equally spaced bins with $Δ a = \frac{a_{max} - a_{min}}{N_{a}}$ , a_max and a_min being the maximum and minimum fiber size. At a SERVE boundary node at x_i, the discrete perturbed displacement components in equation (7) are evaluated for an applied strain $ε_{ij}^{0}$ as

\begin{matrix} u_{i}^{*} (x) = [\sum_{γ = 1}^{N_{a}} \frac{2 π (R - γ Δ a)}{N_{r} N_{θ}} \sum_{α = 1}^{N_{r}} \\ \sum_{β = 1}^{N_{θ}} α L_{imn} (x - (α Δ r, β Δ θ, γ Δ a)) \\ \times A_{mnkl} (x - (α Δ r, β Δ θ)) S_{2} (x - (α Δ r, β Δ θ)) \\ \times PDF (γ Δ a)] ε_{ij}^{0} \end{matrix}

(8)

The ESBCs on the boundary nodes are computed and applied as

u_{i}^{ESBC} (x) = u_{i}^{A} (x) + u_{i}^{*} (x)

(9)

Selection of a candidate P-SERVE from stiffness convergence

Candidate P-SERVEs are extracted from the SEMVE domain in Figure 3(c) for simulations leading to the homogenized stiffness evaluation. Figure 5(a) shows a set of five concentric cross-sections (i–v) with increasing number of fibers in the SEMVE that can be used as candidate P-SERVEs. The P-SERVE boundaries coincide with the Voronoi cell boundaries at the edges of the MVE, and hence does not intersect any fiber. The thickness of the composite domai $Γ^{serve}$ n is considered to be $10 μ m$ . The FE model is discretized into 4-noded tetrahedral elements with 10 elements in the z-direction.

Figure 5.

(a) Concentrically increasing candidate P-SERVE domains in the MVE for micromechanical simulations, and (b) convergence of the normalized homogenized stiffness tensor ${\bar{C}}_{1111} / E^{M}$ with P-SERVE size. MVE: microstructural volume element; P-SERVE: property-based statistically equivalent RVE.

The composite system is assumed to be the weakest in transverse loading, as corroborated by the transverse modulus in Table 1. Hence, the analyses are performed in the transverse direction. The P-SERVEs are subjected to both the ATDBCs or ESBCs corresponding to a far-field unit uniaxial strain in the transverse direction $ε_{11}^{0} = 1$ . The axes 1 (or x) and 2 (or y) in Figure 5(a) correspond to the transverse directions and the fiber axis is in the 1 (or(z) direction. All other strain components are set to zero. 3D finite element simulations of the P-SERVEs are performed and the homogenized stiffness ${\bar{C}}_{ijkl}, i, j, k, l = 1, 2, 3$ is evaluated by averaging the stress and strain fields as described in Ghosh and Kubair³ and Kubair and Ghosh⁴. The analyses represent a tensile loading in the critical transverse direction, for which the tensile strength is the least and the microstructure is most suscept $E_{2 T}$ ible to cracking. Convergence in the homogenized stiffness with increasing P-SERVE size is used as a metric to determine the optimal P-SERVE size. In this study, the critical stiffness component ${\bar{C}}_{1111}$ is used to determine the effectiveness of the applied boundary conditions on the converged P-SERVE size.

The homogenized stiffness component ${\bar{C}}_{1111}$ normalized by the matrix stiffness is plotted as a function of the P-SERVE size L in Figure 5(b). The figure shows that the homogenized modulus obtained with the ESBCs (superscript ESBC) converges at a SERVE size of approximately $L = 70 μ m$ consisting of 179 fibers. In contrast, a much larger P-SERVE size of approximately $L \approx 200$ µm is necessary for convergence when subjected to the ATDBC.

Comparing P-SERVEs with ESBCs with SVEs

SVEs^7,18 are based on the ergodicity hypothesis that the composite microstructure with dispersed heterogeneities is statistically homogeneous. Hence, its volume-averages are identical to the ensemble-averages. The homogenized modulus of the material is expected to be equal the mean of the volume-averaged modulus obtained from a large number of instantiations of a much smaller volume SVEs. The ensemble-average of any spatially varying field quantity $Ψ (x)$ over the SVE domain may thus be expressed in terms of the homogenized value over the MVE as

\bar{Ψ} = \frac{1}{Ω^{mve}} \int_{Ω^{mve}} Ψ (x) d Ω = \frac{1}{N} \sum_{I = 1}^{I = N} (\frac{1}{Ω^{{sve}_{I}}} \int_{Ω^{{sve}_{I}}} Ψ (x) d Ω)

(10)

where

\bar{Ψ}

is the volume-averaged value,

Ω^{{sve}_{I}}

is the domain of the I^th SVE instantiation and N corresponds to the number of sample SVEs in the ensemble. In general, the volume of any SVE is much smaller than that of the SERVE i.e.

Ω^{{sve}_{I}} < < Ω^{serve}

For comparison with the P-SERVE predictions, the SVE problem is set up with individual square SVEs of size L^I= 70 µm, 135 µm and 200 µm; 100 candidate SVEs are chosen for each SVE size. 2D plane-strain analysis of the candidate SVEs is performed subjected to ATDBCs. The ensemble-averaged stiffness components ${\bar{C}}_{ijkl}$ are obtained for the population of N SVEs as: ${\bar{C}}_{ijkl} = (\sum_{I = 1}^{I = N} {\bar{C}}_{ijkl}^{I}) / N$ , where ${\bar{C}}_{ijkl}^{I}$ are the volume-averaged stiffness components for the $I -$ th SVE. With increasing number of instantiations in the ensemble population N, the ensemble averaged stiffness components are expected to converge to their respective homogenized values for the MVE ( ${\bar{C}}_{ijkl}^{\infty}$ ). The convergence criterion is defined in terms of the minimum number of instantiations or SVEs N required in the ensemble to attain a steady-state, invariant value of the homogenized stiffness. Convergence is ascertained from the plot of the cumulative mean (CM) of the normalized stiffness as a function of the ensemble population size N, as shown in Figure 6(a). The cumulative mean of a stiffness component ${\bar{C}}_{ijkl}$ , normalized by the matrix Young’s modulus E^M, is defined as $CM (\frac{{\bar{C}}_{ijkl}}{E^{M}})_{N} = \frac{1}{N} \sum_{I = 1}^{I = N} \frac{{\bar{C}}_{ijkl}^{I}}{E^{M}}$

Figure 6.

(a) Cumulative mean (CM) of the ensemble-averaged stiffness component ${\bar{C}}_{1111}$ as a function of the number of SVEs for different SVE sizes (L^sve); (b) error in CM of the stiffness as a function of the number of SVEs. SVEs: statistical volume elements.

For an ergodic microstructure, the cumulative mean (CM) of the volume-averaged stiffness is expected to converge to that of the entire MVE, i.e. ${\bar{C}}_{ijkl}^{\infty}$ . The CM of ${\bar{C}}_{1111}$ , obtained from the three SVE sizes, is shown in Figure 6(a). For SVE size $L^{sve} = 70$ µm, the ensemble averaged stiffness even with 100 instantiations does not converge to the accurate value obtained from the P-SERVE analysis with ESBCs in a single iteration. The corresponding error in CM, defined as the difference from the stiffness ${\bar{C}}_{1111}^{\infty}$ , is plotted as a function of the number of SVEs in Figure 6(b). This exemplifies the lack of convergence for smaller SVE sizes. The P-SERVE with ESBCs requires only one instantiation of size L^serve 70 µm for convergence, as illustrated in Figure 5.

Comparing the P-SERVE and SVE stiffness with experimental observations

Results of mechanical tests in Clay and Knoth¹⁴ have been discussed in the section on characterization and testing. The dominant stiffness coefficient is evaluated as

{\bar{C}}_{1111}^{exp} = 11.4 \pm 0.3597

GPa. The homogenized stiffness is obtained from finite element simulations of the P-SERVEs subjected to both ATDBCs and ESBCs, and SVEs subjected to ATDBCs. All the cases in their respective converged situations yield

{\bar{C}}_{1111}^{serve} = {\bar{C}}_{1111}^{sve} = 11.9 GPa

, which is within the acceptable range of the experimental error. The conditions for convergence of the different cases are given in Table 2. The P-SERVEs with ESBCs require the smallest volume

Ω^{serve}

for convergence with only one iteration. This is a significant advantage over the other methods proposed in the literature.

Table 2.

Comparison of convergence conditions for different cases.

	Bound. Cond.	Size	$#$ instantiations
P-SERVE	ATDBC	200 µm	1
SVE	ATDBC	200 µm	10
P-SERVE	ESBC	70 µm	1

P-SERVE: property-based statistically equivalent RVE; ATDBC: affine-transformation-based displacement constant stress; ESBC: exterior statistics-based boundary condition.

Summary and conclusions

This applicability of the P-SERVE with ESBCs, introduced in Ghosh and Kubair³ and Kubair and Ghosh⁴ is extended to polydispersed microstructures in this paper. This P-SERVE-ESBC model reduces micromechanical simulations to an optimal SERVE domain through the use of boundary conditions that incorporate the statistics of the domain exterior to the SERVE. A premise of the efficiency of this method is that the micromechanical simulations of a small P-SERVE are efficient. The large exterior region of the microstructures only needs characterization for statistical analysis, which is rather efficiently accomplished for most material systems. An important aspect is the validation of this P-SERVE-ESBC model for a polymer matrix composite system. The validation results clearly show the significant advantage and potential of this method, both in terms of the volume to be modeled for determining effective mechanical properties and the number of iterations. While the advantage of ESBC-based method over other methods for homogenization has been proven for a few composite materials in this paper and in Ghosh and Kubair³ and Kubair and Ghosh,⁴ the choice of the appropriate method, in general, is material dependent and depends on the contrast in properties between the inclusion and matrix. This aspect has been discussed in Ranganathan and Ostoja-Starzewski.¹⁹

The experimental polymer matrix composite, on which microstructural characterization and mechanical tests are conducted, consists of unidirectional carbon fibers of varying radii dispersed in an epoxy matrix. SEMVEs of the material are constructed from equivalence of morphological distribution and correlation functions, e.g. $PDF (a), S_{1}, S_{2} (r, θ)$ , obtained from the actual micrograph images. The P-SERVEs are extracted from these volumes for micromechanical simulations with ESBCs and also the affine transformation-based displacement boundary conditions (ATDBCs). For convergence to the experimentally measured stiffness values, the P-SERVE with ESBCs requires a single simulation of $Ω^{serve} \approx 70 μ m$ , as opposed to $Ω^{serve} \approx 200 μ m$ with ATDBCs. In comparison, the statistical volume element (SVE) method needs a large number of iterations and instantiations with larger sample sizes ( $Ω^{{sve}_{I}} \approx 200 μ m$ ) for convergence. These validation tests show the significance of the P-SERVE with ESBCs for modeling elastic problems of heterogeneous materials. The authors are working on non-trivial extension of this method to nonlinear problems.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been supported through a grant No. FA9550-12-1-0445 to the Center of Excellence on Integrated Materials Modeling (CEIMM) at Johns Hopkins University awarded by the AFOSR/RSL Computational Mathematics Program (Manager Dr. A. Sayir) and AFRL/RX (Monitors Drs. C. Woodward and C. Przybyla). These sponsorships are gratefully acknowledged. Computing support by the Homewood High Performance Compute Cluster (HHPC) and Maryland Advanced Research Computing Center (MARCC) is gratefully acknowledged.

References

Swaminathan

Ghosh

Pagano

. Statistically equivalent representative volume elements for unidirectional composite microstructures: part I-without damage. J Comp Mater 2006; 40: 583–604.

Ghosh

. Micromechanical analysis and multi-scale modeling using the Voronoi cell finite element method, Boca Raton, Fl: CRC Press/Taylor & Francis, 2011.

Ghosh

Kubair

. Exterior statistics based boundary conditions for representative volume elements of elastic composites. J Mech Phys Solids 2016; 96: 1–24.

Kubair

Ghosh

. Exterior statistics based boundary conditions for establishing statistically equivalent representative volume elements of statistically nonhomogeneous elastic microstructures. Int J Solids Struct 2017; 112: 106–121.

Hill

. Elastic properties of reinforced solids: some theoretical principles. J Mech Phys Solids 1963; 11: 357–372.

Hashin

Shtrikman

. A variational approach to the theory of the elastic behaviour of multiphase materials. J Mech Phys Solids 1963; 11: 127–140.

McDowell

Ghosh

Kalidindi

. Representation and computational structure-property relations of random media. J Miner Metals Mater Soc 2011; 63: 45–51.

Kanit

Forest

Galliet

et al.

Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solids Struct 2003; 40: 3647–3679.

Trovalusci

Ostoja-Starsewski

De Bellis

et al.

Scale-dependent homogenization of random composites as micropolar continua. Eur J Mech 2015; 49: 396–407.

10.

Recciaa

De Bellis

Trovalusci

et al.

Sensitivity to material contrast in homogenization of random particle composites as micropolar continua. Compos Part B 2018; 136: 39–45.

11.

Hazanov

Huet

. Order relationships for boundary conditions effect in heterogeneous bodies smaller than the representative volume. J Mech Phys Solids 1994; 42: 1995–2011.

12.

Zohdi

Wriggers

. An introduction to computational micromechanics, Berlin, Germany: Springer, 2004.

13.

Ostoja-Starzewski

. Microstructural randomness and scaling in mechanics of materials, Boca Raton, Fl: Chapman and Hall/CRC, 2007.

14.

Clay

Knoth

. Experimental results of quasi-static testing for calibration and validation of composite progressive damage analysis methods. J Compos Mater 2016; 10: 1333–1353.

15.

HEXCEL. Composite materials and structures, www.hexcel.com (2017, accessed 6 January 2018).

16.

Montgomery C and Sottos N. Experiments for properties of composites. Unpublished work (private communication), 2017.

17.

Jiao

Stillinger

Torquato

. Modeling heterogeneous materials via two-point correlation functions: basic principles. Phys Rev E 2007; 76: 031110–031110.

18.

Yin

McVeigh

et al.

Statistical volume element method for predicting microstructureconstitutive property relations. Comp Meth Appl Mech Eng 2008; 197: 3516–3529.

19.

Ranganathan

Ostoja-Starzewski

. Scaling function, anisotropy and the size of RVE in elastic random polycrystals. J Mech Phys Solids 2008; 56(2773): 2791–2791.