materials
ArticleAnalysis of Strain Field Heterogeneity at theMicrostructure Level and Inverse Identification ofComposite Constituents by Means of Digital
Image Correlation
Witold Ogierman *and Grzegorz Kokot

Department of Computational Mechanics and Engineering, Faculty of Mechanical Engineering,
Silesian University of Technology, Konarskiego 18a Street, 44-100 Gliwice, Poland; grzegorz.kokot@polsl.pl* Correspondence: witold.ogierman@polsl.pl; Tel.: +48-32-237-2446
Received: 15 November 2019; Accepted: 6 January 2020; Published: 8 January 2020


Abstract: The present paper is devoted to the theoretical study on the estimation of the full-fieldstrain at the microstructural level of composite materials by means of Digital Image Correlation (DIC).
The main aim of the paper is to investigate the influence of speckle size on the accuracy of the strain
field measurement at the microscale. The DIC analysis was conducted based on artificial speckle
patterns generated numerically and the deformation behavior of the composites was simulated byusing the finite element method (FEM). This approach gives the opportunity to compare the results of
the DIC in terms of speckle size with the reference FEM solution. Moreover, the paper focuses on the
inverse identification of the material constants of the composite constituents by using information
associated with the measured strain field. The inverse problem is solved by using a novel two-step
optimization procedure, which reduces the problem complexity. The feasibility and accuracy of the
proposed approach are presented by analysis of two exemplary microgeometries representing the
microstructures of fiber reinforced composites.
Keywords: Digital Image Correlation (DIC); composite materials; micromechanics; inverse problem;
finite element method (FEM); mechanical properties

1. Introduction
Effective mechanical properties of composite materials depend on the material properties of
the constituents, the shape and volume fraction of the reinforcement as well as the interactions
between the reinforcement and the matrix. The strain field at the microscopic level is inhomogeneous
due to the contrast between the matrix and reinforcement material properties. The strain field
heterogeneity can be analyzed by using various analytical, numerical and experimental methods.
Analytical micromechanical models based on Eshelby’s solution [1] such as Mori–Tanaka [2,3] orself-consistent [4,5] methods allow the estimation of the strain partitioning between the matrix andreinforcement by referring to the per phase average strains. On the other hand, when analyticalmethods are used, it is difficult to consider strain concentrations which are effects of fibers clusteringand nonuniform spatial distribution. The other way is to use numerical methods which can take into
consideration complex geometrical models describing the material microstructure. The most versatile
and widely used numerical method for analysis of the deformation behavior of composites at the
microstructural level is the finite element method (FEM) [6–9]. The other numerical methods which canbe used for this purpose are the boundary element method (BEM) [10,11] and meshless methods [12,13].
Recently also experimental optical methods have been applied for measuring the deformation behavior
at the microstructural level. In this case the most popular method is digital image correlation (DIC)
Materials 2020, 13, 287; doi:10.3390/ma13020287

www.mdpi.com/journal/materials

Materials 2020, 13, 287

2 of 19

which is a computer-based procedure that allows full-field displacement information to be obtainedby recording the motion of a speckle pattern on a specimen surface during the deformation [14–17].
The numerical procedures of the DIC are based on pseudo-affine transformation where series ofpictures collected during the deformation of the specimen are correlated and results are presented
in the form of full field color maps of the displacements or strains [18]. The DIC can be applied for
different length scales since the method has no inherent length scale [14]. However, the accuracy of
the DIC method depends on the quality of the speckle pattern which must be placed on the specimen
surface [19,20]. Literature discusses several different approaches to the preparation of the speckle
pattern for DIC testing at the microscale. The feasibility of the DIC technique for the measurementon micro- and nanoscale was demonstrated by Berfield et al. [14], the random speckle pattern was
created by a fine point airbrush (for microscale measurement) and by solution deposition of fluorescentsilica nanoparticles (for nanoscale measurement). Ghadbeigi et al. measured the strain field at the
microstructure level of dual-phase steel [21] and interstitial-free steel [22], the authors used the natural
pattern of the microstructure in order to carry out the DIC analysis. Anzelotti et al. [23] studied the
behavior of a twill-weave carbon fiber reinforced epoxy lamina, the surface was prepared for the DIC
measurement by using a white ink spray. Canal et al. [24] and Mehdikhani et al. [25] analyzed the
behavior of fiber-reinforced composites, the speckle patterns were created by depositing submicron
alumina particles. Joo et al. [26] measured local strain in a dual phase steel by DIC by basing on the
nanodot patterns.
The present paper is devoted to theoretical study on the influence of the speckle size on theaccuracy of the strain field measurement at the microstructural scale by means of DIC. This is a very
important issue from the practical point view because before capturing the images of undeformed
and deformed specimen a speckle pattern that guarantees exposing the strain field heterogeneity due
to the underlying microstructure is required. Theoretical estimation of the required speckle size canprovide information on how the speckle pattern should be prepared, or which preparation methodshould be used. On the other hand, the error of measurement caused by the non-optimal speckle
pattern can be quantified. The paper presents a numerical procedure of the generation of artificial
speckle patterns and the technique of mapping of these patterns on a virtually deformed specimen.
This work concentrates on unidirectional long fiber reinforced composites. The accuracy of DIC
analysis with respect to different sizes of the speckles is investigated by analyzing microgeometries
containing several fibers. Moreover, the paper presents a novel approach of inverse identification of
the properties of the composite constituents. The idea is to perform an inverse identification of theelastic constants of the composite constituents by basing on the full strain field measurement. This way
of indirect estimation of material properties may be very useful for quantification of the properties ofindividual phases of composites which are difficult to capture during standard experimental testing,for instance in the case of in situ composites. There are several works which are devoted to usage offull field measurement for the purpose of inverse identification of elastic properties of materials at themacroscopic level (constituents are not distinguished in this case) [27–31]. One of the most widely used
identification approaches is the finite element model updating method which is based on minimizing
the discrepancy between the displacements measured experimentally and the displacements computedby the FEM depending on the material constants (variables). Rahmani et al. [32] used the full fieldmeasurement data for identification of the elastic constants of composite constituents, the authors
improved the finite element model updating method by adding a regularization term which assumes
that besides the displacements, the results provided by the micromechanical model are also fitted
to the experimental data. The present paper introduces a novel approach of inverse identification
which is based on minimizing the discrepancy between the per phase average strains measured
experimentally and computed by using FEM. The inverse problem is solved in a framework of new
two-step optimization procedure, which reduces the problem complexity. The feasibility and accuracyof the proposed approach are presented by analysis of two exemplary microgeometries representing
the microstructures of fiber reinforced composites.

