Johnson Cook Fortran Code

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

a Centre for Advanced Structural Analysis (CASA), NTNU, Norwegian University of Science and Technology, NO-7491
Trondheim, Norway
b Structural Impact Laboratory (SIMLab), Department of Structural Engineering, NTNU, NO-7491 Trondheim, Norway
c European 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 flow rule, nonlinear
isotropic strain hardening, strain-rate hardening, temperature softening due to adiabatic heating, and duc-
tile failure. The performance of the implementation is verified and validated in EPX. First, a verification
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 finite 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 difference 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
∗ Correspondingauthor
Email address: [email protected] (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 influenced by the
accuracy, robustness and effectiveness 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
flow 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 fit 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 defined 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 defined 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 finally verified and validated through numerical simulations
in EPX. This is first done by a verification 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 infinitesimal 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 infinitesimal deformations,
the strain rate tensor ε̇ may be decomposed by an additive split
ε̇ = ε̇e + ε̇p (1)

e p
where ε̇ is the elastic and ε̇ is the plastic strain rate tensor, respectively. Using the generalized Hooke’s
law, the rate of the stress tensor is defined as
σ̇ = C : ε̇e (2)
where C is the fourth order elasticity tensor given by
C = λel I ⊗ I + 2µel I (3)

where I is the second order identity tensor, and λel and µel are the Lamè constants given as
νE E
λel = and µel = (4)
(1 + ν)(1 − 2ν) 2(1 + ν)
The elasticity tensor C is assumed to be constant, E is the Young’s modulus and ν the Poisson’s ratio.
In general, assuming only isotropic hardening, the yield criterion is defined as
f (σ, R) = ϕ(σ) − σy (σ0 , R) = σeq − (σ0 + R) (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 flow 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
ϕ(aσ) = aϕ(σ) (6)

where a could be any positive scalar. Then, Euler’s theorem for homogeneous functions gives [10]
∂ϕ
σ: =ϕ (7)
∂σ
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
Dp = σ ε̇p ≥ 0 (8)
where Dp expresses the plastic dissipation per unit volume, which is assumed to be dissipated as heat.
The restriction in Eq.(8) is fulfilled by defining the plastic strain rate tensor ε̇p in a way that ensures
non-negative dissipation. This is done by the plastic flow rule, which in the general case is expressed as
∂g
ε̇p = λ̇ (9)
∂σ
3
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
≥0
σ: (10)
∂σ
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 fulfills the dissipation inequality in
Eq.(8) as
∂g
Dp = σ ε̇p = λ̇σ : = λ̇g ≥ 0 (11)
∂σ
By further assuming that the plastic potential function g is defined by the yield function f , Eq.(9) is expressed
as
∂f
ε̇p = λ̇ (12)
∂σ
and it is seen that the plastic potential function is associated with the yield function f . Eq.(12) is therefore
called the associated flow rule, and ensures non-negative dissipation through Eqs.(5), (7) and (8) as
∂f ∂ϕ
Dp = σ ε̇p = λ̇σ : = λ̇σ : = λ̇ϕ = λ̇σeq ≥ 0 (13)
∂σ ∂σ
This implies, that in the case of the associated flow rule, the equivalent plastic strain rate equals the plastic
parameter
ṗ = λ̇ (14)
where the equivalent plastic strain p is defined 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
Dp = σ ε̇p = σeq ṗ (15)

p p
Hence, assuming the associated flow rule implies that the incremental plastic strain vector dε = ε̇ dt is
parallel to the gradient of the yield surface ∇f at σ, and therefore directed along the outward normal n of
the yield surface (see Figure 1). Moreover, the magnitude of ε̇p is given by the plastic multiplier λ̇.
An important remark to this argumentation is that the shape of the yield surface also determines the direc-
tion of the plastic flow, as well as the stress state at which yielding initiates.
The associated flow 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
(σ − σ ∗ ) ε̇p ≥ 0 (16)
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 flow, and in addition it is seen that the yield surface must be convex [11]. Thus, the use of an associ-
ated flow rule ensures that the plastic strain increments dεp are vectors perpendicular to the potential surface.
To distinguish between elastic loading/unloading and plastic loading, the Kuhn-Tucker conditions are used
to determine when plastic flow may occur. These conditions determines the evolution of the yield surface
and are given as [11]
4
σ2 ∇f = ∂f
∂σ e
σ
f (σ) = 0
σ1
∂f
dεp = dλ ∂σ = dλn
Figure 1: Geometric representation of the associated flow 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 first 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]
λ̇f˙ = 0 (18)
The consistency condition ensures that the stress must be located on the yield surface during continuous
plastic flow, and is, in the theory of rate-independent plasticity, used to determined the plastic multiplier λ̇.
To account for viscoplasticity and strain-rate sensitivity, the theory of plasticity has to be extended. The
main difference 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 ⇐⇒ f (σ) = σv (p, ṗ) = (υp − 1)σy > 0 (20)
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

0 for f ≤ 0
ṗ = λ̇ = (21)
φ(σ, R) for f > 0
5
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 flow depends explicitly
on time during the evolution of λ̇.
σ3
σeq
π-plane
σv (ṗ, p) σ
R(p) f =0 f >0
σ0
σ1 σ2
p
(a) Equivalent stress vs equivalent plastic strain (b) The deviatoric plane (π-plane)
Figure 2: Viscoplastic effects due to strain-rate sensitivity.
To account for large displacement and large strain, the constitutive model is defined in a co-rotational
coordinate system in which the basis system rotates with the material [7]. Typically, there are two different
approaches of dealing with finite deformations, i.e., the hypoelastic-plastic formulations and hyperelastic-
plastic formulations. In contrast to the hyperelastic formulation of finite 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 infinitesimal
deformations, and the rate-of-deformation tensor D is decomposed additively as [1]
D = De + Dp (22)
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.
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
σ̂ ∇ = C : De = λel (trDe ) I + 2µel De (23)

where σ is the Cauchy stress tensor, λel and µel are the Lamè constants given in Eq.(4), and σ̂ ∇ is an
objective (or frame-invariant) stress rate given as
σ̂ ∇ = σ̇ − Aσ + σA (24)
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 + λel (trDe ) I + 2µel De (25)

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
6
1
σ̂ ∇J = σ̇ − W σ + σW , L − LT

A=W = (26)
2
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 = σ̇ − ṘRT σ + σ ṘRT = σ̇ − Ωσ + σΩ, A = Ω = ṘRT = −RṘT (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.
x̂2
x̂1
x̂2
x̂1
x̂2
x̂1
x2 U
x1
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 (x̂1 , x̂2 ) represents the local (or
co-rotated) coordinate system.
From Eqs. (26) and (27) it is seen that the Green-Naghdi rate differs from the Jaumann rate only in using
a different measure of the rotation of the material. Despite the difference 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]
1
W = ṘRT − R U̇U−1 + U−T U̇ RT (28)
2
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 = ṘRT = Ω (29)
7
whereas in for instance shear one should be aware of the different formulation and how this may influence
the results of the FE simulation [1].
As concerns the plastic deformation rate tensor Dp this is defined by means of the associated flow rule as
∂f
Dp = λ̇ (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 defines a surface in the stress space which is defined 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.
1
e= √ [1 1 1]
σ3 3
yield surface σ3
Hydrostatic axis
π-plane
σ1 = σ2 = σ3
yield curve
·r Rv
e
v
r
R
P
σ2
q
π-plane 2
3 σy
(Deviatoric plane)
σ1 + σ2 + σ3 = 0 σ2
σ1 σ1
(a) Principal stress space (b) The deviatoric plane (π-plane)
Figure 4: The von Mises yield criterion.

q
2
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]
r · r − (e · r) (e · r) = R2v (31)
By splitting the stress in deviatoric and hydrostatic parts
0 1 σ1 + σ2 + σ3
σ = σ − σH I and σH = tr(σ) = (32)
3 3
0
where σ is the deviatoric stress and σH is the hydrostatic stress, the von Mises yield criterion could be
expressed as
2
1 0 0 0 0 0 0 2 2
R2v = σ: σ − √ tr(σ) 2
= (σij + σH δij )(σij + σH δij ) − 3σH = σij σij = σ : σ = σ (33)
3 3 y
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
r
3 0 p
f (σ) = σeq − σy = σ : σ 0 − σy = 3J2 − σy (34)
2
where the von Mises equivalent stress σeq , in terms of the stress components, reads
q
2 + σ2 + σ2 − σ σ 2 2 2
σeq = σ11 22 33 11 22 − σ22 σ33 − σ33 σ11 + 3(σ12 + σ23 + σ31 ) (35)
and the associated flow rule in Eq.(30) is then given as
0
∂f 3 σ
Dp = λ̇ = λ̇ (36)
∂σ 2 σeq
The rate of the accumulated plastic strain ṗ, in the particular case of associated flow and von Mises isotropic
behaviour, is defined as
r
2 p
ṗ = D : Dp (37)
3
2.2. The modified 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
σy = (A + Bpn ) (1 + Cln(ṗ∗ )) (1 − T ∗ m ) (38)

where A, B, C, n and m are material constants determined from material tests, and
ṗ T − Tr
ṗ∗ = and T∗ = (39)
ṗ0 Tm − Tr
represents the dimensionless strain rate ṗ∗ and the non-dimensional homologous temperature T ∗ , respec-
tively. The parameter ṗ0 is a user-defined reference strain rate, while Tr and Tm are defined as the room
and melting temperature, respectively. The first term in Eq.(38) represents the strain hardening, the sec-
ond term describes the influence 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 flow 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 effects when ṗ∗ < 1, the strain-rate-sensitivity term in the Johnson-Cook
model can be adjusted [13]. This leads to the modified Johnson-Cook model, given by [14]
C
σy = (A + Bpn ) (1 + ṗ∗ ) (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
normalized stress σy / [(A + Bpn ) (1 − T ∗ m )]. Note that, since the formulations in Eqs.(38) and (40) are
slightly different, the calibration of the viscoplastic material parameter C may be somewhat different in the
two versions of the Johnson-Cook model.
The first 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 infinite value of the hardening
modulus HR = Bnpn−1 when plasticity occurs. This could be avoided by defining the initial plastic strain
9
1.2
1.1
Normalized stress
1
0.9
Original Johnson−Cook
Modified Johnson−Cook
0.8 −5 −3 0 3 5
10 10 10 10 10
Dimensionless strain rate ṗ∗
Figure 5: Original vs modified Johnson-Cook model, when CJC = 0.014, CM JC = 0.013, and ṗ0 = 5e−4 .
as an infinite small value. However, it is desirable to formulate the material model as general as possible.
Therefore, the first term of Eq.(40), representing the strain hardening, is replaced by another hardening
rule than that proposed by Johnson and Cook. By replacing the first term of Eq.(40) with a two terms
saturation type of hardening rule proposed by Voce [15], the constitutive model could be expressed as
2
C C
σy = (A + R(p)) (1 + ṗ∗ ) (1 − T ∗ m ) = (A + Qk (1 − exp (−Ck p))) (1 + ṗ∗ ) (1 − T ∗ m )
X
(41)
k=1
where the constant A is the initial yield strength of the material, and Qk and Ck (k = 1, 2) define the
strain hardening. As before, the rate sensitivity is governed by the constant C, while m defines the thermal
softening behaviour. The hardening modulus HR is found from the rate of change in R(p) given as
2
dR dR dp dp dR X
Ṙ = = = HR =⇒ HR = = Qk Ck e(−Ck p) (42)
dt dp dt dt dp
k=1
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

 0
C1 for f ≤ 0
ṗ = λ̇ = ϕ(σ) (43)
 ṗ0 σy (σ0 ,R,T ) −1 for f > 0
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 flow stress. This is
already introduced in the last term of Eq.(41), which accounts for the thermal softening of the flow 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
Dp = χσ : Dp = χσeq ṗ = ρCp Ṫ (44)
where ρ is the material density, Cp is the specific heat capacity of the solid material and χ is the Taylor-
Quinney coefficient. The temperature increase due to adiabatic heating can be computed as
Dp
Ṫ = (45)
ρCp
From Eq.(44) it is seen that the Taylor-Quinney coefficient χ 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 ”fields” 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 defined as the first 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.,
Z p
W 1
D= = hσ1 i dp (46)
Wc Wc 0
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-
sion hσ1 i is equivalent to the function max(0, σ1 ), which implies that only positive values of the maximum
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.
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.,
!
3 + µσ ∗ 3 + µσ
σ1 = σH + p σeq = σ + p σeq (47)
3 3 + µ2σ 3 3 + µ2σ
where the stress triaxiality σ ∗ and the Lode parameter µσ is defined by
11
σH 2σ2 − σ1 − σ3
σ∗ = , µσ = (48)
σeq σ3 − σ1
Remember that σ1 ≥ σ2 ≥ σ3 are the ordered principal stresses. Eq.(46) can then be expressed as
Z p* +
1 ∗ 3 + µσ
D= σ + p σeq dp (49)
Wc 0 3 3 + µ2σ
Since material failure occurs when D = Dc = 1.0 and p = pf , Eq.(49) can be written as
Z pf * +
1 ∗ 3 + µσ
1= σ + p σeq dp (50)
Wc 0 3 3 + µ2σ
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
Wc 3 + µσ
σ∗ = P2 Qk
− p (51)
Apf + k=1 (Qk pf + Ck [exp (−Ck pf ) − 1]) 3 3 + µ2σ
The influence 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 different combinations of the stress
triaxiality and the Lode parameter. Some typical examples are (µσ , σ ∗ ) = (+1, 1/3) which implies uniaxial
tension (σ1 ≥ σ2 = σ3 ), (µσ , σ ∗ ) = (0, 0) corresponding to pure shear (2σ2 = σ1 + σ3 ) and (µσ , σ ∗ ) =
(−1, −1/3) representing uniaxial compression (σ1 = σ2 ≥ σ3 ).
2
7 µσ = −1
µσ = 0
1.75 µσ = +1
6
Plastic failure strain, pf
Plastic failure strain, pf
5 1.5
4 1.25
3 1
2 0.75
1
0.5
0
−1/3 −1 0.25
0 1/2
Stres 0
s tria 1/3 1/2 r, µσ
mete 0
xialit
y, σ ∗ 2/3 1
o de para −2/3 −1/3 0 1/3 2/3
L Stress triaxiality, σ ∗
(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]
0
which uses the principal deviatoric stresses σi to find the principal stresses σi . The principal deviatoric
0
stresses are determined from the eigenvalue problem for the deviatoric stress σ in Eq. (32), i.e.,
12
0 0
det(σ I − σ ) = 0 (52)
where det(·) is the determinant of the matrix argument (·). The result of this operation is a cubic equation
called the characteristic equation, given as
0 0 0 0
(σ )3 − I1 (σ )2 − J2 σ − J3 = 0 (53)
where the principal invariants are given respectively as
0 0 0 1 0 0 1 0 0 0 0
I1 = tr(σ ) = σii = 0, J2 = −I2 = tr(σ 2 ) = σij σij , J3 = I3 = det(σ ) (54)
2 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
r r r
0 J2 θ 0 J2 θ 2π 0 0 J2 θ 2π 0
σ1 = 2 cos , σ2 = 2 cos − ≤ σ1 , σ3 = 2 cos + ≤ σ2 (55)
3 3 3 3 3 3 3 3
where the angle θ is given as
J3
cos θ = p 3 , 0≤θ≤π (56)
2 J2 /27
and the Lode angle is recognized as θL = θ/3. Finally, the principal stresses σi of the stress tensor σ are
determined by
0 0 1
σi = σi + σH = σi + tr(σ) (57)
3
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 defined 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
    
ε11 1 −ν −ν 0 0 0 σ11
 ε22  −ν 1 −ν 0 0 0  σ22 
    
 ε33  1  −ν −ν 1 0 0 0  0 
ε=   = Sσ =     (58)
2ε
 12 
 E  0 0 0 2(1 + ν) 0 0  σ12 
 
 0   0 0 0 0 2(1 + ν) 0  0 
0 0 0 0 0 0 2(1 + ν) 0
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 1 −ν 0 σ11
1
ε =  ε22  = Sσ = −ν 1 0  σ22  (59)
E
2ε12 0 0 2(1 + ν) σ12
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 1 ν 0 ε11
E ν 1
σ = σ22  = Cε = 0   ε22  (60)
1 − ν2 1−ν
σ12 0 0 2 2ε 12
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
s 2
σ11 + σ22 σ11 − σ22 2
σ1 = + + σ12
2 2
s (61)
2
σ11 + σ22 σ11 − σ22 2
σ2 = − + σ12
2 2
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.,
1 2 E
εv = εkk = − P, K = λel + µel = (62)
K 3 3(1 − 2ν)
where P = − 31 σkk is the mean stress and K is the bulk modulus. Note that the hydrostatic pressure P is
equal to the hydrostatic stress σH but has the opposite sign.
The physical meaning of εkk is the relative change of volume

dV
εkk = ε11 + ε22 + ε33 = (63)
V
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 defining 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., εpkk = 0) the plastic strain increment through
the thickness (due to plastic thinning) is given as
∆εp33 = −(∆εp11 + ∆εp22 ) (65)

where ∆εp11 and ∆εp22 could be expressed using strain decomposition
∆εp11 = ∆ε11 − ∆εe11

(66)
∆εp22 = ∆ε22 − ∆εe22
and the elastic strain components are given by Hooke’s law in Eq.(59)
1
∆εe11 =
(∆σ11 − ν∆σ22 )
E (67)
1
∆εe22 = (∆σ22 − ν∆σ11 )
E
By combining Eqs.(65) to (67), the plastic strain increment through the thickness is given by
∆εp33 = −(∆ε11 − ∆εe11 ) − (∆ε22 − ∆εe22 )

= ∆εe11 + ∆εe22 − ∆ε11 − ∆ε22
1 (68)
= (∆σ11 − ν∆σ22 + ∆σ22 − ν∆σ11 ) − (∆ε11 + ∆ε22 )
E
1−ν
= (∆σ11 + ∆σ22 ) − (∆ε11 + ∆ε22 )
E
Then, when the plastic strain increment ∆εp33 is known the total through thickness strain increment ∆ε33
is finally found by strain decomposition
∆ε33 = ∆εe33 + ∆εp33

ν 1−ν
= − (∆σ11 + ∆σ22 ) + (∆σ11 + ∆σ22 ) − (∆ε11 + ∆ε22 ) (69)
E E
1 − 2ν
= (∆σ11 + ∆σ22 ) − (∆ε11 + ∆ε22 )
E
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
∆εp11 = ∆ε11 − ∆εe11

∆εp22 = ∆ε22 − ∆εe22 (70)
∆εp33 = ∆ε33 − ∆εe33
and that the plastic part of the deformation occurs at constant volume (dV = 0), due to plastic incompress-
ibility
∆εp11 + ∆εp22 + ∆εp33 = 0 (71)
15
In the general 3D case, the elastic strain increments are related to the total stress increments by Hooke’s
law
1
∆εe11 = [∆σ11 − ν (∆σ22 + ∆σ33 )]
E
1
∆εe22 = [∆σ22 − ν (∆σ11 + ∆σ33 )] (72)
E
1
∆εe33 = [∆σ33 − ν (∆σ11 + ∆σ22 )]
E
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
1
∆εe11 = ∆σ11
E
ν
∆εe22 = − ∆σ11 = −ν∆εe11 (73)
E
ν
∆εe33 = − ∆σ11 = −ν∆εe11
E
which implies that
∆εe22 = ∆εe33 = −ν∆εe11 (74)

and by utilizing the symmetry it is also assume that
∆εp22 = ∆εp33 (75)

and therefore
∆ε22 = ∆ε33 (76)
ε22 = ε33 (77)

Then, by replacing Eq. (75) into Eq. (71) one obtains
2∆εp22 = 2∆εp33 = −∆εp11 (78)
1 1
∆εp22 = ∆εp33 = − ∆εp11 = − (∆ε11 − ∆εe11 ) (79)
2 2
Using Eq. (73) this now becomes

1 1 1 1
∆εp22 = ∆εp33 = − ∆ε11 − ∆σ11 = − ∆ε11 + ∆σ11 (80)
2 E 2 2E
Finally, the total lateral strain increment ∆ε22 (and ∆ε33 ) is obtained by strain decomposition
∆ε22 = ∆ε33 = ∆εe22 + ∆εp22 = ∆εe33 + ∆εp33

1 1
= − ∆ε11 + ∆σ11 − ν∆εe11
2 2E
1 ν 1 (81)
= ∆σ11 − ∆ε11 − ∆ε11
2E E 2
1 − 2ν 1
= ∆σ11 − ∆ε11
2E 2
16
5. Numerical Return Mapping
Wilkins [23] was one of the first to introduce a return mapping procedure, and a variety of schemes have
been proposed following this work (see e.g. [2, 5, 4, 3, 24]). In general, it is possible to categorize the
stress integration methods in two different 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 flow n = ∂f /∂σ. The forward Euler method uses the known stress state at the
previous or trial configuration to determine the plastic flow 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.
tr ial ial
σ n+1 σtr
n+1
σ1n+1 σ1n+1
2
σn+1 σ2n+1
σn+1 σn+1
σn σn
fn+1 = 0 fn+1 = 0
(a) Explicit stress update using forward Euler (b) Fully implicit stress update using backward Euler
(cutting plane method) (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 fit 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 first 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 flow 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 justified 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
∆σ = C : (∆ε − ∆εp ) = ∆σ e + ∆σ p (82)

where the stress update at time tn+1 is expressed as
trial
σn+1 = σn + ∆σn+1 = σn + C : ∆εn+1 − ∆λn+1 C : nn+1 = σn+1 − ∆σ p (83)
in which the viscoplastic corrector ∆σ p returns the stress to the yield surface.
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 first
order differential equations then reads
εn+1 = εn + ∆εn+1
0
3 σ
εpn+1 = εpn + ∆εpn+1 = εpn
+ ∆λn+1 nn+1 = εpn + ∆λn+1 n+1
2 σeq,n+1
" 1 #
ϕn+1 C
λn+1 = λn + ∆λn+1 = λn + ṗ0 ∆tn+1 −1
σy,n+1
σn+1 = σn + ∆σn+1 = σn + C : (∆εn+1 − ∆εpn+1 ) = σn+1
trial
− ∆λn+1 C : nn+1
2
(84)
X
Rn+1 = Rn + ∆Rn+1 = Rn + ∆λn+1 HR,n+1 = Rn + ∆λn+1 Qk Ck exp(−Ck λn+1 )
k=1
σy,n+1 = σ0 + Rn+1
χσeq,n+1 ∆λn+1
Tn+1 = Tn + ∆Tn+1 = Tn +
ρCp
Dn+1 = Dn + ∆Dn+1 = Dn +hσ1,n+1 i ∆λn+1
C m
Tn+1 − Tr

∆λn+1
fv,n+1 = ϕn+1 − σy,n+1 1 + 1−
ṗ0 ∆tn+1 Tm − Tr
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
trial
gradient nn+1 iteratively, and the viscoplastic corrector is applied from the elastic predictor σ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 configuration n to determine the new converged stress state at step n + 1 (Figure 7a).
By using the normal nin+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.,
i i
i+1 i ∂f ∂f
fn+1 = fn+1 + ∆σn+1 + ∆Xn+1 (85)
∂σ n+1 ∂X n+1
which may be written in terms of the internal variables X as
i+1 i ∂f ∂f ∂f
fn+1 = fn+1 + δσn+1 + δσy,n+1 + δυp,n+1 = 0 (86)
∂σn+1 ∂σy,n+1 ∂υp,n+1
where the flow stress σy,n+1 and the following expressions are chosen as the internal variables X in the
Taylor-expansion
C m
Tn − Tr

∆λn+1
υp (∆λn+1 ) = 1+ Γp (Tn ) = 1− (87)
Introducing the partial derivatives of the internal variables in Eq. (86) this now reads
i+1 ∂f ∂f ∂υp,n+1 i+1

fn+1 i
= fn+1 − C δλi+1 i+1
n+1 − υp,n+1 Γp,n HR,n+1 δλn+1 − σy,n+1 Γp,n δ λ̇n+1 = 0 (88)
∂σn+1 ∂σn+1 ∂ λ̇n+1
Moreover, using ṗ = λ̇ = ∆λ
∆t the internal variable controlling the viscoplasticity υp (ṗ) may be expressed as
υp (∆λ) and Eq. (88) may be written as
i+1 ∂f ∂f ∂υp,n+1 i+1

fn+1 i
= fn+1 − C δλi+1 − υp,n+1 Γp,n HR,n+1 δλi+1
n+1 − σy,n+1 Γp,n δλ =0 (89)
∂σn+1 ∂σn+1 n+1 ∂∆λn+1 n+1
Finally, the incremental change δλ in the incremental plastic multiplier ∆λ is then expressed as
i
fn+1
δλi+1
n+1 = ∂f ∂f p,n+1∂υ
(90)
∂σn+1 C ∂σn+1 + υp,n+1 Γp,n HR,n+1 + σy,n+1 Γp,n ∂∆λ n+1
where
C−1
∂υp (∆λn+1 ) C ∆λn+1
= 1+ (91)
∂∆λ ṗ0 ∆tn+1 ṗ0 ∆tn+1
2
dRn+1 X
HR,n+1 = = Qk Ck exp(−Ck λn+1 ) (92)
dpn+1
k=1
T " 2 2 2 #
∂f ∂f ∂f ∂f ∂f ∂f ∂f
[C] = Cijkl = C1111 + +
∂σ ∂σ ∂σij ∂σkl ∂σ11 ∂σ22 ∂σ33

∂f ∂f ∂f ∂f ∂f ∂f
+ 2C1122 + + (93)
∂σ11 ∂σ22 ∂σ22 ∂σ33 ∂σ33 ∂σ11
" 2 2 2 #
∂f ∂f ∂f
+ 4C1212 + +
∂σ12 ∂σ23 ∂σ31
19
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 effect 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 simplifies the calculations leading to Eq. (90) without any significant
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 nin+1 instead of ni+1
n+1 violates the associated flow 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 effect
on the stress state.
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)
 
C1 C4 C4 0 0 0
C4 C1 C4 0 0 0
 
C4 C4 C1 0 0 0
C = λel I ⊗ I + 2µel I = 
0

 0 0 C7 0 0 
0 0 0 0 C7 0
0 0 0 0 0 C7
(1 − ν)E νC1 1 E
C1 = , C4 = , C7 =
(1 + ν)(1 − 2ν) 1−ν 21+ν
ELSEIF a 2D stress state is used THEN use Eq. (60)

 
C1 C4 0
E 1 E
C = C4 C1 0  , C1 = , C4 = νC1 , C7 =
1 − ν2 21+ν
0 0 C7
ELSE a 1D stress state is used THEN use
 
C1 0 0
C=0 0 0 , C1 = E
0 0 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)
r
trial trial 3 0 trial 0 trial
σeq,n+1 = ϕ(σn+1 ) = σ : σn+1
2 n+1
20
5. Check if temperature softening is included in the material input (i.e., m 6= 0).
6. Set the incremental plastic multiplier ∆λn+1 equal to zero.
trial
7. Check for plastic admissibility, i.e., IF f (σ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 . 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
configuration.
0
3 σn+1
nn+1 =
2 σeq
(ii) Compute hardening modulus HR,n+1 according to Eq. (42)
i
HR,n+1 = Q1 C1 exp(−C1 pin+1 ) + Q2 C2 exp(−C2 pin+1 )
(iii) Compute denominator in Eq. (90) using Eqs. (91), (92), (93) and (98)
i
∂f ∂f i i i
∂υp,n+1
DENOM = i
C i
+ υ Γ H
p,n+1 p,n R,n+1 + σ Γ
y,n+1 p,n
∂σn+1 ∂σn+1 ∂∆λin+1
(iv) Compute incremental change in the plastic multiplier δλi+1

n+1 in Eq. (90) and update variable
i+1
∆λn+1
i
fn+1
δλi+1
n+1 = ∂f ∂f p,n+1 ∂υ
∆λi+1 i i+1
n+1 = ∆λn+1 + δλn+1
(v) Update internal variables dependent on ∆λi+1

n+1
pi+1 i i+1
n+1 = pn+1 + δλn+1
i+1
Rn+1 i
= Rn + HR,n+1 ∆λi+1
n+1
i+1 i+1
σy,n+1 = A + HR,n+1 ∆λi+1
n+1
(vi) Update stress components using the elastic Hooke’s law according to Eq. (83)
i+1
σn+1 trial
= σn+1 − ∆λi+1 i+1
n+1 C : nn+1
Note that the material routine uses a a fully vectorized representation of the stress and strain
vectors, i.e.,
σ = [σ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)
r
i+1 i+1 3 0 ,i+1 0 ,i+1
σeq,n+1 = ϕ(σn+1 ) = σ : σn+1
2 n+1
21
(viii) Update the yield function in Eq. (84) according to the new stress state computed in (vi) and (vii)
!C m
∆λi+1 Tn − Tr

i+1 i+1 i+1 n+1
fv,n+1 = ϕn+1 − σy,n+1 1 + 1−
(viiii) Check for convergence, i.e., IF

f (∆λi+1 )
v n+1
< TOL
σy,n+1 υp Γp

THEN continue to 7.
ELSE i = i + 1 and GO TO (i)
8. Update the temperature according to the low coupling approach using Eq. (45) and check IF T ≥ Tm
χσeq,n+1 ∆λn+1
Tn+1 = Tn +
ρCT
9. Compute the hydrostatic and deviatoric stress respectively according to Eq. (32)
1
σH,n+1 = tr(σn+1 )
3
0
σn+1 = σn+1 − σH,n+1 I
0
10. Computing the principal stresses σi,n+1 based on the principal deviatoric stresses σi,n+1 as shown in
Eqs. (55) and (57)
r
0 J2,n+1 θ
σ1,n+1 = 2 cos
3 3
r
0 J2,n+1 θ 2π 0
σ2,n+1 = 2 cos − ≤ σ1,n+1
3 3 3
r
0 J2,n+1 θ 2π 0
σ3,n+1 = 2 cos + ≤ σ2,n+1
3 3 3
where
J3,n+1
cos θ = q , 0≤θ≤π
3
2 J2,n+1 /27
0 0 1
σi,n+1 = σi,n+1 + σH,n+1 = σi,n+1 + tr(σn+1 )
3
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
Wn+1
Dn+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 stiffness 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
1 − 2ν
∆ε33,n+1 = (∆σ11,n+1 + ∆σ22,n+1 ) − (∆ε11,n+1 + ∆ε22,n+1 )
E
ε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 find the minimum distance from
the trial state to the final state at step n + 1 (see Figure 7b). In the particular case of a 3D stress state, von
Mises yield criterion with associated flow 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 final stress state nn+1 (Figure 8). The return to
the yield surface from the elastic trial state is radial in the deviatoric plane, which simplifies 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)
0 0
trial E
σeq = ϕ(σn+1 ) = ϕ(σn+1 ) − 3G∆λn+1 , G= (94)
2(1 + ν)
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
C m
Tn+1 − Tr

