User materials

Other material laws can be defined by the user by means of the *USER MATERIAL keyword card. More information and examples can be found in section .

A user material card requires the coding of a user material subroutine, either of the CalculiX type or the Abaqus type.

For a user material routine of the CalculiX type the routine umat_user.f can be taken as example. The input primarily exists of the mechanical Lagrange strain tensor at the beginning and end of the increment, requested output is the Piola-Kirchhoff stress tensor of the second kind and the derivative of the Piola-Kirchhoff stress of the second kind with respect to the mechanical Lagrange strain tensor, both at the end of the increment. Details of other input fields is given in the section referenced above.

For a user material routine of the Abaqus type there are two possibilities: either the user wants to apply the routine to linear geometric calculations only. Then, the kind of strain and stress going in or out of the routine is not important and the previous paragraph applies. However, if the user would like to apply the routine to geometrically nonlinear calculations, the strain entering the routine is the corotational mechanical logarithmic strain and the required stress is the corotational Cauchy stress. So CalculiX has to perform some conversions before and after calling the Abaqus user material routine. For a prototype example of a Abaqus user material routine the user is referred to umat.f, for limitations on the use of the Abaqus interface fields section should be consulted.

Here, some information is given on how the fields used in CalculiX (mechanical Lagrange strain and Piola-Kirchhoff stress of the second kind) are converted into the fields required by Abaqus (corotational mechanical logarithmic strain and corotational Cauchy stress) and vice versa. The conversions is coded in umat_abaqusnl.f.

The mechanical logarithmic strain satisfies

$\displaystyle e_{\ln}^{\mathcal{M}} = \sum_{i} \lambda_i^{\mathcal{M}} \boldsymbol{n}_i \otimes \boldsymbol{n}_i,$ (486)

where $ \lambda_i^{\mathcal{M}}$ are the eigenvalues of the mechanical deformation gradient $ \boldsymbol{F} ^{\mathcal{M}}$. In calculiX, the mechanical deformation gradient is obtained by subtracting the thermal expansion from the total deformation gradient, i.e.

$\displaystyle \boldsymbol{F} ^{\mathcal{M}}= \boldsymbol{F} - \alpha \Delta T \boldsymbol{I}$ (487)

for isotropic expansion. The rotation tensor $ \boldsymbol{R}= \boldsymbol{n}_{\underline{i}}
\otimes \boldsymbol{N}_{\underline{i}}$ rotates the principal vectors $ \boldsymbol{N}_i$ of the motion in material coordinates in spatial coordinates $ \boldsymbol{n}_i$:

$\displaystyle \boldsymbol{n}_i = \boldsymbol{R} \cdot \boldsymbol{N}_i.$ (488)

It is a two-point tensor with one leg in the material frame of reference and one leg in the spatial frame of reference. The corotational mechanical logarithmic strain then amounts to (recall that $ \boldsymbol{R}$ is orthogonal, i.e. $ \boldsymbol{R}^{-1} = \boldsymbol{R}^T$) :

$\displaystyle \hat{\boldsymbol{e}}_{\ln}^ {\mathcal{M}} = \boldsymbol{R}^T \cdo...
...} = \sum_{i} \lambda_i^{\mathcal{M}} \boldsymbol{N}_i \otimes \boldsymbol{N}_i.$ (489)

$ \lambda_i^{\mathcal{M}}$ can be obtained from the eigenvalues of the mechanical Lagrange tensor $ \boldsymbol{E}^ {\mathcal{M}}= (\boldsymbol{F}^ {\mathcal{M},T} \cdot
\boldsymbol{F}^ {\mathcal{M}} - \boldsymbol{I})/2 $ and its eigenvectors are $ \boldsymbol{N}_i$.

The Cauchy tensor satisfies

$\displaystyle \boldsymbol{\sigma } = J^{-1} \boldsymbol{R} \cdot \boldsymbol{U} \cdot \boldsymbol{S} \cdot \boldsymbol{U}^T \cdot \boldsymbol{R}^T,$ (490)

and consequently the corotational Cauchy tensor amounts to

