Comptes Rendus Mécanique

Experimentally validated combined sti ﬀ ness expression for


Introduction
The presence of porosities in solids punctuates many technical materials both natural and synthetic alike. In most practical cases, porosity presents itself in randomly shaped and distributed inclusions or voids. Quantifying the combined effect of such inclusions has been at the heart of mechanics in dealing with how to "homogenize" such materials. Much of the reported work on composites containing single and multiple inclusions of various types is based on the original work to determine the elastic field of ellipsoidal inclusion of Eshelby [1] and on the later work of Mori and Tanaka [2]. While the interior Eshelby tensor is constant with respect to the inclusion position, the interior Eshelby tensor is function of the matrix Poisson's ratio, shape, and orientation. The exterior Eshelby tensor is function of the matrix Poisson's ratio and the inclusion's shape, orientation, and position. Reported works (e.g., [3][4][5][6][7], and [8]) have found solutions for interior and exterior components of the Eshelby tensor for a single (mostly ellipsoidal) inclusion taking into account the effect of its volume fraction, orientation, and shape based on the work of Mori-Tanaka [2]. Such works do not account for the effects of the multiple inclusions' positional distribution within the medium nor for the number of these inclusions. Rather than consider only one inclusion's volume fraction or shape at a time as done in conventional homogenization analyses, this work accounts for the composite stiffness of a finite domain with multiple ellipsoidal inclusions.
This work is devoted to validating a new homogenization method for linear elastic composite materials by comparing the resulting estimate of the overall stiffness to estimates produced by finite element computations and also to strain measurements utilizing 3D-printed microstructures. This homogenization scheme of multi inclusions is an extension of recent work by the authors (Hage and Hamade, [9]) that calculates the stiffness of a domain containing a single inclusion with its corresponding shape, volume fraction, and orientation. The work is based in part on the superposition of strain field principle as described by Li et al. [10] that considers one inclusion in finite and infinite domains (without accounting for location) and also expands on the works of Ju and Sun [11][12][13]. The concept of superposition of strain field was used by Askari et al. [14] for modeling of inclusions and precipitation hardening in metal matrix composites. This work advances a combined formulation for the effective stiffness tensor of a finite domain while considering the combined effects of multiple ellipsoidal inclusion's attributes: volume fraction, VF, aspect ratio (accounted for via the Eshelby tensor expression), AR, location (position), r , and 3D orientation (via angles θ, ϕ, γ, where the inclusions are voids. The scheme relies on the average stress consistency condition introduced by (Li et al. [10]) that does not require a balanced stress whenever the Eshelby tensor becomes non-uniform on the ellipsoidal inclusion. This situation occurs as soon as a finite domain is considered in Eshelby's problem, and this situation is at the heart of the present study aiming at accounting for boundary effects. Indeed, in the case of an infinite matrix, the new estimate is perfectly identical to the one of Mori-Tanaka, which is also observed by comparing the new estimate to the one of Mori-Tanaka on configurations taken from the literature. Since the formulation of the stiffness tensor in Li et al. [10] takes into account only one inclusion and only one parameter (i.e., volume fraction), the was further verified here against similar cases to those of Mori-Tanaka's and taking only this one parameter. Also considered in this work the effect of multiple inclusions within the same representative volume element (RVE) by estimating the combined effect of each individual inclusion on the strain field for the composite as the summed effect of all inclusions (interior and exterior strains) within the RVE is calculated. An expression is developed based on a modified Mori-Tanaka method combined with the dual-eigenstrain (interior and exterior eigenstrain) method. From the dual eigenstrain method, an Eshelby tensor, dependent on the inclusion's shape and position, is deduced. Employing the strain field superposition principle, this is applied for multiple inclusions considering all the above variables into one formulation. A global coordinate system for the composite is adopted at its center, and the position of each inclusion is formulated and re-oriented using 3D transformation matrices with respect to this global system. The orientation of each inclusion is taken into account by multiplying the strain concentration tensors with a transformation tensor function of the 3D transformation matrices. The resulting formulation yields a compact expression so that effective composite stiffness may be determined for any number of inclusions with combinations of geometric (ellipsoidal) attributes of a composite containing multiple inclusions. These formulations are reported in Section 2.
A numerical validation section (Section 3) is included as a numerical validation of GSF stiffness solutions as compared with literature solutions. Thus, this work findings are applicable to combinations of variables that are not covered in other techniques based on Mori-Tanaka. Previous works are applicable in cases where one should fix one of the variables (position, shape, volume fraction, orientation, or number of inclusions). However, this work takes into account the combined effect of all these variables. This formulation's numerical estimates of stiffness are compared to simpler cases reported in the literature: cases with specific shapes and orientations of single inclusion in which one variable (e.g. inclusions position, shape, volume fraction, orientation, and number of inclusions) is fixed. Simple inclusion cases are considered that account only for one variable at a time and omitting the effect of location and number of inclusions, this work's estimates mimic those of the simplified Mori-Tanaka scheme. This accomplishment represents verification to literature-reported results using Mori-Tanaka scheme. The section contains one verified case of stiffness variation with inclusions volume fraction and another with shape (aspect ratio).
To its credit, the scheme accounts for the combined effect of inclusions' positions, shapes, volume fractions, orientations, and number. The formulation is also customizable for finite or infinite domains. In addition, it is able to handle multiple numbers of inclusions with specified locations, orientations, volume fractions and aspect ratios. This aspect is crucial in the field of composites where, not only simple cases, but randomly located cases of inclusions with various aspect ratios must be dealt with.

