Télécharger amade2.eso

Retour à la liste

amade2
C AMADE2    SOURCE    CHAT      05/01/12    21:21:02     5004
C AMADE2.eso     SOURCE     bald     MAR 28/03/95     00:00:00
      SUBROUTINE AMADE2(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI,
     &     DEPST,IFOURB,XMAT,NMATT,IVAL,SIGF,DEFP,VARF,KERRE)
C
C-----------------------------------------------------------------------
C
C         PLASTICITE MODELE 2D AMADEI-SAEB ELEMENTS JOINTS
C
C  ENTREES
C
C     IB              = NUMERO DE L'ELEMENT
C     IGAU            = NUMERO DU POINT DE GAUSS
C     NSTRS           = NOMBRE DE COMPOSANTES DE CONTRAINTES
C     SIG0(NSTRS)     = CONTRAINTES INITIALES (AU PAS PRECEDENT)
C     EPIN0(NSTRS)    = DEFORMATIONS INITIALES INEL. (AU PAS PRECEDENT)
C     VAR0(NVARI)     = VARIABLES INTERNES DEBUT
C      VAR0(1)        = DEFORMATION PLASTIQUE EQUIVALENTE: EPSE
C      VAR0(2)        = INCR. DEF. PLAST. DUE A L'OUVERTURE SEULE: EPOU
C      VAR0(3)        = ETAT DU JOINT: STAT
C      VAR0(4)        = U EQUIVALENTE: UEQU
C      VAR0(5)        = DEFORMATION TOTALE: EPS1
C      VAR0(6)        = DEFORMATION TOTALE: EPS2
C     NVARI           = NOMBRE DE VARIABLES INTERNES
C     DEPST(NSTRS)    = INCREMENTS DE DEFORMATION TOTALE
C     XMAT(NCOMAT)    = COMPOSANTES DE MATERIAU
C     NMATT           = NOMBRE DE COMPOSANTES DE MATERIAU
C     IVAL(NCOMAT)    = INDICE DES COMPOSANTES DE MATERIAU
C
C  SORTIES
C
C     SIGF(NSTRS)     = CONTRAINTES FINALES
C     DEFP(NSTRS)     = INCREMENTS DE DEFORMATIONS INELASTIQUES FINALES
C     VARF(NVARI)     = VARIABLES INTERNES FINALES
C     KERRE           = INDICE D'ERREUR
C
C
C-----------------------------------------------------------------------
c
c    déclaration des variables
c
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
      REAL*8 I0, K1, K2, KNI
      REAL*8 KNN, KNT, KTN, KTT
c
      DIMENSION SIG0(*), EPIN0(*), VAR0(*), DEPST(*), XMAT(*)
      DIMENSION SIGF(*), DEFP(*), VARF(*)
      DIMENSION RI0E(2,2), RI0P(2,2), DCON(2)
c
c    paramètres pour le calcul
c
      FATT=XPI/180.D0
c
      NN=100
c
      K1=1.5D0
      K2=4.0D0
c
c    tolérances pour définir la condition de plasticité selon sigma
c    et la condition de nullité des contraintes
c
      EPSS=1.D-6
      EPSI=1.D-6
      TOL1=1.D-20
      TOL2=1.D-15
c
c    extraction des caractéristiques du materiau
c
      FIMU=XMAT(5)*FATT
      SGMT=XMAT(6)
      I0=XMAT(7)*FATT
      S0=XMAT(8)
      B0=XMAT(9)
      UR=XMAT(10)
      UP=XMAT(11)
      KNI=XMAT(12)
      FI0=XMAT(13)*FATT
      VM=XMAT(14)
c
c    extraction de la variable EPOU
c
      EPOU=VAR0(2)
c
c    extraction de la variable UEQU
c
      UEQU=VAR0(4)
c
c    extraction des déformations totales
c
      EPS1=VAR0(5)
      EPS2=VAR0(6)
c
c    calcul des petits pas
c
      DU=DEPST(1)/NN
      DV=DEPST(2)/NN
c
c    incréments de déformation plastique finale
c
      DEFP(1)=0.D0
      DEFP(2)=0.D0
c
c    contraintes finales
c
      SIGF(1)=0.D0
      SIGF(2)=0.D0
