MESO-SCALE HOMOGENIZATION METHODS FOR MOVING WINDOW ANALYSIS

Computational methods are developed to capture the macro-scale behavior of a heterogeneous composite material when the macro-scale is significantly larger than the scale of the heterogeneity. In this case, it may be computationally costly to model the material microstructure directly in a finite element analysis. Often a representative volume element (RVE) can be used with success; however, there are two main disadvantages of an RVE approach. First, when nonlinear behavior or damage initiation is considered, too much information is lost regarding the local variation in the microstructure. Second, an RVE approach effectively eliminates randomness in the meso-scale material representation, which is needed to statistically quantify structural reliability. As an alternative to the RVE approach, a moving windowing homogenization approach has been developed. This captures a degree of small scale material heterogeneity and randomness, and is suited to a priori implementation in a reduced-mesh large scale finite element analysis. However, non-unique solutions may exist for the homogenization of a window (or statistical volume element, SVE) that is below the scale of an RVE. At a scale smaller than the RVE, the inverse of the compliance tensor obtained by a statically uniform boundary condition (SUBC) test does not equal the stiffness tensor obtained by a kinematically uniform boundary condition (KUBC) test; these are lower and upper bonds on the effective behavior of the element, respectively. Mixed uniform boundary conditions and periodic boundary conditions, such as those used in the Generalized Method of Cells (GMC) predict apparent properties between the KUBC and SUBC upper and lower bounds. This work explores the assumption that accurately capturing effective behavior of a composite at the meso-scale and accurately predicting macro-scale behavior in a moving window analysis are equivalent considerations. A composite beam model with randomly varying inclusions is considered, and the deflection at midspan is used as a metric for macro-scale behavior. Homogenization methods including those described above are implemented. Results are compared to determine not only how well each homogenization technique approximates effective behavior, but also by how each predicts the macro-scale behavior of a structure in a finite element analysis.


INTRODUCTION
Digital imaging is a common way to obtain a representation of a composite material.The microstructural morphology is captured in a pixelized array of data.For example, the image shown in Figure 1 shows a basaltic magma obtained in geological studies [22].The full 3D data set describing this material is a 360 3 binary array, which describes a 1cm 3 volume.To enter this data into a finite element analysis with each voxel represented by an element would require a stiffness matrix with degrees of freedom on the order of 10 9 .Even with advanced computational resources, at a large scale it is prohibitively expensive to run finite element models which directly represent the material microstructure.The range of applications where this problem is encountered is diverse.In another example [31], a 3D representation of an austenitic steel microstructure is input into an image based FEM model.In this work, a lower resolution of the original image is used in order to reduce the computational cost of the FEM analysis.
A moving window (MW) homogenization approach may be thought of as an advanced alternative to a low resolution analysis.Windowing decreases the size of the stiffness matrix in the FEM model by using homogenized parameters to describe each element in the reduced FEM mesh; see Figure 1.This type of moving window analysis requires an underlying homogenization algorithm.Though homogenization may be a computationally expensive process, it is much less expensive to model smaller windows separately, and this process can easily be parallelized [38].The larger the window size, the lower the overall computational cost of modeling a large scale structure.The accuracy of models with large window sizes can be improved by tailoring the homogenization technique to maximize accuracy of the macroscale response.
Figure 1: A schematic showing the procedure used in moving window homogenization.The square shown in red is a "window."Each window is homogenized to determine its apparent stiffness where C is a 4 th order tensor mapped to row i and column j on a new grid (right).This grid is then used as a finite element mesh in a macro-scale analysis that is less computationally expensive than modeling the original binary microstructure (left).
Computational expense in the FEM analysis decreases with increasing window size.
When used as input into a finite element analysis of a large scale structure, meso-scale material property fields have the potential to make the structural scale problem tractable.However, as the window size increases, a decrease in accuracy of the prediction of macroscale behavior is often seen.Methods under development have the potential to advance both the accuracy and the efficiency of meso-scale homogenization by windowing.Homogenization techniques have often been studied from the perspective of improving the prediction of effective properties.The goal of this research is to improve the accuracy with which moving window homogenization predicts the macro-scale response of a structure under consideration.
There has been a great deal of research on homogenization of composite materials.The approaches are numerous; one useful broad distinction is made by Hollister and Kikuchi [24] between asymptotic expansion homogenization methods (see [5,19,32]) and standard mechanics analysis.In the latter method, a chosen representative volume element (RVE) is analyzed using either uniform traction or uniform displacement boundary conditions to relate average strain to average stress (to give one example, see [29]).More recently, methods that build upon the asymptotic expansion homogenization method include the Voronoi cell finite element method (VCFEM) [14][15][16], and multi-scale modeling approaches developed by Fish [12,13]; for a review of these and other representative works see [7,11].In both asymptotic homogenization and standard mechanics approaches, an RVE is typically used to decouple the analysis of the composite at the local and global levels to address the multi-scale problem.Another approach involves substructuring large scale problems without the use of an RVE [36,38].This strategy involves a posteriori adaptive local reaveraging to reduce error.The Generalized Method of Cells GMC [1] is a homogenization method requiring relatively low computational cost, employing the use of a periodically repeating unit cell (RUC), rather than an RVE (a comprehensive distinction between the two is made in [8]).Overall, common drawbacks of many of the homogenization methods currently in use are high computational cost and the assumption of the existence an appropriate RVE or RUC.
Often it is difficult to determine the size of an RVE [9], or even prove that one exists, especially in the case of a material exhibiting nonlinear behavior [17].For functionally graded materials, a representative volume element does not exist and non-RVE based numerical methods must be used [10].In certain cases, the scale of the RVE is prohibitively large to be determined by experimental testing or the geometry of a structure does not allow for the incorporation of a full scale RVE element [26].There are two main disadvantages of an RVE approach, which apply to all materials.First, when nonlinear behavior or damage initiation is considered, too much information is lost regarding the local variation in the microstructure to carry out an accurate analysis.Second, an RVE approach effectively eliminates randomness in the meso-scale material representation, which is needed to statistically quantify structural reliability.
As an alternative to the RVE approach, a windowing approach has been used with success in nonlinear analysis [2-4, 28, 30].Meso-scale random fields generated by windowing have been used as a basis for stochastic simulation [18].For a review of efforts in these directions see [6].However, fundamental questions must be addressed so that a windowing method combined with a relatively simple homogenization scheme may ultimately reach its full potential for widespread use in practice.As noted by Yuan and Fish [35], adoption of new computational homogenization methods by industry is in an "embryonic stage," in part due to the fact that new architecture is needed in commercial FEM codes to transfer data between scales.While windowing methods proposed in this work do not provide iterative data transference between scales, the fact that homogenization is performed a priori means that no modification to the FEM code itself need be performed.This is an advantage which promotes its potential use in a wide range of applications.

