Proceedings of the ASME 2011 International Design Engineering Technical Conferences &

Computers and Information in Engineering Conference

IDETC/CIE 2011
August 29-31, 2011, Washington, DC, USA

DETC2011/CIE-47740

ON THE CONSTITUTIVE RESPONSE CHARACTERIZATION FOR COMPOSITE

MATERIALS VIA DATA-DRIVEN DESIGN OPTIMIZATION

John G. Michopoulos John G. Hermanson Athanasios Iliopoulos

Computational Multiphysics Systems Lab. Wood Based Materials & Structures Science Applications Int. Corp.
Center of Computational Material Science Forest Products Laboratory resident at Naval Research Laboratory
Naval Research Laboratory USDA Forest Service Washington DC 20375, USA
Washington, DC, USA Madison, WI, USA

Samuel Lambrakos Tomonari Furukawa

Computational Multiphysics Systems Lab. Computational Multiphysics Lab.
Center of Computational Material Science Dpt. of Mechanical Engineering
Naval Research Laboratory Virginia Polytechnic Institute and State University
Washington, DC, USA Danville, VA, USA

ABSTRACT thetic and actual data demonstrate the successful application of

design optimization for constitutive characterization.
In the present paper we focus on demonstrating the use
of design optimization for the constitutive characterization of
anisotropic material systems such as polymer matrix compos- INTRODUCTION
ites, with or without damage. All approaches are based on the Design optimization as a topic of research relative to en-
availability of experimental data originating from mechatronic gineering applications and product development has been pop-
material testing systems that can expose specimens to multidi- ular within the context of optimal shape determination but not
mensional loading paths and can automate the acquisition of as popular within the context of material characterization. In an
data representing the excitation and response behavior of the attempt to fill this gap, the main objective of the present paper
specimens involved. Material characterization is achieved by is to describe design optimization efforts in the less popular ap-
minimizing the difference between experimentally measured and plication area of the data-driven constitutive characterization of
analytically computed system responses as described by strain anisotropic material systems.
fields and surface strain energy densities. A one dimensional Composite materials are clearly the most widely used
model is presented first to elucidate the design optimization for anisotropic materials for various application areas [1, 2]. Their
the general non-linear constitutive response. Small and large constitutive characterization has been an important topic of in-
strain formulations based on strain energy density decomposi- terest for structural design, material certification and qualifica-
tions are developed and utilized for determining the constitutive tion practitioners. Such characterization has been traditionally
behavior of composite materials. Examples based on both syn- achieved through conventional uniaxial tests and used for deter-

1
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 kind and has continued through the present [8–10]. The most re-
cent prototype of these machines, which is currently under veri-
fication, is shown in Fig. 1.
The energy-based approach associated with mechatronic
testing, although it enables multiaxial loading and inhomoge-
neous states of strain, still requires multiple specimens. It is
significant to state however, that these specimens are tested in a
automated manner with high throughput of specimens per hour,
which have reached values of 30 specimens per hour.
The recent development of flexible full-field displacement
and strain measurements methods has afforded the opportunity
of alternative characterization methodologies [11–14]. Full-field
optical techniques, such as Moire and Speckle Interferometry,
Digital Image Correlation (DIC), and Meshless Random Grid
Method (MRGM), which measure displacement and strain fields
during mechanical tests, have been used mostly for elastic char-
acterization of various materials [14–17]. The resulting measure-
ments are used for identification of constitutive model constants,
via the solution of an appropriately formed inverse problem, with
the help of various computational techniques.
Arguably, the most popular methodology is the mixed nu-
merical/experimental method that identifies the material’s elas-
tic constants by minimizing an objective function formed by
the difference between the full-experimental measurements and
the corresponding analytical model predictions via an optimiza-
tion method [5–10, 14, 17–19]. However, the repetitive finite
element analysis (FEA) required for each iteration of the op-
timization process, makes the computation considerably costly
FIGURE 1: NRL66.3: Most recent 6-DoF mechatronically auto-
[20]. Alternatively, the so-called virtual field method was de-
mated system for the multi-axial testing of composite materials
veloped [20–22] to identify material parameters by finding vir-
tual fields and inversely solving for parameters by substitution of
full-field/surface measurements. That is to say, the virtual field
mining elastic properties. Typically, extraction of these proper- method effectively characterizes materials without finite element
ties, involve uniaxial tests conducted with specimens mounted analysis, provided that appropriate virtual fields are derivable.
on uniaxial testing machines, where the major orthotropic axis Our focus in the present work is to describe recent efforts
of any given specimen is angled relative to the loading direction. concerning design optimization methodologies for constitutive
In addition, specimens are designed such that a homogeneous material characterization. Our approaches are based mostly on
state of strain is developed over a well defined area, which is energy conservation arguments, and they can be classified ac-
for the purpose of measuring kinematic quantities [3, 4]. Con- cording to computational cost in relation to the iterative use of
sequently, the use of uniaxial testing machines imposes require- FEA or not. It is important to clarify that digitally acquired im-
ments of using multiple specimens, griping fixtures and multiple ages are processed by software [23] that implements the MRGM
experiments. The requirement of a homogeneous state of strain [14, 24–28] and is used to extract the full-field displacement and
frequently imposes restrictions on the sizes and shapes of spec- strain field measurements as well as the boundary displacements
imens to be tested. It follows that these requirements result in required for material characterization. Reaction forces and re-
increased cost and time, and consequently to inefficient charac- dundant boundary displacement data are acquired from displace-
terization processes. ment and force sensors integrated with NRL’s multiaxial loader
To address these issues and to extend characterization to called NRL66.3 [29]. In an effort to address the computational
non-linear regimes, multi-degree of freedom automated mecha- cost of the FEA-in-the-loop approaches, the authors have initi-
tronic testing machines, which are capable of loading specimens ated a dissipated and total strain energy density determination
multiaxially in conjunction with energy-based inverse character- approach that has recently been extended to a framework that is
ization methodologies, were introduced at the Naval Research derived from the total potential energy and the energy conserva-
Laboratory (NRL) [5–7]. This introduction was the first of its tion, which can be applied directly with full field strain measure-

