How to Implement Elastoplasticity in a Model Using External Materials

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}

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

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

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:

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.

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:

or

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:

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:

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.

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:

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:

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:

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:

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.

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:

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.

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.

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

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:
  • External Material Examples, Structural Mechanics
• 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.