TASK QUARTERLY 9 No 2, 157166

NUMERICAL MODELING OF
ELASTO-VISCOPLASTIC CHABOCHE
CONSTITUTIVE EQUATIONS USING
MSC.MARC
ANDRZEJ AMBROZIAK
Department of Structural Mechanics,
Faculty of Civil and Environmental Engineering,
Gdansk University of Technology,
Narutowicza 11/12, 80-952 Gdansk, Poland
[email protected]
(Received 22 October 2004; revised manuscript received 17 January 2005)
Abstract: The aim of the present paper is to propose a special kind of finite element procedure
for dynamic and static analysis of civil engineering structures (e.g. trusses, beams, shells and spatial
structures) including elasto-viscoplastic models. The Chaboche model with damage has been chosen
from the wide range of available elasto-viscoplastic constitutive models for numerical calculations.
The main advantage of the presented approach is the possibility of interference in subroutines and
adjusting them to a particular problem (e.g. altering the form of constitutive equations). A userdefined UVSCPL subroutine has been proposed to introduce the elasto-viscoplastic model into the
MSC.Marc system.
Keywords: elasto-viscoplastic constitutive models, damage, finite element method

1. Introduction
In materials where the viscous properties within the plastic range play an
important role it is necessary to apply the elasto-viscoplastic constitutive relations. We
can cite the following elasto-viscoplastic models for example: Perzyna [1], Chaboche
[2], Aubertin [3], Lehmann-Imatani [4], Miller [5], Bodner-Partom [6], Krempl [7],
Tanimura [8], Krieg-Swearengen-Jones [9], Walker [10], Korhonen-Hannjula-Li [11],
Freed-Virrilli [12]. A detailed description of these models can be found in paper [13].
Their relatively large number shows that they are not universal and can be applied
only under certain conditions or for a limited range of materials. From the wide
range of available elasto-viscoplastic laws, the Chaboche model with damage has
been selected in the present paper. This constitutive model is suitable for analysis
of isotropic materials (mainly metals at high temperatures) in the range of small
inelastic strains.

tq209b-e/157

11 V 2005

BOP s.c., http://www.bop.com.pl

158

A. Ambroziak

2. Description of the applied programs

Two computer programs have been applied in the present work. One has
been a self-made FORTRAN code (OFC) for elasto-viscoplastic truss analysis. The
theoretical description, FEM coding, details and application limits of this code are
discussed in [14]. The other is the MSC.Marc system developed by MSC Software,
a multipurpose finite element program for advanced engineering simulations. A great
advantage of this system is the possibility to apply over 100 user-subroutines, which
can be modified by the user. The standard MSC.Marc system does not support the
Chaboche model with the damage concept proposed by Lemaitre (see e.g. [15]).
Therefore, a user-defined UVSCPL subroutine [16] has been used to introduce the
Chaboche model with continuous damage into the MSC.Marc system.

3. The Chaboche model

3.1. Basic equations of the Chaboche model with damage
In the first step, we assume an additive decomposition of the strain rates into
the elastic E and inelastic I rate parts in the following form:
= E + I ,

(1)

so that it is necessary to assume small strains. The relation between the stress and
strain rates can be described for an isotropic material as:


(2)
= (1 D) B : E = B : I ,

where B is the tensor of elastic modules.

In the Chaboche model, the inelastic strain rate has the form of:
3
S0 X0
I = p
,
2 J(S0 X0 )

(3)

where p,
S0 and X0 are the rate of the equivalent plastic strain (the dots denote
differentiation with respect to time) and the deviatoric parts of stress and back stress
tensors, respectively. The J(a) invariant is defined by the following equation:
r
r
3
3 ij
J(a) =
a:a=
a aij .
(4)
2
2
We will consider the isotropic damage expressed by the scalar parameter 0 D 1.
The p rate in the Chaboche model with damage is defined in [15] as:
+n
* J(S0 X0 )
Rk
1D
p =
,
(5)
K
where k, D and K, n are the initial yield stress, the damage parameter and material
parameters, respectively. The angle bracket hxi has the following form:
1
hxi = (x + |x|)
2
and is called the McCauley bracket.

