Plasticity: Mathematical Modelling of Plasticity
In this section, the physical behaviour described previously will be modelled theoretically. Most of the equations in this section were initially developed by “curve fitting” to the behaviour and thus, mathematical idealizations are considered to simplify the relationships. The classical mathematical modelling of plasticity is formulated “incrementally” as shown in the following sections.
Yield Stress and Codes of Practice:
The yielding and strain hardening regions of the metal stress strain curve shown in Figure 1 are characterized by an increase in the “yield stress” with an increase in the amount of plastic deformation. While some metals, for example mild steel, exhibit a distinct “yield stress”, others have a more round stress strain curve with a gradual decrease in the slope of the curve between the stress and the strain. The yield stress in general is history dependent, i.e., it depends on the history of loading of the specimen and is not a unique function of the current state of stress. Codes of practice often prescribe three specified minimum values for the properties of a steel to be fit for a specific application. These are normally the “yield stress” the “ultimate stress” and the “elongation” corresponding to the ultimate stress. The ultimate stress and the elongation corresponding to the ultimate stress can easily be obtained from any of the typical curves shown in Figure 2, however, obtaining the exact elastic limit could be practically very difficult, and therefore, codes of practice often define the “yield stress” to be the stress corresponding to a specific value of strain (e.g. 0.005) or to the stress corresponding to the intersection between the 0.002 offset of the linear stress strain portion (Figure 4). The two criteria do not necessarily predict the same value for the yield stress .

Figure 4. Definition of yield stress or given by different codes of practice
Stress Strain Idealization:
The relationship between the “true” stress versus the “true” or logarithmic strain is idealized to simplify the mathematical model. Points A and B are assumed to coincide and the loading and unloading in the yield and strain hardening region is assumed to follow a linear relationship as shown in Figure 5a. In other applications, especially those involving hand calculations or beam analysis, it is customary to assume an elastic-perfectly plastic stress strain curve, where the yield stress and the ultimate stress coincide (Figure 5b).

Figure 5. Different idealizations for the true stress versus true strain curves of metals: a) Unloading follows the dotted lines, b) elastic-perfectly plastic assumption with unloading following the dotted lines
Strain Decomposition and Incremental Behaviour in a Uniaxial Stress State:
Upon removal of the applied load, metals exhibit what is referred to as permanent or plastic deformation. The true longitudinal strain can thus be decomposed into two components: an elastic component corresponding to the applied stress, and a permanent plastic deformation. In a uniaxial stress state where the force is applied along the direction of the basis vector ,
corresponds to the true stress or the Cauchy stress component and
corresponds to the true longitudinal strain or logarithmic longitudinal strain along the direction of the basis vector
. The strain component can be decomposed into two components: an elastic component
and a plastic component
(Figure 6a). The elastic component is related to the stress using the linear elastic constitutive laws:
The plastic strain component can be obtained from the experimental results as follows:
An explicit relationship describing the uniaxial stress component and the plastic strain component
is often required (Figure 6b) for computations and thus a curve fitting model is often used. One of the very successful models is the Ramberg-Osgood model that fits a power law to the relationship between stress and plastic strain by introducing two material constants
and
as follows:
If an increment in the load is such that is decreased, then the elastic strain component is decreased while the plastic strain component stays the same. If an increment in the load is such that
reaches
, then an increase beyond
will cause an increase in the plastic strain
. i.e., unloading and reloading would follow vertical lines in Figure 6b. Thus, the above relationship in fact describes the magnitude of the yield stress at a certain value of a plastic strain parameter
and is more precisely written in the form:
The plastic strain parameter is a measure of the history of plastic deformation such that
and is always increasing. The plastic strain parameter is also referred to as the equivalent plastic strain. The overbar is used to differentiate between the plastic strain parameter and the plastic strain tensor
.

Figure 6. a) Decomposition of the strain component into an elastic and a plastic component. b) the relationship between the plastic strain component and the applied stress
The remaining elastic strain components

