Sources Esope | Cast3M

cendom
C CENDOM    SOURCE    PV        17/12/08    21:15:33     9660           
       SUBROUTINE CENDOM(wrk52,wrk53,wrk54,WRK6,WRK7,WRK8,WRK9,
     1  NVARI, T0,T1,TREF,IB,IGAU,ifourb,ncourb,iecou)
 
C ENDOM     SOURCE    NICOLAS   99/02/10    21:16:58     3476
*      SUBROUTINE ENDOM(WRK0,WR00,WRK1,WRK6,WRK7,WRK8,WRK9,WTRAV,
*     1  NSTRS,NXMATT,ICARA,INPLAS,NVARI,PRECIS,MFR,IFOURB,KERRE,
*     2   NNVARI,NYOG,NYNU,NYRHO,NYALFA,NNKX,NYKX,NCOURB,NYN,NYM,
*     3         NYKK,T0,T1,TREF,ITHHER,CMATE,N2EL,N2PTEL,IB,IGAU,
*     4         EPAIST,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,CRIGI)
C
c NSTRSS &lt;= NSTRS , MFR1 &lt;- MFR , XCARB &lt;- XCAR
C ==================================================================
C CE SOUS-PROGRAMME EST APPELE DANS "PLAST".
C C'EST UN INTEGRATEUR GENERAL DES LOIS DE COMPORTEMENT DES
C MATERIAUX ENDOMMAGEABLES EN REGIME ANISOTHERME.
C
C ENTREES:
C -------
C WRK0,WR00,WRK1,WRK6,WRK7,WRK8,WRK9,WTRAV = SEGMENTS DE TRAVAIL
C NSTRSS          = NBR. DE COMPOSANTES DES CONTR. OU DES DEFORM.
C SIG0(NSTRSS)    = CONTR. AU DEBUT DU PAS D'INTEGRATION
C DEPST(NSTRSS)   = INCREMENT DES DEFORM. ELASTIQUES
C NVARI          = NBR. DE VARIABLES INTERNES
C VAR0(NVARI)    = VARIABLES INTERNES AU DEBUT DU PAS D'INTEGRATION
C                  CE TABLEAU CONTIENT DANS L'ORDRE:
C                - LA DEFORM. PLASTIQUE CUMULEE
C                - LES VARIABLES INTERNES PILOTANT LES EQUATIONS DU
C                  MODELE
C                - LES DEFORM. PLASTIQUES
C MFR1            = INDICE DE LA FORMULATION MECANIQUE; SEULEMENT
C                  MASSIF OU COQUE POUR LES MATERIAUX ENDOMMAGEABLES
* CMATE = NOM DU MATERIAU
* N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
* N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
* MFR1   = NUMERO DE LA FORMULATION
* IFOURB= OPTION DE CALCUL
* IB    = NUMERO DE L ELEMENT COURANT
* IGAU  = NUMERO DU POINT COURANT
* EPAIST= EPAISSEUR
* NBPGAU= NBRE DE POINTS DE GAUSS
* MELE  = NUMERO DE L ELEMENT FINI
* NPINT = NBRE DE POINTS D INTEGRATION
* NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
* NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
* SECT  = SECTION
* LHOOK = TAILLE DE LA MATRICE DE HOOKE
C ICARA          = NBR. DE CARACT. GEOMETRIQUES DES ELEMENTS FINIS
C XCARB(ICARA)    = CARACT. GEOMETRIQUES DES ELEMENTS FINIS
C INPLAS         = INDICE DU MODELE DECRIVANT LE MATERIAU ENDOMMAGEABLE
C NYKX          = DIMENSION DE YKX
C YKX(NYKX)   = TABLEAU CONTENANT LES COURBES DE TRACTION(EPS,SIGMA)
C                  A DES TEMPERATURES DONNEES (Ti);IL EST DE LA FORME
C                  SUIVANTE :
C                  (T1,X1,F(X1,T1),X2,F(X2,T1),..,T2,Y1,F(Y1,T2),Y2,..)
C NNKX           = NBRE DE COURBES DE TRACTION
C NKX            = TABLEAU CONTENANT LE NBRE DE PTS POUR CHAQUE COURBE
C                  DE TRACTION (EX:SI NKX(1)=4,LA COURBE DE TRACTION
C                  A LA TEMPERATURE T1 SERA DEFINIE PAR 4 POINTS)
C PRECIS         = PRECISION DES ITERATIONS INTERNES
C TO             = TEMPERATURE INITIALE
C T1             = TEMPERATURE FINALE
C Tref           = TEMPERATURE DE REFERENCE
C NYOG          = 2*NBRE DE POINTS DE LA COURBE T-->YOG(T)
C YOG(NYOG)    = TABLEAU CONTENANT LES VALEURS CONNUES DU MODULE
C                  D'YOUNG E A DES TEMPERATURES PRECISES T
C NYNU           = 2*NBRE DE POINTS DE LA COURBE T-->YNU(T)
C YNU(NYNU)      = TABLEAU CONTENANT LES VALEURS CONNUES DU COEFF. DE
C                  POISSON YNU A DES TEMPERATURES PRECISES T
C NYALFA          = 2*NBRE DE POINTS DE LA COURBE T-->YALFA(T)
C YALFA(NYALFA)    = TABLEAU CONTENANT LES VALEURS CONNUES DU COEFF. DE
C                  DILAT. THERMIQUE YALFA A DES TEMPERATURES PRECISES T
C NYN          = 2*NBRE DE POINTS DE LA COURBE T-->YN(T)
C YN(NYN)    = TABLEAU CONTENANT LES VALEURS CONNUES DU SEUIL
C                  D'ENDOMT YN A DES TEMPERATURES PRECISES T
C NYM            = 2*NBRE DE POINTS DE LA COURBE T-->YM(T)
C YM(NYM)        = TABLEAU CONTENANT LES VALEURS CONNUES DE L'ENDOMT
C                  CRITIQUE YM A DES TEMPERATURES PRECISES T
C NYKK          = 2*NBRE DE POINTS DE LA COURBE T-->YKK(T)
C YKK(NYKK)    = TABLEAU CONTENANT LES VALEURS CONNUES DE LA DEFORM.
C                  A RUPTURE YKK A DES TEMPERATURES PRECISES T
C
C SORTIES:
C -------
C NCOURB     = DIMENSION DE TRUC
C TRUC       = TABLEAU DE DIMENSION LA PLUS GRANDE PARMI CELLES
C              DEFINIES DANS NKX;IL CONTIENDRA LA COURBE DE
C              TRACTION A LA TEMPERATURE T,OBTENUE GRACE AU SOUS-
C              PROGRAMME TRACTT.
C SIGF(NSTRSS)= CONTR. A LA FIN DU PAS D'INTEGRATION
C VARF(NVARI)= VARIABLES INTERNES A LA FIN DU PAS D'INTEGRATION
C DEFP(NSTRSS)= INCREMENT DES DEFORM. PLASTIQUES A LA FIN DU PAS
C              D'INTEGRATION
C KERRE      = INDICE QUI REGIT LES ERREURS
C            = 77 SI LA DEFORM. PLAST. CUMULEE ENDOMMAGEE (2IEME VAR.
C                 INT.) EST EN DEHORS DE LA COURBE DE TRACTION, DS.
C                 LE CAS DE L'ECROUISSAGE ET DE L'ENDOMM. ISOTROPES.
C                 CECI PEUT SE PRODUIRE SUITE A L'APPEL A "CRIDAM"
C            = 99 SI LA FORMULATION MECANIQUE N'EST PAS DISPONIBLE
C                  POUR LE MODELE CONSIDERE OU S'IL Y A INCOMPATIBILITE
C                 ENTRE MFR1 ET IFOURB
C ==================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC DECHE
      SEGMENT IECOU
