Sources Esope | Cast3M

deriv
C DERIV     SOURCE    CB215821  18/09/21    21:15:34     9930           
      SUBROUTINE DERIV(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD,
     &DSPT,XMAT,FKX,NPTEUR,NSTRS,NVARI,NCOMAT,NKX,NC,DD,DDV,DDINV,
     &YOUNG,NYOUNG,XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,
     &TI,TPOINT,ITHHER,TRUC,NCOURB)
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
C         EPSV(NSTRS) tenseur des deformations
C         VAR(NVARI) tableau contenant les variables internes
C           pilotant les equations
C           VAR(1)=p ; VAR(2)=r ;VAR(3)=D
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)=ALPHA
C           XMAT(7)=BETA ;XMAT(8)=r; XMAT(9)=A ;
C           XMAT(10)=EVOL ;XMAT(11)=ALPHAT
C           XMAT(12)=RHO; XMAT(13)=SIGY
C         FKX(NKX) tableau de la fonction k(X,T) contenant les courbes
C           de k a differentes temperatures T1,T2,T3,...
C         NPTEUR(NC) tableau de pointeurs sur les differentes courbes
C           de la fonction k
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         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
      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 FKX(*),NPTEUR(*),XMAT(*),DSPT(*),TRUC(*)
       DIMENSION EPSTHD(*),DD(NSTRS,NSTRS),DDV(NSTRS,NSTRS)
       DIMENSION XX(6),A(3,3),D(3),S(3,3)
C
C----------------------------------------------------------------------
C CALCULATE  J0,J1,J2
C   J0: THE BIGEST PRINCIPAL STRESS
C   J1:  MEAN STRESS
C   J2: CONTR. DE VON MISES
C----------------------------------------------------------------------
*  cas des massif et coques
*      write(6,*) ' entree dans deriv mfr ',mfr
      if(mfr.eq.1.or.mfr.eq.3) then
      AJ1  = SIG(1)+SIG(2)+SIG(3)
       AJ1  = AJ1/3.0D0
       DO 30 I=1,3
30       XX(I) = SIG(I)-AJ1
       DO 40 I=4,NSTRS
40       XX(I) = SIG(I)
       AJ2 =  PROCON(XX,XX,NSTRS)
       AJ2 =  SQRT(1.5D0*AJ2)
      elseif (mfr.eq.5) then
*  cas des coques epaisses &lt;= a priori non utilise dans ce modele !
       alfah = 1.D0
       ep1 = xcar(1)
       AJ1=   (SIG(1)+SIG(2)) / 3.D0
       AJ2 = vonmis0(sig,nstrs,mfr,ifourb,ep1,alfah)
       XX(1)= sig(1) - aj1
       xx(2) = sig(2)-aj1
       DO 41 I=3,NSTRS
41       XX(I) = SIG(I)
      endif
       CALL ZERO(A,3,3)
       if(mfr.eq.1.or.mfr.eq.3) then
       DO 50 I=1,3
50       A(I,I) = SIG(I)
       A(1,2) = SIG(4)
       IF(NSTRS.NE.6) GOTO 54
       A(1,3) = SIG(5)
       A(2,3) = SIG(6)
54     CONTINUE
      elseif(mfr.eq.5) then
         a(1,1)=sig(1)
         a(2,2)=sig(2)
         a(1,2)= sig(3)
         a(1,3)= sig(4)
         a(2,3)= sig(5)
       endif
       A(2,1) = A(1,2)
       A(3,1) = A(1,3)
       A(3,2) = A(2,3)
       IF (NSTRS.EQ.4) CALL JACOB3(A,2,D,S)
*  on appele jacod3 car on n'a pas besoin des vecteurs propres
*       write(6,*) ' entree de jacob3'
       IF (NSTRS.GE.5) CALL JACOD3(A,3,D)
*       write(6,*) ' sortie de jacob3'
C
C---------------------------------------------------------------------
C CALCULATE PRINCIPLE STRESS : D(*)
C---------------------------------------------------------------------
       AJ0 =  MAX(D(1),D(2),D(3))
