Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

54, 3, pp. 731-742, Warsaw 2016
DOI: 10.15632/jtam-pl.54.3.731

FINITE STRAIN FORMULATION OF ELASTO-PLASTICITY WITHOUT

YIELD SURFACE: THEORY, PARAMETER IDENTIFICATION AND

APPLICATIONS

Cyprian Suchocki, Paweł Skoczylas

Warsaw University of Technology, Department of Mechanics and Armament Technology, Warsaw, Poland

e-mail: c.suchocki@imik.wip.pw.edu.pl

In this study, a finite strain formulation of elasto-plasticitywithout the yield surface is ana-
lyzed. It is found that the presentedmodel allows for a good description of the elasto-plastic
response for many metallic materials. Furthermore, it is noted that the model involving a
limited number of material parameters is able to capture such effects as indefinite elastic
limit, strain hardening, strain softening and the Bauschinger effect. An algorithm for the
evaluation of constants is described and the determined values are given for NChigh carbon
tool steel, copper, brass, Inconel alloy, tungsten heavy alloy and aluminum.

Keywords: constitutive equation, plasticity, metals, material constants

1. Introduction

In continuum mechanics, two main concepts of the plasticity theory can be distinguished, i.e.
flow theories and theories of plasticitywithout theyield surface, including endochronic plasticity.
The origins of constitutive theories based on the flow rules date back to the mid-19th century.
The fundamentals of plasticity theories not employing the yield surface were formulated much
later (e.g. Ilyushin, 1963; Pipkin and Rivlin, 1965).

Thefirst steps associatedwith replacing theflow rule andyield criterion-based formulation of
elasto-plasticity with a new one are claimed to Ilyushin (1963) who first introduced the concept
of arc length in the strain space. This approach was later extended by Pipkin andRivlin (1965)
who developed a novel, general theory of plasticity for the case of finite strains and rotations,
however, without proposing a specific form of the constitutive equation.

The arc length calculated in nine-dimensional finite strain space is the basic concept for
the new formulation of elasto-plasticity. This scalar value is utilized to build a functional of
strain history analogous to the history functional used in the theory of viscoelasticity. Pipkin
and Rivlin (1965) considered the most general form of the strain history functional, which in
practice can be assumed in the form of a hereditary integral, where, instead of time, the arc
length in the strain space is assumed as the variable of integration. For a constitutive theory
constructed in such a manner, the arc length plays the role of an intrinsic time analogously to
the time in the theory of viscoelasticity.

Pipkin and Rivlin (1965) postulated that the arc length should be calculated utilizing a
Pythagorean metric. A different approach was introduced by Valanis (1971a,b) who proposed
using a non-Pythagorean metric for calculation of the arc length in the infinitesimal strain
space and additional multiplication by a scaling function. The plasticity theories which employ
a scaling function in their definition of the intrinsic time are called endochronic (cf. Valanis,
1971b). Such concepts were strongly criticized (e.g. Rivlin, 1981) as nonphysical. Some of the
imperfections of the original endochronic theory were eliminated in a later constitutive model



732 C. Suchocki, P. Skoczylas

(cf. Valanis, 1980). However, this was achieved at the price of splitting the strain tensor into
elastic and plastic components and introducing a yield criterion.

A number of models exist which transfer some concepts of the endochronic formulation to
the flow theory. In particular, the flow rule for kinematic hardening may utilize a hereditary
integral to govern the evolution of the back stress (e.g. Backhaus, 1976). This framework was
further extended to take into account such effects as the memory of maximal prestrain and the
amplitude-dependent hardening in the case of cyclic loadings (cf. Mróz and Rodzik, 1996). The
effective infinitesimal plastic strain is used as the intrinsic time in models of this type.

In the present work, a finite strain formulation of plasticity without the yield surface, pre-
viously introduced by Suchocki (2015), is analyzed upon its ability to capture the inelastic
response of metals. In order to clearly lay out basic theory assumptions, the infinitesimal strain
formulation is described first. Subsequently, the novel finite strain generalization is presented
and examined. It is found that the proposed constitutive model, which involves a very limited
number of material parameters, is not only able to accurately fit the stress-strain data but also
to predict such behavior as indefinite elastic limit, strain hardening, strain softening and the
Bauschinger effect. A numerical algorithm used for the evaluation of the material constants is
described. The model parameters are determined for a variety of metallic materials. Finally, a
parameric study is conducted.

