* fichier :  umat04.dgibi
*=======================================================================
*
*      CAS TEST DE VALIDATION DE LA PRISE EN COMPTE D'UNE LOI DE
*
*             COMPORTEMENT MECANIQUE NON LINEAIRE EXTERNE
*
*=======================================================================
*
* Modele 'NON_LINEAIRE' 'UTILISATEUR', integrateur specifique UMAT
*
* Test No4 : Modele VISCOPLASTIQUE GATT_MONERIE
*            Reprise du test gatt_3d.dgibi
*
*=======================================================================




*=======================================================================
*
*                       PROCEDURES DEPENDANTES
*
*=======================================================================




*=======================================================================
*
*                      JEU DE DONNEES PRINCIPAL
*
*=======================================================================

* ===============================
* 0 - Options generales de calcul
* ===============================
opti dime 3 mode tridim ;
opti elem cu20 ;
'OPTI' 'TRAC' 'PSC'                  ;
L = 'MOT' LIST ; F = 'MOT' FIN       ;

GRAPH    = FAUX ;
LISTCOUR = VRAI ;

* ============
* 1 - Maillage
* ============
NA  = 1 ;
NR  = 1 ;
H   = 6.2e-3 ;
D   = 5.e-3 ;
R   = D/2. ;
HS2 = H/2. ;

OO  = 0.0  0.0  0.0 ;
A1  = R    0.0  0.0 ;
A2  = R    0.0  HS2 ;
A3  = 0.0  0.0  HS2 ;

LB  = DROIT NR OO A1 ;
LD  = DROIT NA A1 A2 ;
LH  = DROIT NR A2 A3 ;
LG  = DROIT NA A3 OO ;

SU1 = DALL LB  LD  LH  LG  'PLAN' ;
SU2 = SU1 TOUR 30.0 OO A3 ;
VOL1 = VOLU 1 SU1 SU2 ;
ELIM VOL1 (D/1000.0) ;
VOL1 = REGE VOL1 ;

SUA = POIN (COOR 2 VOL1) 'INFERIEUR' (D/1000.0) ;
SUB = POIN (COOR 3 VOL1) 'INFERIEUR' (D/1000.0) ;
SUH = POIN (COOR 3 VOL1) 'SUPERIEUR' (HS2 - (D/1000.0)) ;
LG = VOL1 'POIN' 'DROIT' OO A3 (D/1000.0) ;

* ========================================================
* 2 - Affectation des parametres de la loi de comportement
* ========================================================
TE1  = 1400.0 + 273.0 ;
PO   = 4.52E-2 ;
DF1  = 1.E18 ;
* --------------------------------------------------------------
* 2.1 - AFFECTATION DES PARAMETRES DE MATERIAU : premiere partie
* --------------------------------------------------------------
* Description du module d'Young
* Degre du polynome
IDEGE = 2 ;
* Coefficient du degre 0, coefficient du degre 1, coefficient du degre 2
LEE = PROG 226930.e+06 -1.5399e+07 -9.3597e+03 ;

* Description du module de cisaillement
* Degre du polynome
IDEGG = 2 ;
* Coefficient du degre 0, coefficient du degre 1, coefficient du degre 2
LGG = PROG 85830.e+06 -5.157e+06 -3.747e+03 ;

* Polynome du coefficient de dilatation thermique Polynome(T)/(T-273)
* pour T < Talpha
* Degre du polynome
IDEGAI = 3 ;
* Coefficient du degre 0, coefficient du degre 1, coefficient du degre 2
* coefficient du degre 3
LALI = PROG -2.66e-03 9.802e-6 -2.705e-10 4.391e-13 ;

* Polynome du coefficient de dilatation thermique Polynome(T)/(T-273)
* pour T > Talpha
* Degre du polynome
IDEGAS = 3 ;
* Coefficient du degre 0, coefficient du degre 1, coefficient du degre 2
* coefficient du degre 3
LALS = PROG -3.28e-03 1.179e-05 -2.429e-09 1.219e-12 ;

* Temperature de transition pour les deux polynomes decrivant
* le coefficient de dilatation thermique (en Kelvin)
TALPHA = 923. ;

* Constructions d'objets a l'aide des donnees precedentes

*-- Liste de temperatures (en Kelvin)
TK = 'PROG' 273.15 'PAS' 10. 2603.15 ;
*-- On insere Talpha dans TK
TK = TK 'ET' ('PROG' TALPHA) ;
TK = 'ORDO' TK ;
NTK = 'DIME' TK ;
L273 = 'PROG' NTK * 273.;

*-- Construction des ordonnees : partie temperature
OEET = 'PROG' NTK*('EXTR' LEE 1) ;
'SI' ('>' ('DIME' LEE) 1) ;
  'REPE' BLOC (('DIME' LEE) - 1) ;
    OEET = OEET + (('EXTR' LEE (&BLOC+1))*(TK**(&BLOC)));
  'FIN' BLOC ;
'FINS' ;
'MESS'
  'sur la plage de temperatures' ('MINI' TK) '< T(K) <' ('MAXI' TK) ;
'MESS' 'et pour une porosite nulle' ;
MIOEET = 'MINI' OEET ;
MAOEET = 'MAXI' OEET ;
'MESS' 'module d Young            compris entre' MIOEET 'et' MAOEET ;
'SI' ('<EG' MIOEET 0) ;
  'MESS' 'ERREUR : module d Young negatif ou nul !' ;
  'FIN' ;
