Creep and Viscoplasticity
About Creep
In the literature, the terms viscoplasticity and creep are often used interchangeably to refer to the class of problems related to rate-dependent plasticity.
Creep is an inelastic time-dependent deformation that occurs when a material is subjected to stress (typically much less than the yield stress) at sufficiently high temperatures.
The creep strain rate, in a general case, depends on stress, temperature, and time, usually in a nonlinear manner:
It is often possible to separate these effects as shown in this equation:
Experimental data shows three types of behavior for the creep strain rate at constant stress as function of time. Researchers normally subdivide the creep curve into three regimes, based on the fact that many different materials show similar responses:
In the initial primary creep regime (also called transient creep) the creep strain rate decreases with time to a minimum steady-state value.
In the secondary creep regime the creep strain rate is almost constant. This is also called steady-state creep.
In the tertiary creep regime the creep strain increases with time until a failure occurs.
When this distinction is assumed, the total creep rate can be additively split into primary, secondary, and tertiary creep rates
In most cases, Fcr1 and Fcr3 depend on stress, temperature and time, while secondary creep, Fcr2, depends only on stress level and temperature. Normally, secondary creep is the dominant process. Tertiary creep is seldom important because it only accounts for a small fraction of the total lifetime of a structure.
Figure 3-16: Uniaxial creep as a function of logarithmic time.
Creep and Viscoplastic Material Models
Despite the fact that the creep response of a given material is related to its atomic structure, a macroscopic (continuum mechanics) description is normally appropriate for modeling scientific and engineering problems.
In COMSOL Multiphysics, there are several creep models. These models can be split into two main groups. One set of models are more general, and you will have to express the creep rate yourself, based on other variables as a stress tensor and temperature. These models are:
In addition to the basic models for creep described above, there are also predefined material models for creep in metals and crystalline solids:
All creep models are contributing subnodes to a basic material model like Linear Elastic Material and they can be combined with any other subnodes, such as Plasticity or Thermal Expansion to create more advanced models. They can also be combined with each other to model several creep mechanisms acting at the same time.
Creep Potential
Some authors use a creep potential to describe the secondary creep rate, so that the creep rate is written in a way similar to the flow rule for plasticity:
and
Here, Qcr is a user defined creep potential, which is normally written in terms of invariants of the stress tensor.
Volumetric creep is obtained when the creep potential depends only on the first invariant of Cauchy stress tensor, I1(σ), since
This is equivalent to that the creep potential would depend on the pressure p = −I1/3.
When the creep potential depends only on the second deviatoric invariant of Cauchy stress tensor, J2(σ), the deviatoric creep model is obtained since
This is equivalent to that the creep potential would depend on the effective stress σe = √3J2.
When (in SI units) the creep potential, Qcr, is given in units of Pa, the rate multiplier η is given in units of 1/s.
Volumetric Creep
The creep strain rate is calculated by solving the rate equation
so that the creep rate tensor is a diagonal tensor. The trace of the creep strain rate tensor, the volumetric creep strain rate, equals the user input
The volumetric creep strain rate usually depends on the first invariant of Cauchy stress I1(σ) or the pressure p = −I1/3, in addition to the temperature and other material parameters.
Volumetric creep is not generally used to model creep in metals, but it is commonly used to model creep in soils or other geological materials.
Deviatoric Creep
The creep strain rate is calculated by solving the rate equation
Here, nD is a deviatoric tensor coaxial to the stress tensor.
The effective creep strain rate , , normally depends on the second deviatoric invariant of the stress J2(σ) or the effective or von Mises (effective) stress σe, in addition to the temperature and other material parameters.
The deviatoric tensor nD is defined as
(3-44)
The resulting creep strain rate tensor is also deviatoric, since trace(nD) = 0
Given the property
the effective creep strain rate equals the absolute value of the user input
Deviatoric creep is very popular to model creep in metals and alloys. For example, Norton’s law is a deviatoric creep model.
User Defined Creep
The creep strain tensor is calculated by time-integration the user defined symmetric creep strain rate tensor .
Norton Law (Power law)
The most common model for secondary creep is the Norton equation where the creep strain rate is proportional to a power of the effective stress, σe:
This is normally true at intermediate to high stress levels and at absolute temperatures of T/Tm > 0.5, where Tm is the melting temperature (that is, the temperature in the solid is at least as high as half the melting temperature Tm). An “Arrhenius type” temperature dependency can also be included. It is defined by
where Q is the activation energy (SI unit: J/mol), R is the gas constant, and T is the absolute temperature (SI unit: K).
Norton creep is a deviatoric temperature-dependent creep model, with a creep rate equation written as
(3-45)
Here, A is the creep rate coefficient (SI unit: 1/s), n is the stress exponent (dimensionless), σref a reference stress level (SI unit: Pa), and nD is a deviatoric tensor coaxial to the stress tensor as defined in Equation 3-44.
See also the description of the Norton material model in the Solid Mechanics interface documentation.
Norton-Bailey Law
A common model for modeling primary and secondary creep together is the so-called Norton-Bailey (or Bailey-Norton) model. Here, the creep strain is proportional to a power of time and to a power of the effective stress
which for the creep strain rate becomes a time hardening formulation of Norton’s law. Differentiating with respect to time will give the rate form.
Norton-Bailey creep is a deviatoric temperature-dependent creep model, furbished with either a time-hardening or a strain-hardening primary creep model. The creep rate equation for the time-hardening model used in COMSOL Multiphysics is written as
(3-46)
where nD is a deviatoric tensor coaxial to the stress tensor as defined in Equation 3-44, and Fcr is expressed as in the Norton model:
(3-47)
Here, A is the creep rate coefficient (SI unit: 1/s), n is the stress exponent (dimensionless), σref is a reference stress level (SI unit: Pa), tref and tshift are the reference and shift times (SI unit: s), and m is the time-hardening exponent (dimensionless).
The strain-hardening variant of this creep law is implemented as
(3-48)
where εcr,e is the effective creep strain, and εshift is the effective creep strain shift.
The time and frequency shifts in Equation 3-46 and Equation 3-48 serve two purposes:
The strain rate expressions actually predicts an infinite creep rate at t = 0, unless a shift is used. This singularity is weak in the sense that the time integral is well defined, but it will cause problems for the numerical solution. You can then add a small shift to overcome this problem.
See also the description of the Norton-Bailey material model in the Solid Mechanics interface documentation.
Garofalo Law (Hyperbolic Sine Law)
At very high stress levels, the creep rate is proportional to the exponential of the effective stress
Garofalo showed (Ref. 8, Ref. 9) that the power-law and exponential creep are limiting cases for the general empirical expression
This equation reduces to a power-law (Norton law) for ασe < 0.8 and approaches exponential creep for ασe > 1.2, where 1 is a reference effective stress level.
Garofalo creep is also a deviatoric creep model with a creep rate proportional to the hyperbolic sine function. It can also be augmented by an “Arrhenius type” temperature dependency such that
where Q is the activation energy (SI unit: J/mol), R is the gas constant, and T is the absolute temperature (SI unit: K). The complete creep rate equation as used in COMSOL Multiphysics then reads
where, A is the creep rate (SI unit: 1/s), n is the stress exponent (dimensionless), and σref a reference effective stress level (SI unit: Pa). nD is a deviatoric tensor coaxial to the stress tensor as defined in Equation 3-44.
See also the description of the Garofalo (hyperbolic sine) material model in the Solid Mechanics interface documentation.
NaBarro-Herring Creep (Diffusional Creep)
At low stress levels and high temperatures, Nabarro and Herring (Ref. 6, Ref. 7) independently derived an expression for the creep rate as a function of atomic diffusion
Here, d is the grain diameter, Dv is the volume diffusivity through the grain interior, b is Burgers vector, kB is the Boltzmann’s constant, and T is the absolute temperature. nD is a deviatoric tensor coaxial to the stress tensor as defined in Equation 3-44.
See also the description of the Nabarro-Herring material model in the Solid Mechanics interface documentation.
Coble Creep (Diffusional Creep)
Coble creep (Ref. 6, Ref. 7) is closely related to Nabarro-Herring creep but takes into account the ionic diffusivity along grain boundaries Dgb
See also the description of the Coble material model in the physics interface documentation.
Weertman Creep (Dislocation Creep)
At intermediate to high stress levels and temperatures T/Tm > 0.5, the creep mechanism is assumed to be diffusion-controlled movements of dislocations in the crystal lattices (Ref. 7)
where and nD is a deviatoric tensor coaxial to the stress tensor as defined in Equation 3-44. Generally, the stress exponent n takes values between 3 and 5.
A general relation between creep rate and several material parameters is the Mukherjee-Bird-Dorn equation (Ref. 6)
Here, T is the temperature, d is the grain size, b is the Burgers vector, D is the self diffusion coefficient, G is the shear modulus, and eQ/RT is an “Arrhenius” type of temperature dependency.
For high temperatures, Mukherjee-Bird-Dorn equation describes Weertman creep when setting p = 0. Setting n = 0 and p = 2 describes Nabarro-Herring, and setting n = 0 and p = 3 describes Coble creep. Harper-Dorn creep is obtained by setting n = 1 and p = 0.
See also the description of the Weertman material model in the Solid Mechanics interface documentation.
Anand Viscoplastic Model
The Anand viscoplasticity (Ref. 9) is a deviatoric creep model suitable for large, isotropic, viscoplastic deformations in combination with small elastic deformations.
The viscoplastic strain rate equation reads
where nD is a deviatoric tensor coaxial to the stress tensor as defined in Equation 3-44, and the creep rate is calculated from
Here, A is the creep rate coefficient (SI unit: s1), Q is the activation energy (SI unit: J/mol), m is the stress sensitivity, ξ is the multiplier of stress, R is the gas constant, and T is the absolute temperature (SI unit: K).
The internal variable, sa, is called deformation resistance (SI unit: Pa) and is calculated from the rate equation
with the initial condition sa(0) = sinit. Here, h0 is the hardening constant (SI unit: Pa), and a is the hardening sensitivity.
The variable sa* is the saturation value of the deformation resistance sa, which is calculated from the expression
where s0 is the deformation resistance saturation coefficient (SI unit: Pa), and n is the deformation resistance sensitivity.
Chaboche Viscoplastic Model
The viscoplastic strain rate tensor is given by
Here, A is the viscoplastic rate coefficient (SI unit: 1/s), n is the stress exponent (dimensionless), σref a reference stress level (SI unit: Pa), and nD is a deviatoric tensor coaxial to the stress tensor. The Macaulay brackets are applied on the yield function, which is defined as done for plasticity
The effective stress φ(σ) is either the von Mises, Tresca, or Hill stress; or a user defined expression, and σys is the yield stress (which may include a linear or nonlinear Isotropic Hardening model). The stress tensor used in the effective stress φ(σ) is shifted by what is usually called the back stress, σback when Kinematic Hardening is included.
The deviatoric tensor nD is computed from the plastic potential Qp
When von Misses effective stress is used, the associated flow rule reads Qp = Fy, and the deviatoric tensor nD is defined as done for deviatoric creep
(3-49)
Given the property
the effective viscoplastic strain rate is equivalent to
Perzyna Viscoplastic Model
Perzyna viscoplastic model is similar to the Chaboche model, with the exception than the stress exponent is set equal to one. The viscoplastic strain rate tensor is then given by
Some authors denote the viscosity as the quotient η = σref/A.
See also the description of Viscoplasticity in the physics interface documentation.
Energy Dissipation
Since creep and viscoplasticity are inelastic processes, the dissipated energy density can be calculated by integrating the creep dissipation rate density (SI unit: W/m3) given by
In case many creep sub-nodes are added, the creep dissipation rate density is calculated from the cumulative creep strain rate tensor .
The total energy dissipated in a given volume can be calculated by a volume integration of the dissipated creep energy density Wc (SI unit: J/m3).
When the Calculate dissipated energy check box is selected, the dissipation rate density due to creep is available under the variable solid.Wcdr and the dissipation rate density due to viscoplasticity is available under the variable solid.Wvpdr. The dissipated energy density due to creep is available under the variable solid.Wc and due to viscoplasticity under the variable solid.Wvp. Here solid denotes the name of the physics interface node.
Converting Between Different Creep Data Representations
The equation forms described in for the different creep models above differ from the forms most commonly found in the literature. The difference lies in the introduction of normalizing reference values such as the reference stress σref and reference time tref. These values are in a sense superfluous, and can in principle be chosen arbitrarily. The choice of reference values will however affect the numerical values to be entered for the material data. This system has two advantages
Norton Law
Material data for a Norton law is often available as the parameters AN and n in the equation
The coefficient AN has a physical dimension which depends on the value of n and the unit have an implicit dependence on the stress and time units. Converting the data to the form used in COMSOL Multiphysics (Equation 3-45) requires the introduction of the reference stress σref. It is here convenient to use the implicit stress unit for which AN is given as reference stress. The creep rate coefficient A will then have the same numerical value as AN, and you do not need to do any conversions.
The physical dimension of A is however (time)-1, whereas the physical dimension of AN is (stress)-n(time)-1.
Another popular way of representing creep data is to supply the stress giving a certain creep rate. As an example, σc7 is the stress at which the creep strain rate is 10-7/h. Data on this form is also easy to enter: You set the reference stress to the value of σc7, and enter the creep rate coefficient as 1e-7[1/h].
Example
 
Assume that a carbon steel has the following two equivalent descriptions of its creep properties at a certain temperature:
σc7 = 70 MPa, and stress exponent n = 4.5.
AN  = 4.98·10-16 with respect to units MPa and hours, and stress exponent n = 4.5.
In the first case, enter:
σref as 70[MPa]
n as 4.5
A as 1e-7[1/h].
In the second case, enter:
σref as 1[MPa]
n as 4.5
A as 4.98e-16[1/h]
These two sets of data describe the same material.
Norton-Bailey law
Material data for a Norton-Bailey law usually is usually written in terms of the creep strain, rather than the creep strain rate, so that the form of the constitutive relation is
In this case, the coefficient ANB has implicit dimension and units which depend on the values of n and m, and on the stress and time units.
Converting the data to the form used in COMSOL Multiphysics (Equation 3-46 and Equation 3-47) requires the introduction of an both an arbitrary reference stress σref and an arbitrary reference time tref. If you use the implicit units for which ANB is given as the reference values, then the constant A will have the same numerical value as ANB.
Garofalo Law
Since the stress inside in the Garofalo law appears as an argument to a sinh() function, it must necessarily be nondimensionalized. Most commonly this is however written as
Comparing with the expression in COMSOL Multiphysics,
it is evident that the reference stress should be chosen as
In this case, there is no arbitrariness in the choice of σref, since α is an actual material parameter.