∆λn+1
fv,n+1 (σ) = ϕtrial
n+1 (σ) − 3G∆λ n+1 − σy,n+1 1 + 1 − =0 (95)
ṗ0 ∆t Tm − Tr
Then, the incremental plastic multiplier ∆λn+1 can be solved by means of a local Newton-Raphson iteration,
i.e.,
fv (∆λin+1 )
∆λi+1 i
n+1 = ∆λn+1 − ∂fv (∆λin+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 effect 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
σ3 0
trial
σn+1

∂f
π-plane
+1 al
∂σ
n+1
) σn 0tri
0
σn+1
a 1
( σ 0t λn+
1 l
0
σn
n+ri
∆
3G
ϕ
surface at tn
σ1 σ2
surface at tn+1
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 ) ∂ h 0 trial i

= ϕ σn+1 − 3G∆λn+1 − σy,n+1 υp (∆λn+1 )Γp (Tn )
∂∆λ ∂∆λ (97)
∂σy,n+1 ∂υp (∆λn+1 )
= −3G − υp (∆λn+1 )Γp (Tn ) − σy,n+1 Γp (Tn )
∂∆λ ∂∆λ
where
C
∆λn+1
υp (∆λn+1 ) = 1+
ṗ0 ∆t
C−1
∂υp (∆λn+1 ) C ∆λn+1 (98)
= 1+
∂∆λ ṗ0 ∆t ṗ0 ∆t
m
Tn − Tr