Homogenization formulation for multiple inclusions
In this section, the application of the strain field superposition principle is delineated. To account for multiple inclusions' 3D orientations and locations within the RVE, the effect of each individual inclusion is estimated w.r.t. the RVE Cartesian coordinate using translation and rotation transformation matrix.

Transformation from local to global coordinates
In order to account for the effect of multiple inclusions and their orientations, a representative volume element (RVE) with its corresponding Cartesian coordinate system is shown in Figure 1. Also shown in the 3D RVE are inclusions distributed with multiple orientations with their corresponding rotation angles identified with respect to the RVE Cartesian coordinate. To move from each inclusion's local coordinate system to the global RVE coordinate system, a transformation matrix is utilized to account for the orientation of each individual inclusion within the RVE (Hage and Hamade, [9]). The composite is considered as domain Ωu having multiple inclusions Ω i I from i = 1, . . . , N (where N is the total number of inclusions). Figure 1 is an illustration of such a finite domain with multi ellipsoid inclusions. Shown is the composite domain (and its global coordinate system) with mix of multiple inclusions and various attributes (with their local coordinate systems). The origin, Ou, of the coordinate system is set at the center of the domain. The ori- 3 )) where (1) refers to the considered inclusion (for the sake of non-confusion "(1)" is used for the considered inclusion and "i " is used for the studied point P). The coordinates (x (1) 1 , x (1) 2 , x (1) 3 ) for the origin of each inclusion are obtained with respect to the origin of the composite domain coordinate system. Using the transformation matrix and considering that r is the distance between point P in the domain and the origin of the considered inclusion, the coordinate position vector is reoriented to the inclusion coordinate system and reformulated in the composite domain (Hage and Hamade, [9]) refers to the coordiantes of any point P in the space In addition, Q i j is a transformation matrix is utilized to account for the orientation of the inclusion within the 3D RVE and is utilized to jump from the inclusion local coordinate system to the global RVE coordinate system [9]. This reformulated coordinate position distance r is used to determine the exterior Eshelby tensor and is function of the orientation and coordinates of the center of each inclusion w.r.