2
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
ment for characterization [30–35]. conclusions are presented.
Two techniques, built upon this framework, have been pro-
posed to identify elastic constants and to develop non-parametric
constitutive models of anisotropic materials. THE CASE OF ONE DIMENSIONAL MATERIAL SYS-
The first identification technique estimates the elastic con- TEM
stants for every set of measurements by equating the variation of To introduce design optimization for material characteriza-
the external work, derived from the boundary displacement/force tion in a fashion of increasing complexity we consider first a one-
measurements, with that of the induced strain energy, derived dimensional system that possesses both linear recoverable and
from the full-field strain measurements, and stochastically cor- non-linear irrecoverable responses. The template approach that
recting the estimation using Kalman filter [35]. This technique we will employ first for this simple case and then for more re-
has been proven to identify the elastic constants of anisotropic alistic cases is that of defining a Strain Energy Density function
materials even under the presence of considerable noise in the (SED) that governs the material behavior and indirectly contains
measurements. the actual constitutive behavior.
The second technique develops non-parametric representa-
tions of constitutive models using artificial neural networks [36].
R I
When we first explored this approach in the early nineties [37], U1D = U1D (C, ε) +U1D (C, βi ; ε) =
it was established that the computational performance of the 
1 2
R 
1 2 I

(1)
approach, as it was implemented on the Aspirin/MIGRAINES = Cε + D(βi ; ε) Cε
2 1D 2 1D
neural net simulator framework [38], was not practical for the
amounts of data generated by NRL’s multidimensional testing R (C, ε),U I (C, β ; ε) are the recoverable (elastic) SED
machines. Subsequently, we have applied it on many other ma- where U1D 1D i

terial characterization applications [39–43]. We have established and the irrecoverable (inelastic or dissipated) SED, respectively.
how improvements in the evolving computational technologies The quantities C, βi , D(C, βi ; ε)are the stiffness constant (modu-
have impacted positively the feasibility of more demanding ap- lus of elasticity), the constants participating in the dissipated en-
plications. In that both ANN implementations and full field mea- ergy coefficient function, and the dissipated energy density co-
surement techniques have matured, we have decided to apply efficient function, itself, respectively. Equation 1 implies that
ANN technologies for material characterization [32]. The er- the irrecoverable or dissipated strain energy density (DSED)
ror between the energy quantities is used to develop the neural has been constructed to be a multiplicative decomposition with
network constitutive model, unlike the conventional techniques weighting D(βi ; ε) for the recoverable SED according to:
where stress data are required for the modeling [44–47]. This
technique allows the nonlinear constitutive relations to be mod-
I R
eled comprehensively without the limitations imposed by the U1D (C, βi ; ε) = D(βi ; ε)U1D (C, ε) (2)
parametric expressions of the conventional material models. Ac-
cordingly, this technique has been applied to model the damage The functional form of D(βi ; ε) should be one that ensures en-
behavior of composite materials [48, 49]. ergy dissipation in a manner that yields a softening of nonlin-
In order to maintain reasonable scope this paper considers ear stress-strain constitutive response. There are many forms
only methodologies that require FEA-in-the-loop because of the that have this property, based on transcendental functions, which
simplicity of their implementation and exhaustive capability to have been used in the past [50]. Here we employ a functional
determine the material parameters. The consideration of other form that can be expanded in a Taylor series. This functional
methodologies is more appropriate for future comparative stud- form provides for a polynomial representation that is necessary
ies. condition for algebraic reducibility and is therefore convenient
In the section that follows we present the case of deter- for the application of transformations. This form of the DSED,
mining the properties of the one-dimensional non-linear system. which is initially negligible and then monotonically increasing,
This is done mainly for instructive purposes, which bare rele- can be represented by the following physically consistent choice
vance to subtleties of subsequent formulations presented. Next, of D(βi ; ε):
we present a small strain formulation (SSF) of the general strain
energy density approach followed by a finite strain formulation  m
(FSF), which is for the case of linear and non-linear constitu- 1
− me ε
εf
tive behavior of composite materials with or without damage. D(βi ; ε) = D(m, ε f ; ε) = 1 − e (3)
The paper continuous with a numerical application of design
optimization implementations based on these two formulations, where β1 = m, β2 = ε f are the two material parameters control-
which are in turn based on both synthetic and actual data. Finally ling the dissipative nature of the material behavior. The exponent

3
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 (a) Distribution of elastic, dissipated (or irrecoverable) and total SEDs for one
FIGURE 2: Difference of 1- and 2-term approximations of nor- dimensional case
malized irrecoverable or dissipated energy density as a function
of strain relative to the exact model

m controls the level of non-linearity involved in the model and

ε f expresses the level of failure strain. The series expansion of
D̄(m, ε f ; ε) defined by Eq. 3 yields

n k  mk
(−1)k

1 ε
D̄(m, ε f ; ε) = ∑ . (4)
k=1 k! me εf

Given Eq. 4 and 1, it follows that Eq. 2 may be expressed by

(b) Stress distributions as a function of strain

k  mk
1 n (−1)k

1 ε
I
U1D (C, βi ; ε) = C ∑ ε2 (5) FIGURE 3: SED and corresponding stress distribution
2 k=1 k! me εf

A plot of the function defined by Eq. 5, normalized by the

