from basics to Advanced

How to add your own material behavior law in Code_Aster?

by 8 Comments

Note: This topic is quite advanced, so I will only highlight some of the things you should know to point you in the right direction… I’ll try to detail more in the future, so make sure you are in the newsletter to get informed about the latest articles ;-)

That said…

There are several ways to add new materials laws in Code_Aster!

As usual, some of those ways are way more complex than others…

First, if you know how to code in Fortran, Python, C and how to implement new code in Code_Aster, then you might consider developping a material law directly inside the DEFI_MATERIAU Operator of Code_Aster…
That’s the most difficult way to do it…and I definitely don’t recommend you to do that

The 2 easier ways are to use other kind of code which possess an interface with Code_Aster and allows you to “simply” plug-in your material law into Code_Aster

The 2 kinds of interfaces you can use are:

  • The UMAT subroutine (If you are familiar with Abaqus, that’s the same thing used here)
  • The generator of material behavior law MFRONT

1 The UMAT subroutine

1.1 What is the UMAT subroutine?

UMAT is a subroutine FORTRAN which is fully customizable to allow you to integrate your own material behavior laws.

Here’s a tutorial I created for you:

You will need to download those files to follow the tutorial:
– The HDF file containing the asterstudy case and the mesh
– The UMAT code
Download the files

If you are familiar with Abaqus, you surely have heard about it.

The header of the UMAT subroutine looks like that:
SUBROUTINE UMAT(
STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
RPL,DDSDDT,DRPLDE,DRPLDT,
STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

Each of those parameters are either input or output parameters and you have to know how exactly what they represent and how to calculate them numerically…

You will find a description of how the UMAT subroutine works in Code_Aster here:
https://www.code-aster.org/V2/doc/default/en/man_u/u2/u2.10.01.pdf

1.2 How to run your integrate your UMAT subroutine in Code_Aster?

Let’s take an example of UMAT subroutine taken from the Abaqus manual:

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C ======================================================================
C      INCLUDE 'ABA_PARAM.INC'
C
C     ROUTINE UMAT EXEMPLE CITEE DANS LA DOCUMENATION ABAQUS
C     COMPORTEMENT VISCOELASTIQUE LINEAIRE
C

      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*80 CMNAME
      REAL*8 STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),
     2 DDSDDT(NTENS),DRPLDE(NTENS),
     3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
      REAL*8 DSTRES(6),D(3,3)
C
C  EVALUATE NEW STRESS TENSOR
C
      EV = 0.D0
      DEV = 0.D0
      DO 1001 K1=1,NDI
         EV = EV + STRAN(K1)
         DEV = DEV + DSTRAN(K1)
 1001 CONTINUE
C
      TERM1 = .5D0*DTIME + PROPS(5)
      TERM1I = 1.D0/TERM1
      TERM2 = (.5D0*DTIME*PROPS(1)+PROPS(3))*TERM1I*DEV
      TERM3 = (DTIME*PROPS(2)+2.D0*PROPS(4))*TERM1I
C
      DO 1002 K1=1,NDI
         DSTRES(K1) = TERM2+TERM3*DSTRAN(K1)
     1     +DTIME*TERM1I*(PROPS(1)*EV
     2     +2.D0*PROPS(2)*STRAN(K1)-STRESS(K1))
         STRESS(K1) = STRESS(K1) + DSTRES(K1)
 1002 CONTINUE
C
      TERM2 = (.5D0*DTIME*PROPS(2) + PROPS(4))*TERM1I
      I1 = NDI
      DO 1003 K1=1,NSHR
         I1 = I1+1
         DSTRES(I1) = TERM2*DSTRAN(I1)+
     1     DTIME*TERM1I*(PROPS(2)*STRAN(I1)-STRESS(I1))
         STRESS(I1) = STRESS(I1)+DSTRES(I1)
 1003 CONTINUE