'FINS' ;

OGGT = 'PROG' NTK*('EXTR' LGG 1) ;
'SI' ('>' (DIME LGG) 1) ;
  'REPE' BLOC (('DIME' LGG) - 1) ;
    OGGT = OGGT + (('EXTR' LGG (&BLOC+1))*(TK**(&BLOC)));
  'FIN' BLOC ;
'FINS' ;
MIOGGT = 'MINI' OGGT ;
MAOGGT = 'MAXI' OGGT ;
'MESS' 'module de cisaillement    compris entre' MIOGGT 'et' MAOGGT ;
'SI' ('<EG' MIOGGT 0) ;
  'MESS' 'ERREUR : module de cisaillement negatif ou nul !' ;
  'FIN' ;
'FINS' ;

OALIT = 'PROG' NTK*('EXTR' LALI 1) ;
'SI' ('>' (DIME LALI) 1) ;
  'REPE' BLOC (('DIME' LALI) - 1) ;
    OALIT = OALIT + (('EXTR' LALI (&BLOC+1))*(TK**(&BLOC)));
  'FIN' BLOC ;
  OALIT = OALIT / (TK - L273);
'FINS' ;
*
OALST = 'PROG' NTK*('EXTR' LALS 1) ;
'SI' ('>' (DIME LALS) 1) ;
  'REPE' BLOC (('DIME' LALS) - 1) ;
    OALST = OALST + (('EXTR' LALS (&BLOC+1))*(TK**(&BLOC)));
  'FIN' BLOC ;
  OALST = OALST / (TK - L273);
'FINS' ;
*
ALTI = IPOL TALPHA TK OALIT ;
ALTS = IPOL TALPHA TK OALST ;
*
NTKI = TK 'MASQ' 'INFERIEUR' TALPHA 'SOMM' ;
OALIT = 'EXTR' OALIT ('LECT' 1 'PAS' 1 NTKI) ;
OALST = 'EXTR' OALST ('LECT' (NTKI+1) 'PAS' 1 NTK) ;
OALT = OALIT 'ET' OALST ;
*
* On tronque ALPHA_moy en dessous de 293K
NT293 = TK 'MASQ' 'INFERIEUR' 293. 'SOMM';
XEALTT = 'EXTR' TK ('LECT' (NT293+1) 'PAS' 1 NTK) ;
YEALTT = 'EXTR' OALT ('LECT' (NT293+1) 'PAS' 1 NTK) ;
* On remplace la partie enlevée par une IPOL entre
* ALPHA_moy(T=0) et ALPHA_moy(T=293)
AL0 = ('EXTR' LALI 1) / (0. - 273);
AL293 =  (('EXTR' LALI 1) +
         (('EXTR' LALI 2)*293.) +
         (('EXTR' LALI 3)*293.*293.) +
         (('EXTR' LALI 4)*293.*293.*293.)) / (293. - 273.);
XEALTN = (PROG 0. 293.) ET XEALTT;
YEALTN = (PROG AL0 AL293) ET YEALTT;
*
MIOALT = 'MINI' YEALTN ;
MAOALT = 'MAXI' YEALTN ;
'MESS' 'coefficient de dilatation compris entre' MIOALT 'et' MAOALT ;
'SI' ('<' MIOALT 0) ;
  'MESS' 'ERREUR : coefficient de dilatation negatif !' ;
'FINS' ;
*
ERRALT=100.*(ALTS-ALTI)/ ALTS ;
'SI' ('>' ERRALT 1.D-02) ;
  'MESS' 'ATTENTION : a la temperature de transition :' TALPHA 'K' ;
  'MESS' 'les deux valeurs du coefficient de dilatation thermique' ;
  'MESS' 'different de' ERRALT '%' ALTI ALTS ;
'FINS' ;
*
*-- Construction des evolutions obtenues pour une porosite nulle
EEET = 'EVOL' 'MANU' 'T' TK     'YOUN'  OEET   ;
EGGT = 'EVOL' 'MANU' 'T' TK     'G'     OGGT   ;
EALT = 'EVOL' 'MANU' 'T' XEALTN 'ALPH'  YEALTN ;
* --------------------------------------------------------------
* 2.2 - AFFECTATION DES PARAMETRES DE MATERIAU : deuxieme partie
* --------------------------------------------------------------
PORY = 2.5 ;
PORG = 2.25 ;
PORA = 1. ;

* Constructions d'objets a l'aide des donnees precedentes

*-- Liste de porosites ( )
PP = 'PROG' 0. 'PAS' 0.001 0.3 ;
NPP = 'DIME' PP ;

*-- Construction des ordonnees : partie porosite
OEEP = ('PROG' NPP*1.) - (PORY*PP) ;
MIOEEP = 'MINI' OEEP ;
'SI' ('<EG' MIOEEP 0) ;
  'MESS' 'ATTENTION !' ;
  'MESS' 'sur la plage de porosite :' ('MINI' PP) '< P <' ('MAXI' PP) ;
  'MESS ' 'module d Young negatif ou nul !' ;
'FINS' ;

OGGP = ('PROG' NPP*1.) - (PORG*PP) ;
MIOGGP = 'MINI' OGGP ;
'SI' ('<EG' MIOGGP 0) ;
  'MESS' 'ATTENTION !' ;
  'MESS' 'sur la plage de porosite :' ('MINI' PP) '< P <' ('MAXI' PP) ;
  'MESS ' 'module de cisaillement negatif ou nul !' ;