$\displaystyle \hat{\boldsymbol{\sigma }} = \boldsymbol{R}^T \cdot \boldsymbol{\...
...dsymbol{R} = J^{-1} \boldsymbol{U} \cdot \boldsymbol{S} \cdot \boldsymbol{U}^T.$ (491)

Here, $ J$ is the Jacobian determinant, $ \boldsymbol{S}$ is the Piola-Kirchhoff stress of the second kind and $ \boldsymbol{U}$ is the right-stretch tensor. The latter is the square root of the Cauchy-Green tensor $ \boldsymbol{C}$ and therefore it can also be derived from the Lagrange tensor. However, in routine umat_abaqusnl.f the mechanical Lagrange tensor is available and not the total one. So a relationships is needed between the total right-stretch tensor and the mechanical right-stretch tensor. Since $ \boldsymbol{U}$ has the same eigenvalues as $ \boldsymbol{F}$ and the multiplicative decomposition of the deformation gradient in a mechanical and thermal deformation gradient leads to $ \lambda = \lambda ^ {\mathcal{M}} (1 +
\alpha \Delta T)$ one obtains:

$\displaystyle \hat{\boldsymbol{\sigma }} = J^{-1} \boldsymbol{U}^ {\mathcal{M}}...
...\boldsymbol{S} \cdot (\boldsymbol{U}^ {\mathcal{M}})^T (1 + \alpha \Delta T)^2.$ (492)

The latter equation assumes that the thermal expansion is isotropic. In routine umat_abaqusnl.f $ \boldsymbol{S}$ is known, $ \boldsymbol{U}^
{\mathcal{M}}$ can be obtained from $ \boldsymbol{E}^ {\mathcal{M}}$ and $ (1 +
\alpha \Delta T)$ results from

$\displaystyle (1 +\alpha \Delta T) = \sqrt{\frac{(\boldsymbol{F}^T \cdot \boldsymbol{F})_{11} }{2 \boldsymbol{E}^ {\mathcal{M}}_{11}+1 }},$ (493)

where the index $ _{11}$ refers to element (1,1) of the corresponding matrix. In this way the corotational Cauchy stress is obtained from the Piola-Kirchhoff stress of the second kind. The inverse relationship is used after returning from the Abaqus material user subroutine:

$\displaystyle \boldsymbol{S} = J \boldsymbol{U}^{-1} \cdot \hat{\boldsymbol{\sigma }} \cdot \boldsymbol{U }^{-T}.$ (494)

Now, $ \partial \boldsymbol{S}/\partial \boldsymbol{E}^{\mathcal{M}}$ is needed in CalculiX (i.e. total Lagrangian approach), however, the Abaqus routine returns $ \partial \hat{\boldsymbol{\sigma}} /\partial \hat{\boldsymbol{e}}_{\ln}^{\mathcal{M}}$. Performing an exact transformation is very tedious and nearly impossible. It would require an enormous computational effort and frequently results in an asymmetric $ \partial \boldsymbol{S}/\partial \boldsymbol{E}^{\mathcal{M}}$ - matrix even if $ \partial \hat{\boldsymbol{\sigma}} /\partial \hat{\boldsymbol{e}}_{\ln}^{\mathcal{M}}$ is symmetric. To avoid longer computational times due to the use of an asymmetric linear equation solver a symmetrization is usually performed leading to an approximation anyway. Therefore, a simplified approach is taken which reduced the computational effort to a minimum while leading to acceptable convergence rates despite the approximations involved.

From the relationship between $ \boldsymbol{S}$ and $ \hat{\boldsymbol{\sigma
}}$ one obtains (in Carthesian coordinates)

$\displaystyle \left ( \frac{\partial \boldsymbol{S}_{IJ}}{\partial \boldsymbol{...
... {\partial \boldsymbol{E}_{KL} ^{\mathcal{M}}} \right) \boldsymbol{U}^{-1}_{NJ}$ (495)

(recall that $ \boldsymbol{U}$ is symmetric). In the above equation the dependence of $ \boldsymbol{U}^{-1}$ on $ \boldsymbol{E}^ {\mathcal{M}}$ has been neglected (first approximation). One can further write:

$\displaystyle \left ( \frac{\partial \boldsymbol{S}_{IJ}}{\partial \boldsymbol{...
... {\partial \boldsymbol{E}_{KL} ^{\mathcal{M}}} \right) \boldsymbol{U}^{-1}_{NJ}$ (496)

or, approximating the corotational mechanical logarithmic strain $ \hat{\boldsymbol{e}}_{\ln} ^{\mathcal{M}}$ by the corotational mechanical Eulerian strain $ \hat{\boldsymbol{e}} ^{\mathcal{M}}$ (second approximation):

$\displaystyle \left(\frac{\partial \boldsymbol{S}_{IJ}}{\partial \boldsymbol{E}...
...{\partial \boldsymbol{E}_{KL} ^{\mathcal{M}}} \right) \boldsymbol{U}^{-1}_{NJ}.$ (497)

The Eulerian strain satisfies:

$\displaystyle 2\boldsymbol{e}= \boldsymbol{F}^{-T} \cdot 2\boldsymbol{E} \cdot \boldsymbol{F}^{-1}.$ (498)

which amounts to:

$\displaystyle \left( \boldsymbol{I} - \frac{{\boldsymbol{F}^{\mathcal{M}}}^{-T}...
...al{M}} (1+\alpha \Delta T)^2 -\boldsymbol{I} \right) \cdot \boldsymbol{F}^{-1}.$ (499)

This can be written as:


$\displaystyle \boldsymbol{I}$ $\displaystyle -$ $\displaystyle {\boldsymbol{F}^{\mathcal{M}}}^{-T} \cdot
{\boldsymbol{F}^{\mathcal{M}}}^{-1} + 2\xi \boldsymbol{I}$  
  $\displaystyle =$ $\displaystyle \boldsymbol{F}^{-T}\cdot \left({\boldsymbol{F}^{\mathcal{M}}}^{T}...
...ymbol{F}^{\mathcal{M}} \right) \cdot \boldsymbol{F}^{-1} (1+\alpha \Delta T)^2,$ (500)

where

$\displaystyle 2\xi :=2 \alpha \Delta T + \alpha^2 (\Delta T)^2.$ (501)

This leads to:

$\displaystyle \boldsymbol{e} ^{\mathcal{M}} + \xi \boldsymbol{I} = {\boldsymbol...
... \xi \boldsymbol{C} ^{\mathcal{M}} ) \cdot {\boldsymbol{F}^{\mathcal{M}}}^{-1},$ (502)

or

$\displaystyle \boldsymbol{e} ^{\mathcal{M}} = {\boldsymbol{F}^{\mathcal{M}}}^{-T} \cdot \boldsymbol{E} ^{\mathcal{M}} \cdot{\boldsymbol{F}^{\mathcal{M}}}^{-1}.$ (503)

Consequently, the corotational mechanical Eulerian strain amounts to:

$\displaystyle \hat{\boldsymbol{e}}^{\mathcal{M}}_{PQ} = \boldsymbol{R}^T_{PR} \...
...PR} \boldsymbol{E}^{\mathcal{M}}_{RS} {\boldsymbol{U}^{\mathcal{M}}}^{-1}_{SQ},$ (504)

leading to (again using the first approximation for the inverse mechanical right-stretch tensor)

$\displaystyle \left( \frac{\partial \hat{\boldsymbol{e}}^{\mathcal{M}}_{PQ} }{\...
...{M}}}^{-1}_{PR} (\mathbb{I}_I)_{RSKL} {\boldsymbol{U}^{\mathcal{M}}}^{-1}_{SQ},$ (505)

where

$\displaystyle (\mathbb{I}_I)_{RSKL} := (\delta_{RK} \delta_{SL}+\delta_{RL}\delta_{SK})/2.$ (506)

Due to the symmetry of the corotational mechanical Eulerian strain tensor and the symmetry of the inverse right-stretch tensor, one can write after substitution of Equation (506) into Equation (505):

$\displaystyle \left( \frac{\partial \hat{\sigma}_{MN}} {\partial \hat{\boldsymb...
...oldsymbol{U}^{\mathcal{M}}}^{-1}_{PK} {\boldsymbol{U}^{\mathcal{M}}}^{-1}_{QL}.$ (507)

One finally obtains:

$\displaystyle \left( \frac{\partial \boldsymbol{S}_{IJ}}{\partial \boldsymbol{E...
...oldsymbol{U}^{\mathcal{M}}}^{-1}_{KP} {\boldsymbol{U}^{\mathcal{M}}}^{-1}_{LQ}.$ (508)

This is the expresson which is used for the material tangent. The approximations do not lead to a false result, they just slow down the convergence. Since the wish to avoid asymmetric matrices precludes quadratic convergence anyway the above approximations seem tolerable. This has been confirmed in practical numerical simulations of e.g. the deformation plasticity model and the Johnson-Cook model in implicit dynamic calculations (in explicit dynamic calculations the material tangent is not needed).