Γp (Tn ) = 1 −
Tm − Tr
The numerical scheme within the time increment ∆tn+1 may then be summarized as follows:
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
trial
σn+1 = σn + C : ∆εn+1
2. Compute the equivalent stress in terms of the elastic trial stress

r
trial 0
trial 3 0 trial 0 trial
σeq = ϕ(σn+1 ) = σ : σn+1
2 n+1
and the radial normal vector for the yield surface
0
trial
3 σn+1
nn+1 = trial
2 σeq
24
0
trial
3. Check for plastic admissibility, i.e., IF f (σn+1 ) > 0 THEN perform local Newton-Raphson iteration
to find ∆λn+1
fv (∆λin+1 )
(i) ∆λi+1 i
n+1 = ∆λn+1 − ∂fv (∆λi )
n+1
∂∆λ

fv (∆λi+1 )
(ii) IF σy,n+1 n+1
υp Γp > 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.,
trial
σn+1 = σn+1 − ∆λn+1 C : nn+1
6. Check for element erosion, i.e., IF D > DC THEN the material point stiffness 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. Verification of Material Subroutine
The performance of the material subroutine is verified and validated in EPX. First, a verification 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 specified in the input file [7]. The material
constants used for the verification 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 verification

An extensive single element verification of the VPJC material model is performed in Reference [28]. The
following single element verification is therefore mainly included for completeness of this study, and it is
referred to [28] for more detailed single element verifications. 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 findings in [28], which concluded that the
25
Table 1: Material parameters used in numerical simulations [7, 21].
Elastic constants Yield stress Strain-rate Damage

and density and strain hardening hardening
E ν ρ A Q1 C1 Q2 C2 ṗ0 C WC DC
[MPa] [-] kg/m3 [MPa] [MPa] [-] [MPa] [-] [s−1 ] [-] [MPa] [-]
210000 0.33 7850 370 236.4 39.3 408.1 4.5 5e-4 0.01 473 1.0
Adiabatic heating and temperature softening Convergence Solution algorithm
Cp χ Tm Tr m TOL MXITER SOLU
[J/kgK] [-] [K] [K] [-] [-] [-] [1: cutting plane or 2: radial return]
452 0.9 1800 293 1.0 10−5 100 1
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 finite 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 sufficient 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 finite rotation. The verification of the material routine is therefore limited to the uniaxial tension
and simple shear tests in Figure 9.
v3 = L̇
∆L
u2 v2 = u̇2
L0 θ
L0
(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
The numerical models consist of a simple input file 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 6= 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.,

e p L ∆L
ε33 = ε33 + ε33 = ln = ln 1 + (99)
L0 L0
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 εe33 of the longitudinal strain, Eq. (99) may be simplified
to

L(t)
ε33 ≈ εp33 = p = ṗt = ln (100)
L0
which may be expressed as
L(t) = L0 exp (ṗt) (101)

where the displacement is given by
d33 (t) = L(t) − L0 = L0 (exp (ṗt) − 1) (102)

with corresponding velocity
v3 (t) = L0 ṗ exp (ṗt) , v3 (t = 0) = v3,0 = L0 ṗ (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 final 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̇
defined in Eq. (37), which for the uniaxial test reads
r r h s
2 p 2 p 2 p 2 p 2
i 2 6 p 2
ṗ = D : Dp = (ε̇11 ) + (ε̇22 ) + (ε̇33 ) = (ε̇33 ) = ε̇p33 (104)
3 3 3 4
when assuming that the material is incompressible during plastic deformation (see Eq. (64)). The corre-
sponding expression in simple shear is
v "
r u 2 2 #
2 p u2 1 1 γ12
˙ γ̇
ṗ = D : Dp = t γ12
˙ + γ21
˙ = √ =√ (105)
3 3 2 2 3 3
√
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 first tested with respect to rate-independent plas-
ticity. Then, the model is tested for its precision in predicting both viscoplastic and thermoviscoplastic
27
Table 2: Loading in numerical test programme.
Loading Uniaxial tension Simple shear

Velocity [m/s] v3,0 = 1.0 0.1 0.01 0.001 v2 = 1.732 0.1732 0.01732 0.001732
Termination [s] 0.002 0.02 0.2 2 0.0021 0.021 0.21 2.1
ṗ [s−1 ] 1000 100 10 1 1000 100 10 1
behaviour, before finally 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 verified 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
threshold was also set to an exaggerated value to avoid element erosion during plastic straining, i.e., DC 1.
The analytical solution then reduces to the first 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 ).
1200
1000
800
Stress [MPa]
600
400 σeq (EPX)

R(p) (EPX)
200 σeq (analytical)
R(p) (analytical)
0
0 0.5 1 1.5 2
Plastic strain p [mm/mm]
Figure 10: Results from strain-rate and temperature independent 1D single element (FUN3 ) test in uniaxial
tension (C = 0, m = 0 and DC 1). Numerical results in EPX compared to analytical solution in Eq. (41).
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
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600
400 σeq (EPX) 400 σeq (EPX)

R(p) (EPX) R(p) (EPX)
200 σeq (analytical) 200 σeq (analytical)
R(p) (analytical) R(p) (analytical)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Plastic strain p [mm/mm] Plastic strain p [mm/mm]
(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 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.
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600

0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
(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 1). Numerical results in EPX compared to analytical solution in Eq. (41).
Finally, the material model is verified 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
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600

0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
(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 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 confirmed 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).
σeq ∆σeq R(p = 2.0) ∆R(p = 2.0)

Element
[MPa] [%] [MPa] [%]
Analytical 1014.5 - 644.5 -
FUN3 (uniaxial) 1015.1 0.06 645.1 0.10
Q42L (uniaxial) 1014.1 -0.04 644.1 -0.06
Q42L (shear) 1013.9 -0.06 643.9 -0.09
Q4GS (uniaxial) 1014.1 -0.04 644.1 -0.06
Q4GS (shear) 1013.9 -0.06 643.9 -0.09
CUBE (uniaxial) 1014.9 0.04 644.9 0.06
CUBE (shear) 1014.0 -0.04 644.0 -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
erosion are excluded from the numerical tests by choosing m = 0 and DC 1. The analytical solution then
reduces to the two first terms in Eq. (41), against which the numerical results will be compared. Finally,
the rate-independent response (ṗ ≤ 5e−4 ) from Section 6.1.1 is also plotted for comparison.
30
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 ).
1200
1000
800
Stress [MPa]
600
σeq (ṗ ≤ 5e−4 , EPX)

400 σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX)
σeq (analytical)
0
0 0.5 1 1.5 2
Figure 14: Results from temperature independent 1D single element (FUN3 ) test in uniaxial tension (m = 0
and DC 1). Numerical results in EPX compared to the first two terms of the analytical solution in
Eq. (41).
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).
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600
σeq (ṗ ≤ 5e−4 , EPX) σeq (ṗ ≤ 5e−4 , EPX)

400 σeq (ṗ = 1, EPX) 400 σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
σeq (analytical) σeq (analytical)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 15: Results from temperature independent 2D single element (Q42L) tests (m = 0 and DC 1).
Numerical results in EPX compared to the first two terms of the analytical solution in Eq. (41).
3D shell element
The same tests are performed for the Q4GS element in uniaxial tension (Figure 16a) and simple shear
31
(Figure 16b), where the numerical results are compared to the first two terms of the analytical solution in
Eq. (41).
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600

400 σeq (ṗ = 1, EPX) 400 σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 16: Results from temperature independent 3D single (shell) element (Q4GS ) tests (m = 0 and
DC 1). Numerical results in EPX compared to the first two terms of the analytical solution in Eq. (41).
Finally, the viscoplastic response of the material model is verified for the CUBE element in uniaxial tension
(Figure 17a) and simple shear (Figure 17b).
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600

400 σeq (ṗ = 1, EPX) 400 σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 17: Results from temperature independent 3D single (continuum) element (CUBE ) tests (m = 0 and
DC 1). Numerical results in EPX compared to the first two terms of the analytical solution in Eq. (41).
It is seen that the equivalent stress σeq from the numerical tests is in excellent agreement with the first
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 confirmed by
comparing the equivalent stress σeq and the hardening variable R(p) at the end of the simulations (p = 2.0)
32
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).
Plastic strain rate σeq ∆σeq R(p = 2.0) ∆R(p = 2.0)

Element
[s−1 ] [MPa] [%] [MPa] [%]
Analytical 1014.5 - 644.5 -
FUN3 (uniaxial) 1015.1 0.06 645.1 0.10
Q42L (uniaxial) 1014.1 -0.04 644.1 -0.06
Q42L (shear) 1013.9 -0.06 643.9 -0.09
≤ 5e−4
Q4GS (uniaxial) 1014.1 -0.04 644.1 -0.06
Q4GS (shear) 1013.9 -0.06 643.9 -0.09
CUBE (uniaxial) 1014.9 0.04 644.9 0.06
CUBE (shear) 1014.0 -0.04 644.0 -0.07
Analytical 1094.6 - 644.4 -
FUN3 (uniaxial) 1094.6 0.0 644.4 0.0
Q42L (uniaxial) 1094.6 0.0 644.4 0.0
Q42L (shear) 1094.6 0.0 644.5 0.0
1
Q4GS (uniaxial) 1094.6 0.0 644.4 0.0
Q4GS (shear) 1094.6 0.0 644.5 0.0
CUBE (uniaxial) 1094.6 0.0 644.5 0.0
CUBE (shear) 1094.6 0.0 644.5 0.0
Analytical 1120.1 - 644.4 -
FUN3 (uniaxial) 1120.1 0.0 644.4 0.0
Q42L (uniaxial) 1120.1 0.0 644.4 0.0
Q42L (shear) 1120.1 0.0 644.5 0.0
10
Q4GS (uniaxial) 1120.1 0.0 644.5 0.0
Q4GS (shear) 1120.1 0.0 644.5 0.0
CUBE (uniaxial) 1120.1 0.0 644.5 0.0
CUBE (shear) 1120.1 0.0 644.5 0.0
Analytical 1146.1 - 644.4 -
FUN3 (uniaxial) 1146.1 -0.01 644.4 -0.02
Q42L (uniaxial) 1146.1 -0.01 644.4 -0.01
Q42L (shear) 1146.1 -0.01 644.4 -0.01
100
Q4GS (uniaxial) 1146.1 -0.01 644.4 -0.01
Q4GS (shear) 1146.1 -0.01 644.4 -0.01
CUBE (uniaxial) 1146.1 -0.01 644.4 -0.01
CUBE (shear) 1146.1 0.0 644.4 -0.01
Analytical 1172.8 - 644.4 -
FUN3 (uniaxial) 1171.9 -0.08 643.6 -0.13
Q42L (uniaxial) 1172.0 -0.07 643.7 -0.11
Q42L (shear) 1172.2 -0.06 643.9 -0.09
1000
Q4GS (uniaxial) 1172.0 -0.07 643.7 -0.11
Q4GS (shear) 1172.2 -0.06 643.9 -0.09
CUBE (uniaxial) 1172.2 -0.05 643.9 -0.08
CUBE (shear) 1172.4 -0.04 644.0 -0.07
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 ).
1200 800
700
1000
600
Temperature [K]
800
Stress [MPa]
500
600 400
σeq (ṗ ≤ 5e−4 , EPX) 300 T (p) (ṗ ≤ 5e−4 , EPX)

