How to Implement Elastoplasticity in a Model Using External Materials

Mats Danielsson November 2, 2017
Share this on Facebook Share this on Twitter Share this on LinkedIn

In structural mechanics, there may be situations when you want to implement your own material model. The COMSOL Multiphysics® software gives you the option to program your own material model in C code. The compiled code can then be called from the program using the External Material feature. Here, we demonstrate how to implement an external material model and then use it in an example analysis.

The External Material Functionality

The Nonlinear Structural Materials Module, an add-on product to COMSOL Multiphysics®, provides a plethora of material models, including models for hyperelasticity, isotropic/kinematic hardening plasticity, viscoelasticity, creep, porous plasticity, soils, and more. These material models cover a vast majority of engineering problems within structural analysis.

However, in some situations, the mechanical behavior of a material is not readily expressed in terms of an existing material model. For instance, suppose you developed a specialized material model for a certain alloy and want to use it to solve a large structural mechanics boundary value problem in COMSOL Multiphysics. What do you do?

As a matter of fact, there are three different ways in which you can define your own material:

  1. Many of the material models have a User Defined option. As an example, you can create your own hyperelastic material by typing in an expression for the strain energy density as a function of the strains directly in the GUI.
  2. The next level, which is useful when there is no suitable material model that you can tune through a User Defined option, is to express your material model as a set of separate PDEs or ODEs, and enter the appropriate expressions using one of the mathematical interfaces.
  3. Finally, you can program your own material model in C, using the External Material feature. If the material model is expressed as an algorithm, rather than as a set of equations, this would be the preferred method.

The implementation of a material model as an external DLL can seem like a complex endeavor, but this blog post demonstrates how to implement an elastoplastic material model in COMSOL Multiphysics using hands-on steps that you can follow.

Implementing an Elastoplastic Material Model in COMSOL Multiphysics®

As a starting point, we need to decide on a material model to implement. We choose an isotropic linear-elastic material with isotropic hardening. This is a simple plasticity model that already exists in COMSOL Multiphysics, but it serves nicely to convey some key points.

First, let’s go over some assumptions, definitions, and nomenclature:

  • Strains are additive and assumed small
  • The elastic response is linear and isotropic, and defined by Young’s modulus E and Poisson’s ratio ν (or equivalently by the bulk and shear moduli, K and G)
  • The plastic response is isotropic with an initial yield stress of \sigma_{ys0}
  • The yield stress of the material is a function of effective plastic strain, \sigma_{ys}=\sigma_{ys0}+\sigma_h\left(\epsilon_{pe}\right), with the initial value \sigma_{h}\left(0\right)=0 (the hardening function)
  • Plastic flow follows J_2 theory
  • An increment in effective plastic strain is given by \Delta\epsilon_{pe}=\left(\frac{2}{3}\Delta\bm\epsilon_p:\Delta\bm\epsilon_p\right)^{1/2}

Diagrams illustrating the uniaxial stress-strain curve and yield surface for an elastoplastic material model.
The example material model: Uniaxial stress-strain curve and yield surface in principal stress space.

Now, let’s discuss approaches for implementing a material model as an external material. There are several different ways of calling user-coded routines for external materials, which we refer to as sockets.

Approach 1: General Stress-Strain Relation

We can use the General stress-strain relation socket of the external material to define a complete material model that includes (possibly) both elastic and inelastic strain contributions. This is the most general of the two modeling approaches discussed here. When we use the General stress-strain relation socket, we are faced with two tasks:

  1. Computation of the stress tensor \bm{\sigma}
  2. Computation of the Jacobian of stress \bm{\sigma} with respect to strain \bm{\epsilon}

Approach 2: Inelastic Residual Strain

We can also use the Inelastic residual strain socket to define a description of an inelastic strain contribution \bm\epsilon_{inel} to the overall material model. An example of this would be if we wanted to add our own creep strain term to the built-in linear elastic material. The Inelastic residual strain socket assumes an additive decomposition of the total (Green-Lagrange) strain into elastic and inelastic parts. Thus, this is an adequate assumption when strains are of the order < 10%. When we use the Inelastic residual strain socket of the external material model, we are faced with two tasks:

  1. Computation of the inelastic strain tensor \bm{\epsilon}_{inel}
  2. Computation of a Jacobian of inelastic strain \bm{\epsilon}_{inel} with respect to strain \bm{\epsilon}

