Available online at www.sciencedirect.

com

Draft
Draft Draft
Draft Draft
Draft Structural Integrity Procedia 00 (2016) 000000
www.elsevier.com/locate/procedia

21st European Conference on Fracture, ECF21, 20-24 June 2016, Catania, Italy

A continuum damage model for creep fracture and fatigue analyses

Petteri Kauppilaa , Reijo Kouhiab,, Juha Ojanperaa , Timo Saksalab , Timo Sorjonena
a Valmet Technologies Oy, P.O. Box 109, FI-33101 Tampere, Finland
b Tampere University of Technology, Department of Mechanical Engineering and Industrial Systems, P.O. Box 589, FI-33101 Tampere, Finland

Abstract
In this paper a thermodynamically consistent formulation for creep and creep-damage modelling is given. The model is developed
for isotropic solids by using proper expressions for the Helmholtz free energy and the complementary form of the dissipation
potential, and can be proven to fullfil the dissipation inequality. Also the coupled energy equation is derived. Continuum damage
model with scalar damage variable is used to facilitate simulations with tertiary creep phase. The complementary dissipation
potential is written in terms of the thermodynamic forces dual to the dissipative variables of creep strain-rate and damage-rate. The
model accounts for the multiaxial stress state and the difference in creep rupture time in shear and axial loading as well as in tensile
and compressive axial stress. In addition, the model is simple and only four to eight material model parameters are required in
addition to the elasticity parameters. A specific version of the proposed model is obtained when constrained to obey the Monkman-
Grant relationship between the minimum creep strain-rate and the creep rupture time. The applicability of the Monkman-Grant
hypothesis in the model development is discussed. The proposed 3D-model is implemented in the ANSYS finite element software
by the USERMAT subroutine. Material parameters have been estimated for the 7CrMoVTiB10-10 steel (T24) for temperatures
ranging from 500 to 600 degrees of celcius. Some test cases with cyclic thermal fatigue analysis are presented.
c 2016 The Authors. Published by Elsevier B.V.
Peer-review under responsibility of the Scientific Committee of ECF21.

Keywords: creep fatigue; damage mechanics; Ansys USERMAT

1. Introduction

Creep is an important deformation mechanism for metal structures operating at temperatures over 30 % of their
absolute melting temperature. Many technical applications require higher operating temperatures. The basic test in
studying creep deformation and fracture is the uniaxial creep test where a smooth tensile bar is subjected to a con-
stant load. As already observed by da Costa Andrade (1910), the experimental creep curve consists of three phases
corresponding to decreasing, constant and increasing strain rate. These phases are termed as primary, secondary and
tertiary creep stages. Strain hardening and thermally activated recovery of the dislocation structure are the main mech-
anisms in the primary and secondary phases, while formation of grain boundary cavities and changes in a dislocation
microstructure can be ascribed to the tertiary phase.

Corresponding author. Tel.: +358-40-849-0561; fax: +0-000-000-0000.

E-mail address: reijo.kouhia@tut.fi

2452-3216 c 2016 The Authors. Published by Elsevier B.V.

Peer-review under responsibility of the Scientific Committee of ECF21.
2 P. Kauppila et al. / Structural Integrity Procedia 00 (2016) 000000

Developing creep models applicable for large stress and temperature range is a complicated task. Continuum dam-
age mechanics provides a promising path to such goal, as can be seen e.g. from Altenbach et al. (2009); Betten (2005);
Gorash (2008); Hayhurst (1972, 1994); Murakami (2012); Naumenko (2006). In this study, a thermodynamically con-
sistent creep damage model for modelling secondary and tertiary creep behaviour is described.

2. Termodynanic formulation

In deriving thermodynamically consistent material models the first and second laws of thermodynamics constitute
the basic formalism that is followed here. The constitutive models can be derived from two potential functions, namely
the Helmholtz free energy and the dissipation potential, see Fremond (2002); Lemaitre and Chaboche (1990); Ottosen
and Ristinmaa (2005).

2.1. Energy balance

The first law of thermodynamics, i.e. the energy balance can be written in the form
d
(E + K) = Pmech + Pheat , (1)
dt
where E, K are the internal- and kinetic energies, which are defined as
Z Z
1
E= e dV, K = v v dV, (2)
V 2 V
where e is the specific internal energy, a state function depending on the specific entropy s and strain and v is the
velocity vector. The power of mechanical external forces and the power of non-mechanical sources, here assumed to
consist only of heat, are given as
Z Z Z Z
Pmech = b v dV + t v dS , Pheat = r dV q n dS , (3)
V S V S

