THE JOURNAL OF CHEMICAL PHYSICS
140
, 074905 (2014)
An optimum approximation of npoint correlation functions of randomheterogeneous material systems
M. Baniassadi,
1,2,a)
M. Safdari,
3
H. Garmestani,
4
S. Ahzi,
2,4
P. H. Geubelle,
3
and Y. Remond
2
1
School of Mechanical Engineering, College of Engineering, University of Tehran,P.O. Box 111554563, Tehran, Iran
2
University of Strasbourg, ICube/CNRS, 2 Rue Boussingault, 67000 Strasbourg, France
3
Aerospace Engineering Department, University of Illinois, 104 S Wright St., Urbana, Illinois 61801, USA
4
School of Materials Science and Engineering, Georgia Institute of Technology, 771 Ferst Drive N.W., Atlanta,Georgia 303320245, USA
(Received 26 December 2013; accepted 5 February 2014; published online 21 February 2014)An approximate solution for npoint correlation functions is developed in this study. In the approximate solution, weight functions are used to connect subsets of (n1)point correlation functionsto estimate the full set of npoint correlation functions. In previous related studies, simple weightfunctions were introduced for the approximation of three and fourpoint correlation functions. Inthis work, the general framework of the weight functions is extended and derived to achieve optimum accuracy for approximate npoint correlation functions. Such approximation can be utilizedto construct global npoint correlation functions for a system when there exist limited informationabout these functions in a subset of space. To verify its accuracy, the new formulation is used toapproximate numerically threepoint correlation functions from the set of twopoint functions directly evaluated from a virtually generated isotropic heterogeneous microstructure representing aparticulate composite system. Similarly, threepoint functions are approximated for an anisotropicglass ﬁber/epoxy composite system and compared to their corresponding reference values calculated from an experimental dataset acquired by computational tomography. Results from both virtual and experimental studies conﬁrm the accuracy of the new approximation. The new formulationcan be utilized to attain a more accurate approximation to global npoint correlation functions forheterogeneous material systems with a hierarchy of length scales.
© 2014 AIP Publishing LLC
.[http://dx.doi.org/10.1063/1.4865966]
I. INTRODUCTION
Statistical description constitutes an essential tool for themathematical representation of heterogeneous materials andtheir reconstruction and homogenization.
1
Statistical continuum approaches have been developed to describe the structural morphology through statistical correlation functions
1–3
and experiments or theoretical models have shown that someeffective properties of heterogeneous materials are stronglydependent on their correlation functions.
1,4–6
Npoint correlation functions can be exploited to characterize the properties of a wide range of heterogeneous systems, includinggalaxies.
1,2,4,5,7,8
For any nphase random homogenous medium, we candeﬁne the following characteristic function for one of thephases, for instance for phase I,
χ
(phase
i
)
(
x
)
=
1 if
x
inphase
i
0 otherwise
,
(1)where
x
is a position vector within the heterogeneous system.The npoint correlation functions are then deﬁned by
9
C
n
(
x
1
,
x
2
,
x
3
, . . . ,
x
n
)
=
χ
α
1
(
x
1
)
χ
α
2
(
x
2
)
. . . .χ
α
n
(
x
n
)
,
(2)
a)
Author to whom correspondence should be addressed. Electronicmail: m.baniassadi@ut.ac.ir. Telephone: (9821) 880204035. Fax: (9821)88013029.
where the superscript
α
i
shows the desired event for the vectorposition
x
i
. In general, for a ddimensional isotropic media,twopoint correlation functions can be extracted from an mdimensional subspace (m
=
1,2,..., d1). For example, twopoint correlation functions can be extracted from a 2d or 1dcut section of the isotropic sample.
10
The twopoint correlations of eigenmicrostructures arecalculated using Fast Fourier Transform.
11,12
In fact, thesquare of the amplitude of the Fast Fourier Transform (FFT)of the characteristic function is equal to the FFT of any autocorrelation of respective microstructure function. Higherorder correlation functions can be computed using higherorder spectra.
8
In the present work, our previously developed approximate solution for npoint correlation functions
13
is modiﬁedto derive a comprehensive equation of the Npoint correlation functions of heterogeneous materials. The main objectiveof this work is to propose a new set of weight functions forthe previously developed formulation.
13
However, our currentwork is not limited to threeand fourpoint correlation functions, as in our previous work, but is extended to higherordercorrelation functions. This approximation can be optimizedfor each microstructure (or heterogeneous system) using neural network or other optimization techniques when thereis limited information available about these functions. To
00219606/2014/140(7)/074905/6/$30.00 © 2014 AIP Publishing LLC
140
, 0749051
0749052 Baniassadi
et al.
J. Chem. Phys.
140
, 074905 (2014)
illustrate this point, the optimum approximation is implemented in the current study to approximate threepoint correlation functions from twopoint function in two differentrandom heterogeneous systems; the ﬁrst one comprising ahardcore platelets/polymer matrix system generated virtually,andthesecond one composed ofaglassﬁber/epoxy composite system processed and assessed experimentally.
II. APPROXIMATION OF NPOINT CORRELATIONFUNCTIONS
For noneigen microstructures, using npoint correlations(n
>
2) of the structure can reveal more morphological details of the features and their distribution. Theoretically, noneigen structures can be uniquely deﬁned by exploiting aninﬁniteorder correlation functions (
n
→ ∞
). In statisticalcontinuum theory, approximating higherorder statistical correlation functions provides a good precision and efﬁciencyfor the homogenization and reconstruction of heterogeneousmicrostructures.In our previous research work,
13
npoint correlation functions have been approximated with the aid of mpoint correlation functions, where 0
<
m
<
n
. In that approximation,npoint correlation functions were divided into n subsets of (n1)point correlation functions. The key idea of this approximation was that the correlation functions can be approximated by a weighted sum of the products of n subsets of (n1)point correlation functions, which are linked using theconditional probabilities. In addition, the weight functionswere calculated using the boundary limiting conditions. Thegeneral formulation of the approximation of npoint correlation functions (
n
>
3) was derived as
13
C
n
(
x
1
,x
2
,x
3
,....,x
n
)
=
n
i
=
1
W
ni
n
−
1
n
−
2
l
=
1
C
(
n
−
1)
(
x
i
,....,x
(
n
−
1)
)
n
−
1
n
−
3
l
=
1
C
(
n
−
2)
(
x
i
,....,x
(
n
−
2)
)
.
n
−
1
n
−
4
l
=
1
C
(
n
−
3)
(
x
i
,....,x
(
n
−
3)
)
n
−
1
n
−
5
l
=
1
C
(
n
−
4)
(
x
i
,....,x
(
n
−
4)
)
...
,
(3)where
W
ni
are the dependency weight functions. In the formulation above, (
x
m
,...
x
m
...,
x
p
) is deﬁned as the subset of
(n
−
1
)
points that include
x
i
as a member of the subset. The weight functions
W
nm
can be calculated using the boundary limits. The ﬁrstboundary condition is given for each
i
aslim
x
(
i
)
→∞
C
n
(
x
1
,x
2
,x
3
,... ,x
(
i
−
1)
,x
(
i
)
,x
(
i
+
1)
....,x
n
)
=
C
1
(
x
i
)
C
n
−
1
(
x
1
,x
2
,x
3
,... ,x
(
i
−
1)
,x
(
i
+
1)
....,x
n
)
.
(4)Here,
C
1
represents the onepoint correlation function. This boundary limit condition can be written aslim
x
(
i
)
→∞
n
i
=
1
W
ni
n
−
1
n
−
2
l
=
1
C
(
n
−
1)
(
x
i
,....,x
(
n
−
1)
)
n
−
1
n
−
3
l
=
1
C
(
n
−
2)
(
x
i
,....,x
(
n
−
2)
)
n
−
1
n
−
4
l
=
1
C
(
n
−
3)
(
x
i
,....,x
(
n
−
3)
)
n
−
1
n
−
5
l
=
1
C
(
n
−
4)
(
x
i
,....,x
(
n
−
4)
)
...
=
C
1
(
x
i
)
C
n
−
1
(
x
1
,x
2
,x
3
,...,x
(
i
−
1)
,x
(
i
+
1)
....,x
n
)
..
(5)Applying this boundary limit, we get
x
i
→∞
:
W
ni
=
0
,W
nj
=
0 for
j
=
i.
(6)The second boundary condition is given bylim
x
1
→∞
...
x
n
→∞
C
n
(
x
1
,x
2
,x
3
,...,x
n
)
=
n
i
=
1
C
1
(
x
i
)
.
(7)This equality condition leads to
x
i
→∞
(
i
=
1
,...,n
)
,
n
i
=
1
W
ni
=
1
.
(8)The third boundary limit is expressed aslim
x
(
i
)
→
x
(
j
)
C
n
(
x
1
,x
2
,x
3
,...,x
(
i
)
,...,x
(
j
)
,....,x
n
)
=
C
n
−
1
(
x
1
,x
2
,x
3
,...,x
(
j
)
,....,x
(
n
−
1)
)
.
(9)From thisboundary condition for compatible events and usingEq. (3), we get
x
i
→
x
j
(for
i
=
j
)
, W
nk
=
0
,k
=
i
and
k
=
j.
(10)Therefore, the necessary conditions for the weight functionsare summarized as follows:
x
i
→∞
W
ni
=
0
, W
nj
=
0 for
j
=
i,
(11)
x
i
→∞
(
i
=
1
,...,n
)
,
n
i
=
1
W
ni
=
1
,
(12)
x
i
→
x
j
(for
i
=
j
)
, W
nk
=
0
, k
=
i
and
k
=
j.
(13)In the proposed approximation, a unique solution does notexist for the weight functions. Therefore, any chosen set of the weight functions that satisfy the necessary boundary limitconditions is useful for this approximation. For example, for
0749053 Baniassadi
et al.
J. Chem. Phys.
140
, 074905 (2014)
the approximation of threepoint correlation functions, a simple choice for the weight functions has been proposed in ourprevious work
13
as
W
3
m
=
x
k
x
l