'FINS' ;

OALP = ('PROG' NPP*1.) - (PORA*PP) ;
MIOALP = 'MINI' OALP ;
'SI' ('<EG' MIOALP 0) ;
  'MESS' 'ATTENTION !' ;
  'MESS' 'sur la plage de porosite :' ('MINI' PP) '< P <' ('MAXI' PP) ;
  'MESS ' 'coefficient de dilatation negatif ou nul !' ;
'FINS' ;

*-- Construction des nuages
PP1 = 'EXTR' PP 1 ;

EE = 'NUAG' 'PORO'*'FLOTTANT' 'YOUN'*'EVOLUTION'
            ('EXTR' PP 1)     (EEET*('EXTR' OEEP 1)) ;
'REPE' BLOC (NPP-1) ;
  EEI = 'NUAG' 'PORO'*'FLOTTANT'     'YOUN'*'EVOLUTION'
               ('EXTR' PP (&BLOC+1)) (EEET*('EXTR' OEEP (&BLOC+1))) ;
  EE = EE 'ET' EEI ;
'FIN' BLOC ;

ETT1 = 'EVOL' 'MANU' 'T' TK 'NU' ('PROG' ('DIME'TK)*1.) ;
NU = 'NUAG' 'PORO'*'FLOTTANT' 'NU'*'EVOLUTION'
            ('EXTR' PP 1)
        ((-1.*ETT1)+(0.5*EEET*('EXTR' OEEP 1)/(EGGT*('EXTR' OGGP 1)))) ;
'REPE' BLOC (NPP-1) ;
  NUI = 'NUAG' 'PORO'*'FLOTTANT'     'NU'*'EVOLUTION'
               ('EXTR' PP (&BLOC+1))
        ((-1.*ETT1)+(0.5*(EEET*('EXTR' OEEP (&BLOC+1)))/
          (EGGT*('EXTR' OGGP (&BLOC+1))))) ;
  NU = NU 'ET' NUI ;
'FIN' BLOC ;

AL = 'NUAG' 'PORO'*'FLOTTANT' 'ALPH'*'EVOLUTION'
            ('EXTR' PP 1)     (EALT*('EXTR' OALP 1)) ;
'REPE' BLOC (NPP-1) ;
  ALI = 'NUAG' 'PORO'*'FLOTTANT'     'ALPH'*'EVOLUTION'
               ('EXTR' PP (&BLOC+1)) (EALT*('EXTR' OALP (&BLOC+1))) ;
  AL = AL 'ET' ALI ;
'FIN' BLOC ;

*-- Construction de la table contenant les informations
TA = 'TABL' ;
TA . 'YOUN' = EE ;
TA . 'NU'   = NU ;
TA . 'ALPH' = AL ;
* ---------------------------------------------------------------
* 2.3 - AFFECTATION DES PARAMETRES DE MATERIAU : troisieme partie
* ---------------------------------------------------------------
  LCOMP = 'MOTS'  'RHO'   'R'     'DG0'   'DG'    'DGCR'
                  'K0'    'EXKP'  'KPAF'  'AP'    'BP'  'QP'
                  'K1'    'M1'    'Q1'    'N1'    'K2'  'M2'  'Q2'  'N2'
                  'OMEG'  'OMCR'  'Q'     'H'     'BETA'
                  'DYN1'  'DYCR'  'DYN2'  'DYN3'
                  'K'     'A'     'Q3'    'N3'
                  'ADEN'  'KGON'  'POR0'  'BUMI'  'EFIS'
                  'CR1'   'CR2'   'CR3'   'CR';

LDONN = 'PROG' 10950.    8.314   7e-6   8.2e-6 56e-6
   1e-9      3.48     1.8e-27  4.277    1.212    500000.
7.42e-14       -2.     377000.    1.    1.e-43    2.  462000. 8.
 47350.4    47350.     0.189    600.    1.
2.06e-11    50000.      1615.    40.
2.44e-19  5.19e-36      8000.     1.
               -0.21e-02 7.75463e-15 4.52e-02 1.728e9 3.2e-11
               400. 0.094E-2 0.01E-2 0.   ;

'REPE' BLOC ('DIME' LCOMP) ;
  XX = 'EXTR' LDONN &BLOC ;
  MOTI = 'EXTR' LCOMP &BLOC ;
  TA . MOTI = XX ;
'FIN' BLOC ;
*
*-- Calcul de KP pour fluage primaire UO2
 TA.'KPUO' = TA.'K0' * (TA.'DG' ** TA.'EXKP') ;
* fluage primaire rendu significatif pour UO2 standard
 TA.'KPUO'   = 1E3 * TA.'KPUO' ;

* ==================================================================
* 3 - Definition du modele et affectation des proprietes de materiau
* ==================================================================
LCMAT07 = 'MOTS' 'YOUN' 'NU  ' 'RHO ' 'ALPH'
 'R   '   'DG  '
 'K1  '   'M1  '   'Q1  '   'N1  '   'K2  '   'M2  '   'Q2  '   'N2  '
 'OMEG'   'Q   '   'H   '   'BETA'   'K   '   'A   '   'Q3  '   'N3  '
 'ADEN'   'KGON'   'POR0'   'BUMI'   'EFIS'   'DG0 '
 'CR1 '   'CR2 '   'CR3 '   'CR  '   'KP  '   'AP  '   'BP  '   'QP  '
 'VISQ'   'TYPE'   'COMP'   'DYN ' ;