where r is the internal heat source per unit mass and q heat-flux vector. After some manipulations, the energy balance
can be written as
Z Z
e dV = : gradv + r divq dV.


(4)
V V

By using the definition of infinitesimal strain = 12 [ gradu + ( gradu)T ], the energy balance (4) can be written in the
local form as

e = : + r divq. (5)

The specific internal energy e is a state function depending on the specific entropy s and strain . Utilising the
partial Legendre transformation, the specific Helmholtz free energy = e sT is obtained, which is a state function
depending of measurable state variables, absolute temperature T and strain . The local form of the energy balance
(5) can be transformed into the form

( + sT + sT ) = : + r divq. (6)

2.2. Entropy inequality

The second law of thermodynamics imposes restrictions on the process. The entropy inequality which is also known
as the Clausius-Duhem inequality is
r
Z Z Z
d qn
s dV dV dS , (7)
dt V V S
 P. Kauppila et al. / Structural Integrity Procedia 00 (2016) 000000 3

where s is the specific entropy. After some manipulations, the dissipation power of the system can be written as
= ( + sT ) + : T 1 gradT q, (8)
and the entropy inequality is simply
0. (9)
In the geometrically linear theory, the strain tensor can be additively divided into the elastic e , inelastic c and
thermal strain th components as
= e + c + th . (10)
The specific Helmholtz free energy depends on temperature T , integrity , which describes the degradation of the
material and thermoelastic strain te = c i.e. = (T, , te ). Therefore, the dissipation power can be expressed
as

! !
= s + + :te + :c T 1 gradT q. (11)
T te
The integrity and damage D variables are related by = 1 D. Integrity has the value 1 at the undamaged initial
state and the value 0 at the completely damaged state, whereas damage evolves from zero to one.
Dissipative mechanisms of the system are described by the complementary dissipation potential =
(Y, q, ; T, ), which is a monotonous function with respect to all of its arguments Y, q and , giving the defini-
tion for the dissipation power as

= q + : + Y. (12)
q Y
By defining the thermodynamic force Y = / dual to the integrity rate and equating the dissipation power (11)
with the definition (12), results in equation
gradT
! ! ! ! !
s + T + : te + c : + Y + q = 0. (13)
T te Y T q
Since this equation has to be fulfilled with all possible thermodynamically admissible processes T , te , , Y and q, the
following constitutive equations are obtained

s = , = , c = , = , T 1 gradT = . (14)
T te Y q
Substituting these general constitutive equations into the local form of the energy equation (6), the following form
is obtained
2 2
!
c T = divq + r + :c + te + Y , (15)
te T T
where the heat capacity c is defined as
2
c = T . (16)
T 2
Equation (15) is the thermomechanically coupled heat equation, where the heat input due to the mechanical pro-
cesses is described by the three last terms on the right-hand-side, which describe the heat input due to viscoplastic,
thermoelastic and damage processes.
If the complementary dissipation potential is convex with respect to all of its arguments, which are the thermody-
namic forces, the Clausius-Duhem inequality (CDI) (9) will be automatically satisfied. Convexity of the complemen-
tary dissipation potential is not necessary. A sufficient condition which guarantees the satisfaction of the CDI is that
the potential function is monotonous with respect to all of its arguments.
4 P. Kauppila et al. / Structural Integrity Procedia 00 (2016) 000000

3. The specific model

3.1. Potential functions

