Sources Esope | Cast3M

umapli
C UMAPLI    SOURCE    CHAT      05/01/13    03:58:05     5004
      SUBROUTINE UMAPLI ( STRESS, STATEV, DDSDDE,
     &                  STRAN, DSTRAN, DTIME,
     &                  TEMP, DTEMP, PREDEF, DPRED,
     &                  NTENS, NSTATV,
     &                  PROPS, NPROPS, COORDS,
     &                  NOEL, NPT, KINC )
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C
C     DESCRIPTION FONCTIONNELLE :
C     -------------------------
C
C     Dependance du module Utilisateur UMAT
C     Appelee dans le cas suivant :
C
C     ILOI =  5 : Modele 'PLASTIQUE' 'ISOTROPE' de CAST3M, INPLAS=5
C                 RESTRICTION A LA FORMULATION MASSIVE
C
C     Composantes de materiau
C     -----------------------
C     PROPS ( 1) : 'YOUN'
C     PROPS ( 2) : 'NU  '
C
C     Courbe de traction : pour 1 =&lt; I =&lt; LTRAC
C
C     PROPS ( 2+2*I-1) : 'Eiii'   iii formatage I3 de la valeur de I
C     PROPS ( 2+2*I  ) : 'Siii'
C
C     Variables internes
C     ------------------
C     STATEV( 1) : 'EPSE'
C
C     STATEV( 2) : 'ENXX'
C     STATEV( 3) : 'ENYY'
C     STATEV( 4) : 'ENZZ'
C     STATEV( 5) : 'GNXY'
C     STATEV( 6) : 'GNXZ'
C     STATEV( 7) : 'GNYZ'
C
C-----------------------------------------------------------------------
C     Arguments de l'interface
C
      INTEGER       NTENS, NSTATV, NPROPS,
     &              NOEL, NPT, KINC
C
      REAL*8        STRESS(NTENS), STATEV(*),
     &              DDSDDE(NTENS,NTENS),
     &              STRAN(NTENS), DSTRAN(NTENS),
     &              DTIME,
     &              TEMP, DTEMP, PREDEF(*), DPRED(*),
     &              PROPS(NPROPS),
     &              COORDS(3)
C
C     Variables locales
C
C     Tableaux de travail
C
      REAL*8     SIG(130), EPS(130)
      REAL*8     S(8), SX(8), DS(8), DSIG(8), SPHER1(8), DSPHER1(8),
     &           DDEPSE(8), F(8), W1(8), W2(8), SIGB(8), Z1(8)
C
      REAL*8     TRAC(260), DSIGT(8), WORK(1)
C
C     Variables scalaires
C
      INTEGER    MFR1, IFOURB, I, J, LTRAC, IBI, NITER, ITRY
C
      REAL*8     YUNG, XNU, ALFAH, EP1, EPST, AUX0, AUX, AUX1, GEGE,
     &           SI, SEQ, PHI, PETI, PHI0, SI0, RR,
     &           COEF, PENTE, RP, C, DENOM, DELTA, DMU, DP, DR,
     &           APHI, APHI0, TEST
C
      REAL*8     UN, DEUX, UNDEMI
      PARAMETER (UN=1.D0,DEUX=2.D0,UNDEMI=.5D0)
C
      REAL*8     PRECIS
      DATA       PRECIS/1.e-8/
C
      REAL*8     VONMIS0
      EXTERNAL   VONMIS0
C
C------------------- Debut du code executable --------------------------
C
C=======================================================================
C 1 - INITIALISATIONS
C=======================================================================
C
C     Indicateurs de modelisation pour les dependances issues de CAST3M
C     + constantes geometriques (initialisations arbitraires)
C
      MFR1   = 1
      IFOURB = 2
C
      ALFAH=1.D0
      EP1=0.D0
C
C     Tableaux de travail
C
      DO 100 I=1,NTENS
         S(I)=0.D0
         SPHER1(I)=0.D0
 100  CONTINUE
C
C     Recuperation des constantes de materiau
C
      YUNG=PROPS(1)
      XNU=PROPS(2)
C
C     Recuperation de la courbe de traction
C
      CALL CRBTRA(PROPS,NPROPS,TRAC,LTRAC,KINC)
      IF (KINC.LT.0) RETURN
      DO 101 I=1,LTRAC
         SIG(I)=TRAC(2*I-1)
         EPS(I)=TRAC(2*I)
 101  CONTINUE
C
C     Recuperation des variables internes (deformations plastiques)
C
      EPST=STATEV(1)
C
C=======================================================================
C 2 - Prediction elastique des contraintes
C     Calcul de l'increment de contraintes conforme a CALSIG
C     Calcul des coefficients de la matrice de Hooke conforme a HOOKIS
C=======================================================================
C
      AUX0 = UN/((UN+XNU)*(UN-DEUX*XNU))
      AUX  = YUNG*AUX0*(UN-XNU)
      AUX1 = YUNG*AUX0*XNU
      GEGE = YUNG*UNDEMI/(UN+XNU)
C
      DO 200 J=1,NTENS
         DO 201 I=1,NTENS
            DDSDDE(I,J)=0.0D0
 201     CONTINUE
 200  CONTINUE
C
      DDSDDE(1,1) = AUX
      DDSDDE(1,2) = AUX1
      DDSDDE(1,3) = AUX1
C
      DDSDDE(2,1) = AUX1
      DDSDDE(2,2) = AUX
      DDSDDE(2,3) = AUX1
