View metadata, citation and similar papers at core.ac.

uk brought to you by CORE

provided by Elsevier - Publisher Connector

International Journal of Solids and Structures 51 (2014) 2687–2695

Contents lists available at ScienceDirect

International Journal of Solids and Structures

journal homepage: www.elsevier.com/locate/ijsolstr

Prediction of the effective elastic moduli of materials

with irregularly-shaped pores based on the pore projected areas
Borys Drach a,⇑, Andrew Drach b, Igor Tsukrov c
a
Department of Mechanical & Aerospace Engineering, New Mexico State University, Las Cruces, NM, USA
b
Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, TX, USA
c
Mechanical Engineering Department, University of New Hampshire, Durham, NH, USA

a r t i c l e i n f o a b s t r a c t

Article history: Statistical modeling is used to correlate geometric parameters of pores with their contributions to the
Received 7 October 2013 overall Young’s moduli of linearly elastic solids. The statistical model is based on individual pore
Received in revised form 21 March 2014 contribution parameters evaluated by finite element simulations for a small pore subset selected using
Available online 12 April 2014
the design of experiments approach, so there is no need to solve the elasticity problem for all pores in
the material. A polynomial relating pore geometric parameters to the contribution parameters is then
Keywords: fitted to the results of the simulations. We found a good correlation between normalized projected areas
Stochastic approach
of the pores on three coordinate planes and their contributions to the corresponding effective Young’s
Effective elastic moduli
Irregular pores
moduli. The model is applied and validated for two large sets of pore geometries obtained by X-ray
Micromechanical modeling microcomputed tomography of a carbon/carbon and a 3D woven carbon/epoxy composite specimens.
Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction superspherical shapes (Onaka, 2001; Hashemi et al., 2009;

Sevostianov and Giraud, 2012) To evaluate contributions of pores
Effective stiffness of linear elastic solids is reduced by the pres- having irregular shapes, when no analytical solution for a given
ence of pores. Historically, this reduction has been quantified in shape is available, various numerical techniques, such as finite
terms of the relative volume fraction of pores (porosity) assuming element or boundary element analyses, can be used.
all pores to be spherical, see Dewey (1947), Mackenzie (1950). Direct micromechanical modeling is the preferred way to pre-
Later, the pore shapes were included in consideration by utilizing dict effective elastic properties of materials with pores (with other
the famous Eshelby solution (Eshelby, 1957, 1959) for ellipsoidal possibilities being variational bounds on a property as in Hashin
inclusions and modeling non-spherical pores as ellipsoids of and Shtrikman, 1962 and later publications, or finite element sim-
various eccentricities, see Wu (1966) and discussions in ulations for the entire representative volume element with many
Kachanov et al. (1994) and David and Zimmerman (2011). For pores as in Roberts and Garboczi (2000), Arns et al. (2002) and
two-dimensional pores, in addition to circular and elliptical González et al. (2004)). It involves evaluation of contribution of a
shapes, the analytical solutions for various quazi-polygonal holes certain pore type to the effective elastic properties by solving the
have been utilized, see Zimmerman (1986), Tsukrov and elasticity problem for a pore of that type, and then incorporating
Kachanov (1993), Jasiuk et al. (1994), Kachanov et al. (1994), the solution into a certain micromechanical method. For irregular
Ekneligoda and Zimmerman (2006, 2008b), Zou et al. (2010). In pore shapes, the elasticity problem is usually solved numerically.
three-dimensional case, only a limited number of analytical However, actual microstructures may contain large numbers of
solutions for non-ellipsoidal inhomogeneities are available, diverse pore shape types, and solving elasticity problems for all
including cuboids (Lee and Johnson, 1978), polyhedra (Rodin, 1996; of them with reasonable accuracy may become computationally
Nozaki and Taya, 1997; Lubarda and Markenscoff, 1998), and challenging and practically impossible.
It is therefore desirable to be able to evaluate how irregularly
shaped pores influence the effective material properties without
⇑ Corresponding author. Address: Mechanical & Aerospace Engineering, P.O. having to solve the elasticity problem. For this purpose, researchers
Box 30001/MSC 3450, New Mexico State University, Las Cruces, NM 88003-8001, have been trying to identify the appropriate geometric parameters
USA. Tel.: +1 (575) 646 8041; fax: +1 (575) 646 6111.
(in addition to the overall porosity of the material) to be used for
E-mail addresses: borys@nmsu.edu (B. Drach), andrew.drach@utexas.edu
the predictions of effective properties.
(A. Drach), igor.tsukrov@unh.edu (I. Tsukrov).

http://dx.doi.org/10.1016/j.ijsolstr.2014.03.042
0020-7683/Ó 2014 Elsevier Ltd. All rights reserved.
2688 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695

