*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:

- Write your UMAT subroutine into a file that you can call “umat.f” (the name is not important)
- 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,),);

- 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:

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:

- Elasto-visco-plastic laws: http://code-aster.org/doc/default/en/man_v/v1/v1.03.126.pdf
- Concrete and soil behavior laws: http://code-aster.org/doc/default/en/man_v/v1/v1.03.127.pdf
- Cristalline plasicity laws: http://code-aster.org/doc/default/en/man_v/v1/v1.03.128.pdf
- Laws with Mettalurgy: http://code-aster.org/doc/default/en/man_v/v1/v1.03.129.pdf
- Laws considering Damage: http://code-aster.org/doc/default/en/man_v/v1/v1.03.130.pdf

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

Thomas Helfe says

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

Cyprien says

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!

Thomas Helfer says

Hi Cyprien,

Do not hesitate to contact me in january. I will be glad to help you

Best wishes,

Thomas

Thomas Helfer says

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 !