400 (ṗ = 1, EPX)
σeq T (p) (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) 200 T (p) (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) T (p) (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) 100 T (p) (ṗ = 1000, EPX)
σeq (analytical) T (p) (analytical)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
(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 1). Numerical results in EPX compared to the analytical solution in Eq. (41) and Eq. (45).
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
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600

400 σeq (ṗ = 1, EPX) 400 σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
(a) Equivalent stress σeq (uniaxial tension) (b) Equivalent stress σeq (simple shear)
800 800
700 700
600 600
Temperature [K]
Temperature [K]
500 500
400 400
300 T (p) (ṗ ≤ 5e−4 , EPX) 300 T (p) (ṗ ≤ 5e−4 , EPX)
T (p) (ṗ = 1, EPX) T (p) (ṗ = 1, EPX)
200 T (p) (ṗ = 10, EPX) 200 T (p) (ṗ = 10, EPX)
T (p) (ṗ = 100, EPX) T (p) (ṗ = 100, EPX)
100 T (p) (ṗ = 1000, EPX) 100 T (p) (ṗ = 1000, EPX)
T (p) (analytical) T (p) (analytical)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
(c) Temperature T (uniaxial tension) (d) Temperature T (simple shear)
Figure 19: Results from temperature dependent (m = 1.0 and DC 1) 2D single (continuum) element
(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
(Figure 20b).
36
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600

400 σeq (ṗ = 1, EPX) 400 σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
800 800
700 700
600 600
Temperature [K]
Temperature [K]
500 500
400 400
300 T (p) (ṗ ≤ 5e−4 , EPX) 300 T (p) (ṗ ≤ 5e−4 , EPX)
T (p) (ṗ = 1, EPX) T (p) (ṗ = 1, EPX)
200 T (p) (ṗ = 10, EPX) 200 T (p) (ṗ = 10, EPX)
T (p) (ṗ = 100, EPX) T (p) (ṗ = 100, EPX)
100 T (p) (ṗ = 1000, EPX) 100 T (p) (ṗ = 1000, EPX)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 20: Results from temperature dependent (m = 1.0 and DC 1) 3D single (shell) element (Q4GS )
tests in uniaxial tension (left) and simple shear (right). Numerical results in EPX compared to the analytical
solution in Eq. (41) and Eq. (45).
Finally, the thermoviscoplastic response of the material model is verified for the CUBE element in uniaxial
tension (Figure 21a) and simple shear (Figure 21b).
37
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600

400 σeq (ṗ = 1, EPX) 400 σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
800 800
700 700
600 600
Temperature [K]
Temperature [K]
500 500
400 400
300 T (p) (ṗ ≤ 5e−4 , EPX) 300 T (p) (ṗ ≤ 5e−4 , EPX)
T (p) (ṗ = 1, EPX) T (p) (ṗ = 1, EPX)
200 T (p) (ṗ = 10, EPX) 200 T (p) (ṗ = 10, EPX)
T (p) (ṗ = 100, EPX) T (p) (ṗ = 100, EPX)
100 T (p) (ṗ = 1000, EPX) 100 T (p) (ṗ = 1000, EPX)
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
Figure 21: Results from temperature dependent (m = 1.0 and DC 1) 3D single (continuum) element
(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
Table 5: Summary of single element tests for thermo-viscoplastic behaviour behaviour at the end of the
simulations (p = 2.0).
Plastic strain rate σeq ∆σeq T (p = 2.0) ∆T (p = 2.0)

Element
[s−1 ] [MPa] [%] [K] [%]
Analytical 733.0 - 711.3 -
FUN3 (uniaxial) 733.7 0.10 711.2 -0.02
Q42L (uniaxial) 732.9 -0.07 710.5 -0.01
Q42L (shear) 732.6 -0.05 711.1 -0.01
≤ 5e−4
Q4GS (uniaxial) 732.9 -0.07 710.5 -0.01
Q4GS (shear) 732.6 -0.05 711.1 -0.01
CUBE (uniaxial) 733.5 0.06 710.9 -0.01
CUBE (shear) 732.7 -0.04 711.1 -0.01
Analytical 771.1 - 738.4 -
FUN3 (uniaxial) 771.1 0.0 738.3 -0.02
Q42L (uniaxial) 771.1 0.0 738.3 -0.02
Q42L (shear) 770.7 0.0 738.9 -0.02
1
Q4GS (uniaxial) 771.1 0.0 738.2 -0.02
Q4GS (shear) 770.7 0.0 738.9 -0.02
CUBE (uniaxial) 771.1 0.0 738.3 -0.02
CUBE (shear) 770.7 0.0 738.9 -0.02
Analytical 782.7 - 747.0 -
FUN3 (uniaxial) 782.7 0.0 746.9 -0.01
Q42L (uniaxial) 782.7 0.0 746.9 -0.01
Q42L (shear) 782.3 0.0 747.4 -0.02
10
Q4GS (uniaxial) 782.7 0.0 746.9 -0.01
Q4GS (shear) 782.3 0.0 747.4 -0.02
CUBE (uniaxial) 782.7 0.0 746.9 -0.01
CUBE (shear) 782.3 0.0 747.5 -0.02
Analytical 794.3 - 755.8 -
FUN3 (uniaxial) 794.2 -0.01 755.6 -0.01
Q42L (uniaxial) 794.2 -0.01 755.6 -0.01
Q42L (shear) 793.8 -0.01 756.3 -0.02
100
Q4GS (uniaxial) 794.3 0.0 755.6 -0.01
Q4GS (shear) 793.8 -0.01 756.3 -0.02
CUBE (uniaxial) 794.3 0.0 755.6 -0.01
CUBE (shear) 793.8 0.0 756.3 -0.02
Analytical 806.0 - 764.5 -
FUN3 (uniaxial) 805.5 -0.06 764.3 -0.02
Q42L (uniaxial) 805.5 -0.07 764.3 -0.01
Q42L (shear) 805.2 -0.05 764.9 -0.01
1000
Q4GS (uniaxial) 805.6 -0.07 764.2 -0.01
Q4GS (shear) 805.2 -0.05 764.9 -0.01
CUBE (uniaxial) 805.6 -0.05 764.3 -0.01
CUBE (shear) 805.2 -0.04 764.9 -0.01
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 file under the directive defining 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 stiffness 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 sufficiently 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 ).
1200
1000
800
Stress [MPa]
600
400 σeq (ṗ ≤ 5e−4 , EPX)

σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX)
Failure
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
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).
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
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600
400 σeq (ṗ ≤ 5e−4 , EPX) 400 σeq (ṗ ≤ 5e−4 , EPX)
σeq (ṗ = 1, EPX) σeq (ṗ = 1, EPX)
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
Failure Failure
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
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
(Figure 24b).
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
Failure Failure
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
Figure 24: Results from temperature and strain-rate dependent 3D single (shell) element (Q4GS ) tests
including failure (m = 1.0 and DC = 1).
Finally, the models capability in predicting ductile failure is verified for the CUBE element in uniaxial
tension (Figure 25a) and simple shear (Figure 25b).
41
1200 1200
1000 1000
800 800
Stress [MPa]
Stress [MPa]
600 600
σeq (ṗ = 10, EPX) σeq (ṗ = 10, EPX)
200 σeq (ṗ = 100, EPX) 200 σeq (ṗ = 100, EPX)
σeq (ṗ = 1000, EPX) σeq (ṗ = 1000, EPX)
Failure Failure
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
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 flow stress σy (σ0 , p, ṗ). Consequently, the equivalent stress will increase to fulfill
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 difference 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 differences observed in failure strain in Figures 23 - 25
corresponds well with the corresponding difference predicted by Eq. (51) in Figure 6.
6.2. Tension test