2. Small strain formulation

Below, the basic notions of rate-independent plasticity are discussed in the case of infinitesimal
deformations, cf. Rivlin (1981). The following assumptions have been adopted:

• Only isothermal deformation processes are considered.
• Isotropic materials with elasto-plastic properties are analyzed.
• The volumetric response is purely elastic (Fig. 1a).
• The total shear stress is assumed as a sum of the elastic stress and N rate-independent
stress-like internal state variables (Fig. 1b).

• The internal state variables depend on shear deformation only and are governed by diffe-
rential evolution equations.

Fig. 1. Mechanical scheme of the constitutive model; (a) volumetric response, (b) shear response

The total Cauchy stress tensor σ is a sum of the volumetric stress tensor p1 and the stress
deviator s, i.e.

σ= p1+s p=Kǫ s= s0+
N∑

k=1

h̃k s0 =2Ge (2.1)



Finite strain formulation of elasto-plasticity without yield surface... 733

where

ǫ= trε e= ε−
1

3
(trε)1 (2.2)

whereas h̃k (k = 1,2, . . . ,N) are tensorial internal state variables while K and G are the bulk
and shear moduli, respectively. The internal state variable evolution equation is a modification
of the differential equation following from theMaxwell mechanistic model, i.e.

˙̃
hk+

1

D̃kM(|ė|)
h̃k = γkṡ0 M(|ė|)= (tr ė2)−

1

2 k=1,2, . . . ,N (2.3)

where γk is a dimensionless parameter and the product D̃kM(|ė|) is a rate-dependent relaxation
time. The particular choice of M(|ė|) given by (2.3)2 allows one to eliminate time from the
evolution equation. After inserting Eqs (2.1)4 and (2.3)2 into Eq. (2.3)1 andmultiplying it by a
time differential, an incremental equation is obtained

dh̃k+

√
trde2

D̃k
h̃k =2γkG de k=1,2, . . . ,N (2.4)

which states that the evolution of h̃k (k=1,2, . . . ,N) depends solely on the deformation path
in the strain space and is rate-independent.
An intrinsic time ζ is introduced, cf. Valanis (1971a,b), i.e.

ζ̇ =
1

M(|ė|)
ζ(τ)=

τ∫

−∞

dτ ′

M(τ ′)
ζ(t)= z (2.5)

After substituting Eq. (2.3)2 into Eq. (2.5)2

ζ(τ)=

τ∫

−∞

√
de(τ ′) ·de(τ ′) (2.6)

It can be seen from (2.6) that the intrinsic time has an interpretation of the arc length in the
deviatoric strain space. Integration of Eq. (2.3)1 yields

h̃k(t)=

z(t)∫

0

γke
−
z−ζ

D̃k
∂s0(ζ)

∂ζ
dζ k=1,2, . . . ,N (2.7)

which is reminiscent of a similar equation present in both linear and nonlinear viscoelasticity.

3. Finite strain formulation

All the basic assumptions listed in the previous Section apply to the finite strain model aswell.
The total stressS is given as a sumof the volumetric stressSvol0 , elastic isochoric stressS

iso
0 and

the stress-like tensorial state variables H̃k (k=1,2, . . . ,N), i.e.

S=Svol0 +S
iso
0 +

N∑

k=1

H̃k S
vol
0 = J

∂U

∂J
C
−1

S
iso
0 =2J

−
2

3

[∂W
∂C

]
(3.1)

where, C and C are the total and isochoric right Cauchy-Green (C-G) tensors, respectively,
J = detF is the determinant of the deformation gradient tensor F, whereas U = U(J) and



734 C. Suchocki, P. Skoczylas

W = W(Ī1, Ī2) are the elastic stored energy functions with Ī1 and Ī2 being scalar invariants

of C while DEV[•] = [•]− 1
3
([•] ·C)C−1. The evolution of internal state variables is governed

by differential equations

˙̃
Hk+

1

D̃kM(|Ċ|)
H̃k = γkṠ

iso
0 k=1,2, . . . ,N (3.2)

with

M(|Ċ|)= J̄−
1