* MODIF
* Correction 25/11/2003
* Coherence Gibiane umat04.dgibi avec Esope umagat.eso :
* 6 variables internes supplementaires correspondant aux composantes du
* tenseur de deformations inelastiques

************** LCVAR07 = MOTS 'EPSE' 'PORO' 'BU  ' 'TETA' 'EPSD' 'EPSG' ;
LCVAR07 = MOTS 'EPSE' 'PORO' 'BU  ' 'TETA' 'EPSD' 'EPSG'
               'ENXX' 'ENYY' 'ENZZ' 'GNXY' 'GNXZ' 'GNYZ' ;

* MODIF

LCPAR07 = MOTS 'T   ' 'DFIS' ;

mo_util = MODE VOL1 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE'
                    'NON_LINEAIRE' 'UTILISATEUR'
                    'NUME_LOI' 7 'C_MATERIAU' LCMAT07
                    'C_VARINTER' LCVAR07 'PARA_LOI' LCPAR07 ;

MATREE_u = MATE mo_util 'YOUN' (TA.'YOUN') ;
PP_u = MANU 'CHML' mo_util 'T   ' TE1 'PORO' PO 'RIGIDITE' ;
EE_u = VARI 'NUAG' mo_util MATREE_u PP_u ;

MATR11_u = MATE mo_util 'YOUN' EE_u    'NU  ' 0.3
                        'RHO ' 10950.0 'ALPH' 0.0 ;

MATR12_u = MATE mo_util 'R   '  (TA.'R' )     'DG  '  (TA.'DG') ;

MATR13_u = MATE mo_util 'K1  '  (TA.'K1')     'M1  '  (TA.'M1')
                        'Q1  '  (TA.'Q1')     'N1  '  (TA.'N1')
                        'K2  '  (TA.'K2')     'M2  '  (TA.'M2')
                        'Q2  '  (TA.'Q2')     'N2  '  (TA.'N2') ;

MATR14_u = MATE mo_util 'OMEG'  (TA.'OMEG')   'Q   '  (TA.'Q')
                        'H   '  (TA.'H')      'BETA'  (TA.'BETA')
                        'K   '  (TA.'K')      'A   '  (TA.'A')
                        'Q3  '  (TA.'Q3')     'N3  '  (TA.'N3')  ;

MATR15_u = MATE mo_util 'ADEN'  (TA.'ADEN')   'KGON'  (TA.'KGON')
                        'POR0'  (TA.'POR0')   'BUMI'  (TA.'BUMI')
                        'EFIS'  (TA.'EFIS')                      ;

MATR16_u = MATE mo_util 'DG0 '  (TA.'DG0')
                        'CR1 '  (TA.'CR1')    'CR2 '  (TA.'CR2')
                        'CR3 '  (TA.'CR3')    'CR  '  (TA.'CR')  ;

MATR17_u = MATE mo_util 'KP  '  (TA.'KPUO')   'AP  '  (TA.'AP')
                        'BP  '  (TA.'BP')     'QP  '  (TA.'QP')
'VISQ' 0.  'TYPE'  0.   'COMP'  0.            'DYN '  0.         ;

* TYPE = 0. combustible UO2            sinon combustible AFA3GLAA
* COMP = 0. combustible compressible   sinon incompressible
* DYN  = 0. couplage statique          sinon dynamique

MATR1_u = MATR11_u ET MATR12_u ET MATR13_u ET MATR14_u ET MATR15_u ET
          MATR16_u ET MATR17_u ;

* ========================================================
* 4 - Definition des conditions aux limites et chargements
* ========================================================
CLZB = 'BLOQ' UZ SUB ;
CLZH = 'BLOQ' UZ SUH ;
CLXG = 'BLOQ' UX LG ;
CLYG = 'BLOQ' UY SUA ;
CLT  = CLZB 'ET' CLZH 'ET' CLXG 'ET' CLYG ;

TMIL = 1.E6 ;
TFIN = 2.E6 ;
*
RHO0 = 1.0 - TA.'POR0' ;
AKBU = TA.'EFIS'*270./238./10950./RHO0*DF1 ;
* CE QUI SUIT EST VRAI CAR LE FLUX DE FISSIONS EST CONSTANT
TBU = (TA.'BUMI' / AKBU) - 1. ;
TTBU = ('ENTI' (TBU/10000.)) * 10000. ;
*------------------------
* 4.1 - Instants calcules
*------------------------
'SI' LISTCOUR ;
*--------------
LIST1 = 'PROG' 0   'PAS' 50  300  'PAS' 100  1000 'PAS' 1000  10000
                   'PAS' 10000. TTBU
                   TBU (TTBU+10000.) (TTBU+20000.)
                   'PAS' 10000. TMIL ;
'SINON' ;
*--------
LIST1 = 'PROG' 0   'PAS' 50  300  'PAS' 100  1000 'PAS' 1000  10000
                   'PAS' 10000. TTBU
                   TBU (TTBU+10000.) (TTBU+20000.)
                   'PAS' 10000. TMIL
                   'PAS' 100. (TMIL + ((TFIN - TMIL) / 200.))
                   'PAS' 1000. (TMIL + ((TFIN - TMIL) / 20.))
                   'PAS' 2000. (TMIL + ((TFIN - TMIL) / 4.))
                   'PAS' 10000. TFIN ;