Two related External Material sockets are the General stress-deformation relation and the Inelastic residual deformation. These are more general versions of those discussed above. Instead of defining the deformation in terms of the Green-Lagrange strain tensor, the deformation gradient is provided. Many large-strain elastoplastic material models use a multiplicative decomposition of the deformation gradient into elastic and plastic parts. In these situations, you would likely want to use one of these sockets instead.

Tip: We link to the source file and model file at the bottom of this blog post.

Computation of the Stress Tensor: Stress Update Algorithm

The complexity of computing the stress tensor varies significantly between material models. In practice, the computation of the stress tensor often needs to be formulated as an algorithm. This is often called a stress update algorithm in literature. In essence, the objective of a stress update algorithm for a material model is to compute the stresses, knowing:

  • The total strains corresponding to the current estimate of the displacement field
  • The material state of the previously converged increment

These quantities are provided to the external material as input.

The term “material state” represents any solution-dependent internal variables that are required to describe the material. Examples of such variables are plastic strain tensor components, current yield stress, back stress tensor components, damage parameters, effective plastic strain, etc. The choice of such state variables will depend on the material model. We must ensure that the material state is properly initialized at the start of the analysis, and that it is updated at the end of the increment.

We first need to investigate if there is plastic flow occurring during the increment. We do this by assuming that the elastic strain is equal to the total strain of the current increment, less the (deviatoric) plastic strain of the previously converged increment, ^{old}\bm\epsilon_{p}. This assumption would hold true if there was, indeed, no plastic flow during the increment. The deviatoric stress tensor that is computed this way is aptly called a trial stress deviator and is given by

\mathrm{dev}\left(\bm\sigma_{trial}\right) = 2G\,\mathrm{dev}\left(\bm\epsilon-^{old}\bm\epsilon_{p}\right)

with an effective (von Mises) value \sigma_{trial,e}.

The effective trial stress is compared to the yield stress of the material by assuming that there is no plastic flow during the increment. The yield stress corresponding to the previously converged increment is given by

^{old}\sigma_{ys} = \sigma_{ys0}+\sigma_h\left(^{old}\epsilon_{pe}\right)

Notice that the left-superscripted ^{old}\bm\epsilon_p and ^{old}\epsilon_{pe} in Steps 1 and 2 represent the material state of the previously converged increment, as we discussed earlier.

Now we check if the trial stress causes plastic flow. For another way of expressing this, if the trial stress is inside the yield surface, the response will be purely elastic during the increment. If not, plastic flow will result. The check is performed using the yield condition:

