Formulation and implementation of the VPJC material model in EUROPLEXUS.
Date: 2016-07-20
Institutions: Institute for the Protection and Security of the Citizen
Subjects: EU research policy, accident prevention, materials technology, research report, resistance of materials, simulation, software

Formulation and Implementation of 
the VPJC Material Model in 
EUROPLEXUS 

Vegard Aune 

Folco Casadei 

Georgios Valsamos 

Tore Borvik 
2016 

EUR 27982 EN

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
                                                                    
 
 
 
 
 
Formulation and Implementation of 
the VPJC Material Model in 
EUROPLEXUS 

 
 
 
 
This publication is a Technical report by the Joint Research Centre, the European Commission’s in-house science 
service. It aims to provide evidence-based scientific support to the European policy-making process. The scientific 
output expressed does not imply a policy position of the European Commission. Neither the European 
Commission nor any person acting on behalf of the Commission is responsible for the use which might be made 
of this publication. 

JRC Science Hub 
https://ec.europa.eu/jrc 

JRC102336 

EUR 27982 EN 

ISBN 978-92-79-59746-6  

ISSN 1831-9424  

doi:10.2788/609529  

© European Union, 2016 

Reproduction is authorised provided the source is acknowledged. 

Printed in Italy 

All images © European Union 2016 

How to cite: Authors; title; EUR; doi 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Formulation and Implementation of the VPJC Material Model in
EUROPLEXUS

Vegard Aunea,b,

∗, Folco Casadeic, Georgios Valsamosc, Tore Børvika,b

aCentre for Advanced Structural Analysis (CASA), NTNU, Norwegian University of Science and Technology, NO-7491
Trondheim, Norway
bStructural Impact Laboratory (SIMLab), Department of Structural Engineering, NTNU, NO-7491 Trondheim, Norway
cEuropean Commission, Joint Research Centre (JRC), Institute for the Protection and Security of the Citizen (IPSC),
European Laboratory for Structural Assessment (ELSA), 21027 Ispra, Italy

Abstract

A constitutive model of elastic-thermoviscoplasticity and ductile failure is implemented in EUROPLEXUS
(EPX) using a fully vectorized version of the forward Euler integration algorithm and a special case of the
backward Euler algorithm. The material model is based on the constitutive model originally proposed by
Johnson and Cook, and the ductile failure criterion proposed by Cockcroft and Latham. Constitutive equa-
tions are formulated in terms of linear elasticity, von Mises yield criterion, associated ﬂow rule, nonlinear
isotropic strain hardening, strain-rate hardening, temperature softening due to adiabatic heating, and duc-
tile failure. The performance of the implementation is veriﬁed and validated in EPX. First, a veriﬁcation
is performed using single element simulations in uniaxial tension and simple shear. Then, the model is
validated against experimental data from quasi-static tension tests and airblast experiments. The model is
applicable for one-, two- and three-dimensional stress analysis, and coupled with element deletion options
available in EPX.

Keywords:
Elastic-thermoviscoplasticity, Isotropic hardening, von Mises yield criterion, Adiabatic thermal softening,
Cutting plane return map algorithm, Radial return map algorithm, VPJC, EUROPLEXUS

1. Introduction

In metal plasticity, extreme loading conditions in fast transient dynamics often involve large strains, high
strain rates, temperature softening and ductile failure. A widely used design tool for this class of problems
is the non-linear ﬁnite element (FE) method. In the FE method, the formulation and numerical solution
of non-linear problems in continuum mechanics relies on the weak form of the momentum balance equation
(also known as the principal of virtual work [1]). The integral form of the principal of virtual work is well
suited for direct application to the FE method, and by spatial discretization the solution is integrated in
time.
In the case of fast transient dynamics, it is often used an explicit time integration based on the
central diﬀerence method. Iterations for FE analyses involving non-linear material behaviour can in general
be divided into two levels, hereby denoted the global and local level. The global level involves the explicit
establishment of global equilibrium between internal stresses and external loads, and the local level updates
the corresponding stress state in each FE integration point (for a given strain increment) in terms of the
governing constitutive equations. Constitutive equations (also known as material models) are mathematical
descriptions of material behaviour which gives the stress as a function of the deformation history of the body

∗Corresponding author
Email address: vegard.aune@ntnu.no (Vegard Aune )

Technical report - Formulation and implementation of the VPJC material model

June 24, 2016

[1]. This implies that the local integration of the constitutive equations controls the accuracy and stability
of the global equilibrium, and that the overall solution in the FE model is considerably inﬂuenced by the
accuracy, robustness and eﬀectiveness of the integration algorithm.

A commonly used scheme for the local stress integration is the predictor-corrector method (often called
return mapping). The basis for this return mapping is two successive steps, i.e., the predictor step and the
corrector step. The predictor step is used to estimate a trial stress state, while the corrector step applies a
ﬂow rule by using a return mapping scheme to ensure the consistency condition. That is, ensuring that the
stress state is on the yield surface. The present study uses the explicit cutting plane method proposed by
Ortiz and Simo [2] and the implicit closest point projection method proposed by Simo and Taylor [3]. During
these approaches an elastic predictor uses the total strain increment and the linear elastic relation to obtain
a trial stress state. Then, if the elastic trial stress is outside the plastic domain, the elastically predicted
stresses are corrected to return to the updated yield surface by iterations of the plastic corrector. The basis
for the cutting plane method is to utilize the known stress state at the last converged time increment to
determine the normal to the yield surface, while the implicit closest point projection method utilizes the
cylindrical shape of the yield surface and the deviatoric stress space to scale the trial state to ﬁt a suitably
updated yield surface by iterations of the plastic corrector [4, 5].

This work considers a well-known constitutive relation proposed by Johnson and Cook [6] and von Mises
plasticity is used to update the stresses. A model including elastic-thermoviscoplasticity and ductile failure
is implemented as a new material called VPJC in EUROPLEXUS (EPX) [7]. The implementation uses a
fully vectorized version of the forward Euler integration algorithm and a special case of the backward Euler
algorithm. Large deformations are accounted for by using a hypoelastic-plastic formulation in the constitu-
tive equations and co-rotational formulations in EPX (i.e., the Jaumann stress rates [7]). This implies that
the constitutive equations are deﬁned in a co-rotational coordinate system in which the basis system rotates
with the material. The subroutine interface has been implemented using the Cauchy stress components, and
all stress and strain quantities are deﬁned in terms of the rate of deformation tensor. The model is applicable
for one- (1D), two- (2D) and three-dimensional (3D) stress analyses (i.e., for bar, shell, solid, axisymmetric
and plane strain elements). The model is coupled with the element deletion options available in EPX by
introducing a state variable controlling the element erosion. That is, as the failure variable reaches its critical
value DC the element is removed. The model is ﬁnally veriﬁed and validated through numerical simulations
in EPX. This is ﬁrst done by a veriﬁcation using single element simulations in uniaxial tension and simple
shear. Then, the model is validated against experimental data from quasi-static tension tests and airblast
experiments.

2. Constitutive Equations

There are in general two categories within the plasticity theory, i.e. the mathematical theories and the
physical theories [8]. Mathematical theories are purely phenomenological and aim to represent experimental
observations through mathematical relations. Most of the common material models today are phenomeno-
logical. Physical theories try to describe the underlying mechanisms, and quantify the plastic deformation,
at a microscopical level. A well-known and acknowledged physical theory is the crystal plasticity theory.
This section considers the phenomenological-based formulation of the constitutive equations for the compu-
tational model used in this study, starting with the basic theory for inﬁnitesimal strains and then expanding
the theory to also include large deformations. The proposed computational model includes linear elasticity,
large strains, high strain rates, temperature softening and ductile failure. Assuming negligible elastic strains
compared to the plastic strains, an hypoelastic formulation may be adopted through an additive decom-
position of the rate of deformation tensor. Since the main objective is the numerical implementation of
thermo-viscoplasticity, the constitutive equations are restricted to linear isotropic elastic behavior. Thermo-
dynamic coupling is not included, however, the temperature evolution is assumed to be adiabatic. Ductile
failure is considered by using a simple one-parameter failure criterion suggested by Cockcroft and Latham

2

[9] which is based on plastic work per unit volume.

From a computational point of view, the plastic response is an incremental process which implies that the
constitutive equations need to be formulated on rate form. Starting by assuming inﬁnitesimal deformations,
the strain rate tensor ˙ε may be decomposed by an additive split

where ˙εe is the elastic and ˙εp is the plastic strain rate tensor, respectively. Using the generalized Hooke’s
law, the rate of the stress tensor is deﬁned as

˙ε = ˙εe + ˙εp

(1)

where C is the fourth order elasticity tensor given by

˙σ = C : ˙εe

where I is the second order identity tensor, and λel and µel are the Lam`e constants given as

C = λelI

⊗

I + 2µelI

The elasticity tensor C is assumed to be constant, E is the Young’s modulus and ν the Poisson’s ratio.

λel =

νE
(1 + ν)(1

2ν)

and

µel =

E
2(1 + ν)

−

In general, assuming only isotropic hardening, the yield criterion is deﬁned as

f (σ, R) = ϕ(σ)

σy(σ0, R) = σeq −

−

(σ0 + R)

(2)

(3)

(4)

(5)

where ϕ(σ) = σeq is the equivalent stress, σ0 is the initial yield stress and R = R(p) is the isotropic hardening
variable depending on the equivalent (or accumulated) plastic strain p. The evolution of the yield surface is
therefore governed by the elastic domain, described by the ﬂow stress σy. The function f is assumed to be a
continuous function of the stress tensor σ, and results in negative values when the material is in the elastic
domain. In the case of rate-independent plasticity and isotropic hardening, it is assumed that f (σ) > 0 is
inadmissible and, consequently, that plastic deformation occurs when f (σ) = 0.

Furthermore, ϕ(σ) is assumed to be an isotropic and positive homogeneous function of order one of the
stress tensor. This implies that

where a could be any positive scalar. Then, Euler’s theorem for homogeneous functions gives [10]

ϕ(aσ) = aϕ(σ)

(6)

∂ϕ
∂σ
The plastic deformation process is governed by an important physical restriction, i.e., that work has to be
done to the body all the time for the deformation to continue. This restriction is expressed as

= ϕ

σ :

(7)

where Dp expresses the plastic dissipation per unit volume, which is assumed to be dissipated as heat.

Dp = σ ˙εp

0

≥

(8)

The restriction in Eq.(8) is fulﬁlled by deﬁning the plastic strain rate tensor ˙εp in a way that ensures
non-negative dissipation. This is done by the plastic ﬂow rule, which in the general case is expressed as

˙εp = ˙λ

∂g
∂σ

3

(9)

where g = g(σ)

≥

0 is the plastic potential function and ˙λ

≥

0 is a scalar denoted the plastic multiplier.

Combining Eqs.(8) and (9) result in the following restriction on the plastic potential function

∂g
∂σ ≥
where it is assumed that g is a positive homogeneous function of order one. Then, by using the Euler’s
theorem for homogeneous functions, it is seen that this assumption fulﬁlls the dissipation inequality in
Eq.(8) as

(10)

σ :

0

By further assuming that the plastic potential function g is deﬁned by the yield function f , Eq.(9) is expressed
as

Dp = σ ˙εp = ˙λσ :

∂g
∂σ

= ˙λg

0

≥

(11)

˙εp = ˙λ

∂f
∂σ

(12)

and it is seen that the plastic potential function is associated with the yield function f . Eq.(12) is therefore
called the associated ﬂow rule, and ensures non-negative dissipation through Eqs.(5), (7) and (8) as

Dp = σ ˙εp = ˙λσ :

∂f
∂σ

= ˙λσ :

∂ϕ
∂σ

= ˙λϕ = ˙λσeq ≥

0

(13)

This implies, that in the case of the associated ﬂow rule, the equivalent plastic strain rate equals the plastic
parameter

where the equivalent plastic strain p is deﬁned as the energy conjugate variable to the equivalent stress σeq.
The plastic dissipation could then be expressed in terms of the equivalent stress and the plastic strain rate
as

˙p = ˙λ

(14)

Dp = σ ˙εp = σeq ˙p
Hence, assuming the associated ﬂow rule implies that the incremental plastic strain vector dεp = ˙εpdt is
f at σ, and therefore directed along the outward normal n of
parallel to the gradient of the yield surface
the yield surface (see Figure 1). Moreover, the magnitude of ˙εp is given by the plastic multiplier ˙λ.

(15)

∇

An important remark to this argumentation is that the shape of the yield surface also determines the direc-
tion of the plastic ﬂow, as well as the stress state at which yielding initiates.

The associated ﬂow rule may also be derived from a stronger assumption than the dissipation inequality in
Eq.(8), i.e., the principle of maximum plastic dissipation expressed as

where σ is the stress state at which the plastic dissipation occurs, and σ∗ is an arbitrary stress state inside
or on the yield surface f (σ∗)
0. As for the dissipation inequality this principle states the normality of
plastic ﬂow, and in addition it is seen that the yield surface must be convex [11]. Thus, the use of an associ-
ated ﬂow rule ensures that the plastic strain increments dεp are vectors perpendicular to the potential surface.

≤

σ∗) ˙εp

(σ

−

0

≥

(16)

To distinguish between elastic loading/unloading and plastic loading, the Kuhn-Tucker conditions are used
to determine when plastic ﬂow may occur. These conditions determines the evolution of the yield surface
and are given as [11]

4

σ2

f = ∂f

∂σ e

∇

σ

f (σ) = 0

σ1

dεp = dλ ∂f

∂σ = dλn

Figure 1: Geometric representation of the associated ﬂow rule for an elliptic yield surface in plane stress
space. The orthogonal unit vectors e are pointing in the coordinate directions.

f

0,

˙λ

0,

˙λf = 0

(17)

≥
where the ﬁrst of these conditions indicates that the stress states must lie on or within the yield surface, the
second indicates that the plastic multiplier is non-negative, and the last condition ensures that the stresses
lie on the yield surface during plastic loading. From the Kuhn-Tucker conditions, the consistency condition
can be stated as [10]

≤

The consistency condition ensures that the stress must be located on the yield surface during continuous
plastic ﬂow, and is, in the theory of rate-independent plasticity, used to determined the plastic multiplier ˙λ.

˙λ ˙f = 0

(18)

To account for viscoplasticity and strain-rate sensitivity, the theory of plasticity has to be extended. The
main diﬀerence is that in the viscoplastic domain the stress is allowed to move outside the yield surface, by
assuming that only elastic strains are developed while the stress is inside or on the yield surface (f
0) [11].
This means that plastic strains are assumed to evolve only when the stresses move outside the yield surface
(Figure 2b). By using an additive constitutive relation, this implies that the yield function f is equal to the
viscous stress σv( ˙p), i.e.,

≤

f (σ) = σv( ˙p) > 0

fv(σ) = ϕ(σ)

σy(σ0, R)

σv( ˙p) = 0

(19)

⇐⇒

−

−

From Eq.(19) it is seen that the viscous stress controls the strain-rate sensitivity of the material (Figure
2). Although the additive formulation is illustrative for the viscoplastic material behavior, a multiplicative
formulation may be a preferred choice in some applications. In general, the multiplicative formulation reads

fv(σ) = ϕ(σ)

σy(σ0, R)υp( ˙p) = 0

−

⇐⇒
where υp( ˙p) is determined by the constitutive relation. It should be noted that in both Eq.(19) and Eq.(20)
the viscoplastic constraint fv(σ) = 0 is introduced, which allows for an expansion of the elastic domain due
to the rate sensitivity of the material. The rate sensitivity of the material is controlled by the constitutive
relation for the plastic multiplier, generally expressed as

f (σ) = σv(p, ˙p) = (υp −

1)σy > 0

(20)

(cid:26)

˙p = ˙λ =

0
φ(σ, R)

for
for

f
0
≤
f > 0

5

(21)

where φ(σ, R) depends on the state of the material. From Eq.(21) it is seen that the consistency condition in
Eq.(18) no longer holds in the case of a viscoplastic formulation, and that the plastic ﬂow depends explicitly
on time during the evolution of ˙λ.

σeq

σ0

σv( ˙p, p)

R(p)

π-plane

σ3

σ

f = 0

f > 0

(a) Equivalent stress vs equivalent plastic strain

(b) The deviatoric plane (π-plane)

σ1

p

σ2

Figure 2: Viscoplastic eﬀects due to strain-rate sensitivity.

To account for large displacement and large strain, the constitutive model is deﬁned in a co-rotational
coordinate system in which the basis system rotates with the material [7]. Typically, there are two diﬀerent
approaches of dealing with ﬁnite deformations, i.e., the hypoelastic-plastic formulations and hyperelastic-
plastic formulations.
In contrast to the hyperelastic formulation of ﬁnite plasticity models, hypoelastic
theories do not rely on the assumption of the existence of a strain energy potential to model the reversible
behaviour. Moreover, in the particular case of metal plasticity, the elastic strains are negligible compared to
the plastic strains and it is reasonable to adopt an hypoelastic formulation for large deformations. This is
analogous to the additive decomposition of the strain-rate tensor ˙ε in Eq. (1) in the theory for inﬁnitesimal
deformations, and the rate-of-deformation tensor D is decomposed additively as [1]

where De and Dp are the elastic and plastic part, respectively. The rate-of-deformation tensor D is also
called velocity strain and is a rate measure of the deformation and stretching of the elements.

D = De + Dp

(22)

The starting point for hypoelastic models is the formulation of the constitutive equations in terms of objective
stress rates [5], in which the hypoelastic relation in terms of the objective stress rate is expressed as

where σ is the Cauchy stress tensor, λel and µel are the Lam`e constants given in Eq.(4), and ˆσ∇ is an
objective (or frame-invariant) stress rate given as

ˆσ∇ = C : De = λel (trDe) I + 2µelDe

(23)

−
where ˙σ is the time derivative of the Cauchy stress tensor σ and A is an appropriate angular velocity tensor.
By combining Eqs.(23) and (24) the objective stress rate is obtained as

ˆσ∇ = ˙σ

Aσ + σA

(24)

The choice of an appropriate angular velocity tensor A is typically done by the respective FE codes, where
EPX [7] uses the Jaumann rate of the Cauchy stress ˆσ∇

J given as

˙σ = Aσ

−

σA + λel (trDe) I + 2µelDe

(25)

6

ˆσ∇

J = ˙σ

−

W σ + σW ,

A = W =

(cid:0)L

1
2

−

LT (cid:1)

(26)

where W is the spin tensor indicating the rate of rotation of an element and L is the velocity gradient (see
e.g. [1] for more details). Another frequently used objective stress is obtained by using the Green-Naghdi
stress rate ˆσ∇

G [1] expressed as

ˆσ∇

G = ˙σ

−

˙RRT σ + σ ˙RRT = ˙σ

Ωσ + σΩ,

−

A = Ω = ˙RRT =

R ˙RT

−

(27)

where R is the rotation tensor representing the rigid body motion in the polar decomposition F = RU (see
Figure 3) and Ω is the skew-symmetric rate of the rotation tensor.

Figure 3: Polar decomposition of a typical homogeneous deformation. The Cartesion coordinate system
with coordinate axes (x1, x2) refers to the global reference frame, while (ˆx1, ˆx2) represents the local (or
co-rotated) coordinate system.

From Eqs. (26) and (27) it is seen that the Green-Naghdi rate diﬀers from the Jaumann rate only in using
a diﬀerent measure of the rotation of the material. Despite the diﬀerence in Eqs. (26) and (27), it is worth
noting that when using polar decomposition F = RU the spin tensor W can be expressed as [12]

W = ˙RRT

(cid:16)

R

1
2

−

˙UU−

1 + U−

T ˙U

(cid:17)

RT

(28)

where U is the symmetric right stretch tensor and it is used that ˙F = LF . Thus, for pure rigid motion
( ˙U = 0) the Green-Naghdi and Jaumann stress rates give the same result since

W = ˙RRT = Ω

(29)

7

x1x2Rˆx1ˆx2ˆx2Uˆx1ˆx2ˆx1Fwhereas in for instance shear one should be aware of the diﬀerent formulation and how this may inﬂuence
the results of the FE simulation [1].

As concerns the plastic deformation rate tensor Dp this is deﬁned by means of the associated ﬂow rule as

Dp = ˙λ

∂f
∂σ

(30)

2.1. The von Mises yield criterion

The von Mises yield criterion is an isotropic and pressure-insensitive criterion commonly used for metals and
alloys [11]. The pressure insensitivity in metals and alloys is due to the fact that the plastic deformation
to a large extent takes place by plastic slip, which is a shear driven deformation mode [10]. Geometrically,
the yield criterion f (σ) = 0 deﬁnes a surface in the stress space which is deﬁned by the components of the
stress tensor σ (i.e., the vector space in Figure 4a). The elastic domain is given by the stress space inside
this yield surface, while the plastic domain is the surface itself.

yield surface

σ3

e = 1√
3

[1

1

1]

Hydrostatic axis
σ1 = σ2 = σ3

π-plane

yield curve

r
e·
r

Rv
P

σ2

σ3

Rv

π-plane

(Deviatoric plane)
σ1 + σ2 + σ3 = 0

σ1

σ1

(cid:113) 2

3 σy

σ2

(a) Principal stress space

(b) The deviatoric plane (π-plane)

Figure 4: The von Mises yield criterion.

From Figure 4a it is seen that the von Mises yield surface is a circular cylinder of radius Rv =
3 σy.
The stress state on the yield surface could then be represented in the principal stress space as the vector
r = [σ1 σ2 σ3], which has to satisfy the following equation [12]

(cid:113) 2

·
By splitting the stress in deviatoric and hydrostatic parts

−

·

·

r

r

(e

r) (e

r) = R2
v

(cid:48)

σ

= σ

σH I

−

and

σH =

1
3

tr(σ) =

σ1 + σ2 + σ3
3

(31)

(32)

(cid:48)

where σ
expressed as

is the deviatoric stress and σH is the hydrostatic stress, the von Mises yield criterion could be

R2

v = σ : σ

(cid:18) 1
√3

−

(cid:19)2

tr(σ)

= (σ

(cid:48)

ij + σH δij)(σ

ij + σH δij)

(cid:48)

3σ2

H = σ

ijσ

(cid:48)

−

(cid:48)

ij = σ

(cid:48)

(cid:48)

: σ

=

2
3

σ2
y

(33)

Eq.(33) implies that isotropic (or hydrostatic) states of stress do not result in yielding. Thus, the von Mises
yield criterion is purely deviatoric and additional hydrostatic stress only translates the stress state along the
hydrostatic axis in Figure 4a. Furthermore, the result in Eq.(33) could be presented in terms of the yield
function f from Eq.(5), i.e.,

8

f (σ) = σeq −

σy =

(cid:114)

3
2

σ(cid:48) : σ(cid:48)

σy =

−

(cid:112)

3J2 −

σy

where the von Mises equivalent stress σeq, in terms of the stress components, reads

(cid:113)

σeq =

11 + σ2
σ2

22 + σ2

33 −

σ11σ22 −

σ22σ33 −

σ33σ11 + 3(σ2

12 + σ2

23 + σ2

31)

and the associated ﬂow rule in Eq.(30) is then given as

Dp = ˙λ

∂f
∂σ

= ˙λ

3
2

(cid:48)

σ
σeq

(34)

(35)

(36)

The rate of the accumulated plastic strain ˙p, in the particular case of associated ﬂow and von Mises isotropic
behaviour, is deﬁned as

˙p =

(cid:114)

2
3

Dp : Dp

(37)

2.2. The modiﬁed Johnson-Cook constitutive model

As discussed earlier, engineering structures subjected to fast transient dynamics may experience large strains,
high strain rates and temperature softening. A widely used constitutive relation allowing for these factors is
the constitutive model of Johnson and Cook [6], which is empirical and dates back to 1983. The constitutive
relation has a multiplicative formulation and is given by

where A, B, C, n and m are material constants determined from material tests, and

σy = (A + Bpn) (1 + Cln( ˙p∗)) (1

T ∗

m)

−

˙p∗ =

˙p
˙p0

and

T ∗ =

T
−
Tm −

Tr
Tr

represents the dimensionless strain rate ˙p∗ and the non-dimensional homologous temperature T ∗, respec-
tively. The parameter ˙p0 is a user-deﬁned reference strain rate, while Tr and Tm are deﬁned as the room
and melting temperature, respectively. The ﬁrst term in Eq.(38) represents the strain hardening, the sec-
ond term describes the inﬂuence of the strain rate, and the last term represents the temperature softening
behaviour. The second and third term shift the hardening curve up or down, depending on the strain rate
and the temperature conditions. It is seen from Eq.(38) that the ﬂow stress σy increases with increasing
plastic strain rate (which is a necessary characteristic in viscoplasticity), while it decreases with increasing
temperature. Thus, an increase in the plastic strain rate expands the yield surface whereas an increase in
temperature contracts the yield surface.

To avoid non-physical softening eﬀects when ˙p∗ < 1, the strain-rate-sensitivity term in the Johnson-Cook
model can be adjusted [13]. This leads to the modiﬁed Johnson-Cook model, given by [14]

σy = (A + Bpn) (1 + ˙p∗)C (1

T ∗

m)

−

(40)