*      COMMON/IECOU/NYOG,NYNU,NYALFA,icow1,NYN,NYM,NYKK,
      INTEGER NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
     1  icow8,icow9,icow10,icow11,NYRHO,icow13,NNKX,NYKX,icow16,
*    .  NYALF1,NYBET1,NYR,  NYA,   NYRHO,NSIGY,  NNKX,  NYKX,  IND,
     2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
*     .  NSOM,NINV,NINCMA,NCOMP,JELEM,LEGAUS,INAT,NCXMAT,
     3 icow25,icow26,icow27,icow28,icow29,icow30,icara,
*     .  LTRAC,MFR,IELE,NHRM,NBNN,NBELEM,ICARA,
     4  icow32,icow33, NSTRSS,MFR1,NBGMAT,NELMAT,icow38,
*     .  LW2,   NDEF,  NSTRSS,MFR1,NBGMAT,NELMAT,MSOUPA,
     5  icow39,icow40,icow41,NNVARI,icow43,icow44
*     .  NUMAT1,LENDO,NBBB,NNVARI,KERR1,MELEME
        INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
     .          icow51,icow52,icow53,icow54,icow55,icow56
     .          icow57,icow58
      ENDSEGMENT
 
*
      SEGMENT WRK6
        REAL*8 BB(NSTRS,NNVARI),R(NSTRS),XMU(NSTRS)
        REAL*8 S(NNVARI),QSI(NNVARI),DDR(NSTRS),BBS(NSTRS)
        REAL*8 SIGMA(NSTRS),SIGGD(NSTRS)
      ENDSEGMENT
