Stress equal to zero using UMAT #656

WeisLeDocto · 2024-12-10T17:48:57Z

WeisLeDocto
Dec 10, 2024

Hi !

I have a .dgibi code that I run with CAST3M 24.1, in which I'm using a UMAT generated by mfront (5.0.0-dev).
It is very basic, adapted from the online tutorial on mfront website.
Leaving out the UMAT, it runs and converges with a MECANIQUE ELASTIQUE ISOTROPE model, with consistent results.

However, when including the UMAT using a 'MECANIQUE' 'ELASTIQUE' 'NON_LINEAIRE' 'UTILISATEUR' model, I'm obtaining zero strain, and therefore zero stress. I'm still getting the correct displacement, though.
As mentioned above, the same code runs when replacing the UMAT model with an isotropic one.
I quickly checked in mfront, using std::cout, and the value of eto + deto in the @Integrator block is always zero.

I'm not sure if I'm not accessing the strain the correct way in the .mfront file, or if I again missed something in the .digibi file.
Any clues about this issue ?

Details

System information as returned by uname -a:
Linux LOrdiDeWeis 6.8.0-49-generic #49~22.04.1-Ubuntu SMP PREEMPT_DYNAMIC Wed Nov 6 17:42:15 UTC 2 x86_64 x86_64 x86_64 GNU/Linux

.dgibi file

*
* \file   test.mfront
* \brief  example of how to use the HomogenizedKelvin behaviour law
* in the Cast3M finite element solver
* \author Antoine Weisrock
* \date   08 / 11 / 2024
*

OPTION 'DIMENSION' 3 ;
OPTION 'MODELISER' 'TRID' ;
OPTION 'ERREUR' 'FATALE' ;
OPTION 'ELEM' CUB8 ;

p0 = 'POINT' 0.0 0.0 0.0 'DENS' 1 ;
px = 'POINT' 1.0 0.0 0.0 'DENS' 1 ;
py = 'POINT' 0.0 1.0 0.0 'DENS' 1 ;
pz = 'POINT' 0.0 0.0 1.0 'DENS' 1 ;

vy0 = 'DROIT' 1 p0 py ;
vy1 = vy0 'PLUS' pz ;
syz0 = vy0 'REGLER' 1 vy1 ;
syz1 = syz0 'PLUS' px ;
vol = syz0 'VOLUME' 1 syz1 ;

* trac (vol et (coul syz0 'ROUG') et (coul syz1 'BLEU'));

coel = 'MOTS' 'YOUN' 
              'NU' 
              'RHO' 
              'ALPH'
              'TALP'
              'TREF'
              'LAMH' 
              'LAM1' 
              'LAM2' 
              'LAM3'
              'LAM4' 
              'LAM5' 
              'SIGS' 
              'M_PO' 
              'THET' ;

params = 'MOTS' 'T' ;

mo = 'MODELISER' vol 
                 'MECANIQUE' 'ELASTIQUE' 'NON_LINEAIRE' 'UTILISATEUR' 
                 'LIB_LOI' 'src/libUmatBehaviour.so' 
                 'FCT_LOI' 'umathomogenizedkelvin_malls'
                 'PARA_LOI' params
                 'C_MATERIAU' coel ;

xyoun = 1000000.0 ;
xnu = 0.2 ;
xrho = 1.0 ;
xalph = 1.0 ;
xtalp = 293.0 ;
xtref = 293.0 ;
xlamh = 11.0 ;
xlam1 = 10.0 ;
xlam2 = 1.0 ;
xlam3 = 1.0 ;
xlam4 = 1.0 ;
xlam5 = 1.0 ;
xsigstd = 1.0 ;
xm_po = 1.0 ;
xtheta = 0.0 ;

ma = 'MATERIAU' MO 
                'YOUN' xyoun 
                'NU' xnu 
                'RHO' xrho 
                'ALPH' xalph
                'TALP' xtalp
                'TREF' xtref
                'LAMH' xlamh 
                'LAM1' xlam1 
                'LAM2' xlam2 
                'LAM3' xlam3 
                'LAM4' xlam4 
                'LAM5' xlam5 
                'SIGS' xsigstd 
                'M_PO' xm_po 
                'THET' xtheta ;

tmax = 10.0 ;
npas = 50 ;
umax = 1.1 ;

lit1 = 'PROG' 0.0 'PAS' (tmax / npas)  tmax ;

cl1 = 'BLOQUE' 'DEPL'  syz0 ;
cl2 = 'BLOQUE' 'UX'  syz1 ;

cle1 = cl1 'ET' cl2 ;

li1 = 'PROG' 0.0 tmax ;
li2 = 'PROG' 0.0 1.0 ;
ev = 'EVOL' 'MANU' 'T' li1 'F(T)' li2 ;
dep2 = 'DEPIMPOSE' cl2 umax ;
cha2 = 'CHARGEMENT' 'DIMP' dep2 ev ;