This is illustrated by a comparison of Eqs.(38) and (40) in Figure 5, where the ordinate axis refers to the
m)]. Note that, since the formulations in Eqs.(38) and (40) are
normalized stress σy/ [(A + Bpn) (1
slightly diﬀerent, the calibration of the viscoplastic material parameter C may be somewhat diﬀerent in the
two versions of the Johnson-Cook model.

T ∗

−

The ﬁrst term of the Johnson-Cook model could lead to numerical instabilities when n < 1 due to the initial
value of p, i.e., p equal zero before any plasticity occur. This results in an inﬁnite value of the hardening
1 when plasticity occurs. This could be avoided by deﬁning the initial plastic strain
modulus HR = Bnpn

−

9

(38)

(39)

Figure 5: Original vs modiﬁed Johnson-Cook model, when CJC = 0.014, CM JC = 0.013, and ˙p0 = 5e−

4.

as an inﬁnite small value. However, it is desirable to formulate the material model as general as possible.
Therefore, the ﬁrst term of Eq.(40), representing the strain hardening, is replaced by another hardening
rule than that proposed by Johnson and Cook. By replacing the ﬁrst term of Eq.(40) with a two terms
saturation type of hardening rule proposed by Voce [15], the constitutive model could be expressed as

σy = (A + R(p)) (1 + ˙p∗)C (1

T ∗

m) = (A +

−

2
(cid:88)

k=1

Qk(1

exp (

−

−

Ckp))) (1 + ˙p∗)C (1

T ∗

m)

−

(41)

where the constant A is the initial yield strength of the material, and Qk and Ck (k = 1, 2) deﬁne the
strain hardening. As before, the rate sensitivity is governed by the constant C, while m deﬁnes the thermal
softening behaviour. The hardening modulus HR is found from the rate of change in R(p) given as

˙R =

dR
dt

=

dR
dp

dp
dt

= HR

dp
dt

=

⇒

HR =

dR
dp

=

2
(cid:88)

k=1

QkCke(
−

Ckp)

(42)

It is seen from Eq.(42) that no numerical instabilities are expected for the hardening modulus HR when the
plastic strain equals zero.

Finally, evaluating yielding in the viscoplastic domain and combining Eqs.(34) and (41), the strain-rate
sensitivity of the material is controlled by

˙p = ˙λ =






(cid:20)(cid:16)

˙p0

0

ϕ(σ)
σy(σ0,R,T )

(cid:17) 1

C

(cid:21)

1

−

for

for

f

0

≤
f > 0

(43)

2.3. Adiabatic thermal softening

Assuming that the plastic work is dissipated as heat implies that large plastic deformation results in an
increase in temperature. The result of this increase in temperature is a decrease in ﬂow stress. This is
already introduced in the last term of Eq.(41), which accounts for the thermal softening of the ﬂow stress
at elevated temperatures. However, the evolution of the temperature remains to be established. Since the
plastic response of the material is of a very short time duration in fast transient dynamics, the heat transfer
is modeled by assuming adiabatic conditions. This implies that there is no heat transfer into or out of the
system during the plastic straining. Using the thermal energy balance per unit volume given by [10]

10

10−510−31001031050.80.911.11.2Dimensionlessstrainrate˙p∗Normalizedstress  Original Johnson−CookModified Johnson−Cookwhere ρ is the material density, Cp is the speciﬁc heat capacity of the solid material and χ is the Taylor-
Quinney coeﬃcient. The temperature increase due to adiabatic heating can be computed as

Dp = χσ : Dp = χσeq ˙p = ρCp ˙T

(44)

˙T =

Dp
ρCp

(45)

From Eq.(44) it is seen that the Taylor-Quinney coeﬃcient χ represents the fraction of the plastic power that
is converted into heat. The remaining fraction 1
χ is assumed to be stored in the material due to structural
rearrangements (e.g. elastic ”ﬁelds” around dislocations). A value of 0.9 is commonly found in the literature
for metals, however, a conservative choice would be 1.0 which implies maximum thermal softening. In reality
the χ-value will not be constant, however, advanced experimental techniques are necessary to determine the
evolution of this parameter. Moreover, one should always keep in mind that, while strain hardening (i.e.,
work hardening due to creation of dislocations) implies an increase of strength locally in the material and
distributes the plasticity, softening results in localization of plasticity.

−

2.4. Ductile failure

A material model is not complete without some form of material degradation or failure. The degradation
or failure in a material is usually given in terms of a damage parameter, and failure occurs through damage
evolution. Ductile fracture arises from the nucleation, growth and coalescence of microscopic voids that
initiate at inclusions and second phase particles. The voids around particles grow when the material is
subjected to plastic straining and hydrostatic tension, and fracture occurs when the growing voids reach a
critical size, relative to their spacing, resulting in a local plastic instability between the voids [16]. This work
considers only ductile failure which is deﬁned as the ﬁrst sign of fracture, i.e., the nucleation (or initiation) of
voids. This means that the fracture criterion is uncoupled from the constitutive equations, and is therefore
denoted a failure criterion in the following. The reader is referred to [5, 17] if damage coupling is important
in the constitutive model.

Cockcroft and Latham [9] suggested a simple failure criterion which is based on plastic work per unit volume.
They reasoned that the failure criterion needed to be based on some combination of stress and strain, and
that damage accumulates during plastic straining. To account for hydrostatic tension they based the criterion
on the magnitude of the major principal stress σ1, i.e.,

D =

W
Wc

=

1
Wc

(cid:90) p

σ1(cid:105)

0 (cid:104)

dp

(46)

where Wc is the failure parameter which can be found by integrating the major principal stress in a uniaxial
tension test during the entire equivalent plastic strain path until the plastic strain at failure pf . The expres-
is equivalent to the function max(0, σ1), which implies that only positive values of the maximum
sion
principal stress σ1 (i.e., tensile stress) contribute to the damage evolution. Material failure takes place when
the damage parameter D reaches unity. This failure criterion is therefore easy to implement in a FE code
and it requires only one uniaxial tension test for calibration.

σ1(cid:105)
(cid:104)

The failure criterion by Cockcroft and Latham has proven to give equally good results as more complex
criteria [18, 19, 20], and it is also shown that the criterion indirectly accounts for both deviatoric and
hydrostatic stress states [21]. This is seen by using an alternative expression for the major principal stress
[22], i.e.,

σ1 = σH +

3 + µσ
(cid:112)

3 + µ2
σ

3

(cid:32)

σeq =

σ∗ +

(cid:33)

σeq

3 + µσ
(cid:112)
3

3 + µ2
σ

(47)

where the stress triaxiality σ∗ and the Lode parameter µσ is deﬁned by

11

Remember that σ1 ≥

σ2 ≥

σ∗ =

,

µσ =

σH
σeq

2σ2 −

σ3

σ1 −
σ1

σ3 −
σ3 are the ordered principal stresses. Eq.(46) can then be expressed as

D =

1
Wc

(cid:42)

(cid:90) p

σ∗ +

0

3 + µσ
(cid:112)

3 + µ2
σ

3

(cid:43)

σeqdp

Since material failure occurs when D = Dc = 1.0 and p = pf , Eq.(49) can be written as

1 =

1
Wc

(cid:42)

(cid:90) pf

σ∗ +

0

3 + µσ
(cid:112)

3 + µ2
σ

3

(cid:43)

σeqdp

(48)

(49)

(50)

Assuming that the stress triaxiality σ∗ and Lode parameter µσ remain constant throughout the entire loading
history (i.e., proportional loading) and using the Voce hardening rule from Eq.(41) to express the equivalent
stress σeq through the yield condition in Eq.(34), Eq.(50) now reads

σ∗ =

Apf + (cid:80)2

Wc
k=1(Qkpf + Qk
Ck

[exp (

−

Ckpf )

−

1]) −

3 + µσ
(cid:112)
3

3 + µ2
σ

(51)

The inﬂuence of stress triaxiality σ∗ and Lode parameter µσ on the plastic failure strain pf can then be
illustrated as shown in Figure 6, where the material parameters by Holmen et al. [21] are used to generate
the respective curves. Various loading scenarios are now represented by diﬀerent combinations of the stress
triaxiality and the Lode parameter. Some typical examples are (µσ, σ∗) = (+1, 1/3) which implies uniaxial
σ2 = σ3), (µσ, σ∗) = (0, 0) corresponding to pure shear (2σ2 = σ1 + σ3) and (µσ, σ∗) =
tension (σ1 ≥
(
−
−

1/3) representing uniaxial compression (σ1 = σ2 ≥

σ3).

1,

(a) Failure surface

(b) Selected Lode parameters

Figure 6: Representation of the plastic failure strain pf in the (µσ, σ∗) space when using the Cockcroft-
Latham failure criterion and assuming proportional loading. Note that the necessary material parameters
to be used in Eq. (51) are taken from Table 1.

From Eq.(46) it is seen that, to make the material model as general as possible, it is preferable to compute
the principal stresses directly in the subroutine. This is done by following the procedure proposed by [12]
which uses the principal deviatoric stresses σ
i to ﬁnd the principal stresses σi. The principal deviatoric
stresses are determined from the eigenvalue problem for the deviatoric stress σ

in Eq. (32), i.e.,

(cid:48)

(cid:48)

12

01234567−11/20Lodeparameter,µσ1/212/3Stresstriaxiality,σ∗1/30Plasticfailurestrain,pf−1/300.250.50.7511.251.51.752−2/3−1/301/32/3Stresstriaxiality,σ∗Plasticfailurestrain,pf  µσ=−1µσ=0µσ=+1where det(
called the characteristic equation, given as

) is the determinant of the matrix argument (

·

). The result of this operation is a cubic equation

·

det(σ

(cid:48)

I

−

(cid:48)

σ

) = 0

(52)

−
where the principal invariants are given respectively as

−

(cid:48)

)3

(σ

(cid:48)

I

1(σ

(cid:48)

)2

(cid:48)

J2σ

−

J3 = 0

(53)

(cid:48)

(cid:48)

I

1
2
and the J-symbols for the principal invariants are introduced since these are commonly used in the literature
for metal plasticity. The general solution is then given as

ii = 0, J2 =

ij, J3 = I

3 = det(σ

1 = tr(σ

(cid:48) 2) =

) = σ

ijσ
σ

I
−

tr(σ

(54)

2 =

1
2

)

(cid:48)

(cid:48)

(cid:48)

(cid:48)

(cid:48)

(cid:48)

(cid:114)

(cid:48)

σ

1 = 2

J2
3

cos

θ
3

,

(cid:114)

(cid:48)

2 = 2
σ

J2
3

cos

(cid:18) θ

3 −

(cid:19)

2π
3

(cid:48)

1,

σ

≤

(cid:114)

(cid:48)

3 = 2
σ

J2
3

cos

(cid:19)

(cid:18) θ
3

+

2π
3

(cid:48)
2

σ

≤

(55)

where the angle θ is given as

J3
J 3
2 /27
and the Lode angle is recognized as θL = θ/3. Finally, the principal stresses σi of the stress tensor σ are
determined by

cos θ =

(cid:112)
2

(56)

≤

≤

π

0

θ

,

(cid:48)

(cid:48)

σi = σ

i + σH = σ

i +

1
3

tr(σ)

(57)

where the major principal stress σ1 is used in the Cockcroft-Latham (CL) failure criterion in Eq.(46).

3. Two-dimensional stress analysis

This section contains useful relationships when dealing with the particular case of plane stress. Plane stress
is deﬁned as the state of stress in which the normal stress σ33 and the shear stresses σ13,23 are assumed to
be zero [11]. From a practical point of view this assumption is typically used in the analysis of bodies in
which one of the dimensions (i.e., the thickness) is much smaller than the others and is subjected to loads
that result in dominant stresses perpendicular to the thickness direction [5].

3.1. Two-dimensional theory of elasticity

In the particular case of plane stress, the linear elastic stress-strain compliance relationship for an isotropic
material can be written as a system of equations in matrix form given as

ε =











= Sσ =

1
E











ε11
ε22
ε33
2ε12
0
0











1
ν
−
ν
−
0
0
0

ν
−
1
ν
−
0
0
0

ν
−
ν
−
1
0
0
0

0
0
0
2(1 + ν)
0
0

0
0
0
0
2(1 + ν)
0

0
0
0
0
0
2(1 + ν)





















σ11
σ22
0
σ12
0
0











(58)

The three stress components which is assumed to be zero in the stress vector indicate that their associated
columns in the compliance matrix S can be neglected (i.e., column 3, 5 and 6). Furthermore, if we neglect
the rows associated with the strain components in the thickness direction εi3 (i = 1, 2, 3), the stress-strain
compliance relationship reduces to a 3

3 matrix given as

×

13

ε =





ε11
ε22
2ε12



 = Sσ =

1
E





1
ν
−
0

ν
−
1
0





0
0
2(1 + ν)





σ11
σ22
σ12





(59)

Then, the elasticity matrix C for plane stress is found by inverting the plane stress compliance matrix S
and given by Hooke’s law as

σ =





σ11
σ22
σ12



 = Cε =



1
ν

0

ν
1
0

E

ν2

1

−





0
0

ν

1

−
2





ε11
ε22
2ε12





(60)

It is important to note that the elasticity matrix for plane stress is not found by removing columns and rows
from the general isotropic elasticity matrix in Eq.(3).

Finally, in the particular case of plane stress the principal stresses (σ1, σ2) may be calculated as

σ1 =

σ11 + σ22
2

+

σ2 =

σ11 + σ22
2

−

(cid:115)(cid:18) σ11 −
(cid:115)(cid:18) σ11 −

2

2

(cid:19)2

σ22

+ σ2
12

(cid:19)2

σ22

+ σ2
12

(61)

where the major principal stress σ1 is used in the CL failure criterion in Eq.(46).

3.2. Through thickness strain increment

Although the stress in the thickness direction σ33 is assumed to be zero, this is not the case for the corre-
sponding strain ε33. The expression for the change in thickness during plane stress may be found by using
the relationship between the volumetric strain εv(= εkk) and the hydrostatic pressure P from the theory of
three-dimensional elasticity, i.e.,

εv = εkk =

1
K

−

P, K = λel +

2
3

µel =

3(1

2ν)

(62)

E

−

where P =
equal to the hydrostatic stress σH but has the opposite sign.

1
3 σkk is the mean stress and K is the bulk modulus. Note that the hydrostatic pressure P is

−

The physical meaning of εkk is the relative change of volume

εkk = ε11 + ε22 + ε33 =

dV
V

(63)

where experimental evidence has shown that plastic deformation is volume preserving (dV = 0) even at high
hydrostatic pressures [10]. This means that εkk = 0 and that the material is incompressible during plastic
deformation. Thus, the strain through the thickness for plates and shells is given as

ε33 =

−

(ε11 + ε22)

(64)

where it is seen that the strain in the thickness direction ε33 is a result of the in-plane strains ε11,22. This is
typical for thin-walled structures where there are no constraints for the strain in the thickness direction. The
only inconsistency is that in the constitutive equations for plates and shells the thickness is considered to
be constant while Eq.(64) shows that in reality there will be a small change. This illustrates the importance
of deﬁning the thickness strain increment ∆ε33 in the material routine to predict the actual change in the
thickness of the element.

14

Using the assumption of a plastic incompressible material (i.e., εp
the thickness (due to plastic thinning) is given as

kk = 0) the plastic strain increment through

where ∆εp

11 and ∆εp

22 could be expressed using strain decomposition

∆εp

33 =

−

(∆εp

11 + ∆εp

22)

and the elastic strain components are given by Hooke’s law in Eq.(59)

∆εp
∆εp

11 = ∆ε11 −
22 = ∆ε22 −

∆εe
11
∆εe
22

∆εe

11 =

∆εe

22 =

1
E
1
E

(∆σ11 −
(∆σ22 −

ν∆σ22)

ν∆σ11)

By combining Eqs.(65) to (67), the plastic strain increment through the thickness is given by

∆εp

33 =

−
= ∆εe
1
E
1

=

=

(∆ε11 −
11 + ∆εe

(∆σ11 −
ν
−
E

∆εe

22)

∆εe

11)

(∆ε22 −
∆ε22

−
∆ε11 −
22 −
ν∆σ22 + ∆σ22 −

(∆σ11 + ∆σ22)

(∆ε11 + ∆ε22)

−

ν∆σ11)

−

(∆ε11 + ∆ε22)

(65)

(66)

(67)

(68)

Then, when the plastic strain increment ∆εp
is ﬁnally found by strain decomposition

33 is known the total through thickness strain increment ∆ε33

33 + ∆εp
∆ε33 = ∆εe
ν
E

=

33

(∆σ11 + ∆σ22) +

−
1

2ν

−
E

(∆σ11 + ∆σ22)

(∆ε11 + ∆ε22)

−

1

ν

−
E

(∆σ11 + ∆σ22)

(∆ε11 + ∆ε22)

−

(69)

=

4. Uniaxial stress case

The case of uniaxial stress, which occurs in bars and in some types of beams, can be treated similarly to the
two-dimensional stress analysis case presented in Section 3.1.

In this case one assumes that the two lateral normal stresses are zero (σ22 = σ33 = 0), so that the corre-
sponding strain increments ∆ε22 and ∆ε33 must be computed and returned by the consitutive routine.

We assume that the strain increments can be decomposed into elastic and plastic parts, so that

11 = ∆ε11 −
22 = ∆ε22 −
33 = ∆ε33 −
and that the plastic part of the deformation occurs at constant volume (dV = 0), due to plastic incompress-
ibility

(70)

∆εe
11
∆εe
22
∆εe
33

∆εp
∆εp
∆εp

∆εp

11 + ∆εp

22 + ∆εp

33 = 0

15

(71)

In the general 3D case, the elastic strain increments are related to the total stress increments by Hooke’s
law

∆εe

11 =

∆εe

22 =

∆εe

33 =

1
E
1
E
1
E

[∆σ11 −
[∆σ22 −
[∆σ33 −

ν (∆σ22 + ∆σ33)]

ν (∆σ11 + ∆σ33)]

ν (∆σ11 + ∆σ22)]

(72)

However, in the uniaxial stress case the lateral normal stresses are assumed to vanish (σ22 = σ33 = 0) and
therefore also the stress increments vanish (∆σ22 = ∆σ33 = 0). Eq. (72) then reads

∆εe

11 =

∆εe

22 =

∆εe

33 =

1
E

−

−

which implies that

∆σ11
ν
E
ν
E

∆σ11 =

∆σ11 =

ν∆εe
11

ν∆εe
11

−

−

22 = ∆εe
and by utilizing the symmetry it is also assume that

∆εe

33 =

ν∆εe
11

−

∆εp

22 = ∆εp

33

and therefore

Then, by replacing Eq. (75) into Eq. (71) one obtains

∆ε22 = ∆ε33

ε22 = ε33

2∆εp

22 = 2∆εp

33 =

∆εp
11

−

∆εp

22 = ∆εp

33 =

1
2

−

∆εp

11 =

1
2

−

(∆ε11 −

∆εe

11)

Using Eq. (73) this now becomes

(cid:18)
∆ε11 −
Finally, the total lateral strain increment ∆ε22 (and ∆ε33) is obtained by strain decomposition

22 = ∆εp

∆ε11 +

1
2E

∆σ11

∆σ11

33 =

∆εp

1
E

1
2

1
2

−

=

−

(cid:19)

22 = ∆εe

33 + ∆εp

33

=

∆ε11 +

∆ε22 = ∆ε33 = ∆εe
1
2
−
1
2E
1

∆σ11 −
2ν
−
2E

22 + ∆εp
1
2E
ν
E
∆σ11 −
16

∆σ11 −
1
∆ε11 −
2
1
∆ε11
2

=

=

ν∆εe
11

∆ε11

(73)

(74)

(75)

(76)

(77)

(78)

(79)

(80)

(81)

5. Numerical Return Mapping

Wilkins [23] was one of the ﬁrst to introduce a return mapping procedure, and a variety of schemes have
been proposed following this work (see e.g.
In general, it is possible to categorize the
[2, 5, 4, 3, 24]).
stress integration methods in two diﬀerent approaches, i.e., the forward Euler method and the backward
Euler method. The choice of method is illustrated in Figure 7 and is mainly based on the assumption of
the direction of the plastic ﬂow n = ∂f /∂σ. The forward Euler method uses the known stress state at the
previous or trial conﬁguration to determine the plastic ﬂow direction n [2] (see Figure 7a). This assumption
is therefore only valid for very small time steps, and is therefore suitable for the explicit time integration
FE method. Since the direction is known there is only one unknown during the return mapping procedure,
i.e., the plastic multiplier λ. The backward Euler method evaluates the normal n to the yield surface at the
current unknown stress state [3] (see Figure 7b), which gives good accuracy even for larger time steps making
it appropriate for both the explicit and implicit time integration FE methods. However, the variation of the
normal to the yield surface must be taken into consideration during the return mapping procedure resulting
in a more complex integration, due to the involvement of possible 2nd derivatives of the yield function. In
the particular case of isotropic von Mises plasticity, the accuracy of the forward and backward Euler method
is found to be equal [25]. It should also be emphasized that when considering viscoplasticity the plastic
strain rate will change during the local iterations in the forward Euler method. That is, the plastic strain
rate will increase during the iterative return mapping (due to an increase in plastic strain). This implies
that the return to the yield surface (f = 0) occurs at strain rates that are too low. However, the plastic
strain rate at the converge increment ∆tn+1 should be the same as in the backward Euler method.

σtrial
n+1
σ1
n+1
σ2

n+1

σn+1

σn

fn+1 = 0

σtrial
n+1
σ1
n+1
σ2
n+1
σn+1

σn

fn+1 = 0

(a) Explicit stress update using forward Euler
(cutting plane method)

(b) Fully implicit stress update using backward Euler
(closest point projection method)

Figure 7: The elastic predictor viscoplastic corrector return map algorithms for the von Mises yield criterion
and associated plasticity. Starting from an elastic stress state at the last converged increment tn before
arriving at the new converged increment tn+1.

During both approaches an elastic predictor uses the total strain increment and the linear elastic relation
to obtain a trial stress state. Then, if the elastic trial stress is outside the plastic domain, the elastically
predicted stresses are corrected to ﬁt a suitably updated yield surface by iterations of the plastic corrector.
Thus, the proposed integration algorithms are formulated as an elastic predictor viscoplastic corrector stress
update using a return mapping scheme. The solution is strain controlled in the sense that the total (co-
rotational) strain increment ∆εn+1 = Dn+1∆tn+1, within the time interval ∆tn+1 = tn+1 −
tn, is known
from the global equilibrium and used to update the stresses and internal variables at step n + 1. Then, as in

17

a rate-independent case, it makes sense to ﬁrst compute an elastic trial state by assuming that the material
behaviour is purely elastic within the interval. If the trial state is within the elastic domain (f (σ)
0),
no viscoplastic ﬂow takes place within the considered time step and the trial state is the actual state at
the end of the step. Otherwise, the evolution of stresses and internal variables is computed by means of a
suitable return mapping method. It should be emphasized that the consistency condition (f (σ) = 0) no
longer holds in the viscoplastic case, since the updated stress state at tn+1 generally lies outside the yield
surface (f (σ) > 0). However, the terminology viscoplastic return mapping is justiﬁed in the present case
since the updated stresses are obtained by returning (or moving) the trial stress towards the yield surface.
Hence, the application of the procedure is similar to that in the rate-independent case, and Eq. (23) reads

≤

−
where the stress update at time tn+1 is expressed as

∆σ = C : (∆ε

∆εp) = ∆σe + ∆σp

σn+1 = σn + ∆σn+1 = σn + C : ∆εn+1 −
in which the viscoplastic corrector ∆σp returns the stress to the yield surface.

∆λn+1C : nn+1 = σtrial

n+1 −

(82)

(83)

∆σp

All internal variables and functions in Eqs.(20), (23), (34), (36), (42), (43), (45) and (46) are therefore
adopted to the return mapping scheme for integration, and replaced with the corresponding incremental
values within and at the end of the considered interval ∆tn+1. The incremental form of the system of ﬁrst
order diﬀerential equations then reads

εn+1 = εn + ∆εn+1

εp
n+1 = εp

n + ∆εp

n+1 = εp

λn+1 = λn + ∆λn+1 = λn + ˙p0∆tn+1

n + ∆λn+1nn+1 = εp
n +
(cid:34)(cid:18) ϕn+1
σy,n+1
∆εp

∆λn+1

(cid:48)
σ
n+1
σeq,n+1

3
2
(cid:19) 1

C

(cid:35)

1

−
n+1) = σtrial

σn+1 = σn + ∆σn+1 = σn + C : (∆εn+1 −
Rn+1 = Rn + ∆Rn+1 = Rn + ∆λn+1HR,n+1 = Rn + ∆λn+1

n+1 −
2
(cid:88)

∆λn+1C : nn+1

QkCkexp(

Ckλn+1)

−

(84)

σy,n+1 = σ0 + Rn+1

Tn+1 = Tn + ∆Tn+1 = Tn +

Dn+1 = Dn + ∆Dn+1 = Dn +

fv,n+1 = ϕn+1 −

σy,n+1

1 +

(cid:18)

χσeq,n+1∆λn+1
ρCp
σ1,n+1(cid:105)
(cid:104)
∆λn+1
˙p0∆tn+1

∆λn+1
(cid:19)C (cid:18)

1

−

k=1

(cid:21)m(cid:19)

(cid:20) Tn+1 −
Tr
Tr
Tm −

where ∆λn+1 is the incremental plastic multiplier at time tn+1.

The set of nonlinear algebraic equations in Eq. (84) are solved by updating the normal to the yield surface
gradient nn+1 iteratively, and the viscoplastic corrector is applied from the elastic predictor σtrial
n+1 throughout
the iteration until fv,n+1 = 0.

5.1. Cutting plane method

The present study focus mainly on the cutting plane method originally proposed by Ortiz and Simo [2]. The
basis for this return mapping is an explicit elastic predictor plastic corrector stress update using a forward
Euler scheme. That is, the objective of the numerical integration is to use the known stress state at the

18

previous or trial conﬁguration n to determine the new converged stress state at step n + 1 (Figure 7a).

By using the normal ni
n+1 at the previous stress state (see Figure 7a) and assuming weak coupling of the
temperature, i.e., using Tn instead of Tn+1, there is only one unknown during the return mapping procedure
and that is the plastic multiplier ∆λi+1
n+1. This reduces the fully set of equations in Eq. (84) to a single
equation which can be solved by a linear Taylor-expansion around fn+1(σ, X) = fn+1(σ, σy, υp( ˙p), Γp(Tn)),
i.e.,

f i+1
n+1 = f i

n+1 +

∂f
∂σ

(cid:12)
(cid:12)
(cid:12)
(cid:12)

i

n+1

∆σn+1 +

∂f
∂X

(cid:12)
(cid:12)
(cid:12)
(cid:12)

i

n+1

∆Xn+1

(85)

which may be written in terms of the internal variables X as

f i+1
n+1 = f i

∂f
∂σy,n+1
where the ﬂow stress σy,n+1 and the following expressions are chosen as the internal variables X in the
Taylor-expansion

∂f
∂υp,n+1

∂f
∂σn+1

δυp,n+1 = 0

δσy,n+1 +

δσn+1 +

n+1 +

(86)

(cid:18)

υp(∆λn+1) =

1 +

(cid:19)C

∆λn+1
˙p0∆tn+1

(cid:18)

Γp(Tn) =

1

−

(cid:20) Tn −
Tm −

(cid:21)m(cid:19)

Tr
Tr

(87)

Introducing the partial derivatives of the internal variables in Eq. (86) this now reads

f i+1
n+1 = f i

n+1 −

∂f
∂σn+1

C

∂f
∂σn+1

δλi+1

n+1 −

υp,n+1Γp,nHR,n+1δλi+1

n+1 −

σy,n+1Γp,n

∂υp,n+1
∂ ˙λn+1

δ ˙λi+1

n+1 = 0

(88)

Moreover, using ˙p = ˙λ = ∆λ
υp(∆λ) and Eq. (88) may be written as

∆t the internal variable controlling the viscoplasticity υp( ˙p) may be expressed as

f i+1
n+1 = f i

n+1 −

∂f
∂σn+1

C

∂f
∂σn+1

δλi+1

n+1 −

υp,n+1Γp,nHR,n+1δλi+1

n+1 −

σy,n+1Γp,n

∂υp,n+1
∂∆λn+1

δλi+1

n+1 = 0

(89)

Finally, the incremental change δλ in the incremental plastic multiplier ∆λ is then expressed as

δλi+1

n+1 =

∂f
∂σn+1

C ∂f

∂σn+1

+ υp,n+1Γp,nHR,n+1 + σy,n+1Γp,n

∂υp,n+1
∂∆λn+1

f i
n+1

where

∂υp(∆λn+1)
∂∆λ

=

C
˙p0∆tn+1

(cid:18)

1 +

∆λn+1
˙p0∆tn+1

(cid:19)C

1

−

HR,n+1 =

dRn+1
dpn+1

=

2
(cid:88)

k=1

QkCkexp(

Ckλn+1)

−

(cid:27)T

(cid:26) ∂f
∂σ

[C]

(cid:27)

(cid:26) ∂f
∂σ

=

∂f
∂σij

Cijkl

∂f
∂σkl

= C1111

(cid:34)(cid:18) ∂f
∂σ11

(cid:19)2

(cid:19)2

+

(cid:18) ∂f
∂σ22

+

(cid:18) ∂f
∂σ33

(cid:19)2(cid:35)

+ 2C1122

+ 4C1212

∂f
∂σ22
(cid:19)2

(cid:20) ∂f
∂σ11
(cid:34)(cid:18) ∂f
∂σ12

19

+

+

∂f
∂σ22
(cid:18) ∂f
∂σ23

∂f
∂σ33
(cid:19)2

+

∂f
∂σ33
(cid:18) ∂f
∂σ31

+

(cid:21)

∂f
∂σ11
(cid:19)2(cid:35)

(90)

(91)

(92)

(93)

For more details regarding the linear Taylor-expansion, it is referred to Appendix A. It is also emphasized
that since the time increment ∆tn+1 is often very small during explicit time integration, the change in tem-
perature ∆Tn+1 is assumed to have a negligible eﬀect on the update of the state variables. Therefore, the
temperature is assumed to have no evolution during the time increment ∆tn+1, and a low coupling approach
for the temperature is applied. This implies that Tn is used in the return mapping to update the stress
and internal variables at time tn+1, and then the resulting state variables at tn+1 are used to calculate and
update Tn+1. Thus, this assumption simpliﬁes the calculations leading to Eq. (90) without any signiﬁcant
loss of accuracy in the model. This is also the case for the normal nn+1 = ∂f /∂σn+1 to the yield surface
where the use of the previous normal ni
n+1 violates the associated ﬂow rule in Eqs. (36) and
(84) since a unique normal for the given strain increment ∆εn+1 must be found at each stress state. Thus,
the deformation does not follow the minimum plastic work path when the cutting plane method is used (see
Figure 7). However, when the time increment ∆tn+1 is very small this is assumed to have negligible eﬀect
on the stress state.

n+1 instead of ni+1

The numerical scheme within the time increment ∆tn+1 may then be summarized as follows (see also
Appendix E):

1. Set the initial values of internal variables to the converged values of the previous step at tn, check if

the integration point under evaluation has already failed, and calculate the speed of sound.

2. Compute the elasticity tensor.

That is, IF a 3D stress state is used THEN use Eqs. (3) and (4)

C = λelI

⊗

I + 2µelI =

C1 =

ν)E

(1
(1 + ν)(1

−

2ν)

−













C1 C4 C4
C4 C1 C4
C4 C4 C1
0
0
0
0
0
0

0
0
0
0 C7
0
0
νC1
ν
1

, C4 =

0
0
0
0
0 C7
0









0
0
0
0
0
0 C7
E
1
1 + ν
2

, C7 =

−

ELSEIF a 2D stress state is used THEN use Eq. (60)

C =





C1 C4
C4 C1
0

0
0
0 C7



 , C1 =

E

1

−

ν2 , C4 = νC1, C7 =

1
2

E
1 + ν

ELSE a 1D stress state is used THEN use

C =





C1
0
0

0
0
0


0
0
 , C1 = E
0

3. Use the total stain increment ∆εn+1 from the global equilibrium to compute the elastic predictor

σtrial

n+1 = σn + C : ∆εn+1

4. Compute the von Mises equivalent stress in terms of the elastic trial stress and Eq. (35)

eq,n+1 = ϕ(σtrial
σtrial

n+1 ) =

(cid:114)

3
2

n+1 : σ(cid:48)trial
σ(cid:48)trial
n+1

20

5. Check if temperature softening is included in the material input (i.e., m

= 0).

6. Set the incremental plastic multiplier ∆λn+1 equal to zero.

7. Check for plastic admissibility, i.e., IF f (σtrial

n+1 ) > 0 THEN apply return mapping by using the cutting
plane algorithm (SOLU =1) to ﬁnd the incremental change in the plastic multiplier ∆λn+1. Note that
superscript i denotes the local iteration counter in the following.

(i) Compute the normal nn+1 to the yield surface based on the initial (i.e., trial) or previous stress

conﬁguration.

(ii) Compute hardening modulus HR,n+1 according to Eq. (42)

nn+1 =

(cid:48)
σ
n+1
σeq

3
2

H i

R,n+1 = Q1C1 exp(

C1pi

n+1) + Q2C2 exp(

−

C2pi

n+1)

−

(iii) Compute denominator in Eq. (90) using Eqs. (91), (92), (93) and (98)

DENOM =

∂f
∂σi

n+1

C

∂f
∂σi

n+1

+ υi

p,n+1Γp,nH i

R,n+1 + σi

y,n+1Γp,n

p,n+1

∂υi
∂∆λi

n+1

(iv) Compute incremental change in the plastic multiplier δλi+1

n+1 in Eq. (90) and update variable

∆λi+1
n+1

δλi+1

n+1 =

∂f
∂σn+1

C ∂f

∂σn+1

+ υp,n+1Γp,nHR,n+1 + σy,n+1Γp,n

∂υp,n+1
∂∆λn+1

f i
n+1

∆λi+1

n+1 = ∆λi

n+1 + δλi+1
n+1

(v) Update internal variables dependent on ∆λi+1
n+1
pi+1
n+1 = pi

n+1 + δλi+1
n+1

Ri+1

n+1 = Rn + H i

R,n+1∆λi+1
n+1

y,n+1 = A + H i+1
σi+1

R,n+1∆λi+1

n+1

(vi) Update stress components using the elastic Hooke’s law according to Eq. (83)

Note that the material routine uses a a fully vectorized representation of the stress and strain
vectors, i.e.,

σi+1
n+1 = σtrial

n+1 −

∆λi+1

n+1C : ni+1
n+1

σ = [σ11, σ22, σ33, σ12, σ23, σ31]T

ε = [ε11, ε22, ε33, γ12, γ23, γ31]T = [ε11, ε22, ε33, 2ε12, 2ε23, 2ε31]T

Moreover, the application of objective stress rates (i.e., the Jaumann rate in EPX) to update the
stress tensor σ in plane stress states allow for the possibility that the normal stress σ33 will not
be zero (outside the material routine). The material subroutine should therefore initialize the
normal stress σ33 to zero (σ33 = 0).

(vii) Compute the updated von Mises equivalent stress using Eq. (35)

eq,n+1 = ϕ(σi+1
σi+1

n+1) =

(cid:114)

3
2

21

σ

(cid:48),i+1
n+1 : σ

(cid:48),i+1
n+1

(cid:54)
(viii) Update the yield function in Eq. (84) according to the new stress state computed in (vi) and (vii)

f i+1
v,n+1 = ϕi+1

n+1 −

σi+1
y,n+1

1 +

(cid:32)

(cid:33)C (cid:18)

∆λi+1
n+1
˙p0∆tn+1

(cid:21)m(cid:19)

1

−

(cid:20) Tn −
Tm −

Tr
Tr

(viiii) Check for convergence, i.e., IF

THEN continue to 7.
ELSE i = i + 1 and GO TO (i)

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)

fv(∆λi+1
n+1)
σy,n+1υpΓp

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)