C
C  CREATE NEW JACOBIAN
C
      TERM2 = (DTIME*(.5D0*PROPS(1)+PROPS(2))+PROPS(3)+
     1  2.D0*PROPS(4))*TERM1I
      TERM3 = (.5D0*DTIME*PROPS(1)+PROPS(3))*TERM1I
      DO 1004 K1=1,NTENS
         DO 1005 K2=1,NTENS
            DDSDDE(K2,K1) = 0.D0
 1005    CONTINUE
 1004 CONTINUE
C
      DO 1006 K1=1,NDI
         DDSDDE(K1,K1) = TERM2
 1006 CONTINUE
C
      DO 1007 K1=2,NDI
         N2 = K1-1
         DO 1008 K2=1,N2
            DDSDDE(K2,K1) = TERM3
            DDSDDE(K1,K2) = TERM3
 1008    CONTINUE
 1007 CONTINUE
      TERM2 = (.5D0*DTIME*PROPS(2)+PROPS(4))*TERM1I
      I1 = NDI
      DO 1009 K1=1,NSHR
         I1 = I1+1
         DDSDDE(I1,I1) = TERM2
 1009 CONTINUE
C
C  TOTAL CHANGE IN SPECIFIC ENERGY
C
      TDE = 0.D0
      DO 1010 K1=1,NTENS
         TDE = TDE + (STRESS(K1)+.5D0*DSTRES(K1))*DSTRAN(K1)
 1010 CONTINUE
C
C  CHANGE IN SPECIFIC ELASTIC STRAIN ENERGY
C
      TERM1 = PROPS(1) + 2.D0*PROPS(2)
      DO 1011 K1=1,NDI
         D(K1,K1) = TERM1
 1011 CONTINUE
      DO 1012 K1=2,NDI
         N2 = K1-1
         DO 1013 K2=1,N2
            D(K1,K2) = PROPS(1)
            D(K2,K1) = PROPS(1)
 1013    CONTINUE
 1012 CONTINUE
      DEE = 0.D0
      DO 1014 K1=1,NDI
         TERM1 = 0.D0
         TERM2 = 0.D0
         DO 1015 K2=1,NDI
            TERM1 = TERM1 + D(K1,K2)*STRAN(K2)
            TERM2 = TERM2 + D(K1,K2)*DSTRAN(K2)
 1015    CONTINUE
         DEE = DEE + (TERM1+.5D0*TERM2)*DSTRAN(K1)
 1014 CONTINUE
      I1 = NDI
      DO 1016 K1=1,NSHR
         I1 = I1+1
         DEE = DEE + PROPS(2)*(STRAN(I1)+.5D0*DSTRAN(I1))*DSTRAN(I1)
 1016 CONTINUE
      SSE = SSE + DEE
      SCD = SCD + TDE - DEE

      END

Here’s the steps to follow:

  1. Write your UMAT subroutine into a file that you can call “umat.f” (the name is not important)
  2. Create the dynamic librarie associated to your file by compiling it with as_run

(as_run is the software which allows you to run code_aster from command line)

  • cd path/of/your/study
  • as_run –make_shared -o libumat.so umat.f
  • Specify the following info in ASTK so that it can link to your new material library
  • Give it the /path/of/your/study/libumat.so
  • type “nom”
  • LU = 0 (the unit of the file)
  • Check the “D” because it is an input data

Your library file will be copied in the study directory automatically
4)In the DEFI_MATERIAU command, you need to define all your material field data in the keyword UMAT
Example:

ACIER[0]=DEFI_MATERIAU(ELAS=_F(E=YOUNG_Pa,
                               NU=POISSON,
                               ALPHA=11.8e-6),
            UMAT=_F( LISTE_COEF = (LAMBD_Pa ,MU_Pa ,LAMBB_Pa,MUB_Pa,NUB)),
            ECRO_LINE=_F(D_SIGM_EPSI=pente_Pa, SY=SY_Pa,),);
  1. In the .comm file, you need to find the keyword “COMPORTEMENT” (Behavior) in the STAT_NON_LINE, DYNA_NON_LINE or SIMU_POINT_MAT Operator