Effective stiffness tensor: generalized stiffness formulation
The effective stiffness tensor of a composite with multiple inclusions is calculated based on the superposition principle of strain fields (Li et al., [10]). The effect of the boundary conditions on the inclusion is taken into account by the position factor and the distance between the boundary and the inclusion, which is added as a variable. The authors examined two cases: one in which inclusions are located near the center of the RVE and another case where inclusions are closer to the boundary of the RVE. Although the code allows for ellipsoidal or spherical shaped RVEs, reported in this work is a finite cube-shaped RVE. Cluster of inclusions is inserted into the finite RVE. The Matlab® code checks whether inclusions satisfy the condition of being fully contained within RVE.
At the inner side of the RVE, the elastic field strain is divided into two fields: one from the background (external boundary loads) and another due to the presence of the inclusion field [10]. The new average strain in each phase is the sum of the background strain and the disturbance strain. The averaged interior and exterior Eshelby tensors (being key for determining the elastic field of an ellipsoidal inclusion interiorly and exteriorly) for each inclusion are summed.
In this work, each inclusion is treated as a subdomain within an infinite domain where a stressfree eigenstrain is prescribed in the inclusion and vanishes outside. Since multiple inclusions are considered inside the bounded RVE, the inclusion inside the cube is a finite domain and the exterior Eshelby is calculated between the inclusion and its boundary as can be shown in Figure 2 whereC is the reference solid stiffness tensor and C I i is the stiffness tensor of each inclusion, Ω E is the exterior domain of the inclusion used to calculate the strain field outside the inclusion domain, and the dashed circle corresponds to the shape of the inclusion domain which could be elliptical (or spherical, a special case).
The Eshelby tensor for one inclusion is calculated using the summation of the interior classical Eshelby tensor (adapting the formulation by Mura [15]) and exterior Eshelby function of the position x (using the formulation by Ju et al. [11]). These formulations related to the Eshelby tensor are based on two components: one is uniform but a function of the matrix Poisson's ratio, shape, and orientation and is constant with position and the other is a function of the matrix Poisson's ratio, shape, orientation, and the re-oriented coordinate vector position.
Translating the effective stiffness tensor into the global domain, the strain concentration tensors for each inclusion A E , and A I is calculated in the local domain and equivalent strain concentration tensors are calculated in the global system using a 4th rank transformation tenor, Q i j kl (Koay [16], this tensor is function of the transformation matrix Q i j [9] defined in the previous section). The "Bond transformation" [17] is a tensor that accounts for the inclusion's 3-dimensional attributes.
The derived expression for the effective stiffness tensor for multiple inclusions reported by Li et al. [10] and Hage and Hamade [9] becomes where interior and exterior Eshelby tensors calculated for each inclusion (Mura [15] and Ju and Sun [11]) Q i j kl rotation tensor function of Q i j as defined in [9].
In (3), both A E i and A I i as shown are function of averaged interior and exterior Eshelby tensor S(r ) (these tensors are function of ν m is the Poisson ratio of the matrix, α is the aspect ratio of the inclusion and the position of the inclusions r ) and the stiffness tensors fo the matrix, C E , and inclusions C I i and the reference solid stiffness tensorC i = aC I i + (1 − a)C E for each inclusion thus A E i and its complex relation with all the aforementioned variables could be simplified as a function named herein by "g " as follows , the effective stiffness tensor is function of all 5 considered inclusions attributes: volume fraction ( f ), aspect ratio (α evaluated in S(r ), orientation (evaluated in r and using Q i j kl ), locations (r ), and number of inclusions (averaging over multiple inclusion by using the principle of superposition of strain fields and summing for i = 1 to N ). As "a" is a parameter used to describe the comparison solid, it has to be chosen either to yieldC = C E orC = C I where the method would degenerate either to interior homogenization or the exterior homogenization, respectively. In this work, the parameter "a" is set to 0.01 thus favoring the exterior finite RVE. This value was chosen in order to take into account the effect of all the inclusions on each other inside the exterior cubic RVE without neglecting the effect of each inclusion with respect to its interior domain. It was so chosen to lie in between these two special cases (exterior and interior).
Using Matlab®, the inverse tensors are calculated. Based on the component's inverse matrices, an extension formulation of Li et al. [10] is the aggregate stiffness formulation becomes: Based on strain field superposition (summation) of each inclusion, Equation (4) incorporates the aggregate effects of multi inclusions including location/distribution within the medium. In addition to the 3 classical variables of volume fraction, shape, and orientation, this formulation accounts for the additional effects of 2 additional variables: number count and locations of inclusions.
Most homogenization literature report on calculations of stiffness as function of inclusions' volume fraction and shape considering either uniform orientation or random orientations. To the best of the author's knowledge, no other reported work takes into account the combined effects of all 5 parameters (volume fraction, shape, orientation, location, and number) at the same time in one formulation with a range of variation for each parameter. This GSF work accomplishes this by extending the original formulations by Li et al. [10] to include the effects of inclusions' (1) Orientations: orientations were never taken as a parameter with a level range as in this work (2) Locations: while formulating the exterior point Eshelby's tensor of an ellipsoidal inclusion. To account for this parameter, works in [11] and [15] were combined. Considering the concept of superposition of strain fields (interior and exterior Eshelby and the interaction with multiple inclusions) which has never been reported before (3) Inclusions' number by taking the interaction of the strain field for multiple inclusions. Noting that the superposition assumptions was based on calculating the interior and exterior Eshelby tensors/eigenstrain method where the exterior tensors are function of the position of the inclusions in the domain and there is no restriction regarding the positions of the inclusions in calculating the exterior tensors, the tensor could be calculated for any position, thus it was assumed that the superposition of strain field or in other words overlapping of strain fields is applicable at any position.
Comparing (3) and (4) with the work done in [10], one can see the added terms of (1) the summation of strain fields due to the effect of multiple (N ) inclusions in one matrix thus the effect of location is taken into account in this term inside the new Eshelby tensor dependent on r the inclusions' locations, (2) the effect of the inclusion orientation in the Q tensor. One contribution of this work is to add in one formulation (dubbed GSF) the effect of the inclusions' number, location, orientation, volume fraction and shape, which is a novelty to the field of composites homogenization. This GSF enhanced formulation, takes into consideration anisotropic composites, due to the added terms of inclusions' multiple shapes and orientations which makes the composite anisotropic in nature, and thus the authors will verify this formulation later on in this paper to 3D printed ABS plastics which are not isotropic, and prove the validity of GSF to anisotropic composites.

Solutions
In this section, the numerical solution of (3) and (4) are validated for 2 cases reported in the literature. The homogenization scheme in this work is a generalized method customizable to finite and infinite domains and may be applied to particular cases with same assumptions as in Mori-Tanaka. The comparisons below validate the notion that (3) and (4) yield comparable results to those of Mori-Tanaka in which one of the 5 variables (inclusion position, shape, volume fraction, orientation, and number) is fixed. While classical homogenization methods can be applied for cases where the number of inclusions is larger than 10, GSF could be used for any number of inclusions.

Volume fraction validation
This case compares GSF calculations of normalized Young's modulus values versus volume fraction of porosity values obtained from the homogenization work reported by Dai et al. [18] for 3-phase composite containing a centered spherical inclusion (Zener anisotropic ratio = 1). Each constituent is elastic and isotropic.
The spatial locations of the inclusions are arranged in such way that the composite is isotropically homogeneous (thus E 11 = E 22 = E 33 = E ) with material properties:

Shape validation
Rouhi and Rais-Rohani [19] reported on the general framework of modeling the stiffness properties of carbon nano-fibers enhanced composite materials. A nano-enhanced matrix with CNF as the reinforcement material and vinyl ester resin as the matrix having the following elastic properties: E 0 = 3.5 GPa, E 1 = 450 GPa, and ν 0 = ν 1 = 0.3. Figure 4 represents the axial Young's modulus with respect to the elliptical aspect ratio of the inclusion obtained using the developed homogenization code in comparison with the published results.

Validation of the generalized stiffness formulation: experimental compression testing and fem of domain cubes containing multiple inclusions
In Mori-Tanaka, one inclusion variable (e.g., inclusion position, shape, volume fraction, orientation and number) has to be fixed. The scheme in this work is presented as capable of addressing  any number and combination of inclusion states. In this Section, the scheme will be validated for finite domains where novelty is claimed at the expense of relying on the "average stress consistency condition" of Li et al. [10]. Corroboration is performed by comparing against compression experiments on 3D printed composite specimens containing multiple inclusions and by the finite element method (FEM). This section reports on the design, details, and results of the performed experimental compression tests.

Design of experiments of finite domain with multiple inclusions
Several composites are designed using design of experiment (DOE) SAS/JMP® software to determine the effect on the effective Young's modulus of the composite of inclusions' volume fractions, shapes, orientations, numbers, and locations within the matrix. Full factorial DOE is utilized with 5 considered variables and their value ranges as listed in Table 1. "Agglomerated" locations refer to inclusions are positioned near the center of the cube while and "dispersed" represents inclusions that are far away from the center and nearby the walls of the cube. These positions are chosen to cover 2 extreme cases of inclusions locations within the cubical RVE to reveal whether there is an effect of locations.
In full-factorial experiments, the number of experiments required is 13 * 6 * 14 * 7 * 2 = 15288 cases. Excluding cases where inclusions protrude outside the composite domain of the cube or overlap with other inclusions reduces this number. The cases that satisfy this condition are listed in Table 2. These cases are tested for two schemes for inclusion distribution: 2520 agglomerated inclusions case and 2520 evenly dispersed inclusions cases for a total number of 5040 cases.

3D printing of domain cubes containing multiple inclusions
In order to generate test specimens with predefined inclusion's volume fraction, aspect ratio, orientation, number, and location, 3D printing technology has been adopted using a secondgeneration cube 3D printer (Cubify, Cube 3D Printer 2nd Generation: printer model number: 381000) and using ABS plastic filament (product number CC3D2 GEN-ABS). A subset of the cases listed in Table 2 are 3D-printed, experimentally tested, and modeled using FEM. Taguchi's orthogonal arrays full factorial is utilized with 5 variables: • Shape variable with 2 levels: spheres (AR = 1) and elliptical inclusions (AR = 1.25).
• Position variable with 2 levels: dispersed and agglomerated inclusions. The number of experiments was reduced due to the deletion of repeated or redundant cases. Additional reductions were required due to practical (printing) reasons where the aspect ratios of interest were limited to a range of 1-1.25. Therefore, the scheme is currently limited to spherical and elliptical inclusions but remains untested in the case of flat and elongated inclusions. Consequently, the number of test cases was reduced to the 36 cases listed in Table 3.
The effect of number of inclusions is taken into account via 3 experiments all at VF = 0.1: (1) 2 spheres and 6 ellipses with dispersed positions, (2) 4 spheres and 4 ellipses with dispersed positions, and (3) 6 spheres and 2 ellipses with dispersed positions, with ellipses at 0 • orientation. The effect of inclusions shape (calculated as % circular or % elliptical inclusions) is also considered.
For all cases, the finite domain is taken as uniform with cubic dimensions of 16 * 16 * 16 mm 3 . The composite has a matrix composed of ABS plastic shaped by fusing filament wire in a 3D printing apparatus. Using Cubify Design® software, the cases listed above are designed to generate *.stl files with dimensional tolerance of 0.001 mm. The dimensions of the spheres and ellipses inclusions are designed to have the predefined volume fractions and aspect ratios for all 36 cases.

Compression testing of cubes with multiple inclusions
Compression testing was performed on the cases listed in Table 3 in accordance with ASTM D695-10 Standard Test Method for Compressive Properties of Rigid Plastics. This test method covers the determination of the mechanical properties of unreinforced and reinforced rigid plastics when loaded in compression at relatively low uniform rates of straining or loading and when test specimens of standard shape are employed. Utilized was a Tinius Olsen Universal testing machine (UTM) using load cell of 10 KN with speed of testing of 1 mm/min. Compressive modulus of elasticity values were arrived at from the linear portion of the stress-strain curves. Figure 6 shows such load-stroke (left) and stress-strain (right) curves. This was done by selecting a point along the straight-line portions and dividing the compressive stress by the corresponding strain value. Young's modulus values are obtained experimentally from results reported for all 36 cases and were calculated for all points after 1000 N at the linear part of the curve since before 1000 N recordings correspond to aligning the upper plate grip with the horizontal surface of the cubes this is why those points were disregarded and after 1000 N the alignment of both surfaces was perfect and this is where the modulus of elasticity was calculated. Each experiment was repeated 3 times by printing and testing these 3 different case cubes. Test results obtained from compression tests were almost identical for the same test case and the values reported below are the corresponding average values.

FEM of domain cubes containing multiple inclusions
This section reports on the details and results of the FEM simulations that are run along with the experimental cases tested in the previous Section. Discussed first is finite element simulations description followed by results of the FEM.

FEM model
All cases listed in Table 3 are modeled with FEM. The same *.stl files generated using Cubify Design® software used for 3D printing are used as computer-aided design (CAD) input files imported into the FEM software. The mechanical properties assigned to the defined cube are elastic and isotropic. For the matrix, the strength was tested in compression mode for multiple solid cubes (no inclusions) along the principal directions 1, 2, and 3. The results varied within a range of ±5% which suggests that the matrix is isotropic. An average of the measured modulus values in the three directions has been considered and the ABS printed compressive modulus was found to be 0.398 GPa and its Poisson's ratio of 0.35. This modulus value is significantly lower than reported moduli of solid ABS (typically 2.2 GPa) or that reported for other printed ABS (of 1.5 GPa in Weiss, [20]). No apparent reason for this obvious disparity but given the wide variation in printer vendors and materials, it is possible that the proprietary filament was an ABS blend with lower stiffness or that the fusion process has negatively affected the final mechanical properties. Inclusions are set as voids thus E i = 0 GPa and ν i = 0.

Meshing, simulation control, and boundary conditions
The commercially available DEFORM® 3D version 11 software package (from SFTC) was used. Tetrahedral elements are used to mesh the cube. Finer mesh is set at the top surface resulting in 130,000 elements. Active re-meshing is considered with relative interference ratio of 0.7. The time step is set at 0.001 s/step based on the smallest element size and the displacement rate as recommended by the software guidelines. Boundary conditions (BCs) are set so that −z(−1) surface of the cube has displacement uz = 0; +z(+1) surface has velocity displacement BC of 1.3 mm/min to mimic the upper rigid plate motion during uniaxial compression test. A heat transfer condition for all surfaces with the environment is set.

FEM results
Load versus stroke curves are used in FEM to determine the elastic Young's modulus for all cases considered. Figure

Results and discussion
The effective stiffness tensors are determined numerically for all cases using the developed homogenization methodology. The Young's modulus values normalized with respect to that of the matrix (E 11 /E m ) from (4) are plotted versus volume fraction (VF) and location, average aspect ratio (AR), orientation, and number of inclusions. Overlaid when available, are experimentally and FEM-obtained stiffness values for those cases. homogenization work in this work, VUL: Voigt upper limit (Voigt [21]), and HSLL: Hashin-Shtrikman lower limit (Hashin and Shtrikman [22]). Figure 8 contains data points and fit lines corresponding to these 4 cases:

Normalized longitudinal Young's modulus versus volume fraction (and location/ agglomeration)
• 8 spheres dispersed • 8 spheres agglomerated • 8 ellipses dispersed at 0 • • 8 ellipses agglomerated at 0 • Results suggest that the stiffness calculated values from this homogenization methodology yield good agreement between numerical, experimental, and FEM results. The effective stiffness values of the composite depend strongly on the inclusions' volume fraction. Results of the normalized Young's modulus shows that E 11 /E m decreased with increased volume fraction for all cases considered. These results are in congruence with results reported by Xu et al. [23] where the predicted elastic moduli of particulate cementitious composite decreased with the fraction of voids of microencapsulated phase change materials (MPCMs) and increased with the fraction of quartz grains (hard inclusions). As the fraction of hard inclusions increase, the effective moduli would increase. Chen et al. [24] found that the elasticity of shale rock decreases with increasing volume fraction of multi-inclusions (as voids have a null modulus compared to the matrix similar to the case studied in this work) that is similar to the results found herein. However, their results are the inverse of those found in Xu et al. [25] since the inclusions are now harder than the matrix. The normalized effective Young's moduli of polymer nanocomposites increased with hard interfaces having modulus of elasticity larger than that of the matrix by 12 times.  Results also suggest that clustering of inclusions versus them being dispersed have little effect on the axial Young's modulus E 11 /E m . The difference between stiffness of spherical or elliptical inclusions whether dispersed or agglomerated is relatively insignificant and is estimated at about 2% only. that higher AR values imply rounder or more spherical inclusions). It could also be inferred that normalized modulus exhibits slightly decreasing trend the more the inclusions are circular by 1.676% (at VF = 0.1) and 3.3084% (at VF = 0.2).