< TOL

8. Update the temperature according to the low coupling approach using Eq. (45) and check IF T

Tn+1 = Tn +

χσeq,n+1∆λn+1
ρCT

Tm

≥

9. Compute the hydrostatic and deviatoric stress respectively according to Eq. (32)

σH,n+1 =

1
3

tr(σn+1)

(cid:48)

n+1 = σn+1 −
σ

σH,n+1I

10. Computing the principal stresses σi,n+1 based on the principal deviatoric stresses σ

(cid:48)

i,n+1 as shown in

Eqs. (55) and (57)

(cid:114)

(cid:48)

1,n+1 = 2
σ

J2,n+1
3

cos

(cid:48)

(cid:48)

σ

σ

2,n+1 = 2

3,n+1 = 2

(cid:114)

(cid:114)

J2,n+1
3

cos

(cid:18) θ

3 −

J2,n+1
3

cos

(cid:18) θ
3

+

2π
3

2π
3

θ
3

(cid:19)

(cid:19)

(cid:48)
1,n+1

σ

(cid:48)
2,n+1

σ

≤

≤

where

cos θ =

J3,n+1
J 3
2,n+1/27

,

(cid:113)
2

0

θ

≤

≤

π

σi,n+1 = σ

i,n+1 + σH,n+1 = σ

(cid:48)

(cid:48)

i,n+1 +

1
3

tr(σn+1)

11. Compute the damage parameter Dn+1 by using Eq. (46), i.e., integrating the major principal stress

during the entire equivalent plastic strain path

Wn+1 = Wn + max(σ1,n+1, 0)∆λn+1

Dn+1 =

Wn+1
Wc

12. Update state variables to be returned to the FE code and check for element erosion, i.e., IF Dn+1 > DC

THEN the material point stiﬀness is deleted.

22

13. Update the (actual) total strain ε33,n+1 and incremental strain ∆ε33,n+1 through the thickness accord-

ing to Eq. (69) and return this to the FE code

∆ε33,n+1 =

1

2ν

−
E

(∆σ11,n+1 + ∆σ22,n+1)

(∆ε11,n+1 + ∆ε22,n+1)

−

ε33,n+1 = ε33,n + ∆ε33,n+1

In a FE code, this is typically done for each material integration point in the FE assembly. The stresses are
updated in the material interface using a fully vectorized version of the forward Euler integration algorithm
and a two-state architecture where the initial values at tn are stored in the old arrays and the new values
at tn+1 must be updated and stored in the new arrays returned to the global FE analysis. The VPJC
material routine is valid for 1D, 2D and 3D stress states, i.e., for bar, shell, solid, axisymmetric, plane strain
and plane stress elements. The stresses for 3D elements are stored similar to that of symmetric tensors
as σ = [σ11, σ22, σ33, σ12, σ23, σ31]T , and plane stress, axisymmetric and plane strain elements are stored as
σ = [σ11, σ22, σ33, σ12]T . The deformation gradient and strains are stored similarly to the stresses, however,
one should be aware of that the shear strain is stored as engineering shear strains, e.g. γ12 = 2ε12.

5.2. Radial return method

The objective of the implicit backward Euler integration algorithm is to ﬁnd the minimum distance from
the trial state to the ﬁnal state at step n + 1 (see Figure 7b). In the particular case of a 3D stress state, von
Mises yield criterion with associated ﬂow rule and isotropic hardening, the fully implicit set of equations in
Eq. (84) is reduced to a single equation and the radial return algorithm. This is due to the circular cylindrical
shape of the yield surface in Figure 4a, which implies that the gradient of the yield surface at the trial state
ntrial
n+1 is co-linear with the corresponding gradient at the ﬁnal stress state nn+1 (Figure 8). The return to
the yield surface from the elastic trial state is radial in the deviatoric plane, which simpliﬁes the algorithm
and makes the return mapping exceptionally stable and accurate. This makes the radial return mapping one
of the most widely used integration algorithms for plain strain and 3D J2 (or von Mises) elastoplasticity [26].

By considering the deviatoric stress space, the equivalent stress could then be expressed as (see also Appendix
B)

σeq = ϕ(σ

n+1) = ϕ(σ

(cid:48)

(cid:48)trial
n+1 )

3G∆λn+1,

G =

−

E
2(1 + ν)

(94)

and the system of equations in Eq. (84) is only dependent on the incremental plastic multiplier ∆λn+1,
where the viscoplastic constraint fv,n+1(σ) = 0 reads

fv,n+1(σ) = ϕtrial

n+1 (σ)

3G∆λn+1 −

−

σy,n+1

1 +

(cid:18)

∆λn+1
˙p0∆t

(cid:19)C (cid:18)

1

−

(cid:20) Tn+1 −
Tr
Tr
Tm −

(cid:21)m(cid:19)

= 0

(95)

Then, the incremental plastic multiplier ∆λn+1 can be solved by means of a local Newton-Raphson iteration,
i.e.,

∆λi+1

n+1 = ∆λi

n+1 −

fv(∆λi
∂fv(∆λi
∂∆λ

n+1)
n+1)

(96)

where (again) i is the local iteration counter.

Since the time increment ∆tn+1 is often very small during explicit time integration, the change in temperature
∆Tn+1 is assumed to have a negligible eﬀect on the update of the state variables. Hence, the temperature
is assumed to have no evolution during the time increment ∆tn+1, and a low coupling approach for the
temperature is applied. This implies that Tn is used in the return mapping to update the stress and internal

23

π-plane

∂f
∂σ

σ3

σ

(cid:48)trial
n+1

(cid:12)
(cid:12)
(cid:12)n+1

(cid:48)trial
+1
n
σ
)
(cid:48)trial
∆λn
1
+
3G
n
ϕ(σ

1
+

σ

(cid:48)
n+1

(cid:48)
σ
n

σ1

surface at tn

surface at tn+1

σ2

Figure 8: The implicit elastic predictor viscoplastic corrector return map algorithm for the von Mises yield
criterion and associated plasticity in the deviatoric plane.

variables at time tn+1, and then the resulting state variables at tn+1 are used to calculate Tn+1. This
assumption reduces the number of equations in Eq. (84), and the residual derivative in Eq. (96) is expressed
as

∂fv(∆λn+1)
∂∆λ

=

=

where

(cid:16)

(cid:104)

ϕ

∂
∂∆λ

(cid:17)

σ

(cid:48)trial
n+1

3G

−

−

∂σy,n+1
∂∆λ

−

3G∆λn+1 −
υp(∆λn+1)Γp(Tn)

(cid:105)
σy,n+1υp(∆λn+1)Γp(Tn)

σy,n+1

−

∂υp(∆λn+1)
∂∆λ

Γp(Tn)

(cid:19)C

(cid:18)

υp(∆λn+1) =

1 +

∂υp(∆λn+1)
∂∆λ

=

Γp(Tn) =

C
˙p0∆t
(cid:18)
1

−

∆λn+1
˙p0∆t
(cid:18)

1 +

(cid:20) Tn −
Tm −

(cid:19)C

1

−

∆λn+1
˙p0∆t
Tr
Tr

(cid:21)m(cid:19)

The numerical scheme within the time increment ∆tn+1 may then be summarized as follows:

(97)

(98)

1. Set the initial values of internal variables to the converged values of the previous step at tn, and set
the incremental plastic multiplier ∆λn+1 equal to zero. Use the total stain increment ∆εn+1 from the
global equilibrium to compute the elastic predictor

2. Compute the equivalent stress in terms of the elastic trial stress

σtrial

n+1 = σn + C : ∆εn+1

σtrial
eq = ϕ(σ

(cid:48)trial
n+1 ) =

(cid:114)

3
2

n+1 : σ(cid:48)trial
σ(cid:48)trial
n+1

and the radial normal vector for the yield surface

nn+1 =

3
2

(cid:48)trial
σ
n+1
σtrial
eq

24

3. Check for plastic admissibility, i.e., IF f (σ

(cid:48)trial
n+1 ) > 0 THEN perform local Newton-Raphson iteration

to ﬁnd ∆λn+1

(i)

(ii)

∆λi+1

n+1 = ∆λi

n+1 −

fv(∆λi
∂fv (∆λi
∂∆λ

n+1)
)

n+1

IF

(cid:12)
(cid:12)
(cid:12)
(cid:12)

fv(∆λi+1
n+1)
σy,n+1υpΓp

(cid:12)
(cid:12)
(cid:12)
(cid:12)

> TOL THEN i = i + 1 and GO TO (i)

4. Update internal variables dependent on ∆λn+1. Note that only the deviatoric stress will change during
the radial return, and that the hydrostatic stress could be calculated based on the elastic trial stress.

5. Update stress components using the elastic Hooke’s law, i.e.,

σn+1 = σtrial

n+1 −

∆λn+1C : nn+1

6. Check for element erosion, i.e., IF D > DC THEN the material point stiﬀness is deleted.

7. Update state variables to be returned to the global analysis in EPX.

As for the cutting plane method in Section 5.1, this is done for each material integration point in the FE
assembly and the stresses are updated using a two-state architecture and returned to the global analysis in
EPX. It must be emphasized that the radial return method is only valid in the case of a 3D stress state
(Figure 4), and therefore limited to solid, axisymmetric and plane strain elements. The stresses for three-
dimensional elements are stored similar to that of symmetric tensors as σ = [σ11, σ22, σ33, σ12, σ23, σ31]T , and
plane strain and axisymmetric elements are stored as σ = [σ11, σ22, σ33, σ12]T . The deformation gradient
and strains are stored similarly to the stresses, however, one should be aware of that the shear strain is
stored as engineering shear strains, e.g. γ12 = 2ε12.

6. Veriﬁcation of Material Subroutine

The performance of the material subroutine is veriﬁed and validated in EPX. First, a veriﬁcation is per-
formed using single element simulations in uniaxial tension and simple shear. Then, the model is validated
against experimental data found in the literature for quasi-static tension tests [21] and airblast experiments
[27].

The material model requires the user to specify 20 parameters, where 16 are mandatory while the room (or
ambient) temperature Tr, convergence criterion TOL, maximum number of local iterations MXITER and
preferred solution algorithm SOLU are taken as default if not speciﬁed in the input ﬁle [7]. The material
constants used for the veriﬁcation and validation in the following are taken from the literature for a ductile
structural steel [21] and summarized in Table 1. Previous studies by De Borst and Feenstra [25] have shown
that the choice of return mapping algorithm is more important for orthotropic yield criteria (e.g. the Hill
criterion) than for the isotropic von Mises plasticity theory used in this study. The cutting-plane algorithm
(SOLU = 1) is the default choice and will therefore be used for the return mapping because this is robust
and readily implemented for all stress states (1D, 2D and 3D).

6.1. Singel element veriﬁcation

An extensive single element veriﬁcation of the VPJC material model is performed in Reference [28]. The
following single element veriﬁcation is therefore mainly included for completeness of this study, and it is
referred to [28] for more detailed single element veriﬁcations. Since the material model is applicable for one-
(1D), two- (2D) and three-dimensional (3D) stress states, the scope of this work is limited to one represen-
tative bar element in 1D (FUN3 ), one continuum element in 2D (Q42L), one shell element in 3D (Q4GS )
and one continuum element in 3D (CUBE ). This is based on the ﬁndings in [28], which concluded that the

25

Table 1: Material parameters used in numerical simulations [7, 21].

Elastic constants
and density
ν
[ - ]
0.33

E
[MPa]
210000

ρ
kg/m3
7850

A
[MPa]
370

Yield stress
and strain hardening
Q1
[MPa]
236.4

Q2
[MPa]
408.1

C1
[ - ]
39.3

Adiabatic heating and temperature softening
Cp
[J/kgK]
452

Tm
[K]
1800

m
[ - ]
1.0

χ
[ - ]
0.9

Tr
[K]
293

C2
[ - ]
4.5
Convergence
TOL MXITER
[ - ]
10−5

[ - ]
100

Damage

Strain-rate
hardening
C
˙p0
[s−1]
[ - ]
0.01
5e-4
Solution algorithm
SOLU
[1: cutting plane or 2: radial return]
1

WC
[MPa]
473

DC
[ - ]
1.0

results are more or less independent of element type. Note that the JAUM option must be activated when
using CEA continuum elements (e.g. Q42L and CUBE ) to treat large strains, and that the LMST option
must be included to update the thickness for the Q4GS elements during plastic straining [7].

The numerical tests of each element are performed in uniaxial tension and simple shear for the 2D and 3D
elements, while the 1D element is limited to uniaxial tests only. Note that, in addition to tests in uniaxial
tension and simple shear, it is common to validated a material subroutine in terms of uniaxial compression,
uniaxial tension/compression in oblique directions, and uniaxial tension/compression with a ﬁnite rotation.
However, since the formulation of the material model is isotropic the response is the same in tension and
compression and independent of loading direction. It is therefore considered suﬃcient to only test one of
these loading conditions. Moreover, since the co-rotational formulation in the material routine is taken care
of by EPX outside the material subroutine environment [7], there is no need to test the return map algorithm
in terms of ﬁnite rotation. The veriﬁcation of the material routine is therefore limited to the uniaxial tension
and simple shear tests in Figure 9.

(a) Uniaxial tension

(b) Simple shear

Figure 9: Illustration of the loading cases for the single element tests of the brick element (CUBE ). Note
that γ = tan(θ) = u2/L0.

26

∆LL0v3=˙Lθv2=˙u2u2L0The numerical models consist of a simple input ﬁle containing one element with element lengths of L0 = 1
mm. The tests are performed with a prescribed displacement resulting in an equivalent plastic strain of
approximately 2.0. This is obtained by applying a constant velocity throughout the analysis (Figure 9), and
the results are compared with the analytical solution for the equivalent stress in Eq. (41). In the particular
case of including viscoplasticity (C
= 0), it is convenient to test the elements at constant strain rate since
the analytical solution of the equivalent stress is readily obtained from Eq. (41). Due to the logarithmic
nature of the longitudinal strain in uniaxial tension, i.e.,
(cid:18) L
L0

ε33 = εe

33 + εp

∆L
L0

33 = ln

= ln

(99)

1 +

(cid:18)

(cid:19)

(cid:19)

It is necessary with an exponential displacement-driven loading to obtain a constant strain rate. By consid-
ering large strains and neglecting the elastic part εe
33 of the longitudinal strain, Eq. (99) may be simpliﬁed
to

which may be expressed as

where the displacement is given by

ε33 ≈

εp
33 = p = ˙pt = ln

(cid:19)

(cid:18) L(t)
L0

L(t) = L0 exp ( ˙pt)

with corresponding velocity

d33(t) = L(t)

−

L0 = L0 (exp ( ˙pt)

1)

−

v3(t) = L0 ˙p exp ( ˙pt) ,

v3(t = 0) = v3,0 = L0 ˙p

(100)

(101)

(102)

(103)

Eq. (102) was therefore used to ensure constant strain rate in uniaxial tension, by choosing typical values
of the plastic strain rate ˙p and controlling the respective termination times such that the ﬁnal displacement
was the same in all tests. The corresponding viscoplastic response in the simple shear tests was obtained by
varying the constant velocity (v2 in Figure 9) and termination time to obtain an equivalent plastic strain
of approximately 2.0. To compare the tests in simple shear with the analytical solution in Eq.(41), it was
necessary to scale the loading velocity in simple shear (v2 = constant) with respect to the initial velocity in
uniaxial tension (v3(t = 0) = v3,0). The scaling factor was obtained by considering the plastic strain rate ˙p
deﬁned in Eq. (37), which for the uniaxial test reads

˙p =

(cid:114)

2
3

Dp : Dp =

(cid:114)

(cid:104)

2
3

( ˙εp

11)2 + ( ˙εp

22)2 + ( ˙εp

33)2(cid:105)

(cid:115)

=

(cid:20) 6
4

2
3

(cid:21)

( ˙εp

33)2

= ˙εp
33

(104)

when assuming that the material is incompressible during plastic deformation (see Eq. (64)). The corre-
sponding expression in simple shear is

˙p =

(cid:114)

2
3

Dp : Dp =

(cid:118)
(cid:117)
(cid:117)
(cid:116)

(cid:34)(cid:18) 1
2

2
3

(cid:19)2

˙γ12

+

(cid:18) 1
2

˙γ21

(cid:19)2(cid:35)

=

˙γ12
√3

=

˙γ
√3

(105)

This implies that the loading velocity in simple shear was v2 = √3v3. In addition, the termination time
in the case of simple shear was scaled such that tests in uniaxial tension and simple shear had the same
equivalent plastic strain at the end of the simulation (see Table 2).

To illustrate the capabilities of the model, the model is ﬁrst tested with respect to rate-independent plas-
ticity. Then, the model is tested for its precision in predicting both viscoplastic and thermoviscoplastic

27

(cid:54)
Table 2: Loading in numerical test programme.

Loading
Velocity [m/s]
Termination [s]
˙p [s−

1]

Uniaxial tension
0.01
0.1
0.2
0.02
10
100

v3,0 = 1.0
0.002
1000

0.001
2
1

v2 = 1.732
0.0021
1000

Simple shear
0.1732
0.021
100

0.01732
0.21
10

0.001732
2.1
1

behaviour, before ﬁnally tested with respect to ductile failure. The numerical results were compared to ana-
lytical solutions for the constitutive equations, where this was available, to ensure that the overall behaviour
of the material model is correct.

6.1.1. Rate-independent plasticity
The rate-independent behaviour of the material model is veriﬁed by excluding the strain-rate and temper-
ature sensitivity. This is readily done in the material input by choosing C = 0 and m = 0. The damage
1.
threshold was also set to an exaggerated value to avoid element erosion during plastic straining, i.e., DC (cid:29)
The analytical solution then reduces to the ﬁrst term in Eq. (41), which will be used for comparison to the
numerical solution.

1D bar element

Figure 10 shows the performance of the model in describing the rate-independent behaviour in uniaxial
tension for a bar element (FUN3 ).

Figure 10: Results from strain-rate and temperature independent 1D single element (FUN3 ) test in uniaxial
1). Numerical results in EPX compared to analytical solution in Eq. (41).
tension (C = 0, m = 0 and DC (cid:29)

2D continuum element

Figure 11 shows the capability of the model in describing the rate-independent behaviour in uniaxial tension
and simple shear for a 2D continuum and plane stress condition by using the fully integrated quadrilateral
Q42L element. Note that the JAUM option must be activated when using CEA continuum elements (e.g.
Q42L) to treat large strains.

28

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq (EPX)R(p) (EPX)σeq (analytical)R(p) (analytical)(a) Uniaxial tension (Q42L)

(b) Simple shear (Q42L)

Figure 11: Results from strain-rate and temperature independent 2D single element (Q42L) tests (C = 0,
m = 0 and DC (cid:29)

1). Numerical results in EPX compared to analytical solution in Eq. (41).

3D shell element

The same tests in uniaxial tension and simple shear are performed also for a 3D (quadrilateral) shell element
(Q4GS ), where the numerical results are compared to the analytical solution from Eq. (41) in Figure 12.
Note that the LMST option must be included to update the thickness for the Q4GS elements during plastic
straining.

(a) Uniaxial tension (Q4GS )

(b) Simple shear (Q4GS )

Figure 12: Results from strain-rate and temperature independent 3D single (shell) element (Q4GS ) tests
(C = 0, m = 0 and DC (cid:29)

1). Numerical results in EPX compared to analytical solution in Eq. (41).

3D continuum element

Finally, the material model is veriﬁed for a 3D (hexahedron) continuum element (CUBE ) in uniaxial tension
and simple shear (see Figure 13). Again, note that the JAUM option must be activated when using CEA
continuum elements (e.g. CUBE ) to treat large strains.

29

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq (EPX)R(p) (EPX)σeq (analytical)R(p) (analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq (EPX)R(p) (EPX)σeq (analytical)R(p) (analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq (EPX)R(p) (EPX)σeq (analytical)R(p) (analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq (EPX)R(p) (EPX)σeq (analytical)R(p) (analytical)(a) Uniaxial tension (CUBE )

(b) Simple shear (CUBE )

Figure 13: Results from strain-rate and temperature independent 3D single element (CUBE ) tests (C = 0,
m = 0 and DC (cid:29)

1). Numerical results in EPX compared to analytical solution in Eq. (41).

It is seen that the equivalent stress σeq from the numerical tests is in excellent agreement with the analytical
solution in Eq. (41) for all stress states (Figures 10 - 13), and that the material model is capable of describing
the rate-independent behaviour during plastic straining. This is also conﬁrmed by comparing the equivalent
stress σeq and the hardening variable R(p) at the end of the simulations (p = 2.0) in Table 3. Note that
negative values of ∆σeq and ∆R(p) imply that the analytical value is larger than that in the respective
numerical test.

Table 3: Summary of single element tests for strain rate and temperature independent behaviour at the end
of the simulations (p = 2.0).

Element

Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)

σeq
[MPa]
1014.5
1015.1
1014.1
1013.9
1014.1
1013.9
1014.9
1014.0

∆σeq R(p = 2.0) ∆R(p = 2.0)
[%]
-
0.06
-0.04
-0.06
-0.04
-0.06
0.04
-0.04

[MPa]
644.5
645.1
644.1
643.9
644.1
643.9
644.9
644.0

[%]
-
0.10
-0.06
-0.09
-0.06
-0.09
0.06
-0.07

6.1.2. Viscoplasticity
Similar tests are performed by also including viscoplasticity (C = 0.01) in the model. This will verify if the
model is capable of describing the viscoplastic response at elevated strain rates by varying the plastic strain
rate ˙p according to Table 2. That is, starting with a plastic strain rate ˙p = 1 and increasing this successively
by one order of magnitude until ˙p = 1000. The constant strain rates are obtained in the numerical tests by
using Eq. (102) (or Eq. (103)). As for the rate-independent tests, the temperature sensitivity and element
1. The analytical solution then
erosion are excluded from the numerical tests by choosing m = 0 and DC (cid:29)
reduces to the two ﬁrst terms in Eq. (41), against which the numerical results will be compared. Finally,
5e−
the rate-independent response ( ˙p

4) from Section 6.1.1 is also plotted for comparison.

≤

30

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq (EPX)R(p) (EPX)σeq (analytical)R(p) (analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq (EPX)R(p) (EPX)σeq (analytical)R(p) (analytical)1D bar element

Figure 14 shows the performance of the model in describing the rate-dependent behaviour in uniaxial tension
for a bar element (FUN3 ).

Figure 14: Results from temperature independent 1D single element (FUN3 ) test in uniaxial tension (m = 0
1). Numerical results in EPX compared to the ﬁrst two terms of the analytical solution in
and DC (cid:29)
Eq. (41).

2D continuum element

Figure 15 shows the capabilities of the model in describing the viscoplastic response at elevated strain rates
for the Q42L element in uniaxial tension (Figure 15a) and simple shear (Figure 15b).

(a) Uniaxial tension (Q42L)

(b) Simple shear (Q42L)

Figure 15: Results from temperature independent 2D single element (Q42L) tests (m = 0 and DC (cid:29)
Numerical results in EPX compared to the ﬁrst two terms of the analytical solution in Eq. (41).

1).

3D shell element

The same tests are performed for the Q4GS element in uniaxial tension (Figure 16a) and simple shear

31

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)(Figure 16b), where the numerical results are compared to the ﬁrst two terms of the analytical solution in
Eq. (41).

(a) Uniaxial tension (Q4GS )

(b) Simple shear (Q4GS )

Figure 16: Results from temperature independent 3D single (shell) element (Q4GS ) tests (m = 0 and
1). Numerical results in EPX compared to the ﬁrst two terms of the analytical solution in Eq. (41).
DC (cid:29)

3D continuum element

Finally, the viscoplastic response of the material model is veriﬁed for the CUBE element in uniaxial tension
(Figure 17a) and simple shear (Figure 17b).

(a) Uniaxial tension (CUBE )

(b) Simple shear (CUBE )

Figure 17: Results from temperature independent 3D single (continuum) element (CUBE ) tests (m = 0 and
1). Numerical results in EPX compared to the ﬁrst two terms of the analytical solution in Eq. (41).
DC (cid:29)

It is seen that the equivalent stress σeq from the numerical tests is in excellent agreement with the ﬁrst
two terms of the analytical solution in Eq. (41) for all stress states (Figures 14 - 17), and that the material
model is capable of describing the viscoplastic response at elevated strain rates. This is also conﬁrmed by
comparing the equivalent stress σeq and the hardening variable R(p) at the end of the simulations (p = 2.0)

32

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)in Table 4. By comparing to the results in Table 3 it is seen that the hardening variable R(p) is only
dependent on the plastic straining and independent of strain rate (as expected). Again, note that negative
values of ∆σeq and ∆R(p) imply that the analytical value is larger than that in the respective numerical
test. Also note that the analytical reference in Table 4 is found using the plastic strain and temperature
history from the uniaxial tests with the FUN3 element, while the values given for ∆σeq and ∆T (p = 2.0)
for each element and numerical test are found using the corresponding histories from the respective tests.
That is, the analytical reference used to compute ∆σeq and ∆T (p = 2.0) in Table 4 is found using the plastic
strain and temperature history from the respective test.

33

Table 4: Summary of single element tests for temperature independent behaviour at the end of the simula-
tions (p = 2.0).

∆σeq R(p = 2.0) ∆R(p = 2.0)
[%]
-
0.06
-0.04
-0.06
-0.04
-0.06
0.04
-0.04
-
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
0.0
-
-0.08
-0.07
-0.06
-0.07
-0.06
-0.05
-0.04

[MPa]
644.5
645.1
644.1
643.9
644.1
643.9
644.9
644.0
644.4
644.4
644.4
644.5
644.4
644.5
644.5
644.5
644.4
644.4
644.4
644.5
644.5
644.5
644.5
644.5
644.4
644.4
644.4
644.4
644.4
644.4
644.4
644.4
644.4
643.6
643.7
643.9
643.7
643.9
643.9
644.0

[%]
-
0.10
-0.06
-0.09
-0.06
-0.09
0.06
-0.07
-
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-
-0.02
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
-
-0.13
-0.11
-0.09
-0.11
-0.09
-0.08
-0.07

Plastic strain rate
[s−

1]

4

5e−

≤

1

10

100

1000

Element

Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)

σeq
[MPa]
1014.5
1015.1
1014.1
1013.9
1014.1
1013.9
1014.9
1014.0
1094.6
1094.6
1094.6
1094.6
1094.6
1094.6
1094.6
1094.6
1120.1
1120.1
1120.1
1120.1
1120.1
1120.1
1120.1
1120.1
1146.1
1146.1
1146.1
1146.1
1146.1
1146.1
1146.1
1146.1
1172.8
1171.9
1172.0
1172.2
1172.0
1172.2
1172.2
1172.4

34

6.1.3. Thermoviscoplasticity
Now, thermoviscoplasticity is also included in the model by specifying m = 1.0 in the material input (ac-
cording to Table 1). The temperature T from EPX was used as input in Eq. (41) to obtain the analytical
solution for the equivalent stress σeq, while the incremental plastic strain ∆p and equivalent stress σeq from
the numerical solution is used as input in Eq. (45) to obtain an analytical solution for the temperature
evolution.

1D bar element

Figure 18 shows the performance of the model in describing the thermoviscoplastic behaviour (Figure 18a)
and temperature evolution (Figure 18b) in uniaxial tension for a bar element (FUN3 ).

(a) Equivalent stress σeq

(b) Temperature T

Figure 18: Results from temperature dependent 1D single element (FUN3 ) test in uniaxial tension (m = 1.0
and DC (cid:29)

1). Numerical results in EPX compared to the analytical solution in Eq. (41) and Eq. (45).

2D continuum element

Figure 19 shows the capabilities of the computational model in predicting the thermoviscoplastic response
at elevated strain rates for the Q42L element in uniaxial tension (Figure 19a) and simple shear (Figure 19b).

35

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.520100200300400500600700800Plastic strain p [mm/mm]Temperature [K]  T(p)(˙p≤5e−4,EPX)T(p)(˙p=1,EPX)T(p)(˙p=10,EPX)T(p)(˙p=100,EPX)T(p)(˙p=1000,EPX)T(p)(analytical)(a) Equivalent stress σeq (uniaxial tension)

(b) Equivalent stress σeq (simple shear)

(c) Temperature T (uniaxial tension)

(d) Temperature T (simple shear)

1) 2D single (continuum) element
Figure 19: Results from temperature dependent (m = 1.0 and DC (cid:29)
(Q42L) tests in uniaxial tension (left) and simple shear (right). Numerical results in EPX compared to the
analytical solution in Eq. (41) and Eq. (45).

3D shell element

The same tests are performed for the Q4GS element in uniaxial tension (Figure 20a) and simple shear
(Figure 20b).

36

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.520100200300400500600700800Plastic strain p [mm/mm]Temperature [K]  T(p)(˙p≤5e−4,EPX)T(p)(˙p=1,EPX)T(p)(˙p=10,EPX)T(p)(˙p=100,EPX)T(p)(˙p=1000,EPX)T(p)(analytical)00.511.520100200300400500600700800Plastic strain p [mm/mm]Temperature [K]  T(p)(˙p≤5e−4,EPX)T(p)(˙p=1,EPX)T(p)(˙p=10,EPX)T(p)(˙p=100,EPX)T(p)(˙p=1000,EPX)T(p)(analytical)(a) Equivalent stress σeq (uniaxial tension)

(b) Equivalent stress σeq (simple shear)

(c) Temperature T (uniaxial tension)

(d) Temperature T (simple shear)

1) 3D single (shell) element (Q4GS )
Figure 20: Results from temperature dependent (m = 1.0 and DC (cid:29)
tests in uniaxial tension (left) and simple shear (right). Numerical results in EPX compared to the analytical
solution in Eq. (41) and Eq. (45).

3D continuum element

Finally, the thermoviscoplastic response of the material model is veriﬁed for the CUBE element in uniaxial
tension (Figure 21a) and simple shear (Figure 21b).

37

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.520100200300400500600700800Plastic strain p [mm/mm]Temperature [K]  T(p)(˙p≤5e−4,EPX)T(p)(˙p=1,EPX)T(p)(˙p=10,EPX)T(p)(˙p=100,EPX)T(p)(˙p=1000,EPX)T(p)(analytical)00.511.520100200300400500600700800Plastic strain p [mm/mm]Temperature [K]  T(p)(˙p≤5e−4,EPX)T(p)(˙p=1,EPX)T(p)(˙p=10,EPX)T(p)(˙p=100,EPX)T(p)(˙p=1000,EPX)T(p)(analytical)(a) Equivalent stress σeq (uniaxial tension)

(b) Equivalent stress σeq (simple shear)

(c) Temperature T (uniaxial tension)

(d) Temperature T (simple shear)

1) 3D single (continuum) element
Figure 21: Results from temperature dependent (m = 1.0 and DC (cid:29)
(CUBE ) tests in uniaxial tension (left) and simple shear (right). Numerical results in EPX compared to the
analytical solution in Eq. (41) and Eq. (45).

Again, it is seen that the equivalent stress σeq from the numerical tests is in excellent agreement with the
analytical solution in Eq. (41) for all stress states (Figures 18 - 21), and that the computational model is
capable of describing the temperature softening of the equivalent stress due to adiabatic heating. This is also
supported by comparing the equivalent stress σeq and the temperature T at the end of the simulations (p =
2.0) in Table 5. Note that the analytical reference in Table 5 is found using the plastic strain and temperature
history from the uniaxial tests with the FUN3 element, while the values given for ∆σeq and ∆T (p = 2.0)
for each element and numerical test are found using the corresponding histories from the respective tests.
That is, the analytical reference used to compute ∆σeq and ∆T (p = 2.0) in Table 5 is found using the plastic
strain and temperature history from the respective test. Also note that in a thermoviscoplastic model there
is always a competition between viscoplasticity and thermal softening, where (viscoplastic) strain hardening
implies an increase of strength locally in the material and distributes the plasticity while thermal softening
results in localization of plasticity.

38

00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.52020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)σeq(analytical)00.511.520100200300400500600700800Plastic strain p [mm/mm]Temperature [K]  T(p)(˙p≤5e−4,EPX)T(p)(˙p=1,EPX)T(p)(˙p=10,EPX)T(p)(˙p=100,EPX)T(p)(˙p=1000,EPX)T(p)(analytical)00.511.520100200300400500600700800Plastic strain p [mm/mm]Temperature [K]  T(p)(˙p≤5e−4,EPX)T(p)(˙p=1,EPX)T(p)(˙p=10,EPX)T(p)(˙p=100,EPX)T(p)(˙p=1000,EPX)T(p)(analytical)Table 5: Summary of single element tests for thermo-viscoplastic behaviour behaviour at the end of the
simulations (p = 2.0).

