Sources Esope | Cast3M

cmazar
C CMAZAR    SOURCE    OF166741  25/11/04    21:15:26     12349          
      SUBROUTINE CMAZAR (WRK52,WRK53,WRK54,WRKK2,NSTRS1,NVARI,
     1                   ICARA,JDIM,IFOUR2)
C MAZARS    SOURCE    AM        98/12/23    21:38:30     3409
C
C
C     variables en entree
C
C
C     WRK0     pointeur sur un segment deformation au pas precedent
C
C     WRK1     pointeur sur un segment increment de deformation
C
C     WRKK2     pointeur sur un segment variables internes au pas precedent
C
C     WRK5      pointeur sur un segment de deformations inelastiques
C
C     XMATER     constantes du materiau
C
C     NSTRS1      nombre de composantes dans les vecteurs des contraintes
C                                        et les vecteurs des deformations
C
C     NVARI      nombre de variables internes (doit etre egal a 2)
C
C     NMATT      nombre de constantes du materiau
C
C     ISTEP      flag utilise pour separer les etapes dans un calcul non local
C                ISTEP=0 -----> calcul local
C                ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
C                ISTEP=2 -----> calcul non local etape 2 on continue le calcul
C                               a partir des seuils moyennes
C
C Modif L.Bode - 14/10/92
C Nouveaux parametres en entree
C     JDIM       Dimension de travail
C                ( Coques JDIM =2 , Massifs JDIM = IDIM )
C     IFOUR2     Type de formulation
C                ( Coques IFOUR2 = -2 => contraintes planes ,
C                  Massifs IFOUR2 = IFOUR)
C
C     variables en sortie
C
C     VARINF      variables internes finales
C
C     SIGMAF      contraintes finales
C
C     C. LA BORDERIE MARS 1992
C     declaration des variables
C
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC PPARAM
-INC CCOPTIO
-INC DECHE
 
      SEGMENT WRKK2
        REAL*8 EPSILI(NSTRS1)
        REAL*8 EPSILO(NSTRS1)
      ENDSEGMENT
 
      INTEGER NSTRS1,NVARI
      REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3)
      REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC
      REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA
      INTEGER ISTRS,JSTRS,KCAS,IRTD
      REAL*8 XZERO,UN,DEUX,XPETIT
      REAL*8 DINI,D,DT,DC,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA
      PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
      PARAMETER (YPETIT=1.)
C
C     recuperation des variables initiales dans les tableaux
C
      YOUN = XMAT(1)
      XNU  = XMAT(2)
      EPSD0= XMAT(5)
      ACOM = XMAT(6)
      BCOM = XMAT(7)
      ATRA = XMAT(8)
      BTRA = XMAT(9)
      BETA = XMAT(10)
      DINI = VAR0(2)
C
C     calcul des seuils
C
C        calcul de la deformation totale
C
      DO 100 ISTRS=1,NSTRS1
         EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS)
100   CONTINUE
C
C        calcul des deformations principales
C
C        on reecrit les deformations sous forme matricielle
C
C Modif L.Bode - 14/10/92
C Rajout de IFOUR2 en argument de ENDOCA
*      print*,'on appelle ENDOCB'
      CALL ENDOCB (EPSILO,EPS33,2,IFOUR2)
*      print*,'apres endocb'
C Fin modif L.Bode
C
C        et on diagonalise
C
C Modif L.Bode - 14/10/92
C Pour les elts Coques, on travaille en contraintes planes => JDIM =2
C Pour les elts Massifs JDIM =IDIM
*      print*,'avant JACOB3'
      CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33)
*      print*,'apres JACOB3'
C Fin modif L.Bode
      IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN
C
C       on calcule la matrice de hooke et les contraintes ppales
C
         CMATE = 'ISOTROPE'
         KCAS=1
C Modif L.Bode - 14/10/92
C IFOUR --> IFOUR2 dans appel DOHMAS
*         print*,'avant dohmas'
         CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS1,KCAS,DDHOOK,IRTD)
*         print*,'apres dohmas'
C Fin modif L.Bode
         DO 200 ISTRS=1,3
            r_z = XZERO
            DO 210 JSTRS=1,3
              r_z = r_z + DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
210         CONTINUE
            SIGP(ISTRS)= r_z
200      CONTINUE
      END IF
C
C        on complete la deformation dans le cas des contraintes planes
C
C Modif L.Bode - 14/10/92
C IFOUR remplace par IFOUR2
      IF (IFOUR2.EQ. -2) THEN
         EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU)
      END IF