tek = 293.15;
the1 = MANU 'CHPO' vol 1 'T' tek ;
ev2 = EVOL MANU (PROG 0.0 tmax) (PROG 1.0 1.0) ;
charter = CHAR 'T' the1 ev2 ;

tab1 = 'TABLE' ;
tab1 . 'GRANDS_DEPLACEMENTS' = FAUX ;
tab1 . 'MOVA' = 'MOT' 'RIEN' ;
tab1 . 'TEMPERATURES' = 'TABLE' ;
tab1 . 'TEMPERATURES' . 0 = the1 ;
tab1 . 'VARIABLES_INTERNES' = 'TABLE' ;
tab1 . 'VARIABLES_INTERNES' . 0 = (ZERO mo 'VARINTER');
tab1 . 'BLOCAGES_MECANIQUES' = cle1 ;
tab1 . 'MODELE' = mo ;
tab1 . 'CHARGEMENT' = cha2 'ET' charter ;
tab1 . 'CARACTERISTIQUES' = ma ;
tab1 . 'TEMPS_CALCULES' = lit1 ;
tab1 . 'TEMPS_SAUVES' = lit1 ;
tab1 . 'PRECISION' = 1.e-8;

PASAPAS TAB1 ;

opti donn 5 ;

'REPETER' i ('DIME' tab1 . 'CONTRAINTES') ;
   idx = &i '-' 1 ;
   s = tab1 . 'CONTRAINTES' . idx ;
   d = tab1 . 'DEPLACEMENTS' . idx ;
   'MESSAGE' ('MAXIMUM' ('EXCO' d 'UX')) ' '
             ('MAXIMUM' ('EXCO' d 'UY')) ' '
             ('MAXIMUM' ('EXCO' d 'UZ')) ' '
             ('MAXIMUM' ('EXCO' s 'SMXX')) ' '
             ('MAXIMUM' ('EXCO' s 'SMYY')) ' '
             ('MAXIMUM' ('EXCO' s 'SMZZ')) ' '
             ('MAXIMUM' ('EXCO' s 'SMXY')) ' '
             ('MAXIMUM' ('EXCO' s 'SMXZ')) ' '
             ('MAXIMUM' ('EXCO' s 'SMYZ')) ' ' ;
'FIN' i ;

nb = ('DIME' tab1 . 'CONTRAINTES') '-' 1 ;
psig = 'PROG' ;
peto = 'PROG' ;
'REPETER' i nb ;
   s = tab1 . 'CONTRAINTES' . &i ;
   d = tab1 . 'DEPLACEMENTS' . &i ;
   psig = psig 'ET' ('PROG' ('MAXIMUM' ('EXCO' s 'SMXX'))) ;
   peto = peto 'ET' ('PROG' ('MAXIMUM' ('EXCO' d 'UX'))) ;
'FIN' i ;


evsig = 'EVOL' 'ROUG' 'MANU' 'UX' peto 'SMXX' psig;
'LISTE' evsig ;
'DESSIN' evsig ;

'FIN' ;

.mfront file

// mfront -I $(pwd)/extra-headers -I /usr/include/eigen3 --interface=castem21 --obuild test.mfront

@DSL Default;
@Behaviour HomogenizedKelvin;

@Author Antoine Weisrock;
@Date 08/11/2024;
@Description {This file implements a model based on the homogenization of a 
  stiffness tensor pojected on the Kelvin projectors and expressed as a 
  function of its eigenvalues.
};

/*
@MaterialProperty stress lamh;
@MaterialProperty stress lam1;
@MaterialProperty stress lam2;
@MaterialProperty stress lam3;
@MaterialProperty stress lam4;
@MaterialProperty stress lam5;
*/

@MaterialProperty stress lamh;
@MaterialProperty stress lam1;
@MaterialProperty stress lam2;
@MaterialProperty stress lam3;
@MaterialProperty stress lam4;
@MaterialProperty stress lam5;

/*
@Parameter real sigstd = 1.3;
@Parameter real m = 0.75;
@Parameter real theta = 0.8;
*/

@MaterialProperty real sigstd;
@MaterialProperty real m_po;
@MaterialProperty real theta;