'FINSI' ;
*--------
* -------------------------------
* 4.2 - Chargement en temperature
* -------------------------------
CHTEMP = 'MANU' 'CHPO' VOL1 1 'T' 1.     ;
EVT = 'EVOL' 'MANU' ('PROG' 0.  (2.*TFIN))
                    ('PROG' TE1 TE1) ;
CHARTEMP = 'CHAR' 'T' CHTEMP EVT        ;
* --------------------------------------
* 4.3 - Chargement en densite de fission
* --------------------------------------
CHFISS = 'MANU' 'CHPO' VOL1 1 'DFIS' 1.     ;
EVF = 'EVOL' 'MANU' ('PROG' 0.  (2.*TFIN))
                    ('PROG' DF1 DF1) ;
CHARFISS = 'CHAR' 'DFIS' CHFISS EVF      ;
* --------------------------------------
* 4.4 - Chargement en deplacement impose
* --------------------------------------
DEPZ = 'DEPI' CLZH 1 ;
DELTAZ = HS2/200. ;
EVEPS = 'EVOL' 'MANU'
   ('PROG' 0. TMIL   TFIN (2.*TFIN))
   ('PROG' 0. DELTAZ 0.    0.) ;
CHARMECA = 'CHAR' 'DIMP' DEPZ EVEPS ;
*
CHARTOT  = CHARMECA 'ET' CHARTEMP 'ET' CHARFISS ;

* =========================================
* 5 - Initialisation des variables internes
* =========================================
VAR00_u = ZERO mo_util 'VARINTER' ;
VAR01_u = MANU 'CHML' mo_util 'PORO' PO
                              'TYPE' 'VARIABLES INTERNES' 'STRESSES' ;
VAR0_u  = VAR00_u + VAR01_u ;

* ==========================
* 6 - Resolution par PASAPAS
* ==========================
TABU = TABLE                           ;
TABU.'MODELE' = mo_util                ;
TABU.'CARACTERISTIQUES' = MATR1_u      ;
TABU.'BLOCAGES_MECANIQUES' = CLT       ;
TABU.'CHARGEMENT' = CHARTOT            ;
TABU.'VARIABLES_INTERNES'= TABLE       ;
TABU.'VARIABLES_INTERNES' . 0 = VAR0_u ;
TABU.'TEMPS_CALCULES' = LIST1          ;
TABU.'TEMPS_SAUVES'   = LIST1          ;
*
TMASAU                   = TABLE       ;
TABU . 'MES_SAUVEGARDES' = TMASAU      ;
TMASAU .'DEFTO' = VRAI                 ;
TMASAU .'DEFIN' = VRAI                 ;
*
TEMPS 'ZERO' ;
PASAPAS TABU ;
*
CPUext = TEMPS 'NOEC' ;

* ================================================
* 7 - Post-traitement : verification des resultats
* ================================================
AP    = TA.'AP'        ;
BP    = TA.'BP'        ;
N1    = TA.'N1'        ;
N2    = TA.'N2'        ;

* Facteurs multiplicatifs dus au dopage par le Chrome
WC1 = 'TANH' ( (TA.'CR' - TA.'CR2') / TA.'CR3' ) ;
WC1 = 1. + ( 0.5 * TA.'CR1' * (1. + WC1) )       ;
CV  = 180. / PI                                  ;
WC2 = 1. - ( 'COS' (CV * TA.'DG' / TA.'DG0') )   ;
WC2 = 2. * (TA.'DG0'**TA.'M2') * WC2             ;

BUMI = TA.'BUMI' ;
AKEVD = AKBU*TA.'KGON' ;
ADEN = TA.'ADEN' ;
KGON = TA.'KGON' ;
* --------------------------------------------------------------
* 7.1 - Calcul de la cte AAAA pour le calcul de la def. de dens.
* --------------------------------------------------------------
BUMAX0=60.D0*BUMI ;
CRIT=1.D-10 ;
'REPE' BLOC 100 ;
  BUMAX = BUMI* ('EXP' (1. - (ADEN/(KGON*BUMAX0)))) ;
  BUMAX = (0.2*BUMAX) + (0.8*BUMAX0) ;
  TEST='ABS' ((BUMAX-BUMAX0)/BUMAX0) ;
  'SI' ('<' TEST CRIT) ;
    'QUIT' BLOC ;
  'FINS' ;
  BUMAX0=BUMAX ;
'FIN' BLOC ;
AAAA = (-1.D0*RHO0*(ADEN-(KGON*BUMAX))) /
                ((1.D0+ADEN)*(LOG(BUMAX/BUMI)))  ;
*
* ----------------------------
* 7.2 - Controle des resultats
* ----------------------------
SS_u  = TABU . 'CONTRAINTES'                ;
VV_u  = TABU . 'VARIABLES_INTERNES'         ;
IN_u  = TABU . 'DEFORMATIONS_INELASTIQUES'  ;
NCONT = ('DIME' (TABU . 'CONTRAINTES')) - 1 ;
*
T1     = 0.        ;
SMZZ1  = 0.        ;
FF1    = TA.'POR0' ;
ERMAX1 = 0.        ;
ERMAX2 = 0.        ;
*
LTT_u    = PROG ;
LDFF_u   = PROG ;
LEVFF_u  = PROG ;
LEVP3_u  = PROG ;
LEVTH_u  = PROG ;
LEVCALCu = PROG ;
LSMZZ_u  = PROG ;
*
'REPE' BLOC_u NCONT ;
*
* -------------------------------------
* VITESSE DE DEFORMATION THEORIQUE EVTH
* -------------------------------------
  T2   = 'EXTR' (TABU.'TEMPS_SAUVES') (&BLOC_u + 1);
  EVTH = DELTAZ / HS2 / TMIL                       ;
  'SI' (> T2 (TMIL+0.1))                           ;
    EVTH = -1. * DELTAZ / HS2 / (TFIN - TMIL)      ;
  'FINS'                                           ;