x
1
x
2
+
x
1
x
3
+
x
2
x
3

m
=
lm
=
k ,
(14)where
m
,
k
,and
l
areequalto1,2,or3.Inthecurrentproposedwork, we extend the approximation for the npoint correlationfunctions by ﬁnding the best weight functions for each systemstructure. First, we show that weight functions for the approximation of 3point, 4point, and 5point correlation functionscan be expressed using CayleyMenger
14
determinant as
W
nm
=
(
−
1)
n
−
1
2
n
−
2
((
n
−
2)!)
2
n
−
1
(
{
x
1
x
2

,

x
1
x
3

,...,

x
(
n
−
2)
x
(
n
−
1)
}
m

n
=
m
)
αm
2
n
k
=
1
(
−
1)
n
−
1
2
n
−
2
((
n
−
2)!)
2
n
−
1
(
{
x
1
x
2

,

x
1
x
3

,...,

x
(
n
−
2)
x
(
n
−
1)
}
k

n
=
k
)
αk
2

n <
6
,
(15)where
x
i
are the position vectors and
x
i
x
j
are the correlationvectors (Figure 1). In Eq. (15), {}
ς
represents the subset of vector lengths of correlations which do not include
x
ς
, theexponents
α
m
and
α
k
are optimization parameters, and
n
−
1
isexpressed using CayleyMenger determinant
14
as
n
−
1
(
{
x
1
x
2

,

x
1
x
3

,

,...,

x
(
n
−
2)
x
(
n
−
1)
}
)
=
0 1
... ...
1
sym
0