*
      SEGMENT WRK7
        REAL*8 F(NCOURB,2),W(NCOURB),TRUC(NCOURB)
      ENDSEGMENT
*
      SEGMENT WRK8
        REAL*8 DD(NSTRS,NSTRS),DDV(NSTRS,NSTRS),DDINV(NSTRS,NSTRS)
        REAL*8 DDINVp(NSTRS1,NSTRS1)
      ENDSEGMENT
*
      SEGMENT WRK9
        REAL*8 YOG(NYOG),YNU(NYNU),YALFA(NYALFA),YSMAX(NYSMAX)
        REAL*8 YN(NYN),YM(NYM),YKK(NYKK),YALFA1(NYALF1)
        REAL*8 YBETA1(NYBET1),YR(NYR),YA(NYA),YKX(NYKX),YRHO(NYRHO)
        REAL*8 SIGY(NSIGY)
        INTEGER NKX(NNKX)
      ENDSEGMENT
 
*
      DIMENSION DEPS(6)
*
C
C ------------------------------------------------
C DEFINITION DU NBR. MAXIMUM D'ITERATIONS INTERNES
C ------------------------------------------------
      DATA ITMAX/30/
c
      NXMATT = NMATT
C
      ZEERO=0.D0
      PRECI=PRECIS
      PRECIG=-2.D0*PRECIS
      PRECG=ABS(PRECIG)
      INDIC=0
      IALPH=0
      T=T1
      YUNGV=0.D0
      XNUV=0.D0
C
C ================================================================
C MAJ. POUR LE CAS ELASTIQUE ET TEST POUR SAVOIR SI ON A PLASTIFIE
C ================================================================
      CALL ZDANUL(DEFP,NSTRSS)
      ALPH=1.D0
      ALPH0=0.D0
      ALPH1=1.D0
      X=VAR0(2)
 
      CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,NYN,
     1     YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T0,
     2                  XMAT1,NXMATT,TRUC,NCOURB)
      CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,NYN,
     1     YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T1,
     2                  XMAT2,NXMATT,TRUC,NCOURB)
      ALFA0=XMAT1(4)
      ALFA1=XMAT2(4)
      CTEPS=ALFA0*(T0-TREF)-ALFA1*(T1-TREF)
C
C==================================================
C  CALCUL DE L INCREMENT DE CONTRAINTES
C  DEPST EST L INCREMENT DE DEFORMATIONS ELASTIQUES
C==================================================
C
      CALL ELAST1(1,IFOUR,VAR0,NNVARI,XMAT2,NXMATT,YUNGV,XNUV,
     1          XCARB,ICARA,MFR1,NSTRSS,DD,DDV,KERRE,INDIC,ITHHER)
C
      CALL ZDANUL(DSIGT,NSTRSS)
      DO 13 I=1,NSTRSS
      DO 13 J=1,NSTRSS
             DSIGT(I)=DSIGT(I)+DD(I,J)*DEPST(J)
 13   CONTINUE