Materials 2020, 13, 287
Materials 2020, 13, 287

3 of 19
3 of 19

2. Materials and Methods
2. Materials and Methods
2.1. Artificial Speckle Patterns
2.1. Artificial Speckle Patterns
In order to investigate the influence of the speckle size on the accuracy of the DIC analysis
In orderto investigate
the influenceof thesize ofon thethe patternaccuracygenerationof the DICassumes
analysis artificial
artificial
specklepatterns weregenerated.
Thespeckle
procedurethat thespeckleconsists
patternsofwere
generated.number
The procedureof thepattern
generation
assumes
thatintensity
the image
image
a predefinedof squares(gridof squares)
whose
grayscaleisconsistsofapredefined
numberofsquares(gridofsquares)whosegrayscaleintensityisdetermined
determined in a random manner. Exemplary speckle patterns associated with different sizes of thein a randommanner. Exemplary
specklepatterns
associatedsizesofofartificial
the speckles
arespecklesare presented
in Figure 1.
A morechallenging
issue withis thedifferent
preparation
specklepresentedwhich
in Figure
1.to
Athe
more
challenging
is the preparationof artificial
speckleofpatterns
whichpatterns
relatespecimen
afterissue
deformation,
Figure 2 presentsa scheme
generationofrelateimages.
to the specimen
deformation,
Figure2 presentsa schemeof generationof suchthe
images.
The
such
The first after
step isassociated withcreationof theFE modelwhich describes
geometryfirst constitutive
step is associatedwithof theFE modelwhich
describes
theconditions.
geometry andconstitutive
andbehavior
ofcreation
the material
phasesas wellas theboundary
The resultof thebehaviorof isthea material
phasesas wellas the boundary
conditions.
Theisresultof the FE
is aFE
analysis
displacement
field.
Afterwards
the displacement
field
interpolated
to analysis
the corners
displacement
field.
the displacement
is interpolated
to theof the squares
thatofthe squares
thatAfterwards
built the artificial
speckle field
patterns
representing
thecorners
undeformed
state. The
built the artificial
specklea patterns
the undeformed
state.areThe
coordinatesmodified,describingina
coordinates
describing
locationrepresentingof the cornersof the squares
respectively
location of thethe
cornersof thestatesquaresare respectively
in consequence
state
consequence
deformedis represented
by anmodified,image consistingof a settheof deformed
quadrilaterals.
is represented
by animage consistingof a set of quadrilaterals.
Finally,
representing
Finally,
the images
representing
the undeformed
and deformed state
canthebe images
treated as
input data theforundeformed
and deformed state can be treated as input data for the DIC analysis.
the
DIC analysis.

(a)

(b)

(c)

(d)

Figure 1.
1. Artificial
Artificial speckle
specklepatterns
patternsgenerated
generatedwith
withrespect
respectto
todifferent
different speckle
speckle sizes,
sizes, grids
gridsofofsquares:
squares:(a) 50 ×× 50
50 (b)(b) 100
100 ××100;
100;(c)(c)200200××200;
200;(d)(d)400400××400.
400.

Materials 2020, 13, 287
Materials 2020, 13, 287

4 of 19
4 of 19

Figure 2. Scheme of generation of artificial speckle patters representing the deformed state.
Figure 2. Scheme of generation of artificial speckle patters representing the deformed state.

2.2. Inverse Identification Problem
2.2. Inverse Identification Problem
The aim of the identification procedure is determination of the per phase material properties
The aimidentification
procedure
is determination
phase material
properties
byby basingon ofthethedata
captured duringexperimental
testing ofof thethepercompositematerial.
In the casebasingonphasethe data
captured
duringexperimental
testing moduli
of the composite
In theE caseof aof a two
compositethevariables
are the Young’sof matrix material.
Em and fiber
f and thetwo
phase
compositethevariablesaretheYoung’smoduliofmatrixEm and fiber Ef and the Poisson
Poisson ratios of matrix υm and fiber υf respectively. The method proposed in this paper assumes thatratiosof matrix
υm andoutfiber
υf respectively.
The method
proposedin thisassumesofthat
identification
is carriedby solving
the optimization
problem
associatedwithpaper
minimization
theidentification
is carriedby solving
the optimization
problem
associated withthediscrepancy between
theoutper-phaseaveragestrains measured
experimentally
andminimization
computed by ofusing
discrepancy
betweenper-phaseand computed by using
the theoreticalmodel.theTheobjectiveaveragefunctionstrains
can bemeasured
defined inexperimentallythe following way:
the theoretical model. The objective function can be defined in the following way:
2 X
2 "
2 D ETHE D EDIC
2 #
X
ETHE
D EDIC
2 2 D
THE
DIC 2
THE
DIC 2 min Fmin(Em F(E, Ef ,mυ,mE,f υ, υf )m =−− εεijij+(1)+ εεijij f − −ε ij εij f  , ,, υf ) =   εij εm
(1)
ij m
mm
f
f
=11 j =1 
i=1 ji=

(

) (

)

where
where superscript
superscript DIC
DIC denotes
denotes thethe experimental
experimental results
results andand superscript
superscript THE
THE denotes
denotes thethe results
results
obtained
by
using
the
theoreticalmodel.
The
described
optimization
problemcanbesimplified
obtained by using the theoretical model. The described optimization problem can be simplified bybydividing
dividing thethe optimization
optimization procedure
procedure into
into two
two steps.
steps. In
In thethe casecase ofof twotwo phase
phase composites
composites thethe strainstrainfield
field depends
depends onon thethe υυm,, υυf andandrelationshiprelationshipbetween
betweenthe
theYoung’s
Young’smoduli:
moduli:m

f

E
R =E f ,R = Ef m,
Em

(2)
(2)

therefore,variables cancan bebe reduced
reduced toto three.
three. The
The optimization
optimization constraints
constraints associatedassociatedtherefore, thethe number
number ofof variables
with
the
lower
and
upper
boundsforthevariables
providedinthefollowingway:with the lower and upper bounds for the variables provided in the following way:R ∈ 1,500 and υf , υ m ∈ 0.1, 0.45 .
(3)
R ∈ h1, 500i and υf , υm ∈ h0.1, 0.45i.
(3)should be suitable for most of composite materials. The described problem is ill-posed in theshould be suitable
for most of
materials.
The
described
problem is ill-posedin the Hadamard
Hadamard
sense, therefore
incomposite
order to finda uniquesolution
an appropriate
regularization
term mustsense,
thereforeorder
to find a uniquesolution
an appropriate
must xy
be however
included.
beincluded.
TheinDIC
measurementof strainsat themicroscale isregularization
performed interm
the planeThe
DIC
measurementofstrainsatthemicroscaleisperformed
in
the
plane
xy
howevertheadditional
the additional physical information may be associated with the effective (composite) out-of-plane
properties. For example, Rahmani et al. [32] used out-of-plane stiffness as the additional information.
In order to extend the objective function described by Equation (1) it is more convenient to add