x
1
x
2

2

x
1
x
3

2

x
1
x
4

2
...
sym
0

x
2
x
3

2

x
2
x
4

2
......
sym
...............
sym
...
sym sym sym sym
0
(16)We note that
k
gives a formula for (k1)dimensional volume of convex hull of the points (
x
i
) in terms of the Euclideandistances which are deﬁned using the magnitude of vectorlengths for npoint correlation functions with n
<
6 (seeFigure 1). For example, the second polynomial (or
3
) yieldsthe wellknown Heron’s formula for the area A of a triangle with known edge lengths. These proposed weight functions, using CayleyMenger determinant, are no longer physically meaningful for k
>
4 (or n
>
5) because the CayleyMenger determinant in the Euclidean 3D space becomesequal to zero.
15
We also note that the proposed weight function in the work of Baniassadi
et al.
13
, for n
<
5, can also be
FIG. 1. Schematic of correlation vectors of the npoint correlation functions.
used to derive the above CayleyMenger determinantbasedapproximation.Given the limitation of the above proposed approximation for the weight functions, we propose a new generalizedapproximation based on simple weight functions that satisfyall necessary conditions in Eq. (11)–(13) and which are valid
for npoint correlation functions with n
>
2. These are givenby
W
nm
=
((
{
x
1
x
2

x
1
x
3

...

x
(
n
−
2)
x
(
n
−
1)
}
m

n
=
m
))
α
m
n
k
=
1
(
{
x
1
x
2

x
1
x
3

...

x
(
n
−
2)
x
(
n
−
1)
}
k

n
=
k
)
α
k
.
(17)These weight functions are derived by simply consideringmultiplier correlation lengths of the subsets. Unlike the previous approximation (Eq. (15)), the new generalized approximation is a mathematical description with no particular physical meaning, particularly for n
>
3. In fact, we note that, forn
=
3, the previous (Eq. (15)) and new (Eq. (17)) approxima
tions yield the same result, and that Eq. (17) reduces toW
3
m
=
x
k
x
l