Add the following sub-keywords:

  • RELATION=’UMAT’
  • LIBRAIRIE=’libumat.so’ (you can also specify the absolute path /path/of/your/study/libumat.so)
  • NOM_ROUTINE=’umat_’ (or just ‘umat’)

Example:

RESU=SIMU_POINT_MAT(
         COMPORTEMENT=_F(RELATION=UMAT,
                      NB_VARI=1,
                      LIBRAIRIE='libumat1.so',
                      NOM_ROUTINE='umat_'),
         NEWTON=_F(MATRICE='TANGENTE',
                   REAC_ITER=1,),
         CONVERGENCE=_F(ITER_GLOB_MAXI=10,),
         MATER      = ACIER[imat],
         ARCHIVAGE = _F(LIST_INST = temps_ar),
         INCREMENT=_F(LIST_INST=temps,),
         EPSI_IMPOSE=_F(EPXX=eps_imp[0], EPYY=eps_imp[1], EPZZ=eps_imp[2],
                        EPXY=eps_imp[3], EPXZ=eps_imp[4], EPYZ=eps_imp[5]),
            INFO=2,  );

2 Generating material behavior laws with MFRONT

Whereas UMAT allows you to have full control over your material behavior, you still have to know in depth how coding a material behavior law works
(How to calculate the stress tensor, the internal variables, the tangent operator, the integration points, etc…)

By the way, if you need to learn FORTRAN, check my video tutorials here, I am going from the very beginning on how to run a simple Fortran Code and make a simple program:

[Fortran Tuto 1] Hello World !

 

MFRONT is simplier than that

2.1 What is MFRONT?

It is a generator of material behavior laws which can help you to create simply behavior laws by inputting formulas and equations close to the physical equations without having to bother about how those laws need to be solved numerically.
(You can choose between the algorithms already proposed in MFRONT to integrate your behavior)

If you worry about the validation of this code, then I can assure you that you shouldn’t because MFRONT was developped by the CEA in France for many years to be used in highly sensible nuclear areas so it works for sure
(Doesn’t mean you don’t need to test and validate what YOU developped with MFRONT ;-) )

2.2 How to use MFRONT with Code_Aster?

The following part will give you some “keys” but you’ll definitely need to investigate more in details how it work…

MFRONT is coupled with Code_aster and installed with it like an external software library which can interact with Code_Aster through a special interface.

Here are a few steps you can take:

2.2.1 Write your material behavior law into a “Law.mfront” file

For example let’s take this AnisoLemaitre Law provided in Code_Aster:

@Parser    Implicit;
@Behaviour AnisoLemaitre;
@Author    jmp;
@Date      03/2014;
@Description{
 "Modèle de comportement élasto-visqueux"
 "META_LEMA_ANI avec prise en compte de la "
 "métallurgie pour les tubes de gaine du crayon "
 "combustible."
 "Documentation du Code Aster R4.04.05"
 "http:www.code-aster.org"
 "Modèle utilisé par le département MMC pour"
 "décrire les gaines de crayon combustible en"
 "situation d'APRP."
}

@Includes{
#include<TFEL/Material/Hill.hxx>
#include<TFEL/Material/Lame.hxx>
}

@OrthotropicBehaviour;
@Algorithm NewtonRaphson_NumericalJacobian;
@Theta 1.;
@Epsilon 1.e-10;

@MaterialProperty real young;
young.setGlossaryName("YoungModulus");
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");

@MaterialProperty real a[3];
@MaterialProperty real m[3];
@MaterialProperty real pn[3];
@MaterialProperty real Q[3];
@MaterialProperty real M1[6];
@MaterialProperty real M3[6];

@StateVariable real p;      /* Equivalent viscoplastic strain */

@AuxiliaryStateVariable real seq;
seq.setGlossaryName("HillStress");
@AuxiliaryStateVariable real svi[3];