c
c    variables internes finales
c
      VARF(1)=0.D0
      VARF(2)=0.D0
      VARF(3)=0.D0
      VARF(4)=0.D0
      VARF(5)=0.D0
      VARF(6)=0.D0
c
c **********************************************************************
c ********** test pour verifier la condition initiale du joint *********
c **********************************************************************
c
      IF(EPIN0(2).GT.EPSS)THEN
c
c    le joint est ouvert
c
         IF(DEPST(2).LT.(-EPIN0(2)-EPSI))THEN
c
c       il y a rapprochement des lèvres avec compression (1 cas)
c
            R0=ABS(EPIN0(2)/DEPST(2))
c
            NN=INT(ABS((EPIN0(2)+DEPST(2))/DV))
c
            EPOU=-EPIN0(2)
c
            EPIN0(1)=EPIN0(1)+R0*DEPST(1)
            EPIN0(2)=0.D0
c
            DEFP(1)=R0*DEPST(1)
            DEFP(2)=EPOU
c
            EPS1=EPS1+R0*DEPST(1)
            EPS2=EPS2+R0*DEPST(2)
c
         ELSE IF(DEPST(2).GT.(-EPIN0(2)+EPSS))THEN
c
c       il y a rapprochement des lèvres sans fermeture, ou elles restent
c       bloquées, ou elles s' éloignent (3 cas)
c
            SIGF(1)=0.D0
            SIGF(2)=0.D0
c
            DEFP(1)=DEPST(1)
            DEFP(2)=DEPST(2)
c
            DEDE1=EPIN0(1)+DEPST(1)
            DEDE2=EPIN0(2)+DEPST(2)
c
            VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2))
c
            VARF(2)=DEPST(2)
c
            VARF(3)=1.D0
c
            VARF(4)=VAR0(4)
c
            VARF(5)=EPS1+DEPST(1)
            VARF(6)=EPS2+DEPST(2)
c
            RETURN
c
         ELSE
c
c       il y a fermeture sans compression (1 cas)
c
            SIGF(1)=0.D0
            SIGF(2)=0.D0
c
            DEFP(1)=DEPST(1)
            DEFP(2)=DEPST(2)
c
            DEDE1=EPIN0(1)+DEPST(1)
            DEDE2=EPIN0(2)+DEPST(2)
c
            VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2))
c
            VARF(2)=DEPST(2)
c
            VARF(3)=0.D0
c
            VARF(4)=VAR0(4)
c
            VARF(5)=EPS1+DEPST(1)
            VARF(6)=EPS2+DEPST(2)
c
            RETURN
c
         END IF
c
      ELSE IF((EPIN0(2).LE.EPSS).AND.(SIG0(2).GE.-TOL2))THEN
c
c    le joint est fermé mais sans compression
c
         IF(DEPST(2).GT.EPSS)THEN
c
c       les lèvres s' éloignent (1 cas)
c
            SIGF(1)=0.D0
            SIGF(2)=0.D0
c
            DEFP(1)=DEPST(1)
            DEFP(2)=DEPST(2)
c
            DEDE1=EPIN0(1)+DEPST(1)
            DEDE2=EPIN0(2)+DEPST(2)
c
            VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2))
c
            VARF(2)=DEPST(2)
c
            VARF(3)=1.D0
c
            VARF(4)=VAR0(4)
c
            VARF(5)=EPS1+DEPST(1)
            VARF(6)=EPS2+DEPST(2)
c
            RETURN
c
         ELSE
c
c       les lèvres restent bloquées ou il y a compression (2 cas)
c
         END IF
c
      ELSE
c
c    le joint est fermé avec compression
c
      END IF
c
c **********************************************************************
c ******************* fin du test et debut du calcul *******************
c **********************************************************************
c
      DO 300 KAPPA=1,NN
c
c    definition des deux cas (sigma<SGMT et sigma>=SGMT)
c
         IF(ABS(SIG0(2)).LT.ABS(SGMT))THEN
c
c    cas avec sigma < SGMT
c
c    prédiction élastique
c
            IF(UEQU.LT.UR)THEN
c
c       c' est le cas élastique si U < UR
c
c       definition des quantités necessaires au calcul
c
               VUPU=((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0)
               AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1)
               SR=S0-SIG0(2)*TAN(FI0)
               TP=-SIG0(2)*TAN(FIMU+ATAN(VUPU))*(1.D0-AS)+AS*SR
               TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT)