Normalized longitudinal Young's modulus versus orientation (and location/ agglomeration)
As orientation is redundant in the case of spherical inclusions, the orientation of other shapes (here elliptical) is investigated as for their effect on the stiffness of the composite. Figure 10 shows results for the normalized axial Young's modulus E 11 /E m versus orientation for ellipses having AR = 1.25 at orientations of 0 • , 90 • and random, RAND VF = 0.1 and VF = 0.2. At such a value of aspect ratio, the effect of orientation is noticeable. Plotted, for all orientations, are results for experimental results, EXP, finite element FEM and numerical solution from (4), NUM. For the experimental results, shown in Figure 10 are error bars with respect to the arithmetic mean of the tests done for each case using the t -student distribution at a 95% confidence interval. Analytical, finite element method, and experimental stiffness results show the difference between elliptical inclusions oriented at 0 • with respect to the axial axis x1 and inclusions oriented at 90 • with respect to the axial axis x1 is relatively significant with about 4.5% max relative difference recorded at AR = 1.25 (VF = 0.1) and 9.022% (VF = 0.2). Very good matches between experimental and numerical/FEM results are recorded. Observed variations (error bars) in the experimental results may be related to 3D printing of the tested cubes, UTM results, or human factors may affect the accuracy of experimental data. For the case of 8 ellipses dispersed, normalized Young's modulus values decreased by 2.17% in going from 8 ellipses oriented at 0 • to random orientations (VF = 0.1), 4.14% (VF = 0.2), 1.53% for the case of 8 agglomerated ellipses (VF = 0.1) and by 6.96% (VF = 0.2). It can be inferred that the difference in orientations is more apparent at higher volume fractions. . For larger number of pores, this effect dwindles significantly. This is in agreement with Ju and Sun [12] and with Xu et al. [26] when the number of particles is greater than 49, the average value of elastic moduli of particle-reinforced composites PRCs particles tend to be stable) containing nonspherical. As the number of inclusions increases, their size would decrease for the same volume fraction. To ensure that all of the inclusions are bound within the RVE, investigated only is the case of clustered inclusions (agglomerated) near the RVE center.

