Télécharger deriv1.eso

Retour à la liste

deriv1
C DERIV1    SOURCE    CB215821  17/11/30    21:15:57     9639           
      SUBROUTINE DERIV1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD,
     &DSPT,XMAT,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,YOUNG,NYOUNG,
     &XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI,TPOINT,ITHHER,VARF,
     & IB,IGAU,KERRE)
 
C-----------------------------------------------------------------------
C Objet: Cette subroutine calcule les derivees des variables internes
C        pour un materiau viscoplastique a endommagement et ecrouissage
C        isotrope en regime anisotherme
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
C Entree: TAU pas de temps
C         TI temperature au debut du pas
C         TPOINT  derivee de la temperature
C         SIG(NSTRS) tenseur des contraintes dans le repère général
C         EPSV(NSTRS) tenseur des deformations dans le repère général
C         VAR(NVARI) tableau contenant les variables internes
C           pilotant les equations
C           VAR(1)=p ; VAR(2)=D11 ; VAR(3)=D22;  VAR(4)=D33;
C           VAR(5)=D12;  VAR(6)=D13;  VAR(7)=D23; VAR(8)=DMAX
C         EPSTHD  vitesse de dilatation thermique au debut du pas
C         DSPT(NSTRS,NSTRS) vitesse de deformation totale
C         XMAT(NCOMAT) tableau des parametres scalaires du materiau
C           a une temperature TI donnee
C           XMAT(1)=YOUNG ;XMAT(2)=XNU ;XMAT(3)=N
C           XMAT(4)=M ;XMAT(5)=KK; XMAT(6)=ALF2;
C           XMAT(7)=BET2 ;XMAT(8)=r; XMAT(9)=A ;
C           XMAT(10)=q ;XMAT(11)=ALPHAT
C           XMAT(12)=RHO; XMAT(13)=SIGY
C         MFR indice de la formulation mecanique(seulement massif
C           ou coque pour les materiaux endommageables
C         ICARA nombre de caracteristiques geometriques des elements
C           finis
C         XCAR(ICARA) tableau des caracteristiques geometriques des
C           elements finis
C         IFOURB= -2 EN CONTR. PLANES
C                 -1 EN DEFORM. PLANES
C                  0 EN AXISYMETRIE
C                  1 EN SERIE DE FOURIER
C                  2 EN TRIDIM
C         INDIC=0, 1 OU -1 pour plasticite avec endommagement
C              =2 OU -2 pour viscoplasticite avec endommagement
C         ITHHER = 0 pas de chargement thermique et materiau constant
C                = 1 chargement thermique et materiau constant
C                = 2 chargement thermique et materiau(T)
C------------------------------------------------------------------------
C
C------------------------------------------------------------------------
C Sortie: EPSVD(NSTRS) derivee du tenseur des deformations
C         VARD(NVARI) tableau contenant les derivees des variables
C           internes
C         VARF(NVARI) état des variables internes à la fin du sous programme
C         SIGD(NSTRS) derivee des contraintes
C         DD(NSTRS,NSTRS) matrice de Hooke  au debut du pas
C         DDV(NSTRS,NSTRS) derivee de DD
C         DDINV(NSTRS,NSTRS) inverse de DD
C------------------------------------------------------------------------
C
C
       IMPLICIT INTEGER(I-N)
       IMPLICIT REAL*8(A-H,O-Z)
       DIMENSION SIG(*),EPSV(*),VAR(*),XCAR(*),YOUNG(*),XNU(*)
       DIMENSION SIGD(*),EPSVD(*),VARD(*),DDINV(NSTRS,NSTRS)
       DIMENSION XMAT(*),DSPT(*),EPSTHD(*),VARF(NVARI)
       DIMENSION DDV(NSTRS,NSTRS),DDINVp(NSTRS,NSTRS),DD(NSTRS,NSTRS)
       DIMENSION EPSED(6),XX(6),Di(3),VaP(3),AA(3),BB(6)
       DIMENSION CO(3,3),COM(3,3),DCOT(3,3),COE(3,3),DCOE(3,3)
       DIMENSION AAA(3,3),BBB(3,3),CCC(3,3),UN2(3,3)
       DIMENSION FTH(3,3),ROT(3,3),ROTT(3,3),VcP(3,3),S(3,3)
       DIMENSION D1(3,3),D2(3,3),DV1(3,3),DI1(3,3),DI2(3,3)
       PARAMETER (UN=1.D0, UNM=-1.D0)
       LOGICAL endonul
C
C
       endonul = .FALSE.
       IF ((VAR(2).EQ.0).AND.(VAR(3).EQ.0).AND.(VAR(4).EQ.0).AND.
     &     (VAR(5).EQ.0).AND.(VAR(6).EQ.0).AND.(VAR(7).EQ.0).AND.
     &     (VAR(8).EQ.0)) THEN
         endonul = .TRUE.
       ENDIF
 
C
C-----------------------------------------------------------------
C
C EXPRESSIONS DES TENSEURS UTILISES DANS LE MODELE
C
C-----------------------------------------------------------------
C
C----------------------------------------------
C TENSEUR DES CONTRAINTES (dans repere general)
C----------------------------------------------
C
       CALL ZERO(CO,3,3)
       IF (mfr.eq.1) THEN
       DO I=1,3
         CO(I,I) = SIG(I)
       ENDDO
       CO(1,2) = SIG(4)
       IF(NSTRS.NE.6) GOTO 50
       CO(1,3) = SIG(5)
       CO(2,3) = SIG(6)
50     CONTINUE
       CO(2,1) = CO(1,2)
       CO(3,1) = CO(1,3)
       CO(3,2) = CO(2,3)
       ENDIF
C
C------------------------------------
C TENSEUR D ENDOMMAGEMENT ET
C TENSEUR D ENDOMMAGEMENT DIAGONALISE
C------------------------------------
C
      XD11=VAR(2)
      XD22=VAR(3)
      XD33=VAR(4)
      XD12=VAR(5)
      XD13=VAR(6)
      XD23=VAR(7)
       CALL ZERO(D1,3,3)
       D1(1,1) = XD11
       D1(2,2) = XD22
       D1(3,3) = XD33
       D1(1,2) = XD12
       D1(1,3) = XD13
       D1(2,3) = XD23
       D1(2,1) = D1(1,2)
       D1(3,1) = D1(1,3)
       D1(3,2) = D1(2,3)
 
*
*  AM  14/12/09  TEST SELON NSTRS
*
       IF(IFOURB.LE.0) THEN
           JDIM = 2
           CALL JACOB4(D1,JDIM,Di,ROT)
       ELSE
           JDIM = 3
           CALL JACOB4(D1,JDIM,Di,ROT)
       ENDIF
*
       CALL ZERO(D2,3,3)
       D2(1,1) = Di(1)
       D2(2,2) = Di(2)
       D2(3,3) = Di(3)
       w1 = 1.D0 - Di(1)
       w2 = 1.D0 - Di(2)
       w3 = 1.D0 - Di(3)
       IF (D2(1,1).GE.1.D0) w1=1.D-6
       IF (D2(2,2).GE.1.D0) w2=1.D-6
       IF (D2(3,3).GE.1.D0) w3=1.D-6
 
       ADMAX =  MAX(Di(1),Di(2),Di(3))
C
C------------------------
C TENSEUR UNITE D ORDRE 2
C------------------------
C
       CALL ZERO(UN2,3,3)
       DO I=1,3
         UN2(I,I) = 1.D0
       ENDDO
C
C-------------------------------------------------------
C TENSEUR DES CONTRAINTES MOTRICES (dans repere general)
C-------------------------------------------------------
C
       EQ=XMAT(10)
       DO I=1,3
         DO J=1,3
           ROTT(I,J)=ROT(I,J)
         ENDDO
       ENDDO
       CALL INVALM (ROTT, 3 , 3 , KER, 1.D-12)
C
      IF (endonul) THEN
        DO I=1,3
          DO J=1,3
            COM(I,J) = CO(I,J)
          ENDDO
        ENDDO
      ELSE
C --- Matrice de [(I-D)/det(I-D)]**(q/2) dans le repère principal de D ---
        CALL ZERO(DI1,3,3)
        DO I=1,3
          XNB=(1.D0-Di(I))/(w1*w2*w3)
          IF (XNB.GE.0.D0) THEN
            DI1(I,I)=(XNB)**(EQ/2.D0)
          ELSE
            DI1(I,I)=0.D0
              WRITE(6,*) 'erreur1, IB = ', IB, ' Igau=',igau
              KERRE = 99
              GO TO 100
c          Do j = 1,7
c          write(6,*) 'var(',j,')=', var(j)
c          ENDDO
          ENDIF
        ENDDO
C --- passage dans le repère de général ---
        CALL MULMAT(AAA,ROT,DI1,3,3,3)
        CALL MULMAT(DI2,AAA,ROTT,3,3,3)
 
        CALL MULMAT(AAA,DI2,CO,3,3,3)
        CALL MULMAT(COM,AAA,DI2,3,3,3)
      ENDIF
C
C-----------------------------------
C TENSEUR DES CONTRAINTES EFFECTIVES
C-----------------------------------
C
      IF (endonul) THEN
        DO I=1,3
          DO J=1,3
            COE(I,J) = CO(I,J)
          ENDDO
        ENDDO
      ELSE
C --- Matrice de [I-D]**(-1/2) dans le repère principal de D ---
        CALL ZERO(DI1,3,3)
        CALL ZERO(DI2,3,3)
        DO I=1,3
          XNB=1.D0-Di(I)
          IF (XNB.GT.0.D0) THEN
            DI1(I,I)=SQRT(1.D0/(XNB))
          ELSE
            DI1(I,I)=0.D0
            WRITE(6,*) 'erreur2 IB=', IB, '  IGAU= ', IGAU
            KERRE=99
            GO TO 100
          ENDIF
        ENDDO
C --- passage dans le repère de général ---
        CALL MULMAT(AAA,ROT,DI1,3,3,3)
        CALL MULMAT(DI2,AAA,ROTT,3,3,3)
c        DO i=1,3
c        DO j=1,3
c        WRITE(6,*) 'DI2(',i,',',j,')=',DI2(i,j)
c        enddo
c        enddo
 
        CALL MULMAT(AAA,DI2,CO,3,3,3)
        CALL MULMAT(COE,AAA,DI2,3,3,3)
      ENDIF
C
C----------------------------------
C DEVIATEUR DE CONTRAINTE EFFECTIVE
C----------------------------------
C
        CALL ZERO(BBB,3,3)
        DO I=1,3
          BBB(I,I)=-(COE(1,1)+COE(2,2)+COE(3,3))/(3.D0)
        ENDDO
        CALL PLMATR(COE,UN,BBB,UN,3,3,DCOE)
C --- Calcul de (I-D)**(-1/2).DCOE.(I-D)**(-1/2) ---
      IF (endonul) THEN
        DO I=1,3
          DO J=1,3
            DCOT(I,J) = DCOE(I,J)
          ENDDO
        ENDDO
      ELSE
        CALL MULMAT(AAA,DI2,DCOE,3,3,3)
        CALL MULMAT(DCOT,AAA,DI2,3,3,3)
      ENDIF
C
C----------------------------------------------------------------------
C
C CALCUL DE J1,J2
C   J1: CONTRAINTE MOYENNE
C   J2: CONTRAINTE DE VON MISES
C
C----------------------------------------------------------------------
C
       AJ1 = (SIG(1)+SIG(2)+SIG(3))/(3.D0)
       DO I=1,3
         XX(I) = SIG(I)-AJ1
       ENDDO
       DO I=4,NSTRS
         XX(I) = SIG(I)
       ENDDO
       AJ2 =  PROCON(XX,XX,NSTRS)
       AJ2 =  SQRT(1.5D0*AJ2)
C
C----------------------
C FORCE THERMODYNAMIQUE
C----------------------
C
C --- Etat de triaxialité de la contrainte ---
       IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN
         tri=1.D0
       ELSE
         tri = (AJ2)/(SIG(1)+SIG(2)+SIG(3))
       ENDIF
c           write(6,*) 'AJ2=  ',AJ2
c           write(6,*) 'SIG(1)=  ',SIG(1)
c           write(6,*) 'SIG(2)=  ',SIG(2)
c           write(6,*) 'SIG(3)=  ',SIG(3)
 
c           write(6,*) 'tri=  ',tri
 
C --- Valeurs propres et vecteurs propres du tenseur des contraintes motrices ---
       CALL JACOB3(COM,3,VaP,VcP)
       ER=XMAT(8)
       EA=XMAT(9)
       DO I=1,3
         AA(I) = ((tri*(VaP(I))/(EA))+ABS(tri*(VaP(I))/(EA)))/(2.D0)
         AA(I) = (AA(I))**(ER)
       ENDDO
c       DO I=1,3
c         DO J=1,3
c           write(6,*) 'COM(',i,',',j,')= ',COM(I,J)
c         ENDDO
c       ENDDO
C         DO J=1,3
C           write(6,*) 'VaP(',j,')= ',VaP(J)
C         ENDDO
c         DO J=1,3
c           write(6,*) 'AA(',j,')= ',AA(J)
c         ENDDO
c       DO I=1,3
c         DO J=1,3
c           write(6,*) 'VcP(',i,',',j,')= ',VcP(I,J)
c         ENDDO
c       ENDDO
C
C -- Dans le repère général ---
       CALL ZERO(FTH,3,3)
       FTH(1,1) = AA(1)*(VcP(1,1))**2+AA(2)*(VcP(1,2))**2
       FTH(1,1) = FTH(1,1)+AA(3)*(VcP(1,3))**2
       FTH(2,2) = AA(1)*(VcP(2,1))**2+AA(2)*(VcP(2,2))**2
       FTH(2,2) = FTH(2,2)+AA(3)*(VcP(2,3))**2
       FTH(3,3) = AA(1)*(VcP(3,1))**2+AA(2)*(VcP(3,2))**2
       FTH(3,3) = FTH(3,3)+AA(3)*(VcP(3,3))**2
       FTH(1,2) = AA(1)*VcP(1,1)*VcP(2,1)+AA(2)*VcP(1,2)*VcP(2,2)
       FTH(1,2) = FTH(1,2)+AA(3)*VcP(1,3)*VcP(2,3)
       FTH(1,3) = AA(1)*VcP(1,1)*VcP(3,1)+AA(2)*VcP(1,2)*VcP(3,2)
       FTH(1,3) = FTH(1,3)+AA(3)*VcP(1,3)*VcP(3,3)
       FTH(2,3) = AA(1)*VcP(2,1)*VcP(3,1)+AA(2)*VcP(2,2)*VcP(3,2)
       FTH(2,3) = FTH(1,3)+AA(3)*VcP(2,3)*VcP(3,3)
       FTH(3,1) = FTH(1,3)
       FTH(2,1) = FTH(1,2)
       FTH(3,2) = FTH(2,3)
c       DO I=1,3
c         DO J=1,3
c           write(6,*) 'FTH(',i,',',i,')= ',FTH(I,i)
c         ENDDO
c       ENDDO
C
C---------------------------------------------------------------------
C
C CALCUL DE LA VITESSE DE DEFORMATION PLASTIQUE CUMULEE VARD(1)
C
C---------------------------------------------------------------------
C
C --- Paramètres matériau ---
       EM=XMAT(4)
       EN=XMAT(3)
       EK=XMAT(5)
       ALF2=XMAT(6)
C --- Calcul de la contrainte équivalente effective ---
       AJ3 = (COE(1,1)+COE(2,2)+COE(3,3))/(3.D0)
       DO I=1,3
         XX(I) = COE(I,I)-AJ3
       ENDDO
       XX(4) = COE(1,2)
       XX(5) = COE(1,3)
       XX(6) = COE(2,3)
       AJ4 =  PROCON(XX,XX,NSTRS)
       AJ4 =  SQRT(1.5D0*AJ4)
C
C --- Calcul de la norme euclidienne ---
       DO I=1,3
         DO J=1,3
           AAA(I,J)=(DCOT(I,J)*3.D0)/(2.D0*AJ4)
         ENDDO
       ENDDO
       CALL PLMATR(UN2,UN,D1,UNM,3,3,BBB)
       CALL INVALM (BBB, 3 , 3 , KER, 1.D-12)
       DO I=1,3
         DO J=1,3
           BBB(I,J)=BBB(I,J)*ALF2/3.D0
         ENDDO
       ENDDO
       CALL ZERO(CCC,3,3)
       CALL PLMATR(AAA,UN,BBB,UN,3,3,CCC)
c       DO I=1,3
c       DO j=1,3
c         write(6,*) 'CCC(',i,',',j,')= ', CCC(i,j)
c       ENDDO
c       ENDDO
 
       val = 0.D0
       DO I=1,3
         DO J=1,3
           val = val+(CCC(I,J))**(2.D0)
         ENDDO
       ENDDO
       xnor = SQRT(val)
C       WRITE(6,*) 'xnor=',xnor
C --- Calcul du premier terme du produit ---
       IF ( VAR(1).EQ.0.D0 ) VAR(1)=VAR(1)+1.D-12
       ter = (AJ4+(ALF2*AJ3))/(EK*(VAR(1)**(1.D0/EM)))
       ter = ter + ABS((AJ4+(ALF2*AJ3))/(EK*(VAR(1)**(1.D0/EM))))
       ter = ter/(2.D0)
       ter = ter**(EN)
       ter = ter*SQRT(2.D0/3.D0)
c       WRITE(6,*) 'ter=',ter
C --- Calcul de l evolution de p ---
       IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN
         VARD(1) = 0.D0
       ELSE
         VARD(1) = ter*xnor
       ENDIF
C       write(6,*) 'vard(1)=', vard(1)
C
C---------------------------------------------------------------------
C
C CALCUL DE LA VARIATION DE L ENDOMMAGEMENT VARD(I) (I=2..7)
C
C---------------------------------------------------------------------
 
       BET2=XMAT(7)
C --- Dans le repère principal de D ---
       CALL ZERO(DV1,3,3)
       DV1(1,1) = (BET2*FTH(1,1))+FTH(2,2)+FTH(3,3)
       DV1(2,2) = (BET2*FTH(2,2))+FTH(1,1)+FTH(3,3)
       DV1(3,3) = (BET2*FTH(3,3))+FTH(2,2)+FTH(1,1)
       DV1(1,2) = (BET2-1.D0)*FTH(1,2)
       DV1(1,3) = (BET2-1.D0)*FTH(1,3)
       DV1(2,3) = (BET2-1.D0)*FTH(2,3)
       DV1(2,1) = DV1(1,2)
       DV1(3,1) = DV1(1,3)
       DV1(3,2) = DV1(2,3)
C
C --- Expressions des VARD(I) (I=2..7) ---
       DO I=2,4
         VARD(I) = DV1(I-1,I-1)
       ENDDO
       VARD(5) = DV1(1,2)
       VARD(6) = DV1(1,3)
       VARD(7) = DV1(2,3)
c       DO I=2,4
c         write(6,*) 'vard(2)= ', vard(2)
c       ENDDO
C
C---------------------------------------------------------------------
C
C CALCUL DE LA VARIATION DE DEFORMATION PLASTIQUE
C
C---------------------------------------------------------------------
 
       IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN
         DO I=1,NSTRS
           EPSVD(I) = 0.D0
         ENDDO
       ELSE
         DO I=1,3
           EPSVD(I) = CCC(I,I)*VARD(1)*SQRT(3.D0/2.D0)/xnor
         ENDDO
           EPSVD(4) = CCC(1,2)*VARD(1)*SQRT(3.D0/2.D0)/xnor
*
*  AM  CORRECTION LE 9/12/09
*
****     IF (IFOUR.GT.0) THEN
         IF (IFOURB.GT.0) THEN
           EPSVD(5) = CCC(1,3)*VARD(1)*SQRT(3.D0/2.D0)/xnor
           EPSVD(6) = CCC(2,3)*VARD(1)*SQRT(3.D0/2.D0)/xnor
         ENDIF
       ENDIF
 
 
c       DO I=1,3
c         write(6,*) 'EPSVD(',i,')= ', EPSVD(i)
c       ENDDO
C
C----------------------------------------------------------------------
C
C CALCUL DE LA VARIATION DES CONTRAINTES
C
C----------------------------------------------------------------------
 
      CALL DERTRA(NYOUNG,YOUNG,TI,YUNG,YUNGV,TINF,TSUP)
      CALL DERTRA(NXNU,XNU,TI,ENU,ENUV,TINF,TSUP)
C
C --- Calcul de l inverse de la matrice de Hooke au debut du pas ---
      CALL ELAST4(2,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
     &XCAR,ICARA,MFR,NSTRS,DDINV,DDV,KERRE,INDIC,ITHHER)
      IF(KERRE.NE.0) GO TO 100
C
C --- Calcul de la matrice de Hooke et de sa derivee au debut du pas ---
      CALL ELAST4(1,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
     &XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER)
      IF(KERRE.NE.0) GO TO 100
C
C --- Détermination de la variation de la matrive de Hooke inversée ---
      CALL HOOKP(VAR,VARD,NVARI,XMAT,NCOMAT,MFR,NSTRS,DDINVp)
C
C --- Calcul de la variation de contrainte ---
      CALL ZDANUL(SIGD,NSTRS)
      DO I=1,3
        XX(I) = CO(I,I)
      ENDDO
      XX(4) = CO(1,2)
      XX(5) = CO(1,3)
      XX(6) = CO(2,3)
 
      CALL ZDANUL(EPSED,NSTRS)
      DO I=1,NSTRS
        EPSED(I)=DSPT(I)-EPSVD(I)-EPSTHD(I)
C        write(6,*) 'DSPT(',I,') =', DSPT(I)
c        write(6,*) 'EPSVD(',I,') =', EPSVD(I)
c        write(6,*) 'EPSTHD(',I,') =', EPSTHD(I)
C        write(6,*) 'EPSED(',I,') =', EPSED(I)
      ENDDO
 
      CALL MULMAT(BB,DDINVp,XX,NSTRS,1,NSTRS)
 
      DO I=1,NSTRS
        XX(I)=EPSED(I)-BB(I)
      ENDDO
 
      CALL MULMAT(SIGD,DD,XX,NSTRS,1,NSTRS)
c      DO I=1,3
c       IF (iGAU.EQ.1) THEN
c       write(6,*) 'SIGD(1) =', SIGD(1)
c       ENDIF
c      ENDDO
C
C----------------------------------------------------------------------
C
C MISE A JOUR DES VARIABLES INTERNES
C
C----------------------------------------------------------------------
 
      VARF(8)=ADMAX
 
100   CONTINUE
      RETURN
      END
 
 
 
 
 
 

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