c
               AA1=(SIG0(2)*(1.D0-AS)*K2)/(SGMT*(COS(FIMU+ATAN(VUPU))**2
     $              ))
               AA2=((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0)
               AA3=1.D0+((1.D0-SIG0(2)/SGMT)**(2.D0*K2))*(TAN(I0)**2)
               AA4=K1*SIG0(2)/SGMT
               AA5=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU+ATAN(VUPU))
     $              *AA4
               AA6=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0))
               AA7=AS*TAN(FI0)
               dTPS1=-(1.D0-AS)*TAN(FIMU+ATAN(VUPU))+AA1*AA2/AA3+AA5+AA6
     $              -AA7
               dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT
               AA8=-UEQU*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0
     $              )
               AA8=AA8+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2)
c
c       calcul de la matrice Dt equivalente (2D)
c
               KNN=1.D0/AA8
               KNT=-((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0)*KNN
               KTN=UEQU/UP*KNN*dTPS1
               KTT=UEQU/UP*KNT*dTPS1+TP/UP
c
            ELSE
c
c       c' est le cas élastique si U >= UR
c
c       definition des quantités necessaires au calcul
c
               VUPU=0.D0
               AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1)
               SR=S0-SIG0(2)*TAN(FI0)
               TP=-SIG0(2)*TAN(FIMU)*(1.D0-AS)+AS*SR
               TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT)
c
               AA1=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU)*K1*SIG0(2)
     $              /SGMT
               AA2=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0))
               AA3=AS*TAN(FI0)
               dTPS1=-(1.D0-AS)*TAN(FIMU)+AA1+AA2-AA3
               dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT
               AA4=-UR*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0)
               AA4=AA4+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2)
c
c       calcul de la matrice Dt equivalente (2D)
c
               KNN=1.D0/AA4
               KNT=0.D0
               KTN=UEQU/UP*KNN*dTPS1
               KTT=TP/UP
c
            END IF
c
c    appel à CALC2 pour calculer les incréments de contraintes
c
            CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS,
     $           DELTA,RI0E,DCON)
c
c    calcul de l' incrément de deformation plastique
c
            DEP1=0.D0
            DEP2=0.D0
c
            IF(UEQU.LT.UP)THEN
c
c       il n' y a pas de plasticité
c
               GO TO 200
c
            ELSE IF((UEQU.GE.UP).AND.(UEQU.LT.UR))THEN
c
c       il y a de la plasticité: trait decroissant de la courbe tau
c
c       test pour vérifier si continuer en plasticité ou reprendre le trait
c       lineaire
c
               EFFE=ABS(SIG0(1))-(TR-TP)/(UR-UP)*(UEQU-UP)-TP
c
               IF(ABS(SIG0(1)).LT.TOL1)THEN
c
                  DTEF=ABS(DCON(1))
c
               ELSE
c
                  DTEF=SIG0(1)*DCON(1)/ABS(SIG0(1))
c
               END IF
c
               DTEF=DTEF-(TR-TP)/(UR-UP)*DELTA
               DTEF=DTEF-((UEQU-UP)/(UR-UP)*(dTRS1-dTPS1)+dTPS1)*DCON(2)
c
               IF((EFFE.LT.0.D0).OR.(DTEF.LT.0.D0))GO TO 200
c
c       definition des quantités necessaires au calcul
c
               VUPU=((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0)
               AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1)
               SR=S0-SIG0(2)*TAN(FI0)
               TP=-SIG0(2)*TAN(FIMU+ATAN(VUPU))*(1.D0-AS)+AS*SR
               TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT)
c
               AA1=(SIG0(2)*(1.D0-AS)*K2)/(SGMT*(COS(FIMU+ATAN(VUPU))**2
     $              ))
               AA2=((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0)
               AA3=1.D0+((1.D0-SIG0(2)/SGMT)**(2.D0*K2))*(TAN(I0)**2)
               AA4=K1*SIG0(2)/SGMT
               AA5=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU+ATAN(VUPU))
     $              *AA4
               AA6=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0))
               AA7=AS*TAN(FI0)
               dTPS1=-(1.D0-AS)*TAN(FIMU+ATAN(VUPU))+AA1*AA2/AA3+AA5+AA6
     $              -AA7
               dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT
               AA8=-UEQU*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0
     $              )
               AA8=AA8+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2)