Plastic strain rate
[s−

1]

4

5e−

≤

1

10

100

1000

T (p = 2.0) ∆T (p = 2.0)

[K]
711.3
711.2
710.5
711.1
710.5
711.1
710.9
711.1
738.4
738.3
738.3
738.9
738.2
738.9
738.3
738.9
747.0
746.9
746.9
747.4
746.9
747.4
746.9
747.5
755.8
755.6
755.6
756.3
755.6
756.3
755.6
756.3
764.5
764.3
764.3
764.9
764.2
764.9
764.3
764.9

[%]
-
-0.02
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
-
-0.02
-0.02
-0.02
-0.02
-0.02
-0.02
-0.02
-
-0.01
-0.01
-0.02
-0.01
-0.02
-0.01
-0.02
-
-0.01
-0.01
-0.02
-0.01
-0.02
-0.01
-0.02
-
-0.02
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01

Element

Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)
Analytical
FUN3 (uniaxial)
Q42L (uniaxial)
Q42L (shear)
Q4GS (uniaxial)
Q4GS (shear)
CUBE (uniaxial)
CUBE (shear)

∆σeq
[%]
-
0.10
-0.07
-0.05
-0.07
-0.05
0.06
-0.04
-
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-
-0.01
-0.01
-0.01
0.0
-0.01
0.0
0.0
-
-0.06
-0.07
-0.05
-0.07
-0.05
-0.05
-0.04

σeq
[MPa]
733.0
733.7
732.9
732.6
732.9
732.6
733.5
732.7
771.1
771.1
771.1
770.7
771.1
770.7
771.1
770.7
782.7
782.7
782.7
782.3
782.7
782.3
782.7
782.3
794.3
794.2
794.2
793.8
794.3
793.8
794.3
793.8
806.0
805.5
805.5
805.2
805.6
805.2
805.6
805.2

39

6.1.4. Ductile failure
Finally, ductile failure is included in the model by specifying Dc = 1 in the material input (see Table 1).
To couple the model with the element erosion option available in EPX, it is also necessary to include EROS
in the input ﬁle under the directive deﬁning the type of computation (typically CPLA LAGR EROS 1.0 or
TRID LAGR EROS 1.0 ) [7]. Then, when the damage variable in Eq. (46) reaches its critical value D = DC,
the material (or integration) point stiﬀness is deleted and the element is removed. Note that the respective
material (or integration) point may fail either due to the CL criterion when Dc = 1.0 or when the material
melts at T = Tm. The latter scenario is not included in these single element tests since the material melts
at very large strains for the material used in this study. However, the user should be aware that this type
of material failure may occur in problems with suﬃciently large strains.

1D bar element

Figure 22 shows the performance of the VPJC model in predicting ductile failure during the rate-dependent
behaviour in uniaxial tension for a bar element (FUN3 ).

Figure 22: Results from temperature and strain-rate dependent 1D single element (FUN3 ) test in uniaxial
tension including failure (m = 1.0 and DC = 1).

2D continuum element

Figure 23 shows the capabilities of the model in predicting ductile failure during the rate-dependent behaviour
at elevated strain rates for the Q42L element in uniaxial tension (Figure 23a) and simple shear (Figure 23b).

40

00.10.20.30.40.50.60.7020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)Failure(a) Uniaxial tension (Q42L)

(b) Simple shear (Q42L)

Figure 23: Results from temperature and strain-rate dependent 2D single element (Q42L) tests including
failure (m = 1.0 and DC = 1).

3D shell element

The same tests are performed for the Q4GS element in uniaxial tension (Figure 24a) and simple shear
(Figure 24b).

(a) Uniaxial tension (Q4GS )

(b) Simple shear (Q4GS )

Figure 24: Results from temperature and strain-rate dependent 3D single (shell) element (Q4GS ) tests
including failure (m = 1.0 and DC = 1).

3D continuum element

Finally, the models capability in predicting ductile failure is veriﬁed for the CUBE element in uniaxial
tension (Figure 25a) and simple shear (Figure 25b).

41

00.10.20.30.40.50.60.7020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)Failure00.10.20.30.40.50.60.70.80.911.1020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)Failure00.10.20.30.40.50.60.7020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)Failure00.10.20.30.40.50.60.70.80.911.1020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)Failure(a) Uniaxial tension (CUBE )

(b) Simple shear (CUBE )

Figure 25: Results from temperature and strain-rate dependent 3D single (continuum) element (CUBE )
tests including failure (m = 1.0 and DC = 1).

Figures 22 - 25 show that the VPJC material model accounts for viscoplasticity when the CL criterion is
used to predict ductile failure. That is, elevated strain rates result in reduced elongation to failure due to
the viscoplastic increase in the equivalent stress. This is explained by Eq. (20) where an increased strain
rate results in an increased ﬂow stress σy(σ0, p, ˙p). Consequently, the equivalent stress will increase to fulﬁll
the yield criterion (fv(σ) = 0). This also results in an increase in the major principal stress σ1 and that the
damage parameter Wc in Eq. (46) is reached at smaller elongations. The diﬀerence in elongation at failure
pf in uniaxial tension and simple shear is mainly explained by the dependence of stress triaxiality σ∗ and
Lode parameter µσ in Eqs. (46)-(51). Note that the diﬀerences observed in failure strain in Figures 23 - 25
corresponds well with the corresponding diﬀerence predicted by Eq. (51) in Figure 6.

6.2. Tension test

The material parameters (Q1, C1, Q2 and C2) in the ﬁrst multiplicative term in Eq. (41) and CL parameter
Wc in Eq. (46) were determined based on inverse modelling of the quasi-static tension tests presented in
[21]. These tension tests therefore serve as an excellent basis for validation of the quasi-static part of the
VPJC material model.

The uniaxial tension tests were modelled based on the experimental geometry (see Figure 26a) and pre-
scribing the same boundary conditions as in the material tests (see Figure 26b). The numerical model was
discretized in 3116 4-node quadrilaterals where the initial element size in the gauge area was equal to the
thickness (0.8 mm) in an attempt to capture the diﬀuse and local necking observed in the experiments.
Two simulations were performed using both the Q42L and Q4GS elements, where the former is a plane
stress elements with 2 dofs per node and 4 integration points (for full integration) while the latter is a 3D
shell element with 6 dofs per node and 20 integration points (5 through the thickness) [7]. Mass-scaling
by increasing the density by a factor 109 (ρ = 7850E9) was used to speed up the simulation, because an
explicit solution in time results in signiﬁcant computational times (i.e., a considerable amount of time in-
crements) for quasi-static tests. The quasi-static nature of the test was ensured by comparing the kinetic
energy present in the simulations to the plastic dissipation and found negligible. Moreover, since the plastic
response of the material cannot be assumed to occur within a short time duration in a quasi-static test, the
assumption of an adiabatic condition is no longer valid and the simulation should be run without thermal
softening (m = 0). It was also assumed negligible viscoplasticity in these tests (C = 0), while the remaining
material parameters and physical constants were taken from Table 1. The reader is referred to [21, 27] and

42

00.10.20.30.40.50.60.7020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)Failure00.10.20.30.40.50.60.70.80.911.1020040060080010001200Plastic strain p [mm/mm]Stress [MPa]  σeq(˙p≤5e−4,EPX)σeq(˙p=1,EPX)σeq(˙p=10,EPX)σeq(˙p=100,EPX)σeq(˙p=1000,EPX)Failure[29] for more information about the experimental work and numerical model, respectively, for these tension
tests.

(a) Geometry (experimental)

(b) Numerical model

Figure 26: Geometry and numerical model of the uniaxial tension test.

A comparison of the experimental and numerical results in terms of nominal stress-strain curves is presented
in Figure 27. It is seen that the numerical simulations are in excellent agreement with the experimental data,
and the quasi-static behaviour of the VPJC material model may be assumed to be properly implemented
in EPX.

Figure 27: Nominal stress-strain curves from uniaxial tension tests along three diﬀerent loading directions
for the Docol 600 DL steel. Numerical results from EPX (FEA) for both Q42L and Q4GS elements are
compared to the experimental data in [21]. The red and green dots denote the point of failure in the
determination of Wc.

43

7020012.540φ12.7R1525.0550.05yxvx=2.1mm/minNodalconstraintsNodalloads7000.050.10.150.20.250.30100200300400500600700800Nominal strain [mm/mm]Nominal stress [MPa]  Experiment 0°Experiment 45°Experiment 90°FEA (Q42L)FEA (Q4GS)Wc = 473 MPaWc = 473 MPa6.3. Airblast experiments

A series of airblast experiments on the same steel (Docol 600 DL) material were further used as an validation
of the dynamic part of the VPJC material model. Note that the material tests by Holmen et al. [21] were
taken from the same plates (i.e., the same batch) as those used in this study. Thus, the material data
provided in Table 1 also apply for these airblast experiments.

×

×

0.4

1.0

0.015 m3 and a square opening of 0.3

The experimental setup is shown in Figure 28, and consists of a steel mounting frame ﬁxed to the concrete
0.3 m2 at the centre. The
ﬂoor with outer dimensions 1.0
×
0.0008 m3 were clamped to the rigid mounting frame
square plate specimen with dimensions of 0.4
×
using M12 bolts and a clamping frame in an attempt to achieve ﬁxed boundary conditions. The explosive
mass W was positioned at stand-oﬀ distances R = 0.125 m, 0.250 m and 0.375 m relative to the centre
point of the plate. The explosive material was Composition C-4 with a spherical charge, a mass of 30 g
(equivalent to 40.2 g of TNT) and a diameter of approximately 34.5 mm (Figure ??). Finally, the plates
were painted with a speckle pattern (Figure 28c) and three-dimensional digital image correlation (3D-DIC)
analyses were conducted to measure the transient displacement ﬁeld during the tests. The stereovision setup
of the high-speed cameras (Phantom v1610) is illustrated in Figure 28a. The reader is referred to [27] for
further details regarding the experimental work.

×

(a) Sketch of setup

(b) Spherical charge

(c) DIC speckle pattern

Figure 28: Experimental setup airblast test [27]. The speckle pattern in (c) is seen from the cameras.

44

SteelmountingframeBlastpencil∼25◦TestspecimenClampingframeChargeof30gC4∼3300mmStand-oﬀ-distanceCam1Cam2Figure 29 shows the assembly of the numerical model, where the symmetry of the problem was utilized to
model only one quarter of the experimental setup using symmetric boundary conditions. The plate was
modelled using a Lagrangian discretization with an element size of approximately 2.5 mm (Figure 29b)
and Q4GS shell elements, while the bolts and clamping frames were represented by 8-node brick elements
(CUB8 ) with 8 integration points and the VPJC model with a high elastic limit to ensure elastic behaviour
using the physical constants for steel in Table 1. The steel mounting frame and the bolts were modelled as
one component (Figure 29a) and the bolts were modelled as stress-free. The pre-tensioning of the bolts was
accounted for by applying an external pressure at the contact area (Ac = 175 mm2) between the bolt heads
and the clamping frame (see Figure 29c). Dividing the pre-tensioning force Fp in each bolt by this contact
area resulted in a contact pressure of 527 MPa between each bolt head and the clamping frame. Contact
between the plate, bolts and frames was modelled using a surface-to-surface contact algorithm (GLIS ). This
contact algorithm is based on slave nodes and master surfaces where contact was enforced by Lagrangian
multipliers when a slave node penetrated a master surface. The plate was modelled as the slave and the
static friction coeﬃcient between the plate and clamping frames was set to 0.15, while the dynamic friction
coeﬃcient was equal to 0.10. The blast loading was applied on the plate by using the AIRB last directive in
EPX where the mass and positioning of the charge were used as input [7]. A more detailed description of
the numerical model is given in [29].

(a) Steel mounting frame and bolts

(b) Plate specimen

(c) Complete clamping assembly

Figure 29: Numerical model of the airblast experiments showing (a) steel mounting frame and bolts as one
component (in cyan), (b) plate specimen (in green) and (c) complete assembly including the clamping frame
(also in cyan) and contact area between bolt heads and clamping frame (in magenta) used to model the
eﬀect of the pre-tensioning of the bolts [29].

A comparison of the experimental and numerical results in terms of mid-point deﬂection versus time is pre-
sented in Figure 30. Considering the complexity in these airblast tests, it is observed a very good agreement
between the numerical simulations and experimental observations. In particular, the maximum displacement
corresponds very well with the experimental data for the two largest stand-oﬀ distances and the reversed
snap buckling is captured at the largest stand-oﬀ distance. The representation of the maximum displacement
is also quite good at the closest stand-oﬀ distance, however, the displacement is somewhat underestimated.
This is probably due to poorer representation of the actual loading on the plate, since we are outside the
range of the empirical methods used in the AIRB directive (Z = R/W 1/3 < 0.4 m/kg1/3). However, since
the scope of this report is to validate the implementation of the VPJC model this is not given any more
consideration. Thus, the VPJC material model is found to perform excellent also in fast transient dynamics.

45

Figure 30: Comparison of experimental and numerical results (FEA) for the airblast tests [29].

7. Concluding Remarks

An elastic-thermoviscoplastic Johnson-Cook type constitutive model was implemented in EUROPLEXUS
(EPX) and called VPJC. The VPJC model allows the user to choose between either a fully vectorized ver-
sion of the forward Euler integration (i.e., cutting plane) algorithm or a special case of the backward Euler
(i.e., radial return) algorithm. Ductile failure is also included in the material model through the criterion
proposed by Cockcroft and Latham. This allows the model to be coupled with the element erosion option
available in EPX. The material model is applicable for applicable for one-, two- and three-dimensional stress
states and a wide range of elements in EPX.

The performance of the model was ﬁrst veriﬁed using single element simulations in uniaxial tension and
simple shear. The model was found to respond very well to all the relevant test cases in terms of rate-
independent, viscoplastic and thermoviscoplastic behaviour. The dependence of the stress triaxiality on the
failure strain was also shown by including ductile failure in the single element veriﬁcation. Finally, to ensure
a realistic behaviour, the VPJC model was validated against experimental data from quasi-static tension
tests and airblast experiments. It was found that the numerical simulations were in excellent agreement
with the experimental data, and the VPJC model is therefore found to be properly implemented in EPX.

8. Further work

1. The hardening (Power law) given in Eq. (40) could be included in the model. By introducing a FLAG
in the material card the user could then choose between Voce hardening (e.g. FLAG = 0) or Power
law hardening (e.g. FLAG = 1). This may be done in the plastic corrector by

IF F LAG = 0 THEN use Voce hardening from Eq. (42)
C1p) + Q2C2 exp(

HR = Q1C1 exp(

C2p)

−

−

ELSEIF F LAG = 1 THEN use Power law hardening from Eq. (40)

HR = Bnpn

−

1

ELSE Material routine stops and print message:
’VPJC - INVALID TYPE OF HARDENING, FLAG=0 or 1’

It should be noted that the material card must be extended to include three more parameters, i.e., B,
n and a FLAG.

46

0246810−15−10−505101520253035Time [ms]Mid−point deflection [mm]  ExperimentsFEA (R = 0.125 m)FEA (R = 0.250 m)FEA (R = 0.375 m)2. Extend the radial return method to include plane stress. This could be done by constraining the
original three-dimensional constitutive equation to obtain a plane stress state. That is, prescribe only
the in-plane strains in the three-dimensional equation and use plane stress constraints as an additional
condition to determine the in-plane stresses and out-of-plane strains. The reader is referred to [5]
(Chapter 9) and [26] for more information on this approach.

3. A substepping scheme could be implemented to increase the stability and accuracy of the return map

algorithms.

4. Implement and add an option of using the Johnson-Cook failure criterion [30]. This is done by adding

the following during the update of variables depending on DLAMBDA

C COMPUTE FAILURE STRAIN BASED ON JOHNSON-COOK
D3 ∗
((1 + DLAM BDA/(P DOT
(1 + D5 ∗

P F = (D1 + D2 ∗
.
∗

SIGM AH/P HI))

∗
T R)))

T R)/(T M

EXP (

DT ))

D4)

