search

jfriedlein / usrmat_LS-Dyna_Fortran

Basics to implement user-defined materials (usrmat, umat) in LS-Dyna with Fortran
53 stars 28 forks source link
documentation fortran90 ls-dyna material-model ttb tutorial umat utan

readme

usrmat_LS-Dyna_Fortran

Basics to implement user-defined materials (usrmat, umat, utan) in LS-Dyna with Fortran

What is this all about?

LS-Dyna offers the interfaces and solvers to, among many other things, simulate mechanical systems and the related material behaviour. To reproduce realistic material responses, we need to utilise adequate material models. In case, standard available models cannot generate valid results, new user-defined material models can be implemented. The latter typically requires a stress-strain routine (umat, computes the stress for a given strain) and for implicit analyses also the stress-strain tangent (utan, how the strain changes the stress). This guide introduces the basics to implement user-defined material models in LS-Dyna using the standard Fortran interface.

Software requirements and suggestions

To implement and apply UMAT-routines we recommend the following software and tools:

@todo Add more infos on MPI (currently it only works for R>=10.2) Install MPI via "msmpisdk.msi" https://www.researchgate.net/post/How_to_solve_LINK_fatal_error_LNK1181_cannot_open_input_file_msmpilib_during_Abaqus_and_fortran_linking, so that the desired files 'msmpi.lib' and 'msmpifec.lib' are created. Then link them inside the makefile, for instance as MPICH_LIB = /libpath:"C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x64" msmpi.lib msmpifec.lib and MPICH_INC_PATH = "C:\Program Files (x86)\Microsoft SDKs\MPI\Include".

For the coding:

@todo Add an open source option (e.g. gfortran?, Notepad (Syntax: "Fortran with fixed format")), check under Linux

Test setup and compiler

To check whether everything works as desired, we simply compile the original LS-Dyna object files.

  1. Unzip the object version (*_lib.zip) into your working directory {working directory}.

  2. Start the Intel Fortran Compiler, e.g. via start menu -> Intel Parallel Studio XE 2017 -> Compiler 17.0 Update 6 for Intel 64 Visual Studio 2017 environment. A command window should open up.

  3. Change the directory of the command window to your working directory.

  4. Run the command nmake (Starts the nmake.exe) to compile the Fortran files and create the lsdyna.exe

  5. Wait a bit until all files are compiled. When you run nmake for the first time, it usually takes a bit longer, because every file is compiled. Later on, when you implement material models, you typically only change two files, hence nmake only recompiles those two, which is much faster.

  6. After the compilation finished successfully, you should have a file lsdyna.exe in your working directory. This standalone executable contains LS-Dyna together with your material models.

@todo Test possible issues when passing the "standalone" exe around.

  1. Start an LS-Dyna simulation, for instance from LS-Run, where you choose the created lsdyna.exe as the executable (e.g. instead of [...]/ls-dyna_smp_d_R11.1_winx64_ifort160.exe you choose {working directory}/lsdyna.exe).
  2. The simulation should run just the same as using the LSTC version 'ls-dyna_smp_d_R11.1_winx64_ifort160.exe`.

@todo Check whether there is indeed a noticable performance loss by using the self-compiled exe.

  1. Now we can start implementing our own material models into the files of the working directory and compile and run it using the above steps.
  2. For some reason, when calling nmake, even when you changed the Fortran files, the message 'lsdyna.exe' is up-to-date is shown and the files aren't compiled. As a simple workaround, just delete the lsdyna.exe everytime before you call nmake. This can be accomplished using the start_nmake.bat available in this repository. First, you should enter the paths to your working directory inside the batch file. Then just execute the .bat file and enter "del lsdyna.exe & nmake", which deletes the possibly already existing lsdyna.exe and execute nmake to create a new lsdyna.exe.

Implementation

Before we take a closer look at the files that need to be extended by our material models and how we do this. A few notes on the used programming language FORTRAN.

A few notes on FORTRAN