c
c       calcul de la matrice Dt equivalente (2D)
c
               KNN=1.D0/AA8
               KNT=-((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0)*KNN
               KTN=KNN/(UP-UR)*(dTPS1*(UEQU-UR)+(UP-UEQU)*dTRS1)
               KTT=(TP-TR)/(UP-UR)
               KTT=KTT+KNT/(UP-UR)*(dTPS1*(UEQU-UR)+(UP-UEQU)*dTRS1)
c
c       appel à CALC2 pour calculer les incréments de contraintes
c
               CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS,
     $              DELTA,RI0P,DCON)
c
c       appel à RISPL2 pour calculer les incréments de déformation plastique
c
               CALL RISPL2(RI0E,SIG0,DCON,DU,NSTRS,SGMT,DUp)
c
               DEP1=DUp
               DEP2=0.D0
c
            ELSE IF(UEQU.GE.UR)THEN
c
c       il y a de la plasticité parfaite: trait horizontal de la courbe tau
c
c       test pour vérifier si continuer en plasticité ou reprendre le trait
c       lineaire
c
               EFFE=ABS(SIG0(1))-TR
c
               IF(ABS(SIG0(1)).LT.TOL1)THEN
c
                  DTEF=ABS(DCON(1))
c
               ELSE
c
                  DTEF=SIG0(1)*DCON(1)/ABS(SIG0(1))
c
               END IF
c
               DTEF=DTEF-dTRS1*DCON(2)
c
               IF((EFFE.LT.0.D0).OR.(DTEF.LT.0.D0))GO TO 200
c
c       definition des quantités necessaires au calcul
c
               VUPU=0.D0
               AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1)
               SR=S0-SIG0(2)*TAN(FI0)
               TP=-SIG0(2)*TAN(FIMU)*(1.D0-AS)+AS*SR
               TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT)
c
               AA1=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU)*K1*SIG0(2)
     $              /SGMT
               AA2=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0))
               AA3=AS*TAN(FI0)
               dTPS1=-(1.D0-AS)*TAN(FIMU)+AA1+AA2-AA3
               dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT
               AA4=-UR*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0)
               AA4=AA4+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2)
c
c       calcul de la matrice Dt equivalente (2D)
c
               KNN=1.D0/AA4
               KNT=0.D0
               KTN=KNN*dTRS1
               KTT=0.D0
c
c       appel à CALC2 pour calculer les incréments de contraintes
c
               CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS,
     $              DELTA,RI0P,DCON)
c
c       calcul de l' incrément de deformation plastique
c
               DEP1=DU
               DEP2=0.D0
c
            END IF
c
         ELSE IF(ABS(SIG0(2)).GE.ABS(SGMT))THEN
c
c    cas avec sigma >= SGMT
c
c    il n' y a pas de plasticité
c
c    prédiction élastique
c
c    definition des quantités necessaires au calcul
c
            VUPU=0.D0
            AS=1.D0
            SR=S0-SIG0(2)*TAN(FI0)
            TP=SR
            TR=TP
            dTPS1=-TAN(FI0)
            AA1=((VM**2)*KNI)/((KNI*VM-SIG0(2))**2)
c
c    calcul de la matrice Dt equivalente (2D)
c
            KNN=1.D0/AA1
            KNT=0.D0
            KTN=UEQU/UP*KNN*dTPS1
            KTT=TP/UP
c
c    appel à CALC2 pour calculer les incréments de contraintes
c
            CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS,
     $           DELTA,RI0E,DCON)
c
c    calcul de l' incrément de deformation plastique
c
            DEP1=0.D0
            DEP2=0.D0
c
            IF(UEQU.LT.UP)GO TO 200
c
c       il y a de la plasticité parfaite: trait horizontal de la courbe tau
c
c       test pour vérifier si continuer en plasticité ou reprendre le trait
c       lineaire
c
            EFFE=ABS(SIG0(1))-TP
c
            IF(ABS(SIG0(1)).LT.TOL1)THEN
 
               DTEF=ABS(DCON(1))
c
            ELSE
c
               DTEF=SIG0(1)*DCON(1)/ABS(SIG0(1))