In the above idealization, three functions were implicitly presented to fully describe the constitutive law of the plasticity of metals. These are the yield function, the hardening rule and the flow rule. As plastic deformation is dependent on the history of the loading, these functions can be conveniently written in terms of one parameter that represents the loading history. This one parameter is chosen to be the value of The following section describes these three rules.
Yield Function:
The yield function is the function that describes whether the material is about to exhibit plastic deformation or not. In a uniaxial state of stress, if is less than
then, the behaviour is elastic. Once
reaches
then the material is about to exhibit plastic deformation.
can never be higher than
since as
increases,
increases as well. If the stress matrix of the uniaxial stress state is denoted by
, then the yield function can be defined as:
If , then the material is exhibiting elastic deformations, while if
then, the material is about to exhibit plastic strain. As the material exhibits plastic strain, the change in the applied stress has to be such that the new stress state mathematically lies on the yield surface, i.e.,
In other words, the rate of change of the stress components is related to the rate of change of the equivalent plastic strain
. The above relationship is called the consistency condition. Traditionally, the consistency condition is written in the following incremental form:
The incremental form relates the increment in the increase in the applied stress to the incremental change of the yield stress of the material, or the increment in the applied stress to the increment in the plastic strain parameter. The incremental form represents a numerical integration scheme by which the increments in can be calculated as the applied stress is changed. The backward Euler numerical integration scheme is usually employed as it results in an approximation that is more accurate than the traditional forward Euler numerical integration scheme.
Hardening Rule:
The hardening rule is the function that gives the value of the yield stress as a function of the history of the loading. A typical relationship obtained from experiments is shown in Figure 6b and can conveniently be represented by the Ramberg-Osgood approximation:
This relationship is termed the hardening rule. Given the current amount of the accumulated plastic strain parameter in the material, the value of the yield stress
is obtained.
Flow Rule:
The flow rule is the function that dictates the amount of plastic strain that will accumulate during an increment in the loading. For a uniaxial state of stress the rate of change of the plastic strain aligned with the direction of the applied stress is equal to the rate of change of the plastic strain parameter (the equivalent plastic strain) . The rate of change of the applied stress is equal to the rate of change of the yield stress, i.e,.
. Thus, the rates of change of
and
are related as follows:
In addition, the isochoric nature of the plastic deformation leads to the following relationship between the rates of change of the plastic strain components:
The above relationships can be written in incremental form such that and
, then the increments in the plastic strain components are given by:
The total value of the plastic strain components at the end of each increment is equal to the sum of all the plastic strain components increments developed throughout the loading. The backward Euler numerical integration scheme is usually employed for better accuracy.
Strain Decomposition and Incremental Behaviour in Multi-axial Stress States with Isotropic Hardening (No Bauschinger Effect):
The three functions that were presented in the previous section can be extended to cover a multi-axial stress state.
Multi-axial von Mises Yield Function:
The introduction of the multi-axial yield function dates back to the early years of the twentieth century with various researchers conducting experiments on the yielding of metals. As it was observed that metals are, in general, linear elastic isotropic materials and that yielding occurs independently of the hydrostatic stress, it was proposed that the material would start yielding once the deviatoric component of the elastic strain energy stored inside the material reaches a critical state. The expression for the deviatoric strain energy is equal to:
Where, the von Mises stress, denoted by , is an invariant of the deviatoric stress tensor. The “critical” value at which the material starts to yield can be obtained from a uniaxial state of stress experiment. In the uniaxial state of stress,
is equal to
and yielding will occur once
. Thus, the critical strain energy at the onset of yield is equal to:
Thus, in a state of stress described by a Cauchy stress matrix , the yield function can be defined as follows:
In this relationship, describes the state of the stress inside the material while
gives the maximum value that
can attain before yielding.
is presented as
in the stress based failure criteria section. Similar to the uniaxial case, if
, then the material is exhibiting elastic deformations while if
then, the material is about to exhibit plastic strain. The three dimensional version of the consistency condition presented above has the following form:
In other words, the rates of change of the stress components are related to the rate of change of the equivalent plastic strain
. The incremental form of the consistency condition is such that
. This implies that as the stress changes, the yield function should not change and should also be equal to 0. Thus:
The incremental form of the consistency condition in a multi-axial state of stress relates the increment in the increase in the applied stress represented by to the incremental change of the yield stress of the material, or the increment in the applied stress components
to the increment in the accumulated plastic strain parameter
. The backward Euler integration scheme is recommended for better accuracy of the results.
As shown in the stress based failure criteria section, if a coordinate system is chosen such that stress matrix is diagonal and the only possible nonzero components are ,
, and
, then the von Mises yield function is represented in the stress space by a cylinder whose axis is the line representing the hydrostatic stress states (
). To show this, we will investigate the the relationship:
after a coordinate transformation into the new coordinate system ,
, and
. Obviously,
represents the line of the hydrostatic stress while
and
are perpendicular to that line. The coordinate transformation matrix into the new coordinate system has the form:
A point with coordinates will have coordinates
in the new coordinate system such that:
By replacing ,
, and
in the yield surface relationship with the coordinates
,
, and
:
Suggesting that in this new coordinate system, the relationship represents a cylinder whose axis is the vector and whose radius is
.
View Mathematica Code
Q = {{1/Sqrt[3], 1/Sqrt[3], 1/Sqrt[3]}, {1/Sqrt[2], -1/Sqrt[2], 0}, {1/Sqrt[6], 1/Sqrt[6], -2/Sqrt[6]}}; a = {a1, a2, a3}; s = Transpose[Q].a; s // MatrixForm FullSimplify[(s[[1]] - s[[2]])^2 + (s[[2]] - s[[3]])^2 + (s[[1]] - s[[3]])^2]
Another method to find the radius of the cylinder is by realizing that the stress state represented by the vector in Figure 7 with the components
,
lies on the surface of the cylinder. The unit vector representing the hydrostatic direction can be obtained by normalizing the vector
as follows:
The radius of the cylinder is given by the norm of the vector as follows:

Figure 7. von Mises yield function in a multi-axial stress state represented in the stress space of the principal stresses
Multi-axial Tresca Yield Function:
Another yield function that predicts similar results is that given by the maximum shear stress as described in the stress based failure criteria section. The Tresca Yield function predicts that a material will yield when the maximum shear stress reaches a critical value :
The “critical” value can be obtained from a uniaxial state of stress experiment. In the uniaxial state of stress,
is equal to
. Thus, in a state of stress described by a Cauchy stress matrix
, the Tresca yield function can be defined as follows:
This yield function is represented graphically in the stress space of the principal stresses by a hexagon. Computational algorithms that use the von Mises yield criterion and those that use the Tresca yield criterion predict almost the same behaviour. While the von Mises yield criterion has a smoother curvature, the Tresca yield criterion may be easier for hand calculations. In general, when viewed on the plane perpendicular to the hydrostatic stress, the differences are practically negligible (Figure 8).

Figure 8. View perpendicular to the hydrostatic stress showing the difference between the von Mises (circle) and the Tresca (hexagon) yield surfaces.
Isotropic Hardening Rule:
In a multi-axial state of loading, the isotropic hardening rule predicts that the radius of the von Mises circle will change according to the relationship given between and
. Usually, additional plastic strain (or plastic energy dissipation) applied to the material results in an increase in the radius of the von Mises circle (Figure 9).