Materials 2020, 13, 287

5 of 19

physical information may be associated with the effective (composite) out-of-plane properties. For
example, Rahmani et al. [32] used out-of-plane stiffness as the additional information. In order toMaterials
2020,
13, 287 function described by Equation (1) it is more convenient to add information related
5 of 19extend theobjective
to the transverse Poisson ratio υc . It is worth mentioning that it could be determined straightforwardly
information related to the transverse Poisson ratio υc. It is worth mentioning that it could be
for instance by using a standard tensile test including two extensometers [33] or DIC (at the macroscopic
determined straightforwardly for instance by using a standard tensile test including twolevel of composite) [33,34]. The discrepancy between the transverse Poisson ratio for the composite
extensometers [33] or DIC (at the macroscopic level of composite) [33,34]. The discrepancy between
determined experimentally and theoretically can be expressed in the form of an additional term
the transverse Poisson ratio for the composite determined experimentally and theoretically can be
included in the objective function. In consequence the objective function for the first optimization step
expressed in the form of an additional term included in the objective function. In consequence thecan be formulated as follows:
objective function for the first optimization step can be formulated as follows:
2 X
2 "D E
2 2 D ETHE D EDIC 22 # X
THETHE D EDIC
2 2 
DIC
THE
DIC

2 2
THE EXP
EXP




THE
min Fmin(R, υF(R,
ε
−
ε
+
ε
−
ε
+
υ
−
υ
(4)
m , υfυ) =
ij
ij
ij
ij
c
c
, υf ) =   εmij
. . (4)− εijm  +  εij f− εij f   + υc − υc
m
m
m 
f
f

 
i =1 j = 1 
i=1 j=


(

)

The
first optimization
optimization step
step isis illustrated
illustrated schematicallyschematically inin Figure
Figure 3.
3.
The first

Figure 3. Scheme of the first optimization step.
Figure 3. Scheme of the first optimization step.

To solve the described optimization problem a global optimization method like for example a
To algorithm
solve the described
optimization
problem
globaloptimization
methodlike foragenetic
can be applied[35–39]. In
the casea ofthe theoreticalmodel,
naturally
theexample
FE modelgenetic
algorithm
canbeapplied[35–39].
In
the
case
of
the
theoreticalmodel,
naturally
the
FE
model
which reflects the actual microstructure can be used. Here, the strains associated with the matrix andwhich
reflects
theobtained
actual microstructure
can beHere,thematrixand:inclusion
can beby averaging overtheused.
volumeof thethe strains
domainassociatedoccupied withby thematrixdΩminclusion can be obtained by averaging over the volume of the domain occupied by the matrix dΩm:
Z
D EFEMεij ε FEM
= = εεij ((EE , E , υ , υ )dΩ ,(5)mij
 ij mm ,Eff , υmm , υff )dΩmm ,
(5)
m

ΩmΩm

and the volume of the domain occupied by the fiber dΩf:
εij

FEM
f

=  εij (E m ,E f , υ m , υf )dΩ f .
Ωf(6)

The FE based modeling of the deformation behavior of the microstructure can provide very good
accuracy however it should be noted that it may lead to very long times of computations due to many

Materials 2020, 13, 287

6 of 19

and the volume of the domain occupied by the fiber dΩf :D EFEM Z
εij=
εij (Em , Ef , υm , υf )dΩf .

(6)

f

Ωf

The FE based modeling of the deformation behavior of the microstructure can provide very good
accuracy however it should be noted that it may lead to very long times of computations due to many
objective function evaluations which are required by the global optimization algorithm. Moreover a
two-dimensional FE model used during the optimization is not able to provide the transverse Poisson
ratio thus a three-dimensional FE model should be applied which in turn may lead to prohibitive timesof computations. On the other hand, the proposed objective function formulated in terms of per-phaseaverage fields allows the application of an efficient mean field homogenization model instead of FEM.
During this work the Mori–Tanaka model [2,3,40–44], which found wide popularity due to goodpredictive capabilities, was used. In this case the strain average over the composite is related to thestrain averages over the matrix and fiber in the following way:hεic = hεim (1 − ff ) + hεif ff ,(7)

where ff is a volume fraction of the fibers. A relationship between the strain average over the matrix
and fibers can be expressed in terms of the strain concentration tensor A:
h 

i−1
hεif = S Cm −1 Cf − I + I hεim ≡ Ahεim(8)where S is the Eshelby tensor, Cm and Cf are elastic stiffness tensors for matrix and fiber respectively.
The Eshelby tensor depends on the Poisson ratio of the matrix and shape of the inclusion, componentsof these tensors for continuous fiber can be found for example in the book of Mura [45]. The strain
averages over the fiber and the matrix depend on the strain average over the composite in thefollowing way:
−1
−1
(9)
hεif = A[ff A + (1 − ff )I] hεic , hεim = [ff A + (1 − ff )I] hεic
The effective stiffness tensor for the composite can be expressed in the following form:
C = Cm + ff (Cf − Cm )A[(1 − ff )I + fA]−1(10)After identifying the υm , υi , and R the second optimization step, whose aim is to determine the
Em and Ei must be carried out. There is only one variable during the second optimization step since
R is known. The objective function is determined as the minimization of the stress average over thecomposite domain computed by using the theoretical model and stress determined experimentally:min F(Em ) =σ11 (Em ) THE− hσ11 iEXP
c
c

2

(11)

It is worth mentioning that the experimental determination of stress (in plane xy) can be carried out
by using the force transducer integrated with the loading device (during the microscopic observation).
The other way, which also provides accurate results, is comparison of the longitudinal Young’s
modulus EL determined experimentally and theoretically (experimental data can be captured during
the transverse Poisson ratio measurement carried out for the purpose of regularization):2min F(Em ) = EL (Em )THE − EEXP
.
L

(12)

The problem related to the second step is well-posed and can be solved by using a gradient based
optimization method very efficiently. The second optimization step is illustrated schematically in
Figure 4.

Materials 2020, 13, 287
Materials 2020, 13, 287

7 of 19
7 of 19

Figure 4. Scheme of the second optimization step.
Figure 4. Scheme of the second optimization step.