The material parameters (Q1 , C1 , Q2 and C2 ) in the first 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 diffuse 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 significant 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
[29] for more information about the experimental work and numerical model, respectively, for these tension
tests.
R
50.05
15
12.5
φ1
40
2.7
70 25.05
200
(a) Geometry (experimental)
x vx = 2.1mm/min
70
Nodal constraints
Nodal loads
(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.
800
700
600
Nominal stress [MPa]
500
400 Experiment 0°
Experiment 45°
300 Experiment 90°
FEA (Q42L)
200 FEA (Q4GS)
Wc = 473 MPa
100
Wc = 473 MPa
0
0 0.05 0.1 0.15 0.2 0.25 0.3
Nominal strain [mm/mm]
Figure 27: Nominal stress-strain curves from uniaxial tension tests along three different 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
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.
The experimental setup is shown in Figure 28, and consists of a steel mounting frame fixed to the concrete
floor with outer dimensions 1.0 × 1.0 × 0.015 m3 and a square opening of 0.3 × 0.3 m2 at the centre. The
square plate specimen with dimensions of 0.4 × 0.4 × 0.0008 m3 were clamped to the rigid mounting frame
using M12 bolts and a clamping frame in an attempt to achieve fixed boundary conditions. The explosive
mass W was positioned at stand-off 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 field 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.
2
Cam
Test specimen
Charge of 30 g C4 Blast pencil

∼ 25 ◦
Clamping frame
Cam
Steel mounting frame 1
Stand-off-distance ∼ 3300mm
(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
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 coefficient between the plate and clamping frames was set to 0.15, while the dynamic friction
coefficient 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
effect of the pre-tensioning of the bolts [29].
A comparison of the experimental and numerical results in terms of mid-point deflection 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-off distances and the reversed
snap buckling is captured at the largest stand-off distance. The representation of the maximum displacement
is also quite good at the closest stand-off 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
35
30
25
Mid−point deflection [mm]

20
15
Experiments
10 FEA (R = 0.125 m)
5 FEA (R = 0.250 m)
FEA (R = 0.375 m)
0
−5
−10
−15
0 2 4 6 8 10
Time [ms]
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 first verified 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 verification. 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)

HR = Q1 C1 exp(−C1 p) + Q2 C2 exp(−C2 p)
ELSEIF F LAG = 1 THEN use Power law hardening from Eq. (40)
n−1
HR = Bnp
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
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

P F = (D1 + D2 ∗ EXP (−D3 ∗ SIGM AH/P HI))
. ∗ ((1 + DLAM BDA/(P DOT ∗ DT )) ∗ ∗D4 )
. ∗ (1 + D5 ∗ ((T − T R)/(T M − T R)))
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 ...
IF (D >= DC) THEN
WRITE(59,*) ’Element deleted due to CL-criterion’
IF (DJC >= DC) THEN
WRITE(59,*) ’Element deleted due to JC-criterion’
Note that the material input must be extended to include five 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.54, D4 = −0.015 and D5 = 0.
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 five different 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 influence 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, Influence 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 flow, 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
i i
i+1 i ∂f ∂f
fn+1 = fn+1 + ∆σn+1 + ∆Xn+1 (A.1)
∂σ n+1 ∂X n+1
which may be written in terms of the internal variables X as
i+1 i ∂f ∂f ∂f
fn+1 = fn+1 + δσn+1 + δσy,n+1 + δυp,n+1 = 0 (A.2)
∂σn+1 ∂σy,n+1 ∂υp,n+1
where the flow stress σy,n+1 and the following expressions are chosen as the internal variables X in the
Taylor-expansion
C m
Tn − Tr

∆λn+1
υp (∆λn+1 ) = 1+ Γp (Tn ) = 1− (A.3)
Introducing the partial derivatives of the internal variables Eq. (A.2) reads
i+1 i ∂f ∂f ∂f
fn+1 = fn+1 + δσn+1 + δσy,n+1 + δυp,n+1 =0 (A.4)
∂σn+1 | {z } ∂σy,n+1 | {z } ∂υp,n+1 | {z }
−C ∂f
δλi+1 ∂σy,n+1
δλi+1
∂υp,n+1
δ ṗi+1
| {z } | {z }
∂σn+1 n+1 n+1 n+1
(−1)υp,n+1 Γp,n ∂λn+1 (−1)σy,n+1 Γp,n ∂ ṗn+1
i+1 ∂f ∂f ∂σy,n+1 i+1 ∂υp,n+1 i+1

fn+1 i
= fn+1 − C δλi+1
n+1 − υp,n+1 Γp,n δλn+1 − σy,n+1 Γp,n δ λ̇n+1 = 0 (A.5)
∂σn+1 ∂σn+1 ∂λn+1 ∂ λ̇n+1
i+1 ∂f ∂f ∂υp,n+1 i+1

fn+1 i
= fn+1 − C δλi+1 i+1
n+1 − υp,n+1 Γp,n HR,n+1 δλn+1 − σy,n+1 Γp,n δ λ̇n+1 = 0 (A.6)
∂σn+1 ∂σn+1 ∂ λ̇n+1
Moreover, using ṗ = λ̇ = ∆λ
∆t the internal variable controlling the viscoplasticity υp (ṗ) may be expressed as
υp (∆λ) and Eq. (A.6) may be written as
i+1 ∂f ∂f ∂υp,n+1 i+1

fn+1 i
= fn+1 − C δλi+1 − υp,n+1 Γp,n HR,n+1 δλi+1
n+1 − σy,n+1 Γp,n δλ =0 (A.7)
∂σn+1 ∂σn+1 n+1 ∂∆λn+1 n+1
Finally, the incremental change δλ in the incremental plastic multiplier ∆λ is then expressed as
i
fn+1
δλi+1
n+1 = ∂f ∂f p,n+1 ∂υ
(A.8)
where
C−1
∂υp (∆λn+1 ) C ∆λn+1
= 1+ (A.9)
∂∆λ ṗ0 ∆t ṗ0 ∆t
2
dRn+1 X
HR,n+1 = = Qk Ck e(−Ck λn+1 ) (A.10)
dpn+1
k=1
T
∂f ∂f
[C] (A.11)
∂σ ∂σ
49
 
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
[C] =  C2122 C2133 C2112 C2121 C2123 C2132 C2113 C2131 
 (A.12)
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
 ∂f 
∂σ11
 
 ∂f  f,11
 ∂σ22  f,22 
 ∂f   
 ∂σ33  f,33 
 ∂f   
 ∂σ12  f,12 
∂f  ∂f   
 = f,21 
= ∂σ (A.13)
∂σ  21  
 ∂f  f,23 

 ∂σ23   
 ∂f  f,32 
 ∂σ32   
 ∂f  f,13 
 ∂σ 
∂f
13 f,31
∂σ31
T
∂f

[C] =
∂σ
f
,11 C1111 + f,22 C2211 + f,33 C3311 + f,12 C1211 + f,21 C2111 + f,23 C2311 + f,32 C3211 + f,13 C1311 + f,31 C3111 
f,11 C1122 + f,22 C2222 + f,33 C3322 + f,12 C1222 + f,21 C2122 + f,23 C2322 + f,32 C3222 + f,13 C1322 + f,31 C3122 

f,11 C1112 + f,22 C2212 + f,33 C3312 + f,12 C1212 + f,21 C2112 + f,23 C2312 + f,32 C3212 + f,13 C1312

+ f,31 C3112  (A.14)
 
 
 
f,11 C1113 + f,22 C2213 + f,33 C3313 + f,12 C1213 + f,21 C2113 + f,23 C2313 + f,32 C3213 + f,13 C1313 + f,31 C3113
 
f,11 C1131 + f,22 C2231 + f,33 C3331 + f,12 C1231 + f,21 C2131 + f,23 C2331 + f,32 C3231 + f,13 C1331 + f,31 C3131
T
∂f ∂f

[C] =
∂σ ∂σ
(f,11 C1111 + f,22 C2211 + f,33 C3311 + f,12 C1211 + f,21 C2111 + f,23 C2311 + f,32 C3211 + f,13 C1311 + f,31 C3111 )f,11
+ (f,11 C1122 + f,22 C2222 + f,33 C3322 + f,12 C1222 + f,21 C2122 + f,23 C2322 + f,32 C3222 + f,13 C1322 + f,31 C3122 )f,22
+ (f,11 C1112 + f,22 C2212 + f,33 C3312 + f,12 C1212 + f,21 C2112 + f,23 C2312 + f,32 C3212 + f,13 C1312 + f,31 C3112 )f,12 (A.15)
 
C1111 C1122 C1133 0 0 0 0 0 0
C2211 C2222 C2233 0 0 0 0 0 0 
 
C3311 C3322 C3333 0 0 0 0 0 0 
 
 0 0 0 C1212 C1221 0 0 0 0 
 
 0
[C] =  0 0 C2112 C2121 0 0 0 0 (A.16)


 0 0 0 0 0 C2323 C2332 0 0 
 
 0 0 0 0 0 C3223 C3232 0 0 
 
 0 0 0 0 0 0 0 C1313 C1331 
0 0 0 0 0 0 0 C3113 C3131
50
Due to the symmetry of the stress, σij = σji , it is seen that also
Cijkl = Cjikl (A.17)

Since the strain tensor also is symmetric, εkl = εlk , it is also assumed that
Cijkl = Cijlk (A.18)

and the 9x9 matrix in Eq. (A.16) may be reduced to the following 6x6 matrix
 
C1111 C1122 C1133 0 0 0
C2211 C2222 C2233 0 0 0 
 
C3311 C3322 C3333 0 0 0 
[C] = 
 0
 (A.19)
 0 0 C1212 0 0 

 0 0 0 0 C2323 0 
0 0 0 0 0 C1313
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
T
∂f ∂f
[C] =
∂σ ∂σ
(f,11 C1111 + f,22 C2211 + f,33 C3311 )f,11 + (f,11 C1122 + f,22 C2222 + f,33 C3322 )f,22
+ (f,11 C1133 + f,22 C2233 + f,33 C3333 )f,33
+ (f,12 C1212 + f,21 C2112 )f,12 + (f,12 C1221 + f,21 C2121 )f,21 (A.20)
+ (f,23 C2323 + f,32 C3223 )f,23 + (f,23 C2332 + f,32 C3232 )f,32
+ (f,13 C1313 + f,31 C3113 )f,13 + (f,13 C1331 + f,31 C3131 )f,13
2 2 2

= c1 f,11 + f,22 + f,33 + 2c4 (f,11 f,22 + f,22 f,33 + f,33 f,11 )
2 2 2

+ 4c7 f,12 + f,23 + f,31
However, since the strain tensor is symmetric, εkl = εlk when k 6= l, it is seen that the number of independent
coefficients is reduced to 3
   
C1111 C1122 C1133 0 0 0 c1 c4 c4 0 0 0
C2211 C2222 C2233 0 0 0  c4 c1 c4 0 0 0
   
C3311 C3322 C3333 0 0 0  c4 c4 c1 0 0 0
[C] = 
 0
=  (A.21)
 0 0 C1212 0 0  0
  0 0 c7 0 0 
 0 0 0 0 C2323 0  0 0 0 0 c7 0
0 0 0 0 0 C1313 0 0 0 0 0 c7
where it is used Voigt notation and symmetric stress tensors
σ4 = σ12 = σ21 = σ7
σ5 = σ23 = σ32 = σ8 (A.22)
σ6 = σ31 = σ13 = σ9
51
The symmetry of the stress and strain tensors allows us to express them as column matrices (or vectors) in
the form
   
σ11 ε11
σ22  ε22 
   
σ33  ε33 
[σ] =   and [ε] = 
   (A.23)
σ12 
γ12 
 
σ23  γ23 
σ31 γ31
where γ12 = 2ε12 , γ23 = 2ε23 and γ31 = 2ε31 are the engineering shear strains.
 
c1 c4 c4 0 0 0
c4 c1 c4 0 0 0
 
c4 c4 c1 0 0 0
[C] = 
0
 (A.24)
 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.
T
∂f ∂f ∂f ∂f
[C] = Cij
∂σ ∂σ ∂σi ∂σj
" 2 2 2 #
∂f ∂f ∂f ∂f ∂f ∂f ∂f ∂f ∂f
= c1 + + + 2c4 + + (A.25)
∂σ1 ∂σ2 ∂σ3 ∂σ1 ∂σ2 ∂σ2 ∂σ3 ∂σ3 ∂σ1
" 2 2 2 #
∂f ∂f ∂f
+ 4c7 + +
∂σ4 ∂σ5 ∂σ6
52
Appendix B. Radial Return Algorithm
As discussed in Section 5.2, the system of equations in Eq.(84) can be extensively simplified and the return
mapping for the von Mises model, with associated flow 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 efficient, and to
improve the performance of the overall finite 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 : ∆εn+1 − ∆εpn+1 = σn+1

trial
− C : ∆εpn+1 = σn+1
trial
− ∆σ̃ p

(B.1)
trial
where σn+1 is the elastic predictor moving the stress outside the yield surface, and the plastic corrector
∆σ̃ p returns the stress to the yield surface (Figure B.31).
Introducing the associated flow rule, Eq.(B.1) reads

0
3 σ
σn+1 = trial
σn+1 − ∆λn+1 C : n+1 (B.2)
2 σeq
Using the relation (Appendix C)
0 0
C : σn+1 = 2Gσn+1 (B.3)
Eq.(B.2) is given by
0
σ
σn+1 = trial
σn+1 − 3G n+1 ∆λn+1 (B.4)
σeq
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
flow vector ε̇p is purely deviatoric, and the hydrostatic stress σH,n+1 is unaffected by the plastic flow
trial
(σH,n+1 = σH,n+1 ). This implies that the stress update can be formulated more efficiently 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.
0 0
0 0 σ 0 σn+1
σn+1 = trial
σn+1 − 3G n+1 ∆λn+1 = σn+1
trial
− 3G 0 ∆λn+1 (B.5)
σeq ϕ(σn+1 )
Rearranging Eq.(B.5) results in an important result, i.e.

0 3G∆λn+1 0
trial
σn+1 1 + 0 = σn+1 (B.6)
ϕ(σn+1 )
which implies that the trial and updated deviatoric stresses are co-linear, i.e.
0 0
trial
σn+1 σn+1
0
= 0 trial
(B.7)
ϕ σn+1 ϕ σn+1
53
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 simplification

0 0
trial 3G∆λn+1
σn+1 = σn+1 1− 0 trial (B.8)
ϕ(σn+1 )
0
q
trial 0 trial

where ϕ(σn+1 )= 3J2 σn+1 is the elastic trial von Mises equivalent stress.
0 0
trial
It is also noted that σn+1 is a constant vector in the return mapping, and that the deviatoric stress σ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
0
trial
stress by a factor 1 − 3G∆λn+1 /ϕ(σn+1 ). The geometric representation of this update formula in the
deviatoric plane is illustrated in Figure B.31.
σ3 0
trial
σn+1

∂f
∂σ
n+1
+1 al
π-plane
) σn 0tri
0
σn+1 a 1
( σ 0t λn+
1 l
n+ri
∆
3G
0
σn
ϕ
surface at tn
σ1 σ2
surface at tn+1
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
0
f (∆λ) = ϕ σn+1 − σy,n+1
0 3G∆λn+1

trial
= ϕ σn+1 1− 0 trial − σy,n+1
ϕ(σn+1 ) (B.9)
0
trial
= ϕ σn+1 − 3G∆λn+1 − σy,n+1
0
trial
= ϕ σn+1 − 3G∆λn+1 − (σy,n + HR,n+1 ∆λn+1 ) = 0
54
0
where the incremental plastic multiplier ∆λ is the only unknown. In Eq.(B.9) it is utilized that ϕ(aσ ) =
0 0
aϕ(σ ), since the function ϕ(σ ) is assumed to be a positive homogeneous function of order one of the stress.
Thus,

0 0
trial 3G∆λn+1
ϕ(σn+1 ) = ϕ σn+1 1 − 0 trial
ϕ(σn+1 )
(B.10)
0
trial 3G∆λn+1
= ϕ(σn+1 ) 1 − 0 trial
ϕ(σn+1 )
By using Eq.(B.9), the incremental plastic multiplier at time tn+1 is expressed as
0
trial
ϕ σn+1 − σy,n
∆λn+1 = (B.11)
3G∆λn+1 + HR,n+1
which, in the case of nonlinear hardening, could be solved by a Newton-Raphson method
f (∆λ)
∆λi+1 i
n+1 = ∆λn+1 −
f 0 (∆λ)
0
trial
(B.12)
ϕ σn+1 − 3G∆λin+1 − σy,n+1
= ∆λin+1 −
−3G − HR
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)

0
3 σ
trial
σn+1 = σn+1 − ∆λn+1 C : n+1
2 σeq
0
trial
trial 3 σn+1
= σn+1 − ∆λn+1 C : 0 trial

2 ϕ σn+1 (B.13)
pn+1 = pn + ∆λn+1
0
trial
3 σn+1
εpn+1 = εpn + 0 trial ∆λn+1

2 ϕ σn+1
55
Appendix C. Isotropic elastic material
Assuming isotropic linear elastic behavior of the material reduces the stiffness tensor to only depend on the
Lamè constants, i.e
Cijkl = λel δij δkl + µel (δik δjl + δil δjk ) (C.1)
and the stress is then given by the linear relation
σij = Cijkl εkl = λel δij δkl εkl + µel (δik δjl εkl + δil δjk εkl )
= λel δij εkk + µel (εij + εji ) (C.2)
= λel δij εkk + 2µel (εij )
where the following relations are used
δik δjl εkl = δik εjk = εij

δil δjk εkl = δil εlj = εji (C.3)
εij = εji
By inserting the expression for the Lamè constants, Eq.(C.1) reads
νE E
Cijkl = δij δkl + (δik δjl + δil δjk )
(1 + ν)(1 − 2ν) 2(1 + ν)

E 1 ν
= (δik δjl + δil δjk ) + δij δkl (C.4)
1+ν 2 1 − 2ν

1 ν
= 2G (δik δjl + δil δjk ) + δij δkl
2 1 − 2ν
and by multiplying the isotropic elasticity tensor in Eq.(C.4) with the deviatoric part of the stress tensor,
this product reduces to

0 1 0 0
ν 0
Cijkl σkl = 2G δik δjl σkl + δil δjk σkl + δij δkl σkl
2 1 − 2ν
(C.5)
1 0 0
ν 0
= 2G σij + σji + δij σkk
2 1 − 2ν
0 0
Utilizing the symmetry of the stress tensor, i.e. σij = σji , and the definition of the deviatoric stress
0 1 0 1
σij = σij − σkk δij , σii = σii − σkk δii = 0 (C.6)
3 3
Eq.(C.5) reduces to the following expression
0 0
Cijkl σkl = 2Gσij (C.7)
which is a useful relation in the radial return map algorithm.
56
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 defined 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
    
ε11 1 −ν −ν 0 0 0 σ11
 ε22  −ν 1 −ν 0 0 0  σ22 
    
 = Sσ = 1 −ν −ν
 ε33   1 0 0 0  0 

2ε12 
  (D.1)
 0 0 0 2(1 + ν) 0 0
E  σ12 
   
2ε23   0 0 0 0 2(1 + ν) 0  σ23 
2ε31 0 0 0 0 0 2(1 + ν) σ31
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
    
ε11 1 −ν 0 0 0 σ11
 ε22  −ν 1 0 0 0  σ22 
2ε12  = Sσ = 1  0
    
0 2(1 + ν) 0 0  σ12  (D.2)
 
2ε23  E 0
 
0 0 2(1 + ν) 0  σ23 
2ε31 0 0 0 0 2(1 + ν) σ31
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 1 ν 0 0 0 ε11
σ22  ν 1 0 0 0   ε22 
σ12  = Cε = E 1−ν
    
0 0 0 0  2ε12  (D.3)
2
 
σ23  1 − ν2 
0 1−ν
 
 2ε23 
0 0 2 0
1−ν
σ31 0 0 0 0 2
2ε31
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 finding the principal stresses σi of the stress tensor σ are now
determined by using Eq. (57), i.e.
0 0 1
σi = σi + σH = σi + tr(σ) (D.4)
3
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

Define useful parameters related to the theory of SUBROUTINE VPJC (ITAU, IPLANC, PROPS, STRAININC, STRESSOLD,
> STRAIN, STATE, STRESSNEW, CSON, FAIL_FLAG)
elasticity * sound speed
T1 = R1 + POISS
T2 = R1 - R2*POISS
T 1 = 1 + ν, T 2 = 1 − 2ν, T3 = 1 − ν T3 = R1 - POISS
*
SELECT CASE (IPLANC)
Compute the speed of sound based on the CASE (0) ! NO PLANE STRESS CONDITION
CSON = SQRT (YOUNG*T3 / (RHO*T1*T2))
stress state, i.e., CASE (1) ! PLANE STRESS CONDITION (SIG(3)=0)
IF 3D stress state THEN CSON = SQRT (YOUNG / (RHO*T1*T3))
CASE (2) ! UNIAXIAL STRESS CONDITION (SIG(2)=SIG(3)=0)
s CSON = SQRT (YOUNG / RHO)
E(1 − ν) CASE DEFAULT ! CHECK THIS ONLY HERE
CSON = CALL ERRMSS (’VPJC’, ’WRONG IPLANC’)
ρ(1 + ν)(1 − 2ν) STOP ’VPJC - WRONG IPLANC’
END SELECT
ELSEIF plane stress THEN *
s
E
CSON =
ρ(1 − ν 2 )
ELSEIF uniaxial stress THEN