Figure 9. Isotropic hardening represented by a change (usually an increase) in the radius of the circle with the increase of the plastic strain parameter
Multi-axial Flow Rule:
The last ingredient of incremental plasticity is to find the rate of change of the plastic strain components associated with the rate of change of the stress components that is causing plastic flow. Experimental evidence suggest that the rates of the change of the plastic strain components are directly proportional to the applied deviatoric stress components as follows:
i.e.,
The incremental form predicts that the increments in the plastic strain components are directly proportional to the deviatoric stress components in a material as follows:
i.e.,
This observation conforms with the isochoric nature of the plastic deformation (why?). Notice that there is an inherent difference between the calculations of the plastic strain components and the elastic strain components. The increments of the elastic strain components are directly proportional to the increments in the stress components, while the plastic strain components, as observed from experiments, are directly proportional to the components of the deviatoric stress tensor.
In order to be able to compute the plastic strain components associated with a plastic flow instant, the constant of proportionality or
needs to be computed and then multiplied by the respective deviatoric stress component. This constant of proportionality can be obtained from energy considerations by assuming that the energy that goes into plastic deformation in a general state of stress is equal to the energy that goes into plastic deformation in a uniaxial state of stress. In a uniaxial state of the stress, the rate of the energy that is dissipated during the plastification of the material is equal to
which is equal to
. The energy dissipated during the plasticification of the material in a general state of stress is equal to
. Assuming that the energy dissipated in the plastification process in a general state of stress is equal to that in a uniaxial state of stress results in the following:
Replacing the plastic strain components with the deviatoric stress components and the constant of proportionality results in the following expression:
However, since and since
, the stress components
can be replaced with the deviatoric components as follows:
Note that by definition, the von Mises stress is equal to:
Thus, the relationship between and
is as follows:
Therefore, the plastic strain components can be calculated by integrating the following differential equations functions in :
The equivalent plastic strain is always positive and is a measure of the history of loading of the material. The term “equivalent plastic strain” is similar to the term “equivalent stress” that is used to denote . In addition, the rate of the “equivalent plastic strain” can be related to the rates of the plastic strain components with a relationship similar in form to that between
and the deviatoric stress components:
Therefore:
The above relationship can be written in the following incremental form:
Associated Flow Rule:
The assumption that the incremental plastic strain components are directly proportional to the deviatoric stress components resulted in the final expression (as shown above):
It has been suggested in the early twentieth century that, similar to the existence of an “Elastic Potential Energy Function” from which the stress components and the elastic strain components are derived, the plastic strain components can also be derived from a “Plastic Potential Function” where:
Careful investigation of the above expression in comparison with the previous expression for the incremental plastic strain components suggests that can be taken to be the yield function
. In this case:
Bonus Exercise: Show the following relation:
The term “associated flow rule” is used to signify that the plastic strain components are derived from the yield function and thus, the flow rule is associated with the yield function. Associated flow rules are usually used for metal plasticity, while for plasticity of other materials, the incremental plastic strain components could be derived from different functions that are not associated with the yield function. In the latter case, the flow rule is termed: “non-associated flow rule”. Associated flow rules in metals that are derived from the von Mises or the Tresca yield surfaces predict zero volumetric plastic strain, however, in soils and many other materials, plasticity is associated with a change in volume and different “plastic potential” functions are utilized.
The Normality Rule:
Consider a smooth function such that
. With 9 independent components,
can also be viewed as an element of a 9 dimensional vector space. The function
can then be considered as a smooth scalar valued function
. The gradient of
with respect to its variables is denoted by
and is defined as a vector field in
with components
(See the section on the gradient of scalar fields). The direction pointing towards constant
in
can be obtained as the unit vector
along which the directional derivative is equal to zero:
Thus, the vector is normal to the vectors pointing in the direction of constant
. If
is taken to be the yield function defined on the nine dimensional Euclidean vector space of the nine stress components (termed the stress space), then the vector
is perpendicular to the direction of constant
in the stress space. Since
is constant (
) on the yield surface, then
is normal to the yield surface. The incremental plastic strain components can also be considered to be elements of a nine dimensional Euclidean vector space and the associated flow rule asserts that the nine dimensional vector
with components
is linearly dependent on (parallel to)
as:
Thus, the nine dimensional vector can also be viewed as normal (perpendicular) to the yield surface and thus, the associated flow rule leads to the normality rule. In a three dimensional vector space of the principal stresses of
, the components of
are perpendicular to the yield surface (Figure 10).
See Example 1 in the examples and exercises page for a detailed example of an isotropic hardening material model.