@PhysicalBounds lamh in [0.0: *[;
@PhysicalBounds lam1 in [0.0: *[;
@PhysicalBounds lam2 in [0.0: *[;
@PhysicalBounds lam3 in [0.0: *[;
@PhysicalBounds lam4 in [0.0: *[;
@PhysicalBounds lam5 in [0.0: *[;

@PhysicalBounds sigstd in [0.0: *[;
@Bounds m_po in [-1.0: 1.0];
@Bounds theta in [-3.15: 3.15];

@ModellingHypothesis Tridimensional;
@StrainMeasure Hencky;
// @StrainMeasure GreenLagrange;

@LocalVariable st2tost2<3u, real> stiffness;

@Includes{
  #include "TFEL/Math/KelvinDecomposition.hxx"
}

@InitLocalVariables {
  st2tost2<3u, real> integrated = kelvin_integrated_tensor(
    lamh, lam1, lam2, lam3, lam4, lam5, sigstd, m_po);
  
  st2tost2<3u, real> homogenized = pow_1_over_m(integrated, m_po);
  
  stiffness = rotate(homogenized, theta);
  
}

@Integrator{
  
  sig = stiffness * (eto + deto);
  
}

@TangentOperator {
  
  Dt = stiffness;
  
}

working isotropic version

diff --git a/HomogenizedKelvin.dgibi b/HomogenizedKelvinlinear.dgibi
index a56223e94..f65e5f8d7 100644
--- a/HomogenizedKelvin.dgibi
+++ b/HomogenizedKelvinlinear.dgibi
@@ -41,13 +41,8 @@ coel = 'MOTS' 'YOUN'
               'THET' ;
 
 params = 'MOTS' 'T' ;
-
-mo = 'MODELISER' vol 
-                 'MECANIQUE' 'ELASTIQUE' 'NON_LINEAIRE' 'UTILISATEUR' 
-                 'LIB_LOI' 'src/libUmatBehaviour.so' 
-                 'FCT_LOI' 'umathomogenizedkelvin_malls'
-                 'PARA_LOI' params
-                 'C_MATERIAU' coel ;
+                 
+mo = MODE vol MECANIQUE ELASTIQUE ISOTROPE ;
 
 xyoun = 1000000.0 ;
 xnu = 0.2 ;
@@ -63,24 +58,9 @@ xlam4 = 1.0 ;
 xlam5 = 1.0 ;
 xsigstd = 1.0 ;
 xm_po = 1.0 ;
-xtheta = 0.0 ;
-
-ma = 'MATERIAU' MO 
-                'YOUN' xyoun 
-                'NU' xnu 
-                'RHO' xrho 
-                'ALPH' xalph
-                'TALP' xtalp
-                'TREF' xtref
-                'LAMH' xlamh 
-                'LAM1' xlam1 
-                'LAM2' xlam2 
-                'LAM3' xlam3 
-                'LAM4' xlam4 
-                'LAM5' xlam5 
-                'SIGS' xsigstd 
-                'M_PO' xm_po 
-                'THET' xtheta ;
+xtheta = 0.0 ;            
+
+ma = MATE mo 'YOUN' xyoun 'NU' xnu ;
 
 tmax = 10.0 ;
 npas = 50 ;
@@ -107,13 +87,11 @@ charter = CHAR 'T' the1 ev2 ;
 tab1 = 'TABLE' ;
 tab1 . 'GRANDS_DEPLACEMENTS' = FAUX ;
 tab1 . 'MOVA' = 'MOT' 'RIEN' ;
-tab1 . 'TEMPERATURES' = 'TABLE' ;
-tab1 . 'TEMPERATURES' . 0 = the1 ;
 tab1 . 'VARIABLES_INTERNES' = 'TABLE' ;
 tab1 . 'VARIABLES_INTERNES' . 0 = (ZERO mo 'VARINTER');
 tab1 . 'BLOCAGES_MECANIQUES' = cle1 ;
 tab1 . 'MODELE' = mo ;
-tab1 . 'CHARGEMENT' = cha2 'ET' charter ;
+tab1 . 'CHARGEMENT' = cha2;
 tab1 . 'CARACTERISTIQUES' = ma ;
 tab1 . 'TEMPS_CALCULES' = lit1 ;
 tab1 . 'TEMPS_SAUVES' = lit1 ;

Answered by thelfer

Dec 10, 2024

If you set GRANDS_DEPLACEMENTS to FAUX, I suspect that Cast3M sends the identity as the deformation gradient to the behaviour, hence 0 logarithmic strain, hence 0 stress

View full answer

thelfer · 2024-12-10T20:15:44Z

thelfer
Dec 10, 2024
Maintainer

If you set GRANDS_DEPLACEMENTS to FAUX, I suspect that Cast3M sends the identity as the deformation gradient to the behaviour, hence 0 logarithmic strain, hence 0 stress

1 reply

WeisLeDocto Dec 13, 2024
Author

Thank you, indeed that was the problem

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Stress equal to zero using UMAT #656

Uh oh!

{{title}}

Uh oh!

Replies: 1 comment 1 reply

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

Stress equal to zero using UMAT #656

Uh oh!

WeisLeDocto Dec 10, 2024

Replies: 1 comment · 1 reply

Uh oh!

thelfer Dec 10, 2024 Maintainer

Uh oh!

WeisLeDocto Dec 13, 2024 Author

WeisLeDocto
Dec 10, 2024

Replies: 1 comment 1 reply

thelfer
Dec 10, 2024
Maintainer

WeisLeDocto Dec 13, 2024
Author