Zimmerman (1991), for example, showed that compressibility use the design of experiments (DoE) approach (Ryan, 2007;
of 2D pores correlates well with the undimensionalized ratio of Myers et al., 2009) to identify the list of input parameters (in our
pore perimeter P to its area A defined as P2 =A (a three-dimen- case, pore shape geometric parameters) that corresponds to the
sional analog of this parameter, the surface-to-volume ratio, was specified level of model prediction variance across the whole
used in the stochastic model of Drach et al., 2013). Zimmerman’s parameter space. We design the ‘‘experiment’’ (the combination
observation was later supported by numerical studies on several of irregular shapes for FE simulations) using an I-optimal design
irregular 2D shapes (Tsukrov and Novak, 2002). However, such a module in JMP software (SAS Institute Inc, 2010). The optimization
parameter cannot be considered universal. For example, substitu- criterion for I-optimal designs is minimization of the integrated
tion of smooth pore boundaries by the very jagged ones (as illus- variance of the model predictions over the entire design space. This
trated in Fig. 1) will significantly increase the surface-to-volume means that the expected width of confidence intervals for model
ratio without making big changes in the overall response. This predictions will be approximately the same across the full param-
immediately follows from the comparison (or ‘‘auxiliary’’) theorem eter space.
of Hill (1965) as interpreted by Kachanov and Sevostianov (2005). We apply this method to two sets of pore shapes obtained by
According to their interpretation, contributions of two pores with microcomputed tomography of two different material samples:
very different surface-to-volume ratios shown by solid lines in carbon/carbon and 3D woven carbon/epoxy composites. The exam-
Fig. 1 are bounded by the same contributions of pores shown with ples of pore shapes are given in Fig. 2. For both sets we compare
dashed lines. The influence of corrugated boundaries on the overall contributions of the pores to effective elastic properties deter-
compressibility of pores is also discussed in Ekneligoda and mined by the direct FEA simulations for each pore with the ones
Zimmerman (2008a). based on its projected areas. Note that we use these sets of pore
Another possible set of geometric parameters to characterize shapes only to check the hypothesis that statistical estimates based
contribution of a pore to the effective elastic properties is its prin- on the projected areas can be used instead of elasticity solutions.
cipal moments of inertia (PMIs). PMIs are defined as the eigen- Our predictions cannot be directly used for carbon/carbon and
values of the matrix of moments of inertia given by 3D woven carbon/epoxy composites, because in our analysis we
Iij ¼ V qðrÞðr 2 dij  xi xj ÞdV, where qðrÞ is the material density of
R
assume the pores to be placed in an isotropic homogeneous mate-
the body for which PMIs are calculated (pores can be treated as rial which is obviously not the case for these composites. However,
homogeneous bodies with qðrÞ ¼ q ¼ 1), r is the distance to the we feel that the considered pore shape sets provide good data to
axis around which the moment is calculated, dij is the Kronecker test the hypothesis, as the pores in these two cases result from very
delta, xi are the coordinates x1 ; x2 ; x3 . PMIs have been previously different mechanical processes (porosity from incomplete filling by
used to approximate pores by the equivalent ellipsoids (Li et al., chemical vapor deposition of pyrolytic carbon vs. damage due to
1999; Drach et al., 2013). In Drach et al. (2013) we utilized statis- chemical and thermal shrinkage of epoxy) and thus have different
tical approaches to correlate pore principal moments of inertia to morphologies and shape distributions.
their contribution to the effective elastic moduli.
In this paper we investigate the hypothesis that contribution 2. Micromechanical modeling approach utilizing statistical
of irregularly-shaped pores to the effective Young’s moduli of analysis
porous material can be evaluated based on their projected areas
(‘‘shadows’’) by statistical analysis of effective compliance Characterization of individual pore contributions to the effec-
parameters. We begin by selecting a certain subset of pores from tive material response is based on the pore compliance contribu-
the full dataset. Then, we perform finite element (FE) simulations tion tensor – H -tensor proposed in Kachanov et al. (1994). Note
for the pores in the subset to quantify their individual contribu- that a similar tensor in the context of a microcrack was earlier
tions to the effective compliance of the material. This data is given by Horii and Nemat-Nasser (1983). The fourth rank tensor
then used in regression analysis to construct an approximating H is defined as a set of proportionality coefficients between remo-
polynomial model for the dependence of pore’s contribution to tely applied homogeneous stress field r0 and the additional strain
the effective compliance on its projected area. To verify that the De generated in the material due to the presence of a cavity:
model works for the full dataset, the constructed model is then
checked against the FE simulations on the new subset of randomly De ¼ H : r 0 ð1Þ
chosen irregular shapes not used in the model fitting. where ‘‘:’’ denotes the contraction over two indices. For a material
Because the accuracy of the model predictions (width of the with a large number of pores, a proper representative volume ele-
error bars) depends on the number of direct simulations used in ment (RVE) (Hill, 1963; Nemat-Nasser and Hori, 1999) can be
the data fitting, it is important to choose an ‘‘optimal’’ subset selected, and the effective compliance tensor of the material is
which corresponds to the desired model prediction variance. We given by

(a) (b)
Fig. 1. Two pores having similar contributions to the overall elastic properties but significantly different surface-to-volume ratios. Dashed lines show upper and lower limits
of their contributions.
 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695 2689

Fig. 2. Examples of pores extracted by microcomputed tomography from (a) carbon/carbon composite sample; (b) 3D woven carbon/epoxy composite sample.