s
E
CSON =
ρ
2. Compute the elasticity tensor.
Goverening equations Implementation in EUROPLEXUS

IF a 3D stress state is used THEN use Eqs. (3)-(4) * compute the elasticity tensor
CASE (0) ! NO PLANE STRESS CONDITION
C
1 C4 C5 0 0 0  C1 = YOUNG*T3 / (T1*T2)
C5 C2 C4 0 0 0  C4 = POISS*C1 / T3
C

C6 C3 0 0 0 
 CASE (1) ! PLANE STRESS CONDITION (SIG(3)=0)
C = λel I ⊗ I + 2µel I =  5 C1 = YOUNG / (T1*T3)
 0 0 0 C7 0 0 

 0 0 0 0 C8 0  C4 = POISS*C1
0 0 0 0 0 C9 CASE (2) ! UNIAXIAL STRESS CONDITION (SIG(2)=SIG(3)=0)
C1 = YOUNG
(1 − ν)E νC1 C4 = R0 ! THESE VALUES HAVE TO BE CHECKED
C1 = C2 = C3 = , C4 = C5 = C6 = ,
(1 + ν)(1 − 2ν) 1−ν END SELECT
*
1 E C2 = C1
C7 = C8 = C9 =
C3 = C1
2 1+ν
C5 = C4
C6 = C4
ELSEIF a 2D stress state is used THEN use C7 = R5*YOUNG / T1
Eq. (60) C8 = C7
C9 = C7
  *
C1 C4 0
C = C5 C2 0 
0 0 C7
E 1 E
C1 = C2 = , C4 = C5 = νC1 , C7 =
1 − ν2 2 1+ν
ELSE a 1D stress state is used THEN use
 
C1 0 0
C =  0 0 0 , C1 = E
0 0 0
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 .

* compute the elastic trial stress and the trial equivalent stress
*
trial STRESS(1) = STRESSOLD(1) + C1*STRAININC(1) + C4*STRAININC(2)
σ11,n+1 = σ11,n + C1 ∆ε11,n+1 + C4 ∆ε22,n+1
> + C5*STRAININC(3)
+ C5 ∆ε33,n+1 *
SELECT CASE (IPLANC) ! TREAT NORMAL COMPONENTS FIRST
IF a 3D stress state is used THEN STRESS(2) = STRESSOLD(2) + C4*STRAININC(1) + C2*STRAININC(2)
> + C6*STRAININC(3)
trial STRESS(3) = STRESSOLD(3) + C5*STRAININC(1) + C6*STRAININC(2)
σ22,n+1 = σ22,n + C4 ∆ε11,n+1 + C2 ∆ε22,n+1 > + C3*STRAININC(3)
PHITRIAL = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)
+ C6 ∆ε33,n+1 > + STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
trial > - STRESS(2)*STRESS(3) - STRESS(3)*STRESS(1)
σ33,n+1 = σ33,n + C5 ∆ε11,n+1 + C6 ∆ε22,n+1 CASE (1) ! PLANE STRESS CONDITION (STRESS(3)=0)
STRESS(2) = STRESSOLD(2) + C4*STRAININC(1) + C2*STRAININC(2)
+ C3 ∆ε33,n+1
> + C6*STRAININC(3)
trial STRESS(3) = R0
σ12,n+1 = σ12,n + C7 ∆ε12,n+1 PHITRIAL = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)
trial > - STRESS(1)*STRESS(2)
σ23,n+1 = σ23,n + C8 ∆ε23,n+1 CASE (2) ! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)
STRESS(2) = R0
trial
σ31,n+1 = σ31,n + C9 ∆ε31,n+1 STRESS(3) = R0
PHITRIAL = STRESS(1)*STRESS(1)
3 0 trial 0 trial trial 2 trial 2 END SELECT
σn+1 : σn+1 = (σ11,n+1 ) + (σ22,n+1 )
2 *
SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
trial 2 trial trial
+ (σ33,n+1 ) − σ11,n+1 σ22,n+1 CASE (0) ! NO SHEAR STRESS IS PRESENT
CASE (1) ! ONLY TAU_XY (STRESS(4)) IS PRESENT
trial trial trial trial
+ σ22,n+1 σ33,n+1 − σ33,n+1 σ11,n+1 STRESS(4) = STRESSOLD(4) + C7*STRAININC(4) ! THE GAMMAS ARE USED!
PHITRIAL = PHITRIAL + R3*(STRESS(4)*STRESS(4))
trial 2 trial 2 CASE (3) ! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT
+ 3 (σ12,n+1 ) + (σ23,n+1 )
STRESS(4) = STRESSOLD(4) + C7*STRAININC(4) ! SHEAR STRAINS
trial 2 STRESS(5) = STRESSOLD(5) + C8*STRAININC(5) ! ARE THE GAMMAS, NOT
+ (σ31,n+1 )
STRESS(6) = STRESSOLD(6) + C9*STRAININC(6) ! THE EPSILONS
PHITRIAL = PHITRIAL + R3*(STRESS(4)*STRESS(4)
> + STRESS(5)*STRESS(5)
ELSEIF a 2D plane stress state is used THEN > + STRESS(6)*STRESS(6))
CASE DEFAULT ! WE CHECK THIS ONLY IN THE FIRST SELECT CASE ON ITAU
trial CALL ERRMSS (’VPJC’, ’WRONG ITAU’)
σ22,n+1 = σ22,n + C4 ∆ε11,n+1 + C2 ∆ε22,n+1 STOP ’VPJC - WRONG ITAU’
END SELECT
+ C6 ∆ε33,n+1
*
trial PHITRIAL = SQRT (PHITRIAL)
σ33,n+1 = 0 *
trial
σ12,n+1 = σ12,n + C7 ∆ε12,n+1
3 0 trial 0 trial trial 2 trial 2

σn+1 : σn+1 = (σ11,n+1 ) + (σ22,n+1 )
2
trial trial trial 2
− σ11,n+1 σ22,n+1 + 3 (σ12,n+1 )
ELSE a 1D stress state is used and

trial
σ22,n+1 = 0
trial
σ33,n+1 = 0
3 0 trial 0 trial trial 2

σn+1 : σn+1 = (σ11,n+1 )
2
END IF
s
trial trial 3 0 0
σeq,n+1 = ϕ(σn+1 ) = σ trial : σ trial
n+1 n+1
2
59
5. Define initial values of internal variables, to be used in return mapping algorithm, to the converged
values of the previous step at tn . Check if temperature softening is included by the user (i.e., m 6= 0
trial
in material input), compute a value for the yield function f (σn+1 ) and set the incremental plastic
multiplier ∆λn+1 equal to zero.

*--------------------------------------------------------------
PN = STATE(3)
pn = STATE(3) RP = STATE(6)
SIGMAY = SY + RP
Rn = STATE(6) *
σy,n = A + Rn * import old value damage variables
WE = STATE(14)
Wn = STATE(14) *
IF (TEMP_SOFT) THEN
TSTAR = (T-TR) / (TM-TR)
TSOFT = R1 - TSTAR**M
IF temperature softening is specified in material F = PHITRIAL - SIGMAY*TSOFT
ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1
input (m 6= 0) THEN F = PHITRIAL - SIGMAY
ENDIF
Tn − Tr DLAMBDA = R0
∗
T = DDLAMBDA = R0
Tm − Tr IF (PHITRIAL == R0) THEN

∗m
PHI = R1
Γp (Tn ) = 1 − T ELSE
PHI = PHITRIAL
trial
fn+1 = σeq,n+1 − σy,n Γp (Tn ) ENDIF
*
ELSE m = 0 and the CPU time can be reduced

by skipping the softening term
trial
fn+1 = σeq,n+1 − σy,n
ENDIF
∆λn+1 = 0
δλn+1 = 0
IF ϕtrial
n+1 = 0 THEN
ϕn+1 = 1
ELSE
trial trial
ϕn+1 = ϕn+1 = σeq,n+1
ENDIF
60
6. Check for plastic admissibility.

trial ) > 0 THEN apply return mapping by using the
IF f (σn+1 IF (F > R0) THEN
*
cutting plane algorithm (SOLU = 1) to find the incremental * plastic step: apply return mapping
* the cutting plane method uses values from the previous
change in the plastic multiplier ∆λn+1 . * iteration to update the plastic parameter, dlambda.
* thus, all subsequent values (until dlambda) are related to
Superscript i denotes the local iteration counter. * 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
configuration. Note that Eqs.(35) and (34) are used to calculate the gradients, i.e., σ12 = σ21 ,
σ23 = σ32 and σ13 = σ31 .

IF a 3D stress state is used THEN * compute the respective gradients of the von mises yield surface
RNV(2:6) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
n11 = ∂f /∂σ11 = (σ11 − 0.5(σ22 + σ33 ))/σeq CASE (0) ! NO PLANE STRESS CONDITION
RNV(1) = (STRESS(1) - R5*(STRESS(2) + STRESS(3))) / PHI
n22 = ∂f /∂σ22 = (σ22 − 0.5(σ33 + σ11 ))/σeq RNV(2) = (STRESS(2) - R5*(STRESS(3) + STRESS(1))) / PHI
n33 = ∂f /∂σ33 = (σ33 − 0.5(σ11 + σ33 ))/σeq RNV(3) = (STRESS(3) - R5*(STRESS(1) + STRESS(2))) / PHI
CASE (1) ! PLANE STRESS CONDITION (STRESS(3)=0)
n12 = ∂f /∂σ12 = 3σ12 /σeq RNV(1) = (STRESS(1) - R5*STRESS(2)) / PHI
RNV(2) = (STRESS(2) - R5*STRESS(1)) / PHI
n23 = ∂f /∂σ23 = 3σ23 /σeq CASE (2) ! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)
n31 = ∂f /∂σ31 = 3σ31 /σeq RNV(1) = STRESS(1) / PHI ! TO BE CHECKED ...
END SELECT
*
ELSEIF a plane stress state is used THEN SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
CASE (0) ! NO SHEAR STRESS IS PRESENT
n11 = ∂f /∂σ11 = (σ11 − 0.5(σ22 ))/σeq RNV(4) = R3*STRESS(4) / PHI
CASE (3) ! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT
n22 = ∂f /∂σ22 = (σ22 − 0.5(σ11 ))/σeq RNV(4) = R3*STRESS(4) / PHI
RNV(5) = R3*STRESS(5) / PHI
n33 = 0
n12 = ∂f /∂σ12 = 3σ12 /σeq END SELECT
n11 = ∂f /∂σ11 = σ11 /σeq
(ii) Compute hardening modulus HR,n+1 according to Eq. (42)

* compute hardening modulus
HR = QR1*CR1*EXP(-CR1*P) + QR2*CR2*EXP(-CR2*P)
i i i
HR,n+1 = Q1 C1 exp(−C1 pn+1 )+Q2 C2 exp(−C2 pn+1 )
61
(iii) Compute denominator in Eq. (90) using Eqs. (91), (92), (93) and (98)

IF a 3D stress state is used THEN * compute denominator in the residual derivative of dlambda
T 2 2 2
h i
n Cn = C1 (n11 ) + (n22 ) + (n33 ) 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))
+ 2C4 [n11 n22 + n22 n33 + n33 n11 ] CASE (1) ! PLANE STRESS CONDITION (RNV(3)=0)
h
2 2 2
i DENOMCP = C1*(RNV(1)*RNV(1) + RNV(2)*RNV(2))
+ 4C7 (n12 ) + (n23 ) + (n31 ) > + 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 ...
ELSEIF a plane stress state is used THEN END SELECT
*
SELECT CASE (ITAU) ! TREAT SHEAR COMPONENTS
T 2 2
h i
n Cn = C1 (n11 ) + (n22 ) CASE (0) ! NO SHEAR STRESS IS PRESENT
CASE (1) ! ONLY TAU_XY (RNV(4)) IS PRESENT
2
h i
+ 2C4 [n11 n22 ] + 4C7 (n12 ) 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)
ELSE a 1D stress state is used and > + RNV(6)*RNV(6))
END SELECT
*
T 2
n Cn = C1 (n11 ) 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)
C DENOM = DENOMCP + HR*DENOMV*TSOFT + SIGMAY*DENOMDV*TSOFT
∆λi

i n+1  ELSE ! M = 0 --> THE CODE ASSUMES TSTAR=0, SO THAT TSOFT=1
υp,n+1 = 1 +
ṗ0 ∆tn+1 DENOM = DENOMCP + HR*DENOMV + SIGMAY*DENOMDV
ENDIF
C−1
i ∆λi

∂υp,n+1 C n+1
= 1 + 
∂∆λi ṗ0 ∆tn+1 ṗ0 ∆tn+1
n+1
IF temperature softening is activated (m 6= 0)

THEN
" #m !
Tn − Tr
Γp (Tn ) = 1−
Tm − Tr
T i i
DENOM = n Cn + υp,n+1 Γp,n HR,n+1
i
∂υp,n+1
i
+ σy,n+1 Γp,n
∂∆λi
n+1
ELSE m = 0 and skip the softening term
T i i
DENOM = n Cn + υp,n+1 HR,n+1
i
∂υp,n+1
i
+ σy,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

* compute residual increment ddlambda and update variable dlambda
DDLAMBDA = F / DENOM
i
fn+1 DLAMBDA = DLAMBDA + DDLAMBDA
i+1
δλ =
n+1
DENOM
i+1 i i+1
∆λ = ∆λn+1 + δλ
n+1 n+1
62
(v) Update internal variables dependent on ∆λi+1
n+1

* update internal variables dependent on dlambda
P = P + DDLAMBDA
i+1 i i+1 R = RP + HR*DLAMBDA
p = pn+1 + δλ
n+1 n+1
SIGMAY = SY + R
i+1 i i+1
R = Rn + HR,n+1 ∆λ
n+1 n+1
i+1 i+1
σ = A+R
y,n+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