C
      DDSDDE(3,1) = AUX1
      DDSDDE(3,2) = AUX1
      DDSDDE(3,3) = AUX
C
      DDSDDE(4,4) = GEGE
      DDSDDE(5,5) = GEGE
      DDSDDE(6,6) = GEGE
C
      DO 210 I=1,NTENS
         DSIGT(I)=0.0D0
 210  CONTINUE
      DO 211 J=1,NTENS
         DO 212 I=1,NTENS
            DSIGT(I)=DSIGT(I)+DDSDDE(I,J)*DSTRAN(J)
 212     CONTINUE
 211  CONTINUE
C
      DO 220 I=1,NTENS
         S(I)=STRESS(I)+DSIGT(I)
         SIGB(I)=S(I)
         SX(I)=S(I)-SPHER1(I)
 220  CONTINUE
C
C=======================================================================
C 3 - Calcul de la limite elastique SI
C=======================================================================
C
      CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI)
      IF (IBI.EQ.1) THEN
         KINC = -511
         RETURN
      ENDIF
C
C=======================================================================
C 4 - Calcul de la contrainte equivalente SEQ
C=======================================================================
C
      SEQ=VONMIS0(SX,NTENS,MFR1,IFOURB,EP1,ALFAH)
C
C=======================================================================
C 5 - Verification du critere de plasticite
C     Demarrage des iterations internes le cas echeant
C=======================================================================
C
      PHI=SEQ-SI
      NITER=0
      PETI=1.1D0*0.5D0*PRECIS*SEQ
      CALL EPSPRE(SEQ,SI,PETI,ITRY)
C
C     SI ON NE PLASTIFIE PAS avec la prediction elastique => ON SORT
C     L'etat de contraintes est determine par la prediction elastique
C     Les variables internes (deformations plastiques) sont inchangees
C
      IF ( (ITRY.EQ.1).OR.(SEQ.LE.SI) ) THEN
         DO 500 I=1,NTENS
            STRESS(I)=S(I)
 500     CONTINUE
         RETURN
      ENDIF
C
C     SINON : ON PLASTIFIE DONC DEMARRAGE DES ITERATIONS INTERNES
C
C     Initialisations complementaires
C
      PHI0=PHI
      SI0=SI
      RR=0.D0
C
      DO 510 I=1,NTENS
         DSIG(I)=0.D0
         DSPHER1(I)=0.D0
 510  CONTINUE
C
C        _______________________________________________________________
C        Demarrage des iterations internes (bloc REPETER)
C
  10     NITER=NITER+1
C
C        Calcul de W1=DF/D(SIGMA)
C
         F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
         F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
         F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
         DO 1010 I=4,NTENS
            F(I)=SX(I)
1010     CONTINUE
         DO 1011 I=1,3
            W1(I)=1.5D0*F(I)/SEQ
            Z1(I)=W1(I)
1011     CONTINUE
         DO 1012 I=4,NTENS
            W1(I)=3.D0*F(I)/SEQ
            Z1(I)=1.5D0*F(I)/SEQ
1012     CONTINUE
C
         DO 1020 I=1,NTENS
            W2(I)=0.D0
            DO 1021 J=1,NTENS
               W2(I)=W2(I)+DDSDDE(I,J)*W1(J)
1021        CONTINUE
1020     CONTINUE
C
         COEF=0.D0
         DO 1030 I=1,NTENS
            COEF=COEF+W1(I)*W2(I)
1030     CONTINUE
C
         CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC,1,IBI)
         IF (IBI.EQ.1) THEN
            KINC = -512
            RETURN
         ENDIF
C
         RP=PENTE
         C=0.D0
         DENOM=COEF+C+RP
         DELTA=PHI/DENOM
         DMU=C*DELTA/SEQ
C
         DO 1040 I=1,NTENS
            DSIG(I)=-DELTA*W2(I)
            DSPHER1(I)=DMU*SX(I)
1040     CONTINUE
         DP=DELTA
         DR=RP*DP
C
         RR=RR+DR
         EPST=EPST+DP
         DO 1050 I=1,NTENS
            S(I)=S(I)+DSIG(I)
            SPHER1(I)=SPHER1(I)+DSPHER1(I)
            SX(I)=S(I)-SPHER1(I)
1050     CONTINUE
         SEQ=VONMIS0(SX,NTENS,MFR1,IFOURB,EP1,ALFAH)
C
C        Test de convergence
C        Calcul de la nouvelle valeur de PHI
C
         CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI)
         PHI=SEQ-SI
C
         PETI=1.E-7
         APHI=ABS(PHI)
         APHI0=ABS(PHI0)
         TEST=PETI*APHI0
C
         IF (NITER.GT.50) THEN
            KINC = -513
            RETURN
         ENDIF
C
C        SI ON A CONVERGE : on calcule l'etat final, on stocke les
C        resultats et on sort
C
         IF (APHI.LE.TEST) THEN
C
            DO 1060 I=1,NTENS
               DS(I)=S(I)-SIGB(I)
1060        CONTINUE
            CALL EPSIG0(DS,DDEPSE,MFR1,IFOURB,YUNG,XNU,WORK,NTENS)
            STATEV(1)=EPST
            DO 1070 I=1,NTENS
               STRESS(I)=S(I)
               STATEV(1+I)=STATEV(1+I)-DDEPSE(I)
1070        CONTINUE
C
C        SINON : on itere
C
         ELSE
            GOTO 10
         ENDIF
C
C        _______________________________________________________________
C        Fin des iterations internes
C
      RETURN
      END