The reversible prosesses are described by the specific Helmholtz free energy = (T, te , ), which depends on
absolute temperature T , thermoelastic strains te , and integrity . In this study, the following form of the specific
Helmholtz free energy is chosen
!
T 1
= c T T ln + (te th ): Ce : (te th ), (17)
Tr 2
where th = (T T r ), is the thermal strain, Ce the elasticity tensor, a second order tensor containing the thermal
expansion coefficients and T r is an arbitrary reference temperature. For an isotropic solid the elasticity and the thermal
expansion tensor have the forms
E E
Ce = II+ I, and = I, (18)
(1 + )(1 2) 1+
where I and I are the second- and fourth order identity tensors, respectively. All the material coefficients, the Youngs
modulus E, the Poissons ratio and the linear coefficient of thermal expansion can depend on temperature.
The complementary dissipation potential can be additively decomposed into thermal, damaging and viscoplastic
parts as
(Y, q, ; T, ) = th (q; T ) + d (Y; T, ) + c (; T, ), (19)
where the thermal part is
1 1 1
th (q; T ) = T q q. (20)
2
For an isotropic solid the thermal conductivity tensor is simply = I, where the thermal conductivity can depend
on temperature. Damage affects also to the thermal conductivity, and thus the coefficient of thermal conductivity can
also depend on the integrity. Here this effect is neglected.
For creep the following Norton type potential function is adopted
! p+1
hc (T ) rc
c (; T, ) = , (21)
p + 1 tc rc

where = 3J2 is the von Mises effective stress (J2 is the second invariant of the deviatoric stress), tc is a charac-
teristic time for creep and it is directly related to the relaxation time, hc Arrhenius-type thermal activation function
hc (T ) = exp(Qc /RT ), where Qc is the creep activation energy and R is the universal gas constant. Choice of the
reference stress rc (also known as a drag stress) is discussed later.

3.2. Monkman-Grant hypothesis

From the creep tests of many metals and alloys it is concluded that the product of minimum creep strain-rate and
rupture time is a constant (Riedel, 1987; Nabarro and Villers, 1995), i.e.
cmin trup = constant. (22)
This constant is almost independent of temperature and stress and it is known as the Monkman-Grant parameter. It is
often expressend in a modified form
CMG = (cmin )m trup , (23)
where m is in the range 0.8 1. The value of the Monkman-Grant parameter equals roughly to the rupture strain
CMG / rup . (24)
 P. Kauppila et al. / Structural Integrity Procedia 00 (2016) 000000 5

3.3. Damage potentials

In this study two variants for the damage potential are compared. Both of them are of Kachanov-Rabotnov type
(Kachanov, 1958, 1986):
!r+1
hd (T ) Yr Y
d (Y; T, ) = , (model 1) (25)
r + 1 td k Yr
! 1 p+1
hc (T ) Yr Y 2
d (Y; T, ) = 1 , (model 2) (26)
( 2 p + 1)(1 + k + p) td k Yr
where td is a characteristic time for damage evolution, hd is an Arrhenius-type thermal activation function for damage
processes hd (T ) = exp(Qd /RT ), where Qd is the damage activation energy and R is the universal gas constant as in
the creep activation function. The reference value of the thermodynamic force, Yr , is chosen by Yr = rd 2 /(2E), where
rd is a reference stress for the damage process.
The first version (25) results in a more general model, whereas the second one (26) is restricted to satisfy the
Monkman-Grant hypothesis (with m = 1) exactly and therefore its parameters are coupled to the creep model. For
these two models the Monkman-Grant parameter have the values
! p2r
1 td hc td
CMG = min trup =
c
(model 1) and CMG = (model 2). (27)
1 + k + 2r tc hd r tc
To account for different damage evolution in tensile and compressive regions the damage potentials could be
changed to the form
d (Y; T, , e ) (Y/Yr + tr e )r+1 , (28)
where is an extra material parameter. Due to the lack of material data available to the authors this form has not been
used.

3.4. Constitutive equations

From the general constitutive equations (14) and with the specific choices (17), (20) and (21) the following consti-
tutive equations are obtained
!p
hc
= Ce : e , c = , q = gradT. (29)
tc rc
The integrity rates resulting from the considered two models (25) and (26) are
!r !1 p
hd Y hc Y 2
= , (model 1) and = (model 2). (30)
td k Yr (1 + k + p)d k Yr

4. Response in uniaxial creep test

The material parameters can be determined from the uniaxial creep tests at different temperatures and under several
constant stress values. In this simplified loading case, the constitutive equations can be integrated in a closed form,
the resulting integrity development and the creep rupture time for the model (25) are
!2r 1/(1+k+2r) !2r
t 1
= 1 (1 + k + 2r)hd , trup = td .

(31)
r (1 + k + 2r)hd r

td
The yield stress is selected as the reference value in the creep and damage evolution equations, i.e. rd = rc =
y0 (T ) r . For the total strain, the following formula is obtained
! p2r ! p #(1+k+2rp)/(1+k+2r)
hc
"
1 t
= + 1 1 (1 + k + 2r)hd ,