IF a 3D stress state is used THEN * update stress components using the elastic (hooke’s) law
i+1 i i+1
σ = σ11,n+1 − δλ (C1 n11 + C4 n22 + C5 n33 ) STRESS(1) = STRESS(1) -
11,n+1 n+1
> DDLAMBDA*(C1*RNV(1) + C4*RNV(2) + C5*RNV(3))
i+1 i i+1
σ = σ22,n+1 − δλ (C4 n11 + C2 n22 + C6 n33 ) STRESS(2) = STRESS(2) -
22,n+1 n+1
i+1 i i+1 STRESS(3) = STRESS(3) -
σ = σ33,n+1 − δλ (C5 n11 + C6 n22 + C3 n33 )
33,n+1 n+1 > DDLAMBDA*(C5*RNV(1) + C6*RNV(2) + C3*RNV(3))
i+1 i i+1 CASE (1) ! PLANE STRESS CONDITION (RNV(3)=0)
σ = σ12,n+1 − δλ C n
12,n+1 n+1 7 12 STRESS(1) = STRESS(1) -
i+1 i i+1 > DDLAMBDA*(C1*RNV(1) + C4*RNV(2))
σ = σ23,n+1 − δλ C n
23,n+1 n+1 8 23 STRESS(2) = STRESS(2) -
> DDLAMBDA*(C4*RNV(1) + C2*RNV(2))
i+1 i i+1
σ = σ31,n+1 − δλ C n CASE (2) ! UNIAXIAL STRESS CONDITION (RNV(2)=RNV(3)=0)
31,n+1 n+1 9 31
STRESS(1) = STRESS(1) - DDLAMBDA*C1*RNV(1) ! TO BE CHECKED ...
END SELECT
ELSEIF a plane stress state is used THEN *
i+1 i i+1 CASE (1) ! ONLY TAU_XY (RNV(4)) IS PRESENT
σ = σ11,n+1 − δλ (C1 n11 + C4 n22 + C5 n33 )
11,n+1 n+1 STRESS(4) = STRESS(4) - DDLAMBDA*C7*RNV(4)
i+1 i i+1 CASE (3) ! ALL THREE SHEAR STRESSES (RNV(4:6)) ARE PRESENT
σ = σ22,n+1 − δλ (C4 n11 + C2 n22 + C6 n33 )
22,n+1 n+1 STRESS(4) = STRESS(4) - DDLAMBDA*C7*RNV(4)
i+1 STRESS(5) = STRESS(5) - DDLAMBDA*C8*RNV(5)
σ = 0 STRESS(6) = STRESS(6) - DDLAMBDA*C9*RNV(6)
33,n+1
END SELECT
i+1 i i+1
σ = σ12,n+1 − δλ C n *
12,n+1 n+1 7 12

i+1 i i+1
σ = σ11,n+1 − δλ C n
11,n+1 n+1 1 11
63
(vii) Compute the updated von Mises equivalent stress using Eq. (35)

IF a 3D stress state is used THEN * compute updated equivalent stress
3 0 ,i+1 0 ,i+1 i+1 2 i+1 2 PHI = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)
σ : σ = (σ ) + (σ )
n+1 n+1 11,n+1 22,n+1 > + STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
2
> - STRESS(2)*STRESS(3) - STRESS(3)*STRESS(1)
i+1 2 i+1 i+1
+ (σ ) −σ σ CASE (1) ! PLANE STRESS CONDITION (STRESS(3)=0)
33,n+1 11,n+1 22,n+1
PHI = STRESS(1)*STRESS(1) + STRESS(2)*STRESS(2)
i+1 i+1 i+1 i+1 > - STRESS(1)*STRESS(2)
+σ σ −σ σ
22,n+1 33,n+1 33,n+1 11,n+1 CASE (2) ! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)
i+1 2 i+1 2 PHI = STRESS(1)*STRESS(1)
+ 3 (σ ) + (σ )
12,n+1 23,n+1 END SELECT
i+1 2 *
+ (σ )
31,n+1 SELECT CASE (ITAU) ! ADD SHEAR COMPONENTS IF PRESENT
ELSEIF a plane stress state is used THEN PHI = PHI + R3*(STRESS(4)*STRESS(4))
PHI = PHI + R3*(STRESS(4)*STRESS(4)
3 0 ,i+1 0 ,i+1 i+1 2 i+1 2
σ : σ = (σ ) + (σ ) > + STRESS(5)*STRESS(5)
n+1 n+1 11,n+1 22,n+1
2 > + STRESS(6)*STRESS(6))
i+1 i+1 i+1 2 END SELECT
−σ σ + 3 (σ ) *
11,n+1 22,n+1 12,n+1
PHI = SQRT (PHI)
3 0 ,i+1 0 ,i+1 i+1 2

σ : σ = (σ )
n+1 n+1 11,n+1
2
ENDIF s
i+1 i+1 3 0 ,i+1 0 ,i+1
σ = ϕ(σ ) = σ : σ
eq,n+1 n+1 n+1 n+1
2
(viii) Update the yield function in Eq. (84) according to the new stress state computed in (vi) and
(vii). Check for convergence.

IF temperature softening (m 6= 0) THEN * compute yield function and new residual
IF (TEMP_SOFT) THEN
TSTAR = (T - TR) / (TM - TR) ! THIS RECALCULATION OF TSTAR
i+1 i+1 i+1

TSOFT = R1 - TSTAR**M ! AND TSOFT SEEMS REDUNDANT HERE
f = ϕ −σ •
v,n+1 n+1 y,n+1 F = PHI - SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT
C *
i+1

∆λ Tn − Tr m
" # !

n+1  * check convergence
• = 1 + 1−

 RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT)
ṗ0 ∆tn+1 Tm − Tr RESNOR = ABS (RESNOR)
i+1
fv (∆λ )

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

RESNOR =

* check convergence
σ i+1 i+1

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

RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C))
RESNOR = ABS (RESNOR)
ENDIF
ELSE m = 0 and skip the softening term *
IF (RESNOR <= TOL) THEN
i+1 i+1 i+1

f = ϕ −σ •
v,n+1 n+1 y,n+1
C
i+1

∆λ
n+1 

• = 1 +


ṗ0 ∆tn+1

fv (∆λi+1 )

n+1
RESNOR =

σ i+1 i+1

y,n+1 υp,n+1

ENDIF
IF RESNOR < TOL THEN continue to 7.
ELSE i = i + 1 and GO TO (i)
64
7. Update the temperature according to the low coupling approach using Eq. (45) and erode element if
T > Tm

χσeq,n+1 ∆λn+1 * update temperature (only if temperature softening is activated)
Tn+1 = Tn + IF (TEMP_SOFT) THEN
ρCT 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)

* compute hydrostatic stress based on new stress state,
SIGMAH = STRESS(1)
1 1 SELECT CASE (IPLANC)
σH,n+1 = tr(σn+1 ) = (σ11,n+1 +σ22,n+1 +σ33,n+1 ) CASE (0) ! NO PLANE STRESS CONDITION
3 3
SIGMAH = SIGMAH + STRESS(2) + STRESS(3)
SIGMAH = SIGMAH + STRESS(2)
CASE (2) ! UNIAXIAL STRESS CONDITION (STRESS(2)=STRESS(3)=0)
END SELECT
*
SIGMAH = SIGMAH / R3
0
9. Computing the principal stresses σi,n+1 based on the principal deviatoric stresses σi,n+1 as shown in
Eqs. (55) and (57)

* 3d stress state and compute deviatoric stress
S1 = STRESS(1) - SIGMAH
0 S2 = STRESS(2) - SIGMAH
σ11,n+1 = σ11,n+1 − σH,n+1
0 *
σ22,n+1 = σ22,n+1 − σH,n+1 * finding principal stress based on the deviatoric stress and invariants
A1 = R5*(S1*S1 + S2*S2 + S3*S3)
0
σ33,n+1 = σ33,n+1 − σH,n+1 *
1 0 2 0 2 0 2

J2,n+1 = (σ11,n+1 ) + (σ22,n+1 ) + (σ33,n+1 ) A2 = -S1*S2*S3 ! DETERMINANT OF THE DEVIATORIC STRESS TENSOR
2 CASE (1) ! ONLY TAU_XY (STRESS(4)) IS PRESENT
2 2 2 A1 = A1 + STRESS(4)*STRESS(4)
+ (σ12,n+1 ) + (σ23,n+1 ) + (σ31,n+1 )
A2 = S3*STRESS(4)*STRESS(4) - S1*S2*S3
0 2 0 2 CASE (3) ! ALL THREE SHEAR STRESSES (STRESS(4:6)) ARE PRESENT
J3,n+1 = σ11,n+1 (σ23,n+1 ) + σ22,n+1 (σ31,n+1 )
A1 = A1 +
0 0 0 0 > STRESS(4)*STRESS(4) + STRESS(5)*STRESS(5) +
2
+σ33,n+1 (σ12,n+1 ) − σ11,n+1 σ22,n+1 + σ33,n+1 > STRESS(6)*STRESS(6)
A2 = S1*STRESS(5)*STRESS(5) +
− 2σ12,n+1 σ23,n+1 σ31,n+1 > S2*STRESS(6)*STRESS(6) +
  > S3*STRESS(4)*STRESS(4) -
 J3,n+1  > S1*S2*S3 - R2*STRESS(4)*STRESS(5)*STRESS(6)
θ = arccos 
 r
, 0 ≤ θ ≤ π END SELECT
2 J3

/27 *
2,n+1 A3 = -SQRT (R27 / A1)*A2*R5 / A1
s A3 = SIGN (MIN (ABS (A3), R1), A3)
0 J2,n+1 θ
σ1,n+1 = 2 cos A4 = ACOS (A3) / R3
3 3 *
s ! * compute (chronological) principal stresses
0 J2,n+1 θ 2π AUX = SQRT (A1 / R3)
σ2,n+1 = 2 cos −
3 3 3 S1 = SIGMAH + R2*AUX*COS(A4)
s S2 = SIGMAH + R2*AUX*COS(A4-R2*PI/R3) ! only for postprocessing
!
0 J2,n+1 θ 2π S3 = SIGMAH + R2*AUX*COS(A4+R2*PI/R3) ! only for postprocessing
σ3,n+1 = 2 cos + *
3 3 3
0
σ1,n+1 = σH,n+1 + σ1,n+1
0
σ2,n+1 = σH,n+1 + σ2,n+1
0
σ3,n+1 = σH,n+1 + σ3,n+1
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

*
WE = WE + MAX (S1, R0)*DLAMBDA
Wn+1 = Wn + max(σ1,n+1 , 0)∆λn+1 D = WE / WC
*
Wn+1
Dn+1 =
Wc
11. Update state variables to be returned to the element routine in EUROPLEXUS and check for element
erosion.

* update state variables before exit
STATE(1) = SIGMAH
ST AT E(1) = σH,n+1 STATE(2) = PHI
STATE(3) = P
ST AT E(2) = σeq,n+1 STATE(4) = PHITRIAL
STATE(5) = F
ST AT E(3) = pn+1 STATE(6) = R
trial STATE(7) = DDLAMBDA
ST AT E(4) = σeq,n+1 STATE(8) = DLAMBDA
STATE(9) = NRITER
ST AT E(5) = fn+1 STATE(10) = DLAMBDA / DT
STATE(11) = D
ST AT E(6) = Rn+1
*
ST AT E(7) = δλn+1 STATE(13) = T
STATE(14) = WE
ST AT E(8) = ∆λn+1 STATE(15) = CSON
*
ST AT E(9) = N RIT ER STATE(16) = S1
ST AT E(10) = ṗn+1 = ∆λn+1 /∆tn+1 STATE(17) = S2
STATE(18) = S3
ST AT E(11) = Dn+1 *
* decide if material point stiffness is zero, i.e. failed
ST AT E(13) = Tn+1 *--------------------------------------------------------------
990 IF (D >= DC) THEN
ST AT E(14) = Wn+1 *
ST AT E(15) = CSONn+1 * the current gp was already failed or is failing now :
* state(1:10,12,13), the relevant components
ST AT E(16) = σ1,n+1 * of stressnew and cson are set to 0.
* the total strains (strain) are not set to 0 so they continue to be
ST AT E(17) = σ2,n+1 * updated as long as the element is not eroded, i.e. as long
* as the element has some virgin gauss points.
ST AT E(18) = σ3,n+1 *
SELECT CASE (ITAU)
IF Dn+1 > DC THEN
CASE (0)
failure and the material point stiffness is deleted. STRESSNEW(1:3) = R0
CASE (1)
STATE(12) = 1 STRESSNEW(1:4) = R0
ELSE CASE (2)
STRESSNEW(1:5) = R0
no failure and the material point stiffness should CASE (3)
not be deleted. STRESSNEW(1:6) = R0
END SELECT
STATE(12) = 0 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