C
C===============================================
C  CALCUL DE L INCREMENT DE DEFORMATIONS TOTALES
C===============================================
C
      CALL ELAST1(2,IFOUR,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,XNUV,
     1          XCARB,ICARA,MFR1,NSTRSS,DD,DDV,KERRE,INDIC,ITHHER)
      CALL ELAST1(2,IFOUR,VAR0,NNVARI,XMAT2,NXMATT,YUNGV,XNUV,
     1          XCARB,ICARA,MFR1,NSTRSS,DDINV,DDV,KERRE,INDIC,ITHHER)
           CALL ZDANUL(DEPS,NSTRSS)
           AA=1.D0
           DO 45 I=1,NSTRSS
             DEPS(I)=DEPST(I)
             IF (I.GT.3) AA=0.D0
             DO 46 J=1,NSTRSS
               DEPS(I)=DEPS(I)+(DDINV(I,J)*SIG0(J))
               DEPS(I)=DEPS(I)-(DD(I,J)*SIG0(J))
46           CONTINUE
             DEPS(I)=DEPS(I)-(AA*CTEPS)
45         CONTINUE
C
C ------
C
      CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,NYN,
     1     YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T,
     2     XMAT,NXMATT,TRUC,NCOURB)
 
    5 DO 10 I=1,NSTRSS
                      SIGF(I)=SIG0(I)+DSIGT(I)
                      SIGMA(I)=SIG0(I)
   10 CONTINUE
      DO 20 I=1,NVARI
                      VARF(I)=VAR0(I)
   20 CONTINUE
C
      IF (ITHHER.EQ.2.AND.IALPH.GT.0) GOTO 32
   30 CALL CRIDAM(SIGF,NSTRSS,VARF,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1NYKX,YKX,NNKX,NKX,NNVARI,MFR1,XF,KERRE,SIGGD,YY,T,TRUC,
     2                                                    NNFUS)
      IF (XF.LT.PRECI) GOTO 1000
      CALL CRIDOM(SIG0,NSTRSS,VAR0,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1NNVARI,MFR1,XG,KERRE,SIGGD,Y)
C
C............ Y DESIGNE LE SEUIL D'ENDOMMAGEMENT .................
      IF (Y.EQ.ZEERO) THEN
                           INDIC=1
                           GOTO 32
      ENDIF
      TTEST=XG/Y
      IF (TTEST.GT.PRECIG) INDIC=1
C
C
C ==============
C ON A PLASTIFIE
C ==============
C ------------------------------------------------------------------
C DEFINITION DE DD, BB, R, S, XMU, ET QSI:
C  .          .         .
C SIG = DD*EPSELAS + BB*Q   AVEC:
C         - SIG     = CONTR.
C         - EPSELAS = DEFORM. ELASTIQUES
C         - DD      = MATR. DE HOOK ENDOMMAGEE
C         - Q       = VARIABLES INTERNES PILOTANT LES LOIS DU MODELE
C R   = LOI D'ECOULEMENT PLASTIQUE
C S   = LOIS D'EVOLUTION DES VARIABLES Q
C XMU = DERIVEE DU CRITERE PAR RAPPORT AUX CONTR.
C QSI = DERIVEE DU CRITERE PAR RAPPORT AUX VARIABLES Q
C ------------------------------------------------------------------
   32 ICO=0
C
C
C
C ------------------------------------------
C DEBUT DE LA BOUCLE DES ITERATIONS INTERNES
C --------------------------------------------------------------------
C A CHAQUE ITERATION, CORRESPONDANT A SIGF ET VARF CALCULES AU RANG I,
C "ENDOM1" CALCULE DD, BB, R, S, XMU ET QSI EN FONCTION DU MODELE DU
C MATERIAU ENDOMMAGEABLE.
C --------------------------------------------------------------------
C
   35 CALL TRACTT(YKX,NKX,TRUC,NYKX,NNKX,NNFUS,T)
      NFUST=2*NNFUS
      CALL DERTRA(NYOG,YOG,T,TYUNG,YUNGV,TINF,TSUP)
      CALL DERTRA(NYNU,YNU,T,TNU,XNUV,TINF,TSUP)
      CALL ENDOM1(SIGF,NSTRSS,VARF,XMAT,NXMATT,XCARB,ICARA,
     1NNFUS,TRUC,NNVARI,MFR1,IFOURB,DD,DDINV,
     2BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3YUNGV,XNUV)
      IF (KERRE.NE.0) GOTO 1000