((T

∗ ∗

−

−

−

C COMPUTE AND UPDATE DAMAGE VARIABLE

∗

.

DJC = DJC + DLAM BDA/P F

∗
...continue as is ...

∗
990

IF (D >= DC OR DJC >= DC) THEN

the current gp was already failed or is failing now ...

∗

∗
...continue as is ...

IF (D >= DC) THEN

WRITE(59,*) ’Element deleted due to CL-criterion’

IF (DJC >= DC) THEN

WRITE(59,*) ’Element deleted due to JC-criterion’

...continue as is ...

Note that the material input must be extended to include ﬁve more parameters, i.e., D1, D2, D3, D4
and D5. Material constants for a Weldox steel can be found in [14], i.e., D1 = 0.0705, D2 = 1.732,
D3 =

0.015 and D5 = 0.

0.54, D4 =

−

−

References

[1] T. Belytschko, W. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons (2000).
[2] M. Ortiz, J. Simo, An analysis of a new class of integration algorithms for elastoplastic relations, International Journal of

Numerical Methods in Engineering 23 (1986) 353–366.

[3] J. Simo, R. L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Computer Methods in Applied

Mechanics and Engineering 48 (3) (1985) 101–118.

[4] R. D. Krieg, D. B. Krieg, Accuracies of numerical solution methods for elastic-perfectly plastic model, ASME Journal of

Pressure Vessel Technology 99 (4) (1977) 510–515.

[5] E. A. de Souza Neto, D. Peric, D. R. J. Owen, Computational methods for plasticity : theory and applications, Wiley,

UK (2008).

[6] G. R. Johnson, W. H. Cook, A Constitutive Model and Data for Metals Subjected to Large Strains, High Strain Rates

and High Temperatures, Proceedings of the 7th International Symposium on Ballistics, Hague (1983) 541–547.

[7] Europlexus, User’s manual, Joint Research Centre, http://europlexus.jrc.ec.europa.eu/public/manual html/index.html

[cited 07.06.16].

47

[8] S. A. Khan, S. Huang, Continuum Theory of Plasticity, Wiley-Interscience, 1995.
[9] M. G. Cockcroft, D. J. Latham, Ductility and the workability of metals, Journal of the institute of metals 96 (1968) 33–39.
[10] O. S. Hopperstad, T. Børvik, Lecture Notes KT8306 Plasticity Theory, Department of Structural Engineering, SIMLab

(2013).

[11] J. Lubliner, Plasticity theory, Dover (2008).
[12] F. Irgens, Continuum Mechanics, Springer-Verlag (2008).
[13] M. Ortiz, G. T. Camacho, Adaptive Lagrangian modelling of ballistic penetration metallic targets, Comput. Methods

Appl. Mech. Engrg. 142 (1997) 269–301.

[14] T. Børvik, O. S. Hopperstad, T. Berstad, M. Langseth, A computational model of viscoplasticity and ductile damage for

impact and penetration, European Journal of Mechanics - A/Solids 20 (2001) 685–712.

[15] E. Voce, The relationship between stress and strain for homogenous deformation, Journal of the Institute of Metals 74

(1948) 536–562.

[16] T. L. Anderson, Fracture Mechanics - Fundamentals and Applications, Taylor and Francis Group, third edition (2005).
[17] J. Lemaitre, J.-L. Chabocke, Mechanics of Solid Materials, Cambridge University Press.
[18] T. Børvik, S. Dey, A. H. Clausen, Perforation resistance of ﬁve diﬀerent high-strength plates subjected to small-arms

projectiles, International Journal of Impact Engineering 36 (2009) 948–964.

[19] A. Kane, T. Børvik, A. Benallal, O. S. Hopperstad, Failure criteria with unilateral conditions for simulation of plate

perforation, European Journal of Mechanics - A/Solids 30 (2011) 468–476.

[20] S. Dey, T. Børvik, O. S. Hopperstad, M. Langseth, On the inﬂuence of fracture criterion in projectile impact of steel

plates, Computational Materials Science 38 (2006) 176–191.

[21] J. K. Holmen, O. S. Hopperstad, T. Børvik, Low-velocity impact on multi-layered dual-phase steel plates, International

Journal of Impact Engineering 78 (2015) 161–177.

[22] J. K. Holmen, J. Johnsen, O. S. Hopperstad, T. Børvik, Inﬂuence of fragmentation on the capacity of aluminium alloy

plates subjected to ballistic impact, European Journal of Mechanics A/Solids 55 (2016) 221–233.

[23] M. L. Wilkins, Calculation of elastic-plastic ﬂow, Technical report UCRL-7322, Lawrence radiation Laboratory, University

of California. Livermore, California April 19.

[24] R. P. R. Cardoso, J. W. Yoon, Stress integration method for a nonlinear kinematic/isotropic hardening model and its

characterization based on polycrystal plasticity, International Journal of Plasticity 25 (2009) 1684–1710.

[25] R. de Borst, P. H. Feenstra, Studies in anisotropic plasticity with reference to the Hill criterion, International Journal for

Numerical Methods in Engineering 29 (1990) 315–336.

[26] J. C. Simo, R. L. Taylor, A Return Mapping Algorithm for Plane Stress Elastoplasticity, International Journal for Nu-

merical Methods in Engineering 22 (1986) 649–670.

[27] V. Aune, E. Fagerholt, K. Hauge, M. Langseth, T. Børvik, Experimental study on the response of aluminium and
steel plates subjected to free airblast loading, International Journal of Impact Engineering 90 (2016) 106–121. doi:
10.1016/j.ijimpeng.2015.11.017.

[28] F. Casadei, V. Aune, G. Valsamos, M. Larcher, Testing of the Johnson-Cook material model VPJC in EUROPLEXUS,
European Commission, Joint Research Centre, Technical Report, JRC 98848 EUR 27594 (OPOCE LB-NA-27594-EN-N)
(2015) 1–62. doi:10.2788/02760.

[29] V. Aune, G. Valsamos, F. Casadei, M. Larcher, M. Langseth, T. Børvik, Numerical study on the structural response of

blast-loaded thin aluminium and steel plates, Submitted for journal publication (2016).

[30] G. R. Johnson, W. H. Cook, Fracture characteristics of three metals subjected to various strains, strain rates, temperatures

and pressures, Engineering Fracture Mechanics 21 (I) (1985) 31–48.

48

Appendix A. Cutting plane algorithm by using a linear Taylor-expansion

A linear Taylor-expansion around fn+1(σ, X) = fn+1(σ, σy, υp( ˙p), Γp(Tn)) reads

f i+1
n+1 = f i

n+1 +

∂f
∂σ

(cid:12)
(cid:12)
(cid:12)
(cid:12)

i

n+1

∆σn+1 +

∂f
∂X

(cid:12)
(cid:12)
(cid:12)
(cid:12)

i

n+1

∆Xn+1

(A.1)

which may be written in terms of the internal variables X as

f i+1
n+1 = f i

∂f
∂σy,n+1
where the ﬂow stress σy,n+1 and the following expressions are chosen as the internal variables X in the
Taylor-expansion

∂f
∂υp,n+1

∂f
∂σn+1

δυp,n+1 = 0

δσy,n+1 +

δσn+1 +

n+1 +

(A.2)

(cid:18)

υp(∆λn+1) =

1 +

(cid:19)C

∆λn+1
˙p0∆tn+1

(cid:18)

Γp(Tn) =

1

−

Introducing the partial derivatives of the internal variables Eq. (A.2) reads

(cid:21)m(cid:19)

(cid:20) Tn −
Tm −

Tr
Tr

(A.3)

f i+1
n+1 = f i

n+1 +

∂f
∂σn+1

δσn+1
(cid:124) (cid:123)(cid:122) (cid:125)
∂σn+1

C ∂f

−

δλi+1
n+1

+

∂f
∂σy,n+1
(cid:124)
(cid:125)
(cid:123)(cid:122)
1)υp,n+1Γp,n

(
−

δσy,n+1
(cid:124) (cid:123)(cid:122) (cid:125)

∂σy,n+1
∂λn+1

δλi+1
n+1

+

(

∂f
∂υp,n+1
(cid:124) (cid:123)(cid:122) (cid:125)
1)σy,n+1Γp,n

−

δυp,n+1
(cid:124) (cid:123)(cid:122) (cid:125)

∂υp,n+1
∂ ˙pn+1

δ ˙pi+1
n+1

= 0

(A.4)

f i+1
n+1 = f i

n+1 −

∂f
∂σn+1

C

∂f
∂σn+1

δλi+1

n+1 −

υp,n+1Γp,n

∂σy,n+1
∂λn+1

δλi+1

n+1 −

σy,n+1Γp,n

∂υp,n+1
∂ ˙λn+1

δ ˙λi+1

n+1 = 0

(A.5)

f i+1
n+1 = f i

n+1 −

∂f
∂σn+1

C

∂f
∂σn+1

δλi+1

n+1 −

υp,n+1Γp,nHR,n+1δλi+1

n+1 −

σy,n+1Γp,n

∂υp,n+1
∂ ˙λn+1

δ ˙λi+1

n+1 = 0

(A.6)

Moreover, using ˙p = ˙λ = ∆λ
υp(∆λ) and Eq. (A.6) may be written as

∆t the internal variable controlling the viscoplasticity υp( ˙p) may be expressed as

f i+1
n+1 = f i

n+1 −

∂f
∂σn+1

C

∂f
∂σn+1

δλi+1

n+1 −

υp,n+1Γp,nHR,n+1δλi+1

n+1 −

σy,n+1Γp,n

∂υp,n+1
∂∆λn+1

δλi+1

n+1 = 0

(A.7)

Finally, the incremental change δλ in the incremental plastic multiplier ∆λ is then expressed as

δλi+1

n+1 =

∂f
∂σn+1

C ∂f

∂σn+1

+ υp,n+1Γp,nHR,n+1 + σy,n+1Γp,n

∂υp,n+1
∂∆λn+1

f i
n+1

where

∂υp(∆λn+1)
∂∆λ

=

C
˙p0∆t

(cid:18)

1 +

∆λn+1
˙p0∆t

(cid:19)C

1

−

HR,n+1 =

dRn+1
dpn+1

=

2
(cid:88)

k=1

QkCke(
−

Ckλn+1)

(cid:27)T

(cid:26) ∂f
∂σ

[C]

(cid:27)

(cid:26) ∂f
∂σ

49

(A.8)

(A.9)

(A.10)

(A.11)

[C] =

















C1111 C1122 C1133 C1112 C1121 C1123 C1132 C1113 C1131
C2211 C2222 C2233 C2212 C2221 C2223 C2232 C2213 C2231
C3311 C3322 C3333 C3312 C3321 C3323 C3332 C3313 C3331
C1211 C1222 C1233 C1212 C1221 C1223 C1232 C1213 C1231
C2111 C2122 C2133 C2112 C2121 C2123 C2132 C2113 C2131
C2311 C2322 C2333 C2312 C2321 C2323 C2332 C2313 C2331
C3211 C3222 C3233 C3212 C3221 C3223 C3232 C3213 C3231
C1311 C1322 C1333 C1312 C1321 C1323 C1332 C1313 C1331
C3111 C3122 C3133 C3112 C3121 C3123 C3132 C3113 C3131

















(cid:27)

(cid:26) ∂f
∂σ

=





































∂f
∂σ11
∂f
∂σ22
∂f
∂σ33
∂f
∂σ12
∂f
∂σ21
∂f
∂σ23
∂f
∂σ32
∂f
∂σ13
∂f
∂σ31

=

































f,11
f,22
f,33
f,12
f,21
f,23
f,32
f,13
f,31

(cid:27)T

(cid:26) ∂f
∂σ

[C] =















f,11C1111 + f,22C2211 + f,33C3311 + f,12C1211 + f,21C2111 + f,23C2311 + f,32C3211 + f,13C1311 + f,31C3111
f,11C1122 + f,22C2222 + f,33C3322 + f,12C1222 + f,21C2122 + f,23C2322 + f,32C3222 + f,13C1322 + f,31C3122
f,11C1133 + f,22C2233 + f,33C3333 + f,12C1233 + f,21C2133 + f,23C2333 + f,32C3233 + f,13C1333 + f,31C3133
f,11C1112 + f,22C2212 + f,33C3312 + f,12C1212 + f,21C2112 + f,23C2312 + f,32C3212 + f,13C1312 + f,31C3112
f,11C1121 + f,22C2221 + f,33C3321 + f,12C1221 + f,21C2121 + f,23C2321 + f,32C3221 + f,13C1321 + f,31C3121
f,11C1123 + f,22C2223 + f,33C3323 + f,12C1223 + f,21C2123 + f,23C2323 + f,32C3223 + f,13C1323 + f,31C3123
f,11C1132 + f,22C2232 + f,33C3332 + f,12C1232 + f,21C2132 + f,23C2332 + f,32C3232 + f,13C1332 + f,31C3132
f,11C1113 + f,22C2213 + f,33C3313 + f,12C1213 + f,21C2113 + f,23C2313 + f,32C3213 + f,13C1313 + f,31C3113
f,11C1131 + f,22C2231 + f,33C3331 + f,12C1231 + f,21C2131 + f,23C2331 + f,32C3231 + f,13C1331 + f,31C3131















(cid:27)T

(cid:26) ∂f
∂σ

[C]

(cid:27)

(cid:26) ∂f
∂σ

=

(f,11C1111 + f,22C2211 + f,33C3311 + f,12C1211 + f,21C2111 + f,23C2311 + f,32C3211 + f,13C1311 + f,31C3111)f,11
+ (f,11C1122 + f,22C2222 + f,33C3322 + f,12C1222 + f,21C2122 + f,23C2322 + f,32C3222 + f,13C1322 + f,31C3122)f,22
+ (f,11C1133 + f,22C2233 + f,33C3333 + f,12C1233 + f,21C2133 + f,23C2333 + f,32C3233 + f,13C1333 + f,31C3133)f,33
+ (f,11C1112 + f,22C2212 + f,33C3312 + f,12C1212 + f,21C2112 + f,23C2312 + f,32C3212 + f,13C1312 + f,31C3112)f,12
+ (f,11C1121 + f,22C2221 + f,33C3321 + f,12C1221 + f,21C2121 + f,23C2321 + f,32C3221 + f,13C1321 + f,31C3121)f,21
+ (f,11C1123 + f,22C2223 + f,33C3323 + f,12C1223 + f,21C2123 + f,23C2323 + f,32C3223 + f,13C1323 + f,31C3123)f,23
+ (f,11C1132 + f,22C2232 + f,33C3332 + f,12C1232 + f,21C2132 + f,23C2332 + f,32C3232 + f,13C1332 + f,31C3132)f,32
+ (f,11C1113 + f,22C2213 + f,33C3313 + f,12C1213 + f,21C2113 + f,23C2313 + f,32C3213 + f,13C1313 + f,31C3113)f,13
+ (f,11C1131 + f,22C2231 + f,33C3331 + f,12C1231 + f,21C2131 + f,23C2331 + f,32C3231 + f,13C1331 + f,31C3131)f,13

[C] =

















C1111 C1122 C1133
C2211 C2222 C2233
C3311 C3322 C3333
0
0
0
0
0
0

0
0
0
0
0
0

0
0
0
0
0
0

0
0
0

0
0
0

C1212 C1221
C2112 C2121

0
0
0
0
0

0
0
0
0
0

0
0
0
0
0
0
0

0
0
0
0
0
0
0

















C2323 C2332
C3223 C3232

0
0

0
0

C1313 C1331
C3113 C3131

0
0
0
0

0
0
0
0

50

(A.12)

(A.13)

(A.14)

(A.15)

(A.16)

Due to the symmetry of the stress, σij = σji, it is seen that also

Since the strain tensor also is symmetric, εkl = εlk, it is also assumed that

Cijkl = Cjikl

and the 9x9 matrix in Eq. (A.16) may be reduced to the following 6x6 matrix

Cijkl = Cijlk

[C] =











C1111 C1122 C1133
C2211 C2222 C2233
C3311 C3322 C3333
0
0
0

0
0
0

0
0
0

0
0
0
C1212
0
0

0
0
0
0
C2323
0











0
0
0
0
0
C1313

(A.17)

(A.18)

(A.19)

Adding major symmetry, i.e. Cijkl = Cklij, to Eq. (A.19) and assuming an isotropic material where all axes
are symmetry axes, i.e. c1 = C1111 = C2222 = C3333, c4 = C1122 = C2211 = C1133 = C3311 = C2233 = C3322
and c7 = C1212 = C2323 = C1313. Then, Eq. (A.15) reduces to

(cid:27)T

(cid:26) ∂f
∂σ

[C]

(cid:27)

(cid:26) ∂f
∂σ

=

(f,11C1111 + f,22C2211 + f,33C3311)f,11 + (f,11C1122 + f,22C2222 + f,33C3322)f,22
+ (f,11C1133 + f,22C2233 + f,33C3333)f,33
+ (f,12C1212 + f,21C2112)f,12 + (f,12C1221 + f,21C2121)f,21
+ (f,23C2323 + f,32C3223)f,23 + (f,23C2332 + f,32C3232)f,32
+ (f,13C1313 + f,31C3113)f,13 + (f,13C1331 + f,31C3131)f,13
= c1 (cid:0)f 2
(cid:1) + 2c4 (f,11f,22 + f,22f,33 + f,33f,11)
(cid:1)
+ 4c7 (cid:0)f 2

,22 + f 2
,33
,23 + f 2
,31

,11 + f 2
,12 + f 2

(A.20)

However, since the strain tensor is symmetric, εkl = εlk when k
coeﬃcients is reduced to 3

= l, it is seen that the number of independent











c1
c4
c4
0
0
0

[C] =











C1111 C1122 C1133
C2211 C2222 C2233
C3311 C3322 C3333
0
0
0

0
0
0

0
0
0

0
0
0
C1212
0
0

0
0
0
0
C2323
0











0
0
0
0
0
C1313

=

where it is used Voigt notation and symmetric stress tensors

σ4 = σ12 = σ21 = σ7
σ5 = σ23 = σ32 = σ8
σ6 = σ31 = σ13 = σ9

51

c4
c1
c4
0
0
0

c4
c4
c1
0
0
0

0
0
0
c7
0
0

0
0
0
0
c7
0











0
0
0
0
0
c7

(A.21)

(A.22)

(cid:54)
The symmetry of the stress and strain tensors allows us to express them as column matrices (or vectors) in
the form

[σ] =











and [ε] =











σ11
σ22
σ33
σ12
σ23
σ31





















ε11
ε22
ε33
γ12
γ23
γ31

where γ12 = 2ε12, γ23 = 2ε23 and γ31 = 2ε31 are the engineering shear strains.

[C] =











c1
c4
c4
0
0
0

c4
c1
c4
0
0
0

c4
c4
c1
0
0
0

0
0
0
c7
0
0

0
0
0
0
c7
0











0
0
0
0
0
c7

which implies that Eq. (A.15) may be written on a more compact form, i.e.

(A.23)

(A.24)

(cid:27)T

(cid:26) ∂f
∂σ

[C]

(cid:27)

(cid:26) ∂f
∂σ

=

∂f
∂σi

∂f
∂σj

Cij
(cid:34)(cid:18) ∂f
∂σ1
(cid:34)(cid:18) ∂f
∂σ4

= c1

+ 4c7

(cid:19)2

+

(cid:19)2

(cid:19)2

+

(cid:19)2

(cid:18) ∂f
∂σ2
(cid:18) ∂f
∂σ5

(cid:18) ∂f
∂σ3
(cid:18) ∂f
∂σ6

+

+

(cid:19)2(cid:35)

(cid:19)2(cid:35)

+ 2c4

(cid:20) ∂f
∂σ1

∂f
∂σ2

+

∂f
∂σ2

∂f
∂σ3

+

∂f
∂σ3

∂f
∂σ1

(cid:21)

(A.25)

52

Appendix B. Radial Return Algorithm

As discussed in Section 5.2, the system of equations in Eq.(84) can be extensively simpliﬁed and the return
mapping for the von Mises model, with associated ﬂow rule and isotropic hardening, was reduced to a single
nonlinear equation with the incremental plastic multiplier ∆λ as the only unknown. The objective of this
section is to give a theoretical understanding of the result ϕ = ϕtrial
3G∆λ, which leads to this special
case called the radial return algorithm. It is readily seen that the radial return algorithm reduce the number
of equations and is important to make the state-update procedure more computationally eﬃcient, and to
improve the performance of the overall ﬁnite element scheme. For simplicity only the isotropic strain hard-
ening is considered in the constitutive equations.

−

Starting with the general elastic predictor plastic corrector approach, the stress state at time tn+1 may be
expressed as

σn+1 = σn + C : (cid:0)∆εn+1 −

∆εp

n+1

(cid:1) = σtrial

n+1 −

C : ∆εp

n+1 = σtrial

n+1 −

∆ ˜σp

(B.1)

where σtrial
∆ ˜σp returns the stress to the yield surface (Figure B.31).

n+1 is the elastic predictor moving the stress outside the yield surface, and the plastic corrector

Introducing the associated ﬂow rule, Eq.(B.1) reads

σn+1 = σtrial

n+1 −

3
2

∆λn+1C :

(cid:48)
σ
n+1
σeq

Using the relation (Appendix C)

Eq.(B.2) is given by

(cid:48)

C : σ

n+1 = 2Gσ

(cid:48)
n+1

σn+1 = σtrial

n+1 −

3G

(cid:48)
σ
n+1
σeq

∆λn+1

(B.2)

(B.3)

(B.4)

In the case of a three-dimensional stress state, the von Mises yield criterion is a circular cylinder in the
principal stress space with the deviatoric plane normal to the hydrostatic axis (Figure 4). The von Mises
ﬂow vector ˙εp is purely deviatoric, and the hydrostatic stress σH,n+1 is unaﬀected by the plastic ﬂow
(σtrial
H,n+1 = σH,n+1). This implies that the stress update can be formulated more eﬃciently be splitting the
stress into deviatoric and hydrostatic parts, where only the deviatoric stress component need to be included
in the return mapping, i.e.

(cid:48)

σ

n+1 = σ

(cid:48)trial
n+1 −

3G

(cid:48)
σ
n+1
σeq

∆λn+1 = σ

(cid:48)trial
n+1 −

3G

σ
ϕ(σ(cid:48)

(cid:48)
n+1
n+1)

∆λn+1

Rearranging Eq.(B.5) results in an important result, i.e.

(cid:18)

1 +

σ

(cid:48)
n+1

(cid:19)

3G∆λn+1
ϕ(σ(cid:48)
n+1)

= σ

(cid:48)trial
n+1

which implies that the trial and updated deviatoric stresses are co-linear, i.e.

(cid:48)
n+1

σ
ϕ (cid:0)σ(cid:48)

n+1

σ

(cid:48)trial
n+1
ϕ (cid:0)σ(cid:48)trial

n+1

(cid:1)

(cid:1) =

53

(B.5)

(B.6)

(B.7)

This means that the gradient of the yield surface coincide at the trial and updated stress states (Figure
B.31). Substituting the result in Eq.(B.7) into (B.5) results in the following simpliﬁcation

(cid:48)

σ

n+1 = σ

(cid:48)trial
n+1

(cid:18)

1

−

3G∆λn+1
ϕ(σ(cid:48)trial
n+1 )

(cid:19)

(B.8)

where ϕ(σ

(cid:48)trial
n+1 ) =

(cid:113)

3J2

(cid:0)σ(cid:48)trial

n+1

(cid:1) is the elastic trial von Mises equivalent stress.

(cid:48)trial
n+1

is a constant vector in the return mapping, and that the deviatoric stress σ

(cid:48)
It is also noted that σ
n+1
is a function only dependent of ∆λ in Eq.(B.8). Furthermore, Eq.(B.8) implies that, in the radial return
algorithm with von Mises yielding, the updated deviatoric stress is obtained by scaling the trial deviatoric
(cid:48)trial
stress by a factor 1
n+1 ). The geometric representation of this update formula in the
−
deviatoric plane is illustrated in Figure B.31.

3G∆λn+1/ϕ(σ

σ3

∂f
∂σ

(cid:12)
(cid:12)
(cid:12)n+1

σ

(cid:48)
n+1

σ

(cid:48)trial
n+1

(cid:48)trial
+1
n
σ
)
(cid:48)trial
∆λn
1
+
3G
n
ϕ(σ

1
+

π-plane

(cid:48)
σ
n

σ1

surface at tn

surface at tn+1

σ2

Figure B.31: The implicit elastic predictor plastic corrector radial return map algorithm for von Mises
yielding in the deviatoric plane.

Finally, by inserting Eq.(B.8) into the yield criterion, the following nonlinear scalar equation is obtained

f (∆λ) = ϕ

= ϕ

= ϕ

= ϕ

(cid:16)

(cid:16)

(cid:16)

(cid:16)

(cid:17)

σ

(cid:48)
n+1

σ

(cid:48)trial
n+1

−
(cid:17) (cid:18)

(cid:17)

(cid:17)

σ

(cid:48)trial
n+1

σ

(cid:48)trial
n+1

−

−

σy,n+1

1

−

3G∆λn+1
ϕ(σ(cid:48)trial
n+1 )
3G∆λn+1 −
3G∆λn+1 −

σy,n+1

54

(cid:19)

σy,n+1

−

(B.9)

(σy,n + HR,n+1∆λn+1) = 0

(cid:48)

where the incremental plastic multiplier ∆λ is the only unknown. In Eq.(B.9) it is utilized that ϕ(aσ
aϕ(σ
Thus,

) =
) is assumed to be a positive homogeneous function of order one of the stress.

), since the function ϕ(σ

(cid:48)

(cid:48)

(cid:48)

ϕ(σ

n+1) = ϕ

(cid:18)

σ

(cid:48)trial
n+1

= ϕ(σ

(cid:48)trial
n+1 )

(cid:20)
1

(cid:18)

−

1

−

3G∆λn+1
ϕ(σ(cid:48)trial
n+1 )
3G∆λn+1
ϕ(σ(cid:48)trial
n+1 )

(cid:21)(cid:19)

(cid:19)

By using Eq.(B.9), the incremental plastic multiplier at time tn+1 is expressed as

∆λn+1 =

(cid:16)

ϕ

σ

(cid:48)trial
n+1

(cid:17)

σy,n

−

3G∆λn+1 + HR,n+1

which, in the case of nonlinear hardening, could be solved by a Newton-Raphson method

∆λi+1

n+1 = ∆λi

n+1 −

= ∆λi

n+1 −

f (∆λ)
f (cid:48)(∆λ)
(cid:16)

ϕ

σ

(cid:48)trial
n+1

(cid:17)

−

3G∆λi

n+1 −

σy,n+1

3G

HR

−

−

(B.10)

(B.11)

(B.12)

where i is the local iteration counter. For linear hardening ∆λ is found directly and no local iterations are
necessary.

Finally, the stresses and strains are updated according to Eq.(B.2)

(cid:1)

(B.13)

σn+1 = σtrial

n+1 −

= σtrial

n+1 −

3
2

3
2

∆λn+1C :

∆λn+1C :

pn+1 = pn + ∆λn+1

(cid:48)
σ
n+1
σeq
σ

(cid:48)trial
n+1
ϕ (cid:0)σ(cid:48)trial

n+1

εp
n+1 = εp

n +

3
2

σ

(cid:48)trial
n+1
ϕ (cid:0)σ(cid:48)trial

n+1

(cid:1) ∆λn+1

55

Appendix C. Isotropic elastic material

Assuming isotropic linear elastic behavior of the material reduces the stiﬀness tensor to only depend on the
Lam`e constants, i.e

and the stress is then given by the linear relation

Cijkl = λelδijδkl + µel (δikδjl + δilδjk)

σij = Cijklεkl = λelδijδklεkl + µel (δikδjlεkl + δilδjkεkl)
= λelδijεkk + µel (εij + εji)
= λelδijεkk + 2µel (εij)

where the following relations are used

δikδjlεkl = δikεjk = εij
δilδjkεkl = δilεlj = εji

εij = εji

By inserting the expression for the Lam`e constants, Eq.(C.1) reads

Cijkl =

=

νE
(1 + ν)(1
(cid:20) 1
2

E
1 + ν
(cid:20) 1
2

δijδkl +

2ν)

−
(δikδjl + δilδjk) +

E
2(1 + ν)
ν

(δikδjl + δilδjk)

(cid:21)

δijδkl

(cid:21)

2ν

−

1
ν

δijδkl

2ν

−

= 2G

(δikδjl + δilδjk) +

1

(C.1)

(C.2)

(C.3)

(C.4)

and by multiplying the isotropic elasticity tensor in Eq.(C.4) with the deviatoric part of the stress tensor,
this product reduces to

Cijklσ

(cid:48)

kl = 2G

= 2G

(cid:16)

(cid:16)

(cid:20) 1
2
(cid:20) 1
2

δikδjlσ

(cid:48)

kl + δilδjkσ

(cid:17)

(cid:48)
kl

+

1

(cid:48)

ij + σ

(cid:48)
ji

σ

(cid:17)

+

ν

2ν

1

−

δijσ

(cid:48)
kk

(cid:21)

(cid:48)
kl

δijδklσ

2ν

ν

−
(cid:21)

Utilizing the symmetry of the stress tensor, i.e. σ

(cid:48)

ij = σ

(cid:48)

ji, and the deﬁnition of the deviatoric stress

σ

(cid:48)

ij = σij −

1
3

σkkδij,

(cid:48)

ii = σii −
σ

1
3

σkkδii = 0

Eq.(C.5) reduces to the following expression

Cijklσ

(cid:48)

kl = 2Gσ

(cid:48)
ij

which is a useful relation in the radial return map algorithm.

56

(C.5)

(C.6)

(C.7)

Appendix D. Two-dimensional stress analysis

This section is an extened version of Section 3 and shows useful relationships when dealing with the particular
case of plane stress and three shear stress components. Plane stress is already deﬁned in Section 3 as the
state of stress in which the normal stress σ33 and the shear stresses σ13,23 are assumed to be zero [11]. From
a practical point of view this assumption is typically introduced in the analysis of bodies in which one of
the dimensions (i.e. the thickness) is much smaller than the others and is subjected to loads that result in
dominant stresses perpendicular to the thickness direction [5]. However, the same assumption may also be
used when only the normal stress σ33 is assumed to be zero. Thus, all the shear components are considered
in the calculation.

Appendix D.1. Two-dimensional theory of elasticity
Recall that, in the particular case of plane stress, the linear elastic stress-strain compliance relationship for
an isotropic material can be written as a system of equations in matrix form given as











= Sσ =

1
E











ε11
ε22
ε33
2ε12
2ε23
2ε31











1
ν
−
ν
−
0
0
0

ν
−
1
ν
−
0
0
0

ν
−
ν
−
1
0
0
0

0
0
0
2(1 + ν)
0
0

0
0
0
0
2(1 + ν)
0

0
0
0
0
0
2(1 + ν)





















σ11
σ22
0
σ12
σ23
σ31











(D.1)

The normal stress σ33 which is assumed to be zero in the stress vector indicates that its associated column
in the compliance matrix S can be neglected (i.e. column 3). Furthermore, if we neglect the row associated
with the strain components in the thickness direction εi3 (i = 1, 2, 3), the stress-strain compliance relation-
ship reduces to a 5

5 matrix

×









= Sσ =

1
E









ε11
ε22
2ε12
2ε23
2ε31









1
ν
−
0
0
0

ν
−
1
0
0
0

0
0
2(1 + ν)
0
0

0
0
0
2(1 + ν)
0









0
0
0
0
2(1 + ν)









σ11
σ22
σ12
σ23
σ31









(D.2)

Then, the elasticity matrix C for plane stress is found by inverting the plane stress compliance matrix S
and given by Hooke’s law as









σ11
σ22
σ12
σ23
σ31









= Cε =

E

ν2

1

−



1
ν


0


0

0

ν
1
0
0
0

0
0

ν

1

−
2
0
0

0
0
0

ν

1

−
2
0









0
0
0
0

ν

1

−
2









ε11
ε22
2ε12
2ε23
2ε31









(D.3)

It is important to note that the elasticity matrix for plane stress is not found by removing columns and rows
from the general isotropic elasticity matrix in Eq. (C.1).

It should also be noted that the procedure in ﬁnding the principal stresses σi of the stress tensor σ are now
determined by using Eq. (57), i.e.

(cid:48)

(cid:48)

σi = σ

i + σH = σ

i +

1
3

tr(σ)

(D.4)

where the major principal stress σ1 is used in the Cockcroft-Latham failure criterion in Eq. (46). Thus, the
procedure in Section 3.1 is not applicable when using all the shear components.

57

Appendix E. Implementation of Return Map Algorithm in EUROPLEXUS

This appendix shows a somewhat more detailed version of the numerical scheme given in Section 5.1. Within
the time increment ∆tn+1 the scheme may then be explained as follows:

1. Set the initial values of internal variables to the converged values of the previous step at tn, and
calculate the speed of sound. Note that the number of shear stress components (ITAU ) and plane
stress conditions (IPLANC ) (provided by the element routine) are used to SELECT the appropriate
stress state in the material routine.

Governing equations

Implementation in EUROPLEXUS

Deﬁne useful parameters related to the theory of
elasticity

SUBROUTINE VPJC (ITAU, IPLANC, PROPS, STRAININC, STRESSOLD,
>
STRAIN, STATE, STRESSNEW, CSON, FAIL_FLAG)
* sound speed

T 1 = 1 + ν, T 2 = 1

2ν, T 3 = 1

ν

−

−

the

Compute
stress state, i.e.,
IF 3D stress state THEN

speed of

(cid:115)

CSON =

E(1
−
ρ(1 + ν)(1

2ν)

−
ELSEIF plane stress THEN

sound based on the

ν)

T1 = R1 + POISS
T2 = R1 - R2*POISS
T3 = R1 - POISS

SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

CSON = SQRT (YOUNG*T3 / (RHO*T1*T2))

CASE (1)

! PLANE STRESS CONDITION (SIG(3)=0)

CSON = SQRT (YOUNG / (RHO*T1*T3))

CASE (2)

! UNIAXIAL STRESS CONDITION (SIG(2)=SIG(3)=0)

CSON = SQRT (YOUNG / RHO)

CASE DEFAULT ! CHECK THIS ONLY HERE

CALL ERRMSS (’VPJC’, ’WRONG IPLANC’)
STOP ’VPJC - WRONG IPLANC’

END SELECT

*

*

CSON =

(cid:115)

E

−

ν2)

ρ(1

ELSEIF uniaxial stress THEN

CSON =

(cid:115)

E
ρ

2. Compute the elasticity tensor.

Goverening equations

Implementation in EUROPLEXUS

IF a 3D stress state is used THEN use Eqs. (3)-(4)

C = λel I ⊗ I + 2µel I =










C1
C5
C5
0
0
0

C4
C2
C6
0
0
0

C5
C4
C3
0
0
0

0
0
0
C7
0
0

0
0
0
0
C8
0

C1 = C2 = C3 =

(1 − ν)E

(1 + ν)(1 − 2ν)

, C4 = C5 = C6 =

C7 = C8 = C9 =

1

E

2

1 + ν










,

0
0
0
0
0
C9
νC1
1 − ν

ELSEIF a 2D stress state is used THEN use

Eq. (60)

C =





C1
C5
0

C4
C2
0





0
0
C7

C1 = C2 =

E

1 − ν2

, C4 = C5 = νC1 , C7 =

1

E

2

1 + ν

ELSE a 1D stress state is used THEN use

C =





C1
0
0

0
0
0



0
0
0

 , C1 = E

* compute the elasticity tensor
SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

C1 = YOUNG*T3 / (T1*T2)
C4 = POISS*C1 / T3

CASE (1)

! PLANE STRESS CONDITION (SIG(3)=0)

C1 = YOUNG / (T1*T3)
C4 = POISS*C1

CASE (2)

! UNIAXIAL STRESS CONDITION (SIG(2)=SIG(3)=0)

C1 = YOUNG
C4 = R0
END SELECT

! THESE VALUES HAVE TO BE CHECKED

*

*

C2 = C1
C3 = C1
C5 = C4
C6 = C4
C7 = R5*YOUNG / T1
C8 = C7
C9 = C7

58

3. & 4. Use the total stain increment ∆εn+1 from the global equilibrium to compute the elastic predictor.
Then, compute the von Mises equivalent stress in terms of the elastic trial stress and Eq. (35). Note
that Eqs.(35) and (34) are used to calculate the gradients, i.e., σ12 = σ21, σ23 = σ32 and σ13 = σ31.

Goverening equations

Implementation in EUROPLEXUS

* compute the elastic trial stress and the trial equivalent stress
*

STRESS(1) = STRESSOLD(1) + C1*STRAININC(1) + C4*STRAININC(2)

>

*

+ C5*STRAININC(3)

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

STRESS(2) = STRESSOLD(2) + C4*STRAININC(1) + C2*STRAININC(2)

>

>

>
>

+ C6*STRAININC(3)

STRESS(3) = STRESSOLD(3) + C5*STRAININC(1) + C6*STRAININC(2)

+ C3*STRAININC(3)

PHITRIAL = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)
+ STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
- STRESS(2)*STRESS(3) - STRESS(3)*STRESS(1)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

STRESS(2) = STRESSOLD(2) + C4*STRAININC(1) + C2*STRAININC(2)

>

>

+ C6*STRAININC(3)

STRESS(3) = R0
PHITRIAL = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)