α

x
k
x
l

α
+
x
m
x
k

β
+
x
m
x
l

γ
m
=
lm
=
k,
(18)where
α
,
β
, and
γ
are nonzero positive real numbers and
l
,
m
, and
k
are equal to 1, 2, and 3.
FIG. 2. Twopoint correlation functions (TPCF) via correlation length
r
/
d
,where
d
is the diameter of the platelets, for the twophase heterogeneoussystems (shown in inset image).
0749054 Baniassadi
et al.
J. Chem. Phys.
140
, 074905 (2014)
FIG. 3. Dispersion contour of the error for optimization parameters entering Eq. (18): (a)
α
or Xdirection, (b)
β
or Ydirection, and (c)
γ
or Zdirection.
III. COMPUTATIONAL VERIFICATION
To show the validity of the proposed weight functionsapproximation, we use a virtual system based on random microstructures and consider the 3point correlation functionsonly. The virtual microstructure consists of a representativevolume element (RVE) ﬁlled with 4.5% inclusions. Isotropicrealizations of randomly distributed hardcore platelet inclusions, with three different aspect ratios, are generated andused to calculate the optimum threepoint correlation functions of the generated virtual microstructures. In this study, asa ﬁrst step, threepoint correlation functions are approximatedusing Eq. (18). In the next step, optimum values of the optimization parameters (
α
,
β
,and
γ
) are calculated using neuralnetworks. To calculate the twopoint correlation functions, weuse MonteCarlo simulations
2
for the virtual noneigen samples. The obtained results for the twopoint correlation functions are then used to get the approximate solution for thethreepoint correlation functions. Then, we optimize the approximation using Monte Carlo results for the 3point correlation points. For this purpose, we use a three layer neuralnetwork with one input and one output. The network has twohidden layers with 32 and 10 neurons and one output layerwith one neuron and uses LevenbergMarquardt back propagation algorithm for training.The platelet geometries of inclusions are deﬁned by thecorresponding radius center and the normal vector for eachinclusion surface. The size of the RVE is chosen large enoughto produce convergence for the twopoint correlation functions. In Figure 2, we show the calculated twopoint corre
lation functions (TPCF) of the three realized microstructureswith different inclusion aspect ratios.Threepoint correlation functions (platelet/platelet/ platelet) are calculated using Monte Carlo simulation, and theresults are used to optimize the weight functions in Eq. (18).
TABLE I. Optimum values of optimization parameters.Aspect ratio
α β γ
Average of minimum error10 2.2 2.2 2.2 0.1215 2.1 2.1 2.1 0.1320 2 2 2 0.16
In the simulation, more than 1000 threepoint correlation points with different magnitude of vector lengths havebeen computed and used to optimize the threepoint correlation functions (THPCF) given by Eq. (3). In this optimization
process, we deﬁne an error function asError
=
THPCF(
r
)
simulation
−
THPCF(
r
)
approximation