c
            END IF
c
            DTEF=DTEF-dTPS1*DCON(2)
c
            IF((EFFE.LT.0.D0).OR.(DTEF.LT.0.D0))GO TO 200
c
c       definition des quantités necessaires au calcul
c
            VUPU=0.D0
            AS=1.D0
            SR=S0-SIG0(2)*TAN(FI0)
            TP=SR
            TR=TP
            dTPS1=-TAN(FI0)
            AA1=((VM**2)*KNI)/((KNI*VM-SIG0(2))**2)
c
c       calcul de la matrice Dt equivalente (2D)
c
            KNN=1.D0/AA1
            KNT=0.D0
            KTN=KNN*dTPS1
            KTT=0.D0
c
c       appel à CALC2 pour calculer les incréments de contraintes
c
            CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS,
     $           DELTA,RI0P,DCON)
c
c       calcul de l' incrément de deformation plastique
c
            DEP1=DU
            DEP2=0.D0
c
         END IF
c
 200     CONTINUE
c
c    chargement de la matrice de rigidité au debut du pas
c
         IF(KAPPA.EQ.1)THEN
c
            RTT=KTT
            RTN=KTN
c
         END IF
c
c    chargement des valeurs à la fin des sub-pas
c
         SIG0(1)=SIG0(1)+DCON(1)
         SIG0(2)=SIG0(2)+DCON(2)
c
         EPS1=EPS1+DU
         EPS2=EPS2+DV
c
         IF(SIG0(2).GE.-TOL2) THEN
c
c       on a plasticité selon la composante normale: le joint se ouvre
c       ou il reste bloqué
c
            SIG0(1)=0.D0
            SIG0(2)=0.D0
c
            UEQU=UEQU
c
            DEFP(1)=DEFP(1)+DU
            DEFP(2)=DEFP(2)+DV
c
            EPIN0(1)=EPIN0(1)+DU
            EPIN0(2)=EPIN0(2)+DV
c
         ELSE
c
c       on n' a pas de plasticité selon la composante normale: le joint
c       reste fermé
c
            UEQU=ABS(UEQU+DELTA)
c
            DEFP(1)=DEFP(1)+DEP1
            DEFP(2)=DEFP(2)+DEP2
c
            EPIN0(1)=EPIN0(1)+DEP1
            EPIN0(2)=EPIN0(2)+DEP2
c
            DPP1=DEFP(1)
            DPP2=DEFP(2)
c
            EPP1=EPIN0(1)
            EPP2=EPIN0(2)
c
         END IF
c
 300  CONTINUE
c
c **********************************************************************
c ******** fin de la boucle du calcul et preparation des sorties *******
c **********************************************************************
c
c    test pour vérifier si le deux tau sont presque nulles
c
      TOL3=1.D-1*ABS(RTT*(UEQU-VAR0(4))+RTN*DEPST(2))
c
      IF(ABS(SIG0(1)).LT.TOL3)THEN
c
         SIG0(1)=0.D0
c
      END IF
c
c    chargement des valeurs des contraintes finales
c
      SIGF(1)=SIG0(1)
      SIGF(2)=SIG0(2)
c
c    état du joint
c
      IF(EPIN0(2).GT.EPSS)THEN
c
c       le joint est ouvert
c
         VARF(3)=1.D0
c
      ELSE IF((EPIN0(2).LE.EPSS).AND.(SIG0(2).GE.-TOL2))THEN
c
c       le joint est fermé mais sans compression
c
         DEFP(1)=DPP1
         DEFP(2)=DPP2
c
         EPIN0(1)=EPP1
         EPIN0(2)=EPP2
c
         VARF(3)=0.D0
c
      ELSE
c
c       le joint est fermé avec compression
c
         VARF(3)=2.D0
c
      END IF
c
c    chargement de la variable EPSE
c
      VARF(1)=SQRT(2.D0/3.D0*(EPIN0(1)**2+EPIN0(2)**2))
c
c    chargement de la variable EPOU
c
      VARF(2)=DEFP(2)
c
c    chargement de la variable UEQU
c
      VARF(4)=UEQU
c
c    chargement des déformations totales
c
      VARF(5)=EPS1
      VARF(6)=EPS2
c
      RETURN
c
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales