Elastoplastic Material Models
In this section:
Introduction
Many materials have a distinct elastic regime, in which the deformations are recoverable and path independent. When the stresses exceed a certain level, the yield limit, permanent plastic strains will appear.
The elastic part of the constitutive relation can be described by either a Linear Elastic Material, Nonlinear Elastic Material Models, or by Hyperelastic Material Models.
Elastoplastic material models are common, both when modeling metals and soils.
In geotechnical applications it is common to define compressive stresses as having positive signs. In COMSOL Multiphysics, the convention is however to always use positive signs for tensile stresses.
Defining the Yield Surface
A yield criterion serves to define the stress condition under which plastic deformation occurs. Stress paths inside the yield surface result in purely recoverable deformations (elastic behavior), while paths intersecting the yield surface produces both recoverable and permanent deformations (plastic strains).
In general, the yield surface can be described as
where fc can be a constant value (for perfectly plastic materials), or a variable for strain-hardening materials. The yield surface F is a surface in the space of principal stresses, in which the elastic regime (F ≤ 0) is enclosed.
For brittle materials, the yield surface represents a failure surface, which is a stress level at which the material collapses instead of deforms plastically.
Some authors define the yield criterion as f (σ) =  fc, while the yield surface is an isosurface in the space of principal stresses F = 0, which can be chosen for numerical purposes as .
Isotropic Plasticity
For isotropic plasticity, the plastic potential Qp is written in terms of at most three invariants of Cauchy’s stress tensor
where the invariants of the stress tensor are
so that the increment of the plastic strain tensor can be decomposed into
The increment in the plastic strain tensor includes in a general case both deviatoric and volumetric parts. The tensor is symmetric given the following properties
(3-21)
A common measure of inelastic deformation is the effective plastic strain rate, which is defined as
(3-22)
The trace of the incremental plastic strain tensor, which is called the volumetric plastic strain rate , is only a result of dependence of the plastic potential on the first invariant I1(σ), sinceJ2/∂σ and J3/∂σ are deviatoric tensors
For metal plasticity under the von Mises or Tresca criteria, the volumetric plastic strain rate is always zero because the plastic potential is independent of the invariant I1(σ). This is known as J2 plasticity.
Incompressible plastic deformation is experimentally observed in metals, but it not the case for most materials used in geotechnical applications. For instance, a nonzero volumetric plastic strain is explicitly used in the Cam-Clay material.
Yield Function
When an associated flow rule is applied, the yield function must be smooth, that is, continuously differentiable with respect to the stress. In COMSOL Multiphysics, the following form is used:
where σys is the yield stress. The scalar function φ(σ) is called effective stress. The default form of the effective stress is the von Mises stress, which is often used in metal plasticity:
Other expressions can be defined, such as Tresca stress, Hill orthotropic plasticity, or a user-defined expression.
The Tresca effective stress is calculated from the difference between the largest and the smallest principal stress
A user-defined yield function can by expressed in terms of invariants of the stress tensor such as the pressure (volumetric stress)
the effective (von Mises) stress σmises, or other invariants, principal stresses, or stress tensor components.
von Mises Criterion
The von Mises criterion suggests that the yielding of the material begins when the second deviatoric stress invariant J2 reaches a critical value. This criterion can be written in terms of the elements of Cauchy’s stress tensor (Ref. 1)
or equivalently .
The von Mises criterion is implemented as
where σys is the yield stress level (yield stress in uniaxial tension).
The effective or von Mises stress () is available in the variable solid.mises, where solid is the name of the physics interface node.
Tresca Criterion
The Tresca yield surface is normally expressed in terms of the principal stress components
The Tresca criterion is a hexagonal prism with its axis equally inclined to the three principal stress axes. When the principal stresses fulfill σ1 ≥ σ2 ≥ σ3, this criterion is written as
By using the representation of principal stresses in term of the invariants J2 and the Lode angle 0 ≤ θ ≤ π/3, this criterion can alternatively be written as
or equivalently
The maximum shear stress is reached at the meridians (θ = 0 or θ = π/3). The Tresca criterion can be circumscribed by setting the Lode angle θ = 0, or equivalently, by a von Mises criterion
The minimum shear is reached at θ = π/6, so the Tresca criterion can be inscribed by setting a von Mises criterion
When dealing with soils, the parameter k is also called undrained shear strength.
Tresca criterion can be used with either an associated or non-associated flow rule, in which case von Mises stress is applied in order to get better numerical performance.
The Tresca effective stress, σtresca = σ1 − σ3, is implemented in the variable solid.tresca, where solid is the name of the physics interface node.
Figure 3-4: The upper and lower limits of the Tresca criterion.
Figure 3-5: Classical yield criteria for metals. Tresca criterion (left) and von Mises criterion (right).
The von Mises and Tresca criteria are independent of the first stress invariant I1 and are mainly used for the analysis of plastic deformation in metals and ductile materials, though some researchers also use these criteria for describing fully saturated cohesive soils under undrained conditions. The von Mises and Tresca criteria belongs to what researchers call volume preserving or J2 plasticity, as the plastic flow is independent on the mean pressure.
Porous Plasticity
The description of plastic deformation in soils, porous metals and aggregates has a main difference with respect to traditional metal plasticity; as the yield function and plastic potential are not only defined in terms of the deviatoric stress tensor component (or the deviatoric stress invariant J2), but they also include dependencies on the hydrostatic pressure.
Shima-Oyane Criterion
Shima and Oyane (Ref. 14) proposed a plastic yield surface for modeling the compaction of porous metallic structures fabricated by sintering. This criterion can been applied for powder compaction at both low and high temperatures. The yield function and associated plastic potential is defined by an ellipsoid in the stress space. The plastic potential Qp is written in terms of both von Mises effective stress and mean pressure, and it also considers isotropic hardening due to changes in porosity. The plastic potential is defined by
here, σe is the effective stress, σ0 is the yield stress, pm is the pressure, and ρrel is the relative density, which is related to the porosity φ by ρrel = 1 − φ. The material parameters α, γ, and m are obtained from curve fitting experimental data. Typical material parameter values for copper powder are α = 6.2, γ = 1.03, and m = 5.
Gurson Criterion
Gurson criterion (Ref. 15) consists in a pressure dependent yield function to describe the constitutive response of porous metals, this yield function is derived from the analytical expression of an isolated void immersed in a continuum medium. The void volume fraction, or porosity φ is chosen as main variable. The yield function and associated plastic potential is not an ellipse in the stress space, as in Shima-Oyane Criterion, but it is defined in terms of the hyperbolic cosine function. The plastic potential for Gurson criterion reads
here, σe is the effective stress, σ0 is the initial yield stress, pm is the pressure, and φ is the porosity.
Gurson-Tvergaard-Needleman Criterion
Tvergaard and Needleman modified Gurson Criterion for porous plasticity to include parameters to better fit experimental data (Ref. 16-17). The resulting criterion is called in the literature Gurson-Tvergaard-Needleman (GTN) criterion. The plastic potential for GTN criterion reads
here, σe is the effective stress, σ0 is the initial yield stress, pm is the pressure, and φe is the effective void volume fraction (effective porosity). Typical correction parameter values are q1 = 1.5, q2 = 1.03, and q3 = q12.
The effective void value fraction (or effective porosity) φe used in the plastic potential is a function of the current porosity φ and other material parameters:
here φc is the critical void volume fraction (critical porosity) at which void coalescence begins, and φf is the void volume fraction at failure. When the porosity increases up to value of failure, the effective porosity takes a maximum value of φm; at this high porosity, the porous material loses the capacity to carry stresses. This maximum porosity value is derived from other parameters
since typically the correction parameters are related as q3 = q12, then φm = 1/q1.
Fleck-Kuhn-McMeeking Criterion
The Fleck-Kuhn-McMeeking criterion (Ref. 19), also called FKM criterion, was developed to model the plastic yielding of metal aggregates of high porosity. The yield function and associated plastic potential is derived from expressions for randomly distributed particles. The criterion is considered relevant for aggregates with porosity between 10% and 35%. The plastic potential for FKM criterion reads
here, σe is the effective stress and pm is the pressure.The flow strength of the material under hydrostatic loading, pf, is computed from
here, σ0 is the initial yield stress, and φ is the void volume fraction (porosity). The maximum void volume fraction φm typically takes the value of 36%, the limit of dense random packing of sintered powder.
FKM-GTN Criterion
The FKM-GTN criterion is a combination of the Fleck-Kuhn-McMeeking Criterion and Gurson-Tvergaard-Needleman Criterion, intended to cover a wider range of porosities (Ref. 20-21). For low void volume fractions (porosity lower than 10%), the GTN model is used and for void volume fractions higher than 25%, the FKM criterion is used. In the transition zone, a linear combination of both criteria is used.
Soil Plasticity
Mohr-Coulomb Criterion
The Mohr-Coulomb criterion is the most popular criterion in soil mechanics. It was developed by Coulomb before the Tresca and von Mises criteria for metals, and it was the first criterion to account for the hydrostatic pressure. The criterion states that failure occurs when the shear stress and the normal stress acting on any element in the material satisfy the equation
here, τ is the shear stress, c the cohesion, and denotes the angle of internal friction.
With the help of Mohr’s circle, this criterion can be written as
The Mohr-Coulomb criterion defines an irregular hexagonal pyramid in the space of principal stresses, which generates singularities in the derivatives of the yield function.
Figure 3-6: The Mohr-Coulomb criterion. The cone opens toward the compressive axis.
The Mohr-Coulomb criterion can be written in terms of the invariants I1 and J2 and the Lode angle 0 ≤ θ ≤ π/3 (Ref. 1, Ref. 9) when the principal stresses are sorted as σ1 ≥ σ2 ≥ σ3. The yield function then reads
The tensile meridian is defined when θ = 0 and the compressive meridian when θ = π/3.
Rearranging terms, the Mohr-Coulomb criterion reads
where
, , and
In the special case of frictionless material, (, α = 0, k = c), the Mohr-Coulomb criterion reduces to a Tresca’s maximum shear stress criterion, 1 − σ3= 2k or equivalently
See also the description of the Mohr-Coulomb material model in the Solid Mechanics interface documentation.
Drucker-Prager Criterion
The Mohr-Coulomb criterion causes numerical difficulties when treating the plastic flow at the corners of the yield surface. The Drucker-Prager model neglects the influence of the invariant J3 (introduced by the Lode angle) on the cross-sectional shape of the yield surface. It can be considered as the first attempt to approximate the Mohr-Coulomb criterion by a smooth function based on the invariants I1 and J2 together with two material constants (which can be related to Mohr-Coulomb’s coefficients)
This is sometimes also called the extended von Mises criterion, since it is equivalent to the von Mises criterion for metals when setting α = 0.
Figure 3-7: The Drucker-Prager criterion. The cone opens toward the compressive axis.
The coefficients in the Drucker-Prager model can be matched to the coefficients in the Mohr-Coulomb criterion by
and
The symbol ± is related to either matching the tensile meridian (positive sign) or the compressive meridian (negative sign) of Mohr-Coulomb’s pyramid.
The matching at the tensile meridian (θ = 0) comes from setting
in the Mohr-Coulomb criterion, and the matching at the compressive meridian (θ = π/3) from setting
Figure 3-8: The Drucker-Prager criterion showing the tensile and compressive meridians (inner and outer circles), and the Lode angle compared to the cross section of Mohr-Coulomb criterion in the π-plane.
In the special case of frictionless material, (, α = 0, ), the Drucker -Prager criterion reduces to the von Mises criterion
When matching Drucker-Prager criterion to Mohr-Coulomb criterion in 2D plane-strain applications, the parameters are
and
and when matching both criteria in 2D plane-stress applications, the matching parameters are:
and
Dilatation Angle
The Mohr-Coulomb yield criterion is sometimes used with a non-associated plastic potential. This plastic potential could be either a Drucker-Prager criterion, or the same Mohr-Coulomb yield function but with a different slope with respect to the hydrostatic axis, in which case the angle of internal friction is replaced by the dilatation angle, which is normally smaller (Ref. 5).
Also, when using a Drucker-Prager criterion matched to a Mohr-Coulomb criterion, the plastic potential could also be non-associated, in which case the difference between the dilatation angle and the angle of internal friction would result in a yield surface and plastic potential portrayed by two cones with different angles with respect to the hydrostatic angle.
Elliptic Cap
The Mohr-Coulomb and Drucker-Prager criteria portray a conic yield surface which opens in the hydrostatic axis direction. Normally, these soil models are not accurate above a given limit pressure because real-life materials cannot bear infinite loads and still behave elastically. A simple way to overcome this problem is to add an elliptical end-cap on the compressive side to these models.
The elliptic cap is an elliptic yield surface of semi-axes as shown in Figure 3-9. The initial pressure pa (SI units: Pa) denotes the pressure at which the elastic range circumscribed by either a Mohr-Coulomb pyramid or a Drucker-Prager cone is not valid any longer, so a cap surface is added. The limit pressure pb gives the curvature of the ellipse, and denotes the maximum admissible hydrostatic pressure for which the material starts deforming plastically. Pressures higher than pb are not allowed
Figure 3-9: Elliptic cap model in Haigh–Westergaard coordinate system.
Note that the sign convention for the pressure is taken from the Structural Mechanics Module: positive sign under compression, so pa and pb are positive parameters. Figure 3-9 shows the cap in terms of the variables p and q
and
In terms of these variables, the equation for the elliptic cap reads
the point (pa, qa) in the Haigh–Westergaard coordinate system is where the elliptic cap intersects either the Mohr-Coulomb or the Drucker-Prager cone.
Elliptic Cap With Hardening
It is also possible to add isotropic hardening to the cap surface. In this case, the center of the ellipse is shifted as the volumetric plastic strain increases, also, the size of the ellipse’s semi-axes grow as hardening evolves. The intersection of the elliptic cap with the pressure axis is given by
here, pb0 is the initial value for the limit pressure pb, Kiso is the isotropic hardening modulus, εpvol the volumetric plastic strain, and εpvol,max the maximum volumetric plastic strain. Instead of providing the value for the initial pressure pa (SI units: Pa), the ellipse’s aspect ratio R is entered.
Note that the volumetric plastic strain εpvol is negative in compression, so the limit pressure pb is increased from pb0 as hardening evolves.
See also the description of the Drucker-Prager material model in the Solid Mechanics interface documentation.
Matsuoka-Nakai Criterion
Matsuoka and Nakai (Ref. 3) discovered that the sliding of soil particles occurs in the plane in which the ratio of shear stress to normal stress has its maximum value, which they called the mobilized plane. They defined the yield surface as
where the parameter μ = (τ/σn)STP equals the maximum ratio between shear stress and normal stress in the spatially mobilized plane (STP-plane), and the invariants are applied over the effective stress tensor (this is the Cauchy stress tensor minus the fluid pore pressure).
The Matsuoka-Nakai criterion circumscribes the Mohr-Coulomb criterion in dry soils, when
and denotes the angle of internal friction in Mohr-Coulomb criterion.
Figure 3-10: The Matsuoka-Nakai criterion and Mohr-Coulomb criterion in the principal stress space.
See also the description of the Matsuoka-Nakai material model in the Solid Mechanics interface documentation.
Lade-Duncan Criterion
The Lade-Duncan criterion was originally developed to model a large volume of laboratory sample test data of cohesionless soils. This criterion is defined as
where I1 and I3 are the first and third stress invariants respectively, and k is a parameter related to the direction of the plastic strain increment in the triaxial plane. The parameter k can vary from 27 for hydrostatic stress conditions (σ1 = σ2 = σ3), up to a critical value kc at failure. In terms of the invariants I1, J2,and J3, this criterion can be written as
The Lade-Duncan criterion can be fitted to the compressive meridian of the Mohr-Coulomb surface by choosing
with as the angle of internal friction in Mohr-Coulomb criterion
Figure 3-11: Comparing the Mohr-Coulomb, Matsuoka-Nakai, and Lade-Duncan criteria when matching the tensile meridian.
See also the description of the Lade-Duncan material model in the Solid Mechanics interface documentation.
Tension Cut-Off
It appears that the Mohr-Coulomb and Drucker-Prager criteria predict tensile strengths larger than the experimental measurements on soil samples. This discordance can be mended by the introduction of the Rankine or tension cut-off criterion.
The Rankine criterion states that a material stops deforming elastically when the biggest principal stress σ1 reaches a maximum tensile stress, also called tension cut-off limit σt.
In terms of the principal stress, Rankine criterion reads
For soils and clays, the maximum tensile stress can be estimated from the material parameters, such as the cohesion c and the friction angle . For instance, the tip of the cone in Mohr-Coulomb criterion is reached when
therefore, the tension cut-off should be chosen such as
The Mohr-Coulomb criterion together with a tension cut-off is sometimes called modified Mohr-Coulomb criterion (Ref. 19).
Hill Orthotropic Plasticity
Hill (Ref. 12, Ref. 13) proposed a quadratic yield function (and associated plastic potential) in a local coordinate system given by the principal axes of orthotropy ai
(3-23)
The six parameters F, G, H, L, M, and N are related to the state of anisotropy. As with isotropic plasticity, the elastic region Qp < 0 is bounded by the yield surface Qp = 0.
Hill demonstrated that this type of anisotropic plasticity is volume preserving, this is, given the associated flow rule
the trace of the plastic strain rate tensor is zero, which follows from the expressions for the diagonal elements of
so the plastic volumetric strain rate is zero
Hill plasticity is an extension of J2 (von Mises) plasticity, in the sense that it is volume preserving. Due to this assumption, six parameters are needed to define orthotropic plasticity, as opposed to orthotropic elasticity, where nine elastic coefficients are needed.
Expressions for the Coefficients F, G, H, L, M, N
Hill noticed that the parameters L, M, and N are related to the yield stress in shear with respect to the axes of orthotropy ai, thus they are positive parameters
, ,
Here, σysij represents the yield stress in shear on the plane ij.
The material parameters σys1, σys2, and σys3 represent the tensile yield stress in the direction, a1, a2, and a3, and they are related to Hill’s parameters F, G, and H as
or equivalently
Note that at most one of the three coefficients F, G, and H can be negative.
In order to define a yield function and plastic potential suitable for isotropic or kinematic hardening, the average initial yield stress σys0 is calculated from the Hill’s parameters F, G, and H (this is equivalent to the initial yield stress σys0 in von Mises plasticity)
(3-24)
Defining Hill’s effective stress as (Ref. 13)
makes it possible to write the plastic potential in a similar way to von Mises plasticity.
Isotropic hardening is then applied on the average yield stress variable σys0, by using the plastic potential
Here, the average yield stress
now depends on the initial yield stress σys0, the hardening function σh, and the effective plastic strain εpe.
Isotropic Hardening
Plasticity implements seven different kinds of isotropic hardening models for elastoplastic materials:
Perfect (or Ideal) Plasticity
For perfect elastoplastic materials, the yield surface is fixed in the space of principal stresses, and therefore, plastic deformations occur only when the stress path moves on the yield surface (the regime inside the yield surface is elastic, and stress paths beyond the yield surface are not allowed).
In this case the plasticity algorithm solves either the associated or non-associated flow rule for the plastic potential Qp
with the yield function
In the settings for plasticity you specify the effective stress φ(σ) for the yield function from von Mises stress, Tresca stress, Hill effective stress, or a user-defined expression; and σys0 is the initial yield stress that defines the onset of plastic deformation.
When Large plastic strain is selected as the plasticity model for the Plasticity node, either the associate or non-associated flow rule is applied as written in Equation 3-30.
Linear Isotropic Hardening
In this case the plasticity algorithm solves either the associated or non-associated flow rule for the plastic potential Qp
with the yield function
where the yield stress σyspe) now depends on the effective plastic strain εpe.
The yield stress σyspe) is then a function of the effective plastic strain and the initial yield stress σys0
here, the isotropic hardening modulus Eiso is calculated from
For linear isotropic hardening, the isotropic tangent modulus ETiso is defined as (stress increment / total strain increment). A value for ETiso is entered in the isotropic tangent modulus section for the Plasticity node. The Young’s modulus E is taken from the parent material (Linear Elastic, Nonlinear Elastic or Hyperelastic material model). For orthotropic and anisotropic elastic materials, E represents an effective Young’s modulus.
Ludwik
In Ludwik model for nonlinear isotropic hardening, the yield stress σyspe) is defined by a nonlinear function of the effective plastic strain. Ludwik equation (also called Ludwik-Hollomon equation) for isotropic hardening is given by the power-law
here, k is the strength coefficient and n is the hardening exponent. Setting n = 1 would result in linear isotropic hardening.
Swift
For non-saturating materials, the Swift power-law equation relates the initial yield stress σys0 and the isotropic hardening σh, to the effective plastic strain as
here, k is the strength coefficient, n is the hardening exponent and ε0 is a reference strain. Noting that at zero plastic strain the initial yield stress is related to the strength coefficient and hardening exponent as
the yield stress is written as
Voce
Voce rule for nonlinear isotropic hardening is intended for materials that exhibit a saturating evolution of hardening. The isotropic hardening σh is exponentially related to the effective plastic strain as
The yield stress σyspe) is then defined as
The value of the saturation exponent parameter β determines the saturation rate of the hysteresis loop for cyclic loading. The saturation flow stress σsat characterizes the maximum distance by which the yield surface can expand in the stress space. For values εpe >> 1/β, the yield stress saturates to
Hockett-Sherby
Hockett-Sherby rule for nonlinear isotropic hardening is also intended for materials which yield stress saturates as effective plastic strain increases. It is similar to Voce rule, but it includes an exponential dependency of the form
where σ is the steady-state flow stress, m the saturation coefficient and n the saturation exponent. For values mεpen >> 1 the yield stress saturates to
User defined
The yield stress versus the effective plastic strain can be specified with the help of a hardening curve that could also depend on other variables, such as temperature.
In this case, define the (usually nonlinear) hardening function σhpe) such that the yields stress reads
The internal variable for the effective plastic strain is named solid.epe. The effective plastic strain evaluated at Gauss points is named solid.epeGp.
When Large plastic strain is selected as the plasticity model for the Plasticity node, either the associate or non-associated flow rule is applied as written in Equation 3-30.
Kinematic Hardening
There are few options for computing either linear or nonlinear kinematic hardening for plasticity:
For any of the kinematic hardening models, the algorithm solves either the associated or non-associated flow rule for the plastic potential Qp
with the yield function defined as
Here, σys is the yield stress (which may include a linear or nonlinear isotropic hardening model), and the effective stress φ(σ) is either the von Mises, Tresca, or Hill stress; or a user-defined expression. The stress tensor used in the yield function is shifted by what is usually called the back stress, σback.
Linear Kinematic Hardening
The back stress is generally not only a function of the current plastic strain but also of its history. In the case of linear kinematic hardening, the back stress σback is a linear function of the plastic strain tensor εp, this is also known as Prager’s hardening rule.
The implementation of linear kinematic hardening assumes a linear evolution of the back stress tensor with respect to the plastic strain tensor:
where the kinematic hardening modulus Ck is calculated from
The value for Ek is entered in the kinematic tangent modulus section and the Young’s modulus E is taken from the linear or nonlinear elastic material model. For orthotropic and anisotropic elastic materials, E represents an averaged Young’s modulus. Note that some authors define the kinematic hardening modulus as Hk = 2/3Ck.
Armstrong-Frederik Hardening Model
Armstrong and Frederick (Ref. 7) added memory to Prager’s linear kinematic hardening model. This nonlinear kinematic hardening model allows to capture the Bauschinger effect and nonlinear behavior by non symmetrical tension-compression loading.
The nonlinear evolution of the back stress σback is governed by the rate
here, Ck is the kinematic hardening modulus, γk is a kinematic hardening parameter, and εpe the effective plastic strain. Setting γk = 0 recovers Prager’s rule for linear kinematic hardening.
To solve this rate, internal degrees of freedom are added to account for the back stress components. In order to have the same units as used for the plastic strain, the algorithm solves for the back strain εback as proposed in (Ref. 8), which is related to the back stress as
The nonlinear evolution for the back strain reads
Chaboche Hardening Model
Chaboche (Ref. 8) proposed a nonlinear kinematic hardening model based on the superposition of N back stress tensors
each of these back stress tensors σback,i follows a nonlinear Frederick-Armstrong kinematic hardening rule
Practitioners would normally select γk = 0 for one of the back stress equations, thus recovering Prager’s linear rule for that branch
The back stress tensor σback is then defined by the superposition of N back stress tensors
As done for Armstrong-Frederik kinematic hardening, the algorithm solves for the back strain tensors εback,i instead of the back stress tensors. The change of variables is
and the nonlinear evolution for the back strain tensors reads
Introduction to Small and Large Plastic Strains
There are two implementations of plasticity available in COMSOL Multiphysics. One is based on the additive decomposition of strains, which is the most suitable approach in the case of small strains, and the other one is based on the multiplicative decomposition of the deformation gradient, which is more suitable when large plastic strains occurs. The additive and multiplicative decomposition of strains is described in Inelastic Strain Contributions.
When small plastic strain is selected as the plasticity model, an additive decomposition is used. If the elastic or plastic strains are large, the additive decomposition might produce incorrect results. As an example, the volume preservation, which is an important assumption in metal plasticity, will no longer be respected. The additive decomposition of strains is however widely used both for metal and soil plasticity.
When Small plastic strain is selected as the plasticity model for the Plasticity node, and the Include geometric nonlinearity check box is selected on the study Settings window, a Cauchy stress tensor is used to evaluate the yield function and plastic potential. The components of this stress tensor are oriented along the material directions, so it can be viewed as a scaled second Piola-Kirchhoff stress tensor. The additive decomposition of strains is understood as the summation of Green-Lagrange strains.
When large plastic strain is selected as the plasticity model, the total deformation gradient tensor is multiplicatively decomposed into an elastic deformation gradient and a plastic deformation gradient.
Plastic Flow for Small Strains
The flow rule defines the relationship between the increment of the plastic strain tensor  and the current state of stress, σ, for a yielded material subject to further loading. When Small plastic strain is selected as the plasticity model for the Plasticity node, the direction of the plastic strain increment is defined by
Here, λ is a positive multiplier (also called the consistency parameter or plastic multiplier), which depends on the current state of stress and the load history, and Qp is the plastic potential.
The “dot” (for ) means the rate at which the plastic strain tensor changes with respect to Qp/∂σ. It does not represent a true time derivative. Some authors call this formulation rate independent plasticity.
The direction of the plastic strain increment is perpendicular to the surface (in the hyperspace spanned by the stress tensor components) defined by the plastic potential Qp.
The plastic multiplier λ is determined by the complementarity or Kuhn-Tucker conditions
(3-25), and
where Fy is the yield function. The yield surface encloses the elastic region defined by Fy < 0. Plastic flow occurs when Fy = 0.
If the plastic potential and the yield surface coincide with each other (Qp = Fy), the flow rule is called associated, and the rate in Equation 3-26 is solved together with the conditions in Equation 3-25.
(3-26)
For a non-associated flow rule, the yield function does not coincide with the plastic potential, and together with the conditions in Equation 3-25, the rate in Equation 3-27 is solved for the plastic potential Qp (often, a smoothed version of Fy).
(3-27)
The evolution of the plastic strain tensor (with either Equation 3-26 or Equation 3-27, plus the conditions in Equation 3-25) is implemented at Gauss points in the plastic element elplastic.
Plastic Flow for Large Strains
When Large plastic strain is selected as the plasticity model for the Plasticity node, a multiplicative decomposition of deformation (Ref. 9, Ref. 10, and Ref. 11) is used, and the associated plastic flow rule can be written as the Lie derivative of the elastic left Cauchy-Green deformation tensor Bel:
(3-28)
The plastic multiplier λ and the yield function Φ (written in terms of the Kirchhoff stress tensor τ) satisfy the Kuhn-Tucker condition, as done for infinitesimal strain plasticity
, and
The yield function Φ in Ref. 9 and Ref. 10 was written in terms of Kirchhoff stress τ and not Cauchy stress σ because the authors defined the plastic dissipation with the conjugate energy pair τ and d, where d is the rate of strain tensor.
The Lie derivative of Bel is then written in terms of the plastic right Cauchy-Green rate
(3-29)
By using Equation 3-28 and Equation 3-29, the either associated or non-associated plastic flow rule for large strains is written as (Ref. 10)
(3-30)
together with the Kuhn-Tucker conditions for the plastic multiplier λ and the yield function Fy
(3-31), and
For the associated flow rule, the plastic potential and the yield surface coincide with each other (Qp = Fy), and for the non-associated case, the yield function does not coincide with the plastic potential.
In COMSOL Multiphysics, the elastic left Cauchy-Green tensor is written in terms of the deformation gradient and the right Cauchy-Green tensor, so Bel = FCp1FT. The flow rule then reads
(3-32)
The plastic flow rule is then solved at Gauss points in the plastic element elplastic for the inverse of the plastic deformation gradient Fp1, so that the variables in Equation 3-32 are replaced by
, and
After integrating the flow rule in Equation 3-32, the plastic Green-Lagrange strain tensor is computed from the plastic deformation tensor
and the elastic Green-Lagrange strain tensor is computed from the elastic deformation gradient tensor Fel = FFp1
When Large plastic strain is selected as the plasticity model for the Plasticity node, the effective plastic strain variable is computed as the true effective plastic strain (also called Hencky or logarithmic plastic strain).
When either Large plastic strain or Small plastic strain is selected as the plasticity model for the Plasticity node, the out-of-plane shear strain components are not computed in 2D, neither for plane stress nor plane strain assumption.
Numerical Solution of the Elastoplastic Conditions
A backward Euler discretization of the pseudo-time derivative is used in the plastic flow rule. For small plastic strains, this gives
where “old” denotes the previous time step and Λ = λΔt, where Δt is the pseudo-time step length.
For large plastic strains, Equation 3-32 is numerically solved with the so-called exponential mapping technique
where and .
For each Gauss point, the plastic state variables (εp or Fp1, depending on whether small strain or large strain plasticity is selected) and the plastic multiplier, Λ, are computed by solving either of the above time-discretized flow rules together with the complementarity conditions
This is done as follows (Ref. 9):
1
Elastic-predictor: Try the elastic solution εp = εp,old (or ) and Λ = 0. If this satisfies Fy ≤ 0 it is done.
2
Plastic-corrector: If the elastic solution does not work (this is Fy > 0), solve the nonlinear system consisting of the flow rule and the equation Fy = 0 using a damped Newton method..
The numerical tolerance to fulfill the condition Fy = 0 is given in SI units of Pascals, and it depends on the initial yield stress (in case of plasticity and porous plasticity) or it is defined in terms of other material parameter (for soil plasticity). This numerical tolerance is 0.1% the value defined in the variable item.tol., where item is the name of the node.
Energy Dissipation
Since plasticity is an inelastic process, the dissipated energy density can be calculated by integrating the pseudo-rate given by
As plasticity is rate independent, the plastic dissipation density Wp is obtained after integrating an extra variable in the plastic flow rule.
The total energy dissipated by plasticity in a given volume can be calculated by a volume integration of the plastic dissipation density Wp.
When the Calculate dissipated energy check box is selected, the plastic dissipation density is available under the variable solid.Wp and the total plastic dissipation under the variable solid.Wp_tot.