C Fin modif L.Bode
C
C        on calcule le epsilontild
C
      EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
     1       MAX( XZERO , EPSIPP(2) )**2 +
     2       MAX( XZERO , EPSIPP(3) )**2
      EPSTIL=SQRT (EPSTIL)
      epstil=max(xpetit,epstil)
      IF (ISTEP .EQ. 0) THEN
         EPSTIM=EPSTIL
         VARF(1)=EPSTIL
      ELSE IF (ISTEP .EQ. 1) THEN
         VARF(2)=DINI
         VARF(1)=EPSTIL
      ELSE IF (ISTEP .EQ. 2) THEN
         EPSTIM=VAR0(1)
         VARF(1)=EPSTIM
      ELSE
         PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP
      END IF
      IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN
C
C      on calcule l'endommagement et les contraintes
C
C        on calcule ALFAT ALFAC DT et DC puis D
C        dans le cas ou le seuil initial est depasse
C
         IF ( EPSTIM .GT. EPSD0) THEN
C
C      calcul de l'endommagement
C
C        on calcule le signe des contraintes elastiques
C
            DO 300 ISTRS=1,3
               IF (SIGP(ISTRS).LT. XZERO) THEN
                  SIGPC(ISTRS)=SIGP(ISTRS)
                  SIGPT(ISTRS)=XZERO
               ELSE
                  SIGPT(ISTRS)=SIGP(ISTRS)
                  SIGPC(ISTRS)=XZERO
               END IF
300         CONTINUE
            TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
            TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
C
C        on calcule les deformations dues aux contraintes positives
C
            DO 400 ISTRS=1,3
               EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
400         CONTINUE
C
C        on en deduit ALFAT et ALFAC
C
            ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
     1              MAX(XZERO,EPSIPP(2))*EPSILT(2) +
     2              MAX(XZERO,EPSIPP(3))*EPSILT(3)
            ALFAT = ALFAT/(EPSTIL*EPSTIL)
            ALFAC = UN - ALFAT
C
C       modification pour la bi ou tricompression
C
            IF (TRSIGC.LT. -YPETIT .AND. TRSIGT.LT.YPETIT) THEN
               GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
     1                                 SIGPC(3)*SIGPC(3)
               GAMMA=-SQRT(GAMMA)/TRSIGC
               EPSTIM=EPSTIM*GAMMA
            END IF
C
C        amelioration de la reponse en cisaillement pour beta > 1.
C
            IF (BETA .GT. UN) THEN
               IF ( ALFAT .GT. XPETIT ) THEN
                  ALFAT=ALFAT**BETA
               END IF
               IF ( ALFAC .GT. XPETIT ) THEN
                  ALFAC=ALFAC**BETA
               END IF
            END IF
C
C        on calcule DT et DC puis D
C        dans le cas ou le seuil initial est depasse
C
C        on est oblige de verifier car on a pu multiplier par gamma
C
C            IF (EPSTIM .GT. EPSD0) THEN
C               DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
C     1            ATRA*EXP(-BTRA*(EPSTIM-EPSD0))
C               DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
C     1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
C            ELSE
C               DT=XZERO
C               DC=XZERO
C            END IF
CLB 2010/01/20
C            trois lois d'évolution:
C             ATRA > 0 : loi de mazars classique
C           -10&lt;ATRA&lt;0 : loi d'évolution exponentielle modifiée pour le GF
C             ATRA &lt; -10 : loi d'évolution linéaire  --> BTRA est alors
C                  la déformation pour laquelle la contrainte s'annule C
C
            IF (EPSTIM .GT. EPSD0) THEN
               IF (ATRA .GT. 0.D0) THEN
                 DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
     1              ATRA*EXP(BTRA*(EPSD0-EPSTIM))
               ELSE IF (ATRA . GT. -10.) THEN
                 DT=UN - epsd0/epstim*EXP(BTRA*(EPSD0-EPSTIM))
               ELSE
                 IF(EPSTIM .LT. BTRA) THEN
                   DT=UN - EPSD0*(BTRA - EPSTIM)/EPSTIM/(BTRA - EPSD0)
                 ELSE
                   DT=1.D0
                 ENDIF
               END IF
               DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
     1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
            ELSE
 
               DT=XZERO
               DC=XZERO
            END IF
            D = ALFAT*DT + ALFAC*DC
C
C           on borne la valeur de D a 0.99999999 pour avoir une raideur residuelle evitant des points flottants
C
*pv         D=MIN ( D , UN-1.D-20 )
            D=MIN ( D , UN-1.D-8 )
         ELSE
            D=XZERO
         END IF
C
C           on teste la croissance de D
C
         D=MAX ( D , DINI )
C
C           on le stocke dans les variables internes finales
C
         VARF(2)= D
C
C           on calcule les contraintes finales
C
         CALL MATVE1 (DDHOOK,EPSILO,NSTRS1,NSTRS1,SIGF,2)
         DO 500 ISTRS=1,NSTRS1
            SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
500      CONTINUE
C
C         et les deformations inelastiques finales
C
         DO 600 ISTRS=1,NSTRS1
            EPINF(ISTRS)=EPSILO(ISTRS)*D
600      CONTINUE
 
      ENDIF
 
      RETURN
      END