tq209b-e/158

11 V 2005

BOP s.c., http://www.bop.com.pl

(6)

Numerical Modeling of Elasto-Viscoplastic Chaboche Constitutive Equations. . .

159

The isotropic hardening, R, and kinematic hardening, X, are defined by the

following expressions:

(7)
R = b(R1 R) p,
2
= a I cX p,
X

(8)
3
where a, c, b and R1 are material parameters.
According to the damage model proposed by Kachanov [17], the damage
evolution D used in Equation (5) is expressed by the following equation (see [15]):
s

Y

D=
p,

(9)
S
where S, s are the material parameters of the damage and the Y function is given as:


2
1
2
2
(1 + )eq + 3(1 2)H ,
(10)
Y =
2(1 D)2 E 3

where E is Youngs modulus, is Poissons ratio, and eq and H are the HuberMisses equivalent stress and hydrostatic stress.
In the case of an uniaxial tension or compression state, the inelastic strain rate
can be written as:
+n
* |X|

k
1D
I = p sign( X) =
sign( X)
(11)
K
while kinematic hardening has a simple form, [2]:
X = a I cX | I |.

(12)

It is worth noting that the stress tensor, S, has one non-vanishing component,
, while the S0 and X0 tensors have three non-zero components:

2
2

0
0
X
0
0
0 0

3
3

S = 0 0 0 , S0 = 0 13
0 , X0 = 0
0 . (13)
13 X
1
1
0
0
3
0
0
3X
0 0 0

3.2. UVSCPL the calculation algorithm

The UVSCPL [16] routine allows very general material laws to be entered into
MSC.Marc. This user subroutine has been used in the present work for computing
the inelastic strain increment for the elastic-viscoplastic Chaboche material model.
In general, the inelastic strain rate and the stress increment must be defined in the
procedure. The main part of the algorithm used in the computation is presented in
Figure 1 (where j is the number of time steps). As in expressions for step j, the
components from the same step are used and iterations are necessary at each step.
Additionally, the PLOTV [16] routine has been used to visualize the evolution
of the damage parameter. In this procedure, it is simple to define the variable to be
plotted and written to the post file.

4. Numerical examples
Numerical calculations have been performed for simple solid (Figure 2) and
truss (Figure 3) structures. Such tests are usually performed in a laboratory during

tq209b-e/159

11 V 2005

BOP s.c., http://www.bop.com.pl

160

A. Ambroziak

t 
j ,
Xj1 + X
2

t 

[step 1]
R = 2 Rj1 + Rj ,


t 
j ,
Dj1 + D
D =
2

[step 2] S0j , X0j , J(S0j X0j ), J(S0j ),

X =

"

"

S0j X0j
3
p j
2 J(S0j X0j )

[step 3] p j =

[step 4] Ij =

Xj = Xj1 + X

Rj = Rj1 + R

Dj = Dj1 + D
tr(Sj )

J(S0j X0j )/(1 Dj ) Rj k

K

n #



j = 2 a I cXj p j
[step 5] X
3 j
i
h
j = b(R1 Rj ) p j
[step 6] R


s 
j = Y p j
[step 7] D
S


[step 8] Bj = (1 Dj )B



[step 9] Ij = Ij tj



[step 10] Sj = Bj (j Ij )

Figure 1. Flow graph of the UVSCPL subroutine

identification of the Chaboche model material parameters [18]. The identifications

results can be immediately applied in 3D or 2D analysis, on the assumptions of small
strains and isotropy of the material.
Three-dimensional eight-node isoparametric solid (Element 7, see [19]) and
three-dimensional two-node straight truss (Element 9, see [19]) elements were applied
in the MSC.Marc calculations. Trilinear interpolation functions were used in the former

Figure 2. Geometry of the solid structure

tq209b-e/160

11 V 2005

BOP s.c., http://www.bop.com.pl

Numerical Modeling of Elasto-Viscoplastic Chaboche Constitutive Equations. . .