MOVING WINDOW HOMOGENIZATION
In this work, statistical volume elements (SVE), also called windows, are considered, which are below the scale of the representative volume.These windows are homogenized; the material property fields generated are then used as input into a finite element analysis in the manner shown in Figure 1.The first question that arises is how to generate so-called "apparent properties" on each meso-scale SVE.These are distinguished from effective properties on an RVE.
The following definition of the RVE was given by Hill [23]: "the apparent overall moduli" must be "effectively independent of the surface values of traction and displacement, so long as these values are macroscopically uniform."It has been proven by Huet [25] and further established in subsequent work (see [20,21]) that a test of the material under kinematically uniform boundary conditions (KUBC) provides an upper bound and a test under statically uniform boundary conditions (SUBC) provides a lower bound for the actual effective modulus of a material when the size of the material is smaller than a representative volume.As the size of the material increases, the two bounds converge.For a representative volume, these two values are identical, as stated by Hill.It has further been proven [20,21] that under mixed loading conditions (i.e. a combination of traction and displacement boundary conditions) the prediction of the effective modulus is bounded by the KUBC and SUBC apparent properties.Similarly, a periodic boundary condition has been shown to produce a prediction of the effective modulus that is also bounded by the KUBC and SUBC apparent properties [27].The order relationships described above lead to the so-called hierarchy of bounds [34,37].This hierarchy can be described by the following series of inequalities: Where is the effective stiffness, and are the Reuss and Voigt bounds, respectively, is the apparent compliance determined by an SUBC test of the material, is the apparent stiffness determined by a KUBC test of the material, is the size of the window of the material under consideration, and denotes averaging.It is noted that implies that for any tensor .The difference between and unity at a given length scale may be thought of as a measure of the departure of the material from its effective representation, [34]. ( Thus a fundamental problem when working with an SVE is that the inverse of the compliance tensor obtained in an SUBC test does not necessarily equal the stiffness tensor obtained in a KUBC test.A unique value for the constitutive tensor at the meso-scale is not readily defined.In the following analysis, the goal is to most accurately quantify the mesoscale constitutive tensor of each window in order to best approximate macro-scale behavior.

Beam Model
The macro-scale behavior of a simple beam model is evaluated where different mesoscale modeling techniques are used.Low resolution, binary representations of the two-phase material are compared against homogenized representations.Relatively large and small window sizes are used.Each representation of the original microstructure is used as input into a finite element analysis, and the results of these analyses are compared.The effect of discretization error in the models is separated from the effect of modeling error by comparing fully homogenized models with different levels of mesh refinement.Consider the beam in Figure 2 with dimensions and loading as shown.The inclusions are uniformly randomly distributed and do not overlap; the radius of each inclusion is 0.033m.Two contrast ratios are considered.In case 1, the inclusions have a Young's modulus 0 GPa and for the matrix GPa.In case 2, 0 GPa and GPa.Poisson's ration for both materials is taken to be zero.
A single microstructure is generated and discretized using a refined triangular mesh (Figure 2b).A pixelized mesh is also generated as shown in Figure 2c.The pixelized mesh is composed of 20x400 elements, each with dimensions of 0.1 x 0.1 m 2 .Both benchmark models are analyzed by finite elements to determine the deflected shape of the beam under load.
First, "fully homogenized" models are considered, meaning that a single homogenized constitutive tensor is calculated to describe the entire beam.This is performed by dividing the beam into 2x2 windows, homogenizing each window, then averaging the constitutive properties into one effective constitutive tensor.Also, 5x5 windows are used and the same averaging procedure is followed.The homogenization approach used in each case is a test using kinematic uniform boundary conditions.While the KUBC test is easily performed using a finite element approach, the SUBC test is not.As detailed in [37], the traction boundary conditions produce rigid body motion, and statically admissible fields are not generated by a standard FEM approach.Therefore, in the current models, the KUBC test is used at the meso-scale.Homogenization by KUBC over predicts the effective stiffness (Equation 1), and the fully homogenized 5x5 model is expected to provide a more accurate prediction than the 2x2 model according to the established hierarchy of bounds.As a way to evaluate the discretization error associated with a less refined mesh, both sets of fully homogenized effective properties are used in a "fine" mesh, where the beam is represented by a 20x400 element mesh as in the binary pixelized original representation.
A second class of models, which will be termed "moving window (MW) homogenized," is compared against the benchmark and fully homogenized models.In the MW homogenized models, the beam is again partitioned into 2x2 and 5x5 windows, and each window is homogenized using the KUBC mechanical analysis.In these models, the constitutive tensor of each window is directly input into a finite element analysis (see Figure 1).In the windowed models, a randomly varying meso-scale field is used to describe the constitutive behavior of the beam.
Finally, a set of "low resolution" models are generated.In these models, no homogenization is used; the representation of the two-phase material is binary.In the low resolution 2x2 model, every second pixel is sampled and the element assigned these properties; in the 5x5 model, sampling is performed on every fifth pixel.This is the simplest algorithm for reducing the computational cost of the macro-scale finite element analysis.Table 1 gives a summary and comparison of all the material models used in the analysis.

RESULTS AND DISCUSSION
Results of the beam models are shown in Figures 3-5 where beam deflection is plotted as a function of position along the length of the beam.In Figure 3, for both contrast ratio 2:1 (Figure 3a) and 10:1 (Figure 3b), the fully homogenized models are less accurate than the MW homogenized model when compared to the benchmark pixelized and refined mesh models.This difference is more pronounced with the high contrast ratio material (Figure 3b).The discretization error between a coarse and a fine mesh is also captured in these plots, and is shown to be small compared to the material modeling error.The difference between the pixelized and refined mesh models is comparable to the discretization error, and also smaller than the difference between the moving windowed and fully homogenized model results.Figure 4: Beam deflection as a function of position along the length of the beam for case 1 (material with 2:1 contrast ratio).For different window sizes, fully homogenized and MW homogenized models are compared against the benchmark pixelized model.For reference, legend numbers refer to model numbers listed in Table 1.
Figure 5: Beam deflection as a function of position along the length of the beam for case 1 (material with 2:1 contrast ratio).The plot shows a zoomed-in view of the region of maximum deflection for the most accurate models.For reference, legend numbers refer to model numbers listed in Table 1.
Figure 4 shows the fully homogenized material models and the MW homogenized models compared to the benchmark pixelized material representation for the 2:1 material contrast ratio case.Based on figures 3 and 4, moving window homogenized models show more accuracy than fully homogenized models.This is true for both contrast ratios (Figure 3) and for both the 5x5 and 2x2 cases (Figure 4).The MW homogenized 2x2 model is the most accurate, followed by the MW homogenized 5x5 model.Based on Figure 5, which isolates the most accurate material models for comparison, it is also observed that low resolution models perform better than windowed models, where mesh size is the same.In Figure 5, the 2x2 low resolution model is more accurate than the MW homogenized 2x2 model (note that both have the same mesh, so discretization error is not a factor).Similarly, the 5x5 low resolution model is more accurate than the 5x5 MW homogenized model.
Typically, when more computational effort is spent in characterizing a parameter, more accuracy is expected.However, the opposite is true in this case.The use of a static uniform boundary condition (SUBC) would provide a lower bound on the effective behavior of each window; the properties obtained by KUBC homogenization are known to be over-predictive of the effective stiffness of the material.But if the over-predictiveness of the KUBC homogenization were the main cause of the inaccuracy in the windowed models, the 5x5 MW model would be expected to be more accurate than the 2x2 MW model (see Equation1); this is not the case.
Finally, it is observed that the small window (2x2) MW homogenized model is more accurate than the large window (5x5) MW homogenized model, a result which cannot fully be explained by discretization error.In the fully homogenized models (Figure 4), the 5x5 model is more accurate than the 2x2, even in the case where the 5x5 coarse mesh model (4 x 80 elements) is compared to the 2x2 fine mesh model (20 x 400 elements).This suggests that the relatively large error in the MW 5x5 model compared to the MW 2x2 model in Figure 5 is largely due to modeling as opposed to discretization error.

CONCLUSION
A distinction is made between predicting effective stiffness and predicting macro-scale behavior.Accuracy in predicting effective stiffness refers to the nearness of the predicted apparent property to the effective stiffness bounded in Equation 1 by SUBC and KUBC apparent properties.It is shown in comparing the fully homogenized models based on ensemble averaging of 2x2 elements and of 5x5 elements (Figures 3 and 4) that a larger window size provides greater accuracy in this respect.Accuracy in predicting macro-scale behavior refers, in this analysis, to predicting the beam deflection.Often the implicit assumption is made that these two types of accuracy are the same.However, in the case of meso-scale moving window homogenization, it is shown that the larger 5x5 window size is less accurate.
It is seen that the models with more stress concentrations in the macro-scale FEM analysis are more accurate.Material property mismatch creates stress concentration; the higher the mismatch, the higher the stress concentrations that develop and the harder it is to model the effective composite material properties.A moving window homogenized model has greater material property mismatch than a fully homogenized model; the moving window homogenized models are shown to be more accurate.The low resolution (binary) models have greater material property mismatch than the MW homogenized models; the low resolution models are shown to be more accurate.Finally, the MW homogenized models with smaller windows have more stress concentrations which develop at the interfaces between unlike elements; the smaller window size models show greater accuracy than larger window size models.
These results suggest that accurately representing the material property mismatch (and associated stress concentrations) is of primary importance in modeling the macro-scale behavior of a composite.This is a different consideration than accurately representing effective material properties.This analysis motivates the development of an optimal homogenization scheme specific to moving window meso-scale analysis.
The models homogenized by moving window analysis are shown to be a dramatic improvement over fully homogenized models.It is realistic to expect that using a meso-scale homogenization technique optimized to predict macro-scale structural behavior, the moving window homogenized 5x5 results (shown in Figure 5) can be an improvement on the simple low resolution 5x5 prediction.In addition to the computational savings gained by using moving window analysis, meso-scale material property fields may be used as a basis for stochstic simulation.Using these fields, nominally identical random microstructures may be generated, which may then be used in a structural scale reliability analysis.To this end as well, it is important to develop optimal moving window homogenization methods.

Figure 2 :
Figure 2: (a) A schematic showing a composite beam under uniform load (not to scale).The fiber shape is shown using (b) a refined triangular mesh; and (c) a pixelized mesh.

Figure 3 :
Figure 3: Beam deflection as a function of position along the length of the beam for (a) case 1(material with 2:1 contrast ratio between phases); (b) case 2 (material with 10:1 contrast ratio between phases).For reference, legend numbers refer to model numbers listed in Table1.
Figure 3: Beam deflection as a function of position along the length of the beam for (a) case 1(material with 2:1 contrast ratio between phases); (b) case 2 (material with 10:1 contrast ratio between phases).For reference, legend numbers refer to model numbers listed in Table1.

Table 1 :
Material models used in the beam problem.