* ----------------------------------------------------
* SMZZ EST PREPONDERANTE DEVANT LES AUTRES CONTRAINTES
* ----------------------------------------------------
  SMZZ = 'MAXI' ('EXCO' (SS_u . &BLOC_u) 'SMZZ') ;
  SSM = SMZZ / 3. ;
  SSEQ = 'ABS' SMZZ ;
  SSPRIM1 = -1. * SMZZ / 3 ;
  SSPRIM2 = -1. * SMZZ / 3 ;
  SSPRIM3 =  2. * SMZZ / 3 ;
* ---------------------------------------------
* DEFORMATION VISCO-PLASTIQUE EQUIVALENTE EVPEQ
* ---------------------------------------------
  EVXX  = 'MAXI' ('EXCO' VV_u.&BLOC_u 'ENXX')                       ;
  EVYY  = 'MAXI' ('EXCO' VV_u.&BLOC_u 'ENYY')                       ;
  EVZZ  = 'MAXI' ('EXCO' VV_u.&BLOC_u 'ENZZ')                       ;
  EVH   = (EVXX+EVYY+EVZZ) / 3.                                     ;
  EVPEQ = ((EVXX - EVH)**2) + ((EVYY - EVH)**2) + ((EVZZ - EVH)**2) ;
  EVPEQ = ((2./3.) * EVPEQ)**0.5                                    ;
* ------------------------------------------
* VITESSE DE DEFORMATION VISCO-PLASTIQUE EVP
* ------------------------------------------
*-- fluage primaire
*
'SI' ('<' EVPEQ 1E-10) ; EVPEQ = 1E-10 ; 'FINSI' ;
*
  EVP01 = 1.5 * TA.'KPUO' *
  ('EXP' (-1.*TA.'QP'/(TA.'R'*TE1))) *
  ( (SSEQ**(AP-1)) / (EVPEQ**BP) ) * SSPRIM1 ;
*
  EVP02 = 1.5 * TA.'KPUO' *
  ('EXP' (-1.*TA.'QP'/(TA.'R'*TE1))) *
  ( (SSEQ**(AP-1)) / (EVPEQ**BP) ) * SSPRIM2 ;
*
  EVP03 = 1.5 * TA.'KPUO' *
  ('EXP' (-1.*TA.'QP'/(TA.'R'*TE1))) *
  ( (SSEQ**(AP-1)) / (EVPEQ**BP) ) * SSPRIM3 ;
*
*-- fluage secondaire (2 mecanismes)
*
  PP = 'MAXI' ('EXCO' VV_u . &BLOC_u 'PORO') ;

  A1  = (N1* (PP**(-1./N1) - 1.)**(-2.*N1/(N1+1.)))      ;
  B1  = (1. + (2.*PP/3.))/ ((1. - PP) **(2.*N1/(N1+1.))) ;
  AB1 = (9. * A1) / ( (4.*B1) + A1 )                     ;
  BA1 = B1 / ( B1 + (A1/4.) )                            ;
*
  GSSEQ = ((BA1 * (SSEQ**2)) + (AB1 * (SSM**2))) ** 0.5  ;

  A2 = (N2* (PP**(-1./N2) - 1.)**(-2.*N2/(N2+1.)))       ;
  B2 = (1. + (2.*PP/3.))/ ((1. - PP) **(2.*N2/(N2+1.)))  ;
*
  EVP11 = 0.5 * (TA.'K1' * (TA.'DG'**TA.'M1') * WC1 *
  ('EXP' (-1.*TA.'Q1'/(TA.'R'*TE1))) *
  (((A1*((1.5*SSM)**2)) + (B1*(SSEQ ** 2)))**((N1-1.)/2.)) *
  ((A1*1.5*SSM) + (3.*B1*SSPRIM1))) ;
*
  EVP21 = 0.5 * (TA.'K2' * WC2 * ('EXP' (-1.*TA.'Q2'/(TA.'R'*TE1))) *
  (((A2*((1.5*SSM)**2)) + (B2*(SSEQ ** 2)))**((N2-1.)/2.)) *
  ((A2*1.5*SSM) + (3.*B2*SSPRIM1))) ;
*
  EVP12 = 0.5 * (TA.'K1' * (TA.'DG'**TA.'M1') * WC1 *
  ('EXP' (-1.*TA.'Q1'/(TA.'R'*TE1))) *
  (((A1*((1.5*SSM)**2)) + (B1*(SSEQ ** 2)))**((N1-1.)/2.)) *
  ((A1*1.5*SSM) + (3.*B1*SSPRIM2))) ;