- STRESS(1)*STRESS(2)

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

STRESS(2) = R0
STRESS(3) = R0
PHITRIAL = STRESS(1)*STRESS(1)

END SELECT

*

*

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)
! ONLY TAU_XY (STRESS(4)) IS PRESENT
CASE (1)

STRESS(4) = STRESSOLD(4) + C7*STRAININC(4) ! THE GAMMAS ARE USED!
PHITRIAL = PHITRIAL + R3*(STRESS(4)*STRESS(4))

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

STRESS(4) = STRESSOLD(4) + C7*STRAININC(4) ! SHEAR STRAINS
STRESS(5) = STRESSOLD(5) + C8*STRAININC(5) ! ARE THE GAMMAS, NOT
STRESS(6) = STRESSOLD(6) + C9*STRAININC(6) ! THE EPSILONS
PHITRIAL = PHITRIAL + R3*(STRESS(4)*STRESS(4)
+ STRESS(5)*STRESS(5)
+ STRESS(6)*STRESS(6))

>
>

CASE DEFAULT ! WE CHECK THIS ONLY IN THE FIRST SELECT CASE ON ITAU

CALL ERRMSS (’VPJC’, ’WRONG ITAU’)
STOP ’VPJC - WRONG ITAU’

END SELECT

PHITRIAL = SQRT (PHITRIAL)

σ

trial
11,n+1 = σ11,n + C1∆ε11,n+1 + C4∆ε22,n+1

+ C5∆ε33,n+1

IF a 3D stress state is used THEN

σ

trial
22,n+1 = σ22,n + C4∆ε11,n+1 + C2∆ε22,n+1

+ C6∆ε33,n+1

σ

trial
33,n+1 = σ33,n + C5∆ε11,n+1 + C6∆ε22,n+1

+ C3∆ε33,n+1

σ

σ

σ

trial
12,n+1 = σ12,n + C7∆ε12,n+1
trial
23,n+1 = σ23,n + C8∆ε23,n+1
trial
31,n+1 = σ31,n + C9∆ε31,n+1

(cid:48) trial
n+1 : σ

(cid:48) trial
n+1 = (σ

σ

2
trial
11,n+1)

+ (σ

2
trial
22,n+1)

3

2

2 − σ

trial
11,n+1σ

trial
22,n+1

trial
33,n+1 − σ

trial
33,n+1σ

trial
11,n+1

+ (σ

2
trial
23,n+1)

+ (σ

+ σ

trial
33,n+1)
trial
22,n+1σ
2
trial
12,n+1)
2(cid:1)

trial
31,n+1)

+ 3(cid:0)(σ

+ (σ

ELSEIF a 2D plane stress state is used THEN

σ

trial
22,n+1 = σ22,n + C4∆ε11,n+1 + C2∆ε22,n+1

+ C6∆ε33,n+1

σ

σ

trial
33,n+1 = 0
trial
12,n+1 = σ12,n + C7∆ε12,n+1

(cid:48) trial
n+1 : σ

(cid:48) trial
n+1 = (σ

σ

2
trial
11,n+1)

+ (σ

2
trial
22,n+1)

3

2

− σ

trial
11,n+1σ

trial
22,n+1 + 3(cid:0)(σ

trial
12,n+1)

2(cid:1)

ELSE a 1D stress state is used and

σ

σ

trial
22,n+1 = 0
trial
33,n+1 = 0

(cid:48) trial
n+1 : σ

(cid:48)trial
n+1 = (σ

2
trial
11,n+1)

3

2

σ

END IF

eq,n+1 = ϕ(σtrial
trial
σ

n+1 ) =

(cid:115)

3

2

σ

(cid:48) trial
n+1

: σ

(cid:48) trial
n+1

59

5. Deﬁne initial values of internal variables, to be used in return mapping algorithm, to the converged
= 0
values of the previous step at tn. Check if temperature softening is included by the user (i.e., m
in material input), compute a value for the yield function f (σtrial
n+1 ) and set the incremental plastic
multiplier ∆λn+1 equal to zero.

Goverening equations

Implementation in EUROPLEXUS

*--------------------------------------------------------------

PN = STATE(3)
RP = STATE(6)
SIGMAY = SY + RP

*
* import old value damage variables

WE = STATE(14)

*

*

IF (TEMP_SOFT) THEN

TSTAR = (T-TR) / (TM-TR)
TSOFT = R1 - TSTAR**M
F = PHITRIAL - SIGMAY*TSOFT

ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1

F = PHITRIAL - SIGMAY

ENDIF
DLAMBDA = R0
DDLAMBDA = R0
IF (PHITRIAL == R0) THEN

PHI = R1

ELSE

PHI = PHITRIAL

ENDIF

pn = STATE(3)
Rn = STATE(6)

σy,n = A + Rn

Wn = STATE(14)

IF temperature softening is speciﬁed in material

input (m

= 0) THEN

∗

T

=

Γp (Tn ) =

Tn − Tr

Tm − Tr
(cid:16)

1 − T

∗m(cid:17)

fn+1 = σ

trial
eq,n+1 − σy,nΓp(Tn)

ELSE m = 0 and the CPU time can be reduced

by skipping the softening term

fn+1 = σ

trial
eq,n+1 − σy,n

ENDIF

∆λn+1 = 0

δλn+1 = 0

IF ϕtrial

n+1 = 0 THEN

ϕn+1 = 1

ELSE

ENDIF

ϕn+1 = ϕ

trial
n+1 = σ

trial
eq,n+1

60

(cid:54)
(cid:54)
6. Check for plastic admissibility.

Goverening equations

Implementation in EUROPLEXUS

IF f (σtrial

n+1 ) > 0 THEN apply return mapping by using the

cutting plane algorithm (SOLU = 1) to find the incremental

change in the plastic multiplier ∆λn+1.

Superscript i denotes the local iteration counter.

IF (F > R0) THEN

*
* plastic step: apply return mapping
* the cutting plane method uses values from the previous
* iteration to update the plastic parameter, dlambda.
* thus, all subsequent values (until dlambda) are related to
* the previous iteration, "i".
*--------------------------------------------------------------

P = PN
LOOP_NR : DO NRITER = 1, MXITER

(i) Compute the normal nn+1 to the yield surface based on the initial (i.e., trial) or previous stress
conﬁguration. Note that Eqs.(35) and (34) are used to calculate the gradients, i.e., σ12 = σ21,
σ23 = σ32 and σ13 = σ31.

Goverening equations

Implementation in EUROPLEXUS

IF a 3D stress state is used THEN

n11 = ∂f /∂σ11 = (σ11 − 0.5(σ22 + σ33))/σeq

n22 = ∂f /∂σ22 = (σ22 − 0.5(σ33 + σ11))/σeq

n33 = ∂f /∂σ33 = (σ33 − 0.5(σ11 + σ33))/σeq

n12 = ∂f /∂σ12 = 3σ12 /σeq

n23 = ∂f /∂σ23 = 3σ23 /σeq

n31 = ∂f /∂σ31 = 3σ31 /σeq

ELSEIF a plane stress state is used THEN

n11 = ∂f /∂σ11 = (σ11 − 0.5(σ22))/σeq

n22 = ∂f /∂σ22 = (σ22 − 0.5(σ11))/σeq

n33 = 0
n12 = ∂f /∂σ12 = 3σ12 /σeq

ELSE a 1D stress state is used and

n11 = ∂f /∂σ11 = σ11/σeq

* compute the respective gradients of the von mises yield surface

RNV(2:6) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

RNV(1) = (STRESS(1) - R5*(STRESS(2) + STRESS(3))) / PHI
RNV(2) = (STRESS(2) - R5*(STRESS(3) + STRESS(1))) / PHI
RNV(3) = (STRESS(3) - R5*(STRESS(1) + STRESS(2))) / PHI

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

RNV(1) = (STRESS(1) - R5*STRESS(2)) / PHI
RNV(2) = (STRESS(2) - R5*STRESS(1)) / PHI

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

RNV(1) = STRESS(1) / PHI ! TO BE CHECKED ...

END SELECT

*

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (STRESS(4)) IS PRESENT

RNV(4) = R3*STRESS(4) / PHI

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

RNV(4) = R3*STRESS(4) / PHI
RNV(5) = R3*STRESS(5) / PHI
RNV(6) = R3*STRESS(6) / PHI

END SELECT

(ii) Compute hardening modulus HR,n+1 according to Eq. (42)

Goverening equations

Implementation in EUROPLEXUS

H

i
R,n+1 = Q1C1 exp(−C1p

i
n+1)+Q2C2 exp(−C2p

i
n+1)

* compute hardening modulus

HR = QR1*CR1*EXP(-CR1*P) + QR2*CR2*EXP(-CR2*P)

61

(iii) Compute denominator in Eq. (90) using Eqs. (91), (92), (93) and (98)

Goverening equations

Implementation in EUROPLEXUS

IF a 3D stress state is used THEN

nT Cn = C1

(cid:104)
(n11)

2

+ (n22)

2

+ (n33 )

2(cid:105)

+ 2C4 [n11n22 + n22n33 + n33n11 ]

+ 4C7

(cid:104)
(n12)

2

+ (n23)

2

+ (n31)

2 (cid:105)

ELSEIF a plane stress state is used THEN

nT Cn = C1

(cid:104)

2

(n11)

+ (n22 )

2(cid:105)

+ 2C4[n11 n22] + 4C7

(cid:104)
(n12 )

2(cid:105)

ELSE a 1D stress state is used and

2
nT Cn = C1 (n11)



υ

i
p,n+1 =

1 +

∆λi

n+1
˙p0∆tn+1


1 +

C



C



∆λi

n+1
˙p0∆tn+1



C−1



˙p0∆tn+1

∂υi

p,n+1

∂∆λi

n+1

=

IF temperature softening is activated (m
THEN

= 0)

* compute denominator in the residual derivative of dlambda
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION

CASE (0)
DENOMCP = C1*(RNV(1)*RNV(1) + RNV(2)*RNV(2) + RNV(3)*RNV(3))

>

>

*

>

*

+ R2*C4*(RNV(1)*RNV(2) + RNV(2)*RNV(3) + RNV(3)*RNV(1))

CASE (1)
DENOMCP = C1*(RNV(1)*RNV(1) + RNV(2)*RNV(2))

! PLANE STRESS CONDITION (RNV(3)=0)

+ R2*C4*(RNV(1)*RNV(2))

CASE (2)
DENOMCP = C1*(RNV(1)*RNV(1)) ! TO BE CHECKED ...

! UNIAXIAL STRESS CONDITION (RNV(2)=RNV(3)=0)

END SELECT

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS

CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (RNV(4)) IS PRESENT

DENOMCP = DENOMCP + R4*C7*(RNV(4)*RNV(4))

CASE (3)

! ALL THREE SHEAR STRESSES (RNV(4:6)) ARE PRESENT

DENOMCP = DENOMCP + R4*C7*(RNV(4)*RNV(4) + RNV(5)*RNV(5)

+ RNV(6)*RNV(6))

END SELECT

DENOMV = (R1 + DLAMBDA / (PDOT*DT))**C
DENOMDV = (C / (PDOT*DT))*(R1 + DLAMBDA / (PDOT*DT))**(C - R1)
IF (TEMP_SOFT) THEN

TSTAR = (T - TR) / (TM - TR)
TSOFT = R1 - TSTAR**M
DENOM = DENOMCP + HR*DENOMV*TSOFT + SIGMAY*DENOMDV*TSOFT
ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1

DENOM = DENOMCP + HR*DENOMV + SIGMAY*DENOMDV

ENDIF

(cid:32)

Γp(Tn) =

1 −

DENOM = nT Cn + υ

(cid:35)m(cid:33)

(cid:34) Tn − Tr
Tm − Tr
i
p,n+1Γp,nH

i
R,n+1

+ σ

i
y,n+1Γp,n

∂υi

p,n+1

∂∆λi

n+1

ELSE m = 0 and skip the softening term

DENOM = nT Cn + υ

i
p,n+1H

i
R,n+1

+ σ

i
y,n+1

∂υi

p,n+1

∂∆λi

n+1

ENDIF

(iv) Compute incremental change in the plastic multiplier δλi+1

n+1 in Eq. (90) and update variable

∆λi+1
n+1

Goverening equations

Implementation in EUROPLEXUS

δλ

i+1
n+1

=

f i
n+1

DENOM

∆λ

i+1
n+1

= ∆λ

i
n+1 + δλ

i+1
n+1

* compute residual increment ddlambda and update variable dlambda

DDLAMBDA = F / DENOM
DLAMBDA = DLAMBDA + DDLAMBDA

62

(cid:54)
(v) Update internal variables dependent on ∆λi+1
n+1

Goverening equations

Implementation in EUROPLEXUS

p

i+1
n+1

= p

i
n+1 + δλ

i+1
n+1

* update internal variables dependent on dlambda

P = P + DDLAMBDA
R = RP + HR*DLAMBDA
SIGMAY = SY + R

R

i+1
n+1

σ

i+1
y,n+1

= Rn + H

i
R,n+1∆λ

i+1
n+1

= A + R

i+1
n+1

(vi) Update stress components using the elastic Hooke’s law according to Eq. (83). Note that the
material routine uses a a fully vectorized representation of the stress and strain vector, i.e., σ =
[σ11, σ22, σ33, σ12, σ23, σ31]T and ε = [ε11, ε22, ε33, γ12, γ23, γ31]T = [ε11, ε22, ε33, 2ε12, 2ε23, 2ε31]T

Goverening equations

Implementation in EUROPLEXUS

(C1 n11 + C4n22 + C5n33 )

STRESS(1) = STRESS(1) -

* update stress components using the elastic (hooke’s) law

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION

CASE (0)

DDLAMBDA*(C1*RNV(1) + C4*RNV(2) + C5*RNV(3))

STRESS(2) = STRESS(2) -

DDLAMBDA*(C4*RNV(1) + C2*RNV(2) + C6*RNV(3))

STRESS(3) = STRESS(3) -

DDLAMBDA*(C5*RNV(1) + C6*RNV(2) + C3*RNV(3))

CASE (1)

! PLANE STRESS CONDITION (RNV(3)=0)

STRESS(1) = STRESS(1) -

DDLAMBDA*(C1*RNV(1) + C4*RNV(2))

STRESS(2) = STRESS(2) -

DDLAMBDA*(C4*RNV(1) + C2*RNV(2))

CASE (2)

! UNIAXIAL STRESS CONDITION (RNV(2)=RNV(3)=0)

STRESS(1) = STRESS(1) - DDLAMBDA*C1*RNV(1) ! TO BE CHECKED ...

END SELECT

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS

CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (RNV(4)) IS PRESENT

STRESS(4) = STRESS(4) - DDLAMBDA*C7*RNV(4)

CASE (3)

! ALL THREE SHEAR STRESSES (RNV(4:6)) ARE PRESENT

STRESS(4) = STRESS(4) - DDLAMBDA*C7*RNV(4)
STRESS(5) = STRESS(5) - DDLAMBDA*C8*RNV(5)
STRESS(6) = STRESS(6) - DDLAMBDA*C9*RNV(6)

END SELECT

>

>

>

>

>

*

*

IF a 3D stress state is used THEN

σ

i+1
11,n+1

= σ

i
11,n+1 − δλ

σ

i+1
22,n+1

= σ

i
22,n+1 − δλ

σ

i+1
33,n+1

= σ

i
33,n+1 − δλ

σ

i+1
12,n+1

= σ

i
12,n+1 − δλ

σ

i+1
23,n+1

= σ

i
23,n+1 − δλ

σ

i+1
31,n+1

= σ

i
31,n+1 − δλ

i+1
n+1

i+1
n+1

i+1
n+1

i+1
n+1

i+1
n+1

i+1
n+1

(C4 n11 + C2n22 + C6n33 )

(C5 n11 + C6n22 + C3n33 )

C7n12

C8n23

C9n31

ELSEIF a plane stress state is used THEN

σ

i+1
11,n+1

= σ

i
11,n+1 − δλ

σ

i+1
22,n+1

= σ

i
22,n+1 − δλ

i+1
n+1

i+1
n+1

(C1 n11 + C4n22 + C5n33 )

(C4 n11 + C2n22 + C6n33 )

σ

i+1
33,n+1

= 0

σ

i+1
12,n+1

= σ

i
12,n+1 − δλ

i+1
n+1

C7n12

ELSE a 1D stress state is used and

σ

i+1
11,n+1

= σ

i
11,n+1 − δλ

i+1
n+1

C1n11

63

(vii) Compute the updated von Mises equivalent stress using Eq. (35)

Goverening equations

Implementation in EUROPLEXUS

IF a 3D stress state is used THEN

(cid:48) ,i+1
n+1

σ

: σ

(cid:48),i+1
n+1

3

2

= (σ

i+1
11,n+1

2

)

+ (σ

i+1
22,n+1

2

)

+ (σ

i+1
33,n+1

)

2 − σ

i+1
11,n+1

σ

i+1
22,n+1

+ σ

i+1
22,n+1

σ

i+1
33,n+1

− σ

i+1
33,n+1

σ

i+1
11,n+1

+ 3(cid:0)(σ

i+1
12,n+1

2

)

+ (σ

i+1
23,n+1

2

)

+ (σ

i+1
31,n+1

2(cid:1)

)

ELSEIF a plane stress state is used THEN

(cid:48) ,i+1
n+1

σ

: σ

(cid:48),i+1
n+1

3

2

= (σ

i+1
11,n+1

2

)

+ (σ

i+1
22,n+1

2

)

− σ

i+1
11,n+1

σ

i+1
22,n+1

+ 3(cid:0)(σ

i+1
12,n+1

2(cid:1)

)

ELSE a 1D stress state is used and

(cid:48) ,i+1
n+1

σ

: σ

(cid:48) ,i+1
n+1

= (σ

i+1
11,n+1

2

)

3

2

ENDIF

σ

i+1
eq,n+1

= ϕ(σi+1
n+1

) =

(cid:115)

3

2

σ

(cid:48),i+1
n+1

: σ

(cid:48) ,i+1
n+1

* compute updated equivalent stress

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION

CASE (0)

PHI = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)

+ STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
- STRESS(2)*STRESS(3) - STRESS(3)*STRESS(1)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

PHI = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)

- STRESS(1)*STRESS(2)

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

PHI = STRESS(1)*STRESS(1)

END SELECT

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)
! ONLY TAU_XY (STRESS(4)) IS PRESENT
CASE (1)

PHI = PHI + R3*(STRESS(4)*STRESS(4))

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

PHI = PHI + R3*(STRESS(4)*STRESS(4)

+ STRESS(5)*STRESS(5)
+ STRESS(6)*STRESS(6))

END SELECT

PHI = SQRT (PHI)

>
>

>

*

>
>

*

(viii) Update the yield function in Eq. (84) according to the new stress state computed in (vi) and

(vii). Check for convergence.

Goverening equations

Implementation in EUROPLEXUS

* compute yield function and new residual

IF (TEMP_SOFT) THEN

TSTAR = (T - TR) / (TM - TR) ! THIS RECALCULATION OF TSTAR
TSOFT = R1 - TSTAR**M
F = PHI - SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT

! AND TSOFT SEEMS REDUNDANT HERE

*
* check convergence

RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT)
RESNOR = ABS (RESNOR)

ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1

F = PHI - SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)

* check convergence

RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C))
RESNOR = ABS (RESNOR)

ENDIF

*

IF (RESNOR <= TOL) THEN

IF temperature softening (m

= 0) THEN

f

i+1
v,n+1

= ϕ

i+1
n+1

− σ

i+1
y,n+1

(cid:16)

(cid:17)

•



(cid:16)

(cid:17)

•

=


1 +

∆λ

i+1
n+1
˙p0∆tn+1



C

(cid:32)




1 −

(cid:34) Tn − Tr
Tm − Tr

(cid:35)m(cid:33)

RESNOR =

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)

fv (∆λ

σ

i+1
y,n+1

υ

)

i+1
n+1
i+1
p,n+1

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)

Γp

ELSE m = 0 and skip the softening term

f

i+1
v,n+1

= ϕ

i+1
n+1

− σ

i+1
y,n+1

(cid:16)

(cid:17)

•



(cid:16)

(cid:17)

•

=


1 +

∆λ

i+1
n+1
˙p0∆tn+1



C




RESNOR =

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)

σ

ENDIF

fv (∆λ
i+1
y,n+1

)

i+1
n+1
i+1
p,n+1

υ

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)

IF

RESNOR < TOL THEN continue to 7.

ELSE i = i + 1 and GO TO (i)

64

(cid:54)
7. Update the temperature according to the low coupling approach using Eq. (45) and erode element if

T > Tm

Goverening equations

Implementation in EUROPLEXUS

Tn+1 = Tn +

χσeq,n+1∆λn+1
ρCT

* update temperature (only if temperature softening is activated)

IF (TEMP_SOFT) THEN

T = T + TQ*PHI*DLAMBDA / (RHO*CP)
IF (T >= TM) THEN ! CHECK IF MELTING OCCURS

STATE(13) = TM
D = DC ! this will cause the GP to fail (see 990)
GO TO 990

ENDIF

ENDIF

8. Compute the hydrostatic stress according to Eq. (32)

Goverening equations

Implementation in EUROPLEXUS

σH,n+1 =

1

3

tr(σn+1) =

1

3

(σ11,n+1+σ22,n+1+σ33,n+1)

* compute hydrostatic stress based on new stress state,

SIGMAH = STRESS(1)
SELECT CASE (IPLANC)

CASE (0)

! NO PLANE STRESS CONDITION

SIGMAH = SIGMAH + STRESS(2) + STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

SIGMAH = SIGMAH + STRESS(2)

CASE (2)
END SELECT

*

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

SIGMAH = SIGMAH / R3

9. Computing the principal stresses σi,n+1 based on the principal deviatoric stresses σ

(cid:48)

i,n+1 as shown in

Eqs. (55) and (57)

Goverening equations

Implementation in EUROPLEXUS

(cid:48)
11,n+1 = σ11,n+1 − σH,n+1
(cid:48)
22,n+1 = σ22,n+1 − σH,n+1
(cid:48)
33,n+1 = σ33,n+1 − σH,n+1

σ

σ

σ

J2,n+1 =

1

2

(cid:18)

(cid:48)
2
11,n+1)

(σ

+ (σ

(cid:48)
2
22,n+1)

+ (σ

(cid:48)
2
33,n+1)

(cid:19)

+ (σ12,n+1)

2

J3,n+1 = σ

(cid:48)
2
11,n+1(σ23,n+1)

+σ

(cid:48)
33,n+1(σ12,n+1)

2

2

+ (σ31,n+1)

+ (σ23,n+1)
(cid:48)
2
22,n+1(σ31,n+1)
(cid:48)
33,n+1

(cid:48)
22,n+1 + σ

(cid:48)
11,n+1σ

2 − σ

+ σ

− 2σ12,n+1σ23,n+1σ31,n+1

θ = arccos







(cid:114)

2

J3,n+1
J3

2,n+1

/27







,

0 ≤ θ ≤ π

(cid:115)

(cid:48)
1,n+1 = 2

(cid:115)

(cid:48)
2,n+1 = 2

(cid:115)

(cid:48)
3,n+1 = 2

σ

σ

σ

J2,n+1
3

J2,n+1
3

J2,n+1
3

θ

3
(cid:32) θ

3
(cid:32) θ

cos

cos

cos

(cid:33)

2π

3

(cid:33)

2π

−

+

3

3

σ1,n+1 = σH,n+1 + σ

σ2,n+1 = σH,n+1 + σ

σ3,n+1 = σH,n+1 + σ

(cid:48)
1,n+1
(cid:48)
2,n+1
(cid:48)
3,n+1

* 3d stress state and compute deviatoric stress

S1 = STRESS(1) - SIGMAH
S2 = STRESS(2) - SIGMAH
S3 = STRESS(3) - SIGMAH

*
* finding principal stress based on the deviatoric stress and invariants

A1 = R5*(S1*S1 + S2*S2 + S3*S3)

*

>
>

>
>
>

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)

A2 = -S1*S2*S3 ! DETERMINANT OF THE DEVIATORIC STRESS TENSOR

CASE (1)

! ONLY TAU_XY (STRESS(4)) IS PRESENT

A1 = A1 + STRESS(4)*STRESS(4)
A2 = S3*STRESS(4)*STRESS(4) - S1*S2*S3

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

A1 = A1 +

STRESS(4)*STRESS(4) + STRESS(5)*STRESS(5) +
STRESS(6)*STRESS(6)

A2 = S1*STRESS(5)*STRESS(5) +

S2*STRESS(6)*STRESS(6) +
S3*STRESS(4)*STRESS(4) -
S1*S2*S3 - R2*STRESS(4)*STRESS(5)*STRESS(6)

END SELECT

A3 = -SQRT (R27 / A1)*A2*R5 / A1
A3 = SIGN (MIN (ABS (A3), R1), A3)
A4 = ACOS (A3) / R3

*
* compute (chronological) principal stresses

AUX = SQRT (A1 / R3)
S1 = SIGMAH + R2*AUX*COS(A4)
S2 = SIGMAH + R2*AUX*COS(A4-R2*PI/R3) ! only for postprocessing
S3 = SIGMAH + R2*AUX*COS(A4+R2*PI/R3) ! only for postprocessing

*

65

10. Compute the damage parameter Dn+1 by using Eq. (46), i.e.

integrating the major principal stress

during the entire equivalent plastic strain path

Goverening equations

Implementation in EUROPLEXUS

Wn+1 = Wn + max(σ1,n+1, 0)∆λn+1

Dn+1 =

Wn+1
Wc

*

*

WE = WE + MAX (S1, R0)*DLAMBDA
D = WE / WC

11. Update state variables to be returned to the element routine in EUROPLEXUS and check for element

erosion.

Goverening equations

Implementation in EUROPLEXUS

ST AT E(1) = σH,n+1

ST AT E(2) = σeq,n+1

ST AT E(3) = pn+1

ST AT E(4) = σ

trial
eq,n+1

ST AT E(5) = fn+1

ST AT E(6) = Rn+1

ST AT E(7) = δλn+1

ST AT E(8) = ∆λn+1

ST AT E(9) = N RIT ER

ST AT E(10) = ˙pn+1 = ∆λn+1/∆tn+1

ST AT E(11) = Dn+1

ST AT E(13) = Tn+1

ST AT E(14) = Wn+1

ST AT E(15) = CSONn+1

ST AT E(16) = σ1,n+1

ST AT E(17) = σ2,n+1

ST AT E(18) = σ3,n+1

IF

Dn+1 > DC

THEN

failure and the material point stiffness is deleted.

STATE(12) = 1

ELSE

no failure and the material point stiffness should

not be deleted.

STATE(12) = 0

* update state variables before exit

STATE(1) = SIGMAH
STATE(2) = PHI
STATE(3) = P
STATE(4) = PHITRIAL
STATE(5) = F
STATE(6) = R
STATE(7) = DDLAMBDA
STATE(8) = DLAMBDA
STATE(9) = NRITER
STATE(10) = DLAMBDA / DT
STATE(11) = D

STATE(13) = T
STATE(14) = WE
STATE(15) = CSON