161

Figure 3. Geometry of the truss structure

and thus the strains tend to be constant throughout the element. Linear interpolation
was used in the latter for coordinates and displacements. At the same time, only the
straight two-nodes truss element with simple linear interpolation was used in the
OFC program.
The geometrical parameters of l = 1.0m (length) and A = 0.001m2 (square
cross section area for the truss and solid elements) were assumed. Due to the purely
numerical character of the work the values of the material parameters were taken
from the literature. The author is aware that sometimes different sets of parameters
can be found in the literature for the same material under the same conditions due
to the strong physical nonlinearity of the problem (see [20]). In practical applications
such sets must be carefully verified.

4.1. Example 1 (without damage)

The following material parameters were assumed in the first step of the numerical calculations (INCO718 nickel-based alloy at 650C from [21]): E = 159.0GPa,
= 0.3, k = 514.21MPa, b = 60.0, R1 = 194.39MPa, a = 170.0GPa, c = 500.0, n = 4.0,
K = 1023.5 (MPa s)1/n , S = 0.0MPa, s = 0.0.
The stress versus strain diagrams for various strain rates (107 101 s1 )
are presented in Figure 4. The result of the creep test is shown in Figure 5 (the
constant stress condition applied at the level of = 1000MPa). There is no difference

Figure 4. Constant stress rate tests

tq209b-e/161

11 V 2005

BOP s.c., http://www.bop.com.pl

162

A. Ambroziak

Figure 5. Creep test

Figure 6. Constant stress rate tests

between the results obtained from the two computer programs and the types of finite
elements used in the numerical calculations.

4.2. Example 2 (with damage)

The following material parameters were taken for the Chaboche model coupled
to damage (INCO718 alloy at 627C from [18]): E = 162.0GPa, = 0.3, k = 501.0MPa,
b = 15.0, R1 = 165.4MPa, a = 80.0GPa, c = 200.0, n = 2.4, K = 12790.0 (MPas)1/n ,
S = 4.48MPa, s = 3.0. The reference calculations for this case were conducted with
damage constants S = 0.0MPa, s = 0.0.
The stress versus strain graphs for the two strain rates, = 0.001s1 and
= 0.01s1 , are shown in Figure 6, where the influence of damage is apparent. In

tq209b-e/162

11 V 2005

BOP s.c., http://www.bop.com.pl

Numerical Modeling of Elasto-Viscoplastic Chaboche Constitutive Equations. . .

163

Figure 7. Change of the damage parameter in the constant stress rate tests

Figure 8. Creep test

the analysed range of strain and stress, damage of the structure and a break in the
numerical calculations was observed in the case of the = 0.01s1 computation. The
evolution of the damage parameter is given in Figure 7; for the = 0.01s1 calculation
damage exceeds the value limit. The simulations of creep tests are presented in
Figures 8 and 9 and for the stress values of = 1000MPa and = 1500MPa. In
this case, the damage has great influence for = 2000MPa. Finally, a diagram of the
cyclic load test for = 0.001s1 is shown in Figure 10. The cyclic test was carried out
in the range of = 0.03. Like in the first example, there is no significant difference
between the results obtained from the two computer programs. It is worth pointing out
the additive character of the damage parameter. Its value increases at each successive
step, expect for the elastic range of deformation, where its value is constant due to

tq209b-e/163

11 V 2005

BOP s.c., http://www.bop.com.pl

164

A. Ambroziak

Figure 9. Change of the damage parameter in the creep tests

Figure 10. Cyclic test

the lack of inelastic effects (see Figures 11 and 12). The failure of the specimen after
a single cycle of loading is attributable to the relatively high strain rate. An error in
the determination of material constants is another possible explanation. For example,
the material parameter K has a very high value compared to other references in the
literature (cf. e.g. [22, 23]).

5. Conclusions
The following conclusion and remarks may now be formulated:
(1) Thanks to the possibility of including users constitutive model subroutines
into the MSC.Marc system, the Chaboche model with damage can be directly
applied in calculations.