2

2 J̄2 = trĊ
2

(3.3)

This specific choice of M(|Ċ|) ensures that the evolution of H̃k (k = 1,2, . . . ,N) is rate-
independent. Integration of (3.2) yields

H̃k(t)=

z(t)∫

0

γke
−
z−ζ

D̃k
∂Siso0 (ζ)

∂ζ
dζ k=1,2, . . . ,N (3.4)

where the following equalities are used

ζ̇ =
1

M(|Ċ|)
ζ(τ)=

τ∫

−∞

√
dC(τ ′) ·dC(τ ′) ζ(t)= z (3.5)

Equation (3.4) is a special case of the general, functional stress-strain relation formulated by
Pipkin andRivlin (1965), and the intrinsic time introduced in Eqs (3.5) is amodification of the
fictitious timemeasure proposed therin.

4. One-dimensional tension-compression

A uniaxial tension along theX1-direction (Fig. 2b) is defined by the following set of equations

x1 =λ1X1 x2 =λX2 x3 =λX3 (4.1)

where the principal stretches λ1 > 1 and λ < 1 for tension, whereas λ1 < 1 and λ > 1 for
compression processes. Applying the boundary condition of the stress-free lateral surface, the
scalar equations can be written as follows

S11 =S
vol
0(11)+S

iso
0(11)+

N∑

k=1

H̃k(11) 0=S
vol
0(22)+S

iso
0(22)+

N∑

k=1

H̃k(22) (4.2)

Is is assumed thatW(C)=W(Ī1), which yields

S
vol
0 =JpC

−1 p=
∂U

∂J
S
iso
0 =2J

−
2

3
∂W

∂Ī1

(
1− 1
3
I1C

−1
)

(4.3)

In the considered case I1 = λ
2
1 +2λ

2, which leads to the following formulas for the isochoric
elastic stress components

Siso0(11) =
4

3
J−

2

3
∂W

∂Ī1

[
1−
( λ
λ1

)2]
Siso0(22) =−

1

2

(λ1
λ

)2
Siso0(11) (4.4)



Finite strain formulation of elasto-plasticity without yield surface... 735

By utilizing Eq. (4.3)1 in Eqs (4.2), one obtains

S11 = Jp
1

λ21
+Siso0(11)+

N∑

k=1

H̃k(11) 0= Jp
1

λ2
+Siso0(22)+

N∑

k=1

H̃k(22) (4.5)

After eliminating the Jp factor, and using (4.4)2, the following form of the process equation is
found

S11 =
3

2
Siso0(11)+

N∑

k=1

[
H̃k(11)−

( λ
λ1

)2
H̃k(22)

]
(4.6)

In this study, two specific forms of W are considered, i.e. neo-Hooke and the one proposed
by Knowles

W(Ī1)=
µ

2
(Ī1−3) W(Ī1)=

µ

2b

{[
1+
b

κ
(Ī1−3)

]κ
−1
}

(4.7)

with µ, b and κ being the elasticity constants. After inserting (4.7)1 into (4.4)1, the following
formula is found

Siso0(11) =
2

3
J−

2

3µ
[
1−
( λ
λ1

)2]
(4.8)

whereas substituting (4.7)2 into (4.4)1 yields

Siso0(11) =
2

3
J−

2

3µ
[
1+
b

κ
(Ī1−3)

]κ−1[
1−
( λ
λ1

)2]
(4.9)

In the case of material incompressibility J = 1 and λ = λ
−
1

2

1 , whereas Ī1 = I1 = λ1 +2λ
−
1

2

1 .
Thus, Eqs (4.8) and (4.9) take respectively forms

Siso0(11) =
2

3
µ
(
1− 1
λ31

)
Siso0(11) =

2

3
µ
[
1+
b

κ
(I1−3)

]κ−1(
1− 1
λ31

)
(4.10)

By applying the incompressibility condition to Eq. (4.6), the one-dimensional process equation
is found, i.e.

S11 =
3

2
Siso0(11)+

N∑

k=1

(
H̃k(11)− H̃k(22)

1

λ31

)
(4.11)

whereas the engineering and Cuachy stresses are calculated according to the formulas

T11 =λ1S11 σ11 =λ1T11 (4.12)