STATE(16) = S1
STATE(17) = S2
STATE(18) = S3

*

*

*
* decide if material point stiffness is zero, i.e. failed
*--------------------------------------------------------------

990 IF (D >= DC) THEN

*
* the current gp was already failed or is failing now :
* state(1:10,12,13), the relevant components
* of stressnew and cson are set to 0.
* the total strains (strain) are not set to 0 so they continue to be
* updated as long as the element is not eroded, i.e. as long
* as the element has some virgin gauss points.
*

SELECT CASE (ITAU)
CASE (0)

STRESSNEW(1:3) = R0

CASE (1)

STRESSNEW(1:4) = R0

CASE (2)

STRESSNEW(1:5) = R0

CASE (3)

STRESSNEW(1:6) = R0

END SELECT
STATE(1:10) = R0

* state(11) (damage d) keeps last value it had before failure

STATE(12) = R0

* state(13) (absolute temperature t) keeps last value before failure
* state(14) (damage accumulation we) keeps last value before failure

CSON = R0
STATE(15) = CSON

ELSE

*
* no failure
*

STATE(12) = R1

*

66

12. Update the actual total strain ε33,n+1 and incremental strain ∆ε33,n+1 through the thickness according

to Eq. (69) and return this to the element routine in EUROPLEXUS

Goverening equations

Implementation in EUROPLEXUS

IF plane stress THEN use Eq.(69) to compute ac-
tual strain trough thickness

∆ε33,n+1 =

1 − 2ν

E

(∆σ11,n+1 + ∆σ22,n+1)

− (∆ε11,n+1 + ∆ε22,n+1)

ε33,n+1 = ε33,n + ∆ε33,n+1

* strain increment(s) and total strain(s) corresponding to
* the imposed zero stress condition(s)
*

SELECT CASE (IPLANC)
CASE (0)
CASE (1)

! NO PLANE STRESS CONDITION
! PLANE STRESS CONDITION (STRESS(3)=0)

* store provisional strain increment

AUX3 = STRAININC(3)

* compute and update the actual strain increment

FACTOR = T2 / YOUNG
STRAININC(3) = FACTOR*(STRESSNEW(1) - STRESSOLD(1)

>
>
* update the total strain

+ STRESSNEW(2) - STRESSOLD(2))
- STRAININC(1) - STRAININC(2)

STRAIN(3) = STRAIN(3) - AUX3 + STRAININC(3)

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

AUX2 = STRAININC(2)
AUX3 = STRAININC(3)
FACTOR = T2 / YOUNG
STRAININC(2) = 0.5D0*(FACTOR*(STRESSNEW(1) - STRESSOLD(1))

>

- STRAININC(1))

STRAININC(3) = STRAININC(2)
STRAIN(2) = STRAIN(2) - AUX2 + STRAININC(2)
STRAIN(3) = STRAIN(3) - AUX3 + STRAININC(3)

END SELECT

ENDIF
*
* erosion (optional)
*

IF (L_EROSION) THEN

IF (D >= DC) FAIL_FLAG = .TRUE.

ENDIF

*
END SUBROUTINE VPJC

In EPX, this is done for each material integration point in the FE assembly. The stresses are updated
in the material interface using a fully vectorized version of the forward Euler integration algorithm and
a two-state architecture where the initial values at tn are stored in the old arrays and the new values at
tn+1 must be updated and stored in the new arrays returned to the global FE analysis. The VPJC ma-
terial routine is valid for 1D, 2D and 3D stress states, i.e., for bar, shell, solid, axisymmetric, plane strain
and plane stress elements. The stresses for 3D elements are stored similar to that of symmetric tensors
as σ = [σ11, σ22, σ33, σ12, σ23, σ31]T , and plane stress, axisymmetric and plane strain elements are stored as
σ = [σ11, σ22, σ33, σ12]T . The deformation gradient and strains are stored similarly to the stresses, however,
one should be aware of that the shear strain is stored as engineering shear strains, e.g. γ12 = 2ε12.

Appendix F. Material routine in EUROPLEXUS

The constitutive model of elastic-thermoviscoplasticity and ductile failure is implemented in EUROPLEXUS
(EPX) [7] using a fully vectorized version of the forward Euler integration algorithm and a special case of the
backward Euler integration algorithm as presented in this report. The FORTRAN implementation of the
numerical scheme in Appendix E is presented in this section. The model is applicable for 1D, 2D and 3D
stress states, i.e., for 1D, 2D and 3D stress states, i.e., for bar, shell, solid, axisymmetric, plane strain and
plane stress elements. The model is coupled with an element deletion option available in EPX by introducing
a state variable controlling the element erosion, i.e., as the damage variable reaches its critical value DC the
stiﬀness in the material point is removed. An element will be removed from the mesh when a user-deﬁned
fraction of the material points in the element are deleted. This is speciﬁed in the input ﬁle for EPX using
the EROS directive where the user should specify a (proportional) number between 0 and 1. The default
value is 1 and indicates that an element is eroded when all its Gauss points have failed, while the value 0
indicate that an element is eroded when any one of its Gauss points fails.

The procedure described in Section 5.1 and Appendix E is formulated in terms of the material subroutine
VPJC, which communicates with the element routine (or even a higher-order routine) in EPX. During an

67

FE analysis the material subroutine VPJC is called by the chosen element routine for each material point
at each increment. When the material subroutine is called, it is provided with the state at the start of
the increment (stress and solution-dependent state variables at tn). Based on the geometrical calculations
performed by the element routine it is also provided with the strain increment ∆ε and total strains ε at the
beginning and the end of the increment. Thus, the VPJC uses a two-state architecture for the stress tensor
σ where the initial values at tn are stored in the old arrays and the new values at tn+1 must be updated and
stored in the new arrays returned to the global analysis in EPX. The stresses for 3D elements are stored
similar to that of symmetric tensors as σ = [σ11, σ22, σ33, σ12, σ23, σ31]T , while plane stress, plane strain
and axisymmetric elements are stored as σ = [σ11, σ22, σ33, σ12]T . The deformation gradient and strains are
stored similarly to the stresses, however, one should be aware of that the shear strain is stored as tensor
shear strains, e.g. ε12 = 1
2 γ12. That is, ε = [ε11, ε22, ε33, 2ε12, 2ε23, 2ε31]T and ε = [ε11, ε22, ε33, 2ε12]T in
the 3D and 2D cases, respectively.

The coding of material routines in EUROPLEXUS follow FORTRAN 90 conventions. In addition to the
vectorized interface, it is essential that all variables are deﬁned and initialized properly using a routine that
reads the material constants from the input data set. Moreover, a conversion routine is used and called by
the element routine just before and after the material routine, to make the arrangement of stress and strain
in VPJC compatible with the one used in the element routine. Thus, the conversation routine ensures that
all tensors have the same number of components and the same organization for a given element type.

The FORTRAN source of the material routine for EPX is listed below.

C M_MATERIAL_VPJC

TOUS GVALSAMOS 16/06/06 20:09:13 #3077

MODULE M_MATERIAL_VPJC

*
* material of type vpjc
*
* the vpjc subroutine can be used for all stress states, by providing
* appropriate values of itau and iplanc
*

USE M_MATERIALS_ARRAY

*

*

*

IMPLICIT NONE
SAVE

PRIVATE
PUBLIC ::

> MAVPJC, ! PROPERTIES READING ROUTINE
> VPJC

! MATERIAL ROUTINE

CONTAINS

*===============================================================================

SUBROUTINE MAVPJC (INUMLDC)

*------------------------------------------------------------------
*
vegard aune, folco casadei 11-14
* material "vpjc" (visco-plastic johnson-cook)
*------------------------------------------------------------------
* mat : index of material read in the input file
*
* organization of newmat%matval: (material properties)
* ----------------------------------------------------
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*

- material density
- elastic modulus
- possion ratio
- initial yield stress
- tolerance newton iteration
- maximum number of iterations
- asymptote first term voce hardening
- hardening parameter first term voce
- asymptote 2nd term voce hardening
- hardening parameter 2nd term voce
- reference strain-rate
- hardening parameter visco term
- taylor-quinney coefficient
- specific heat capacity solid material
- room temperature

matval(1) - rho
matval(2) - young
matval(3) - nu
matval(4) - elas
matval(5) - tol
matval(6) - mxit
matval(7) - qr1
matval(8) - cr1
matval(9) - qr2
matval(10)- cr2
matval(11)- pdot
matval(12)- c
matval(13)- tq
matval(14)- cp
matval(15)- tr

68

matval(16)- tm
matval(17)- m
matval(18)- dc
matval(19)- wc
matval(20)- solu

- melting temperature
- hardening parameter temperature term
- upper limit of d when softening
- cockcroft-latham failure criterion
- 1=cutting plane (default) explicit

requires N-R iterations,

2=radial return (implicit) requires N-R iterations

INCLUDE ’NONE.INC’

INCLUDE ’CARMA.INC’
INCLUDE ’CONSTA.INC’
INCLUDE ’CUNIT.INC’
INCLUDE ’REDCOM.INC’

*
*
*
*
*
*
*
*

*

*

INTEGER, INTENT(IN) :: INUMLDC

*
* local variables
*

INTEGER, PARAMETER :: IFIX = 4

INTEGER, PARAMETER :: NMOT = 20

CHARACTER*4 :: MOT(NMOT)
INTEGER :: KOPT(NMOT), I

*

*

*

DATA MOT /’RO ’,’YOUN’,’NU ’,’ELAS’,’TOL ’,
’MXIT’,’QR1 ’,’CR1 ’,’QR2 ’,’CR2 ’,
’,’TQ ’,’CP ’,’TR ’,
’PDOT’,’C
’,’DC ’,’WC ’,’SOLU’/
’TM ’,’M

>
>
>

*
* the syntax of this material is particularly simple, because
* it contains only a sequence of keywords, each followed by a value.
* the first 4 keywords are mandatory, the others are optional
*
* create new material structure newmat of the needed length
*

CALL CREATE_MATERIAL (NMOT, 0, 0)
NEWMAT%NAME = ’VPJC : VISCO-PLASTIC JOHNSON-COOK’
NEWMAT%TYPE = 120 ! TYPE OF THE VPJC MATERIAL
NEWMAT%NUMLDC = INUMLDC
NEWMAT%LGECR = LGECR(NEWMAT%TYPE)

*

CALL LIRVAL (NMOT, MOT, NEWMAT%MATVAL, KOPT)

*
* input data complete ? (all keywords are mandatory except tol, mxit, tr, solu)
*

DO I = 1, NMOT

SELECT CASE (I)
CASE (5, 6, 15, 20) ! TOL, MXIT, TR, SOLU ARE OPTIONAL
CASE DEFAULT

! ALL OTHERS ARE MANDATORY

IF (KOPT(I) == 0) GO TO 9

END SELECT

END DO

*
* set default values in the optional properties which have not been read
*

IF (KOPT(5) == 0) NEWMAT%MATVAL(5) = 1.E-4
IF (KOPT(6) == 0) NEWMAT%MATVAL(6) = 20
IF (KOPT(15) == 0) NEWMAT%MATVAL(15) = 293.D0 ! TR
IF (KOPT(20) == 0) NEWMAT%MATVAL(20) = 1.D0

! TOL
! MXIT

! SOLU

*
* printouts
*

WRITE(TAPOUT,1000) NEWMAT%NUMLDC

*

WRITE(TAPOUT,1001) NEWMAT%MATVAL(1)
WRITE(TAPOUT,1002) NEWMAT%MATVAL(2)
WRITE(TAPOUT,1003) NEWMAT%MATVAL(3)
WRITE(TAPOUT,1004) NEWMAT%MATVAL(4)
WRITE(TAPOUT,1005) NEWMAT%MATVAL(5)
WRITE(TAPOUT,1006) NINT (NEWMAT%MATVAL(6))
WRITE(TAPOUT,1007) NEWMAT%MATVAL(7)
WRITE(TAPOUT,1008) NEWMAT%MATVAL(8)
WRITE(TAPOUT,1009) NEWMAT%MATVAL(9)
WRITE(TAPOUT,1010) NEWMAT%MATVAL(10)
WRITE(TAPOUT,1011) NEWMAT%MATVAL(11)

69

WRITE(TAPOUT,1012) NEWMAT%MATVAL(12)
WRITE(TAPOUT,1013) NEWMAT%MATVAL(13)
WRITE(TAPOUT,1014) NEWMAT%MATVAL(14)
WRITE(TAPOUT,1015) NEWMAT%MATVAL(15)
WRITE(TAPOUT,1016) NEWMAT%MATVAL(16)
WRITE(TAPOUT,1017) NEWMAT%MATVAL(17)
WRITE(TAPOUT,1018) NEWMAT%MATVAL(18)
WRITE(TAPOUT,1019) NEWMAT%MATVAL(19)
WRITE(TAPOUT,1020) NINT (NEWMAT%MATVAL(20))

*

RETURN

*
* errors
*

9 CALL ERRMSS (’MAVPJC’,’DIRECTIVE "VPJC" INCOMPLETE’)

STOP ’MAVPJC’

*

1000 FORMAT(///,10X,’MATERIAL NUMBER ’,I5,’ ( VPJC )’,/)
1001 FORMAT(30X,’DENSITY’,T65,’= ’,1PE12.5)
1002 FORMAT(30X,’ELASTIC MODULUS’,T65,’= ’,1PE12.5)
1003 FORMAT(30X,’POISSON RATIO’,T65,’= ’,1PE12.5)
1004 FORMAT(30X,’INITIAL YIELD STRESS’,T65,’= ’,1PE12.5)
1005 FORMAT(30X,’TOLERANCE FOR N-R ITERATIONS’,T65,’= ’,1PE12.5)
1006 FORMAT(30X,’MAX NUMBER OF N-R ITERATIONS’,T65,’= ’,I6)
1007 FORMAT(30X,’ASYMPTOTE FIRST TERM VOCE HARDENING’,T65,’= ’,1PE12.5)
1008 FORMAT(30X,’HARDENING PARAMETER FIRST TERM VOCE’,T65,’= ’,1PE12.5)
1009 FORMAT(30X,’ASYMPTOTE 2ND TERM VOCE HARDENING’,T65,’= ’,1PE12.5)
1010 FORMAT(30X,’HARDENING PARAMETER 2ND TERM VOCE’,T65,’= ’,1PE12.5)
1011 FORMAT(30X,’REFERENCE STRAIN RATE’,T65,’= ’,1PE12.5)
1012 FORMAT(30X,’HARDENING PARAMETER VISCOUS TERM’,T65,’= ’,1PE12.5)
1013 FORMAT(30X,’TAYLOR-QUINNEY COEFFICIENT’,T65,’= ’,1PE12.5)
1014 FORMAT(30X,’SPECIFIC HEAT CAPACITY OF THE SOLID’,T65,’= ’,1PE12.5)
1015 FORMAT(30X,’ROOM TEMPERATURE’,T65,’= ’,1PE12.5)
1016 FORMAT(30X,’MELTING TEMPERATURE’,T65,’= ’,1PE12.5)
1017 FORMAT(30X,’HARDENING PARAM. TEMPERATURE TERM’,T65,’= ’,1PE12.5)
1018 FORMAT(30X,’UPPER LIMIT OF D WHEN SOFTENING’,T65,’= ’,1PE12.5)
1019 FORMAT(30X,’COCKCROFT-LATHAM FAILURE CRITERION’,T65,’= ’,1PE12.5)
1020 FORMAT(30X,’SOLUTION (1=CUT PLANE, 2=RAD.RET.)’,T65,’= ’,I6)

*

END SUBROUTINE MAVPJC

*-------------------------------------------------------------------------------

SUBROUTINE VPJC (ITAU, IPLANC, PROPS, STRAININC, STRESSOLD,
STRAIN, STATE, STRESSNEW, CSON, FAIL_FLAG)

>

*
*
*
* added optional radial return solution strategy (solu == 2)
*
*erosion

vegard aune, folco casadei 12-14

02-16

USE M_FAILED_ELEMENTS ! for l_erosion

USE M_PILOT
USE M_KAAPI_DATA ! FOR ELEMENT PRINTOUTS

*

*

IMPLICIT NONE

*
* args

INTEGER, INTENT(IN) :: ITAU, IPLANC
LOGICAL, INTENT(INOUT) :: FAIL_FLAG
REAL(8), INTENT(IN) :: PROPS(*), STRESSOLD(*)
REAL(8), INTENT(OUT) :: STRESSNEW(*), CSON
REAL(8), INTENT(INOUT) :: STRAININC(*), STRAIN(*), STATE(*)

*
* locals

REAL(8) :: STRESS(6), RNV(6)

*

REAL(8) :: DT, RHO, YOUNG, POISS, R3G
REAL(8) :: C1, C2, C3, C4, C5, C6, C7, C8, C9
REAL(8) :: SY, HR, RP, R, QR1, CR1, QR2, CR2
REAL(8) :: TOL, PDOT, C, TQ, CP, TR, TM, M
REAL(8) :: PN, SIGMAY, DENOMV, DENOMDV, SIGMAH
REAL(8) :: P, DENOM, DDLAMBDA, DLAMBDA, RESNOR, FACTOR
REAL(8) :: PHITRIAL, PHI, F, D, T, DC, WC, AUX, WE, AUX3,

>
>

S1, S2, S3, A1, A2, A3, A4, DENOMCP, AUX2, T1, T2, T3,
PHIEFF

REAL(8), PARAMETER :: R0 = 0.D0, R1 = 1.D0, R2 = 2.D0,

>

R3 = 3.D0, R4 = 4.D0, R5 = 0.5D0,

70

>

R27=27.D0

* defining the number pi

REAL(8), PARAMETER :: PI = 3.1415926535D0

*

INTEGER :: MXITER, NRITER, SOLU
LOGICAL :: TEMP_SOFT ! IS TEMPERATURE SOFTENING ACTIVATED (IS M /= 0?)
REAL(8) :: TSTAR, TSOFT

: shear components

: normal components

sigma_z=stress(3)=0)

(sigma_z=stress(3)=0)

: number of shear stress components (stress(4),stress(5),stress(6))

1 = plane stress condition
2 = uniaxial stress condition (sigma_y=stress(2)=0 and

- material density
- elastic modulus
- possion ratio
- initial yield stress
- tolerance newton iteration

0 = no shear stress is present
1 = only tau_xy (stress(4)) is present
3 = all three shear stresses (stress(4:6)) are present

organization of stress and strain arrays
----------------------------------------
stress(1)- xx
stress(2)- yy
stress(3)- zz
stress(4)- xy
stress(5)- yz
stress(6)- xz

list of material properties sent to the subroutine
--------------------------------------------------
props(1) - rho
props(2) - young
props(3) - poiss
props(4) - sy
props(5) - tol
props(6) - mxiter - maximum number of iterations
props(7) - qr1
props(8) - cr1
props(9) - qr2
props(10)- cr2
props(11)- pdot
props(12)- c
props(13)- tq
props(14)- cp
props(15)- tr
props(16)- tm
props(17)- m
props(18)- dc
props(19)- wc
props(20)- solu

*
* input parameters to generalize the present material routine:
* itau
*
*
*
* iplanc : 0 = no plane stress condition
*
*
*
*
*
*
*
*
*
*
*
*
*
* attention: the shear components of strain and straininc arrays are the
* "engineering" values (gamma), i.e. twice the "tensor" values (epsilon)!
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
* list of state variables (old and new)
* ----------------------------------------------------------------
*
intent
* ----------------------------------------------------------------
new
* state(1) - hydrostatic pressure (1/3 sigma_kk)
new
* state(2) - von mises equivalent stress
p old+new
* state(3) - equivalent plastic strain
new
* state(4) - elastic trial equivalent stress
phitrial
f
* state(5) - yield function (should be near zero)
new
r old+new
* state(6) - total hardening of the material
new
* state(7) - change of dlambda over the step
new
* state(8) - incremental plastic multiplier
new
* state(9) - n of iterations for convergence
new
* state(10)- plastic strain rate
d old+new
* state(11)- damage
new
* state(12)- failure: 1 = virgin g.p., 0 = failed g.p. fa
t old+new
* state(13)- absolute temperature
we old+new
* state(14)- cockcroft-latham damage accumulation
new
* state(15)- sound speed
new
* state(16)- first principal stress
new
* state(17)- second principal stress
new
* state(18)- third principal stress

- asymptote first term voce hardening
- hardening parameter first term voce
- asymptote 2nd term voce hardening
- hardening parameter 2nd term voce
- reference strain-rate
- hardening parameter visco term
- taylor-quinney coefficient
- specific heat capacity solid material
- room temperature
- melting temperature
- hardening parameter temperature term
- upper limit of d when softening
- cockcroft-latham failure criterion
- 1=cutting plane (default) requires N-R iterations,

ddlambda
dlambda
nriter
dlambda / dt

(not available in plane or uniaxial stress!)

2=radial return (direct, no iterations)

cson
s1
s2
s3

sigmah
phi

description

name

71

in input it is always .false.
in output it is .true. if:

- the current Gauss point is failed, **and**
- erosion has been activated by the user

*
* fail_flag : output only!
*
*
*
*
*
*
* NOTE:
* - the material routine is called also for an already failed gp, if the
*
* - however, the material routine is not called at all for an eroded element!
*
! Thread variable

element to which the gp belongs has not been eroded yet!

(i.e. l_erosion == .true.)

INTEGER :: NTH

*

DC = PROPS(18)
D = STATE(11)
T = STATE(13) ! T IS INITIALIZED TO TR IN ECROU.FF

*
* check previous failure of current gp

IF (D >= DC) THEN

* current gauss point was previously failed, so let’s just skip it ...

GO TO 990

ENDIF

*
* current gp was not previously failed, but it might fail now ...
*

CALL PILOT_GET_THREAD_NUM (NTH)

! End thread variable
*

CALL ELEM_DT1 (IEL(NTH), DT)

*

*

RHO
YOUNG
POISS
SY
TOL
MXITER
QR1
CR1
QR2
CR2
PDOT
C
TQ
CP
TR
TM
M
WC
SOLU

= PROPS(1)
= PROPS(2)
= PROPS(3)
= PROPS(4)
= PROPS(5)
= NINT (PROPS(6))
= PROPS(7)
= PROPS(8)
= PROPS(9)
= PROPS(10)
= PROPS(11)
= PROPS(12)
= PROPS(13)
= PROPS(14)
= PROPS(15)
= PROPS(16)
= PROPS(17)
= PROPS(19)
= NINT (PROPS(20))

IF (SOLU == 2) THEN

* radial return is not available in either plane stress or uniaxial stress

IF (IPLANC > 0) THEN

CALL ERRMSS (’VPJC’, ’SOLU 2 N/A IN PLANE/UNIAXIAL STRESS’)
STOP ’VPJC - SOLU 2 N/A IN PLANE/UNIAXIAL STRESS’

ENDIF

ENDIF

*

TEMP_SOFT = (M /= 0.D0) ! TEMPERATURE SOFTENING IS ACTIVATED ONLY

! IF M /= 0 (EXPONENTIALS ARE EXPENSIVE ...)

*
* sound speed

T1 = R1 + POISS
T2 = R1 - R2*POISS
T3 = R1 - POISS

*

SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

CSON = SQRT (YOUNG*T3 / (RHO*T1*T2))

CASE (1)

! PLANE STRESS CONDITION (SIG(3)=0)

CSON = SQRT (YOUNG / (RHO*T1*T3))

CASE (2)

! UNIAXIAL STRESS CONDITION (SIG(2)=SIG(3)=0)

CSON = SQRT (YOUNG / RHO)

CASE DEFAULT ! WE CHECK THIS ONLY IN THE FIRST SELECT CASE ON IPLANC

72

CALL ERRMSS (’VPJC’, ’WRONG IPLANC’)
STOP ’VPJC - WRONG IPLANC’

END SELECT

*
*--------------------------------------------------------------
* state update procedure for the von mises elasto-
* thermoviscoplastic material model with non-linear isotropic
* hardening goverened by the modified johnson-cook model:
* explicit elastic predictor/return mapping algorithm.
*--------------------------------------------------------------
*
* elastic predictor: compute elastic trial stress state
*--------------------------------------------------------------
*
* compute the elasticity tensor
SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

C1 = YOUNG*T3 / (T1*T2)
C4 = POISS*C1 / T3

CASE (1)

! PLANE STRESS CONDITION (SIG(3)=0)

C1 = YOUNG / (T1*T3)
C4 = POISS*C1

CASE (2)

! UNIAXIAL STRESS CONDITION (SIG(2)=SIG(3)=0)

C1 = YOUNG
C4 = R0
END SELECT

! THESE VALUES HAVE TO BE CHECKED

*

C2 = C1
C3 = C1
C5 = C4
C6 = C4
C7 = R5*YOUNG / T1
C8 = C7
C9 = C7
R3G = R3*C7

*
* compute the elastic trial stress and the elastic trial equivalent stress
*

STRESS(1) = STRESSOLD(1) + C1*STRAININC(1) + C4*STRAININC(2)

>

*

+ C5*STRAININC(3)

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

STRESS(2) = STRESSOLD(2) + C4*STRAININC(1) + C2*STRAININC(2)

>

>

>
>

+ C6*STRAININC(3)

STRESS(3) = STRESSOLD(3) + C5*STRAININC(1) + C6*STRAININC(2)

+ C3*STRAININC(3)

PHITRIAL = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)
+ STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
- STRESS(2)*STRESS(3) - STRESS(3)*STRESS(1)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

STRESS(2) = STRESSOLD(2) + C4*STRAININC(1) + C2*STRAININC(2)

>

>

+ C6*STRAININC(3)

STRESS(3) = R0
PHITRIAL = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)

- STRESS(1)*STRESS(2)

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

STRESS(2) = R0
STRESS(3) = R0
PHITRIAL = STRESS(1)*STRESS(1)

END SELECT

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)
! ONLY TAU_XY (STRESS(4)) IS PRESENT
CASE (1)

STRESS(4) = STRESSOLD(4) + C7*STRAININC(4) ! THE GAMMAS ARE USED!
PHITRIAL = PHITRIAL + R3*(STRESS(4)*STRESS(4))

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT
STRESS(4) = STRESSOLD(4) + C7*STRAININC(4) ! SHEAR STRAIN AND STRAININC
STRESS(5) = STRESSOLD(5) + C8*STRAININC(5) ! ARE THE GAMMAS, NOT
STRESS(6) = STRESSOLD(6) + C9*STRAININC(6) ! THE EPSILONS
PHITRIAL = PHITRIAL + R3*(STRESS(4)*STRESS(4)
+ STRESS(5)*STRESS(5)
+ STRESS(6)*STRESS(6))

>
>

CASE DEFAULT ! WE CHECK THIS ONLY IN THE FIRST SELECT CASE ON ITAU

CALL ERRMSS (’VPJC’, ’WRONG ITAU’)
STOP ’VPJC - WRONG ITAU’

73

END SELECT

*

PHITRIAL = SQRT (PHITRIAL)

*
* check for plastic admissibility
*--------------------------------------------------------------

PN = STATE(3)
RP = STATE(6)
SIGMAY = SY + RP

*
* import old value damage variables

WE = STATE(14)

*

*

IF (TEMP_SOFT) THEN

TSTAR = (T-TR) / (TM-TR)
TSOFT = R1 - TSTAR**M
F = PHITRIAL - SIGMAY*TSOFT

ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1

F = PHITRIAL - SIGMAY

ENDIF
DLAMBDA = R0
DDLAMBDA = R0
IF (PHITRIAL == R0) THEN

PHI = R1

ELSE

PHI = PHITRIAL

ENDIF

IF (SOLU == 2) GO TO 200

*
*=======================================================================
* solu == 1 : cutting plane solution strategy (default)
*=======================================================================
*

IF (F > R0) THEN

*
* plastic step: apply return mapping
* the cutting plane method uses values from the previous
* iteration to update the plastic parameter, dlambda.
* thus, all subsequent values (until dlambda) are related to
* the previous iteration, "i".
*--------------------------------------------------------------

P = PN
LOOP_NR : DO NRITER = 1, MXITER

*
* compute the respective gradients of the von mises yield surface

RNV(2:6) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

RNV(1) = (STRESS(1) - R5*(STRESS(2) + STRESS(3))) / PHI
RNV(2) = (STRESS(2) - R5*(STRESS(3) + STRESS(1))) / PHI
RNV(3) = (STRESS(3) - R5*(STRESS(1) + STRESS(2))) / PHI

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

RNV(1) = (STRESS(1) - R5*STRESS(2)) / PHI
RNV(2) = (STRESS(2) - R5*STRESS(1)) / PHI

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

RNV(1) = STRESS(1) / PHI ! TO BE CHECKED ...

END SELECT

*

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (STRESS(4)) IS PRESENT

RNV(4) = R3*STRESS(4) / PHI

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

RNV(4) = R3*STRESS(4) / PHI
RNV(5) = R3*STRESS(5) / PHI
RNV(6) = R3*STRESS(6) / PHI

END SELECT

*
* compute hardening modulus

HR = QR1*CR1*EXP(-CR1*P) + QR2*CR2*EXP(-CR2*P)

*
* compute denominator in the residual derivative of dlambda

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

