Elastoplastic Materials
In this section:
Plasticity
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 a Linear Elastic Material Nonlinear Elastic Materials, or Hyperelastic Materials.
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 Criterion
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 .
Plastic Flow
The evolution of plastic deformation is typically referred as plastic flow. There are two main implementations for plastic flow available in COMSOL Multiphysics depending on what strain decomposition is selected. An additive decomposition of strains, is the most suitable approach in the case of small strains, and a multiplicative decomposition of the deformation gradient is more suitable when strains and/or rotations are expected to be large. The additive and multiplicative decomposition of strains is described in Inelastic Strain Contributions.
When an additive decomposition is used while the elastic or plastic strains are large, it 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
The type of strain decomposition used can be controlled in the Geometric Nonlinearity section of the material model.
If the Include geometric nonlinearity check box is selected in Settings window of the study, the default is to use a multiplicative decomposition.
Plastic Flow for Additive Decomposition
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. The direction of the plastic flow 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-49), 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-50 is solved together with the conditions in Equation 3-49.
(3-50)
For a nonassociated flow rule, the yield function does not coincide with the plastic potential, and together with the conditions in Equation 3-49, the rate in Equation 3-51 is solved using the plastic potential Qp (often, a smoothed version of Fy).
(3-51)
When additive plasticity is used and the Include geometric nonlinearity check box is selected on the study Settings window, the Cauchy stress tensor is used to evaluate the yield function and plastic potential. Plastic flow is, however, assumed in the direction the second Piola–Kirchhoff stress tensor, and the computed plastic strain tensor can be interpreted as a Green–Lagrange strain.
Plastic Flow for Multiplicative Decomposition
For a large strain formulation with multiplicative decomposition of deformation (Ref. 3, Ref. 4, and Ref. 5), the plastic flow rule can be written as the Lie derivative of the elastic left Cauchy–Green deformation tensor Bel:
(3-52)
The plastic multiplier λ and the plastic potential Qp (written in terms of the Cauchy stress tensor σ) satisfy the Kuhn–Tucker condition, as done for infinitesimal strain plasticity
, and
The yield function Fy in Ref. 3 and Ref. 4 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-53)
By using Equation 3-52 and Equation 3-53, the plastic flow rule is written as (Ref. 4)
(3-54)
For the associated flow rule, the plastic potential and the yield surface coincide with each other (Qp = Fy), and for the nonassociated 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 F and the plastic right Cauchy–Green tensor Cp, so Bel = FCp1FT. The flow rule then reads
(3-55)
The plastic flow rule is then solved at Gauss points to find the inverse of the plastic deformation gradient Fp1, by replacing the variables in Equation 3-55 with
, , and
After integrating the flow rule in Equation 3-55, 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
The equivalent plastic strain variable is computed as the true equivalent plastic strain (also called Hencky or logarithmic plastic strain) when a multiplicative decomposition is used.
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-56)
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.
A common measure of inelastic deformation is the equivalent plastic strain rate, (or equivalent von Mises strain rate) which is defined from the plastic shear components
(3-57)
Some authors define the equivalent plastic strain rate directly from the yield surface and the yield stress σys
(3-58)
and it is also possible to define it from a user defined expression hp
(3-59)
Incompressible plastic deformation is experimentally observed in metals, but this is not the case for most materials used in geotechnical applications. For instance, a nonzero increment in volumetric plastic strain, , is explicitly used in Porous Plasticity and Elastoplastic Soil Models.
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 equivalent stress. The default form of the equivalent stress is the von Mises stress, which is often used in metal plasticity:
Other expressions can be defined, such as Tresca stress, Hill equivalent stress, or a user defined expression.
The Tresca equivalent 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), p = −I1(σ)/3, the equivalent von Mises stress σmises, or other invariants, principal stresses, or stress tensor components.
The 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)
The von Mises criterion is implemented with the yield surface
where σys is the yield stress level (yield stress in uniaxial tension) and is the equivalent von Mises stress.
The equivalent von Mises stress σmises is available in the variable solid.mises, where solid is the name of the physics interface node.
The 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 are sorted as σ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, the Tresca 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.
The Tresca criterion can be used with either an associated or nonassociated flow rule, in which case von Mises stress is applied in order to get better numerical performance.
The Tresca equivalent stress, σtresca = σ1 − σ3, is implemented in the variable solid.tresca, where solid is the name of the physics interface node.
Figure 3-14: The upper and lower limits of the Tresca criterion.
Figure 3-15: Classical yield criteria for metals. The Tresca criterion (left) and the 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 belong to what researchers call volume preserving or J2 plasticity, as the plastic flow is independent on the mean pressure.
Orthotropic Plasticity
Hill (Ref. 6, Ref. 7) proposed a quadratic yield function (and associated plastic potential) in a local coordinate system given by the principal axes of orthotropy ai
(3-60)
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-61)
Defining Hill’s equivalent stress as (Ref. 7)
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 equivalent plastic strain εpe.
Isotropic Hardening
In the case of isotropic hardening, the plasticity problem is defined with either associated or nonassociated plastic flow and with the yield surface
where the yield stress σyspe) is a function of the equivalent plastic strain εpe.
There are several different kinds of isotropic hardening models for elastoplastic and viscoplastic materials:
Perfect Plasticity (no isotropic hardening)
The internal variable for the equivalent plastic strain is named solid.epe. The equivalent plastic strain evaluated at Gauss points is named solid.epeGp.
Perfect Plasticity
For perfect (or ideal) 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).
Here the yield stress is directly given by the initial yield stress σys0 such that
In the settings for plasticity you specify the equivalent stress σe(σ) for the yield function; and σys0 that defines the onset of plastic deformation.
Linear Isotropic Hardening
The yield stress σyspe) is for linear isotropic hardening a function of the equivalent plastic strain εpe 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 equivalent Young’s modulus.
Ludwik
In the Ludwik model for nonlinear isotropic hardening, the yield stress σyspe) is defined by a nonlinear function of the equivalent plastic strain. The Ludwik equation (also called the 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.
Johnson–Cook
Similar to the Ludwik hardening model, the yield stress σyspe) in the Johnson–Cook model is defined as a nonlinear function of the equivalent plastic strain εpe given by a power law. The difference is that effects of the strain rate are included in the model. The yield stress and hardening function is given by
(3-62)
The Johnson–Cook model also adds the possibility to include thermal softening by including a function which depends on the normalized homologous temperature Th. It is a function of the current temperature T, the melting temperature of the metal Tm, and a reference temperature Tref
As a default, this dependency follows a power law, . The softening function f(Th) should have the properties f(Th < 0) = 0 and f(Th1) = 1.
The natural logarithm is used in Equation 3-62. If material data have been determined using a base 10 logarithm, the value of C must be adjusted accordingly.
The strain rate strength coefficient C and the reference strain rate in the Johnson–Cook model are obtained from experimental data. The flow curve properties are normally investigated by conducting uniaxial tension tests in a temperature range around half the melting temperature Tm, and subjected to different strain rates like 0.1/s, 1/s, and 10/s.
It should be noted that the intent of the strain rate dependent term in Equation 3-62 is to capture the hardening at high strain rates. If extrapolated to very low strain rates, a low or even negative yield stress will be predicted. Low strain rates can however be expected in some regions in a general finite element model. To resolve this, the material model as implemented, never evaluates the logarithmic term to a value smaller than zero. This means that the reference strain rate must be considered the quasistatic limit. Below this strain rate, the hardening function will be a strain rate independent, but temperature dependent, Ludwik law:
Swift
For nonsaturating materials, the Swift power-law equation relates the initial yield stress σys0 and the isotropic hardening σh, to the equivalent 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
for
Voce
The 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 equivalent 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
The Hockett–Sherby rule for nonlinear isotropic hardening is also intended for materials which yield stress saturates as equivalent 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
Hardening Function
The yield stress versus the equivalent 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
Kinematic Hardening
For kinematic hardening, the plasticity problem is defined with either associated or nonassociated plastic flow and with the yield surface
(3-63)
Here, σys is the yield stress (which may include a linear or nonlinear isotropic hardening model), and σe (σ) is the equivalent stress. The stress tensor used in the yield function is shifted by what is usually called the back stress, σb to impose kinematic hardening.
There are few options for computing either linear or nonlinear kinematic hardening for plasticity:
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.
When kinematic hardening is added, both the plastic potential and the yield surface are calculated with effective invariants, that is, the invariants of the tensor defined by the difference between the stress tensor minus the back-stress, σeff = σ - σb. The invariant of effective deviatoric tensor is named solid.II2sEff, which is used when a von Mises, Tresca, or Hill orthotropic plasticity is computed together with kinematic hardening.
Armstrong–Frederick Hardening Model
Armstrong and Frederick (Ref. 1) added memory to Prager’s linear kinematic hardening model. This nonlinear kinematic hardening model makes it possible to capture the Bauschinger effect and nonlinear behavior by nonsymmetric tension-compression loading.
The nonlinear evolution of the back stress σb is governed by the rate
here, Ck is the kinematic hardening modulus, γk is a kinematic hardening parameter, and εpe the equivalent 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 εb as proposed in (Ref. 2), which is related to the back stress as
The nonlinear evolution for the back strain reads
Chaboche Hardening Model
Chaboche (Ref. 2) proposed a nonlinear kinematic hardening model based on the superposition of N back stress tensors
each of these back stress tensors σb,i follows a nonlinear Armstrong–Frederick kinematic hardening rule
(3-64)
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 σb is then defined by the superposition of N back stress tensors
As done for Armstrong–Frederick kinematic hardening, the algorithm solves for the back strain tensors εb,i instead of the back stress tensors. The change of variables is
and the nonlinear evolution for the back strain tensors in Equation 3-64 reads
Multiplicative Decomposition
For multiplicative decomposition, the total deformation gradient F is decomposed into an elastic Fel and an inelastic (plastic) Fpl deformation
as described in Inelastic Strain Contributions. When kinematic hardening is used with multiplicative decomposition, Fpl is similarly split multiplicatively as suggested in Ref. 8, such that
Here Fple is referred to as the “elastic” part of Fpl, and Fplv as its “viscous” part of the plastic deformation. Conceptually, Fple corresponds to a spring and Fplv to a dashpot that are connect in series in a rheological model. Now new kinematic quantities can be defined in a similar fashion as for Fel. For example, the tensor
which is analogous to the elastic right Cauchy–Green deformation tensor Cel. As for the elastic part of the model, a free energy function Wkin(Cpl,e) can be defined for the spring of the kinematic hardening part. From Wkin different stress tensors can be derived, for example
which represents that back stress in the material frame and is analogous to the second Piola-Kirchhoff stress of the elastic part. By applying a push-forward operation by the total deformation gradient F, the back stress in the spatial frame is
Alternatively using the plastic deformation gradient leads
where is a representation of the back stress in an intermediate frame defined by Fpl. The yield function in Equation 3-63 is defined in the spatial frame for large strains, thus to complete the formulation we need to find σb, or Cple, which depends on the kinematic hardening model used.
It is in COMSOL Multiphysics assumed that Wkin has the same format as the elastic strain energy function and that it is always isotropic. Hence for a linear elastic material
where εple = 0.5(Cple-1) and Ck is the kinematic hardening modulus. For hyperelastic materials, a compressible Neo-Hookean strain energy function is used
The flow rule suggested in Plastic Flow for Multiplicative Decomposition leads to an arbitrary plastic rotation, which is undesirable with the above formulation for large strain kinematic hardening. Hence, when kinematic hardening is used, the plastic flow rule is modified to
such that the plastic rotation is uniquely determined.
Linear Kinematic Hardening
To obtain a linear kinematic hardening model the inelastic part of the plastic deformation is omitted by setting Fplv = 1, and thus Cple = Cpl.
Armstrong–Frederick Kinematic Hardening
To obtain an Armstrong–Frederick kinematic hardening model, the inelastic part of the plastic deformation Fplv is considered as an additional unknown variable determined by the flow rule
where γk is a kinematic hardening parameter. The actual unknown used by the plasticity algorithm is the symmetric part of the inverse of Fplv.
Chaboche Kinematic Hardening
The Chaboche kinematic hardening model can be viewed as a combination of linear kinematic and Armstrong–Frederick kinematic hardening by adding N branches. Formally each branch here corresponds to a unique multiplicative split of Fpl such that for branch i we have
Conceptually, each branch corresponds to a spring and a dashpot in series and the branches are then connected parallel. The first branch always corresponds to linear kinematic hardening such that Fplv,1 = 1. Each consecutive branch adds an Armstrong–Frederick kinematic hardening model, where Fplv,i is an additional unknown with its own flow rule. The total back stress is then the sum of the stress in all branches
Soil Plasticity
The 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-16: 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. 11) 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, α = 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.
The 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 (Ref. 1, Ref. 2) 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)
The criterion 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-17: The Drucker–Prager criterion. The cone opens toward the compressive axis.
The coefficients in the Drucker–Prager model are related to the cohesion c and angle of internal friction ϕ in the Mohr–Coulomb criterion by the relation
The coefficients in the Drucker–Prager model can be matched to the coefficients in the Mohr–Coulomb criterion by (Ref. 1)
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-18: The Drucker–Prager criterion showing the tensile and compressive meridians (inner and outer circles), and the Lode angle compared to the cross section of the Mohr–Coulomb criterion in the π-plane.
In the special case of frictionless material, (ϕ = 0, α = 0, ), the Drucker–Prager criterion reduces to the von Mises criterion
When the Drucker–Prager criterion is matched to the Mohr–Coulomb criterion in 2D plane-strain applications, the parameters are
and
when matching the criterion in 2D generalized plane-strain, the matching parameters are:
and
and when matching the criterion in 2D plane-stress, the matching parameters are:
and
Dilatation Angle
The Mohr–Coulomb criterion is sometimes used with a nonassociated 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. 7).
Also, when using a Drucker–Prager criterion matched to the Mohr–Coulomb criterion, the plastic potential could also be nonassociated, in which case the difference between the dilatation angle and the angle of internal friction results 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-19. 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-19: Elliptic cap model in Haigh–Westergaard coordinate system.
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-19 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.
When the Mohr–Coulomb or the Drucker–Prager equations are represented in  − I1 plane, the equation for the elliptic cap reads
where Ia = −3pa, Ib = −3pb, and Ja is the ellipse semiaxis in the  − I1 plane, given by the intersection with the Drucker–Prager criterion
This makes it easier to compute the ellipse’s semiaxis Ja in terms of Drucker–Prager parameters α and k.
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. The volumetric plastic strain εpvol is negative in compression, so the limit pressure pb is increased from pb0 as hardening evolves.
Instead of providing the value for the initial pressure pa (SI units: Pa), the ellipse’s aspect ratio R is entered. The ellipse aspect ratio R is given by
where
The ellipse semiaxis in the  − I1 plane is given by the aspect ratio R
In order to have a smooth transition, the cap intersects the failure surface at the point of tangency (Ref. 3). The intersection coordinates (Ic, Jc) are given by
and
 