For the case of generalized hyperelasticity as special case of
Elastic constant C, is presented in Fig. 2 for m = 4 and ε f =
which is elasticity, the corresponding constitutive law will be
0.0008. Referring to Fig. 2 it can be seen that for more than two
given by
terms the series expression gives essentially identical results to
those of the exact evaluation of DSED. It is interesting to note
that the single term expression is also in good agreement with
R (C, ε) +U I (C, β ; ε)]
∂ [U1D
the exact form. Therefore, we can truncate all but the first term σ1D = 1D i
=
in Eq. 5 to obtain: ∂ε
R I
= σ1D (C, ε) + σ1D (C, βi ; ε) = (8)
!
ε m+1
 
1

1

ε
m 2+m
I
U1D (C, βi ; ε) = − C ε2 (6) = Cε −C
2 me εf 2me εm f

It follows then that, Eq. 1 can be expressed The resulting constitutive law contains the linear elastic part,
as expected, but modified by a non-linear inelastic term. An in-
R I dicative variation of total SED and its components (recoverable
U1D = U1D (C, ε) +U1D (C, βi ; ε) = and irrecoverable SEDs), as described by Eq. 7 are shown in Fig.
   m
1 1 1 ε (7) 3(a). Figure 3(b) shows the corresponding stress distribution de-
= Cε 2 − C ε2
2 2 me εf fined by the constitutive law expressed by Eq. 8.

4
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 In the present formulation, the material parameters to be
identified, or estimated based on experimental data, are the
stiffness parameter and the two dissipated energy parameters I1 = tr(ε) = εii ,
C, m, ε f . 1 1
The results of applying design optimization for determining I2 = tr(ε 2 ) = εi j εi j , (10)
2 2
these material parameters for 20 data points created synthetically 1 3 1
by the material model defined Eq. 7 or 8 are shown in Table 1 I3 = tr(ε ) = εi j ε jk εki
2 3
for various methods.
These results indicate that Nelder-Mead and Differential The invariants are chosen to take advantage of a priori knowledge
evolution Methods are able to determine the properties exactly, that the SED as a scalar quantity, must be invariant under frame
while Simulated Annealing takes much longer and achieves less reference translations and rotations. This follows in that the SED
accurate parameter determinations. should be objective (i.e. independent of the observer’s frame of
reference).
An additive decomposition of this expression in terms of a
recoverable and an irrecoverable SED can be expressed by
COMPOSITE MATERIAL SYSTEM
For the general case of a composite material system we con- R I
USSF = USSF (S; εi j ) +USSF (D; εi j ). (11)
sider that a modified anisotropic hyperelastic strain energy den-
sity function can be constructed to encapsulate both the elastic
and the inelastic responses of the material. However, certain Clearly, all the second order monomials of strain components
classes of composite materials reach failure after small strains R (S; ε ) and the higher
will be forming the recoverable part USSF ij
and some under large strains. For this reason we give two ex- order monomials will be responsible for the irrecoverable part
amples, one involving a small (infinitesimal) strain formulation I (D; ε ). The resulting constitutive law is given by
USSF ij
(SSF) and another involving an finite (large) strain formulation
(FSF).
R I
σi j = ∂USSF /∂ εi j = ∂ (USSF (S; εi j ) +USSF (D; εi j ))/∂ εi j =
∂USSF ∂ I1 ∂USSF ∂ I2 ∂USSF ∂ I3
Small Strain Formulation For the SSF we introduce a SED = + +
∂ I1 ∂ εi j ∂ I2 ∂ εi j ∂ I3 ∂ εi j
function that, in its most general form, can be represented as a
scaled Taylor expansion of the Helmholtz free energy of a de- (12)
formable body, which is in terms of small strain invariants of the
form A general expression which provides a strain dependent version
of Eq. 11, is given by

R I
1 1 USSF = USSF (S; εi j ) +USSF (D; εi j ) =
USSF (I1 , I2 , I3 ) = s0 + s11 I1 + s12 I12 + s13 I13 +
2 3! 1 (13)
1 1 1 = si jkl εi j εkl + di jkl (εi j )εi j εkl
+ s14 I14 + . . . + s21 I2 + s22 I22 + s23 I23 + 2
4! 2 3!
1 1 1 1 where si jkl are the components of the elastic stiffness ten-
+ s24 I24 + . . . + s13 I32 + s13 I33 + s14 I34 +
4! 2 3! 4! (9) sor (Hooke’s tensor) and di jkl (εi j ) are strain-dependent damage
1 1 1 functions, which fully define irrecoverable or dissipated strain
. . . + s1121 I1 I2 + s1221 I12 I2 + s1122 I1 I22 + . . .
2 3! 3! energy density given by enforcing the dissipative nature of en-
1 1 1 ergy density. The quantity di jkl (εi j ) can be defined in a manner
+ s1131 I1 I3 + s1331 I12 I3 + s1133 I1 I32 + . . . +
2 3! 3! analogous to that employed for the 1D system described above
1 1 1 and is given by
s1321 I13 I2 + s1221 I12 I22 + s1123 I1 I23 + . . .
4! 4! 4!

ε p
(− qi j ) i j /(epi j )
where the invariants are defined by di jkl (εi j ) = si jkl (1 − e ij ) (14)

5
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 TABLE 1: Design optimization results for a non-linear 1D material system

Method. Objective Function C [GPa] m εf

Target values 0.00 130.00 4. 0.0008
Nelder Mead -0.0000116825 130.00 4. 0.0008
Differential Evolution -0.00000429153 130.00 4. 0.0008
Simulated Annealing 134234 131.21 3.1 0.000856

As in the 1D case, we perform an equivalent series expan-

sion and subsequently drop all terms except the first, in that it
captures almost all of the dissipative behavior. Accordingly,     
σxx s˘xx sxy sxz 0 0 0 εxx
 σyy   sxy
   s˘yy syz 0 0 0   εyy 
 
 σzz   sxz syz s˘zz 0 0 0    εzz 
 
 σxz  =  0 (18)
  
mp 0 0 s˘xz 0 0 
 εxz 

εi j i j
n  m   
1  σyz   0 0 0 0 s˘yz 0 εyz 
∑ (−1)m
 
di jkl (εi j ) = si jkl mp =
m=1 epi j m!qi j i j σxy 0 0 0 0 0 s˘xy εxy
p 2pi j
εi ji j
   2
1 1 εi j
= si jkl (− p + 2p
− · · · +) ' (15) where:
epi j qi ji j epi j 2qi j i j
  pi j
1 εi j
s˘i j = si j 1 − d¯i j


' −si jkl (19)
epi j qipji j

and
Thus the irrecoverable part of the energy in Eq. 11 becomes:
1 1+pi j
d¯i j = pi j εi j (20)
epi j qi j
I I
USSF (D; εi j ) = USSF (si jkl , pi j , qi j ; εi j ) =
1 1+pi j (16) All terms that are not shown in expression 18 are zero due to
= −si jkl pi j εi j εkl the orthotropic symmetry requirements. Therefore, the material
e(2 + pi j )pi j qi j
parameters are the 9 elastic si j constants and 6 × 2 = 12 damage
constants pi j , qi j for a total of 21 parameters. Clearly, when the
quantities d¯i j do not depend on the strains and they are constants,
Next, substituting Eq. 14 into Eq. 11 yields
Eq. 16 reduces to most of the continuous damage theories given
by various investigators in the past [47,51–53]. For a transversely
isotropic material the number of material parameters drops to
R I
5+10=15 for a 3D state of strain and to 4+8=12 for a plane stress
USSF = USSF (S; εi j ) +USSF (D; εi j ) = state.
1 1 1+pi j (17)
= si jkl εi j εkl − si jkl pi j εi j εkl
2 e(2 + pi j )pi j qi j Finite Strain Formulation The FSF can be written in a double
additive decomposition manner. The first being the decomposi-
tion of the recoverable and irrecoverable SED, and the second
Applying Eq. 10 on Eq. 15, and employing Voight [3] notation being the decomposition between the volumetric (or dilatational)
for the case of a general orthotropic material, yields the constitu- Wv and the distortional (or isochoric) Wd parts of the total SED.
tive relation This decomposition is expressed by:

6
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 The corresponding constitutive behavior is given by the second
Piola-Kirchhoff stress tensor according to [51]
R I
USSF = USSF (αi ; J, C̄) +USSF (αi , βi ; J, C̄) =
= [Wv (J) +Wd (C̄, A ⊗ A, B ⊗ B)]− (21)
∂UFSF
−[dvWv (J) + dd Wd (C̄, A ⊗ A, B ⊗ B)] S=2 (26)
∂C
where αi , βi are the elastic and inelastic material parameters of
or the usual Cauchy stress tensor according to
the system, respectively. A rearrangement of these decompo-
sitions, such as the volumetric vs. distortional decomposition,
which appears on the highest expression level, leads to an ex-
pression introduced in [54], i.e., 2 ∂UFSF T
σFSF = F · ·F . (27)
J ∂C

UFSF = (1 − dv )Wv (J) + (1 − dd )Wd (C̄, A ⊗ A, B ⊗ B), (22) Under the FSF formulation the material characterization problem
involves determining the 36 coefficients (at most) of all mono-
mials when the sums in the expression of distortional SED are
with the damage parameters dk ∈ [0, 1], k ∈ [v, d] defined as expanded in Eq. 24, in addition to the compressibility constant d
and the 4 parameters used in Eq. 23. It follows that potentially
 
a (t)
 there can be a total of 41 material constants.
∞ − k
dk = dka 1 − e ηka (23)

NUMERICAL RESULTS
where ak (t) = max Wko (s) is the maximum energy component For the purpose of demonstrating numerically the afore-
s∈[0,t] mentioned concepts, the material selected for generating the
reached so far, and dka ∞ ,η
are two pairs of parameters control-
ka necessary simulated experimental data is a typical laminate
ling the energy dissipation characteristics of the two components constructed from an epoxy resin/fiber laminae system of type
of SED. In this formulation, J = det F is the Deformation Gra- AS4/3506-1. The elastic moduli of this material are listed in
dient, C̄ = F T F is the right Cauchy Green (Green deformation) Table 2 according to several sources [3, 55–57].
tensor, A, B are constitutive material directions in the undeformed Clearly, what is considered to be a set of material constants
configuration, and A ⊗ A, B ⊗ B are microstructure structural ten- varies widely as it really depends on the fiber volume fraction,
sors expressing fiber directions. Each of the two components of the fiber coating, the manufacturing process of the fiber, resin and
SED are defined as composite and the quality of the experimental procedure over-
all. As can be seen from the bottom of the entries of the table
where we added some statistical observations, % deviation ob-
1
Wv (J) = (J − 1)2Wd (C̄, A ⊗ A, B ⊗ B) = served varies from 11.1 % to 78.2 %. It is therefore important
d to identify the set of elastic material properties before and af-
3 3 6
ter a batch of new material is manufactured or before a material
= ∑ ai (I¯1 − 3)i + ∑ b j (I¯2 − 3) j + ∑ ck (I¯4 − 1)k +
i=1 j=1 k=1 system is used for design, material qualification or material cer-
6 6 6 (24) tification. To demonstrate the usage of the SSF in conjunction
+ ∑ dl (I¯5 − 1)l + ∑ em (I¯6 − 1)m + ∑ fn (I¯7 − 1)n + with design optimization we present here an example of using
l=2 m=2 n=2 real data form a multiaxially loaded specimen from a test con-
6 ducted by utilizing NRL66.3. The model characteristics of the
+ ∑ go (I¯8 − (A · B)2 )o specimen used are presented in Fig. 4 where in Fig. 4(a) the
o=2 discretization model and potential boundary conditions are de-
picted, and in Fig. 4(b), a detail at the area of the left notch
where the strain invariants are defined as follows: shows a stacking of [+60,-60]16 with each lamina made out of
AS4/3506-1.
Two objective functions were constructed. Both utilized the
I¯1 = trC̄, I¯2 = 21 (tr2C̄ − trC̄2 ) fact that through the REMDIS-3D software, developed by our
I¯4 = A · C̄B, I¯5 = A · C̄2 B (25) group, one can obtain full field measurements of the displace-
I¯6 = B · C̄B, I¯7 = B · C̄2 B, I¯8 = (A · B)A · C̄B ment and strain fields over any deformable body [14, 23–28].

7
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 TABLE 2: Engineering properties of AS4/3506-1 laminae

Ref. E11 [GPa] E22 [GPa] E33 [GPa] ν12 ν23 ν13 G12 [GPa] G23 [GPa] G13 [GPa]
[3] 147.0 10.3 10.3 0.27 0.28 0.27 7.0 4.04 7.0
[55] 142.0 9.8 9.8 0.30 0.34 0.30 6.0 3.77 6.0
[56] 135.0 9.0 9.0 0.28 0.28 6.9 6.9
[57] 138.0 9.7 9.7 0.30 0.49 0.30 5.24 3.24 5.24
[58] 139.3 11.1 11.1 0.30 0.40 0.30 6.0 3.964 6.0
[59] 150.0 8.0 8.0 0.30 0.30 5.0 5.0
[60] 147.0 10.3 10.3 0.27 0.27 6.89 6.89
[61] 142.0 10.3 10.3 0.27 0.27 7.2 7.2
Min 135.0 8.0 8.0 0.27 0.28 0.27 5.0 3.24 5.0
Avg. 142.5 9.8 9.8 0.29 0.30 0.29 6.279 3.753 6.279
Max. 150.0 11.1 11.1 0.30 0.49 0.30 7.2 4.039 7.2
Deviation [%] 11.1 38.8 38.8 11.1 78.2 11.1 44.0 24.7 44.0
Avg. Cons. 142.5 9.81 9.81 0.29 0.30 0.29 6.279 3.769 6.279

Thus, our experimental measurements for the formation of the TABLE 3: Engineering properties of AS4/3506-1 laminae
objective functions were chosen to be the strains at the nodal
points of the discretization shown Fig. 4(a). The first objective
function chosen was based entirely on strains and is given by Ref E11 [GPa] E22 [GPa] ν12 ν23 G12 [GPa]
Daniels [3] 147.0 10.3 0.27 0.28 7.0
Present 125.0 10.8 0.27 0.32 7.96
i 2
!
N 2
h 2i h
exp
Jε = ∑ ∑ ∑ εi j − εifjem k
,
k
(28)
k=1 i=1 j=i

strains respectively. An implementation of both objective func-

the second objective function is given in terms of surface strain tions was applied by utilizing an implementation of the DIRECT
energy density according to global optimizer [62], which is available within Matlab [63], and
a custom developed Monte-Carlo optimizer, also implemented in
I
2 Matlab. The loading conditions applied for the moveable edge of
JU ≈ U exp −U f em dS ≈
∂Ω the specimen were ux = 0[m], uy = 0.0005[m], uz = 0.001[m], rx =
I 2 2 2 2   2
! 0[rad], ry = 0[rad], rz = 0[rad]. In Fig.5 we can see the experi-
exp exp f em f em mentally acquired distribution of εyy (left) and the corresponding
≈ ∑∑ ∑ ∑ si jmn εi j εmn − si jmn εi j εmn dS
∂Ω i=1 j=i m=1 n=m distribution produced from FEM analysis for the five identified
(29) unique elastic constants, which are shown in Table 3, in compar-
ison to those of [3].
h
exp
i h i For the case of the FSF a two stage optimization was per-
where εi j , εifjem are the experimentally determined and formed. In the first stage we determined the values of the coef-
k k
the FEM produced components of strain are at node k. The quan- ficients of the strain invariant monomials in Eq. 24 such that the
tities U exp ,U f em are the surface strain energy densities formu- FSF matches the SSF by constructing and minimizing an objec-
lated by using the experimental strains and the FEM produced tive function of the form

8
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 (a) Determined experimentally via MRGM

(a) Discretization and boundary conditions of typical specimen used for ma-
terial characterization

(b) Corresponding matched distribution of the numerical model as it was pro-

duced by FEA

FIGURE 5: Distributions of the vertical component of strain

(εyy ) for the purpose of determining the five elastic constants of
AS4/3506-1 composite lamina material

damage behavior of the FSF were determined through the mini-

mization of the objective function

I
2
1 2
(b) Detail of FEA model near one of the notches
JUF (d, a1 , b1 , c1 , d1 , da∞ , ηa ) ≈ UFSF −UFSF dS ≈
∂Ω
 !2  (31)
I 2 2 2 2
FIGURE 4: Finite Elements Model 1 2


≈ ∑∑ ∑ ∑ UFSF −UFSF dS
∂Ω i=1 j=i m=1 n=m

1 (d,
where UFSF ¯ ā1 , b̄1 , c̄1 , d¯1 , d¯a∞ , η̄a ) is the SED of the FSF for the
I
2 parameters d, ¯ ā1 , b̄1 , c̄1 , d¯1 computed from the previous step of
JURF (d, a1 , b1 , c1 , d1 ) ≈ R
UFSF R
−USSF dS ≈ matching the FSF with the SSF. The parameters d¯a∞ , η̄a were cho-
∂Ω
 !2  (30) sen arbitrarily to be some values that can capture nonlinearity due
I 2 2 2 2
R R

to damage. This complete set of parameters represents the known
≈ ∑∑ ∑ ∑ UFSF −USSF dS
∂Ω i=1 j=i m=1 n=m target model of the material as shown in the second row of Table
4. The unknown constitutive model was chosen to be represented
by another FSF SED, namely UFSF 1 (d, a , b , c , d , d ∞ , η ). The
1 1 1 1 a a
This was done is order to establish the proper parameters of bottom row of Table 4 shows the parameters determined by mini-
the FSF model that match the SSF model. mizing the objective function in Eq. 31 via a Monte-Carlo global
In the second stage the material parameters encoding the optimizer.

9
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 TABLE 4: FSF parameters of AS4/3506-1 laminae with damage

Model d a1 b1 c1 e1 da∞ ηa
Target 0.00007012 1611 713 29212 -458.7 0.9 1.8
Identified 0.00006987 1613 721 29105 -449.3 0.92 1.979

(a) Small strain formulation FEA results

FIGURE 7: Comparison of the load vs. vertical component of

strain (εyy ) at a point in front of the notch, between the target and
the identified model by using the FSF

energy density and full field strain based approaches have been
utilized to incorporate massive full field strain measurements
from specimens loaded by a multiaxial custom-made loading
machine. We have formulated objective functions expressing
the difference between the experimentally observed behavior of
composite materials under various loading conditions, and the
(b) Finite strain formulation FEA results simulated behavior via FEA, which are formulated in terms of
strain energy density functions of a particular structure under
FIGURE 6: FEA results of the vertical component of strain (εyy ) identical loading conditions.
for the case of a specimen loaded both in tension and torsion Two formalisms involving small strains and finite strains
have been utilized in a manner that involves both additive decom-
position of recoverable and irrecoverable strain energy density.
An indication of how well the FSF can capture the SSF, for This was done in order to address both the elastic and inelastic
the case of the recoverable (linear elastic) regime, is shown in response of composite materials due to damage. The finite strain
Fig. 6, where the distribution of εyy is shown for both models formulation further involves a volumetric and distortional energy
for a case of combined torsion and tension of the specimen. In decomposition.
addition, a comparison of load vs. strain evolution, for a point in Demonstrations have been given in terms of numerical ex-
front of the first of two notches, is shown in Fig. 7. amples utilizing both synthetic and actual data in determining
both the elastic and inelastic material parameters.
We are currently working in extending these energy formu-
CONCLUSIONS AND DISCUSSION lations to approaches that achieve two main goals. First, that
We have demonstrated the application of design optimiza- those approaches do not require FEA in the loop and second, that
tion methodologies for the determination of the constitutive re- they incorporate the stochastic nature of the material response
sponse of composite materials with or without damage. Strain and the acquired data in a manner that quantifies uncertainty.

10
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
This quantification should utilize prior knowledge via Bayesian “Towards the robotic characterization of the constitutive
based stochastic formalisms, which enable incremental and re- response of composite materials”. Composite Structures,
cursive algorithms based on Kalman filtering and in general, in- 86(1-3), pp. 154–164.
formation theory. [11] Pierron, F., Vert, G., Burguete, R., Avril, S., Rotinat,
R., and Wisnom, M. R., 2007. “Identification of the or-
thotropic elastic stiffnesses of composites with the virtual
ACKNOWLEDGMENT fields method: Sensitivity study and experimental valida-
The authors acknowledge the support by the Office of Naval tion”. Strain, 43(3), pp. 250–259.
Research, and by NRL’s 6.1 core-program. [12] Cooreman, S., Lecompte, D., Sol, H., Vantomme, J., and
Debruyne, D., 2007. “Elasto-plastic material parameter
identification by inverse methods: Calculation of the sen-
REFERENCES sitivity matrix”. International Journal of Solids and Struc-
[1] Kaczmar, J. W., Pietrzak, K., and Wlosinski, W., 2000. tures, 44(13), pp. 4329 – 4341.
“The production and application of metal matrix compos- [13] Furukawa, T., and Michopoulos, J., 2008. “Computational
ite materials”. Journal of Materials Processing Technology, design of multiaxial tests for anisotropic material charac-
106(1-3), pp. 58 – 67. terization”. International Journal for Numerical Methods
[2] Baker, A., Dutton, S., and Kelly, D., 2004. Composite in Engineering, 74(12), pp. 1872–1895.
Materials for Aircraft Structures, Second Edition, 2nd ed. [14] Michopoulos, J. G., Iliopoulos, A. P., and Furukawa, T.,
AIAA, Oct. 2009. “Accuracy of inverse composite laminate characteri-
[3] Daniel, I. M., and Ishai, O., 2005. Engineering Mechanics zation via the mesh free random grid method”. ASME Con-
of Composite Materials, 2nd ed. Oxford University Press, ference Proceedings, 2009(48999), pp. 367–374.
USA, July. [15] Bruno, L., Furgiuele, F. M., Pagnotta, L., and Poggialini,
[4] Carlsson, L., Adams, D. F., and Pipes, R. B., 1986. Exper- A., 2002. “A full-field approach for the elastic characteri-
imental Characterization of Advanced Composite Materi- zation of anisotropic materials”. Optics and Lasers in En-
als, Third Edition, 3 ed. CRC Press, Dec. gineering, 37(4), pp. 417 – 431.
[5] Mast, P., Beaubien, L., Clifford, M., Mulville, D., Sut- [16] Schmidt, T., Tyson, J., and Galanulis, K., 2003. “Full–
ton, S., Thomas, R., Tirosh, J., and Wolock, I., 1983. field dynamic displacement and strain measurement using
“A semi-automated in-plane loader for materials test- advanced 3d image correlation photogrammetry: Part 1”.
ing”. Experimental Mechanics, 23, pp. 236–241. Experimental Techniques, 27(3), pp. 47–50.
10.1007/BF02320415. [17] Kajberg, J., and Lindkvist, G., 2004. “Characterisation of
[6] Mast, P. W., Nash, G. E., Michopoulos, J. G., Thomas, materials subjected to large strains by inverse modelling
R. W., Badaliance, R., and Wolock, I., 1992. Experimen- based on in-plane displacement fields”. International Jour-
tal determination of dissipated energy density as a measure nal of Solids and Structures, 41(13), pp. 3439 – 3459.
of strain-induced damage in composites. Tech. Rep. Tec. [18] Meuwissen, M. H. H., Oomens, C. W. J., Baaijens, F.
Rpt. NRL/FR/6383–92-9369, Naval Research Laboratory, P. T., Petterson, R., and Janssen, J. D., 1998. “Determi-
Washington, DC, USA. nation of the elasto-plastic properties of aluminium using a
[7] Mast, P., Nash, G., Michopoulos, J., Thomas, R., Wolock, mixed numerical-experimental method”. Journal of Mate-
I., and Badaliance, R., 1995. “Characterization of strain- rials Processing Technology, 75(1-3), pp. 204 – 211.
induced damage in composites based on the dissipated en- [19] Molimard, J., Le Riche, R., Vautrin, A., and Lee, J., 2005.
ergy density Part I. Basic scheme and formulation”. Theo- “Identification of the four orthotropic plate stiffnesses using
retical and Applied Fracture Mechanics, 22(2), pp. 71–96. a single open-hole tensile test”. Experimental Mechanics,
[8] Michopoulos, J., 2003. Recent advances in Composite 45, pp. 404–411. 10.1007/BF02427987.
Materials: In Honor Of S.A. Paipetis. Kluwer Academic [20] Grediac, M., Toussaint, E., and Pierron, F., 2002. “Spe-
Press, ch. Computational and Mechatronic Automation of cial virtual fields for the direct determination of material
Multiphysics Research for Structural and Material Systems, parameters with the virtual fields method. 1–principle and
pp. 9–21. definition”. International Journal of Solids and Structures,
[9] Michopoulos, J., 2004. “Mechatronically automated char- 39(10), pp. 2691 – 2705.
acterization of material constitutive respone”. In Proceed- [21] Grediac, M., and Pierron, F., 2006. “Applying the virtual
ings of the 6th World Congress on Computational Mechan- fields method to the identification of elasto-plastic constitu-
ics (WCCM-VI), Tsinghua University Press and Springer, tive parameters”. International Journal of Plasticity, 22(4),
pp. 486–491. pp. 602 – 627.
[10] Michopoulos, J., Hermanson, J., and Furukawa, T., 2008. [22] Chalal, H., Avril, S., Pierron, F., and Meraghni, F., 2006.

11
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
 “Experimental identification of a nonlinear model for com- Adelaide 2008.
posites using the grid technique coupled to the virtual fields [33] Furukawa, T., and Michopoulos, J., 2008. “Online plan-
method”. Composites Part A: Applied Science and Manu- ning of multiaxial loading paths for elastic material iden-
facturing, 37(2), pp. 315 – 325. CompTest 2004. tification”. Computer Methods in Applied Mechanics and
[23] Michopoulos, J. G., and Iliopoulos, A. P., 2009. “A com- Engineering, 197(9-12), pp. 885–901.
putational workbench for remote full field 2d displacement [34] Michopoulos, J., Furukawa, T., and Lambrakos, S., 2009.
and strain measurements”. ASME Conference Proceedings, “Data-driven characterization of composites based on vir-
2009(48999), pp. 55–63. tual deterministic and noisy multiaxial data”. In 2008
[24] Andrianopoulos, N. P., and Iliopoulos, A. P., 2006. “Dis- ASME International Design Engineering Technical Confer-
placements measurement in irregularly bounded plates us- ences and Computers and Information in Engineering Con-
ing mesh free methods”. In Proc. 16th European Confer- ference, DETC 2008, Vol. 3, pp. 1095–1106.
ence of Fracture, Alexandroupolis, Greece, July 3–7. [35] Furukawa, T., and Pan, J. W., 2010. “Stochastic identifica-
[25] Iliopoulos, A. P., Michopoulos, J. G., and Andrianopoulos, tion of elastic constants for anisotropic materials”. Inter-
N. P., 2008. “Performance sensitivity analysis of the Mesh- national Journal for Numerical Methods in Engineering,
Free Random Grid method for whole field strain measure- 81(4), pp. 429–452.
ments”. ASME Conference Proceedings, 2008(43277), [36] Samarasinghe, S., 2006. Neural Networks for Applied Sci-
pp. 545–555. ences and Engineering: From Fundamentals to Complex
[26] Iliopoulos, A. P., Michopoulos, J. G., and Andrianopoulos, Pattern Recognition, 1 ed. Auerbach Publications, Sept.
N. P., 2010. “Performance analysis of the Mesh-Free Ran- [37] Michopoulos, J., 1992. Constitutive characterization of
dom Grid Method for full-field synthetic strain measure- composites via artificial neural nets. Tech. rep., Code 6383
ments”. Strain, p. In Print. Internal Branch Report, NRL, April.
[27] Iliopoulos, A., and Michopoulos, J., 27–31 Jul 2009. “Sen- [38] Leighton, R., 1994. Aspirin/migraines, v7.0.
sitivity analysis of the mesh-free random grid method for [39] Yagawa, G., Yoshimura, S., Okuda, H., and Furukawa, T.,
measuring deformation fields on composites”. In Pro- 1996. Localized Damage IV: Computer-Aided Assessment
ceedings of the 17th International Conference on Compos- and Control. ch. Innovative approaches for modelling of in-
ite Materials, ICCM-17, The British Composites Society: elastic material behaviors (applications of neural networks
Edinburgh International Convention Centre (EICC), Edin- and evolutionary algorithms), pp. 169–184.
burgh, UK, 2009. [40] Furukawa, T., and Yagawa, G., 1998. “Implicit constitu-
[28] Iliopoulos, A. P., and Michopoulos, J. G., 2009. “Effects tive modelling for viscoplasticity using neural networks”.
of anisotropy on the performance sensitivity of the Mesh- International Journal for Numerical Methods in Engineer-
Free random grid method for whole field strain measure- ing, 43(2), pp. 195–219.
ment”. ASME Conference Proceedings, 2009(48999), Jan., [41] Furukawa, T., and Hoffman, M., 2004. “Accurate cyclic
pp. 65–74. plastic analysis using a neural network material model”.
[29] Michopoulos, J. G., Hermanson, J. C., and Iliopoulos, Engineering Analysis with Boundary Elements, 28(3),
A., 2010. “Towards a recursive hexapod for the multi- pp. 195 – 204. Inverse Problems.
dimensional mechanical testing of composites”. In Pro- [42] Man, H., Furukawa, T., Hoffman, M., Imlao, S., and Lou,
ceedings of the ASME 2010 International Design En- Z. H., 2007. “An indirect implicit model for frequency de-
gineering Technical Conferences & Computers and In- pendent hystereses of piezoelectric ceramics”. Journal of
formation in Engineering Conference IDETC/CIE 2010, the Australian Ceramics Society, 43(2), pp. 169–174.
no. DETC2010/CIE-28699. [43] Man, H., Furukawa, T., Hoffman, M., and Imlao, S., 2008.
[30] Michopoulos, J., and Furukawa, T., 2007. “Multi-level cou- “An indirect implicit technique for modelling piezoelec-
pling of dynamic data-driven experimentation with material tric ceramics”. Computational Materials Science, 43(4),
identification”. Lecture Notes in Computer Science (includ- pp. 629 – 640.
ing subseries Lecture Notes in Artificial Intelligence and [44] Ghaboussi, J., and Sidarta, D., 1998. “New nested adap-
Lecture Notes in Bioinformatics), 4487 LNCS, pp. 1180– tive neural networks (NANN) for constitutive modeling”.
1188. Computers and Geotechnics, 22(1), pp. 29 – 52.
[31] Furukawa, T., Michopoulos, J., and Kelly, D., 2008. “Elas- [45] Shin, H. S., and Pande, G. N., 2003. “Identification of
tic characterization of laminated composites based on mul- elastic constants for orthotropic materials from a structural
tiaxial tests”. Composite Structures, 86(1-3), pp. 269–278. test”. Computers and Geotechnics, 30(7), pp. 571 – 577.
[32] Man, H., Furukawa, T., and Kellermann, D., 2008. “Im- [46] Jung, S., and Ghaboussi, J., 2006. “Neural network con-
plicit constitutive modeling based on the energy princi- stitutive model for rate-dependent materials”. Comput.
ples”. In Proceedings of XXII ICTAM, 25-59 August 2008, Struct., 84, June, pp. 955–963.

12
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.
[47] Haj-Ali, R., and Kim, H.-K., 2007. “Nonlinear constitu- rate-dependent compressive strength of carbon/epoxy com-
tive models for FRP composites using artificial neural net- posites”. Composites Science and Technology, 65(15-16),
works”. Mechanics of Materials, 39(12), pp. 1035 – 1042. pp. 2481 – 2491. 20th Anniversary Special Issue.
[48] Man, H., Furukawa, T., Herszberg, I., and Prusty, G., 18- [62] Jones, D. R., Perttunen, C. D., and Stuckman, B. E., 1993.
21 May 2009, Baltimore, MD 2009. “Implicit modeling of “Lipschitzian optimization without the Lipschitz constant”.
damage behavior for composite materials using an energy- J. Optim. Theory Appl., 79(1), pp. 157–181.
based method”. In SAMPE ’09 Spring Symposium Con- [63] The Mathworks. Matlab. http://www.mathworks.com.
ference Proceedings, p. 54(Paper ID: 5704).
[49] Man, H., Furukawa, T., Pan, J., Iliopoulos, A., Michopou-
los, J. G., Orifici, A., and Hermanson, J. “Experimen-
tally validated neural network constitutive modelling using
energy-based characterisation”. submitted.
[50] Matzenmiller, A., Lubliner, J., and Taylor, R. L., 1995.
“A constitutive model for anisotropic damage in fiber-
composites”. Mechanics of Materials, 20(2), pp. 125 –
152.
[51] Haupt, P., 2002. Continuum Mechanics and Theory of Ma-
terials, 2nd ed. Springer, Apr.
[52] Van Der Meer, F., and Sluys, L., 2009. “Continuum models
for the analysis of progressive failure in composite lami-
nates”. Journal of Composite Materials, 43(20), pp. 2131–
2156.
[53] Ju, J. W., Chaboche, J., and Voyiadjis, G. Z., 1998. Dam-
age Mechanics in Engineering Materials, 1st ed. Elsevier
Science, Jan.
[54] Kaliske, M., 2000. “A formulation of elasticity and vis-
coelasticiy for fibre reinforced material at small and finite
strains”. Comput. Meth. Appl. Mech. Eng., 185, pp. 225–
243.
[55] Ryan, K. F. Dynamic response of
graphite/epoxy plates subjected to impact loading.
http://dspace.mit.edu/handle/1721.1/42982. Thesis
(M.S.)–Massachusetts Institute of Technology, Dept. of
Aeronautics and Astronautics, 1990.
[56] Wegner, P. M., and Adams, D. F., 2000. Verification of
the combined load compression (CLC) test methods. Tech.
Rep. DOT/FAA/AR-00/26, DOT/FAA.
[57] Vaidya, R., and Sun, C., 1997. “Fracture criterion for
notched thin composite laminates”. AIAA J., 35(2),
pp. 311–316.
[58] Minguet, P. J., Dugundji, J., and Lagace, P., 1989. “Post-
buckling behavior of laminated plates using a direct energy-
minimization technique”. AIAA Journal, 27(12), pp. 1785–
1792.
[59] Lessard, L. B., Eilers, O. P., and Shokrieh, M. M., 1995.
“Testing of in-plane shear properties under fatigue load-
ing”. Journal of reinforced plastics and composites, 14(9),
pp. 965–987.
[60] Esp, B., and Chan, W., 2008. “Two parameter approxi-
mation to the orthotropic stress concentration factor”. In
SAMPE 2008 - Long Beach, CA May 18 - 22.
[61] Bing, Q., and Sun, C., 2005. “Modeling and testing strain

13
This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
Approved for public release; distribution is unlimited.