Figure 10. Schematic of the normality rule in the three dimensional vector space of the principal stresses.
Strain Decomposition and Incremental Behaviour in Multi-axial Stress State with Combined Kinematic and Isotropic Hardening (Bauschinger Effect):
In order to include the Bauschinger effect, various modifications to the isotropic hardening models presented in the previous section have been proposed throughout the development of the theory. The most successful and generalized form is that of the combined nonlinear kinematic and isotropic hardening model. In this section, the general form for this material model is presented.
Uniaxial Stress State Representation:
In a uniaxial state of stress, modelling the Bauschinger effect is implemented by assuming that the difference between the yield stress in tension and in compression is due to the shifting of the yield surface by a value that is equal to , which is termed the backstress. The yield stress
is assumed to be the sum of two components:
which is the yield stress with respect to the shifted centre of the yield surface and
which describes the shifting of the yield surface. Both
and
are assumed to be functions of the history of the deformation represented by
and the values of the stress increment. Figure 11 shows the schematics of three different models that can be utilized to model the Bauschinger effect. In the first model (Figure 11a), referred to as the linear kinematic model, the yield stress with respect to the shifted centre
is assumed to be constant while the shifting of the centre of the yield surface is assumed to be linear. This model can be enhanced by assuming that the
varies with the equivalent plastic strain in a nonlinear fashion (Figure 11b). It is also possible to assume that the backstress follows a nonlinear relationship leading to the general nonlinear combined kinematic and isotropic hardening model, as shown in Figure 11c. Obviously, in order to calibrate such material model, numerous experiments have to be conducted to fit material curves for
and for
.

Figure 11. Three different models for the decomposition of the yield stress into a yield stress
with respect to the centre of the yield surface and a backstress
. (a) Kinematic hardening models, (b) combined linear kinematic hardening and nonlinear isotropic hardening model, (c) nonlinear kinematic and isotropic hardening model.
The most successful constitutive law model that describes the evolution of the backstress is the Armstrong-Frederick evolution law which in a uniaxial state of stress has the following form:
or in incremental form:
Where and
are material constants. An earlier version of the above evolution law contained only the first term and was termed “Ziegler’s law”. Ziegler’s law had a limitation that
would be unbounded and keeps increasing with the increase of the plastic strain
. The clever introduction of the term
forces the value of
to be bounded since integration of the above formula leads to:
Thus, at higher levels of the value of
is bounded by the ratio of the material constants
.
Multi-axial Stress State Implementation:
In order to implement the above relationships in a multi-axial stress state, the backstress is assumed to be a tensor with components
. The yield function in this case is given as follows:
where is the von Mises stress but with the difference
substituted instead of the stress tensor
. Figure 12 shows the graphical representation of the yield function in the stress space of the principal stresses. If
is the deviatoric backstress tensor, then, the yield function can be written as:
The Armstrong-Frederick material model describing the evolution of the backstress in a multi-axial stress state is given by the following form:
Or in incremental form:
The bounding limit of is given by the value
and thus, the centre of the yield surface in the stress space is bounded by a cylinder of radius
, while the yield surface is bounded by the cylinder with a radius equal to
(Figure 12).
The associated flow rule can still be used for this material model where the plastic strain components are calculated as follows:
Or in incremental form:
The increment in the plastic strain parameter can be obtained from the consistency condition as follows:
or in incremental form:
With
Bonus exercise:
Deduce the above relationship for
See Example 2 in the examples and exercises page for a detailed example of a nonlinear kinematic hardening material model.

Figure 12. (a) Schematic showing the shifted yield surface in the three dimensional stress space of the principal stresses. The inner dotted circle represents the bounding limits of the centre of the yield surface, while the outer dotted circle represents the bounding limits of the maximum stress. (b) Uniaxial representation