(32)

E 1 + k + 2r p hd r r td





6 P. Kauppila et al. / Structural Integrity Procedia 00 (2016) 000000

Table 1. The calibrated model parameters for the 7CrMoVTiB10-10 steel (T24), qc = Qc /R and qd = Qd /R.

model tc [s] pr td [s] a qc [K] rr qd [K] b

1 3039.9 14.77 37.768 4.804 7137.6 7.545 9350.1 5.201
2 3414.1 14.59 41.26 4.891 7137.6 - - -

400 400
300 300

200 200
[MPa]

[MPa]
100 100

50 50

102 101 100 101 102 103 104 105

c,min [%/103 h] trup [h]

Fig. 1. The calibrated results and the data from Arndt et al. (2000) for the T24-steel based on the minimum creep strain-rate (lhs) and on the creep
strengths (rhs). Results for the model 1 are shown by solid lines and for the model 2 by dashed lines. The sets from top to bottom correspond to
tests at temperatures 500 C, 550 C and 600 C, respectively.

which gives the following expression for the rupture strain

! p2r
1 td hc
rup = + . (33)
E 1 + k + 2r p tc hd r
For the second model (26) the integrity evolution, the creep rupture time and the rupture strain in a constant stress
creep test are
! p #1/(1+k+p) !p
td 1 + k + p td
"
t
= 1 hc , trup = , rup = + . (34)
r td hc r E 1+k tc

5. Determination of the material parameters

At temperatures between 500 C and 600 C the material parameters for the 7CrMoVTiB10-10 steel (T24) have been
determined from the material data of the manufacturer (Arndt et al., 2000). Temperature dependency of the parameters
p and r as well as for the yield stress y0 are assumed to be reasonably well presented by the linear expressions
p(T ) = pr [1 + a(T T r )/T r ] , r(T ) = rr [1 + b(T T r )/T r ] , r = y0 (T ) = cT, (35)
where the values = 1123 MPa and c = 1 MPa/K give a good fit to the data.
In the calibration it is also assumed that k + 2r = p + 2 (Lemaitre, 1992, Chapter 3.3.3), which gives k = 2, if
p = 2r. The model parameters are presented in Table 1 and the model predictions and the material data from Arndt
et al. (2000) are shown in Fig. 1.
It can be seen from the results presented in Table 1 that the material parameters are realistic and the Monkman-
Grant hypothesis is relatively well satisfied also for the model 1 at the reference temperature T r = 773 K. Furthermore
it is seen that the values for the powers p and r are decreasing with increasing temperature, which is typical for most
metals, cf. (Garofalo, 1965, Table 3.1). 1
1
 The displacement is considered to be a result of thermal expansion caused by variable operating temperature and
local temperature differences in a large superheater structure. The analyzed lifetime of the header is 150 creep-
fatigue cycles in which the duration of load changes is one hour and the duration of hold periods at constant load
is 200 hours. P. Kauppila et al. / Structural Integrity Procedia 00 (2016) 000000 7

Displacement
Temperature

Temperature (C)
600

Displacement
0 500
Time

(a) (b) (c)

Fig. 2. The FE mesh, mainly 20 node hexahedral ANSYS SOLID186 elements & some 10 node tetrahedal SOLID187 elements, and the prescribed
displacement history at the end of the tube nozzle. The displacement at the header end and the creep temperature changes periodically and the
Fig. 2.consist
periods The FEM-model mesh
of ramp and holdintimes.
ANSYS (a) and the time-varying displacement at the end of the tube nozzle (b). The displacement at the end of
P. Kauppila
the tube nozzle and the creep temperature changes et al. /and
periodically Structural Integrity
the periods consistProcedia 00hold
of ramp and (2016) 000000
times (c).

0,40 0,010
Model 1 Model 1

Eq. creep strain (-)

0,35
The main results of the analyses, the values of the damage 0,009
Model 2parameter and the equivalent creep strains Modelat2 the
0,008
Damage (-)