Equations (4.10)-(4.12) along with Eq. (4.4)2 and two scalar equations following from Eq.
(3.2) form up a system of nonlinear algebraic and differential equations which can be solved
numerically. It is assumed that the excitation is kinematic with a constant strain rate T , thus,
at the time increment n+1

tn+1 = tn+∆t λ1(tn+1)=1+Ttn+1

In the case of incompressibility J̄2 = J2 = trĊ. Furthermore, an auxiliary quantity is defined as

J̃(tn+1)= J2(tn+1)∆t
2



736 C. Suchocki, P. Skoczylas

Material parameter optimization algorithm

Input:p0, (T̃11)j,wj (j=1,2, . . . ,K),N, T̃ , imax, i=1

1. Calculate the inital value of the weighted error function

F0(p0)=
K∑

j=1

wj

[
(T11(p

0))j − (T̃11)j
]2

2. Assemble the initial matrix of parameter values

M
0
5×(1+2N) =




µ1 γ11 · · · γN1 D̃11 · · · D̃N1
µ2 γ12 · · · γN2 D̃12 · · · D̃N2
...
...
...

...
...

...
...

µ5 γ15 · · · γN5 D̃15 · · · D̃N5




3. Find a parameter combination minimizing the error function

p
i =
[
µopt,γ

opt
1 , . . . ,γ

opt
N
,D̃
opt
1 , . . . ,D̃

opt
N

]T

4. Check the stop criterion, i.e.

∣∣∣
Fi(pi)−Fi−1(pi−1)
Fi−1(pi−1)

∣∣∣< T̃

5. Assemble the matrix of parameter values for the next iteration

(a) Calculate parameter value increments (k=1,2, . . . ,N)

∆p=
[
∆µ,∆γ1, . . . ,∆γN,∆D̃1, . . . ,∆D̃N

]T

∆µ=
|µopt|
5i

∆γk =
|γopt
k
|

5i
∆D̃k =

|D̃opt
k
|

5i

(b) ComputeMi
5×(1+2N)

components (k=1,2, . . . ,5, l=1,2, . . . ,1+2N)

[
M
i
5×(1+2N)

]

kl
=Mikl = p

i
l +
k−3
2
∆pl

6. Go to step3 if i < imax; otherwise display optimalmaterial parameter values

which is determined by the following formula

J̃(tn+1)= 2(T∆t)
2
{2
3

[1
3
(λ41(tn+1)+

2

λ21(tn+1)

)
−2
](
λ31(tn+1)−

1

λ31(tn+1)

)2

+2λ21(tn+1)+
1

λ41(tn+1)

}

According to Eqs (4.10)1 and (4.4)2, the isochoric elastic stresses are given as

Siso0(11)(tn+1)=
2

3
µ
(
1−

1

λ31(tn+1)

)
Siso0(22)(tn+1)=−

1

2
λ31(tn+1)S

iso
0(11)(tn+1)



Finite strain formulation of elasto-plasticity without yield surface... 737

The total second P-K stress is calculated using Eq. (4.11), assumingN =1, i.e.

S11(tn+1)=
3

2
Siso0(11)(tn+1)+ H̃1(11)(tn+1)−

1

λ31(tn+1)
H̃1(22)(tn+1)

whereas the internal state variable components are calculated from the following recurrence-
update formulas

H̃1(11)(tn+1)=

(
1− 1

2D̃1J̃
−

1
2(tn+1)

)
H̃1(11)(tn)+γ1

(
Siso
0(11)
(tn+1)−Siso0(11)(tn)

)

1+ 1
2D̃1J̃

−

1
2(tn+1)

H̃1(22)(tn+1)=

(
1− 1

2D̃1J̃
−

1
2(tn+1)

)
H̃1(22)(tn)+γ1

(
Siso
0(22)
(tn+1)−Siso0(22)(tn)

)

(
1+ 1

2D̃1J̃
−

1
2(tn+1)

)

which have been developed utilizing the central difference method.
Thegiven above set of discretized process equations in the formderived for an incompressible

material is utilized within the parameter evaluation algorithm for computation of theoretical
stress values. The assumption of material incompressibility has been adopted for the sake of
simplicity. The relative error in estimation of the stretch ratio in the lateral direction that
results from this assumption does not exceed 2% for the considered class of materials.

