Esope source files

   1 : C VISCOE    SOURCE    PV        08/09/11    21:17:14     6150           
   2 :       SUBROUTINE VISCOE(CM,EPSIPP,EPSE,DTT,DCT,DT
   3 :      &                      ,DELTAT,DTTDT,DCTDT,DTDT)
   4 : C
   5 : C
   6 :       IMPLICIT REAL*8(A-H,O-Z)
   7 : C
   8 :       DIMENSION EPSIPP(3),EPSILT(3),CM(*),DDHOOK(6,6),
   9 :      1 SIGP(3),SIGPT(3),SIGPC(3)
  10 :       PARAMETER (XUnDemi=0.5D0,XZERO=0.D0,UN=1.D0,DEUX=2.D0)
  11 :       PARAMETER (XPETIT=1.D-12)
  12 : 
  13 : C
  14 : C     recuperation des donnees du modele
  15 : C
  16 : 
  17 :         YOUN=CM(1)
  18 :         XNU=CM(2)
  19 :         XMDt=CM(14)
  20 :         XNDt=CM(15)
  21 :         At=CM(16)
  22 :         Bt=CM(17)
  23 :         XMDc=CM(18)
  24 :         XNDc=CM(19)
  25 :         Ac=CM(20)
  26 :         Bc=CM(21)
  27 :         ED0=CM(22)
  28 : C      
  29 : C     Matrice de hooke (pour le 3D)
  30 : C
  31 :         AUX0=YOUN/((Un+XNU)*(Un-XNU-XNU))
  32 :         AUX=AUX0*(Un-XNU)
  33 :         AUX1=AUX0*XNU
  34 :         GEGE=XUnDemi*YOUN/(Un+XNU)
  35 :         DDHOOK(1,1)=AUX
  36 :         DDHOOK(2,1)=AUX1
  37 :         DDHOOK(3,1)=AUX1
  38 :         DDHOOK(1,2)=AUX1
  39 :         DDHOOK(2,2)=AUX
  40 :         DDHOOK(3,2)=AUX1
  41 :         DDHOOK(1,3)=AUX1
  42 :         DDHOOK(2,3)=AUX1
  43 :         DDHOOK(3,3)=AUX
  44 :         DDHOOK(4,4)=GEGE
  45 :         DDHOOK(5,5)=GEGE
  46 :         DDHOOK(6,6)=GEGE
  47 : C
  48 : C
  49 : C      on calcule l'endommagement et les contraintes
  50 : C
  51 : C
  52 : C        on calcule ALFAT ALFAC DT et DC puis D
  53 : C
  54 : C
  55 :         DO 200 ISTRS=1,3
  56 :            SIGP(ISTRS)= XZERO
  57 :            DO 210 JSTRS=1,3
  58 : 
  59 :              SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  60 : 210        CONTINUE
  61 : 200     CONTINUE
  62 : C
  63 : C        on calcule le signe des contraintes elastiques
  64 : C
  65 :         DO 300 ISTRS=1,3
  66 :          IF (SIGP(ISTRS).LT. XZERO) THEN
  67 :             SIGPC(ISTRS)=SIGP(ISTRS)
  68 :             SIGPT(ISTRS)=XZERO
  69 :          ELSE
  70 :             SIGPT(ISTRS)=SIGP(ISTRS)
  71 :             SIGPC(ISTRS)=XZERO
  72 :          END IF
  73 : 300     CONTINUE
  74 :         TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  75 :         TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  76 : C
  77 : C        on calcule les deformations dues aux contraintes positives
  78 : C
  79 :         DO 400 ISTRS=1,3
  80 :           EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  81 : 400     CONTINUE
  82 : C
  83 : C        on en deduit ALFAT et ALFAC
  84 : C
  85 :         ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
  86 :      1          MAX(XZERO,EPSIPP(2))*EPSILT(2) +
  87 :      2          MAX(XZERO,EPSIPP(3))*EPSILT(3)
  88 :         ALFAT = ALFAT/(EPSE*EPSE)
  89 :         ALFAC = UN - ALFAT
  90 : C
  91 : C       modification pour la bi ou tricompression
  92 : C
  93 :         IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN
  94 :           GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+SIGPC(3)*SIGPC(3)
  95 :           GAMMA=-SQRT(GAMMA)/TRSIGC
  96 :           EPSE=EPSE*GAMMA
  97 :         END IF
  98 : C
  99 : C        on calcule DT et DC puis D
 100 :         dpot=((EPSE-ED0-(UN/At)*(dtt/(UN-dtt))**(UN/Bt))/XMdt)**XNdt
 101 :         dpot=MAX( dpot, XZERO)
 102 :         
 103 :         DTTDT=dtt+dpot*DELTAT
 104 :         DTTDT=MIN ( DTTDT , UN-1.D-05 ) 
 105 :         dpoC=((EPSE-ED0-(UN/Ac)*(dct/(UN-dct))**(UN/Bc))/XMdC)**XNdc
 106 :         dpoC=MAX( dpoC, XZERO)
 107 :         DCTDT=dct+dpoC*DELTAT
 108 :         DCTDT=MIN ( DCTDT , UN-1.D-05 )
 109 : C           
 110 :         DTDT = ALFAT*DTTDT + ALFAC*DCTDT
 111 :                 
 112 : 
 113 : C
 114 : C           on borne la valeur de D a 0.99999
 115 : C
 116 :         DTDT=MIN ( DTDT , UN-1.D-05 )
 117 : C
 118 : C           on teste la croissance de D 
 119 : C
 120 :         DTDT=MAX ( DTDT , DT )
 121 : C
 122 :       RETURN
 123 :       END
 124 :  
 125 :  
 126 :