most critical location are shown in Fig. 3.0,30 The most critical location in these analysis was on the surface of the weld
between the header pipe and the tube 0,25 nozzle in all situations, although it shall be noted 0,007that these models do not
take the special characteristics of a welded joint into account. According to the results, 0,006the values of the damage
0,20 0,005
parameter calculated using the model 1 are slightly higher than the values calculated using the model 2. In any
0,15 0,004
case, the values of the damage parameter calculated using the model 1 are only 37 pp greater than the values
calculated using the model 2, and thus 0,10the results of the both models are relatively 0,003 precise for practical
applications. In practical applications the analyzed 0 0,2
material 0,4
can be0,6 0,8 to be damaged,
considered 0 when 0,2 the0,4 value of 0,6 0,
the damage parameter is over 0.3 and the onset Displacement
of the tertiary(mm)creep phase has occurred. Displacement
Despite the slight (mm)
difference in the values of the damage parameter between the models, the equivalent creep strains at the most
damaged
Fig. 3. Damagelocation are
nearalmost
(a)
distribution the mostequal
criticalbetween theheader.
location of the models The1accumulated
and
(b) 2. damage and the equivalent creep strain at the most critical
(c)
location as functions of the prescribed displacement.

Fig. 3. The accumulated damage of the header at the most critical location (a). The accumulated damage at the most critical locati
6.the accumulated equivalent creep strain at the most damaged location (c) as functions of the length of the displacement in the analys
Implementation