IF plane stress THEN use Eq.(69) to compute ac- * strain increment(s) and total strain(s) corresponding to
* the imposed zero stress condition(s)
tual strain trough thickness *
1 − 2ν CASE (0) ! NO PLANE STRESS CONDITION
∆ε33,n+1 = (∆σ11,n+1 + ∆σ22,n+1 ) CASE (1) ! PLANE STRESS CONDITION (STRESS(3)=0)
E * store provisional strain increment
− (∆ε11,n+1 + ∆ε22,n+1 ) AUX3 = STRAININC(3)
* compute and update the actual strain increment
ε33,n+1 = ε33,n + ∆ε33,n+1 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)
AUX2 = STRAININC(2)
AUX3 = STRAININC(3)
FACTOR = T2 / YOUNG
STRAININC(2) = 0.5D0*(FACTOR*(STRESSNEW(1) - STRESSOLD(1))
> - STRAININC(1))
STRAININC(3) = STRAININC(2)
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
stiffness in the material point is removed. An element will be removed from the mesh when a user-defined
fraction of the material points in the element are deleted. This is specified in the input file 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 = 12 γ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 defined 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)
* ----------------------------------------------------
* matval(1) - rho - material density
* matval(2) - young - elastic modulus
* matval(3) - nu - possion ratio
* matval(4) - elas - initial yield stress
* matval(5) - tol - tolerance newton iteration
* matval(6) - mxit - maximum number of iterations
* matval(7) - qr1 - asymptote first term voce hardening
* matval(8) - cr1 - hardening parameter first term voce
* matval(9) - qr2 - asymptote 2nd term voce hardening
* matval(10)- cr2 - hardening parameter 2nd term voce
* matval(11)- pdot - reference strain-rate
* matval(12)- c - hardening parameter visco term
* matval(13)- tq - taylor-quinney coefficient
* matval(14)- cp - specific heat capacity solid material
* matval(15)- tr - room temperature
68
* matval(16)- tm - melting temperature
* matval(17)- m - hardening parameter temperature term
* matval(18)- dc - upper limit of d when softening
* matval(19)- wc - cockcroft-latham failure criterion
* matval(20)- solu - 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 ’,
> ’PDOT’,’C ’,’TQ ’,’CP ’,’TR ’,
> ’TM ’,’M ’,’DC ’,’WC ’,’SOLU’/
*
* 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 ! TOL
IF (KOPT(6) == 0) NEWMAT%MATVAL(6) = 20 ! MXIT
IF (KOPT(15) == 0) NEWMAT%MATVAL(15) = 293.D0 ! TR
IF (KOPT(20) == 0) NEWMAT%MATVAL(20) = 1.D0 ! SOLU
*
* printouts
*
WRITE(TAPOUT,1000) NEWMAT%NUMLDC
*
WRITE(TAPOUT,1001) NEWMAT%MATVAL(1)
WRITE(TAPOUT,1006) NINT (NEWMAT%MATVAL(6))
69
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)
*
* vegard aune, folco casadei 12-14
*
* added optional radial return solution strategy (solu == 2) 02-16
*
*erosion
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
*
* input parameters to generalize the present material routine:
* itau : number of shear stress components (stress(4),stress(5),stress(6))
* 0 = no shear stress is present
* 1 = only tau_xy (stress(4)) is present
* 3 = all three shear stresses (stress(4:6)) are present
* iplanc : 0 = no plane stress condition
* 1 = plane stress condition (sigma_z=stress(3)=0)
* 2 = uniaxial stress condition (sigma_y=stress(2)=0 and
* sigma_z=stress(3)=0)
*
* organization of stress and strain arrays
* ----------------------------------------
* stress(1)- xx : normal components
* stress(2)- yy
* stress(3)- zz
* stress(4)- xy : shear components
* stress(5)- yz
* stress(6)- xz
*
* attention: the shear components of strain and straininc arrays are the
* "engineering" values (gamma), i.e. twice the "tensor" values (epsilon)!
*
* list of material properties sent to the subroutine
* --------------------------------------------------
* props(1) - rho - material density
* props(2) - young - elastic modulus
* props(3) - poiss - possion ratio
* props(4) - sy - initial yield stress
* props(5) - tol - tolerance newton iteration
* props(6) - mxiter - maximum number of iterations
* props(7) - qr1 - asymptote first term voce hardening
* props(8) - cr1 - hardening parameter first term voce
* props(9) - qr2 - asymptote 2nd term voce hardening
* props(10)- cr2 - hardening parameter 2nd term voce
* props(11)- pdot - reference strain-rate
* props(12)- c - hardening parameter visco term
* props(13)- tq - taylor-quinney coefficient
* props(14)- cp - specific heat capacity solid material
* props(15)- tr - room temperature
* props(16)- tm - melting temperature
* props(17)- m - hardening parameter temperature term
* props(18)- dc - upper limit of d when softening
* props(19)- wc - cockcroft-latham failure criterion
* props(20)- solu - 1=cutting plane (default) requires N-R iterations,
* 2=radial return (direct, no iterations)
* (not available in plane or uniaxial stress!)
*
* list of state variables (old and new)
* ----------------------------------------------------------------
* description name intent
* ----------------------------------------------------------------
* state(1) - hydrostatic pressure (1/3 sigma_kk) sigmah new
* state(2) - von mises equivalent stress phi new
* state(3) - equivalent plastic strain p old+new
* state(4) - elastic trial equivalent stress phitrial new
* state(5) - yield function (should be near zero) f new
* state(6) - total hardening of the material r old+new
* state(7) - change of dlambda over the step ddlambda new
* state(8) - incremental plastic multiplier dlambda new
* state(9) - n of iterations for convergence nriter new
* state(10)- plastic strain rate dlambda / dt new
* state(11)- damage d old+new
* state(12)- failure: 1 = virgin g.p., 0 = failed g.p. fa new
* state(13)- absolute temperature t old+new
* state(14)- cockcroft-latham damage accumulation we old+new
* state(15)- sound speed cson new
* state(16)- first principal stress s1 new
* state(17)- second principal stress s2 new
* state(18)- third principal stress s3 new
71
*
* fail_flag : output only!
* 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
* (i.e. l_erosion == .true.)
*
* NOTE:
* - the material routine is called also for an already failed gp, if the
* element to which the gp belongs has not been eroded yet!
* - however, the material routine is not called at all for an eroded element!
*
! Thread variable
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 = PROPS(1)
YOUNG = PROPS(2)
POISS = PROPS(3)
SY = PROPS(4)
TOL = PROPS(5)
MXITER = NINT (PROPS(6))
QR1 = PROPS(7)
CR1 = PROPS(8)
QR2 = PROPS(9)
CR2 = PROPS(10)
PDOT = PROPS(11)
C = PROPS(12)
TQ = PROPS(13)
CP = PROPS(14)
TR = PROPS(15)
TM = PROPS(16)
M = PROPS(17)
WC = PROPS(19)
SOLU = 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
*
CSON = SQRT (YOUNG*T3 / (RHO*T1*T2))
CASE (1) ! PLANE STRESS CONDITION (SIG(3)=0)
CSON = SQRT (YOUNG / (RHO*T1*T3))
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
C1 = YOUNG*T3 / (T1*T2)
C4 = POISS*C1 / T3
CASE (1) ! PLANE STRESS CONDITION (SIG(3)=0)
C1 = YOUNG / (T1*T3)
C4 = POISS*C1
C1 = YOUNG
C4 = R0 ! THESE VALUES HAVE TO BE CHECKED
END SELECT
*
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
*
> + C5*STRAININC(3)
*
> + C6*STRAININC(3)
> + C3*STRAININC(3)
> + STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
> + C6*STRAININC(3)
STRESS(3) = R0
> - STRESS(1)*STRESS(2)
STRESS(2) = R0
STRESS(3) = R0
PHITRIAL = STRESS(1)*STRESS(1)
END SELECT
*
STRESS(4) = STRESSOLD(4) + C7*STRAININC(4) ! THE GAMMAS ARE USED!
PHITRIAL = PHITRIAL + R3*(STRESS(4)*STRESS(4))
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(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)
F = PHITRIAL - SIGMAY*TSOFT
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(1) = STRESS(1) / PHI ! TO BE CHECKED ...
END SELECT
*
END SELECT
*
* compute hardening modulus
*
* compute denominator in the residual derivative of dlambda
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))
DENOMCP = C1*(RNV(1)*RNV(1)) ! TO BE CHECKED ...
END SELECT
*
DENOMCP = DENOMCP + R4*C7*(RNV(4)*RNV(4))
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
DENOM = DENOMCP + HR*DENOMV*TSOFT + SIGMAY*DENOMDV*TSOFT
DENOM = DENOMCP + HR*DENOMV + SIGMAY*DENOMDV
ENDIF
*
* compute residual increment ddlambda and update variable dlambda
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
STRESS(1) = STRESS(1) -
STRESS(1) = STRESS(1) - DDLAMBDA*C1*RNV(1) ! TO BE CHECKED ...
END SELECT
*
STRESS(4) = STRESS(4) - DDLAMBDA*C7*RNV(4)
END SELECT
*
* compute updated equivalent stress
> + STRESS(3)*STRESS(3) - STRESS(1)*STRESS(2)
> - STRESS(1)*STRESS(2)
PHI = STRESS(1)*STRESS(1)
END SELECT
*
75
PHI = PHI + R3*(STRESS(4)*STRESS(4))
PHI = PHI + R3*(STRESS(4)*STRESS(4)
> + 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 ! AND TSOFT SEEMS REDUNDANT HERE
F = PHI - SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT
*
* check convergence
RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)*TSOFT)
F = PHI - SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C)
* check convergence
RESNOR = F / (SIGMAY*((R1 + DLAMBDA / (PDOT*DT))**C))
ENDIF
*
*
* update stress components based on plastic step
*-----------------------------------------------
STRESSNEW(1) = STRESS(1)
STRESSNEW(2:4) = R0 ! INITIALIZE TO ZERO "OPTIONAL" COMPONENTS
END SELECT
*
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)
STATE(13) = TM
GO TO 990
ENDIF
ENDIF
*
* compute hydrostatic stress based on new stress state,
SIGMAH = STRESS(1)
END SELECT
*
*
76
*
* finding principal stress based on the deviatoric stress and invariants
A1 = R5*(S1*S1 + S2*S2 + S3*S3)
*
A2 = -S1*S2*S3 ! DETERMINANT OF THE DEVIATORIC STRESS TENSOR
A1 = A1 + STRESS(4)*STRESS(4)
A1 = A1 +
> STRESS(4)*STRESS(4) + STRESS(5)*STRESS(5) +
> STRESS(6)*STRESS(6)
> 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
*
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
77
END SELECT
*
END SELECT
*
SIGMAH = STRESS(1)
END SELECT
*
*
*
A1 = R5*(S1*S1 + S2*S2 + S3*S3)
*
A1 = A1 +
END SELECT
*
A3 = -SQRT (R27 / A1)*A2*R5 / A1
A4 = ACOS (A3) / R3
*
*
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
* THIS CASE IS ILLEGAL FOR SOLU == 2 (SEE TEST AT BEGINNING OF ROUTINE)
RNV(1) = STRESS(1) / PHI ! TO BE CHECKED ...
END SELECT
*
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
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 = R0
TSOFT = R1
ENDIF
DENOM = R3G + HR*DENOMV*TSOFT + SIGMAY*DENOMDV*TSOFT
C COMPUTE NEWTON-RAPHSON INCREMENT AND UPDATE VARIABLE DLAMBDA
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)
*
*
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)
STATE(13) = TM
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)
END SELECT
C COMPUTE DEVIATORIC STRESS
C FINDING PRINCIPAL STRESS BASED ON THE DEVIATORIC STRESS AND INVARIANTS
A1 = R5*(S1*S1 + S2*S2 + S3*S3)
A1 = A1 +
END SELECT
A3 = -SQRT (R27 / A1)*A2*R5 / A1
A4 = ACOS (A3) / R3
C COMPUTE (CHRONOLOGICAL) PRINCIPAL STRESSES
D = WE / WC
*
C UPDATE STRESS COMPONENTS USING THE ELASTIC (HOOKE’S) LAW
C--------------------------------------------------------------
STRESSNEW(1) = STRESS(1) -
> DLAMBDA*(C1*RNV(1) + C4*RNV(2) + C5*RNV(3))
> DLAMBDA*(C1*RNV(1) + C4*RNV(2))
> DLAMBDA*(C4*RNV(1) + C2*RNV(2))
> DLAMBDA*C1*RNV(1) ! TO BE CHECKED ...
END SELECT
80
*
STRESSNEW(4) = STRESS(4) - DLAMBDA*C7*RNV(4)
END SELECT
*
C COMPUTE UPDATED EQUIVALENT STRESS
PHI = STRESSNEW(1)*STRESSNEW(1) +
> STRESSNEW(2)*STRESSNEW(2) +
> STRESSNEW(3)*STRESSNEW(3) -
> STRESSNEW(3)*STRESSNEW(1)
PHI = STRESSNEW(1)*STRESSNEW(1) +
> STRESSNEW(1)*STRESSNEW(2)
PHI = STRESSNEW(1)*STRESSNEW(1)
END SELECT
*
CASE (1) ! ONLY TAU_XY (STRESSNEW(4)) IS PRESENT
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
81
END SELECT
*
END SELECT
*
SIGMAH = STRESS(1)
END SELECT
*
*
A1 = R5*(S1*S1 + S2*S2 + S3*S3)
*
A1 = A1 +
END SELECT
*
A3 = -SQRT (R27 / A1)*A2*R5 / A1
A4 = ACOS (A3) / R3
*
*
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)
*
* 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
AUX2 = STRAININC(2)
AUX3 = STRAININC(3)
FACTOR = T2 / YOUNG
STRAININC(2) = 0.5D0*(FACTOR*(STRESSNEW(1) - STRESSOLD(1))
> - STRAININC(1))
STRAININC(3) = STRAININC(2)
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.
LB-NA-27982-EN-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

Johnson Cook Fortran Code

Uploaded by

Copyright:

Available Formats

Johnson Cook Fortran Code

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Johnson Cook Fortran Code

Uploaded by

Copyright:

Available Formats

What is the document about?

What is the document about?

What material model is implemented in EUROPLEXUS?

What material model is implemented in EUROPLEXUS?

Formulation and Implementation of

the VPJC Material Model in

JRC Science Hub

© European Union, 2016

Reproduction is authorised provided the source is acknowledged.

All images © European Union 2016

How to cite: Authors; title; EUR; doi

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

European Laboratory for Structural Assessment (ELSA), 21027 Ispra, Italy

ε̇ = ε̇e + ε̇p (1)

C = λel I ⊗ I + 2µel I (3)

In general, assuming only isotropic hardening, the yield criterion is defined as

f (σ, R) = ϕ(σ) − σy (σ0 , R) = σeq − (σ0 + R) (5)

ϕ(aσ) = aϕ(σ) (6)

Dp = σ ε̇p = σeq ṗ (15)

f (σ) = σv (ṗ) > 0 ⇐⇒ fv (σ) = ϕ(σ) − σy (σ0 , R) − σv (ṗ) = 0 (19)

Figure 2: Viscoplastic effects due to strain-rate sensitivity.

σ̂ ∇ = C : De = λel (trDe ) I + 2µel De (23)

σ̇ = Aσ − σA + λel (trDe ) I + 2µel De (25)

σ̂ ∇G = σ̇ − ṘRT σ + σ ṘRT = σ̇ − Ωσ + σΩ, A = Ω = ṘRT = −RṘT (27)

(a) Principal stress space (b) The deviatoric plane (π-plane)

Figure 4: The von Mises yield criterion.

2.2. The modified Johnson-Cook constitutive model

σy = (A + Bpn ) (1 + Cln(ṗ∗ )) (1 − T ∗ m ) (38)

2.3. Adiabatic thermal softening

2.4. Ductile failure

Plastic failure strain, pf

3. Two-dimensional stress analysis

3.1. Two-dimensional theory of elasticity

3.2. Through thickness strain increment

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

ε33 = − (ε11 + ε22 ) (64)

∆εp33 = −(∆εp11 + ∆εp22 ) (65)

∆εp11 = ∆ε11 − ∆εe11

∆εp33 = −(∆ε11 − ∆εe11 ) − (∆ε22 − ∆εe22 )

∆ε33 = ∆εe33 + ∆εp33

4. Uniaxial stress case

∆εp11 = ∆ε11 − ∆εe11

∆εp11 + ∆εp22 + ∆εp33 = 0 (71)

∆εe22 = ∆εe33 = −ν∆εe11 (74)

∆εp22 = ∆εp33 (75)

∆ε22 = ∆ε33 (76)

ε22 = ε33 (77)

2∆εp22 = 2∆εp33 = −∆εp11 (78)

∆ε22 = ∆ε33 = ∆εe22 + ∆εp22 = ∆εe33 + ∆εp33

∆σ = C : (∆ε − ∆εp ) = ∆σ e + ∆σ p (82)

5.1. Cutting plane method

i+1 ∂f ∂f ∂υp,n+1 i+1

i+1 ∂f ∂f ∂υp,n+1 i+1

2. Compute the elasticity tensor.

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

ELSE a 1D stress state is used THEN use

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

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

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

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

∂fv (∆λn+1 ) ∂ h 0 trial i