3. Results and Discussion
3. Results and Discussion
3.1. Numerical Modeling of the Deformation Behavior
3.1. Numerical Modeling of the Deformation Behavior
For the purpose of the virtual tests, two different microgeometries representing the microstructures
For
the purposeof thevirtualtwo different
microgeometries
representing
the microstructuresof the composites
were
taken
intotests,
account(Figure 5).
Two-dimensional,
plane strain
finite elementofthe composites
were taken
into account(Figure5). computations
Two-dimensional,
planestrain finite
element
modelswere considered.
In orderto performthe FE
ANSYS
software(version
18.2,
modelswere
considered.
In
order
to
performtheFE
computations
ANSYS
software(version
18.2,
Canonsburg, PA, USA) was used. The DIC can be applied for different length scales thus, no particularCanonsburg,
USA) wasThe DICcan be appliedwas
forintroduced,
different lengthscalesthus,
unit of length PA,
associatedwith used.
the modelsof microgeometrya lengthof thesidenoofparticularof length
the modelsof microgeometry
introduced,
a lengthofthe squareunit
relatedto the associatedgeometricalwith
modelswas considered
as unitary.was
A finefinite elementmesh
theside ofofthesquare related
to thegeometricalwas
considered
as unitary.
A fineequals
finiteconsisting
quadrilateralelements
wasgenerated, models
the global
lengthof the edgeof the element
element
mesh6).
consisting
quadrilateralelements
generated,
global lengthof the edge of
0.002 (Figure
Materialof
properties
describing
thewasbehaviorof thetheconstituents
corresponding
to thetheelement
0.002 (Figure6). Material
properties
describing
thebehaviorof the0.35.
constituents
analyzedequals
microgeometries
are collected
in Table
1, the volume
fractionof thefibers equals
Periodiccorresponding
to the[46,47]
analyzed
microgeometries
are εcollected
in Table 1, the volume fraction of theboundary conditions
enforcinguniaxial strain
11 = 0.01 were applied for both microgeometries.
fibersequals strain
0.35. Periodic
enforcinguniaxial
strainε11material= 0.01 were
The obtained
fields areboundarypresented conditionsin Figure 7 [46,47](for material
1) andin Figure8 (for2).
applied for both microgeometries. The obtained strain fields are presented in Figure 7 (for material
Table 1. 2).
Material Properties of the Composite Constituents.
1) and in Figure 8 (for material

Material 1
Material 2

Em , MPa

Ei , MPa

υm

υi

3500
2600

72,000115,000

0.35
0.39

0.20
0.28

Materials 2020, 13, 287

8 of 19

Materials 2020, 13, 287

8 of 19

Materials
Materials 2020,
2020, 13,
13, 287
287

88 of
of 19
19

(a)

(b)

Figure 5. Geometrical models of the composite microstructure: (a) Material 1; (b) Material 2.

(a)
(a)

(b)
(b)

Figure5.
modelscomposite
microstructure:(a)2.
Figure5. Geometrical
Geometricalmodels
models ofofthecomposite
microstructure:(a) Material
Material
1; (b)(b) Material
Material
2. Material 2.
Figure5. Geometrical
ofthe
the
composite
microstructure:(a) 1;
Material
1; (b)

Figure 6. Finite element mesh discretizing the microgeometry.
Table 1. Material Properties of the Composite Constituents.

Em, MPaEi, MPaυmυiMaterial 1
3500
72,0000.35 0.20
Material
2element
2600
115,000the
0.39
0.28
6.
elementmeshmicrogeometry.
Figure6. Finite
Finiteelementmesh discretizing
discretizing
thethemicrogeometry.
FigureFigure
6.
Finitemesh
discretizing
microgeometry.
Table
Table 1.
1. Material
Material Properties
Properties ofof thethe Composite
Composite Constituents.
Constituents.

Material
Material 11
Material
Material 22

EEmm,, MPa
MPa
3500
3500
2600
2600

EEii,, MPa
MPa
72,00072,000115,000
115,000

υυmm
υυii
0.35
0.35 0.20
0.20
0.39
0.39 0.28
0.28

Materials 2020, 13, 287

9 of 19

(a)

(b)

(a)
(a)

(b)
(b)

(c)
Figure 7.fields
Strainobtained
fields obtainedby using
finiteelement
element method(FEM)
for Material
1 (a) ε111; (b)Figure 7. Strain
by using
thethefinite
method(FEM)
for Material(a)ε22ε;11 ; (b) ε22 ;(c) ε12.
(c) ε12 .

(c)Materials
2020,7.13,
287 fields obtained by using the finite element method (FEM) for Material 1 (a) ε11; (b) ε22;9 of 19
FigureStrain(c) ε12.

(b)

(a)

(c)
Figure8.8.Strain
Strainfields
fieldsobtained
obtainedby
byusingusingthetheFEM
FEMforforMaterial
Material22(a)(a)εε1111;; (b)(b) εε2222; (c) ε12.
Figure12 .

The measurement of the transverse Poisson ratio (required as the regularization factor) was also
simulated by using the FEM. In this case, a tensile test of the composite was modeled by assuming athree-dimensional representative
(RVE)
whose
geometry
is presented
in Figure9. Thethree-dimensional
representativevolume
volumeelement
element(RVE)
whosegeometry
is presented
in Figure9.
unit uniaxial strain ε33 = 1 was enforced by prescribing mixed boundary conditions, after solving theboundary value problem the effective transverse Poisson ratio can be determined in the following way:Zυc = −ε11 (Em , Ef , υm , υf )dΩc .

(13)

Ωc

During the tensile test also the longitudinal Young’s modulus can be determined:Z
σ33 (Em , Ef , υm , υf )dΩc .

Ec =
Ωc

(14)

Table 2. Elastic constants determined during the virtual tensile test.

Materials 2020, 13, 287

Material 1
Material 2υc0.28810.3458Ec, MPa
27,734.2
42,335.2

10 of 19

Figure 9. Virtual
Virtual tensile
tensile test
test along
along thethe Z axis, representative volume
volume element
element (RVE)(RVE) representing the
microstructure of the composite.

The Young’s
moduliandCorrelation
Poisson ratios obtained after the FE simulations for both analyzed
3.2. Resultsof the Digital
Imagematerials are collected in Table 2.
The DIC analysis was carried out by using ISTRA 4D software (Dantec Dynamics GmbH, Ulm,
Germany). Four different
patterns determined
were consideredin virtual
accordance
Table 2.speckle
Elastic constants
during thetensilewith
test. those presented in
Figure 1. The coarsest speckle pattern consisting of 50 × 50 quadrilaterals provides a ratio of the fiberυcc , MPa the finest speckle pattern
diameter to the single speckle size equaling approximately
10,Ewhile
Material 1 provides a0.2881consisting of 400 × 400 quadrilaterals
ratio eight times 27,734.2
higher. The results of the analysis
Material
20.345842,335.2
are presented in the form of color maps illustrating the strain distributions. Figures 10–12 present the
distributions of ε11, ε22 and ε12 strains for Material 1. Figures 13–15 present the distributions of ε11, ε223.2. Resultsof theImage
Correlation Moreover, the results are summarized in Tables 3 and 4
andε12 strains
forDigital
Material
2, respectively.
wheretheDICaveragestrains
matrix<εby11>m, <ε22>m and fibers <ε11>f, <ε22>f as well as errors related toThe
analysis
wasforcarriedoutusing ISTRA 4D software (Dantec Dynamics GmbH, Ulm,thereferenceFEsolutionarecollected.
Germany). Four different speckle patterns were considered in accordance with those presented in
Figure 1. The coarsest speckle pattern consisting of 50 × 50 quadrilaterals provides a ratio of thefiber diameter to the single speckle size equaling approximately 10, while the finest speckle patternconsisting of 400 × 400 quadrilaterals provides a ratio eight times higher. The results of the analysis
are presented in the form of color maps illustrating the strain distributions. Figures 10–12 present the
distributions of ε11 , ε22 and ε12 strains for Material 1. Figures 13–15 present the distributions of ε11 , ε22
and ε12 strains for Material 2, respectively. Moreover, the results are summarized in Tables 3 and 4
where the average strains for matrix <ε11 >m , <ε22 >m and fibers <ε11 >f , <ε22 >f as well as errors related
to the reference FE solution are collected.

Materials 2020, 13, 287
Materials 2020,
2020, 13,
13, 287287
Materials

11 of 19
11 of
of 19
19
11

(a)
(a)

(b)
(b)

(c)
(c)

(d)
(d)

Figure10.
Distributionofstrain
11 obtained
obtained forMaterial
withconsiderationofdifferent
speckleFigureFigure10.
10. Distribution
Distributionof
ofstrainstrainεεε11
forMaterial
Material111withwithconsiderationconsiderationof
ofdifferent
differentspeckle
speckle11 obtainedfor
sizes,
gridofsquares:
(a)
50
×
50,
(b)
100
×
100,
(c)
150
×
150,
(d)
200
×
200.
sizes,
150,(d)(d)
200
× 200.
sizes,grid
gridofofsquares:
squares:(a)(a)50
50××50,50,(b)(b)100100× ×100,
100,(c)(c)150
150× ×150,
200
× 200.

(b)
(b)

(a)
(a)
Figure 11. Cont.

Materials 2020, 13, 287
Materials
Materials 2020,
2020, 13,
13, 287287

12 of 19
12
12 of
of 19
19

(c)
(c)

(d)
(d)

FigureFigure11. Distribution
Distributionofstrainobtainedfor
Material
withconsiderationofdifferent
speckleFigure 11.
Distributionof
ofstrainstrainεεε222222obtained
obtainedfor
forMaterial
Material111withwithconsiderationconsiderationof
ofdifferent
differentspeckle
specklesizes,
gridofsquares:
(a)
50
×
50,
(b)
100
×
100,
(c)
150
×
150,
(d)
200
×
200.
sizes,
gridofsquares:
(a)
50
×
50,
(b)
100
×
100,
(c)
150
×
150,
(d)
200
×
200.
sizes, grid of squares: (a) 50 × 50, (b) 100 × 100, (c) 150 × 150, (d) 200 × 200.

(a)
(a)

(b)
(b)

(c)
(c)

(d)
(d)

Figure 12.
12. Distribution
Distributionof
ofstrainstrainεεε121212obtained
obtainedfor
forMaterial
Material111withwithconsiderationconsiderationof
ofdifferent
differentspeckle
speckleFigureFigure12.
Distributionofstrainobtainedfor
Material
withconsiderationofdifferent
specklesizes,grid
gridofofsquares:
squares:(a)(a)50
50×××50,
50,(b)(b)100100×××100,
100,(c)(c)150150×××
150,
(d)
200
×
200.
sizes,
150,(d)
200
×× 200.
sizes,
gridofsquares:
(a)
50
50,
(b)
100
100,
(c)
150
150,
(d)
200
200.

Materials 2020, 13, 287
Materials
Materials 2020,
2020, 13,
13, 287287

13 of 19
13
13 of
of 19
19

(a)
(a)

(b)
(b)

(c)
(c)

(d)
(d)

FigureFigure13.
Distributionofstrain
11 obtainedfor
Material
withconsiderationofdifferent
speckleFigure13.
13. Distribution
Distributionof
ofstrainstrainεεε11
obtainedfor
forMaterial
Material222withwithconsiderationconsiderationof
ofdifferent
differentspeckle
speckle11obtained
sizes,
gridofsquares:
(a)
50
×
50,
(b)
100
×
100,
(c)
150
×
150,
(d)
200
×
200.
sizes,
gridofsquares:
(a)
50
×
50,
(b)
100
×
100,
(c)
150
×
150,
(d)
200
×
200.
sizes, grid of squares: (a) 50 × 50, (b) 100 × 100, (c) 150 × 150, (d) 200 × 200.

(a)
(a)

(b)
(b)
Figure 14. Cont.

Materials 2020, 13, 287
Materials 2020, 13, 287
Materials 2020, 13, 287

14 of 19
14 of 19
14 of 19

(c)
(c)

(d)
(d)

(a)
(a)

(b)
(b)

(c)
(c)

(d)
(d)

FigureFigure 14. Distribution
Distribution ofof strainstrain εε2222obtained
obtainedfor
forMaterial
Material22 with
with considerationconsideration ofof different
different speckle
speckleFigure14.
Distributionofstrainε22obtainedfor
Material2
withconsiderationofdifferent
specklesizes,
gridofsquares:
(a)
50
×
50,
(b)
100
×
100,
(c)
150
×
150,
(d)
200
×
200.
sizes, grid of squares: (a) 50 × 50, (b) 100 × 100, (c) 150 × 150, (d) 200 × 200.
sizes, grid of squares: (a) 50 × 50, (b) 100 × 100, (c) 150 × 150, (d) 200 × 200.

Figure 15. Distribution of strain ε12 obtained for Material 2 with consideration of different speckle
Figure15.
Distribution
εε1212obtained
withofFiguregrid
15. ofDistributionofstrainobtained
forMaterial
Material
withconsiderationof different
different speckle
specklesizes,
squares: (a)of
50strain× 50, (b)100
× 100,for(c)
150
× 150,22(d)200consideration
× 200.
sizes,
200.
sizes,grid
gridofof squares:
squares: (a)(a)50
50××50,50,(b)(b)100100××100,
100,(c)(c)150150× ×150,
150,(d)(d)200200× ×200.

Materials 2020, 13, 287

15 of 19

Table 3. Summary of results of the DIC (Digital Image Correlation) analysis involving different speckle
patterns obtained for Material 1.

FEM (Finite Element Method),reference Solution
DIC, grid 400 × 400
DIC, grid 200 × 200
DIC, grid 100 × 100
DIC, grid 50 × 50-

<ε11 >m<ε11 >f<ε22 >m<ε22 >f0.014710.00124−0.00727−0.00040

0.014193.7%
0.01393
5.3%
0.013488.4%
0.0116920.6%

0.001216.0%
0.001336.9%
0.0013912.1%
0.00319156.3%−0.00697
4.3%−0.00682
6.3%−0.00654
10.0%−0.00577
20.6%−0.00046
16.1%−0.00056
41.9%−0.00052
30.7%−0.00174
340.7%

Table 4. Summary of results of the DIC Analysis involving different speckle patterns obtained forMaterial 2.

FEM, reference solution
DIC, grid 400 × 400
DIC, grid 200 × 200
DIC, grid 100 × 100
DIC, grid 50 × 50-

<ε11 >m<ε11 >f<ε22 >m<ε22 >f0.01505
0.014424.2%
0.01377
8.5%
0.0134310.8%
0.0115523.2%

0.000630.000612.5%
0.00060
4.1%
0.000595.8%
0.00178183.2%−0.00839−0.00796
5.1%−0.00762
9.2%−0.00740
11.7%−0.00650
22.5%−0.00033−0.00031
5.9%−0.00037
13.7%−0.00037
13.6%−0.00180
454.4%

Similar conclusions can be drawn for both analyzed materials. The coarsest speckle pattern does
not allow the strain distribution to be captured properly although it indicates that the strain distribution
is inhomogeneous, the positions of the fibers are exposed. Decreasing the speckle size twice (grid
100 × 100) leads to dramatically better results, in this case the strain distribution is quite similar to thereference one (FE solution) however the peak strains are significantly underestimated and the results
near the fiber-matrix interface are highly averaged. Application of the finest speckle pattern provides a
strain field in good agreement with the reference one, peak strains are only slightly underestimated.
3.3. Results of Inverse Identification
The inverse problem described in Section 2.2 was solved by using the proposed two-step
optimization procedure. The solution of the first optimization step was carried out by using the genetic
algorithm with consideration of a population consisting of 100 chromosomes and 30 generations which
in turn led to 3000 evaluations of the objective function. The second optimization step was solved byusing the gradient based algorithm, in this case 10 evaluations of the objective function formulated by
Equation (12) were enough to find the optimum. The Mori–Tanaka method was used as theoreticalmodel during evaluation of the objective function. At first, in order to verify the proposed identification
procedure, the input data for identification was selected from the results obtained by using the FEM.
Afterwards, the input data was taken from the results obtained by the DIC with consideration ofdifferent sizes of the speckles. The identified elastic constants are collected in Tables 5 and 6 for Material
1 and Material 2 respectively.

Materials 2020, 13, 287

16 of 19

Table 5. Results of the Identification of Constituents Properties for Material 1.

FEM-reference solution
FEM-as input
DIC, grid 400 × 400
DIC, grid 200 × 200
DIC, grid 100 × 100
-

υm

υf

R

Em , MPa

Ef , MPa

0.35000.35932.7%
0.35230.7%
0.35090.3%
0.34491.5%

0.20000.18238.8%
0.1935
3.3%
0.19492.5%
0.2045
2.3%

20.6
21.8
6.0%
18.4
10.4%
15.8
23.4%
13.5
34.5%

3500.03345.64.4%
3902.611.5%
4493.428.4%5162.747.5%

72,000.0
72,936.8
1.3%
71,908.1
0.1%
70,801.7
1.7%
69,565.6
3.4%

Table 6. Results of the identification of constituents properties for Material 2.

FEM-reference solution
FEM-as input
DIC, grid 400 × 400
DIC, grid 200 × 200
DIC, grid 100 × 100
-

υm

υf

R

Em , MPa

Ef , MPa

0.39000.38281.9%0.37483.9%
0.36925.3%
0.36506.4%

0.28000.28973.5%
0.30278.1%
0.310711.0%
0.3175
13.4%

44.2
50.5
14.1%
35.1
20.6%
25.9
41.4%
22.9
48.2%

2600.02312.111.1%
3270.625.8%
4354.867.5%
4886.888.0%

115,000.0
116,646.7
1.4%
114,869.3
0.1%
112,857.5
1.9%
111,872.9
2.7%

An interesting finding is that the Young’s modulus of the fiber was identified with relatively highaccuracy for all considered speckle sizes, the accuracy decreases only slightly with increasing size ofthe speckles. On the other hand, in order to obtain reasonable accuracy of identification of Young’s
modulus of the matrix a fine speckle pattern must be applied, the accuracy decreases significantly withincreasing size of the speckles. The Poisson ratios of matrix and fiber were identified with reasonableaccuracy, in this case accuracy of the identification decreases with increasing size of the speckles moresignificantly than in case of the fiber Young’s modulus. However, the accuracy is much less prone tochanges than in the case of the Young’s modulus of the matrix. It is worth mentioning that in some cases
the results of identification obtained on the basis of data related to the DIC analysis are surprisingly
more accurate than in the case of the input data related to the FEM (“exact” input data). This fact iscaused by differences between the output of the Mori–Tanaka model and the FE analysis based on the
microstructure geometry. In other words, discrepancy between the strain averages obtained by using
the Mori–Tanaka model and the FE analysis may be larger than the discrepancy between the strain
averages obtained by using the Mori–Tanaka model and the DIC analysis (which uses the FE results
as input). Nonetheless, the proposed method of identification provides reasonable accuracy and atime-efficient solution.
4. Concluding Remarks
This paper presents a study on the influence of the speckle size on the accuracy of the strain
field measurement at the microscale by means of DIC. Generation and analysis of artificial images
representing the undeformed and deformed material states gave an opportunity to compare theresults of the DIC in terms of the speckle size with the reference FE solution. Such analysis led to theestablishment of an approximate size of the speckles which is required to observe the heterogenous
strain field at the microstructure level with the desired accuracy. On the other hand, it could be a useful
tool for estimation of the error level of the DIC analysis when a nonoptimal speckle pattern is used.
Moreover, a novel method of inverse identification of the elastic properties of the composite constituentswas proposed. Identification was performed by solving the optimization problem associated with

Materials 2020, 13, 287

17 of 19

minimization of the discrepancy between the per phase average strains measured experimentallyand computed by using the theoretical model. The inverse problem was solved by using a novel
two-step optimization procedure, which reduces the problem complexity. It is worth mentioning thatthe proposed formulation of the objective function allows the use of a time-efficient Mori–Tanaka
model as theoretical model for computing the per-phase average strains which in turn leads to a very
efficient solution of the inverse problem. The measured strain field at the microstructure level does notprovide enough information for identification, therefore an additional regularization term, like forexample the transverse Poisson ratio, must be added. The accuracy of the identification was analyzed
in terms of the size of the speckles; in general the proposed method can provide reasonable accuracyof identification.
Further work will be concerned with the experimental testing of the deformation behavior ofcomposites at the microstructural level. The method of assessment of the required speckle pattern size
proposed in the present paper should be applied in order to perform reliable measurements. In thiscase a digital representation of a real speckle pattern covering a specimen’s surface can be consideredinstead of the artificial speckle pattern. Then the displacements obtained from the FE analysis of thematerial microstructure can be mapped on the digital representation of the real speckle pattern. Finally,
the discrepancy between the reference FE solution and DIC analysis based on the digital representationof the real speckle pattern may be a useful indicator of the accuracy of the DIC analysis carried out on
the basis of images captured during the experiment.
Author Contributions: Conceptualization, W.O.; methodology, W.O. and G.K.; software, W.O.; validation, W.O.
and G.K.; formal analysis, W.O. and G.K.; investigation, W.O. and G.K.; resources, W.O. and G.K.; data curation,
W.O. and G.K.; writing-original draft preparation, W.O.; writing-review and editing, G.K.; visualization, W.O. andG.K.; funding acquisition, W.O. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Science Centre, Poland, grant number 2016/21/N/ST8/01119
and the Faculty of Mechanical Engineering, Silesian University of Technology, grant number 10/990/BK_19/0139.
Conflicts of Interest: The authors declare no conflict of interest.

References
1.
2.
3.
4.
5.
6.
7.
8.

9.
10.
11.

Eshelby, J.D. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R.
Soc. Lond. Ser. A. Math. Phys. Sci. 1957, 241, 376–396.
Benveniste, Y. A new approach to the application of Mori-Tanaka’s theory in composite materials. Mech.
Mater. 1987, 6, 147–157. [CrossRef]
Mori, K.; Tanaka, T. Average stress in matrix and average elastic energy of materials with misfitting inclusions
Title. Acta Metall. 1973, 21, 571–574. [CrossRef]
Benveniste, Y. Revisiting the generalized self-consistent scheme in composites: Clarification of some aspects
and a new formulation. J. Mech. Phys. Solids 2008, 56, 2984–3002. [CrossRef]
Nemat-Nasser, S.; Hori, M. Micromechanics: Overall Properties of Heterogeneous Materials, 2nd ed.; Elsevier
Science Publishers: Amsterdam, The Netherlands, 1999.
Kouznetsova, V.; Brekelmans, W.A.M.; Baaijens, F.P.T. Approach to micro-macro modeling of heterogeneousmaterials. Comput. Mech. 2001, 27, 37–48. [CrossRef]
Segurado, J.; Llorca, J. A numerical approximation to the elastic properties of sphere-reinforced composites.
J. Mech. Phys. Solids 2002, 50, 2107–2121. [CrossRef]
Ganesh, V.V.; Chawla, N. Effect of particle orientation anisotropy on the tensile behavior of metal matrix
composites: Experiments and microstructure-based simulation. Mater. Sci. Eng. A 2005, 391, 342–353.
[CrossRef]
Ogierman, W.; Kokot, G. Generation of representative volume elements of heterogeneous materials with
distributed orientations of inclusions. Compos. Struct. 2018, 201, 636–646. [CrossRef]
Kamiński, M. Boundary element method homogenization of the periodic linear elastic fiber composites. Eng.
Anal. Bound. Elem. 1999, 23, 815–823. [CrossRef]
Fedeliński, P.; Górski, G.; Czyż, T.; Dziatkiewicz, G.; Ptaszny, J. Analysis of effective properties of materials
by using the boundary element method. Arch. Mech. 2014, 66, 19–35.

Materials 2020, 13, 287

12.
13.
14.

15.
16.
17.
18.
19.
20.
21.
22.

23.
24.

25.

26.

27.
28.
29.
30.

31.
32.

33.
34.

18 of 19

Rodrigues, D.E.S.; Belinha, J.; Pires, F.M.A.; Dinis, L.M.J.S.; Jorge, R.M.N. Material homogenization techniquefor composites: A meshless formulation. Sci. Technol. Mater. 2018, 30, 50–59. [CrossRef]
Ahmadi, I.; Aghdam, M.M. Micromechanics of fibrous composites subjected to combined shear and thermalloading using a truly meshless method. Comput. Mech. 2010, 46, 387–398. [CrossRef]
Berfield, T.A.; Patel, J.K.; Shimmin, R.G.; Braun, P.V.; Lambros, J.; Sottos, N.R. Micro-and nanoscale
deformation measurement of surface and internal planes via digital image correlation. Exp. Mech. 2007, 47,
51–62. [CrossRef]
Grédiac, M.; Hild, F.; Pineau, A. Full-Field Measurements and Identification in Solid Mechanics; John Wiley andSons: Hoboken, NJ, USA, 2012.
Barile, C.; Casavola, C.; Pappalettera, G. Digital image correlation comparison of damaged and undamaged
aeronautical CFRPs during compression tests. Materials 2019, 12, 249. [CrossRef] [PubMed]
Kokot, G.; Skalski, K.; Makuch, A.; Ogierman, W. Digital image correlation and nanoindentation in evaluationof material parameters of cancellous bone microstructure. Arch. Mater. Sci. Eng. 2017, 83, 10–16. [CrossRef]
Chu, T.C.; Ranson, W.F.; Sutton, M.A. Applications of digital-image-correlation techniques to experimental
mechanics. Exp. Mech. 1985, 25, 232–244. [CrossRef]
Schreier, H.; Braasch, J.; Sutton, M. Systematic errors in digital image correlation caused by intensity
interpolation. Opt. Eng. 2000, 39, 2915–2921. [CrossRef]
Technical Note, T-Q-400-Accuracy-3DCORR-003-EN. In Error Estimations of 3D Digital Image Correlation
Measurements; Dantec Dynamics GmBH: Ulm, Germany, 2006.
Ghadbeigi, H.; Pinna, C.; Celotto, S.; Yates, J.R. Local plastic strain evolution in a high strength dual-phasesteel. Mater. Sci. Eng. A 2010, 527, 5026–5032. [CrossRef]
Ghadbeigi, H.; Pinna, C.; Celotto, S. Quantitative Strain Analysis of the Large Deformation at the Scale ofMicrostructure: Comparison between Digital Image Correlation and Microgrid Techniques. Exp. Mech. 2012,
52, 1483–1492. [CrossRef]
Anzelotti, G.; Nicoletto, G.; Riva, E. Mesomechanic strain analysis of twill-weave composite lamina underunidirectional in-plane tension. Compos. Part A Appl. Sci. Manuf. 2008, 39, 1294–1301. [CrossRef]
Canal, L.P.; González, C.; Molina-Aldareguía, J.M.; Segurado, J.; Llorca, J. Application of digital image
correlation at the microscale in fiber-reinforced composites. Compos. Part A Appl. Sci. Manuf. 2012, 43,
1630–1638. [CrossRef]
Mehdikhani, M.; Aravand, M.; Sabuncuoglu, B.; Callens, M.G.; Lomov, S.V.; Gorbatikh, L. Full-field strain
measurements at the micro-scale in fiber-reinforced composites using digital image correlation. Compos.
Struct. 2016, 140, 192–201. [CrossRef]
Joo, S.H.; Lee, J.K.; Koo, J.M.; Lee, S.; Suh, D.W.; Kim, H.S. Method for measuring nanoscale local strain
in a dual phase steel using digital image correlation with nanodot patterns. Scr. Mater. 2013, 68, 245–248.
[CrossRef]
Genovese, K.; Lamberti, L.; Pappalettere, C. A new hybrid technique for in-plane characterization oforthotropic materials. Exp. Mech. 2004, 44, 584–592. [CrossRef]
Conceição António, C.; Rasheed, S. A displacement field approach based on FEM-ANN and experiments for
identification of elastic properties of composites. Int. J. Adv. Manuf. Technol. 2018, 95, 4279–4291. [CrossRef]
Lauwagie, T.; Sol, H.; Roebben, G.; Heylen, W.; Shi, Y.; Van der Biest, O. Mixed numerical-experimental
identification of elastic properties of orthotropic metal plates. NDT E Int. 2003, 36, 487–495. [CrossRef]
Avril, S.; Bonnet, M.; Bretelle, A.S.; Grédiac, M.; Hild, F.; Ienny, P.; Latourte, F.; Lemosse, D.; Pagano, S.;
Pagnacco, E.; et al. Overview of identification methods of mechanical parameters based on full-fieldmeasurements. Exp. Mech. 2008, 48, 381–402. [CrossRef]
Barile, C.; Casavola, C.; Pappalettera, G.; Pappalettere, C. Innovative mechanical characterization of materials
by combining ESPI and numerical modelling. Int. J. Mech. 2016, 10, 115–123.
Rahmani, B.; Mortazavi, F.; Villemure, I.; Levesque, M. A new approach to inverse identification of
mechanical properties of composite materials: Regularized model updating. Compos. Struct. 2013, 105,
116–125. [CrossRef]
Carlsson, L.; Adams, D.; Pipes, B. Experimental Characterization of Advanced Composite Materials, 4th ed.; CRC
Press: Boca Raton, FL, USA, 2014.
Pritchard, R.; Lava, P.; Debruyne, D.; Terentjev, E. Precise determination of the Poisson ratio in soft materials
with 2D digital image correlation. Soft Matter 2013, 9, 6037. [CrossRef]

Materials 2020, 13, 287

35.

36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.

47.

19 of 19

Kokot, G.; Orantek, P. The Topology Optimization Using Evolutionary Algorithms. In IUTAM Symposium on
Evolutionary Methods in Mechanics; Burczyński, T., Osyczka, A., Eds.; Springer: Dordrecht, The Netherlands,
2004; pp. 173–186.
Ogierman, W. Inverse Identification of Elastic Properties of Constituents of Discontinuously Reinforced
Composites. Materials 2018, 11, 2332. [CrossRef] [PubMed]
Kuś, W.; Burczyński, T. Bioinspired Algorithms in Multiscale Optimization. In Computer Methods in Mechanics.
Advanced Structured Materials; Kuczma, M., Wilmanski, K., Eds.; Springer: Berlin/Heidelberg, Germany, 2010.
Makowski, P.; Kuś, W. Optimization of bone scaffold structures using experimental and numerical data.
Acta Mech. 2016, 227, 139–149. [CrossRef]
Beluch, W.; Długosz, A. Multiobjective and multiscale optimization of composite materials by means ofevolutionary computations. J. Theor. Appl. Mech. 2016, 54, 397–409. [CrossRef]
Tucker, I.C.; Liang, E. Stiffness predictions for unidirectional short-fibre composites: Review and evaluation.
Compos. Sci. Technol. 1999, 59, 655–671. [CrossRef]
Liu, L.; Huang, Z. A Note on mori-tanaka’s method. Acta Mech. Solida Sin. 2014, 27, 234–244. [CrossRef]
Schjødt-Thomsen, J.; Pyrz, R. The Mori-Tanaka stiffness tensor: Diagonal symmetry, complex fibre orientations
and non-dilute volume fractions. Mech. Mater. 2001, 33, 531–544. [CrossRef]
Böhm, H. (Ed.) Mechanics of Microstructured Materials; Springer: Vienna, Austria, 2004.
Perdahcioǧlu, E.S.; Geijselaers, H.J.M. Constitutive modeling of two phase materials using the mean field
method for homogenization. Int. J. Mater. Form. 2011, 4, 93–102. [CrossRef]
Mura, T. Micromechanics of Defects in Solids; Martinus Nijhoff Publishers: Dordrecht, The Netherlands, 1987.
Liu, H.; Zeng, D.; Li, Y.; Jiang, L. Development of RVE-embedded solid elements model for predicting
effective elastic constants of discontinuous fiber reinforced composites. Mech. Mater. 2016, 93, 109–123.
[CrossRef]
Omairey, S.L.; Dunning, P.D.; Sriramula, S. Development of an ABAQUS plugin tool for periodic RVE
homogenisation. Eng. Comput. 2019, 35, 567–577. [CrossRef]© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open accessarticle distributed under the terms and conditions of the Creative Commons Attribution(CC BY) license (http://creativecommons.org/licenses/by/4.0/).