The described model has been implemented in a structural FE-code ANSYS as a user-defined material subroutine
From aIntegration
USERMAT. practicalofpoint ofdependent
the rate view, themodelbothisdeveloped
performed by models
using theyield relatively
implicit backward equal
Eulerresults
method,within
which the tem
range
has 500600
proven C. However,
to be accurate especially the model
for large 2 is accurate
practically usable timeonlystep
in sizes
relatively highrate-dependent
in solving creep temperatures
problems and it b
(Kouhia et al., 2005), although it does not completely inherit such a nice property when combined with damage, see
slightly inaccurate as the temperature lowers, which is a result of the assumed Monkman-Grant relations
e.g. Wallin and Ristinmaa (2001).
model 1 aistypical
A part of moresteam
accurate
boiler also in lowheader
superheater creep( temperatures
355.6 mm 36 mm mainly because
with two symmetryof its greater
planes) madeflexibility
of T24 due t
number of calibration parameters.
steel and operating within a temperature range from 500 to 600
C has been analysed under creep-fatigue loading.
The loading consist of a static internal pressure of 14.0 MPa(g), time-varying temperature and displacement at the end
ofAcknowledgements.
the tube nozzle 44.5 mm This work
6.3 mm. was
Thecarried out isinconsidered
displacement the researchto be aprogram Flexible
result of thermal Energycaused
expansion Systems (FLE
by variable operating temperature and local temperature differences in a large superheater structure. The Youngs
supported by Tekes and the Finnish Funding Agency for Innovation. The aim of FLEXe is to creat
modulus has values 175 MPa at 500 C, 168 MPa at 550 C and 163 MPa at 600 C, and it is linearly interpolated
technological
between and(Arndt
these values business
et al., concepts
2000). Theenhancing theis radical
Poissons ratio assumedtransition from the
to be independent current energy
of temparature and the systems
sustainable
value = 0.3 is systems. FLEXe consortium
used. The analysed lifetime of theconsists
header is of15017 industrialcycles
creep-fatigue partners andthe10
in which research
duration organisatio
of load
programme is coordinated by CLIC Innovation Ltd. www.clicinnovation.fi
changes is one hour and the duration of hold periods at constant load is 200 hours.
The main results of the analyses, damage field (D = 1 ) and the values of the damage and the equivalent
creep strains at the most critical location are shown in Fig. 3. The most critical location in these analysis is on the
References

Altenbach, H., Gorash, Y., Naumenko, K., 2009. Steady-state creep of a pressurized thick cylinder in both the linear and the power law
Acta Mechanica 195, 263274.
Arndt, J., Haarmann, K., Kottmann, G., Vaillant, J., Bendick, W., Kubla, G., Arbab, A., Deshayes, F., 2000. The T23/T24 Book. 2nd ed., Va
8 P. Kauppila et al. / Structural Integrity Procedia 00 (2016) 000000

surface of the weld between the header pipe and the tube nozzle in all situations, although it should be noted that these
models do not take the special characteristics of a welded joint into account. According to the results, the values of
the damage parameter calculated using the model 1 are slightly higher than the values calculated using the model 2. In
any case, the values of the damage parameters calculated using the model 1 are only 3-7 percentage points greater than
the values calculated using the model 2, and thus the results of the both models are relatively accurate for practical
applications. A material can be considered to be essentially fully damaged when the value of the damage parameter
exceeds the value 0.3 and the onset of tertiary creep phase has occurred. Despite the slight difference in the values
of the damage parameter for the models, the equivalent creep strains at the most damaged location are almost equal
between models 1 and 2.
From an engineering point of view, the both developed models yield practically equal results within the temperature
range 500-600 C. However, the model 2 is accurate only in relatively high creep temperatures and it becomes slightly
inaccurate at lower temperatures, which is a result of the assumed Monkman-Grant relationship. The model 1 is more
accurate also in low creep temperatures mainly because of its greater flexibility due to higher number of calibration
parameters.

Acknowledgements. This work was carried out in the research program Flexible Energy Systems (FLEXe) and sup-
ported by Tekes and the Finnish Funding Agency for Innovation. The aim of FLEXe is to create novel technological
and business concepts enhancing the radical transition from the current energy systems towards sustainable systems.
FLEXe consortium consists of 17 industrial partners and 10 research organisations. The programme is coordinated by
CLIC Innovation Ltd. www.clicinnovation.fi

References

Altenbach, H., Gorash, Y., Naumenko, K., 2009. Steady-state creep of a pressurized thick cylinder in both the linear and the power law ranges.
Acta Mechanica 195, 263274.
Arndt, J., Haarmann, K., Kottmann, G., Vaillant, J., Bendick, W., Kubla, G., Arbab, A., Deshayes, F., 2000. The T23/T24 Book. 2nd ed., Vallourec
& Mannesmann Tubes.
Betten, J., 2005. Creep Mechanics. Springer-Verlag, Berlin.
da Costa Andrade, E., 1910. On the viscous flow in metals, and allied phenomena. Proceedings of the Royal Society A 84, 112.
Fremond, M., 2002. Non-Smooth Thermomechanics. Springer, Berlin.
Garofalo, F., 1965. Fundamentals of Creep and Creep-Rupture in Metals. Macmillan series in Materials Science, Macmillan, New York.
Gorash, Y., 2008. Development of a creep-damage model for non-isothermal long-term strength analysis of high-temperature components operating
in a wide stress range. Ph.D. thesis. Martin-Luther-Universitat. Halle-Wittenberg, Germany.
Hayhurst, D., 1972. Creep rupture under multiaxial states of stress. Journal of Mechanics and Physics of Solids 20, 381390.
Hayhurst, D., 1994. The use of continuum damage mechanics in creep analysis for design. Journal of Strain Analysis 25, 233241.
Kachanov, L., 1958. On the creep fracture time. Iz. An SSSR Ofd. Techn. Nauk. , 2631(in Russian).
Kachanov, L., 1986. Introduction to continuum damage mechanics. volume 10 of Mechanics of Elastic Stability. Martinus Nijhoff Publishers.
Kouhia, R., Marjamaki, P., Kivilahti, J., 2005. On the implicit integration of inelastic constitutive equations. International Journal for Numerical
Methods in Engineering 62, 18321856.
Lemaitre, J., 1992. A Course on Damage Mechanics. Springer-Verlag, Berlin.
Lemaitre, J., Chaboche, J.L., 1990. Mechanics of Solid Materials. Cambridge University Press.
Murakami, S., 2012. Continuum Damage Mechanics. volume 185 of Solid Mechanics and Its Applications. Springer Netherlands.
Nabarro, F., Villers, H., 1995. The Physics of Creep and Creep Resistant Alloys. Taylor & Francis Ltd.
Naumenko, K., 2006. Modeling of high temperature creep for structural analysis applications. Ph.D. thesis. Martin-Luther-Universitat. Halle-
Wittenberg, Germany.
Ottosen, N., Ristinmaa, M., 2005. The Mechanics of Constitutive Modeling. Elsevier.
Riedel, H., 1987. Fracture at High Temperatures. MRE Materials Reserch and Engineering, Springer-Verlag, Berlin, Heidelberg.
Wallin, M., Ristinmaa, M., 2001. Accurate stress updating algorithm based on constant strain rate assumption. Computer Methods in Applied
Mechanics and Engineering 190, 55835601.