5. Experimetal setup

Anumber of experiments have been performed in order to obtain the data which allowed one to
asses the model’s capability of capturing the mechanical response of metals. The Instron 1115
material testing machine was used for the experiments (Fig. 2a).

Fig. 2. (a) Experimental set-up. (b) Deformation process and coordinate systems. (c) Dogbone
specimens

The conducted tests included one-dimensional tension (Fig. 2b) of copper, Inconel No. 1
(supersaturated), Inconel No. 2 (supersaturated and aged), brass, NC tool steel, aluminum and



738 C. Suchocki, P. Skoczylas

tungsten heavy alloy (THA). In each case, unloading was performed shortly before the failure.
The specimens (Fig. 2c) were machined according to the Polish standard PN-91 H-04310. Each
tension test was repeated for at least three times in order to make sure that the registered
stress-strain curveswere reproducible. The specimenswere tested at room temperature of 20◦C.

6. Curve fitting

A pattern search algorithm has been used for the evaluation of the model parameters. It is a
modification of the scheme proposed by Zalewski (2009). In each iteration, amatrixMi

5×(1+2N)
containing several sets of material constants is assembled. Every parameter set is subsequently
used to compute theoretical stress values. A combination of constant values pi+1, which results
in the best approximation of the stress-stretch curve, is allowed to proceed to the next iteration
i+1. The performed calculations have been listed in an outline, which assumed making use
of neo-Hooke stored-energy function. The curve fitting results are presented in Figs. 3 and 4,
whereas the evaluated values of the material parameters are gathered in Table 1.

Table 1.Determined values of material parameters

Material µ [MPa] b [–] κ [–] γ1 [–] D̃1 [–]

copper 367.33 22.41 0.68 132.80 0.0053

Inconel No. 1 1114.68 - - 83.43 0.0034

Inconel No. 2 1152.58 - - 74.12 0.013

brass 681.68 - - 35.67 0.013

steel NC 6 1065.19 - - 47.06 0.0081

aluminum 546.70 - - 30.72 0.0098

THA 1226.94 - - 88.48 0.0060

In each case, a single internal state variable has been assumed. In the case of copper, the
version of the constitutive equation employing theKnowles stored-energy functionhas beenused
in order to reproduce the strain-softening phenomenon (Fig. 3a). For all the remaining metals,
the neo-Hooke function has been utilized.

An excellent agreement between the theoretical and experimental curves has been achieved
for THA (Fig. 4a). In some cases, fitting monotonic stress-stretch data leads to obtaining a
better approximation as it can be seen in the case of brass (Fig. 4b).

7. Parametric studies

Several simulations have been carried out in order to asses the model behavior, its sensitivity
and the influence of the experimental data used for material parameter evaluation.

In Fig. 5a, the results of loading-unloading simulations are presented which have been per-
formed for five different values of the maximum stretch. The material parameter set has been
evaluated from the brass monotonic loading data (model 1). In Fig. 5b, the results of perfor-
ming the same simulations are shown where the set of constants has been determined from the
brass hysteresis data (model 2). In Fig. 5c, the hysteresis loops obtained for the two different
material parameter sets are compared. It can be seen that utilizing themonotonic loading data
for the determination of material constants instead of using the hysteresis data, does not pro-
duce any negative consequences. What is more, an excellent approximation is achieved in the
initial part of the stress-stretch curve (Fig. 4b). In Fig. 5d, the result of one-dimensional, cyclic
loading-unloading simulation is shown (model 1).



Finite strain formulation of elasto-plasticity without yield surface... 739

Fig. 3. Fitting of constitutive model predictions to experimental data: (a) copper, (b) Inconel No. 1,
(c) Inconel No. 2, (d) brass, (e) NC steel, (f) aluminum

Fig. 4. Curve fitting results: (a) THA hysteresis data, (b) brass monotonic tension data



740 C. Suchocki, P. Skoczylas

Fig. 5. Results for the brass deformation process simulation; (a) loading-unloading simulations for
model 1, (b) loading-unloading simulations for model 2, (c) comparison of simulation results,
(d) hysteresis loop for brass, (e) cyclic unloading-loading with small stretch amplitude, (f) cyclic

unloading-loading with large stretch amplitude