Normalized Young's modulus versus dispersion/agglomeration
Some of the data reported in Figure 8 is re-examined here with focus on agglomeration and plotted as bar charts. For the below listed cases, Figure 12 (VF = 0 and VF = 1 not shown) shows the normalized axial Young's modulus E 11 /E m numerical values found from (4) for both dispersed and agglomerated cases: It can be inferred that location has little effect on the modulus where maximum percentage difference recorded is about 2% for the range of volume fractions examined. This may be due to the small size of the RVE considered where the location of the inclusions does not appear to significantly affect the elastic properties of the composite.

Significance of the findings
For all 36 cases tested experimentally, the normalized stiffness values determined from the numerical solution of (4) are plotted in Figure 13 against the experimentally measured normalized values. These Young's modulus values are determined at identical states of volume fractions, orientations, aspect ratios, locations, and numbers of inclusions. The figure shows a cloud of points, segregated by volume fraction values of 0.1, 0.15, 0.2, and 0.3, located close to the y = x line which shows that experimental and numerical results are quite comparable. This is indicative of the quality of the formulations advanced in this work.

Conclusions
The generalized stiffness formulation (GSF) scheme is a combined expression (4) for the effective stiffness tensor of a multiple-inclusion composite. Inclusions states have five geometric and other attributes of volume fraction (size), shape, orientation, location, and number within the solid matrix. The formulation is first validated by comparing its numerical solutions with two literature-reported stiffness values of relatively simple configurations. With a DOE-developed test matrix accounting for wide-ranging values of these attributes, this formulation is further validated for a wide range of multiple inclusions composites using experimental compression tests and FEM simulations.
Examining the values of normalized axial modulus E 11 /E m versus the 5 inclusions' attributes, the following conclusions are drawn: • Volume fraction of inclusions in matrix (porosity, %) is very significant. E 11 /E m values decrease (increase) significantly for increasing (decreasing) porosity. • Average aspect ratio, AR (circularity) of the inclusions. Trends of E 11 /E m with % circular inclusions (average aspect ratio) are inconclusive to make a definitive determination on their effect on composite stiffness. A slight decreasing trend of the modulus is shown with increasing number of spherical inclusions. • Orientation of the inclusions is fairly significant in its effect on E 11 /E m . For the case of 8 ellipses dispersed, normalized Young's modulus values decreased by 2.17% in going from 8 ellipses oriented at 0 • to random orientations and by 1.53% for the case of 8 agglomerated ellipses at VF = 0.1, 4.14% at VF = 0.2, 1.53% for the case of 8 agglomerated ellipses at VF = 0.1, and 6.96% at VF = 0.2. It can be inferred as well that the difference in orientations is more pronounce at higher volume fractions. • Number of inclusions is significant (up to a point). The normalized axial effective Young's modulus E 11 /E m increases with the number of circular inclusions up to N = about 10 after which the effect is minimal. For the same volume fraction, as the number of inclusions increases their size would decrease. Therefore, in order to ensure that all of the inclusions are bound within the RVE, the case of the inclusions clustered near the center of the RVE is investigated. • Location of the inclusions in the matrix (agglomerated versus dispersed) appears to have little effect (4.306%) on the normalized axial modulus E 11 /E m .