* fichier : damage_tc_3d.dgibi
************************************************************************
* Section : Mecanique Endommagement
* Section : Mecanique Elastique
************************************************************************
*******************************************************
*  Cas test de l'implantation numérique du modele
*  DAMAGE_TC 3D
*******************************************************
* Développé par :  
*     Benjamin Richard 
* Contact :  
*     Benjamin.Richard@lmt.ens-cachan.fr 
* Institution : 
*     CEA\DEN\DANS\DM2S\SEMT\EMSI
*******************************************************
* Les cas de charges sont entrés : 
* - 1 : Traction monotone
* - 2 : Compression monotone
* - 3 : Traction cyclique
* - 4 : Compression cyclique
* - 5 : Traction/compression cyclique
*******************************************************
* Choix du cas de charge
ncas  = 1;
graph=mot 'N';

*******************************************************

*
*-------------- Options de calcul ---------------------
*
OPTION DIME 3 ELEM CUB8;

*
*-------------- Definition de la geometrie ------------
*
P1 = 0. 0. 0.;
P2 = 1. 0. 0.;
P3 = 1. 1. 0.;
P4 = 0. 1. 0.;
P5 = 0. 0. 1. ;

L1 = P1 DROIT 1 P2 ;
L2 = P2 DROIT 1 P3 ;
L3 = P3 DROIT 1 P4 ;
L4 = P4 DROIT 1 P1 ;

LTOT = L1 ET L2 ET L3 ET L4 ;
SURF1 = SURF LTOT PLANE;
VOLTO =VOLU SURF1 1 TRANS P5;
SURF2= FACE VOLTO 2;

VOLTOT = VOLTO;

*
*------- Définition des conditions aux limites --------
*----------- et des déplacements imposés --------------
*
CL  = BLOQ SURF1 UZ;
CLL = BLOQ P1 'DEPL';
CLL = CLL ET (bloq p2 'UY');
CL1 = BLOQ SURF2 Uz;

D1 = DEPI CL1 1;

*
*----- Définition du modèle ---------------------------
*
MOD1   = MODE VOLTOT MECANIQUE ELASTIQUE ISOTROPE 
         ENDOMMAGEMENT DAMAGE_TC;

*
*----- Paramètres matériaux ---------------------------
* 

* Young's modulus
YG = 36.e9 ;

* Poisson ratio
NUB = 0.20 ;

* Material density
RHOM = 2210 ;

* Fracture energy
Gvalm = 300. ;

* Tension limit strength
ftulm = 3.3e6 ;

* Drop Factor for Peak Tensile Stress
redcm = 1.7e6 ;

* Elastic Limit Compressive Stress
fc01m = -25.e6 ;

* Equi-biaxial Compressive Ratio
rt45m = 1.18 ;

* Compressive Peak Stress
fcu1m = -42.5e6 ;

* Ultimate Limit Strain
extum = -0.015;

* Reference Stress for Plastic Parameter
strpm = -22.e6 ;

* Reference Strain for Plastic Parameter
extpm = -0.001 ;

* Fitting Point 1 (Strain)
ext1m = -0.006 ;

* Fitting Point 1 (Stress)
str1m = -35.e6 ;

* Fitting Point 2 (Strain)
ext2m = -0.008 ;

* Fitting Point 2 (Stress)
str2m = -22.e6 ;

* NCRIT - 1 exponential softening
*       - 2 linear softening
ncrim = 1 ;

* Computation of the element size field
jaco1 = jaco mod1 ;
jaco2 = chan 'RIGIDITE' mod1 jaco1 ;
hlenm = jaco2**(1./3.) ; 


mat1 = mate mod1  YOUN  yg NU nub RHO rhom HLEN hlenm
                   GVAL Gvalm FTUL ftulm REDC redcm FC01 fc01m
                   RT45 rt45m FCU1 fcu1m STRU extum EXTP extpm
                   STRP strpm EXT1 ext1m STR1 str1m EXT2 ext2m
                   STR2 str2m NCRI ncrim ;
*
*----- Définition des cas de charge -------------------
*

SI (EGA ncas 1) ;
LI1  = PROG 0. 1.;
LI2  = PROG 0. 4.0e-4;
LIS1 = PROG 0. PAS 0.02 1.;
FINSI;

SI (EGA ncas 2) ;
LI1  = PROG 0.  1.;
LI2  = PROG 0. -8.0e-3;
LIS1 = PROG 0. PAS 0.02 1.;
FINSI;

SI (EGA ncas 3) ;
LI1  = PROG 0. 1. 2. 3. 4. 5.;
LI2  = PROG 0. 1.5e-4 9.5e-6 2.0E-4 3.0E-5 2.5E-4;
LIS1 = PROG 0. PAS 0.02 5.;
FINSI;

SI (EGA ncas 4) ;
LI1  = PROG 0. 1. 2. 3. 4. 5.;
LI2  = PROG 0. -3.0e-3 -1.5e-4 -5.0E-3 -3.5E-4 -8.0E-3;
LIS1 = PROG 0. PAS 0.02 5.;
FINSI;

SI (EGA ncas 5) ;
LI1  = PROG 0. 1. 2. 3. 4. 5. 6.;
LI2  = PROG 0. 5.0e-4 -3.0e-3 -1.5E-4 -5.0E-3 -3.5E-4 -8.0E-3;
LIS1 = PROG 0. PAS 0.02 2;
FINSI;

EV = EVOL MANU LI1 LI2 ;
CHA1 = CHAR 'DIMP' D1 EV ;

*
*----------- Calcul par l'operateur PASAPAS ------------
*

TAB1 = TABLE ;
TAB1.'BLOCAGES_MECANIQUES' = CL ET CLL ET CL1;
TAB1.'MODELE' = MOD1;
TAB1.'MOVA' = 'RIEN';
TAB1.'CHARGEMENT' = CHA1;
TAB1.'CARACTERISTIQUES' = MAT1;
TAB1.'TEMPS_CALCULES' = LIS1;
TAB1.'MAXITERATION' = 10;

PASAPAS TAB1 ;
 
*
*----------- Courbe effort-deplacement -----------------
*

Si (ega graph 'O');
ev2=@global tab1 CL1 EV fz;
DESS EV2;
finsi;

nn=dime tab1.deplacements;
nn=nn-1;
ee=tab1.variables_internes.nn;

aa = extr ee 'DAMT' 1 1 1 ;
list aa;
err =  abs (aa  - 0.97626) / 0.97626  ;
message ' erreur relative ';
list err;
si (err > 1.e-3);
erreur (5);
finsi; 

fin;