*
  EVP22 = 0.5 * (TA.'K2' * WC2 * ('EXP' (-1.*TA.'Q2'/(TA.'R'*TE1))) *
  (((A2*((1.5*SSM)**2)) + (B2*(SSEQ ** 2)))**((N2-1.)/2.)) *
  ((A2*1.5*SSM) + (3.*B2*SSPRIM2))) ;
*
  EVP13 = 0.5 * (TA.'K1' * (TA.'DG'**TA.'M1') * WC1 *
  ('EXP' (-1.*TA.'Q1'/(TA.'R'*TE1))) *
  (((A1*((1.5*SSM)**2)) + (B1*(SSEQ ** 2)))**((N1-1.)/2.)) *
  ((A1*1.5*SSM) + (3.*B1*SSPRIM3))) ;
*
  EVP23 = 0.5 * (TA.'K2' * WC2 * ('EXP' (-1.*TA.'Q2'/(TA.'R'*TE1))) *
  (((A2*((1.5*SSM)**2)) + (B2*(SSEQ ** 2)))**((N2-1.)/2.)) *
  ((A2*1.5*SSM) + (3.*B2*SSPRIM3))) ;
*
  PSI1 = TA.'K1' * WC1 * (TA.'DG'**TA.'M1')/(N1+1.) *
  ('EXP' (-1.*TA.'Q1'/(TA.'R'*TE1))) *
  (((A1*((1.5*SSM)**2)) + (B1*(SSEQ ** 2)))**((N1+1.)/2.)) ;
*
  PSI2 = TA.'K2' * WC2/(N2+1.) * ('EXP' (-1.*TA.'Q2'/(TA.'R'*TE1))) *
  (((A2*((1.5*SSM)**2)) + (B2*(SSEQ ** 2)))**((N2+1.)/2.)) ;
*
  EVPC1 = 0.5 * TA.'BETA' / TA.'H' * TA.'Q' * TA.'OMEG' *
  (GSSEQ**(-1.*TA.'Q' - 2.)) *
      (1. -
  (('TANH' ((TE1 - (TA.'OMEG' * (GSSEQ**(-1.*TA.'Q'))))  / TA.'H'))**2))
   * (PSI2 - PSI1) * ( ((AB1*SSM)/3.) +  (1.5*BA1*SSPRIM1) ) ;
*
  EVPC2 = 0.5 * TA.'BETA' / TA.'H' * TA.'Q' * TA.'OMEG' *
  (GSSEQ**(-1.*TA.'Q' - 2.)) *
      (1. -
  (('TANH' ((TE1 - (TA.'OMEG' * (GSSEQ**(-1.*TA.'Q'))))  / TA.'H'))**2))
   * (PSI2 - PSI1) * ( ((AB1*SSM)/3.) +  (1.5*BA1*SSPRIM2) ) ;
*
  EVPC3 = 0.5 * TA.'BETA' / TA.'H' * TA.'Q' * TA.'OMEG' *
  (GSSEQ**(-1.*TA.'Q' - 2.)) *
      (1. -
  (('TANH' ((TE1 - (TA.'OMEG' * (GSSEQ**(-1.*TA.'Q'))))  / TA.'H'))**2))
   * (PSI2 - PSI1) * ( ((AB1*SSM)/3.) +  (1.5*BA1*SSPRIM3) ) ;
*
  TETA = 0.5 * TA.'BETA' *
      (1. +
      ('TANH' ((TE1 - (TA.'OMEG' * (GSSEQ**(-1.*TA.'Q')))) / TA.'H'))) ;
*
  EVP1 = EVP01 + ((1. - TETA) * EVP11) + (TETA * EVP21) + EVPC1 ;
  EVP2 = EVP02 + ((1. - TETA) * EVP12) + (TETA * EVP22) + EVPC2 ;
  EVP3 = EVP03 + ((1. - TETA) * EVP13) + (TETA * EVP23) + EVPC3 ;
* -------------------------------------
* VITESSE DE DEFORMATION ELASTIQUE EVEL
* -------------------------------------
  DSMZZ = (SMZZ - SMZZ1) / (T2 - T1) ;
  EVEL  = DSMZZ / ('MAXI' EE_u)      ;
* ----------------------------------
* VITESSE DE DEFORMATION IRRADIATION
* ----------------------------------
  EVIR1 = 1.5 * TA.'A' * DF1 * (SSEQ **(TA.'N3' - 1.)) *
       ('EXP' (-1.*TA.'Q3'/(TA.'R'*TE1))) * SSPRIM1 ;
  EVIR2 = 1.5 * TA.'A' * DF1 * (SSEQ **(TA.'N3' - 1.)) *
       ('EXP' (-1.*TA.'Q3'/(TA.'R'*TE1))) * SSPRIM2 ;
  EVIR3 = 1.5 * TA.'A' * DF1 * (SSEQ **(TA.'N3' - 1.)) *
       ('EXP' (-1.*TA.'Q3'/(TA.'R'*TE1))) * SSPRIM3 ;
* ------------------------------------
* VITESSE DE DEFORMATION DENSIFICATION
* ------------------------------------
* CE QUI SUIT EST VRAI CAR LE FLUX DE FISSIONS EST CONSTANT
  BU = AKBU * T2 ;
  'SI' (BU '>' BUMI) ;
    NUMER = (RHO0 * RHO0 * AKEVD) - (RHO0 * AAAA * AKEVD) +
            (RHO0 * AAAA * AKEVD * ('LOG' (BU / BUMI))) -
            (RHO0 * AAAA / T2) ;
    DENOM = (RHO0 + (AAAA * ('LOG' (BU / BUMI)))) *
            (RHO0 + (AAAA * ('LOG' (BU / BUMI)))) ;
    EVD = -1. * (AKEVD - (NUMER/DENOM)) / 3. ;
  'SINO' ;
    EVD = 0. ;
  'FINS' ;