>

DENOMCP = C1*(RNV(1)*RNV(1) + RNV(2)*RNV(2) + RNV(3)*RNV(3))
+ R2*C4*(RNV(1)*RNV(2) + RNV(2)*RNV(3) + RNV(3)*RNV(1))

74

*

*

CASE (1)

! PLANE STRESS CONDITION (RNV(3)=0)

DENOMCP = C1*(RNV(1)*RNV(1) + RNV(2)*RNV(2))

>

+ R2*C4*(RNV(1)*RNV(2))

CASE (2)

! UNIAXIAL STRESS CONDITION (RNV(2)=RNV(3)=0)

DENOMCP = C1*(RNV(1)*RNV(1)) ! TO BE CHECKED ...

END SELECT

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (RNV(4)) IS PRESENT

DENOMCP = DENOMCP + R4*C7*(RNV(4)*RNV(4))

CASE (3)

! ALL THREE SHEAR STRESSES (RNV(4:6)) ARE PRESENT

DENOMCP = DENOMCP + R4*C7*(RNV(4)*RNV(4) + RNV(5)*RNV(5)

>

END SELECT

+ RNV(6)*RNV(6))

DENOMV = (R1 + DLAMBDA / (PDOT*DT))**C
DENOMDV = (C / (PDOT*DT))*(R1 + DLAMBDA / (PDOT*DT))**(C - R1)
IF (TEMP_SOFT) THEN

TSTAR = (T - TR) / (TM - TR)
TSOFT = R1 - TSTAR**M
DENOM = DENOMCP + HR*DENOMV*TSOFT + SIGMAY*DENOMDV*TSOFT
ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1

DENOM = DENOMCP + HR*DENOMV + SIGMAY*DENOMDV

ENDIF

*
* compute residual increment ddlambda and update variable dlambda

DDLAMBDA = F / DENOM
DLAMBDA = DLAMBDA + DDLAMBDA

*
* update internal variables dependent on dlambda

P = P + DDLAMBDA
R = RP + HR*DLAMBDA
SIGMAY = SY + R

*
* update stress components using the elastic (hooke’s) law

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

STRESS(1) = STRESS(1) -

DDLAMBDA*(C1*RNV(1) + C4*RNV(2) + C5*RNV(3))

STRESS(2) = STRESS(2) -

DDLAMBDA*(C4*RNV(1) + C2*RNV(2) + C6*RNV(3))

STRESS(3) = STRESS(3) -

DDLAMBDA*(C5*RNV(1) + C6*RNV(2) + C3*RNV(3))

CASE (1)

! PLANE STRESS CONDITION (RNV(3)=0)

STRESS(1) = STRESS(1) -

DDLAMBDA*(C1*RNV(1) + C4*RNV(2))

STRESS(2) = STRESS(2) -

DDLAMBDA*(C4*RNV(1) + C2*RNV(2))

>

>

>

>

>

CASE (2)

! UNIAXIAL STRESS CONDITION (RNV(2)=RNV(3)=0)

STRESS(1) = STRESS(1) - DDLAMBDA*C1*RNV(1) ! TO BE CHECKED ...

END SELECT

*

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (RNV(4)) IS PRESENT

STRESS(4) = STRESS(4) - DDLAMBDA*C7*RNV(4)

CASE (3)

! ALL THREE SHEAR STRESSES (RNV(4:6)) ARE PRESENT

STRESS(4) = STRESS(4) - DDLAMBDA*C7*RNV(4)
STRESS(5) = STRESS(5) - DDLAMBDA*C8*RNV(5)
STRESS(6) = STRESS(6) - DDLAMBDA*C9*RNV(6)

END SELECT

*
* compute updated equivalent stress

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)
PHI = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)
+ STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
- STRESS(2)*STRESS(3) - STRESS(3)*STRESS(1)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

PHI = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)

- STRESS(1)*STRESS(2)

>
>

>

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

PHI = STRESS(1)*STRESS(1)

END SELECT

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT

75

CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (STRESS(4)) IS PRESENT

PHI = PHI + R3*(STRESS(4)*STRESS(4))

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

>
>

PHI = PHI + R3*(STRESS(4)*STRESS(4)
+ STRESS(5)*STRESS(5)
+ STRESS(6)*STRESS(6))

*

END SELECT

PHI = SQRT (PHI)

*
* compute yield function and new residual

IF (TEMP_SOFT) THEN

TSTAR = (T - TR) / (TM - TR) ! THIS RECALCULATION OF TSTAR
TSOFT = R1 - TSTAR**M
F = PHI - SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT

! AND TSOFT SEEMS REDUNDANT HERE

*
* check convergence

RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT)
RESNOR = ABS (RESNOR)

ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1

F = PHI - SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)

* check convergence

RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C))
RESNOR = ABS (RESNOR)

ENDIF

*

IF (RESNOR <= TOL) THEN

*
* update stress components based on plastic step
*-----------------------------------------------
STRESSNEW(1) = STRESS(1)
STRESSNEW(2:4) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

STRESSNEW(2) = STRESS(2)
STRESSNEW(3) = STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

STRESSNEW(2) = STRESS(2)

CASE (2)
END SELECT

*

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)
! ONLY TAU_XY (STRESS(4)) IS PRESENT
CASE (1)

STRESSNEW(4) = STRESS(4)

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

STRESSNEW(4) = STRESS(4)
STRESSNEW(5) = STRESS(5)
STRESSNEW(6) = STRESS(6)

END SELECT

*
* update remaining variables dependent of dlambda
*------------------------------------------------
* update temperature (only if temperature softening is activated)

IF (TEMP_SOFT) THEN

T = T + TQ*PHI*DLAMBDA / (RHO*CP)
IF (T >= TM) THEN ! CHECK IF MELTING OCCURS

STATE(13) = TM
D = DC ! this will cause the GP to fail (see 990)
GO TO 990

ENDIF

ENDIF

*
* compute hydrostatic stress based on new stress state,

SIGMAH = STRESS(1)
SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

SIGMAH = SIGMAH + STRESS(2) + STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

SIGMAH = SIGMAH + STRESS(2)

CASE (2)
END SELECT

*

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

SIGMAH = SIGMAH / R3

*
* 3d stress state and compute deviatoric stress

76

S1 = STRESS(1) - SIGMAH
S2 = STRESS(2) - SIGMAH
S3 = STRESS(3) - SIGMAH

*
* finding principal stress based on the deviatoric stress and invariants

A1 = R5*(S1*S1 + S2*S2 + S3*S3)

*

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)

A2 = -S1*S2*S3 ! DETERMINANT OF THE DEVIATORIC STRESS TENSOR

CASE (1)

! ONLY TAU_XY (STRESS(4)) IS PRESENT

A1 = A1 + STRESS(4)*STRESS(4)
A2 = S3*STRESS(4)*STRESS(4) - S1*S2*S3

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

>
>

>
>
>

A1 = A1 +

STRESS(4)*STRESS(4) + STRESS(5)*STRESS(5) +
STRESS(6)*STRESS(6)
A2 = S1*STRESS(5)*STRESS(5) +
S2*STRESS(6)*STRESS(6) +
S3*STRESS(4)*STRESS(4) -
S1*S2*S3 - R2*STRESS(4)*STRESS(5)*STRESS(6)

END SELECT

A3 = -SQRT (R27 / A1)*A2*R5 / A1
A3 = SIGN (MIN (ABS (A3), R1), A3)
A4 = ACOS (A3) / R3

*
* compute (chronological) principal stresses
AUX = SQRT (A1 / R3)
S1 = SIGMAH + R2*AUX*COS(A4)
S2 = SIGMAH + R2*AUX*COS(A4-R2*PI/R3) ! only for postprocessing
S3 = SIGMAH + R2*AUX*COS(A4+R2*PI/R3) ! only for postprocessing

*

WE = WE + MAX (S1, R0)*DLAMBDA
D = WE / WC

*
* update state variables before exit

STATE(1) = SIGMAH
STATE(2) = PHI
STATE(3) = P
STATE(4) = PHITRIAL
STATE(5) = F
STATE(6) = R
STATE(7) = DDLAMBDA
STATE(8) = DLAMBDA
STATE(9) = NRITER
STATE(10) = DLAMBDA / DT
STATE(11) = D

STATE(13) = T
STATE(14) = WE
STATE(15) = CSON

STATE(16) = S1
STATE(17) = S2
STATE(18) = S3

EXIT LOOP_NR

ELSE ! (RESNOR >= TOL)

IF (NRITER == MXITER) THEN

CALL ERRMSS (’VPJC’,’MAX NO. OF N-R ITERATIONS REACHED’)
STOP ’VPJC - MAX NO. OF N-R ITERATIONS REACHED’

ENDIF

ENDIF ! (RESNOR < TOL)

END DO LOOP_NR

ELSE ! (F <= R0)

*

*

*

*

*

*
* if not plasticity, elastic update

STRESSNEW(1) = STRESS(1)
STRESSNEW(2:4) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

STRESSNEW(2) = STRESS(2)
STRESSNEW(3) = STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

77

STRESSNEW(2) = STRESS(2)

CASE (2)
END SELECT

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)
! ONLY TAU_XY (STRESS(4)) IS PRESENT
CASE (1)

STRESSNEW(4) = STRESS(4)

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

STRESSNEW(4) = STRESS(4)
STRESSNEW(5) = STRESS(5)
STRESSNEW(6) = STRESS(6)

END SELECT

SIGMAH = STRESS(1)
SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

SIGMAH = SIGMAH + STRESS(2) + STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

SIGMAH = SIGMAH + STRESS(2)

CASE (2)
END SELECT

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

*

*

*

SIGMAH = SIGMAH / R3

*
* 3d stress state and compute deviatoric stress

S1 = STRESS(1) - SIGMAH
S2 = STRESS(2) - SIGMAH
S3 = STRESS(3) - SIGMAH

*
* finding principal stress based on the deviatoric stress and invariants

A1 = R5*(S1*S1 + S2*S2 + S3*S3)

*

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)

A2 = -S1*S2*S3 ! DETERMINANT OF THE DEVIATORIC STRESS TENSOR

CASE (1)

! ONLY TAU_XY (STRESS(4)) IS PRESENT

A1 = A1 + STRESS(4)*STRESS(4)
A2 = S3*STRESS(4)*STRESS(4) - S1*S2*S3

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

A1 = A1 +

>
>

>
>
>

STRESS(4)*STRESS(4) + STRESS(5)*STRESS(5) +
STRESS(6)*STRESS(6)
A2 = S1*STRESS(5)*STRESS(5) +
S2*STRESS(6)*STRESS(6) +
S3*STRESS(4)*STRESS(4) -
S1*S2*S3 - R2*STRESS(4)*STRESS(5)*STRESS(6)

END SELECT

A3 = -SQRT (R27 / A1)*A2*R5 / A1
A3 = SIGN (MIN (ABS (A3), R1), A3)
A4 = ACOS (A3) / R3

*
* compute (chronological) principal stresses

AUX = SQRT (A1 / R3)
S1 = SIGMAH + R2*AUX*COS(A4)
S2 = SIGMAH + R2*AUX*COS(A4-R2*PI/R3) ! only for postprocessing
S3 = SIGMAH + R2*AUX*COS(A4+R2*PI/R3) ! only for postprocessing

*

*

STATE(1) = SIGMAH ! COMPUTE SIGMAH AND INSERT IT HERE
STATE(2) = PHITRIAL ! COMPUTE PHI AND INSERT IT HERE
STATE(3) = PN ! COMPUTE P AND INSERT IT HERE
STATE(4) = PHITRIAL
STATE(5) = F
STATE(6) = RP
STATE(7) = R0
STATE(8) = R0
STATE(9) = R0
STATE(10) = R0
STATE(11) = D
STATE(12) = R1
STATE(13) = T
STATE(14) = WE
STATE(15) = CSON

STATE(16) = S1
STATE(17) = S2

78

*

*

*

STATE(18) = S3

ENDIF ! (F > R0)

GO TO 990

200 CONTINUE

*
*=======================================================================
* solu == 2 : radial return solution strategy
*=======================================================================
*
C COMPUTE THE GRADIENT OF THE VON MISES YIELD SURFACE

RNV(2:6) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

RNV(1) = (STRESS(1) - R5*(STRESS(2) + STRESS(3))) / PHI
RNV(2) = (STRESS(2) - R5*(STRESS(3) + STRESS(1))) / PHI
RNV(3) = (STRESS(3) - R5*(STRESS(1) + STRESS(2))) / PHI

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

RNV(1) = (STRESS(1) - R5*STRESS(2)) / PHI
RNV(2) = (STRESS(2) - R5*STRESS(1)) / PHI

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

RNV(1) = STRESS(1) / PHI ! TO BE CHECKED ...

END SELECT

*

*

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (STRESS(4)) IS PRESENT

RNV(4) = R3*STRESS(4) / PHI

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

RNV(4) = R3*STRESS(4) / PHI
RNV(5) = R3*STRESS(5) / PHI
RNV(6) = R3*STRESS(6) / PHI

END SELECT

IF (F > R0) THEN

*
C PLASTIC STEP: APPLY RADIAL RETURN MAPPING - USE NEWTON-RAPHSON
C ALGORITHM TO SOLVE THE RETURN MAPPING EQUATION
C--------------------------------------------------------------

P = PN
LOOP_NR2 : DO NRITER = 1, MXITER

*
C COMPUTE HARDENING MODULUS

HR = QR1*CR1*EXP(-CR1*P) + QR2*CR2*EXP(-CR2*P)

C COMPUTE RESIDUAL DERIVATIVE

DENOMV = (R1 + DLAMBDA/(PDOT*DT))**C
DENOMDV = (C/(PDOT*DT))*(R1+DLAMBDA/(PDOT*DT))**(C-R1)
IF (TEMP_SOFT) THEN

TSTAR = (T - TR) / (TM - TR)
TSOFT = R1 - TSTAR**M

ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1

TSTAR = R0
TSOFT = R1

ENDIF
DENOM = R3G + HR*DENOMV*TSOFT + SIGMAY*DENOMDV*TSOFT
C COMPUTE NEWTON-RAPHSON INCREMENT AND UPDATE VARIABLE DLAMBDA

DDLAMBDA = F / DENOM
DLAMBDA = DLAMBDA + DDLAMBDA

C COMPUTE NEW RESIDUAL

P = P + DDLAMBDA
R = RP + HR*DLAMBDA
SIGMAY = SY + R
DENOMV = (R1 + DLAMBDA/(PDOT*DT))**C
F = PHITRIAL - R3G*DLAMBDA - SIGMAY*DENOMV*TSOFT

C CHECK CONVERGENCE

RESNOR = F / (SIGMAY*DENOMV*TSOFT)
RESNOR = ABS (RESNOR)

*

IF (RESNOR <= TOL) THEN

*
C UPDATE VARIABLES DEPENDENT OF DLAMBDA
C--------------------------------------------------------------
C UPDATE TEMPERATURE (ONLY IF TEMPERATURE SOFTENING IS ACTIVATED)

79

IF (TEMP_SOFT) THEN

PHIEFF = PHITRIAL - R3G*DLAMBDA ! EQUIVALENT (EFFECTIVE) STRESS
T = T + TQ*PHIEFF*DLAMBDA / (RHO*CP)
IF (T >= TM) THEN ! CHECK IF MELTING OCCURS

STATE(13) = TM
D = DC ! this will cause the GP to fail (see 990)
GO TO 990

ENDIF

ENDIF

C COMPUTE HYDROSTATIC STRESS BASED ON TRIAL STRESS STATE,
C I.E. UTILIZE THAT THE HYDROSTATIC STRESS IS THE SAME
C FOR BOTH TRIAL STATE AND FINAL STATE. ONLY THE DEVIATORIC
C STRESS STATE WILL CHANGE.

SIGMAH = STRESS(1)
SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

SIGMAH = SIGMAH + STRESS(2) + STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

SIGMAH = SIGMAH + STRESS(2)

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

END SELECT
SIGMAH = SIGMAH / R3

C COMPUTE DEVIATORIC STRESS

S1 = STRESS(1) - SIGMAH
S2 = STRESS(2) - SIGMAH
S3 = STRESS(3) - SIGMAH

C FINDING PRINCIPAL STRESS BASED ON THE DEVIATORIC STRESS AND INVARIANTS

A1 = R5*(S1*S1 + S2*S2 + S3*S3)
SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)

A2 = -S1*S2*S3 ! DETERMINANT OF THE DEVIATORIC STRESS TENSOR

CASE (1)

! ONLY TAU_XY (STRESS(4)) IS PRESENT

A1 = A1 + STRESS(4)*STRESS(4)
A2 = S3*STRESS(4)*STRESS(4) - S1*S2*S3

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

>
>

>
>
>

A1 = A1 +

STRESS(4)*STRESS(4) + STRESS(5)*STRESS(5) +
STRESS(6)*STRESS(6)
A2 = S1*STRESS(5)*STRESS(5) +
S2*STRESS(6)*STRESS(6) +
S3*STRESS(4)*STRESS(4) -
S1*S2*S3 - R2*STRESS(4)*STRESS(5)*STRESS(6)

END SELECT
A3 = -SQRT (R27 / A1)*A2*R5 / A1
A3 = SIGN (MIN (ABS (A3), R1), A3)
A4 = ACOS (A3) / R3
C COMPUTE (CHRONOLOGICAL) PRINCIPAL STRESSES
AUX = SQRT (A1 / R3)
S1 = SIGMAH + R2*AUX*COS(A4)
S2 = SIGMAH + R2*AUX*COS(A4-R2*PI/R3) ! only for postprocessing
S3 = SIGMAH + R2*AUX*COS(A4+R2*PI/R3) ! only for postprocessing
WE = WE + MAX (S1, R0)*DLAMBDA
D = WE / WC

*
C UPDATE STRESS COMPONENTS USING THE ELASTIC (HOOKE’S) LAW
C--------------------------------------------------------------

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

STRESSNEW(1) = STRESS(1) -

DLAMBDA*(C1*RNV(1) + C4*RNV(2) + C5*RNV(3))

STRESSNEW(2) = STRESS(2) -

DLAMBDA*(C4*RNV(1) + C2*RNV(2) + C6*RNV(3))

STRESSNEW(3) = STRESS(3) -

CASE (1)

! PLANE STRESS CONDITION (RNV(3)=0)

DLAMBDA*(C5*RNV(1) + C6*RNV(2) + C3*RNV(3))

>

>

>

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

>

>

STRESSNEW(1) = STRESS(1) -

DLAMBDA*(C1*RNV(1) + C4*RNV(2))

STRESSNEW(2) = STRESS(2) -

DLAMBDA*(C4*RNV(1) + C2*RNV(2))

CASE (2)

! UNIAXIAL STRESS CONDITION (RNV(2)=RNV(3)=0)

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

STRESSNEW(1) = STRESS(1) -

>

DLAMBDA*C1*RNV(1) ! TO BE CHECKED ...

END SELECT

80

*

SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
CASE (0)
CASE (1)

! NO SHEAR STRESS IS PRESENT
! ONLY TAU_XY (RNV(4)) IS PRESENT

STRESSNEW(4) = STRESS(4) - DLAMBDA*C7*RNV(4)

CASE (3)

! ALL THREE SHEAR STRESSES (RNV(4:6)) ARE PRESENT

STRESSNEW(4) = STRESS(4) - DLAMBDA*C7*RNV(4)
STRESSNEW(5) = STRESS(5) - DLAMBDA*C8*RNV(5)
STRESSNEW(6) = STRESS(6) - DLAMBDA*C9*RNV(6)

END SELECT

*
C COMPUTE UPDATED EQUIVALENT STRESS

SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
! NO PLANE STRESS CONDITION
CASE (0)

>
>
>
>
>

>
>

>
>

*

*

PHI = STRESSNEW(1)*STRESSNEW(1) +
STRESSNEW(2)*STRESSNEW(2) +
STRESSNEW(3)*STRESSNEW(3) -
STRESSNEW(1)*STRESSNEW(2) -
STRESSNEW(2)*STRESSNEW(3) -
STRESSNEW(3)*STRESSNEW(1)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

PHI = STRESSNEW(1)*STRESSNEW(1) +
STRESSNEW(2)*STRESSNEW(2) -
STRESSNEW(1)*STRESSNEW(2)

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

PHI = STRESSNEW(1)*STRESSNEW(1)

END SELECT

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)
! ONLY TAU_XY (STRESSNEW(4)) IS PRESENT
CASE (1)

PHI = PHI + R3*(STRESSNEW(4)*STRESSNEW(4))

CASE (3)

! ALL THREE SHEAR STRESSES (STRESSNEW(4:6)) ARE PRESENT

PHI = PHI + R3*(STRESSNEW(4)*STRESSNEW(4)
+ STRESSNEW(5)*STRESSNEW(5)
+ STRESSNEW(6)*STRESSNEW(6))

END SELECT

PHI = SQRT (PHI)

*
C UPDATE STATE VARIABLES (SDV) BEFORE EXIT
STATE(1) = SIGMAH
STATE(2) = PHI
STATE(3) = P
STATE(4) = PHITRIAL
STATE(5) = F
STATE(6) = R
STATE(7) = DDLAMBDA
STATE(8) = DLAMBDA
STATE(9) = NRITER
STATE(10) = DLAMBDA / DT
STATE(11) = D

*

*

*

*

*

STATE(13) = T
STATE(14) = WE
STATE(15) = CSON

STATE(16) = S1
STATE(17) = S2
STATE(18) = S3

EXIT LOOP_NR2

ELSE ! (RESNOR >= TOL)

IF (NRITER == MXITER) THEN

CALL ERRMSS (’VPJC’,’MAX NO. OF N-R ITERATIONS REACHED’)
STOP ’VPJC - MAX NO. OF N-R ITERATIONS REACHED’

ENDIF

ENDIF ! (RESNOR < TOL)

END DO LOOP_NR2

ELSE ! (F <= R0)

*
C IF NOT PLASTICITY, ELASTIC UPDATE

STRESSNEW(1) = STRESS(1)
STRESSNEW(2:4) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST

81

CASE (0)

! NO PLANE STRESS CONDITION

STRESSNEW(2) = STRESS(2)
STRESSNEW(3) = STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

STRESSNEW(2) = STRESS(2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)
* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

CASE (2)

END SELECT

*

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)
! ONLY TAU_XY (STRESS(4)) IS PRESENT
CASE (1)

STRESSNEW(4) = STRESS(4)

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

STRESSNEW(4) = STRESS(4)
STRESSNEW(5) = STRESS(5)
STRESSNEW(6) = STRESS(6)

END SELECT

SIGMAH = STRESS(1)
SELECT CASE (IPLANC)
CASE (0)

! NO PLANE STRESS CONDITION

SIGMAH = SIGMAH + STRESS(2) + STRESS(3)

CASE (1)

! PLANE STRESS CONDITION (STRESS(3)=0)

* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

SIGMAH = SIGMAH + STRESS(2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)
* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)

CASE (2)

END SELECT
SIGMAH = SIGMAH / R3

*
* 3d stress state and compute deviatoric stress

S1 = STRESS(1) - SIGMAH
S2 = STRESS(2) - SIGMAH
S3 = STRESS(3) - SIGMAH

*
* finding principal stress based on the deviatoric stress and invariants

A1 = R5*(S1*S1 + S2*S2 + S3*S3)

*

*

SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
! NO SHEAR STRESS IS PRESENT
CASE (0)

A2 = -S1*S2*S3 ! DETERMINANT OF THE DEVIATORIC STRESS TENSOR

CASE (1)

! ONLY TAU_XY (STRESS(4)) IS PRESENT

A1 = A1 + STRESS(4)*STRESS(4)
A2 = S3*STRESS(4)*STRESS(4) - S1*S2*S3

CASE (3)

! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT

A1 = A1 +

>
>

>
>
>

STRESS(4)*STRESS(4) + STRESS(5)*STRESS(5) +
STRESS(6)*STRESS(6)
A2 = S1*STRESS(5)*STRESS(5) +
S2*STRESS(6)*STRESS(6) +
S3*STRESS(4)*STRESS(4) -
S1*S2*S3 - R2*STRESS(4)*STRESS(5)*STRESS(6)

END SELECT

A3 = -SQRT (R27 / A1)*A2*R5 / A1
A3 = SIGN (MIN (ABS (A3), R1), A3)
A4 = ACOS (A3) / R3

*
* compute (chronological) principal stresses

AUX = SQRT (A1 / R3)
S1 = SIGMAH + R2*AUX*COS(A4)
S2 = SIGMAH + R2*AUX*COS(A4-R2*PI/R3) ! only for postprocessing
S3 = SIGMAH + R2*AUX*COS(A4+R2*PI/R3) ! only for postprocessing

*

STATE(1) = SIGMAH ! compute SIGMAH and insert it here
STATE(2) = PHITRIAL ! compute PHI and insert it here
STATE(3) = PN ! compute P and insert it here
STATE(4) = PHITRIAL
STATE(5) = F
STATE(6) = RP
STATE(7) = R0
STATE(8) = R0
STATE(9) = R0
STATE(10) = R0
STATE(11) = D

82

STATE(12) = R1
STATE(13) = T
STATE(14) = WE
STATE(15) = CSON

STATE(16) = S1
STATE(17) = S2
STATE(18) = S3

ENDIF ! (F > R0)

*

*

*
* decide if material point stiffness is zero, i.e. failed
*--------------------------------------------------------------

990 IF (D >= DC) THEN

*
* the current gp was already failed or is failing now :
* state(1:10,12,13), the relevant components
* of stressnew and cson are set to 0.
* the total strains (strain) are not set to 0 so they continue to be
* updated as long as the element is not eroded, i.e. as long
* as the element has some virgin gauss points.
*

SELECT CASE (ITAU)
CASE (0)

STRESSNEW(1:3) = R0

CASE (1)

STRESSNEW(1:4) = R0

CASE (2)

STRESSNEW(1:5) = R0

CASE (3)

STRESSNEW(1:6) = R0

END SELECT
STATE(1:10) = R0

* state(11) (damage d) keeps last value it had before failure

STATE(12) = R0

* state(13) (absolute temperature t) keeps last value it had before failure
* state(14) (damage accumulation we) keeps last value it had before failure

CSON = R0
STATE(15) = CSON

ELSE

*
* no failure
*

STATE(12) = R1

*
* strain increment(s) and total strain(s) corresponding to
* the imposed zero stress condition(s)
*

SELECT CASE (IPLANC)
CASE (0)
CASE (1)

! NO PLANE STRESS CONDITION
! PLANE STRESS CONDITION (STRESS(3)=0)

* store provisional strain increment

AUX3 = STRAININC(3)

* compute and update the actual strain increment

FACTOR = T2 / YOUNG
STRAININC(3) = FACTOR*(STRESSNEW(1) - STRESSOLD(1)

>
>

+ STRESSNEW(2) - STRESSOLD(2))
- STRAININC(1) - STRAININC(2)

* update the total strain

STRAIN(3) = STRAIN(3) - AUX3 + STRAININC(3)

CASE (2)

! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)

AUX2 = STRAININC(2)
AUX3 = STRAININC(3)
FACTOR = T2 / YOUNG
STRAININC(2) = 0.5D0*(FACTOR*(STRESSNEW(1) - STRESSOLD(1))

>

- STRAININC(1))

STRAININC(3) = STRAININC(2)
STRAIN(2) = STRAIN(2) - AUX2 + STRAININC(2)
STRAIN(3) = STRAIN(3) - AUX3 + STRAININC(3)

END SELECT

ENDIF

*
* erosion (optional)
*

IF (L_EROSION) THEN

IF (D >= DC) FAIL_FLAG = .TRUE.

ENDIF

83

*

END SUBROUTINE VPJC

*===============================================================================

END MODULE M_MATERIAL_VPJC

84

Europe Direct is a service to help you find answers to your questions about the European Union 
Free phone number (*): 00 800 6 7 8 9 10 11 
(*) Certain mobile telephone operators do not allow access to 00 800 numbers or these calls may be billed. 

A great deal of additional information on the European Union is available on the Internet. 
It can be accessed through the Europa server http://europa.eu 

How to obtain EU publications 

Our publications are available from EU Bookshop (http://bookshop.europa.eu), 
where you can place an order with the sales agent of your choice. 

The Publications Office has a worldwide network of sales agents. 
You can obtain their contact details by sending a fax to (352) 29 29-42758. 

 
 
 
 
 
 
 
 
 
 
 
 
 
L
B
-
N
A
-
2
7
9
8
2
-
E
N
-
N

JRC Mission 

As the Commission’s  
in-house science service,  
the Joint Research Centre’s  
mission is to provide EU  
policies with independent,  
evidence-based scientific  
and technical support  
throughout the whole  
policy cycle. 

Working in close  
cooperation with policy  
Directorates-General,  
the JRC addresses key  
societal challenges while  
stimulating innovation  
through developing  
new methods, tools  
and standards, and sharing  
its know-how with  
the Member States,  
the scientific community  
and international partners. 

Serving society  
Stimulating innovation  
Supporting legislation 

doi:10.2788/609529 

ISBN 978-92-79-59746-6 

2