T. KROUPA et al.: ONE-DIMENSIONAL ELASTO-PLASTIC MATERIAL MODEL WITH DAMAGE ...
213–218
ONE-DIMENSIONAL ELASTO-PLASTIC MATERIAL MODEL
WITH DAMAGE FOR A QUICK IDENTIFICATION OF THE
MATERIAL PROPERTIES
ENODIMENZIJSKI MODEL ELASTOPLASTI^NEGA MATERIALA
S PO[KODBO ZA HITRO UGOTOVITEV LASTNOSTI MATERIALA
Tomá{ Kroupa, Hana Srbová, Jan Klesa
University of West Bohemia, NTIS – New Technologies for the Information Society, Univerzitní 22, 306 14, Plzeò, Czech Republic
kroupa@kme.zcu.cz
Prejem rokopisa – received: 2015-07-01; sprejem za objavo – accepted for publication: 2016-02-10
doi:10.17222/mit.2015.173
Material parameters for elastic, plastic and damage behavior of low-molecular-weight epoxy resin CHS-EPOXY 520 hardened
with CHS-P 11 are identified in the paper. Uniaxial cyclic static tests were performed on specimens in line with the ASTM
standard D638 - 10 using a Zwick Roell/Z050 test machine with a BTC-Exaclbi.001 clip-on biaxial extensometer. Poisson’s
ratio and tensile strength were calculated directly from the test results. Other parameters were identified using a combination of
a one-dimensional material model, which assumes the infinitesimal strain theory and takes into account elastic, plastic and
damage behavior and the optimization method. The model is coded in Python and the values of the material parameters are
identified using the optiSLang optimization software. The optimization process uses simple design improvement and
gradient-based algorithms with the goal to minimize the difference between the force-elongation curves recorded during the
tests and those calculated with the proposed model.
Keywords: epoxy, identification, optimization, elasticity, plasticity, damage
V ~lanku so opredeljeni parametri epoksi smole z majhno molekulsko maso CHS-EPOXY 520, utrjeno z CHS-p 11, pri
elasti~nem, plasti~nem obna{anju in pri po{kodbi. Enoosni stati~ni cikli~ni preizkusi so bili izvedeni na vzorcih, ki ustrezajo
ASTM standardu D638-10, z uporabo Zwick Roell/Z050 preizkusnega stroja z name{~enim dvoosnim ekstenziometrom
BTC-Exaclbi.001. Poissonov koli~nik in natezna trdnost sta bili izra~unani neposredno iz rezultatov preizkusov. Drugi
parametri so ugotovljeni s kombinacijo enodimenzijskega modela materiala, ki predpostavlja infinitezimalno teorijo napetosti in
ki upo{teva obna{anje v elasti~nem, plasti~nem in ob po{kodbi ter metodo optimizacije. Model je kodiran v Pythonu in
vrednosti parametrov materiala so ugotovljene z uporabo optiSLang programske opreme za optimizacijo. Proces optimizacije
uporablja enostavno izbolj{anje oblike in algoritme, ki temeljijo na gradientu, s ciljem po zmanj{anju razlike med krivuljo sila –
raztezek zabele`ene med preizkusom in izra~unane z uporabo predlaganega modela.
Klju~ne besede: epoksi, identifikacija, optimizacija, elasti~nost, plasti~nost, po{kodba
1 INTRODUCTION
A model simulating one-dimensional pure tensile,
compressive or shear tests is proposed. This approach is
able to predict a one-dimensional response of materials
with elastic, plastic and damage behavior. The model is
based on the idea to reduce the computational time
necessary to identify characteristic parameters of the
tested material. The presented model is not necessarily,
or primarily, meant to be used for the final identification
of material parameters; it can only be used for investi-
gating the starting values or boundaries of the parameters
in the subsequent identifications performed using, e.g., a
finite-element model1,2 or other models.3 The model is
written in free-licensed programing language Python
using only widely used and common modules numpy,
matplotlib, os and pickle. Pre-processing of the experi-
mental results for the further identification process is
also performed in Python. The paper presents results for
pure epoxy resin; nevertheless, the model is not limited
to such type of material.
2 SPECIMENS AND THE EXPERIMENT
Nine specimens shaped according to ASTM standard
D638-10 (Figure 1) were cut from a rectangular plate
with dimensions of approx. 140 mm × 250 mm. The
specimens were made of low-molecular-weight epoxy
resin CHS-EPOXY 520 hardened with CHS-P 11. The
Materiali in tehnologije / Materials and technology 51 (2017) 2, 213–218 213
MATERIALI IN TEHNOLOGIJE/MATERIALS AND TECHNOLOGY (1967–2017) – 50 LET/50 YEARS
UDK 66.017:620.1:620.169.2 ISSN 1580-2949
Original scientific article/Izvirni znanstveni ~lanek MTAEC9, 51(2)213(2017)
Figure 1: Dimensions of the specimens
Slika 1: Dimenzije vzorcev
air from the mixture was removed with a vacuum pump
and the resin was then cured at atmospheric pressure,
which increased the quality of the epoxy resin by
reducing the number and volume of bubbles.
A ZWICK ROELL/Z050 test machine with a
BTC-EXACLBI.001 clip-on biaxial extensometer was
used for the experimental tests.
The specimens were subjected to cyclic tensile tests.
During each cycle, unloading was started every time the
elongation between the extensometer arms (the gage
area) exceeded a multiple of l = 0.02 mm. During the
corresponding cycle, loading started when the tensile
force decreased below the multiple of 30 % of the force
at the start of the unloading. All the tests were performed
at room temperature of 22 °C. Figure 2 shows a typi-
cally cracked specimen with a rupture.
3 MATERIAL MODEL
The free-energy-function method was used for the
formulation of the model.4,5 In this paper, the free-energy
function  is proposed in the form of the sum of elastic,
plastic and damage parts:


 	  
	

= − + +
⎡
⎣
⎢
⎤
⎦
⎥∫∫
1 1
2
1 2
00
E D bE( )( ) )P Pd d
P
(1)
where  is the density of the material; E is Young’s
modulus; D is the damage factor; E is the elastic strain;
 P is the equivalent plastic strain; R is the hardening-
associated quantity;  is the damage variable; and B is
the damage-associated quantity. Generalized forces are
described in more detail below. The stress is calculated
in Equation (2) as:

 


= = −
∂
∂ E
EE D( )1 (2)
where stress 
 is considered as the nominal stress acting
in the whole area of the cross-section of a specimen.
Stress D is the stress acting in the undamaged area of
the cross-section of a specimen:4



 
 

D E
E=
−
=
−
=
( ) ( )1 1D D
E
∂
∂
(3)
The damage-associated energy, in this case identi-
cally equal to the strain energy density of the undamaged
material, is
Y
D
E= − =


∂
∂
1
2
2( )E (4)
The hardening-associated quantity, which defines the
plastic material behavior, is:
R( ) 


P
P=
∂
∂
(5)
And the damage-associated quantity, which defines
the damage material behavior, is:
B( ) 


=
∂
∂
(6)
The model is considered as a nonlinear material and
the infinitesimal-strain theory is used in Equation (7):


=
l
l
(7)
where  is the total strain; l is the elongation of the
gage area and l is the original length of the gage area.
The elasto-plastic decomposition of the axial strain  is
also considered:6,7
  = +E P (8)
where P is the plastic strain. The function of plasticity
decides whether the plastic flow will be present. It is
proposed in the following form:
F R RP r
y
r
P= − = − + ≤
 
 
 ( ( ))0 0 (9)
where 
r = 
 D , which implies the assumption that the
plastic flow is present only in the undamaged part of the
cross-section and that the plastic flow will occur in
tension as well as in compression; 
y is the yield stress
and R0 is the initial yield stress. The function of damage
is proposed in the following form:
F Y B BD r= − + ≤( ( ))0 0 (10)
where Yr ≡ Y, therefore, as in the case of plasticity, no
distinction between the tensile and compressive damage
is considered; B0 is the energy necessary for the damage
initiation.
Rates of quantities are expressed below. Time deri-
vatives ("·") are replaced by time step changes ("")
because all the analyses were performed as quasi-static.
The change in the plastic strain P is:
 


 
P P
P
P sign= =
∂
∂
F
( ) (11)
where P is the plastic multiplier. The change in the
equivalent plastic strain  P is:
  P P
P
P= − =
∂
∂
F
R
(12)
The change in the damage factor D is:
  D
F
Y
= − =D
D
D∂
∂
(13)
where D is the damage multiplier. The change in the
damage variable is:
T. KROUPA et al.: ONE-DIMENSIONAL ELASTO-PLASTIC MATERIAL MODEL WITH DAMAGE ...
214 Materiali in tehnologije / Materials and technology 51 (2017) 2, 213–218
MATERIALI IN TEHNOLOGIJE/MATERIALS AND TECHNOLOGY (1967–2017) – 50 LET/50 YEARS
Figure 2: Cracked specimen
Slika 2: Prelomljen vzorec
  = − =D
D
D∂
∂
F
Y
(14)
If inequalities (9) and (10) are fulfilled, the calcu-
lation of responses is simple. Otherwise, the changes of
the unknown variables (P,  P , D, , P, D)
must be calculated. The following section shows the
solution of this problem in the presented model.
4 CALCULATION OF RESPONSES
The calculation of the unknown variables is per-
formed in several steps.
1. The so-called trial state is used. Firstly, the change
in the total strain is considered only as a change in the
elastic strain, which leads to:
  trial
E E= +t (15)
where  trial
E is the trial value of the elastic strain;  t
E is
the elastic strain at the end of the previous time step and
 is the change in the total strain in the present time
step. Trial values of the other quantities are considered
as the values at time t. The values of the plasticity and
damage functions are calculated using (9) and (10). If
the inequalities in the mentioned equations are fulfilled,
no plastic or damage quantities are changed and the trial
values of strain and stress calculated using (2) are con-
sidered as the values in time t+t and the calculations
for the time step are finished. Otherwise the return
mapping algorithm described in step 2 is used.
2. Using relations:
    P trial
E E (16)
and
    P P trial
P (17)
and the fact obvious from relations (13) and (14) that D
=  = D, the set of six equations from (9) to (14) can be
transformed into the set of two equations with two
unknowns (E and D).
The first equation has the following form:
E R K n  E R
P R− =0 0( ) (18)
where KR and nR are the shape parameters of the hard-
ening function. The equivalent plastic strain  P is cal-
culated as:
    P trial
P
trial
P E Esign= + −( ) ( )E (19)
where  trial
P is the trial value of the equivalent plastic
strain. The term in the shape of the classical power law
used in Equations (9) and (18):

 
y nR K= −0 R
P R( ) (20)
is the hardening function.
The second equation is:
1
2
02 0E B B D( ) ( ( ))
E − + = (21)
where term (B0 + B(D)) defines the damage behavior of
the material.
Equations (18) and (21) can be solved separately. The
first equation (18) needs to be solved numerically. In the
present paper, E is obtained with the golden-section
method in the interval   trial
P
trial
P− . Other variables
can be calculated from relations (16) and (19).
The solution of the second equation (21) can be
written in the explicit form. It is not necessary to know
the exact shape of function B = B() = B(D) because it is
only necessary to define inverse function B–1 = B–1(Y) =
B–1(E). Then
D B E E= −⎛⎝
⎜ ⎞
⎠
⎟−1 2 21
2
1
2
( ) ( ) E 0
E (22)
where 0
E is the initial strain for the damage evolution
and B0 is assumed in the following form:
B E0
21
2
= ( )
E (23)
The exact shape of D is, in this paper, proposed in the
following form:
D D K
K
n
=
−⎡
⎣⎢
⎤
⎦⎥
+
⎧
⎨
⎩
+ −
−
u B
E
0
E
f
E
B
E
( ) ( )
( )
( )
( ) (
 

 
2 2
2
2
1
0
E
f
E
)
( )
2
2
⎡
⎣⎢
⎤
⎦⎥
⎫
⎬
⎭
m
(24)
with
n
a
=
⎛
⎝
⎜
⎞
⎠
⎟


r
f
1
(25)
and
m
a
=
⎛
⎝
⎜
⎞
⎠
⎟


f
r
2
(26)
where Du is the ultimate damage factor, for which the
total failure of the material occurs; KB and r are the
parameters of the damage-factor dependency defining
the position of the inflection point of the damage curve;
 f
E is the strain at which failure occurs and a1 and a2 are
the shape parameters.
3. Once the plastic strain, damage factor and all the
other variables including the stress are known, the next
time step can be solved.
The model includes 11 material parameters sum-
marized in Table 1.
The calculation of the whole analysis of the cyclic
tensile test with almost 3000 time steps is done in less
than 1 second on the machine with Intel® Core
I7-4790, a four-core CPU, running at the frequency of
3.60 GHz using Windows 7 operating system with an
SSD disk.
T. KROUPA et al.: ONE-DIMENSIONAL ELASTO-PLASTIC MATERIAL MODEL WITH DAMAGE ...
Materiali in tehnologije / Materials and technology 51 (2017) 2, 213–218 215
MATERIALI IN TEHNOLOGIJE/MATERIALS AND TECHNOLOGY (1967–2017) – 50 LET/50 YEARS
Table 1: Material parameters
Tabela 1: Parametri materiala
Parameter Units Meaning
E Pa Elastic modulus
R0 Pa Initial yield stress
Kr Pa Shape parameter of the hardeningcurve
nR – Shape parameter of the hardeningcurve
Du – Maximum value of damage that thematerial can absorb before failure
KB – Coefficient defining the position ofthe inflection point
0
E – Initial elastic strain for damage
f
E – Strain at which failure occurs
r
E
– Strain defining the position of theinflection point
a1 – Shape parameter of damage curve
a2 – Shape parameter of damage curve
5 IDENTIFICATION
Raw experimental data are pre-processed for the
further identification. Important data points such as the
maxima and minima of the loading cycles and cross
points are found (Figure 3). Averaged values of these
points serve as the target for the further identification.
In order to identify the material parameters, the diffe-
rence between numerical and experimental results is
minimized. Function r defining this difference is the
so-called total residual. It is the sum of five particular
residuals and it is minimized in the identification
process.
The first residual is the difference between tensile
curves rL:
r
N
F l F l
F
i i
i
N
L
L
E F
E
L
=
−⎛
⎝
⎜
⎞
⎠
⎟
=
∑1
1
2
( ) ( )
max
Δ Δ
(27)
where superscripts E and N stand for the experimental
and numerical analysis, respectively; NL is the number
of points, in which the curves are compared (the number
of Pmax points); li is the elongation of the specimen
related to the point Pmax and the denominator Fmax
E is the
maximum force in the target curve used as the weight
coefficient.
According to Figure 3, the target for tensile curves in
the experiment is considered as the connection of Pmax
points. Pure epoxy shows a negligible difference bet-
ween points Pmax and Pint unlike the other materials such
as textile composites, rubber, etc. In a numerical analy-
sis, if a model does not consider viscoelasticity (which is
true of the presented model) these two points coincide.
Figure 4 shows the idea of a comparison of the tensile
curves and the calculation of residual rL.
The second residual rK is the difference between the
tangents of the unload-load cycles of the experimental
data kE and numerical data kN:
r
N
k l k l
k
i i
i
N
K
K
E N
E
K
=
−⎛
⎝
⎜
⎞
⎠
⎟
=
∑1
1
2
( ) ( )
max
Δ Δ
(28)
where NK is the number of the points, in which the
curves are compared; Δli is the elongation of the speci-
men related to the point Pmax and the denominator is the
maximum value of the tangent in the experiment used as
the weight coefficient. The tangents k are the first
derivatives of the straight line connecting points Pmax
and Pmin (Figure 3).
The third residual rmF is the difference between the
values of the maximum forces:
r
F
FmF
N
E= −
⎛
⎝
⎜
⎞
⎠
⎟1 max
max
(30)
where Fmax
N is the maximum force reached in the
numerical analysis and Fmax
E is the value from the target
function.
The fourth residual rmD is the difference between the
values of the elongations related to the points of maximal
forces:
r
l
l
F
F
mD
N
E= −
⎛
⎝
⎜
⎞
⎠
⎟1
Δ
Δ
max
max
(31)
T. KROUPA et al.: ONE-DIMENSIONAL ELASTO-PLASTIC MATERIAL MODEL WITH DAMAGE ...
216 Materiali in tehnologije / Materials and technology 51 (2017) 2, 213–218
MATERIALI IN TEHNOLOGIJE/MATERIALS AND TECHNOLOGY (1967–2017) – 50 LET/50 YEARS
Figure 4: Comparison of the tensile curves
Slika 4: Primerjava nateznih krivulj
Figure 3: One unload-load cycle and important points
Slika 3: En cikel razbremenjeno – obremenjeno in pomembne to~ke
where Δl Fmax
N and Δl Fmax
E are the elongations related to
the points where force-elongation curves reached their
maxima of the force in the numerical analysis and in the
experiment.
Finally, the fifth residual rE is the difference between
the values of the force measured after brittle crack (zero
values) and the forces calculated in the numerical
analysis, subsequent to the point where the curve reaches
its maximum:
r
N
F l
F
i
i
N
E
E
N
E
E
=
⎛
⎝
⎜
⎞
⎠
⎟
=
∑1
1
2
( )
max
Δ
(32)
This particular residual is used to induce the model to
get to the state where the brittle crack occurs.
Once the particular residuals are known, the total
residual can be calculated as:
r r r r r r= + + + +L T mF mD E (33)
Optimization software OptiSLang 3.2.3. was used.
Gradient and simple design improvement algorithms
were used.8
6 RESULTS
Figure 5 shows the force-elongation dependency.
Note that the deviation of the maximum forces in the
experiments is not negligible.
Considering the deviation of the experimental results,
the identified dependency shows sufficient agreement;
nevertheless, it is slightly more curved than the experi-
mental target.
Figure 6 shows the comparison of the results for the
calculated tangents of the unload-load cycles. The first
cycle was ignored because of the influence of the
inaccuracies of the experimental results for low loading
forces. The comparison of the curves is performed within
the interval where all the measured specimens provide
the necessary data. The first specimen reached its
strength limit at an elongation of about 0.25 mm. The
justification for this choice is as follows: after reaching
the mentioned elongation limit, significant jumps in the
dependence of the averaged experimental results occur
because of the disappearing data from each cracked
specimen.
Figure 7 shows the identified hardening function.
Dependency of the damage factor on the elastic strain
is shown in Figure 8. The modulus ED = E(1 – D)
occurs right before the brittle fracture and equals
3.55 GPa and the damage factor is D = 0.47 (see the gray
circle in Figure 8).
The identified material parameters including those
calculated directly from the experiment are given in
Table 2.
T. KROUPA et al.: ONE-DIMENSIONAL ELASTO-PLASTIC MATERIAL MODEL WITH DAMAGE ...
Materiali in tehnologije / Materials and technology 51 (2017) 2, 213–218 217
Figure 7: Hardening function
Slika 7: Funkcija utrjevanja
Figure 5: Force-elongation dependency. Dots represent points Pmax .
Grey bold line is the target experiment; black bold line is the nume-
rical analysis.
Slika 5: Odvisnost sile od raztezka. To~ke predstavljajo Pmax. Siva,
oja~ana linija je cilj preizkusa, ~rna oja~ana linija je numeri~na ana-
liza.
Figure 6: Tangent of unload-load cycle vs. elongation. Grey bold line
is the target experiment; black bold line is the numerical analysis.
Dots represent tangents k for each cycle in each experiment.
Slika 6: Tangenta na cikel obremenjeno – razbremenjeno, raztezek.
Siva oja~ana krivulja je cilj raziskave, ~rna debela krivulja je nume-
Table 2: Identified parameters. Parameters with * were calculated
directly from the experimental results.
Tabela 2: Ugotovljeni parametric. Parametri ozna~eni z *, so bili
izra~unani neposredno iz eksperimentalnih rezultatov
Parameter Value Unit
E 6.64 GPa
R0 0.00 MPa
KR 149.19 GPa
nR 1.10 -
Du 0.76 -
KB 0.62 -
0
E 0.00 %
r
E 0.46 %
f
E 1.17 %
a1 2.48 -
a2 6.95 -
v* 0.29 -

S* 41.67 MPa
7 CONCLUSION
A relatively simple and effective algorithm for the
calculation of the response of a material, showing elastic,
plastic and damage behavior is proposed. The plastic and
damage-flow problems are solved separately. A set of
equations is transformed into a single nonlinear equation
in the case of the plastic-flow problem and the equation
is solved with the golden-section algorithm, which is
extremely effective. The other set of equations that need
to be solved in order to calculate the damage factor is
transformed into one explicit formula for calculating the
damage factor.
An appropriate identification process is proposed and
its effectivity is shown in the case of the cyclic tensile
tests of pure epoxy resin. The function that is minimized
during the identification process is the sum of five inde-
pendent residuals, each quantifying a difference between
experimental and numerical results. These residuals
express the differences between tensile-elongation
dependencies and the values of the tangents of the
unload-load cycles, the distance between the point of the
maximum force reached in force and elongation direc-
tions and the value of the force after the brittle fracture
of the specimen is minimized. The model and the iden-
tification process are robust and time effective.
Acknowledgement
This publication was supported by project LO1506 of
the Czech Ministry of Education, Youth and Sports.
8 REFERENCES
1 T. Kroupa, V. La{, R. Zem~ík: Improved nonlinear stress–strain
relation for carbon–epoxy composites and identification of material
parameters, Journal of Composite Materials, 45 (2011) 9,
1045–1057, doi:10.1177/0021998310380285
2 R. Zem~ík, R. Kottner, V. La{, T. Plundrich, Identification of ma-
terial properties of quasi-unidirectional carbon-epoxy composite
using modal analysis, Mater. Tehnol., 43 (2009) 5, 257–260
3 S. Ogihara, K. L. Reifsnider: Characterization of Nonlinear Behavior
in Woven Composite Laminates, Applied Composite Materials, 9
(2001), 249–263, doi:10.1023/A:1016069220255
4 S. Murakami, Continuum damage mechanics, A continuum mecha-
nics approach to the analysis of damage and fracture, Springer, 2012
5 E. A. de Souza Neto, D. Peri}, D. R. J. Owen, Computational Me-
thods for Plasticity: Theory and Applications, Wiley, 2008
6 V. La{, R. Zem~ík, Progressive Damage of Unidirectional Composite
Panels, Journal of Composite Materials, 42 (2008) 1, 25–44,
doi:10.1177/0021998307086187
7 K. J. Bathe, Finite element procedures, Prentice Hall, Upper Saddle
River, New Jersey 07458, USA, 2007
8 J. S. Arora, Introduction to Optimum Design, Elsevier, 2004
218 Materiali in tehnologije / Materials and technology 51 (2017) 2, 213–218
T. KROUPA et al.: ONE-DIMENSIONAL ELASTO-PLASTIC MATERIAL MODEL WITH DAMAGE ...
MATERIALI IN TEHNOLOGIJE/MATERIALS AND TECHNOLOGY (1967–2017) – 50 LET/50 YEARS
Figure 8: Damage function. Grey circle shows the initiation of brittle
fracture and the final failure of the material.
Slika 8: Funkcija po{kodb. Sivi krogi ka`ejo pri~etke krhkega loma in
kon~no poru{itev materiala