C
      CALL ZDANUL(DDR,NSTRSS)
      DO 50 J=1,NSTRSS
      DO 40 K=1,NSTRSS
                      DDR(J)=DDR(J)+DD(J,K)*R(K)
   40 CONTINUE
   50 CONTINUE
C
      CALL ZDANUL(BBS,NSTRSS)
      DO 70 J=1,NSTRSS
      DO 60 K=1,NNVARI
                      BBS(J)=BBS(J)+BB(J,K)*S(K)
   60 CONTINUE
   70 CONTINUE
C
      DLAMA1=0.D0
      DLAMA2=0.D0
      DO 80 J=1,NSTRSS
                      DLAMA1=DLAMA1+XMU(J)*(DDR(J)-BBS(J))
   80 CONTINUE
      DO 90 J=1,NNVARI
                      DLAMA2=DLAMA2+QSI(J)*S(J)
   90 CONTINUE
C
C -------------------------------------------------------------
C CALCUL DE LA VARIATION DU MULT. PLAST. PAR LE CRIT. DE PLAST.
C -------------------------------------------------------------
      DT=T-T0
      DLAMDA=XF/(DLAMA1-DLAMA2)
C
C ------------------------------------------------------------------
C CALCUL DES CONTR., DES VAR. INT. Q ET DE LA DEFORM. PLAST. CUMULEE
C AU RANG I+1
C ------------------------------------------------------------------
      DO 100 J=1,NSTRSS
         SIGF(J)=SIGF(J)-DLAMDA*(DDR(J)-BBS(J))
  100 CONTINUE
      VARF(2)=VARF(2)+DLAMDA*S(1)
C
      IF (INDIC.EQ.1) THEN
      S2=S(2)
      IF (IFOURB.NE.-2) THEN
      XA=1.D0+VARF(3)
      CALL ENDOM1(SIGF,NSTRSS,VARF,XMAT,NXMATT,XCARB,ICARA,
     1NNFUS,TRUC,NNVARI,MFR1,IFOURB,DD,DDINV,
     2BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3YUNGV,XNUV)
 
      IF (KERRE.NE.0) GOTO 1000
      XB=VARF(3)+DLAMDA*S(2)*(1.D0-VARF(3))
      XDELTA=(XA*XA)-4.D0*XB
            IF (XDELTA.LT.0.D0) THEN
               VARF(3)=VARF(3)+DLAMDA*S2
            ELSE
                XDELTA=SQRT(XDELTA)
                VARF(3)=(XA-XDELTA)/2.D0
            ENDIF
      ELSE
C
C  CAS DES CONTR. PLANES - MAJ DE D
C
      XA=VARF(3)+0.5D0*DLAMDA*S(2)
      CALL ENDOM1(SIGF,NSTRSS,VARF,XMAT,NXMATT,XCARB,ICARA,
     1NNFUS,TRUC,NNVARI,MFR1,IFOURB,DD,DDINV,
     2BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3YUNGV,XNUV)
      XB=0.5D0*DLAMDA*S(2)*(1.D0-VARF(3))
      XDELTA=(1.D0-XA)*(1.D0-XA)-4.D0*XB
            IF (XDELTA.LT.0.D0) THEN
               VARF(3)=VARF(3)+DLAMDA*S2
            ELSE
               XDELTA=SQRT(XDELTA)
               VARF(3)=(1.D0+XA-XDELTA)/2.D0
            ENDIF
      ENDIF
      ENDIF
C
      VARF(1)=VARF(1)+DLAMDA/(1.D0-VARF(3))
C
C CONTR. PLANES - MAJ DES DEFORM. PLAST.
C
      IF (IFOURB.EQ.-2) THEN
        CALL ENDOM1(SIGF,NSTRSS,VARF,XMAT,NXMATT,XCARB,ICARA,
     1  NNFUS,TRUC,NNVARI,MFR1,IFOURB,DD,DDINV,
     2  BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3  YUNGV,XNUV)
        DO 110 J=1,NSTRSS
              VARF(1+NNVARI+J)=VARF(1+NNVARI+J)+DLAMDA*R(J)
  110 CONTINUE
      ENDIF