* ---------------------------------
* VITESSE DE DEFORMATION GONFLEMENT
* ---------------------------------
* CE QUI SUIT EST VRAI CAR LE FLUX DE FISSIONS EST CONSTANT
  EVG = AKEVD/3. ;
* -------------------------------
* VITESSE DE DEFORMATION SELON ZZ
* -------------------------------
  EVCALC = ((1.+(TA.'K'*DF1))*EVP3) + EVEL + EVIR3 +
                       EVD + EVG ;
  ECART1 = (EVCALC - EVTH) / EVTH ;
  ERMAX1 = 'MAXI' ('PROG' ('ABS' ECART1) ERMAX1) ;
* -------------------
* VITESSE DE POROSITE
* -------------------
  FF = 'MAXI' ('EXCO' (VV_u . &BLOC_u) 'PORO') ;
  DFF = (FF - FF1) / (T2 - T1) ;
  EVFF = (3.*EVD) +
         ((1. - FF) * (1.+(TA.'K'*DF1))*(EVP1+EVP2+EVP3)) +
         ((1. - FF) * (EVIR1 + EVIR2 + EVIR3)) ;
  ECART2 = ('ABS' (DFF - EVFF)) / FF ;
  ERMAX2 = 'MAXI' ('PROG' ('ABS' ECART2) ERMAX2) ;

  LTT_u = LTT_u 'ET' ('PROG' T2) ;
  LEVP3_u = LEVP3_u 'ET' ('PROG' EVP3) ;
  LEVCALCu = LEVCALCu 'ET' ('PROG' EVCALC) ;
  LEVTH_u = LEVTH_u 'ET' ('PROG' EVTH) ;
  LDFF_u = LDFF_u 'ET' ('PROG' DFF) ;
  LEVFF_u = LEVFF_u 'ET' ('PROG' EVFF) ;
  LSMZZ_u   = LSMZZ_u   'ET' ('PROG' SMZZ)   ;


  FF1   = FF ;
  SMZZ1 = SMZZ ;
  T1 = T2 ;

'FIN' BLOC_u ;

Message 'CPU avec modele externe = ' CPUext 'centiemes de seconde' ;

'SI' ( ERMAX1 '<EG' 0.05 );
   'ERRE' 0;
'SINO';
   'MESS' 'ERREUR MAXIMALE VITESSE DEFORMATION :'
           ERMAX1 '> 0.05 ' ;
   'ERRE' 5;
'FINS';

'SI' ( ERMAX2 '<EG' 0.05 );
   'ERRE' 0;
'SINO';
   'MESS' 'ERREUR MAXIMALE POROSITE :'
           ERMAX2 '> 0.05 ' ;
   'ERRE' 5;
'FINS';

*-----------------------------------------------------------------------
* 8 - Post-traitement : visualisations graphiques si demandees

'SI' GRAPH ;

  'TITR' 'Tridimensionnel : Vitesses de Deformation'         ;
  EVEVP_u  = 'EVOL' 'MANU' LTT_u 'Temps'  LEVP3_u  'Vitesse' ;
  EVEVC_u  = 'EVOL' 'MANU' LTT_u 'Temps'  LEVCALCu 'Vitesse' ;
  EVEVTH_u = 'EVOL' 'MANU' LTT_u 'Temps'  LEVTH_u  'Vitesse' ;
  TAD_u = TABLE                                              ;
  TAD_u . 1 = 'TIRC'                                         ;
  TAD_u . 3 = 'TIRR'                                         ;
  TAD_u . 'TITRE' = TABLE                                    ;
  TAD_u . 'TITRE' . 1 = 'Viscoplastique Analytique'          ;
  TAD_u . 'TITRE' . 2 = 'Totale Calculee'                    ;
  TAD_u . 'TITRE' . 3 = 'Totale Analytique'                  ;
  'DESS' (EVEVP_u 'ET' EVEVC_u 'ET' EVEVTH_u) TAD_u 'LEGE'   ;
*
  'TITR' 'Tridimensionnel : Vitesses de variation de la porosite' ;
  EVDFF_u  = 'EVOL' 'MANU' LTT_u 'Temps' LDFF_u  'Vitesse'        ;
  EVEVFF_u = 'EVOL' 'MANU' LTT_u 'Temps' LEVFF_u 'Vitesse'        ;
  TAF_u = TABLE                                                   ;
  TAF_u . 2 = 'TIRR'                                              ;
  TAF_u . 'TITRE' = TABLE                                         ;
  TAF_u . 'TITRE' . 1 = 'Calculee'                                ;
  TAF_u . 'TITRE' . 2 =  'Analytique'                             ;
  'DESS' (EVDFF_u 'ET' EVEVFF_u) TAF_u 'LEGE'                     ;
*
  'TITR' 'Tridimensionnel : Contrainte'                         ;
  EVSMZZ_u = 'EVOL' 'MANU' LTT_u 'Temps'  LSMZZ_u 'Contrainte'  ;
  'DESS' EVSMZZ_u                                               ;
'FINS' ;
'FINS' ;

FIN ;