tq209b-e/164

11 V 2005

BOP s.c., http://www.bop.com.pl

Numerical Modeling of Elasto-Viscoplastic Chaboche Constitutive Equations. . .

165

Figure 11. Change of the damage parameter of the cyclic test in the time domain

Figure 12. Change of the damage parameter of the cyclic test in the strain domain

(2) A non-linear analysis of simple solid and truss structures with elastic-viscoplastic physical equations for the Chaboche model has been successfully carried
out.
(3) The MSC.Marc and OFC programs respond almost identically in all the calculated test variants.
(4) The obtained results require a revision for more complex structures, types of
analysis and kinds of load. Especially the influence of layering in the shell
elements should be investigated.
(5) Further research should focus on the development of the UVSCPL routine with
an emphasized influence of temperature.

tq209b-e/165

11 V 2005

BOP s.c., http://www.bop.com.pl

166

A. Ambroziak

Acknowledgements
The research was performed within the Polonium 2005 Polish-French coopera
tion program administrated by KBN and EGIDE
(KBN 5598.II/2004/2005).
Calculations presented in the paper were performed at the Academic Computer
Centre in Gdansk (TASK).
The study was supported by the European Community under the FP5 Programme, key-action City of Tomorrow and Cultural Heritage (Contract No. EVK4CT-2002-80005). This support is acknowledged with gratitude.

References
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]

[22]
[23]

Perzyna P 1978 Thermodynamics of Elastic Materials, PWN, Warsaw (in Polish)

Chaboche J L 1989 Int. J. Plasticity 5 247
Aubertin M, Gill D E and Landanyi B 1991 Mechanics of Materials 11 63
Lehmann T 1984 Cours and Lectures, Vol. 281, Springer, Wien New York
Miller A 1976 J. Engng. Mat. and Technol., ASME 98 97; ibid. 106
Bodner S R and Partom Y 1975 J. Appl. Mech., ASME 42 385
Krempl E, McMahon J J and Yao D 1984 Mechanics of Material 5 35
Tanimura S 1979 J. Int. Engng. Sci. 17 997
Krieg R D, Swearengen J C and Jones W B 1987 Unified Constitutive Equations for Creep
and Plasticity (Miller K, Ed.), Elsevier, pp. 245271
Freed A D and Walker K P 1989 The Winter Annual Meeting of the ASME (Hui D and
Kozik T J, Eds.), San Francisco, pp. 1015
Korhonen M A, Hannjula S P and Li C Y 1987 Unified Constitutive Equations for Creep and
Plasticity (Miller K, Ed.), Elsevier, pp. 89138
Freed A D and Virrilli M J 1988 Proc. of the MECAMAT, Besancon, Vol. I, pp. 2739
Woznica K 1998 Dynamique des structures elasto-viscoplastique, Cahiers de Mechanique, Lille
Kosowski P and Woznica K 2000 Machine Dynamics Problems 24 (3) 33
Lemaitre J 1985 J. Engng. Mat. and Technol. 107 83
Users handbook: MSC.Marc, Volume D: User Subroutines and Special Routines, Version 2003.
MSC.Software Corporation 2003
Kachanov L M 1958 TVZ Akad. Nauk. S.S.R. Otd. Tech. Nauk. 8 26
Ammar G and Dufailly J 1993 European J. Mech., A/Solids 12 197
Users handbook: MSC.Marc, Volume B: Element Library, Version 2003. MSC.Software Corporation 2003
Kosowski P, Schmidt R and Weichert D 1997 Engineering Translations 45 91
Golinval J C 1989 Calcul par elements finis de structures elasto-viscoplastiques somies a
`
des chargements cycliques a
` haute temperature, Collection des publications de la Faculte des
Sciences Appliques, Universitete de Li`ege, No. 126b
Efits J, Abdel Kader M S and Jones D I 1989 Int. J. Plasticity 6 1
Lewandowski M S and Overfelt R A 1999 Acta Mater. 47 4695

tq209b-e/166

11 V 2005

BOP s.c., http://www.bop.com.pl