C
C
C
C -------------------------------------------------------------
C SI LE CRIT. PLAST. EST ATTEINT, ON ARRETE LES ITERATIONS
C SI LE NBR. MAXIMUM D'ITERATIONS EST DEPASSE, IDEM MAIS ERREUR
C -------------------------------------------------------------
  115 CALL CRIDAM(SIGF,NSTRSS,VARF,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1NYKX,YKX,NNKX,NKX,NNVARI,MFR1,XF,KERRE,SIGGD,YY,T,TRUC,
     2                                                   NNFUS)
      TEST=ABS(XF/YY)
      IF (KERRE.NE.0) GOTO 1000
      IF (TEST.LT.PRECI.AND.INDIC.EQ.1) GOTO 120
      IF (TEST.LT.PRECI.AND.INDIC.EQ.0) THEN
      CALL CRIDOM(SIGF,NSTRSS,VARF,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1NNVARI,MFR1,XG,KERRE,SIGGD,Y)
             TTEST=XG/Y
             ATTEST=ABS(TTEST)
             IF (ATTEST.LT.PRECG) GOTO 120
             IF (TTEST.LE.PRECIG) THEN
                IF (ALPH.EQ.1.D0) GOTO 120
                A=1.D0
             ELSE
                A=-1.D0
             ENDIF
             ALPH=ALPH1+0.5D0*A*ABS(ALPH1-ALPH0)
             ALPH0=ALPH1
             ALPH1=ALPH
             IALPH=IALPH+1
             T=(1-ALPH)*T0+ALPH*T1
             X=VARF(2)
             CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,
     1       NYN,YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T,
     2       XMAT,NXMATT,TRUC,NCOURB)
             ALFFA=XMAT(4)
C
             CALL ELAST1(1,IFOURB,VAR0,NNVARI,XMAT,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,MFR1,NSTRSS,DDINV,DDV,KERRE,INDIC,ITHHER)
             CALL ELAST1(2,IFOURB,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,MFR1,NSTRSS,DD,DDV,KERRE,INDIC,ITHHER)
             CALL MULMAT(DDV,DDINV,DD,NSTRSS,NSTRSS,NSTRSS)
C
             CTEPS=ALFA0*(T0-TREF)-ALFFA*(T-TREF)
             CALL ZDANUL(DSIGT,NSTRSS)
C
             DO 8 I=1,NSTRSS
             AA=1.D0
               DO 9 J=1,NSTRSS
                IF (J.GT.3) AA=0.D0
                DSIGT(I)=DSIGT(I)+DDINV(I,J)*((ALPH*DEPS(J))+(AA*CTEPS))
                DSIGT(I)=DSIGT(I)+DDV(I,J)*SIG0(J)
9              CONTINUE
               DSIGT(I)=DSIGT(I)-SIG0(I)
8            CONTINUE
             IF (IALPH.GT.50) GOTO 1000
             GOTO 5
      ENDIF
      ICO=ICO+1
      IF (ICO.GT.ITMAX) THEN
                            KERRE=2
                            GOTO 1000
      ENDIF
      GOTO 35
C ----------------------------------------
C FIN DE LA BOUCLE DES ITERATIONS INTERNES
C ----------------------------------------
C
C
C ====================================================
C FIN DES ITERATIONS INTERNES - LE CRITERE EST VERIFIE
C                 MAJ DES VARIABLES
C ====================================================
C
C ----------------------------------------------------------------
C CALCUL DES DEFORM. ELAST. A LA FIN DU PAS - ELLES SERONT DS. DDR
C ----------------------------------------------------------------
  120 CALL ELAST1(2,IFOURB,VARF,NNVARI,XMAT,NXMATT,YUNGV,XNUV,
     1XCARB,ICARA,MFR1,NSTRSS,DD,DDV,KERRE,INDIC,ITHHER)
      CALL ZDANUL(DDR,NSTRSS)
      DO 140 J=1,NSTRSS
      DO 130 K=1,NSTRSS
                        DDR(J)=DDR(J)+DD(J,K)*SIGF(K)
  130 CONTINUE
  140 CONTINUE