@LocalVariable stress lambda;
@LocalVariable stress mu;
@LocalVariable tfel::math::st2tost2<N,real> H;
@LocalVariable real T_ ;
@LocalVariable real invn[3], f[3], gamma[3], sv[3]  ;

// variables de commande aster
@ExternalStateVariable real SECH,HYDR,IRRA,NEUT1,NEUT2,CORR,ALPHPUR,ALPHBETA;

@IsTangentOperatorSymmetric true;
@TangentOperator{
  using namespace tfel::material::lame;
  if((smt==ELASTIC)||(smt==SECANTOPERATOR)){
    computeElasticStiffness<N,Type>::exe(Dt,lambda,mu);
  } else if(smt==CONSISTENTTANGENTOPERATOR){
    StiffnessTensor Hooke;
    Stensor4 Je;
    computeElasticStiffness<N,Type>::exe(Hooke,lambda,mu);
    getPartialJacobianInvert(Je);
    Dt = Hooke*Je;
  } else {
    return false;
  }
}

@InitLocalVariables{
  using namespace tfel::material::lame;
  lambda = computeLambda(young,nu);
  mu = computeMu(young,nu);
  // proportion en phase alpha en milieu de pas de temps
  const real Z = min(max(ALPHPUR  + theta*dALPHPUR+
			 ALPHBETA + theta*dALPHBETA,0.),1.);
  if (Z >= 0.99) {
    f[0]=1. ;
  } else if (Z >= 0.9) {
    f[0] = (Z-0.9)/0.09 ;
  } else {
    f[0] = 0. ;
  }
  if (Z >= 0.1) {
    f[2]=0. ;
  } else if (Z >= 0.01) {
    f[2] = (0.1-Z)/0.09 ;
  } else {
    f[2] = 1. ;
  }
  if (Z >= 0.99) {
    f[1]=0. ;
  } else if (Z >= 0.9) {
    f[1] = 1.0-(Z-0.9)/0.09 ;
  } else if (Z >= 0.1) {
    f[1] = 1.0 ;
  } else if (Z >= 0.01) {
    f[1] = 1.0-(0.1-Z)/0.09 ;
  } else {
    f[1] = 0. ;
  }

  // Temperature Aster en Celsius
  T_ = 273.0 + T + theta * dT ;
  for(unsigned short i=0;i!=3;++i){
    invn[i] = 1.0 / pn[i] ;
    gamma[i] = a[i]* exp(Q[i]/T_*invn[i]) ;
  }
  // correspondance M aster (repere x,y,z) et H
  real M[6];
  if (Z >= 0.99) {
    for(unsigned short i=0;i!=6;++i){
      M[i]=M1[i] ;
    }
  } else if (Z >= 0.01) {
    for(unsigned short i=0;i!=6;++i){
      M[i]=Z*M1[i]+(1.-Z)*M3[i] ;
    }
  } else {
    for(unsigned short i=0;i!=6;++i){
      M[i]=M3[i] ;
    }
  }
  const real H_F = 0.5*( M[0]+M[1]-M[2]);
  const real H_G = 0.5*(-M[0]+M[1]+M[2]);
  const real H_H = 0.5*( M[0]-M[1]+M[2]);
  const real H_L = 2.0*M[3];
  const real H_M = 2.0*M[4];
  const real H_N = 2.0*M[5];
  H = hillTensor<N,real>(H_F,H_G,H_H,
			 H_L,H_M,H_N);
}

@ComputeStress{
  sig = lambda*trace(eel)*Stensor::Id()+2*mu*eel;
}