f=\left\{\begin{array}{c l}\sigma_{trial,e}-^{old}\sigma_{ys}<0 & \mathrm{Elastic}\\ \sigma_{trial,e}-^{old}\sigma_{ys}\ge0 & \mathrm{Elastoplastic}\end{array}\right.

Diagrams illustrating the yield condition for elastic and elastoplastic computations.
Check of the yield condition to determine elastic or elastoplastic computation.

The stress update algorithm now necessarily branches off into either a purely elastic computation or an elastoplastic computation. We will follow each of these branches, starting with the purely elastic branch.

Elastic Computation

Because we determined that there is no plastic flow during the increment, the trial stress deviator is, in fact, identical to the stress deviator, and the update of the plastic strain tensor and the effective plastic strain is trivial.

\left\{\begin{array}{lcl}\bm\sigma &=& \mathrm{dev}\left(\bm\sigma_{trial}\right) + K\mathrm{trace}\left(\bm\epsilon\right)\mathbf{1}\\\bm\epsilon_p &=& ^{old}\bm\epsilon_p\\\epsilon_{pe} &=& ^{old}\epsilon_{pe}\end{array}\right.

We can directly return the pure elastic stress-strain relation as the Jacobian.

Elastoplastic Computation

The objective of the elastoplastic branch of the stress update algorithm is to compute the stress deviator and update the plastic strains. We begin by again expressing the stress deviator, now knowing that plastic flow takes place during the increment:

\mathrm{dev}\left(\bm\sigma\right) =2G\,\mathrm{dev}\left(\bm\epsilon-^{old}\bm\epsilon_{p}-\Delta\bm\epsilon_p\right)=\mathrm{dev}\left(\bm\sigma_{trial}\right)- 2G\Delta{\lambda}\,\mathrm{dev}\left(\bm{\sigma}\right)


\mathrm{dev}\left(\bm\sigma\right) = \frac{1}{\left(1+2G\Delta\lambda\right)}\mathrm{dev}\left(\bm\sigma_{trial}\right)

In the above equation, we used a discrete form for the flow rule that states that an increment in plastic strain is proportional to the stress deviator through a so-called plastic multiplier \Delta\lambda. Let’s stop for a moment and consider a graphical representation of this equation for stress:

A diagram showing the correction of the trial stress deviator.
Graphical representation of the correction of the trial stress deviator.

If we compute a trial stress deviator that lies outside the yield surface, we need to make a correction so that the stress deviator is returned to the to-be-determined yield surface. The plastic multiplier determines the exact amount by which the trial stress deviator should be scaled back to give the correct stress deviator. If we compute the plastic multiplier, it is straightforward to then compute the stress deviator and the plastic strain increment.

The key steps are to:

  1. Transform the equation for stress into a governing scalar equation
  2. Express all the unknowns in a single parameter
  3. Solve the scalar equation for this parameter

We can relate the plastic multiplier to the effective plastic strain increment \Delta\epsilon_{pe} using the flow rule and then transform the equation for stress into a governing scalar equation:

\sigma_{ys0}+\sigma_h\left(^{old}\epsilon_{pe} + \Delta\epsilon_{pe}\right)+3G\Delta\epsilon_{pe}-\sigma_{trial,e}=0.

This is in general a nonlinear equation, and we need to solve it using a suitable iterative scheme. We are now ready to compute the stress tensor, the plastic strain tensor, and the effective plastic strain.

\left\{ \begin{array}{ccl}\bm\sigma &=& \frac{1}{1+2G\Delta\lambda}\mathrm{dev}\left(\bm\sigma_{trial}\right) + K\mathrm{trace}\left(\bm\epsilon\right)\mathbf{1}\\ \\\bm\epsilon_p &=& ^{old}\bm\epsilon_p + \Delta\bm\epsilon_p\\\epsilon_{pe} &=& ^{old}\epsilon_{pe} + \Delta\epsilon_{pe} \end{array}\right.

The updated plastic strain tensor and effective plastic strain are stored as state variables.

Computing the Jacobian of the Material

We computed the stresses and updated the material state (the state variables) for our material model. Now, we turn our attention to the Jacobian computation. The Jacobian has other names in literature, such as tangent stiffness, tangent modulus, or tangent operator. In the stress update algorithm, we express the deviatoric and hydrostatic parts of the stress tensor as:

\begin{array}{rcl} \mathrm{dev}\left(\bm\sigma\right) &=& \mathrm{dev}\left(\bm\sigma_{trial}\right)-2G\,\Delta\bm\epsilon_p \\\mathrm{trace}\left(\bm\sigma\right)&=&3K\mathrm{trace}\left(\bm\epsilon\right)\end{array}.

The Jacobian that we want to compute is the derivative of the Second Piola-Kirchhoff stress tensor with respect to the Green-Lagrange strain tensor. For our example material, we assume that strains are small. This means that we do not need to distinguish between various measures of stresses and strains, because they are indistinguishable in the small-strain limit. The derivative \mathcal{C} of stress with respect to strain is written as


If we use the equations for the deviatoric and hydrostatic stress and the definition of the trial stress, we can express the Jacobian in the following way:

\mathrm{d}\bm\sigma= \underbrace{ \left[ \underbrace{2G\left(\mathcal{I}-\frac{1}{3}\mathbf{1}\otimes\mathbf{1}\right) + K\left(\mathbf{1}\otimes\mathbf{1}\right)}_{\mathcal{C}^{e}} – 2G\frac{\mathrm{d}\bm\epsilon_p}{\mathrm{d}\bm\epsilon} \right]}_{\mathcal{C}}:\mathrm{d}\bm{\epsilon}

Note that we replaced the increment of the plastic strain tensor by the total plastic strain tensor in the expression above. Their derivatives with respect to strain are the same, by virtue of the additive update of the plastic strains. Recall that our two modeling approaches require differently defined Jacobians. We see immediately how they are related. In the General stress-strain relation, the Jacobian is given by the full expression above. In the Inelastic residual strain, the Jacobian is given by one term in the expression, namely:

\frac{ \mathrm{d}\bm\epsilon_{inel}}{\mathrm{d}\epsilon}=\frac{\mathrm{d}\bm\epsilon_p}{\mathrm{d}\epsilon}

The term \mathcal{C}^e is the elastic Jacobian. For a purely elastic computation, the total Jacobian of the General stress-strain relation equals this quantity, while the Jacobian of the Inelastic residual strain in this case is zero.

If we use the flow rule and the chain rule for differentiation, we arrive at the following expression:

\frac{\mathrm{d}\bm\epsilon_p}{\mathrm{d}\bm\epsilon} = \Delta\lambda\frac{\mathrm{d}\,\mathrm{dev}\left(\bm\sigma\right)}{\mathrm{d}\bm\epsilon}\,+ \,\mathrm{dev}\bm\sigma\otimes\frac{\mathrm{d}\Delta\lambda}{\mathrm{d}\bm\epsilon}

Notice that this expression depends on the plastic multiplier. This suggests that for the current material model, there is little benefit in choosing the Inelastic residual strain over the General stress-strain relation, because both approaches require a full stress update algorithm to compute the plastic multiplier. For other material models, such as creep models, the benefit would be greater. Using the governing scalar equation and the flow rule, we can compute the last derivative in the expression above.

\frac{\mathrm{d}\Delta\lambda}{\mathrm{d}\bm\epsilon} = \left(\frac{3G}{\sigma_{ys}}\right)^2 \, \frac{1}{1+2G\Delta\lambda} \, \frac{1}{3G+\frac{\mathrm{d}\sigma_h}{\mathrm{d}\Delta\epsilon_{pe}}} \, \left(1-\frac{2}{3}\Delta\lambda\frac{\mathrm{d}\sigma_h}{\mathrm{d}\Delta\epsilon_{pe}}\right) \, \mathrm{dev}\left(\bm\sigma\right)

In order to ensure rapid convergence of the global equation solver and ultimately reduce the simulation time, the computed Jacobian should be accurate. Well, what does accurate mean? Simply put, it means that the computed derivative must be consistent with the stress update algorithm that was used to compute stresses. That is, any assumptions or simplifications used in the stress update algorithm should be reflected in the computation of the Jacobian. A derivative based on the stress update algorithm is often called algorithmic or consistent.

In some situations, the Jacobian computation can be cumbersome. It is often possible in these situations to use an approximate Jacobian. Keep in mind that the accuracy of the solution is determined by the stress update algorithm. As long as the Jacobian is not too far off, the global equation solver will still converge to the correct solution, albeit at a lower rate of convergence.

Further Comments on Jacobian Computation

  • The external material uses a reduced vector and matrix form for stresses, strains, and Jacobians. We need to phrase our stress update algorithm and Jacobian computation accordingly.
  • The strain vector that enters into the external material uses tensor components for the shears. The Jacobian is defined assuming tensor components of the shears; therefore, it is generally nonsymmetric.
  • The case of plane stress is handled outside the external material feature, so we only need to make one implementation of our model. In other words, the model we implement can be directly used in 3D, 2D axial symmetry, and 2D plane stress and plane strain analyses.

Application Example: Specialization of the Material Model

In the sections above, we developed a stress update algorithm and outlined how to compute a Jacobian. Now, we will consider a special case for the hardening curve. We assume that the yield stress is a linear function of effective plastic strain. This is usually called linear hardening, and it is defined by a constant “plastic modulus” E_{iso}, which is the constant slope of the hardening curve. As it turns out, linear hardening means that the plastic strain increment can be solved on closed form:

\Delta\epsilon_{pe} = \frac{\sigma_{trial,e}-^{old}\sigma_{ys}}{3G+E_{iso}}

In the example’s source code file, we have made use of this specialization into linear hardening.

Let’s consider an example problem of pulling a plate with a hole.

A schematic of a plate with a hole.
Dimensions, boundary conditions, and loads for a plate with a hole.

The problem has two symmetry planes and only one-quarter of the plate is modeled. We use the following material parameters:

  • E = 70 GPa
  • ν = 0.2
  • \sigma_{ys0} = 243 MPa
  • E_{Tiso} = 2.17 GPa

The problem assumes plane stress. We can compare predictions of the implementation of our material model with the built-in counterpart. We expect differences only within the order of numerical round-off as long as the tolerance when solving the nonlinear equations is tight enough.

A COMSOL Multiphysics model of the effective von Mises stress for the external material model.
A COMSOL Multiphysics model of the effective von Mises stress for the built-in material model.

Computations of the effective von Mises stress in MPa of the external material implementation (left) and built-in material model (right).

A COMSOL Multiphysics model of the effective plastic strain for the external material model.
A COMSOL Multiphysics model of the effective plastic strain for the built-in material model.

Computations of the effective plastic strain for the external material implementation (left) and built-in material model (right).

Some Concluding Remarks

There are actually more scenarios where you can employ the possibility to add external materials.

Consider a situation where you have a source code for a material model that has been verified in another context. You may have created it yourself or found it in a textbook or journal paper. In this case, it may be more efficient to use the external material functionality than casting it into a new form and enter it as set of ODEs. Even when the code is written in, for example, Fortran or C++, it is usually rather straightforward to wrap it into the C interface used by the external material.

A coded implementation may be computationally more efficient than using the User Defined or extra PDE options. The reason is that the detailed knowledge about the material law makes it possible to devise efficient stress updates, using, for example, local time stepping.

You may want to distribute your material model in a compiled form so that the end user cannot access the source code. As a matter of fact, the third-party product PolyUMod is implemented this way.

Next Steps

  • Download the example source code and model files featured in this blog post (linear_hardening.mph, linear_hardening.c, linear_hardening.dll) and get more details on how to use external materials:
  • Read a previous blog post about accessing external material models
  • Learn more about computational plasticity with these highly informative books:
    • J.C. Simo and T.J.R. Hughes, Computational Inelasticity, Interdisciplinary Applied Mathematics, Vol. 7, Springer-Verlag, 1998.
    • M. Kojic and K.J. Bathe, Inelastic Analysis of Solids and Structures, Computational Fluid and Solid Mechanics, Springer-Verlag, 2005.
    • T. Belytschko, W.K. Liu, and B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, 2000.

PolyUMod software is developed by Veryst Engineering LLC. COMSOL AB and its subsidiaries and products are not affiliated with, endorsed by, sponsored by, or supported by Veryst Engineering LLC.


  1. Ijaz Shahid November 11, 2017   8:13 am

    Respected Professor
    I have a question about the elastic-plastic material modeling in COMSOL Multiphysics. Is there is any modeling in COMSOL where i can design my model in elastic plastic region.

  2. Mats Danielsson November 13, 2017   2:44 am


    As mentioned in the blog post, there are numerous material models for elasto-plasticity already in the program. If you are interested in modeling the elasto-plastic behavior of a material not already provided, you can use one of the External Material approaches that are discussed in the blog post. To get started, you can download example files using the link under “Next Steps”.


  3. KAI REN November 15, 2017   10:35 pm

    Dear Dr. Danielsson,

    I have a question about transferring strain fields between two different mph files. Is there a method to transfer strain field in different mph files (when thermal strain, plastic strain are involved), so the later mph file can utilize the result of the former mph file as initial value?
    As there exists “initial strain and stress” tag, I wonder how to add initial plastic strain and thermal strain to let comsol distinguish constitution of different strains.

    Many thanks,

  4. Mats Danielsson November 16, 2017   1:07 am

    Dear Kai,

    This ia a bit off-topic with regard to the External Material feature, so I suggest that you contact our support.


  5. hicham CHAOUKI August 23, 2018   5:35 pm

    Hi Mats,
    Thank you for the detailed informations. It is possible to implement the same constitutive law according to the second approach mentioned above (using the PDE tools with weak or strong formulation) ?
    Best regards,
    Laval University

  6. Mats Danielsson August 24, 2018   2:24 am


    Material models like the one in this blog post use a yield condition. That is why the stress update has to be expressed in terms of an algorithm, where a certain “branch” is traversed depending on the state of the yield condition (elastic vs. elastoplastic). In other words, the elastoplastic behavior cannot be expressed purely in terms of, say, an ODE.

    That said, many other material models, such as creep or viscoelastic models, do not feature a yield condition. In these models, inelastic strains are said to develop at every non-zero level of stress, albeit very slowly at low levels. These types of models can be implemented the way you propose.


  7. Narinder Singh August 24, 2018   2:42 pm

    Hi Mats,

    Thanks for this useful blog. I have a question regarding the total strain increment in the accompanying files.

    For calculating trial stress (or strain), shouldn’t the strain increment be equal to e(n+1)-ep(n) ?

    Your code writes it as eTrial[i] = e[i] – ep[i], where both e and ep correspond to the current values, making it effectively e(n)-ep(n). Where does the ‘increment’ part come in?

    Or if I ask in general, how does one get to access strain increments (for each integration point of each element)? As far I believe, the equation view only shows the current state of quantities.

    I am sure I have missed something somewhere. Please clarify.


  8. Mats Danielsson August 24, 2018   3:12 pm


    The trial stress is the stress that is computed assuming that the increment is purely elastic. Therefore it defines a trial strain as the total strain of the current (global) iteration less the total plastic strain of the previous increment. If you look closer at the code, the plastic strains are using values from the states1 array. This states1 array has not yet been updated when the plastic strains are defined, so the values do indeed represent the converged values of the previous increment. If the trial stress violates the yield condition, the code branches off into the elastoplastic part where it computes an increment in plastic strains, and adds it to states1. And so on.

    It’s not a strain increment that is used, though. This formulation is total, in the sense that it does NOT use strains from the previous increment to compute an increment in stress. It only uses total strains from the current (global) iteration and the plastic state from the previously converged global increment.

    I hope that helps!


  9. Narinder Singh August 24, 2018   4:18 pm

    Thanks for the explanation Mats.

    I agree with you in principle. Does it mean then that the ‘input’ strains in the code, ‘e’ are the current total strains (effectively, e(n+1)) and not e(n) as I was thinking them to be. In that case, it all makes sense.

  10. Mats Danielsson August 24, 2018   4:52 pm

    Hello again.

    Yes, the strains that are fed into the code correspond to the current estimate of the displacement field.


  11. hicham CHAOUKI September 12, 2018   5:42 pm

    Thanks Mats for all informations. I tried the dll file that I found on website and it works well. However, wen I write my own code and I built up the dll file, Comsol does not work. Could you please give us some informations about the compiler that we need and the methodology to follow.
    Best regards,
    Laval University

  12. Mats Danielsson September 13, 2018   2:01 am


    With regard to compilation, I suggest you read the section called “How to Compile and Link an External Material Model” in the documentation.

    Without more detailed information, it is difficult to know what is going wrong in your particular case. I suggest you contact our support.


  13. hicham CHAOUKI October 3, 2018   3:30 pm

    Thank you Mats. Everything is ok with the compilation.
    Best Regards,

  14. hicham CHAOUKI November 13, 2018   2:41 pm

    Hi Mats,
    I have a question about external material routine for mechanical analysis. I’m working on the developpment of a constitutive law where the mechanical behaviour is coupled with thermal and porous behaviours. The constitutive law depends, among others, on the thermal strain and the pore pressure. I’m wondering if it is possible to pass to the external routine values of the pore pressure and the temperature field and then solve the problem by considering a weak coupling between thermal, mehcanical and porous problems.
    Bests regards,

  15. Mats Danielsson November 19, 2018   6:02 am


    Of course you can pass anything to the external material in the parameter list. With regard to the coupled problem, I advise you to contact COMSOL support as it becomes too complex to address in the blog.


  16. hicham CHAOUKI November 19, 2018   4:46 pm

    Hello Mats,
    Thank you for the information.
    Best regards,

Loading Comments...