C
C -----------------------------------------
C CALCUL DES DEFORM. PLAST. A LA FIN DU PAS
C ----------------------------------------------------------------
C ON CALCULE D'ABORD, DS. BBS, LES DEFORM. TOTALES A LA FIN DU PAS
C ----------------------------------------------------------------
      CALL ELAST1(2,IFOURB,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,XNUV,
     1XCARB,ICARA,MFR1,NSTRSS,DD,DDV,KERRE,INDIC,ITHHER)
      CALL ZDANUL(BBS,NSTRSS)
      DO 170 J=1,NSTRSS
      DO 160 K=1,NSTRSS
                      BBS(J)=BBS(J)+DD(J,K)*SIGMA(K)
  160 CONTINUE
  170 CONTINUE
      AA=1.D0
      DO 180 J=1,NSTRSS
      IF (J.GT.3) AA=0.D0
      BBS(J)=BBS(J)+VAR0(1+NNVARI+J)+(ALPH*DEPS(J))+(AA*CTEPS)
  180 CONTINUE
C
C
      DO 190 J=1,NSTRSS
               IF (IFOURB.EQ.-2) GOTO 195
               VARF(1+NNVARI+J)=BBS(J)-DDR(J)
  195          DEFP(J)=VARF(1+NNVARI+J)-VAR0(1+NNVARI+J)
  190 CONTINUE
      IF (ALPH.EQ.1.D0) GOTO 900
      IF (INDIC.EQ.0) THEN
                          INDIC=1
                          ALPH=1.D0-ALPH
                         DO 520 J=1,NVARI
                          VAR0(J)=VARF(J)
  520                    CONTINUE
                         DO 525 J=1,NSTRSS
                          SIGMA(J)=SIGF(J)
  525                    CONTINUE
C
                          T0=T
                          T=T1
                          DO 530 J=1,NXMATT
                          XMAT1(J)=XMAT(J)
                          XMAT(J)=XMAT2(J)
  530                     CONTINUE
                          ALFA0=XMAT1(4)
                          ALFFA=XMAT(4)
             CALL ELAST1(1,IFOURB,VAR0,NNVARI,XMAT,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,MFR1,NSTRSS,DD,DDV,KERRE,INDIC,ITHHER)
             CALL ELAST1(2,IFOURB,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,MFR1,NSTRSS,DDINV,DDV,KERRE,INDIC,
     2       ITHHER)
             CALL MULMAT(DDV,DD,DDINV,NSTRSS,NSTRSS,NSTRSS)
C
             CTEPS=ALFA0*(T0-TREF)-ALFFA*(T-TREF)
             CALL ZDANUL(DSIGT,NSTRSS)
C
             DO 515 I=1,NSTRSS
               AA=1.D0
               DO 510 J=1,NSTRSS
                 IF (J.GT.3) AA=0.D0
                 DSIGT(I)=DSIGT(I)+DD(I,J)*((ALPH*DEPS(J))+(AA*CTEPS))
                 DSIGT(I)=DSIGT(I)+DDV(I,J)*SIGMA(J)
510              CONTINUE
               DSIGT(I)=DSIGT(I)-SIGMA(I)
               SIGF(I)=SIGF(I)+DSIGT(I)
515            CONTINUE
                          GOTO 30
      ENDIF
C
C ==================================================================
C SI LE PT. DE GAUSS EST ROMPU, ON ANNULE PROGRESSIVEMENT SES CONTR.
C ==================================================================
  900 DCT=XMAT(7)
      DCC=1.1D0*DCT
      XD=VARF(3)
      IF (XD.GE.DCC) THEN
                              VARF(3)=1.D0
                              CALL ZDANUL(SIGF,NSTRSS)
      ENDIF
      IF (XD.GT.DCT.AND.XD.LT.DCC) THEN
                              VARF(3)=1.D0
                              BETA=11.D0-10.D0*(XD/DCT)
                              DO 950 J=1,NSTRSS
  950                               SIGF(J)=SIGF(J)*BETA
 
      ENDIF
 1000 RETURN
      END
*******************************************************************