THPCF(
r
)
simulation
.
(19)The average error contours are reported in Figure 3 via theoptimization parameters (
α
,
β
,and
γ
) of the approximation,showing a large dispersion of the error. However, we can ﬁndthe optimum values of these parameters to have the best approximation of the THPCF for the considered microstructure(see Table I).
IV. EXPERIMENTAL VALIDATION
To scrutinize the accuracy of the proposed approximation, a composite specimen composed of 52 vol. % unidirectional glass ﬁbers loaded into an epoxy matrix is analyzed. The internal microstructure of the specimen has beenobtained using a highresolution 3D Xray imaging system
FIG. 4. (a) Sample xray projection image used for reconstruction; (b) 2Dcrosssection view of binary label matrix; and (c) 3D volume rendering of the arrangement of the glass ﬁbers in the unidirectional composite specimen.
0749055 Baniassadi
et al.
J. Chem. Phys.
140
, 074905 (2014)
FIG. 5. Inplane variation of TPCF
01
(a) and TPCF
11
(b), with 1 and 0 respectively denoting the glass and epoxy phases; (c) Variation of the four twopoint correlation functions for a random inplane direction, showing a correlation length
=
10
µ
m for all inplane directions (All dimension are in
µ
m).
(MicroXCT400, Xradia). Figure 4(a) shows a sample 2Dprojection generated by the Xrays passing through the specimen. A number of these Xray projections have been acquiredfrom the specimen in different angles (from
−
170
◦
to 170
◦
around the main axis of the specimen). Using a ﬁltered back projection method, the 3D microstructure of the specimen hasbeenreconstructedfromtheseprojectionimages.Toeliminatenoise and improve quality, a Gaussian smoothing ﬁlter hasbeen applied to the raw data. The binary representation of themicrostructure has been segmented from grayscale data using a threshold ﬁlter. A 2D crosssection of the binary matrixis shown in Figure 4(b). Each voxel of the binary matrix (also
known as label matrix) represents a cubic chunk of the material and a nonzero value is assigned to the each voxel corresponding to the phase occupying the location of the voxel.These operations have been preformed in
MATLAB
using theImageProcessing toolbox. Figure 4(c) also shows a volumetric rendering generated from the acquired data revealing theanisotropic arrangement of unidirectional glass ﬁbers in thecomposite specimen.Twopoint correlation functions (TPCFs) have been evaluated directly from the binary label matrix. Global 3D TPCFhave been constructed by systematically evaluating these twopoint functions in three orthonormal directions over the correlation range. Figures 5(a) and 5(b) show the twopoint corre
TABLE II. Optimum values of optimization parameters.
α β γ
Maximum error1.56 1.77 1.66 0.15FIG. 6. Error distribution in the approximation of THPCF
111
obtained withproposed method, showing a maximum of 15% error.
lationfunctionsTPCF
11
andTPCF
01
with1representingglassphase and 0 the epoxy evaluated along a random inplane direction. Figure 5(c) shows a 2D slice of the global twopointfunctions.Next the threepoint correlation functions (THPCFs) areapproximated from the twopoint functions utilizing Eqs. (3)and (18). For comparison, the reference value of the THPCFs
are also evaluated directly from binary label matrix and compared against their approximate counterparts. An error costfunction is constructed simply by taking absolute differencebetween the approximate and reference values and it is minimized to ﬁnd best values for
α
,
β
, and
γ
. Unique weightingconsidered for all three optimization parameters in the costfunction. Using a collocation approach, a number of pointshave been selected randomly in the correlation range for thespecimen; i.e., about 10
µ
m as shown in Figure 5(c), and the
average value for each parameter are reported in Table II.The distribution of error corresponding to the proposed approximation is depicted in Figure 6 showing a maximum of 15% error.
V. CONCLUSIONS
In the present work, a previously developed approximation for npoint correlation functions based on the conditional probability theory has been modiﬁed. In this study, aset of weight functions has been proposed to obtain an accurate approximation for npoint correlation functions of heterogeneous materials or systems. The approximation can beadapted to different microstructures. Two examples have beenshown to validate our approach. In both examples, threepointcorrelation functions are approximated from the set of twopoint functions using our proposed methodology for optimizing the accuracy. This methodology can be readily extendedto higherorder correlation functions that are needed for applications such as cosmology, biology, and materials science.
1
S. Torquato,
Random Heterogeneous Materials: Microstructure and Macroscopic Properties
(Springer, New York, 2002).
2
M. Baniassadi, A. Laachachi, F. Hassouna, F. Addiego, R. Muller, H.Garmestani, S. Ahzi, V. Toniazzo, and D. Ruch, Compos. Sci. Technol.
71
, 1930 (2011).