The response to cyclic loading-unloading stretch histories with a prestrain has been inve-
stigated (model 1). The stretch was monotonically increased to λmax1 = 1.01. Subsequently,
unloading was performed to λmin1 = 1.0097 which was followed by loading until λ

max
1 [–] was

reached again. The unloading-loading cycle was repeated for 40 times. It can be seen in Fig. 5e
that the unloading slope is larger than in the case of the following reloading, so that the hyste-
resis loop does not form and a stress decrease is observed. This defect is known for the group
of endochronic models and was discussed by Rivlin (1981). The same simulation was repeated
for a larger stretch amplitude, i.e. λmin1 = 1.009 and λ

max
1 = 1.02 which resulted in formation

of a stress-stretch loop (Fig. 5f). It is concluded that for cyclic loadings with a prestrain, the
model’s applicability is limited to the case of large stress/stretch amplitudes.



Finite strain formulation of elasto-plasticity without yield surface... 741

In Fig. 6 the results of two loading-unloading simulations, i.e. in tension and in compression,
are plotted. Thematerial constant values determined for brass (model 1) are utilized. It can be
seen that thematerial responses in tension and compression are almost the same.What ismore,
the Bauschinger effect can be observed.

Fig. 6. Simulation of the Bauschinger effect in brass for tension and compression

Fig. 7. Parametric study results; (a) hysteresis loops for different values of parameter D̃1, (b) hysteresis
loops for different values of the coefficient γ1, (c) the case of a rigid-perfectly plastic solid

The influence of the constant values has been investigated. In Fig. 7a, the uniaxial loading-
-unloading stress-stretch curves obtained for different values of the parameter D̃1 are shown. In
Fig. 7b, the uniaxial loading-unloading curves generated for different values of the coefficient γ1



742 C. Suchocki, P. Skoczylas

are plotted. Some adjustments of the D̃1 value allow one to obtain the response of a rigid-
-perfectly plastic solid (Fig. 7c).

8. Conclusions

In this study, a finite strain formulation of plasticity without a yield surface has been presented
and analyzed. The material parameter values are identified for a number of metallic materials.
A good agreement between the theoretical predictions and the experimental data is achieved.
It is found that the presentedmodel can capture such behavior as indefinite elastic limit, strain
hardening, strain softening and the Bauschinger effect. The model limitations are discussed as
well. The algorithm formaterial constant evaluation has been presented and a parametric study
conducted.

Acknowledgement

This work was supported by project INNOTECH-K2/IN2/27/182101/NCBR/13.

References

1. Backhaus G., 1976, Fließspannungen und Fließbedingung bei zyklischen Verformungen, ZAMM,
56, 8, 337-348

2. Ilyushin A.A., 1963, Plasticity, Foundations of the General Mathematical Theory (in Russian),
Academy of Sciences of the USSR,Moscow

3. Mróz Z., Rodzik P., 1996, On multisurface and integral description of anisotropic hardening
evolution in metals,European Journal of Mechanics – A/Solids, 15, 1, 1-28

4. Pipkin A.C., Rivlin R.S., 1965,Mechanics of rate-independentmaterials, ZAMP, 16, 313-327

5. Rivlin R.S., 1981, Some comments on the endochronic theory of plasticity, International Journal
of Solids and Structures, 17, 231-248

6. Suchocki C., 2015, An internal-state-variable based viscoelastic-plastic model for polymers, Jo-
urnal of Theoretical and Applied Mechanics, 53, 593-604

7. Valanis K.C., 1971a, A theory of viscoplasticity without a yield surface, Part I: General theory,
Archives of Mechanics, 23, 517-534

8. Valanis K.C., 1971b, A theory of viscoplasticity without a yield surface, Part II: Application to
mechanical behavior of metals,Archives of Mechanics, 23, 535-551

9. Valanis K.C., 1980, Fundamental consequences of a new intrinsic time measure, plasticity as a
limit of the endochronic theory,Archives of Mechanics, 32, 171-191

10. Zalewski R., 2009, Numerical method of Chaboche’s model parameters identification basing on
special granular structures experimental data,Modelling in Engineering, 38, 2, 309-319, (in Polish)

Manuscript received July 22, 2015; accepted for print October 28, 2015