S ¼ S0 þ HRVE example Drach et al., 2011), this method may not be practical
when contributions of tens of thousands of different pores are to
where S0 is the compliance tensor of the matrix material and HRVE is
be estimated.
the contribution from all pores present in the RVE.
Our method involves statistical analysis to derive a formula pre-
For dilute distribution of pores, the non-interaction approxima-
dicting contribution of a pore to effective properties based on its
tion can be used, and HRVE is found by direct summation of contri-
geometric parameters. First, we run a relatively small number of
butions from all individual pores in the RVE:
FEA simulations (100–200) for individual pores in a reference
HNI volume. The set of pores is chosen as an ‘‘optimal’’ subset of data
X
RVE ¼ HðiÞ ð2Þ
i using the DoE approach so that the pore shape parameters allow
to obtain the model with a specified level of variance across the
where HðiÞ is the compliance contribution tensor of the ith pore. For whole parameter space. Then a polynomial is fitted to the simula-
higher porosities when the non-interaction approximation is no tion results and estimates for the rest of the pores are inferred from
longer applicable, more advanced micromechanical schemes should the resulting empirical formula. Thus, the statistical analysis
be used. However, most of the first order micromechanical schemes results in a polynomial approximation model, which correlates
can be readily expressed in terms of the non-interaction approxi- pore contributions with their geometric parameters. If m geometric
mation, see Eroshkin and Tsukrov (2005). For example, predictions parameters are chosen, and the degree of approximating polyno-
for the overall pore compliance contribution tensor by the Mori– mial is n, the polynomial has the following general form
Tanaka method (Mori and Tanaka, 1973; Benveniste, 1987) in terms
of HNI
RVE is given by a simple formula Kachanov et al. (1994): PCC ¼ a0 þ a1 F 1 þ a2 F 2 þ    þ am F m þ amþ1 F 21 þ amþ2 F 22 þ   

HMT NI þ anm F nm þ anmþ1 F 1    F n þ    þ aN F mnþ1  F m ð6Þ