The programming language Fortran and especially the older version used in LS-Dyna (with the above setup it's FORTRAN90) has some "features" that might (most certainly) be unknown or unexpected to programmers used to "more modern" languages, such as C++, Matlab, Python, ... So to start with you best consider a few Fortran tutorials or books to get to know the specifics of the language. If you are still reading and skipped the last recommendation, you might be as naive as me. So for the impatient, but experienced C++/Matlab/Python programmers a few remarks to ease the start. Some "features" originate from the time of punchcards (Lochkarten), so maybe it helps to keep this motivation in mind. These notes won't teach you how to code Fortran, but only try to give you a head start. The specific syntax still needs to be figured out by you.

c In a single line:
      i = 1 + 2 + 3 + 4
c Or equivalently split up
      i = 1 + 2
     &      + 3
     &  + 4  ! The horizontal position of the continued line does not matter.

@todo Check use of implicit none and required additions to existing ls-dyna code.

@todo Be careful with "Fortran features", execution order, function/subroutine, use of intent(in), ..., also test in parallel; http://www.cs.rpi.edu/~szymansk/OOF90/bugs.html#2

a = 1./3. or a = 1.0/3.0

Here, the dots that follow the integer, mark them as floating point numbers and hence the variable a is '0.333...'. As a recommendation, always add the dot after coded numbers e.g. also in

a = 2. * 5. * 0.5

because you never know when you might change your mind and rewrite parts of equations using a divison.

Some more basics on Fortran (and Abaqus user interfaces) can be found in this PDF by the Brown University.

Our first user material

  1. Open your working directory (the folder with the unpacked object version, e.g. ls-dyna_smp_d_R11_1_0_139588_winx64_ifort2017vs2017_lib) in Visual Studio.
  2. Implement your material model code (computation of stress, history variables ...) in the file dyn21umats.F (here 'umats' stands for umat-scalar, change the file extension from '.f' to '.F' in case the "#define" and "#include" cause errors), for instance linear elasticity. We code our model in the first unused umat, here it's umat43. Note that we right away start with tensor based models, as noted below, which require the ttb files under the path ttb/... in the working directory.

[dyn21umats.F]:

c Compiler flags and include for the tensor toolbox
#include 'ttb/ttb_library.F'
[...]
      subroutine umat43 (cm,eps,sig,epsp,hsv,dt1,capa,etype,tt,
     1 temper,failel,crv,nnpcrv,cma,qmat,elsiz,idele,reject)
c
c Use the tensor toolbox
c The position of this 'use' is crucial (first entry)
      use Tensor
c Standard LS-Dyna includes and "declarations"
      include 'nlqparm'
      include 'bk06.inc'
      include 'iounits.inc'
      dimension cm(*),eps(*),sig(*),hsv(*),crv(lq1,2,*),cma(*),qmat(3,3)
c Input and output arguments (see the beginning of the 'dyn21umats.F' file)
c cm: contains the material parameters, set in the material card (P1, ...)
c eps: increments (!) in the strain (difference between current strain and last converged load step),
c      stored as eps_11, eps_22, eps_33, eps_12, eps_23, eps_31 (coordinates: 1=x, 2=y, 3=z)
c sig: stresses from the last converged load step, stored in the same order as the strain 'eps'
c hsv: list of history variables (array length: NHV, set in the material card)
      integer nnpcrv(*)
      logical failel,reject
      character*5 etype
      INTEGER8 idele
c Declarations
      real lame_lambda, shearMod_mu, bulkMod_kappa
      type(Tensor2) :: Eye, d_eps, stress, stress_n
c Material parameters, extracted from the material parameter array 'cm'
      lame_lambda = cm(1)
      shearMod_mu = cm(2)
      bulkMod_kappa = lame_lambda + 2./3. * shearMod_mu
c Second order identity tensor
      Eye = identity2(Eye)
c
c The ttb-function 'strain' transforms the incremental strain 'eps' from the vector
c notation to the ttb tensor data type (correct assignment of vector to
c tensor index AND Voigt-factor 0.5).
      d_eps = strain(eps,3,3,6)
c
c The function 'symstore_2sa' stores the stress array 'sig' as the tensor 'stress_n',
      stress_n = symstore_2sa(sig)
c
c Now follows our usual stress equation in tensor notation instead of Voigt, finally
c ... happily, thanks to Andreas Dutzler [https://github.com/adtzlr/ttb]!
c Be aware, that our inputs are the strain increments, so we have to modify our
c usual equations a bit by computing the new stress from the old stress 'stress_n'
c plus a stress update. (Gets irrelevant when you use finite strains and
c the deformation gradient.)
      stress = stress_n + lame_lambda * tr(d_eps) * Eye
     &                  + 2. * shearMod_mu * d_eps
c Transform the stress tensor back into a vector:
      sig(1:6) = asarray(voigt(stress),6)
c
c Everything is done with just a few lines of code ... perfect
c
      return
      end

(Note how compact, slim, insensitive to errors and beautiful tensors can be.)

  1. Following the above test setup: Delete the possibly existing lsdyna.exe in the working directory. Start the Fortran compiler as described in the above steps 2 and 3.
  2. Run nmake.
  3. The code should successfully be compiled and no error message ought to be shown when the lsdyna.exe was created.
  4. Run an LS-Dyna simulation using the freshly compiled lsdyna.exe and select your material model in the material cards (outlined below in the section material card). Also consider the section on the solver settings below.
  5. Often you can test the material model even without having to implement the tangent (note: We haven't implemented the tangent yet, so step 6 will fail for an implicit analysis) by running your simulation as explicit, which does not at all require a tangent (motto of explicit computation: 'Don't look, just go.'). If the material response is as expected, you can again switch to implicit.
  6. Next we implement our tangent. This must be done in the file dyn21utan.F (user tangent) and in the subroutine with the same ID (here 43, so utan43). For linear elasticity the tangent is constant and independent of the input arguments. For more complicated models the umat-utan function split can cause some problems that are discussed later on. Here we just implement the constant tangent, which in Voigt notation fills a 6x6 matrix called es:

[dyn21utan.F]:

[...]
      subroutine utan43(cm,eps,sig,epsp,hsv,dt1,unsym,capa,etype,tt,
     1 temper,es,crv,nnpcrv,failel,cma,qmat)
c
c******************************************************************
c The computations of the tangent in Tensor notation also requires the split of
c the code into three parts.
c 1. Transform the arrays to tensors.
c 2. Compute the quantities with tensors.
c 3. Transform the tensorial results back to arrays.
      use Tensor
      include 'nlqparm'
      dimension cm(*),eps(*),sig(*),hsv(*),crv(lq1,2,*),cma(*)
c cm, eps: as  above
c sig: output stress from the last call to the material model umat (not from the last converged step)
c hsv: output history from the last umat (not from the last converged step)
c For details: see section below
      integer nnpcrv(*)
      dimension es(6,*),qmat(3,3)
      logical failel,unsym
      character*5 etype
c Declarations
      type(Tensor2) :: Eye
      type(Tensor4)  :: tangent_C, IxI, I_dev
c
      real lame_lambda, shearMod_mu, bulkMod_kappa
c Material parameters
      lame_lambda = cm(1)
      shearMod_mu = cm(2)
      bulkMod_kappa = cm(1) + 2./3. * shearMod_mu
c Second order identity tensor
      Eye = identity2(Eye)
c Fourth order tensors
      IxI = Eye.dya.Eye
      I_dev = (Eye.cdya.Eye) - 1./3.*(Eye.dya.Eye)
c Compute the tangent modulus as a fourth order tensor
      tangent_C = bulkMod_kappa * IxI
     &            + 2. * shearMod_mu * I_dev
c Transform tensor 'tangent_C' into the matrix 'es'
      es(1:6,1:6) = asarray(voigt(tangent_C),6,6)
c
      return
      end

Some notes on proper simulation and solver settings for testing umats

The setup of an LS-Dyna keyword file for the comprehensive testing of umats is still "work in progress".

To check the correctness of the umat (first ignoring the utan), you can use an explicit simulation as stated above. To check your implemented tangent, you must use an implicit computation and best get to the bottom of your resulting convergence rate (requires various solver settings, a snippet is given in the following figure). Nader Abedrabbo outlines some very good aspects in his UMAT workshop in the section on "UMAT verification" (unfortunately only for umat, not utan)

When you check the resulting convergence rate (to verify your tangent), you might be surprised that you hardly see quadratic convergence, especially for standard tests with distortion or shear and for every practical example. The culprit here seems to be the LS-Dyna element formulation. From a first glance, you might think that ELFORM=2 is similar to a normal linear Q1 element. It is NOT! LS-Dyna often uses (some kind of) selective reduced integration (SRI, S/R), however in such a way that I cannot understand how the formulation shall be consistent, because they seem to use some averaging. Hence, the convergence rate only approaches 2 for simple undistorted element tests under pure tension (so no shear). Without going into more details here, you might see the approaching dilemma of no decent convergence. work in progress

More details on the setup of these simulations will be given here in the future (small appetiser: "Numerical examples for LS-Dyna").

@todo maybe add one-element test, explain the settings and finish this section

Material models using tensors

In case you like the above equations in Tensor notation and you are not familiar with the in LS-Dyna usual Voigt notation. There is a superb Tensor Toolbox for Fortran (ttb, https://github.com/adtzlr/ttb) with a comprehensive documentation (https://adtzlr.github.io/ttb/) that enables you to use tensors and tensor operations. Regarding the setup of the toolbox, I hand you over to the capable hands of Andreas outlining the steps here https://adtzlr.github.io/ttb/installation.html. You can find a more advanced example in the ttb documentation specific for LS-Dyna (currently e.g. at: "LS-Dyna Tensor Neo-Hooke").

@todo Add the link when the example is added. Refer to some more example files (elasto-plasticity, ...)

Outline of the interface for umat and utan

The following figure shows the typically relevant input/output arguments.

@todo Check if the tangent 'es' can also be used as an input in utan (e.g. computing an incremental tangent, sounds like bs, right?)

The UMAT-subroutine gets the incremental strain eps in Voigt (!) notation. Hence, the shear components contain twice the xy, yz and zx components (all automatically handled by the ttb tensor toolbox).

However, because the input strain is an incremental strain, the typical tensor based models often need to be adapted, cause they are typically based on the total strain tensor. If you want to avoid this, you can set "IHYPER=1" (see below on material cards) and use the deformation gradient to compute the total strain as

@todo Here we need to add some more information on hypo and hyper formulations and Updated Lagrangian, geometric nonlinearity, etc.

The "symmetrize" function is available in the tensor toolbox extension package ttbXkinematics (currently not available online)

@todo Upload the ttbXkinematics package

Here, the index l denotes the values from the current iteration and n the values from the last converged load step. Some background information and context regarding the herein used notation can be found in the following scheme, before we move on.

Secondly, we receive the stresses sig that contain the Cauchy stress from the last converged load step n. The history variables, such as the plastic strain or the hardening for plasticity, are summarised in the list hsv. Material parameters, like the Young's modulus or Poisson's ratio, set in the material card, are stored in the list cm. For the above example the first and second Lame parameters are stored in P1 and P2, respectively. To avoid mixing up the material parameters, you can use the "material parameter manager". Lastly, we can also find the deformation gradient F in the history array after setting the option IHYPER in the material card (see section 'material card').

Now the material model must compute the new Cauchy stress (index tmp above) and update the history variables.

In the UTAN-routine (tangent) we receive the temporary Cauchy stress and history from UMAT. With this we have to compute the tangent. In the world of tensors the latter needs to be the fourth order Eulerian tangent moduli E.

Some considerations on the split of material model (umat) and tangent (utan)

Unfornuately, because LS-Dyna was born "explicitly" (but contains full implicit capabilites, see "LS-Dyna solvers", apart from emerging doubts) the material model routine umat (stress and history update) and the tangent utan (for implicit only) are located in separate subroutines and files. Hence, at first glance it is not straightforward to generate a complicated consistent tangent lacking access and knowledge of quantities computed in the umat routine. For elasto-plasticity as an example, inside the umat routine you determine whether the current step is elastic or plastic and compute the stress accordingly. Obviously, the tangent in utan needs to be set up such that it fits to the chosen elastic or plastic response. But the deciding criterion (here the yield criterion) or quantity (plastic multiplier) is not available in utan (and might be cumbersome to salvage).

So, what options do we have to work around this issue/complication? I came up with the following variants (if you found alternatives, please let me know):

@todo add an example of the "es in hsv" approch and the boilerplate code for utan

Example material card

All of the above was done without even considering LS-Dyna or its pre-/postprocessing (here done via LS-PrePost). In order to apply the material model to a simulation, you need to set up a card for the keyword file. First, create a card *MAT_USER_DEFINED_MATERIAL_MODELS to refer to your user material. The material id in the option "MT" must equal the id in umatXX and utanXX (above it was 43). The value of "NHV" sets the number of used history variables, here (for elasticity) none are used in the material model. The option "IHYPER=1" stores the deformation gradient in the history "hsv" on top of the defined "NHV" (is not necessary for the above small strain elasticity). The parameters "P1" and "P2", for instance, contain the Young's modulus in "cm(1)" and the Poisson ratio in "cm(2)" in the picture, respectively. For our above elasticity material, we would type the values of the first and second Lame parameters into P1 and P2.

@todo Also give the remaining material cards, ideally a good set for proper implicit analysis. E.g. use ILIMIT=1 for full Newton-Raphson, use IGS=1 for a usable linearisation, ...

validation

Elasto-plasticity with linear isotropic hardening can be compared to MAT24. Just be aware that "ETAN" is not the linear hardening modulus "K", but the resulting stiffness ETAN=(E*K)/(E+K) with Young's modulus "E".

2D: plane strain and axial symmetry

Looking for some plane strain examples? Here is how to use plane strain: https://ftp.lstc.com/anonymous/outgoing/support/FAQ/2d_general_condensed

But: For plane strain (considered as shell in LS-Dyna) the deformation gradient that enters the material model is in the corotational system of each element. So, it is not the standard deformation gradient. Be very careful with shell elements (stress/strain expressed in local coordinates system, ...).

@todo What happens for axisym? Transformation to "local material coordinate system" (https://ftp.lstc.com/anonymous/outgoing/support/FAQ/user_defined_materials.faq)

The material models typically utilise the fully 3D state.

Helpful information on the element formulations/integration:

https://ftp.lstc.com/anonymous/outgoing/support/FAQ/2d_general_condensed

tricks

The initialisation step can be detected by checking whether the time step size is zero (e.g. dt1<1e-10 for implicit integration). This can also be helpful when remeshing is used, because LS-Dyna solves an initialisation step at every remeshing (but not all releases use dt1 in this initialisation step, it works e.g. for R13). So, in case of some problems or features that need to be done after remeshing and the solution transfer, you can detect this initialisation step via dt1 and do something special for this step.

Very nice Fortran features: http://www.netlib.org/ccm/page/api/optional.html

Fortran coding standards: http://jules-lsm.github.io/coding_standards/

http://math.hawaii.edu/wordpress/fortran-3/

enumerators via "integer, parameter :: enum_a=1, enum_b=2"

hsv and cm manager

you can open LS-Dyna d3plot files in ParaView and plot values along a line, scale the deformation, ...

FAQ: https://ftp.lstc.com/anonymous/outgoing/support/FAQ/

use "call cstop("E r r o r T e r m i n a t i o n")" to properly end the simulation and free the licence

list of warning and error codes: http://support.moldex3d.com/r15/en/post-processing_errorandwarningmessages_feainterfaceforsolidmodel.html

model checker: https://2021.help.altair.com/2021/hwdesktop/engsol/topics/user_interface/model_checker_ls-dyna_r.htm

use caller environment to speed-up compilation and still use LS-Dyna context

Code design

note

However, this also means that changes in external files, e.g. libraries or outsourced codes, are not automatically detected and need to be compiled manually. The latter can be achieved by simply messing up the path in the include command, calling nmake, correcting the wrong path and calling nmake again. Now this include file is also recompiled.

@todo Find a better and automated way to avoid this silly approach. @todo Check the makefile modifications to compile outsourced/external fortran files

Hypo vs hyper

strain increment (rate-type) vs deformation gradient

for hypo be aware of necessary rotations for tensorial history variables like back stress tensor using the spin

@todo

Some notes on Linux

get the AVX version "grep flags /proc/cpuinfo"

redhat vs centos

eclipse photran

"0.5" vs "5.0d-1" @todo Let Andreas know that in file "libexp.f" line 55 "5.0d-1" is needed instead of 0.5 for Linux build

References/Further reading

Now you are well advised to check out some other resources on this topic, such as:

Older LS-Dyna versions/releases

In R920 everything, so the umat and utan routines, etc. are cramped into the file dyn21.F, which achieves more than 10k lines. Besides that the subroutines are identical so you can implement the code just the same. However, you need a different software setup (Visual Studio, Fortran compiler, e.g. run nmake in "Intel 64 Visual Studio 2010 mode") that fits to the downloaded LS-Dyna object version. Also the interface to UTAN has changed in newer releases, so be aware of that.

todo