@Integrator{
  const real sigeq = sqrt(sig|H*sig);
  real p_=p+theta*dp ;
  real sigv =  0. ;
  real pm[3];
  real dpn[3] ;
  for(unsigned short i=0;i!=3;++i){
    pm[i]  = (p_ > 0.) ? pow(p_,m[i])         : 0.;
    dpn[i] = (dp > 0.) ? pow((dp/dt),invn[i]) : 0. ;
    sv[i]=gamma[i]*pm[i]*dpn[i] ;
    sigv += f[i]*sv[i] ;
  }
  // normale à l'écoulement
  Stensor  n(0.);
  real inv_sigeq(0);
  if(sigeq > 1.e-10*young){
    inv_sigeq = 1/sigeq;
    n  = (H*sig)*inv_sigeq;
  }
  feel += dp*n-deto;
  fp    = (sigeq-sigv)/young;
}

@UpdateAuxiliaryStateVars{
  seq = sqrt(sig|H*sig);
  for(unsigned short i=0;i!=3;++i){
    svi[i]=sv[i] ;
  }
}

2.2.2 Use the Code_Aster Operator CREA_LIB_MFRONT to generate the dynamic library

You can write it like that into your .comm file:

CREA_LIB_MFRONT(UNITE_MFRONT=38,UNITE_LIBRAIRIE=39)

2.2.3 Provide the coefficients required by the law in your DEFI_MATERIAU

MATF=DEFI_MATERIAU(ELAS_META=_F(E = 80000.,
                               NU = 0.35,
                               F_ALPHA = 0.,
                               C_ALPHA = 0.,
                               PHASE_REFE = 'FROID',
                               EPSF_EPSC_TREF = 0.,),
                      MFRONT=_F( LISTE_COEF = (80000.,0.35,2.39,0.22,9.36,0.0,0.77E-04,0.99E-04,4.39,2.96,6.11,19922.8,21023.7,6219.,0.4414,1.,0.714,0.75,0.75,0.75,1.,1.,1.,0.75,0.75,0.75,
                            ),),)
                            
#a_0  =2.39,       ## F1_A
#a_1  =0.22,       ## F2_A
#a_2  =9.36,       ## C_A
#m_0  =0.0,        ## F1_M
#m_1  =0.77E-04,   ## F2_M
#m_2  =0.99E-04,   ## C_M
#pn_0 =4.39,       ## F1_N
#pn_1 =2.96,       ## F2_N
#pn_2 =6.11,       ## C_N
#Q_0 =19922.8,    ## F1_Q
#Q_1 =21023.7,    ## F2_Q,
#Q_2 =6219.,      ## C_Q
#M1_0 =0.4414,     ## F_M11
#M1_1 =1.,         ## F_M22
#M1_2 =0.714,      ## F_M33
#M1_3 =0.75,       ## F_M12
#M1_4 =0.75,       ## F_M13
#M1_5 =0.75,       ## F_M23
#M3_0 =1.,         ## C_M11
#M3_1 =1.,         ## C_M22
#M3_2 =1.,         ## C_M33
#M3_3 =0.75,       ## C_M12
#M3_4 =0.75,       ## C_M13
#M3_5 =0.75,       ## C_M23

2.2.4 Implement the MFRONT law in the COMPORTEMENT KEYWORD OF YOUR STAT_NON_LINE Operator

Example:

VF=STAT_NON_LINE(MODELE=MO,
                CHAM_MATER=CMF,
                EXCIT=(_F(CHARGE=CHPC,FONC_MULT=FMULT,),
                       _F(CHARGE=CH_L,),),

                COMPORTEMENT=_F(RELATION='MFRONT',
                                NOM_ROUTINE='asteranisolemaitre',
                                UNITE_LIBRAIRIE=39,
                                RESI_INTE_MAXI=1e-10
                        ),
               INCREMENT=_F(LIST_INST = L_INST),
               NEWTON=_F(MATRICE = 'TANGENTE',
                         REAC_ITER = 1),
               CONVERGENCE=_F(RESI_GLOB_RELA = 1.E-6,
                              ITER_GLOB_MAXI = 30))

More than 20 material laws have been implemented and tested with Code_Aster, such as:

For the moment, the material behavior laws which are usable in Code_Aster are usable for 3D, 2D, Shell, THM and Joint Elements.