C
C---------------------------------------------------------------------
C CALCULATE VARIABLE X
C---------------------------------------------------------------------
C
       ALFA = XMAT(6)
       BETA = XMAT(7)
       X    = (ALFA*AJ0)+(BETA*AJ1)+((1.0D0-ALFA-BETA)*AJ2)
C
C---------------------------------------------------------------------
C CALCULATE PARAMETER k(X,T)=AK
C---------------------------------------------------------------------
C
       CALL VALPAR(FKX,NPTEUR,NKX,NC,X,AK,AKV,TI,TINF,TSUP,TRUC,NCOURB)
 
C
C---------------------------------------------------------------------
C CALCULATE RATE OF EQUAL EFFECTIVE STRAIN VARD(1)
C           RATE OF REAL EQUAL EFFECTIVE STRAIN VARD(2)
C           WITH VARD(2)=(1-D)*VARD(1)
C---------------------------------------------------------------------
       DEN=1.D0-VAR(3)
       EM=XMAT(4)
       EN=XMAT(3)
       EK=XMAT(5)
       SIGY=XMAT(13)
       IF ( VAR(2).EQ.0.D0 ) VAR(2)=VAR(2)+1.D-12
       VARR=AJ2-(SIGY*DEN)
       D1=DEN*EK*(VAR(2)**(1.D0/EM))
       VARR=VARR/D1
       IF ( VARR.GT.0.D0 ) THEN
           VARD(2)=VARR**EN
         ELSE
           VARD(2)=0.D0
       ENDIF
 
       VARD(1)=VARD(2)/DEN
C
C---------------------------------------------------------------------
C CALCULATE RATE OF DAMAGE
C---------------------------------------------------------------------
       ER=XMAT(8)
       EA=XMAT(9)
       IF (X.LE.0.D0) THEN
           VARD(3)=0.D0
         ELSE
           VARD(3)=((X/EA)**ER)/(DEN**AK)
       ENDIF
C
C---------------------------------------------------------------------
C CALCULATE RATE OF PLASTIC STRAIN
C---------------------------------------------------------------------
       IF ( AJ2.EQ.0.D0 ) THEN
           AA3=0.D0
         ELSE
           AA3=1.5D0*VARD(1)/AJ2
       ENDIF
       if(mfr.eq.1.or.mfr.eq.3) then
       DO 10 I=1,NSTRS
         AGAM=1.D0
         IF (I.GE.4) AGAM=2.D0
         EPSVD(I)=AA3*XX(I)*AGAM
10     CONTINUE
      elseif(mfr.eq.5) then
       DO 11 I=1,NSTRS
         AGAM=1.D0
         IF (I.GE.3) AGAM=2.D0
         EPSVD(I)=AA3*XX(I)*AGAM
11     CONTINUE
      endif
C
C----------------------------------------------------------------------
C CALCULATE RATE OF STRESS
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 ELAST1(2,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
     &XCAR,ICARA,MFR,NSTRS,DDINV,DDV,KERRE,INDIC,ITHHER)
C
C Calcul de la matrice de Hooke et de sa derivee au debut du pas
      CALL ELAST1(1,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
     &XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER)
C
      CALL ZDANUL(SIGD,NSTRS)
      DO I=1,NSTRS
        XFLO=( DSPT(I)-EPSVD(I)-EPSTHD(I) )
        DO J=1,NSTRS
          SIGD(J)=SIGD(J)+ DD(J,I)*XFLO
        ENDDO
      ENDDO
 
      CALL MULMAT(DD,DDV,DDINV,NSTRS,NSTRS,NSTRS)
 
      XFLO1 = VARD(3)/DEN
      DO I=1,NSTRS
        XFLO=SIG(I)
        DO J=1,NSTRS
          SIGD(J)=SIGD(J) + TPOINT*DD(J,I)*XFLO
        ENDDO
        SIGD(I)=SIGD(I) - XFLO1*SIG(I)
      ENDDO
 
      RETURN
      END