RVE ¼ HRVE =ð1  pÞ ð3Þ
where PCC is the predicted contribution to Young’s modulus (E ~1 ; E
~2
where p is the volume fraction of pores. The formulae for the effec-
or E~3 ); F 1 . . . F m are the considered pore geometry parameters
tive Young’s moduli in the next paragraph are given in the assump-
tion of the non-interaction approximation but can be converted to (model factors); a0 . . . aN are the coefficients determined by the
other micromechanical schemes if needed. standard least squares fit of FEA simulation data. In the case of a
To relate the pore compliance contribution tensor to the effec- 3-factor 2nd degree model, the expression is
tive Young’s moduli of the porous material we introduce the
~1 , E
~2 , and E
~3 defined as PCC ¼ a0 þ a1 F 1 þ a2 F 2 þ a3 F 3 þ a4 F 21 þ a5 F 22 þ a6 F 23
dimensionless parameters E
þ a7 F 1 F 2 þ a8 F 2 F 3 þ a9 F 3 F 1 ð7Þ
Ei 1
¼ ð4Þ
~i
E0 1 þ pE Such a statistical model can be used to predict contributions of
any pore with a specific set of geometric parameters.
where E0 is the matrix Young’s modulus (note that in all our consid- The examples of the model are given in Drach et al. (2013),
erations we assume the matrix to be isotropic), and Ei are the effec- where fitting of the FEA predictions is performed using the
tive Young’s moduli in three orthogonal directions. Three response surface methodology (Myers et al., 2009) based on the
parameters E ~i (i = 1, 2, 3) characterize the change in Young’s modu- following geometric parameters (model factors): principal
lus in xi -direction induced by the pores of the same shape in the moments of inertia of each pore and its surface-to-volume ratio.
case when these pores are parallel to each other. These parameters In that work, the set of irregular shapes used to construct the
can be expressed in terms of the H -tensor components for the cor- model was obtained from a microtomography scan of a carbon/car-
responding pore shapes: bon composite sample. The models were constructed for pore con-
tributions to the overall bulk, shear and Young’s moduli based on
~i ¼ E0 Hiiii
E ðno summation over repeating indicesÞ ð5Þ 150 pores selected using the DoE approach. The models were val-
p
idated on another set of 150 randomly chosen pores not included
In most of the previous publications, the pore compliance con- in the model. A reasonably good fit was observed for bulk and
tribution is found deterministically, either analytically for a pore in shear moduli: bulk modulus contribution R2 ¼ 0:98, RMSE ¼ 6:8%
unbounded medium (Eshelby, 1957, 1959; Kachanov et al., 1994) (hereafter the mean error, ME and analogously the root-mean-
or numerically, for example, by finite element analysis (FEA) per- square error, RMSE, values are presented in percentages, e.g.:
formed on a reference volume containing an individual pore, see ME ¼ 100%  MeanError=MeanResponse); shear modulus contribu-
Tsukrov and Novak (2002), Drach et al. (2011). Numerical simula- tion R2 ¼ 0:92, RMSE ¼ 3:6%. However, the regression’s fit values
tions appear to be the only option for the irregularly-shaped pores for Young’s moduli contribution parameters were found to be
for which no analytical solutions are available. Even though the lower (R2 in the range of 0:75 . . . 0:87). Thus, it appears that the
entire numerical procedure can be completely automated (see for principal moments of inertia with surface-to-volume ratio are
2690 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695

not the most suitable combination of geometric parameters to

evaluate the contribution of pores to the effective Young’s moduli
of porous solids.

3. Shape projected areas as the geometric parameters

determining pore contributions to effective Young’s moduli

We consider the projected areas of pores (‘‘shadows’’) as a can-

didate set of geometric parameters for correlation with effective
Young’s moduli. We use the projections on three global coordinate
planes in the original coordinate system of the scanned specimen
(see Fig. 3) and normalize them with respect to the pore volume:
pffiffiffiffiffi
~ i ¼ pA
A ffiffiffiffii ð8Þ
3
V
where Ai and A ~ i are the initial and normalized areas of the pore pro-
jection onto the plane normal to the xi th axis ði ¼ 1; 2; 3Þ, and V is
the volume of considered pore.
Two datasets are considered. The first dataset consists of 8351
pores obtained by processing the microcomputed tomography
(lCT) data for a specimen of carbon/carbon composite, see Drach
et al. (2011, 2013).
The second dataset is based on the lCT scan of a sample of 3D
woven carbon/epoxy composite provided by Albany Engineered
Composites (Rochester, NH). For this dataset, the reconstructed Fig. 4. Reconstructed images of the 3D woven composite sample slices before and
images of the sample slices in the format of 8-bit grayscale VTK after image processing in Fiji.
were processed in the open source image processing software Fiji
(see http://fiji.sc/Fiji) before pore shape extraction. First, a ‘‘Fit all connected objects (in our case, pores) in the input matrix. The
Polynomial’’ filter with parameters ‘‘3  3  3’’ was applied to cor- information included the number of objects and the coordinates
rect for brightness variation within and across the slices. Next, the of the components in the input matrix comprising these objects.
slices were binarized using ‘‘Threshold’’ filter with the value of The criterion by which the components were determined to be
‘‘64’’. In the resulting binary images, white regions represent pores connected to the objects was defined by the connectivity parame-
and black is everything else (matrix and reinforcement), see Fig. 4. ter. In our procedure we used the 6-connected neighborhood con-
The pores smaller than 3 voxels in any dimension were regarded as nectivity, meaning that the matrix component was considered to
noise and not extracted. be a part of the object if it was located in the same row or column
The MATLAB code first described in Drach et al. (2013) was used as one of the object components and is adjacent to it. In the next
for automatic extraction of individual pore meshes and calculation part of the script, individual objects were processed to determine
of pore geometric parameters as follows. VTK file containing series their geometric parameters. As the result of the extraction, a
of black and white images was imported into MATLAB workspace triangular surface mesh was generated and a number of geometric
resulting in a 3D matrix with dimensions equal to the dimensions properties were stored (e.g. linear dimensions, volume, principal
of the VTK file. Components in the matrix have integer values: ‘0’ moments of inertia, projected areas, etc.). The final design
for the black color and ‘1’ for the white color corresponding to space contained 8351 extracted pores for carbon/carbon
the pores. The matrix was then processed using functions from composite and 24,041 extracted pores for 3D woven composite.
MATLAB Image Processing Toolbox to obtain information about The ranges of normalized projected area parameters are given in
Table 1.
The DoE approach was utilized to identify the optimal number
and the combination of pores for two models (carbon/carbon and
3D woven carbon/epoxy pore sets). Using JMP software (SAS
Institute Inc, 2010) we constructed a custom I-optimal experimen-
tal design with the aim of producing statistical models with mini-
mal overall variance across the design space. The number of pores
for the initial design was chosen to be 100, which corresponds to
the average relative prediction variance of 0.05. The complete set
of input variable is provided in Table 2.
For a finite number of pores extracted from the microtomogra-
phy data, it was impossible to find pores with normalized pro-
jected areas that exactly matched the design points produced by
the DoE analysis. Therefore, the pores with the smallest deviation
of the parameters from the designed values were chosen. After
analyzing the available pore data, it was found that the pore
parameter space is bound by a constraint that the sum of the three
parameters must be greater than or equal to 3.15 (note that for a
sphere A ~1 þ A~2 þ A
~ 3 ¼ 3:3 and for a cube with the sides parallel
to the coordinate planes: A ~1 þ A
~2 þ A
~ 3 ¼ 3:0). This constraint is
Fig. 3. Three projected areas (‘‘shadows’’) of a pore. taken into account in the DoE analysis.
 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695 2691

Table 1 of elements on the order of 100,000. Material properties were then

The ranges of the three normalized projected area parameters of the pores from the applied to all matrix elements. Since the results were normalized
carbon/carbon composite and 3D woven carbon/epoxy composite datasets.
by Young’s modulus of the matrix E0 , the value E0 ¼ 1 was used.
Minimum Maximum Average value Midrange value Poisson’s ratio was set to m0 ¼ 0:33 in this study. Six sets of bound-
value value (AVG) (MID) ary conditions in displacements were consecutively applied to sim-
(MIN) (MAX) AVG ¼ MAXþMIN MID ¼ MAXMIN
2 2 ulate three loadcases of uniaxial tension (in the directions of three
Carbon/carbon composite dataset global coordinate axes) and three shear loadcases.
~1
A 1.03 1.79 1.41 0.38 The simulation results’ files were processed via a custom Python
~2
A 1.05 1.75 1.40 0.35
script. The script computed volume averages of all stress compo-
~3
A 0.99 1.59 1.29 0.30 D E
ðkÞ
nents rij for each loadcase k. With the volume averages known,
3D woven carbon/epoxy composite dataset RV
~1
A 0.90 1.35 1.13 0.23 the pore stiffness contribution tensor N ijkl can be calculated: e.g.
~2
A 1.00 1.43 1.22 0.22 from the first loadcase (uniaxial tension in x1 direction):
 D E 
~3 1.20 1.89 1.55 0.35 ð0Þ ð1Þ ð0Þ
A N ij11 ¼ rij  rij =eð0Þ , where rij are the stresses in the
RV
ð1Þ
matrix material subjected to e11 ¼ eð0Þ , hrij iRV are the stress volume
averages in the porous material subjected to e11 ¼ eð0Þ . Pore compli-
Table 2 ance contribution tensor Hijkl was then expressed in terms of N ijkl :
Parameters chosen in the DoE module of JMP software. ð0Þ ð0Þ ð0Þ
Hijkl ¼ Sijmn N mnpq Spqkl , where Sijkl are the components of the compli-
Parameter Value
ance tensor of the matrix material. Finally, parameters E ~2 , E
~1 , E ~3 were
Input variablesa ~1 ; A
A ~2 , A
~3
calculated from the components of H -tensor using formula (5).
Constraints for input variablesb ~1 þ A
A ~2 þ A~ 3 P 3:15
Thus, for each pore characterized by three normalized ‘‘shadow’’
Three model responses ~1 ; E
E ~2 , E
~3
parameters A ~1, A~2 , A
~ 3 , a set of three compliance contribution
Prediction model type 2nd order response surface
~ ~ ~
parameters E1 , E2 , E3 was determined from FEA. These values were
model
Number of experimental runs (design points) 90 entered into JMP software for statistical analysis.
Additional runs (center points) 10
Number of replicates 0 4. Results of the statistical analysis
Number of random starts for search of optimal 1000
designs
The results of the statistical analysis are the sets of coefficients
a
Continuous variables with ranges shown in Table 1. in expression (7) for predictions of E~1 , E
~2 , E
~3 :
b
This constraint is applicable to the considered datasets only.

(a) for carbon/carbon composite dataset

Finite element simulations were performed and results were
^ 1  0:349A
~1 ¼ 2:509 þ 0:844A
E ^ 2  0:195A
^ 3  0:270A ^ 2  0:106A
^1A ^ 3 þ 0:252A
^1 A ^2A^3
processed to determine the components of H -tensor of individual
pores according to the following automated procedure. Surface ^ 1 þ 0:790A
~2 ¼ 2:485  0:362A
E ^ 2  0:201A
^ 3  0:315A ^ 2 þ 0:212A
^1A ^ 3  0:125A
^1 A ^3
^2 A
mesh of a pore was placed in a cube-shaped reference volume with ^ 1  0:205A
~3 ¼ 2:177  0:219A
E ^ 2 þ 0:561A
^ 3 þ 0:209A ^ 2  0:174A
^1A ^ 3  0:163A
^1 A ^3
^2 A
sides five times larger than the largest dimension of the pore, see ð9aÞ
Fig. 5. This setup was auto meshed with second-order tetrahedral
elements. Typical pore shapes yielded FE meshes with the number (b) for 3D woven carbon/epoxy composite dataset
^ 1  0:079A
~1 ¼ 1:873 þ 0:392A
E ^ 2  0:184A
^ 3  0:053A ^ 2  0:101A
^1 A ^ 3 þ 0:064A
^1 A ^3
^2A
~ ^ ^ ^ ^ ^ ^ ^
E2 ¼ 2:086  0:092A1 þ 0:417A2  0:253A3  0:065A1 A2 þ 0:072A1 A3  0:104A2 A3^ ^
^ 1  0:267A
~3 ¼ 3:317  0:226A
E ^ 2 þ 1:148A
^ 3 þ 0:227A
^1 A ^1A
^ 2  0:153A ^2 A
^ 3  0:222A ^3
ð9bÞ
~ AVGðA ~Þ
^i ¼
where A A i i
(see Table 1).
~Þ
MIDðA i

In these models the insignificant factors with p-values (p-value is

the probability that the factor coefficient value is different from zero
due to random error) higher than the level of significance of 0.01 are
excluded sequentially by removing one insignificant parameter at a
time and re-running the model. In both models, the 2nd degree sin-
gle factor parameters turned out to be insignificant, while interac-
tion components were significant in all responses.
Graphically the results for both models are presented in Fig. 6.
In these plots, the straight line represents the perfect fit when
the model predicted values coincide with the actual values
(Actual ¼ Predicted, where the term ‘‘Actual’’ is used for the results
of direct FEA simulations). The two curves below and above the
straight line represent the 95% confidence intervals for the pre-
dicted values. In the case of carbon/carbon composite dataset,
the correlation coefficients and root-mean-square errors for E ~1 ,
~2 , E
E ~3 are a R2 ¼ 0:96; 0:95; 0:96 and RMSE ¼ 5:1%; 5:5%; 3:9%,
correspondingly. In the case of carbon/epoxy dataset,
R2 ¼ 0:96; 0:93; 0:97 and RMSE ¼ 3:2%; 4:3%; 4:4%. Thus, we con-
clude that a good correlation is observed between the normalized
‘‘shadow’’ parameters and Young’s moduli pore contribution
Fig. 5. Reference volume with a pore surface mesh inside. parameters.
2692 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695

4.5 4.5

4 4

3.5 3.5

E1 Actual
E1 Actual

3 3

2.5 2.5

2 2

1.5 1.5
1.5 2 2.5 3 3.5 4 4.5 1.5 2 2.5 3 3.5 4 4.5
E1 Predicted E1 Predicted
RSq=0.96 RMSE=0.1287 RSq=0.96 RMSE=0.0599
4.5 4.5

4 4

3.5 E2 Actual 3.5

E2 Actual

3 3

2.5 2.5

2 2

1.5 1.5
1.5 2 2.5 3 3.5 4 4.5 1.5 2 2.5 3 3.5 4 4.5
E2 Predicted E2 Predicted
RSq=0.95 RMSE=0.1366 RSq=0.93 RMSE=0.0907
4.5 4.5

4 4

3.5 3.5
E3 Actual
E3 Actual

3 3

2.5 2.5

2 2

1.5 1.5
1.5 2 2.5 3 3.5 4 4.5 1.5 2 2.5 3 3.5 4 4.5
E3 Predicted E3 Predicted
RSq=0.96 RMSE=0.084 RSq=0.97 RMSE=0.1493

(a) (b)
~i Actual) vs. model prediction results (E
Fig. 6. The results of response surface fitting presented as plots of direct calculation results (E ~i Predicted) for pore compliance
contribution parameters: (a) carbon/carbon composite dataset; (b) 3D woven carbon/epoxy composite dataset.

For validation of the proposed stochastic 3-factor PCC model Pore geometries were selected from the original datasets with all
given by formulae (9), the estimates of pore compliance contribu- available pores. The results of the validation experiment are shown
tion parameters for each dataset were compared to the direct FEA in Fig. 7. The values of R2 and RMSE for carbon/carbon composite
simulations performed on an independently chosen set of 100 are R2 ¼ 0:96; 0:95; 0:92 and RMSE ¼ 4:6%; 4:7%; 5:3% for E ~1 , E
~2
pores which were not used for the model development. For this and E~3 responses, correspondingly. For 3D woven carbon/epoxy
experiment, the normalized projected area parameters were composite, the values of R2 and RMSE are R2 ¼ 0:96; 0:92; 0:94
selected randomly with uniform distribution over the design space. and RMSE ¼ 2:7%; 4:0%; 4:7% for E ~2 and E
~1 , E ~3 responses,
 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695 2693

4.5 4.5

4 4

3.5 3.5

E1 Actual
E1 Actual

3 3

2.5 2.5

2 2

1.5 1.5
1.5 2 2.5 3 3.5 4 4.5 1.5 2 2.5 3 3.5 4 4.5
E1 Predicted E1 Predicted
R2 = 0.96 RMSE = 0.1152 R2 = 0.96 RMSE = 0.050 3
4.5 4.5

4 4

3.5 3.5
E2 Actual
E2 Actual

3 3

2.5 2.5

2 2

1.5 1.5
1.5 2 2.5 3 3.5 4 4.5 1.5 2 2.5 3 3.5 4 4.5
E2 Predicted E2 Predicted
2
R = 0.95 RMSE = 0.1180 R2 = 0.92 RMSE = 0.0827
4.5 4.5

4 4

3.5 3.5
E3 Actual
E3 Actual

3 3

2.5 2.5

2 2

1.5 1.5
1.5 2 2.5 3 3.5 4 4.5 1.5 2 2.5 3 3.5 4 4.5
E3 Predicted E3 Predicted
R2 = 0.92 RMSE = 0.1165 R2 = 0.94 RMSE = 0.1582

(a) (b)
~1 , E
Fig. 7. Validation of the statistical model for responses E ~2 and E
~3 performed on 100 pores not included in the original model: (a) carbon/carbon composite dataset; (b) 3D
woven carbon/epoxy composite dataset.

correspondingly. These results indicate that the proposed models The mean errors from the validation studies for predictions of E ~1 ,
can be used to estimate with reasonable accuracy the pore compli- ~2 and E
E ~3 are correspondingly 0.1%, 0.2%, 0.4% for 3-factor
ance contribution parameters of the remaining 8251 pores in the model, and 3.5%, 3.0% 1.9% for 4-factor model. Similarly, the
carbon/carbon composite dataset and 23,941 pores in the carbon/ root-mean-square error values are 4.6%, 4.7%, 5.3% for 3-factor
epoxy dataset that were not included in the original experiments. model, and 8.0%, 7.7%, 11.1% for 4-factor model. It can be seen that
The accuracy of the 3-factor model (based on the projected both the bias (mean error) and the scatter (RMSE) are lower for the
areas) for the carbon/carbon composite dataset can be compared model constructed using the normalized projected areas. We con-
to the 4-factor model (based on the principal moments of inertia clude that for the considered carbon/carbon composite dataset, the
and surface-to-volume ratio) presented in Drach et al. (2013). choice of projected areas as the model factors provides more
2694 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695

Table 3
Summary of data on accuracy of the considered statistical models.

Model E1 E2 E3
2 2 2
R ME RMSE R ME RMSE R ME RMSE
4-factor model (Drach et al., 2013) 0.87 0.0000 0.0805 0.75 0.0000 0.1084 0.82 0.0000 0.2836
(0.0%) (5.0%) (0.0%) (5.5%) (0.0%) (10.0%)
4-factor model validation (Drach et al., 2013) 0.69 0.0607 0.1389 0.55 0.0614 0.1579 0.55 0.0536 0.3106
(3.5%) (8.0%) (3.0%) (7.7%) (1.9%) (11.1%)
Proposed 3-factor model (carbon/carbon dataset) 0.96 0.0000 0.1287 0.95 0.0000 0.1366 0.96 0.0000 0.084
(0.0%) (5.1%) (0.0%) (5.5%) (0.0%) (3.9%)
Proposed 3-factor model validation (carbon/carbon dataset) 0.96 0.0030 0.1152 0.95 0.0040 0.1180 0.92 0.0082 0.1165
(0.1%) (4.6%) (0.2%) (4.7%) (0.4%) (5.3%)
Proposed 3-factor model (3D woven carbon/epoxy dataset) 0.96 0.0000 0.0599 0.93 0.0000 0.0907 0.97 0.0000 0.1493
(0.0%) (3.2%) (0.0%) (4.3%) (0.0%) (4.4%)
Proposed 3-factor model validation (3D woven carbon/epoxy dataset) 0.96 0.0044 0.0503 0.92 0.0024 0.0827 0.94 0.0393 0.1582
(0.2%) (2.7%) (0.1%) (4.0%) (1.2%) (4.7%)

accurate predictions of effective Young’s moduli as compared to CMMI-1100409. The microcomputed tomography data was pro-
the model based on the principal moments of inertia and sur- vided by Karlsruhe Institute of Technology (carbon/carbon com-
face-to-volume ratio. Note that another drawback of the 4-factor posite) and Albany Engineered Composites, Inc. (3D woven
model is that orientation of the pore principal axes needs to be carbon/epoxy composite). The New Hampshire Innovation
in the direction of the global axes, otherwise we would need to Research Center is acknowledged for initial support of the collabo-
include orientation vectors in the model. In the model based on ration between the University of New Hampshire and Albany Engi-
‘‘shadows’’, no such adjustment is necessary. Table 3 summarizes neered Composites, Inc.
the data on accuracy of predicting the pore contributions to the The authors would like to thank Philip Ramsey for valuable dis-
overall Young’s moduli based on 4 parameters (principal moments cussions on statistical analysis and design of experiment method-
of inertia and the surface-to-volume ratio) and 3 parameters (nor- ologies. The authors also express their gratitude to Robert
malized projected areas). Zimmerman for his insightful comments on compressibility and
overall elastic behavior of irregular pores. Andrew Drach gratefully
5. Conclusions acknowledges the financial support provided by the ICES Postdoc-
toral Fellowship.
A statistical model based on the normalized projected area
parameters is proposed for correlation of geometry of irregularly References
shaped pores with their contributions to the effective Young’s
moduli of the porous material. The model is applied to large data Arns, C.H., Knackstedt, M.A., Pinczewski, W.V., Garboczi, E.J., 2002. Computation of
sets of pores in two different composite materials. The Young’s linear elastic properties from microtomographic images: methodology and
agreement between theory and experiment. J. Geophys. 67, 1396–1405.
moduli in three orthogonal directions coinciding with the global Benveniste, Y., 1987. A new approach to the application of Mori–Tanaka’s theory in
coordinate system of the specimens are considered. In each case, composite materials. Mech. Mater. 6, 147–157.
one hundred pores were selected to construct the corresponding David, E.C., Zimmerman, R.W., 2011. Elastic moduli of solids containing spheroidal
pores. Int. J. Eng. Sci. 49, 544–560.
models. Design of experiments techniques were utilized to select
Dewey, J.M., 1947. The elastic constants of materials loaded with non-rigid fillers. J.
the pores producing minimal variability of model predictions Appl. Phys. D9, 665–675.
across the entire factor space. The constructed models for three Drach, B., Tsukrov, I., Gross, T., Dietrich, S., Weidenmann, K., Piat, R., Böhlke, T., 2011.
~1 , E
~2 , E
~3 have high values of R2 in both Numerical modeling of carbon/carbon composites with nanotextured matrix
Young’s moduli parameters E
and 3D pores of irregular shapes. Int. J. Solids Struct. 48, 2447–2457.
material systems (0.96, 0.95, 0.96 in carbon/carbon composite Drach, B., Drach, A., Tsukrov, I., 2013. Characterization and statistical modeling of
model, and 0.96, 0.93, 0.97 in 3D woven carbon/epoxy composite irregular porosity in carbon/carbon composites based on X-ray
model), as well as relatively low values of RMSE (5.1%, 5.5% and microtomography data. J. Appl. Math. Mech. 93, 346–366.
Ekneligoda, T.C., Zimmerman, R.W., 2006. Compressibility of two-dimensional
3.9% in carbon/carbon composite model and 3.2%, 4.3% and 4.4% pores having n-fold axes of symmetry. Proc. R. Soc. A Math. Phys. Eng. Sci. 462,
in 3D woven carbon/epoxy composite model). These results indi- 1933–1947.
cate that the proposed models, created based on 100 pore shapes Ekneligoda, T.C., Zimmerman, R.W., 2008a. Boundary perturbation solution for
nearly circular holes and rigid inclusions in an infinite elastic medium. J. Appl.
in each case, can be used to estimate with reasonable accuracy Mech. 75.
the Young’s moduli contribution parameters E ~2 and E
~1 , E ~3 for large Ekneligoda, T.C., Zimmerman, R.W., 2008b. Shear compliance of two-dimensional
data sets of pores, in our case more than 8000 in a sample of car- pores possessing N-fold axis of rotational symmetry. Proc. R. Soc. A Math. Phys.
Eng. Sci. 464, 759–775.
bon/carbon composite and more than 24,000 in a sample of 3D Eroshkin, O., Tsukrov, I., 2005. On micromechanical modeling of particulate
woven carbon/epoxy composite. One of the advantages of the pro- composites with inclusions of various shapes. Int. J. Solids Struct. 42, 409–427.
posed method is that each pore is considered in the global coordi- Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion,
and related problems. Proc. R. Soc. A Math. Phys. Eng. Sci. 241, 376–396.
nate system of the specimen. In contrast to most presently used
Eshelby, J.D., 1959. The elastic field outside an ellipsoidal inclusion. Proc. R. Soc.
micromechanical approaches, there is no need for solving the elas- London Ser. A Math. Phys. Sci. 252, 561–569.
ticity problem in the coordinate system associated with the pore González, C., Segurado, J., LLorca, J., 2004. Numerical simulation of elasto-plastic
deformation of composites: evolution of stress microfields and implications for
followed by the coordinate transformation to the global coordinate
homogenization models. J. Mech. Phys. Solids 52, 1573–1593.
system. Thus, the only information needed to evaluate the pore’s Hashemi, R., Avazmohammadi, R., Shodja, H.M., Weng, G.J., 2009. Composites with
contribution to effective Young’s moduli is its projected areas in superspherical inhomogeneities. Philos. Mag. Lett. 89, 439–451.
the global coordinate system. Hashin, Z., Shtrikman, S., 1962. On some variational principles in anisotropic and
nonhomogeneous elasticity. J. Mech. Phys. Solids 10, 78–96.
Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. J.
Acknowledgements Mech. Phys. Solids 11, 357–372.
Hill, R., 1965. A self-consistent mechanics of composite materials. J. Mech. Phys.
Solids 13, 213–222.
The authors gratefully acknowledge the financial support of the Horii, H., Nemat-Nasser, S., 1983. Overall moduli of solids with microcracks:
National Science Foundation through the grants DMR-0806906 and load-induced anisotropy. J. Mech. Phys. Solids 31, 155–171.
 B. Drach et al. / International Journal of Solids and Structures 51 (2014) 2687–2695 2695

Jasiuk, I., Chen, J., Thorpe, M., 1994. Elastic moduli of two dimensional materials Onaka, S., 2001. Averaged Eshelby tensor and elastic strain energy of a
with polygonal and elliptical holes. Appl. Mech. Rev. 47, 18–28. superspherical inclusion with uniform eigenstrains. Philos. Mag. Lett. 81,
Kachanov, M., Sevostianov, I., 2005. On quantitative characterization of 265–272.
microstructures and effective properties. Int. J. Solids Struct. 42, 309–336. Roberts, A.P., Garboczi, E.J., 2000. Elastic properties of model porous ceramics. J. Am.
Kachanov, M., Tsukrov, I., Shafiro, B., 1994. Effective moduli of solids with cavities of Ceram. Soc. 83, 3041–3048.
various shapes. Appl. Mech. Rev. 47, S151. Rodin, G., 1996. Eshelby’s inclusion problem for polygons and polyhedra. J. Mech.
Lee, J.K., Johnson, W.C., 1978. Calculation of the elastic strain field of a cuboidal Phys. Solids 44, 1977–1995.
precipitate in an anisotropic matrix. Phys. Status Solidi A 46, 267–272. Ryan, T.P., 2007. Modern Experimental Design. John Wiley & Sons, Hoboken, NJ.
Li, M., Ghosh, S., Richmond, O., 1999. An experimental–computational approach to SAS Institute Inc, 2010. JMP 9 Design of Experiments Guide. SAS Institute Inc, Cary,
the investigation of damage evolution in discontinuously reinforced aluminum NC. <www.jmp.com/support/downloads/pdf/jmp902/doe_guide.pdf> (retrieved
matrix composite. Acta Mater. 47, 3515–3532. 28.09.11).
Lubarda, V.A., Markenscoff, X., 1998. On the absence of Eshelby property for non- Sevostianov, I., Giraud, A., 2012. On the compliance contribution tensor for a
ellipsoidal inclusions. Int. J. Solids Struct. 35, 3405–3411. concave superspherical pore. Int. J. Fract. 177, 199–206.
Mackenzie, J.K., 1950. The elastic constants of a solid containing spherical holes. Tsukrov, I., Kachanov, M., 1993. Solids with holes of irregular shapes: effective
Proc. Phys. Soc. (London) 63B, 2–11. moduli and anisotropy. Int. J. Fract. 64, R9–R12.
Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of Tsukrov, I., Novak, J., 2002. Effective elastic properties of solids with defects of
materials with misfitting inclusions. Acta Metall. 21, 571–574. irregular shapes. Int. J. Solids Struct. 39, 1539–1555.
Myers, R.H., Montgomery, D.C., Anderson-Cook, C.M., 2009. Response Surface Wu, T.Te., 1966. The effect of inclusion shape on the elastic moduli of a two-phase
Methodology: Process and Product Optimization using Designed Experiments, material. Int. J. Solids Struct. 2, 1–8.
third ed. John Wiley and Sons, Hoboken, NJ. Zimmerman, R.W., 1986. Compressibility of two-dimensional cavities of various
Nemat-Nasser, S., Hori, M., 1999. Micromechanics: Overall Properties of shapes. J. Appl. Mech. 53, 500–504.
Heterogeneous Solids, second ed. Elsevier Science Publishers. Zimmerman, R.W., 1991. Compressibility of Sandstones. Elsevier.
Nozaki, H., Taya, M., 1997. Elastic fields in a polygon-shaped inclusion with uniform Zou, W., He, Q., Huang, M., Zheng, Q.-S., 2010. Eshelby’s problem of non-elliptical
eigenstrains. J. Appl. Mech. 64, 495–502. inclusions. J. Mech. Phys. Solids 58, 346–372.