2.3 Where can you find MFRONT Tutorials?

Check the documentation here:
http://tfel.sourceforge.net/documentations.html

There is a good tutorial available here, but it’s in French:
http://tfel.sourceforge.net/documents/tutoriel/tutoriel.pdf

Comments

  1. Hi Cyprien,

    Thank your for this kind introduction on how to add non linear behaviours in Code_Aster.

    As the main developer of MFront, I really appreciate your comments about it and particularly on how we strive to provide an industrial strenght product (MFront has now 10 000 “public” unit tests).

    If I may, I would add a few comments.

    First about Code_Aster’ UMAT interface:
    – though strongly inspired by Abaqus/Standard, Code_Aster’ UMAT interface is not totally compatible with Abaqus’ one. You may run into various issues when porting from Abaqus to Code_Aster. One example among the others: implementing orthotropic behaviours is very different.
    – UMAT is limited to small strain behaviours (however your behaviour can be embedded in GROT_GDEP or GDEF_LOG finite strain strategies. The latter works well remarkly well for metals). You may not be able to implement your own hyperelasto-viscoelastic behaviour for example.
    – All in all, I would say that the UMAT interface is great when you do have a reasonnably simple isotropic small strain behaviour that you already have implemented for Abaqus and that you want to port it to Code_Aster (which is a very welcomed move).

    Concerning MFront:
    – You gave a somehow complex example involving phase transition. That’s great to illustrate that one can implement such a complex behaviour quite easily.
    – The given implementation is typical of MFront 2.x versions. In version 3.x, the concept of bricks has been introduced, which strongly reduces the MFront implementation. For version 3.2, which is not yet available in official Code_Aster/Salomé Méca packages, the StandardElastoViscoplastic brick (tfel.sourceforge.net/StandardElastoViscoPlasticityBrick.html) may cover many users’ needs without requiring any C++ code.

    Finally, let me thank you again for your work to advertise and explain the usage of Code_Aster.

    With kind regards,

    Thomas Helfer

    Reply

    • Thank you Thomas for those precisions! I wasn’t aware of the limitations of UMAT with Aster and it was limited to small strains only… which makes the usage of MFRONT even more necessary and beneficial.
      I have still a lot to learn about MFRONT and its correct usage and I was planning to write something more detailed and simple for beginners in the future.
      Very glad to be in contact with you Thomas!

      Reply

        • BTW, I didnot mean that the UMAT interface was limited to isotropic small strain behaviours, but that compatibility with Abaqus was limited to this class of behaviours. You may write orthotropic finite strain behaviours using Code_Aster’UMAT interface (not quite sure if everything is available) but you will have to use some undocumented piece of code. In particular, the definition of the consistent tangent operator is very complicated with the SIMO_MIEHE’ implementation of the finite strain resolution: It took me weeks to set it up properly with the help of Code_Aster developpers !

          Reply

  3. Hi Cyprien
    Thank you for sharing your knowledge about interfacing Umat with Code_Aster.
    Do you know how to pass effectively PNEWDT<1 information to code aster in order to reduce time step?
    Best regards
    Christian

    Reply

  4. I need to develop a UMAT program to introduce FGM properties and to see how the material is mutilated until refraction.

    Reply

  5. Thank you for the walkthrough. I am trying to convert the .f file to a .so file as directed, however the aba_param.inc file can not be opened. I’ve read this is a file from ABAQUS but I do not have access to ABAQUS to link the file. Is there another option to convert my UMAT.f into a .so file?

    the exact error reads: Cannot open included file ‘ABA_PARAM.INC’

    This UMAT file does work in abaqus but I was hoping to use it in code-aster since it is an open-source software. I am fairly new to this and would appreciate your feedback.

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *

View Cyprien Rusu's profile on LinkedIn

Find videos about FEA

Read previous post:
How to determine Weibull Law Parameters with Python? (Part 2)

[latexpage] This article is a second article on the Weibull Law which explains how to use Python to calculate the...

Close