C POUDRA    SOURCE    PV        09/12/22    21:15:02     6589
C
C  PLGDebut : implantation des equations relatives au modele
C  viscoplastique poudre.
C
C  Modification ChD division par des rÅels et non pas des entiers
C  COEC doit etre plus grand que 1.0 et non pas 0.0
C
C  Modification ChD nouvelles equations de c et f selon
C  les equations de DA pour sa these.
C
C  Modification LS nouvelle loi viscoplastique a trois parametres
C  Lois viscoplastiques d Abouaf et 3S sont disponibles ...
C
       SUBROUTINE POUDRA(SIG,EPSV,VAR,EPSVPT,VARPT,XMAT,NSTRS,
     &                  NVARI,NCOMAT,KERRE)
C
      IMPLICIT INTEGER(I-N)
       IMPLICIT REAL*8(A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
       DIMENSION SIG(*),EPSV(*),VAR(*),EPSVPT(*),VARPT(*)
       DIMENSION XMAT(*)
       DIMENSION XX(6)
C
C---------------------------------------------------------------------
C      ====> Connaissant l'etat, c'est a dire SIG, EPSV et VAR on
C      calcule les derivees EPSVPT et VARPT(1) en utilisant les
C      constantes XMAT
C----- calcul des invariants S1=Tr(sigma)
C                             S2C=S2*S2=Tr(Dev(sigma):Dev(sigma))
C                        avec Dev(sigma)=sigma-1/3Tr(sigma)
C
       S1 = SIG(1)+SIG(2)+SIG(3)
C
C      Modif Ch.D du 20/10/94 : Cf DEA Canel page A18
C      assure tr(XX) = 0
       DO I=1,2
          XX(I) = SIG(I)-S1/3.0D0
       END DO
       XX(3) = -XX(1)-XX(2)
C      Fin de modification Ch.D du 20/10/94
C
       DO I=4,NSTRS,1
          XX(I) = SIG(I)
       END DO
       S2C = PROCON(XX,XX,NSTRS)
C
C----- calcul de F et C : fonctions de rho
       RHOR = VAR(NVARI)
       A    = XMAT(6)
       RN   = XMAT(5)
       QSRT = XMAT(7)
       F0 = XMAT(8)
       F1 = XMAT(9)
       F2 = XMAT(10)
       F3 = XMAT(11)
       F4 = XMAT(12)
       F5 = XMAT(13)
       C0 = XMAT(14)
       C1 = XMAT(15)
       C2 = XMAT(16)
       C3 = XMAT(17)
       C4 = XMAT(18)
       C5 = XMAT(19)
       RMIN = F5+1.0D-03
       IF (RHOR.LT.RMIN) THEN
          RHOR1 = RMIN
       ELSE
          RHOR1 = RHOR
       ENDIF
C
C----- Modification L3S du PROGRAMME :
C      Integration de la loi a 3 parametres
C      Le 15/10/1997 par Laurent SANCHEZ ...
C      COEF   : Coefficient F(Rhor)
C      COEC   : Coefficient C(Rhor)
C      COEE   : Coefficient E(Rhor)
C      COSS   : Coefficient S(Rhor)
C
C      COF1   : Coefficient  f1(Rhor)
C      COF2   : Coefficient  f2(Rhor)
C      COF3   : Coefficient  f3(Rhor)
C
       IF ((1.0D0-RHOR1).LT.1.0D-06) THEN
          COEF = 0.0D0
          COEC = 1.0D0
       ELSE
          COEF = F0*(((1.D0-RHOR1)/(RHOR1-F5))**(F1*RHOR1+F2))
          COEC = 1.0D0+C0*(((1.D0-RHOR1)/(RHOR1-F5))**(C1*RHOR1+C2))
       ENDIF
       IF (COEF.GT.(C5)) COEF = C5
       IF (COEC.GT.(C5)) COEC = C5
C
C----- calcul de sigma equivalent :
       SIGEQ = SQRT(COEF*S1*S1+1.5D0*COEC*S2C)
C
C----- calcul de EPSVPT en fonction du parametre A :
C      A                  : INT(F4) = 0
C      Log10(A)           : INT(F4) = 1
C      loi a 3 parametres : INT(F4) = 2
C
*pv pour être sur d'avoir la bonne valeur, on augmente un peu f4
       NF4=INT(F4+0.1)
C
C----- Loi d'Abouaf :
       IF (NF4.EQ.1) THEN
          DO I=1,NSTRS,1
             B = 0.D0
             IF (I.LE.3) B=1.D0
             EPSVPT(I) = 1.D0*(10.D0**A)*(SIGEQ**(RN-1.D0))*RHOR
     &                   *(COEF*S1*B+1.5D0*COEC*XX(I))
          END DO
C
C----- Loi viscoplastique 3S a trois parametres :
       ELSE IF (NF4.EQ.2) THEN
          IF ((1.D0-RHOR1).LT.1.D-6) THEN
             COEE = 0.0D0
             COES = 1.0D0
             COF1 = 0.D0
             COF2 = 1.D0
             COF3 = 0.D0
          ELSE
             COEE = (3.D0*COEF)/(COEF+COEC)
             COSS = (COEF+COEC)**(-(RN+1.D0)/(2.D0*RN))
             COF1 = (2.D0*COEF/3.D0)**((RN+1.D0)/2.D0)
             COF2 = (1.D0-COEE)*(COSS**(-RN))
     &              + 2D0*( (1.5D0)**((RN+1.D0)/2.D0) )*COF1
             COF3 = COEE*(COSS**(-RN))*((2.D0/3.D0)**((RN+5.D0)/4.D0))
     &              - 2.D0*( (1.5D0)**((RN+1.D0)/4.D0) )*COF1
          ENDIF
          DO I=1,NSTRS,1
             B = 0.D0
             IF (I.LE.3) B=1.D0
             EPSVPT(I) = (3.D0/2.D0)**((RN+1.D0)/2.D0)*EXP(-QSRT)
     &          *(10.D0**A)*RHOR *((COF1*S1*(S1**2.D0)**((RN-1.D0)/2.D0)
     &          +0.5D0*COF3*S1*S2C*(S1**2.D0*S2C)**((RN-3.D0)/4.D0))
     &          *B+(COF2*S2C**(RN/2.D0)+0.5D0*COF3*S1**2.D0*S2C**0.5D0
     &          *(S1**2.D0*S2C)**((RN-3.D0)/4.D0))*XX(I))
          END DO
C
C----- Loi viscoplastique par  defaut
       ELSE
          DO I=1,NSTRS,1
             B = 0.D0
             IF (I.LE.3) B=1.D0
             EPSVPT(I) = 1.D0*A*(SIGEQ**(RN-1.D0))*EXP(-QSRT)*RHOR
     &                   *(COEF*S1*B+1.5D0*COEC*XX(I))
          END DO
       ENDIF
C
C----- calcul de VARPT(1) derivee de EPSE
       PPT = SQRT(2.D0/3.D0*PROCON(EPSVPT,EPSVPT,NSTRS))
       VARPT(1) = PPT
C---------------------------------------------------------------------
C
       RETURN
       END