See also the description of the Drucker–Prager material model in the Solid Mechanics interface documentation.
The Matsuoka–Nakai Criterion
Matsuoka and Nakai (Ref. 5) 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 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 the Mohr–Coulomb criterion.
Figure 3-20: The Matsuoka–Nakai and Mohr–Coulomb criteria in the principal stress space.
See also the description of the Matsuoka–Nakai material model in the Solid Mechanics interface documentation.
The 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 the Mohr–Coulomb criterion.
Figure 3-21: 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 Cutoff
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 a tension cutoff criterion.
The Rankine criterion states that a material stops deforming elastically when the largest principal stress σ1 reaches a maximum tensile stress, also called tension cutoff limit σt.
In terms of the principal stress, the 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 cutoff should be chosen such as
The Mohr–Coulomb criterion together with a tension cutoff is sometimes called the modified Mohr–Coulomb criterion (Ref. 21).
There might be numerical problems when implementing Rankine cutoff criterion, as the intersection of the cone in Mohr–Coulomb criterion with the tensile plane produces sharp edges in the stress space. A better approach is to use the mean stress instead of the principal stress, then the mean stress cutoff reads
here, σm is the maximum mean stress and p = −I1/3 is the pressure, which is negative in tension.
Porous Plasticity
The modeling of plastic deformation in soils, porous metals and aggregates has a main difference with respect to traditional metal plasticity; the yield function and plastic potential are not only defined in terms of the deviatoric stress tensor (or the deviatoric stress invariant J2), but also include dependencies on the hydrostatic pressure.
A key concept for porous plasticity models is the evolution of the relative density, which is the solid volume fraction in a porous mixture. The relative density is related to the porosity (or void volume fraction) ϕ by
When compacting a mixture of metal particles, the porosity tends to zero and the relative density tends to one. There are different porous plasticity models to account for the mechanism of compaction and void growth.
The Shima–Oyane Criterion
Shima and Oyane (Ref. 9) proposed a yield surface for modeling the compaction of porous metallic structures fabricated by sintering. The criterion can be 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 equivalent 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 equivalent stress, σys0 is the initial yield stress, pm is the pressure, and ρrel is the relative density. The material parameters α, γ, and m are obtained from curve fitting experimental data. Typical material parameter values for copper aggregates are α = 6.2, γ = 1.03, and m = 5.
The Gurson Criterion
The Gurson criterion (Ref. 10) 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 The Shima–Oyane Criterion, but it is defined in terms of the hyperbolic cosine function. The plastic potential for the Gurson criterion reads
here, σe is the equivalent stress, σys0 is the initial yield stress, pm is the pressure, and ϕ is the porosity.
The Gurson–Tvergaard–Needleman Criterion
Tvergaard and Needleman modified The Gurson Criterion for porous plasticity to include parameters to better fit experimental data (Ref. 11-12). The resulting criterion is called the Gurson–Tvergaard–Needleman (GTN) criterion in the literature. The plastic potential for the GTN criterion reads
here, σe is the equivalent stress, σys0 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. It is often given using a bilinear definition
where ϕ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 the value of failure, the effective porosity takes a maximum value of ϕm and the porous material loses the capacity to carry stresses. The maximum porosity value is derived from other parameters
Since typical values for the parameters are q3 = q12, the maximum porosity value is given by ϕm = 1/q1.
A similar definition for ϕe that gives a smoother response as the material reaches failure can was suggested in Ref. 14. It is based on a modification such that the effective void volume fraction reaches its maxim value asymptotically
where ϕi = 0.75(ϕf − ϕc) + ϕc, and
The Fleck–Kuhn–McMeeking Criterion
The Fleck–Kuhn–McMeeking criterion (Ref. 15), also called the 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 the FKM criterion reads
here, σe is the equivalent stress and pm is the pressure. The flow strength of the material under hydrostatic loading, pf, is computed from
here, σys0 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.
The FKM–GTN Criterion
The FKM–GTN criterion is a combination of The Fleck–Kuhn–McMeeking Criterion and The Gurson–Tvergaard–Needleman Criterion, intended to cover a wider range of porosities (Ref. 1617). The GTN model is used for low void volume fractions (porosity lower than 10%), 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.
Capped Drucker–Prager
See the sections The Drucker–Prager Criterion and Elliptic Cap for details.
Void Nucleation and Growth
For the porous plasticity criteria, the change in relative density is by default computed from the change in plastic volumetric strain
Since the relative density is related to the porosity ϕ by ρrel = 1 − ϕ, the change in porosity is also controlled by the change in plastic volumetric strain
and the change in volumetric plastic strain is given by the porous plasticity model.
Nucleation and shear growth are other mechanism that can trigger an increase of porosity. The increment in porosity based on growth nucleation is given by
for
here, εN is the mean strain for nucleation, fN is the void volume fraction for nucleating particles, and sN is the standard deviation. Typical values for these parameters are εN = 0.04, sN = 0.1, and fN = 0.04. It is assumed that nucleation appears only in tension, and that there is no nucleation in compression.
The other possible mechanism to change the porosity is the so-called void growth in shear
here, kw is a material parameter, ϕ is the current porosity, nD is a deviatoric tensor coaxial to the stress tensor, and is the plastic strain rate, which depends on the porous plasticity model. The weight w is computed from stress invariants as
where θ is the Lode angle. The weight variable has a value of w = 0 at the compressive and tensile meridians (θ = 0 and θ = π/3), and it attains its maximum w = 1 for θ = π/6.
Isotropic Hardening
Porous plasticity implements different kinds of isotropic hardening models to describe the hardening of the porous matrix. When applying isotropic hardening, the average flow stress
now depends on the initial yield stress σys0, the hardening function σh, and the equivalent plastic strain in the porous matrix εpm.
Different isotropic hardening models are implemented for porous plasticity models, which are described for elastoplastic materials:
Perfect Plasticity (no isotropic hardening)
It is also possible to include a power law relation between the equivalent plastic strain in the porous matrix εpm, and the flow stress (yield level) σfm. For uniaxial loading, the strain stress relation is written as (ref. Ref. 12)
for
where σys0 is the initial yield stress, n is the hardening exponent, and the Young’s modulus E is taken from the elastic material properties. By writing the onset of plasticity as ε0 = σys0/E, and noting that ε = ε0 + εpm for σfm > σys0, this reads
for
which is equivalent to Swift isotropic hardening
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-55 is numerically solved with the so-called exponential mapping technique
where and . The purpose of the exponential mapping is to ensure volume preservation of the discretized equations. For plasticity models that are not volume preserving, for example soil plasticity, the computationally leaner backward Euler discretization is used also for large plastic strains.
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. 3):
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).
It is possible to specify the maximum number of iterations and the relative tolerance used to solve the plastic flow rule. Under the Advanced section enter the following settings:
Maximum number of local iterations. To determine the maximum number of iteration when solving the local plasticity equations. The default value is 25 iterations.
Relative tolerance. To check the convergence of the local plasticity equations based on the step size. The final tolerance is computed based on the current solution of the local variable and the entered value. The default value is 1e-6.
To display the Advanced section, click the Show More Options button () and select Advanced Physics Options in the Show More Options dialog box.
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.
Nonlocal Plasticity
The nonlocal plasticity model implemented in COMSOL Multiphysics is based on the so-called implicit gradient method suggested in Ref. 23, and it incorporates some of the generalizations introduced by the micromorphic theory (Ref. 24). The theory assumes that the free energy of the elastoplastic material not only depends on the macroscopic displacement u and macrostate variables, such as the equivalent plastic strain εpe, but also on some micromotion ϕ and its gradient. If the micromotion variable is identified as a nonlocal version of the equivalent plastic strain, εpe,nl, the free energy potential is written as
(3-65)
here Ws is the elastic free energy, Hnl is a model parameter that determines the strength of the coupling between the macro- and micromotions, and lint is the length scale that determines the influence radius of the interaction. For simplicity, Equation 3-65 assumes linear isotropic hardening, but the same concept is applicable for more general isotropic hardening laws and can for example also be extended with kinematic hardening.
Taking the variation of Equation 3-65 with respect to state variables εpe, εpe,nl, and ∇εpe,nl leads to the following set of equations related to the plasticity model:
here, Ω is the domain of the elastoplastic problem. Variation of Equation 3-65 with respect to the remaining state variables leads to the standard elasticity equations. In the nonlocal (or regularized) plasticity problem, εpe,nl is an additional dependent variable that is solved as part of the global problem. The structure of the local plasticity problem is not changed and the numerical solution remains the same as discussed above, but the hardening law now depends on εpe,nl.
For porous plasticity, the same equation applies, but εpe and εpe,nl are replaced by the equivalent plastic strains in the matrix material, εpm and εpm,nl. It is also possible to add nonlocal effects to the void volume fraction ϕ. The porous plasticity model is then extended with the following set of equations
where ϕnl is the nonlocal void volume fraction. In contrast to the equivalent plastic strain, it is assumed that there is a perfect coupling between ϕ and ϕnl, and that the nonlocal void volume fraction ϕnl directly replaces its local counterpart in the model equations.
Resetting Elastoplastic Variables
For some industrial applications, plastic strains are released due to manufacturing processes such as welding or stress relief annealing. To model these type of phenomena, it is possible to reset the internal plasticity variables to zero, or to prescribe them to user defined expressions or values.
Add a Set Variables subnode to the Plasticity, Porous Plasticity, Soil Plasticity, or Elastoplastic Soil Models to set up the plastic strain tensor, plastic deformation gradient, the equivalent plastic strain, and other internal elastoplastic variables.