C GERNAY SOURCE MB234859 23/02/03 21:15:09 11581 SUBROUTINE GERNAY(WRK52,WRK53,WRK54,IB,IGAU,NBPGAU) * *______________________________________________________________________ * * MODELE GERNAY *______________________________________________________________________ * IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC DECHE * PARAMETER(UN=1.D0,DEUX=2.D0,UNDEMI=.5D0) * REAL*8 YOUN, NU, BULK3, G2, LAM, & TETA0, DTETA, TETAF, ALPH, TREF, DYDT, DNDT, & YOUN0, NU0, BULK30, G20, LAM0, & YOUNF, NUF, BULK3F, G2F, LAMF, & DELDSE(6,6), ETHERM(6), DTHERM(6), & YREF, NREF, EPSE0(6), DEPSE(6),DHOOK(6,6), & DSIG2(3,3),DEPST2(3,3),TENS(3,3),DSIG(6), & V1(3),V2(3),V3(3),V(3,3),P(3,3) REAL*8 PARAHOT3(83,1),sigma6v(6), deps6v(6),STRESS(6),DSTRAN(6) REAL*8 STRAN(6),STATEV(32),PROPS(10),DDSDDE(6,6) character*10 cmat logical lerror,lcp * * ON EXPLORE LE DECHE * NPT=1 NOEL=1 KERRE=0 NMAT = XMAT(/1) NVARI = VAR0(/1) NSTRS = SIG0(/1) CALL ZERO(STRESS,6,1) CALL ZERO(STRAN,6,1) CALL ZERO(DSTRAN,6,1) CALL ZERO(STATEV,32,1) CALL ZERO(PARAHOT3,83,1) * * TRANSFER DE L'ETAT INITIAL, DU DEPST, ET DU MATERIAU * DO I=1,NSTRS STRESS(I) = SIG0(I) STRAN(I) = EPST0(I) DSTRAN(I)= DEPST(I) ENDDO * DO I=1,NVARI STATEV(I) = VAR0(I) ENDDO * PROPS( 1) = XMAT( 2) PROPS( 2) = XMAT( 5) PROPS( 3) = XMAT( 6) PROPS( 4) = XMAT( 7) PROPS( 5) = XMAT( 8) PROPS( 6) = XMAT( 9) PROPS( 7) = XMAT(10) PROPS( 8) = XMAT(11) PROPS( 9) = XMAT(12) PROPS(10) = XMAT( 1) C C----------------------------------------------------------------------- C Constantes de materiau C ---------------------- C PROPS(1) : NU Coefficient de Poisson c PROPS(2) : RC Resistance à la compression c PROPS(3) : RT Resistance à la traction c PROPS(4) : RB Resistance à la compression biaxiale c PROPS(5) : EPS1 Deformation au pic en compression c PROPS(6) : DILA Parametre de dilatance alpha_g c PROPS(7) : DAM1 Endommagement de compression au pic c PROPS(8) : XC Paramètre adimensionnel, rapport de l'energie de rupture c en compression dissipée avant le pic sur l'énergie totale c PROPS(9) : GT Densite d'energie de rupture en traction c PROPS(10) : YOUN Young pour le calcul de la rigidite elastique c C----------------------------------------------------------------------- C * print *,' ' * print *, ' POINT DINTEGRATION IGAU = ', IGAU * print *, ' stress', stress * print *, 'dstran', dstran * print *,' STATEV 1 à 4' , (statev(i),i=1,4) * print *,' STATEV 5 à 10 epspl ', (statev(i),i=5,10) * print *,' STATEV 11 à 16 epstr ', (statev(i),i=11,16) * print *,' STATEV 17 à 22 sig,eff,c ', (statev(i),i=17,22) * print *,' STATEV 24 à 29 sig,eff ', (statev(i),i=24,29) * print *,' STATEV 30 à 32 tmax,phi,epsth ',(statev(i),i=30,32) cmat = 'SICO_ETC3D' idimpara3 = 83 c ATTENTION si idimpara3 est modifié il faut modifier aussi la dimension c de parahot3 (plus haut) lerror = .false. **** npttot3=30 npttot3=1 * lcp = .false. * * CAS CONTRAINTES PLANES * if (IFOUR.EQ.-2) then lcp = .true. endif * TEMPERATURES EN DEBUT ET FIN DE PAS TETA1 = ture0(1) TETA2 = turef(1) C Recuperation des constantes de materiau C PARAHOT3(2,NPT) = PROPS(1) PARAHOT3(3,NPT) = PROPS(2) PARAHOT3(4,NPT) = PROPS(3) PARAHOT3(7,NPT) = PROPS(4) PARAHOT3(5,NPT) = PROPS(5) PARAHOT3(8,NPT) = PROPS(6) PARAHOT3(9,NPT) = PROPS(7) PARAHOT3(10,NPT)= PROPS(8) PARAHOT3(11,NPT)= PROPS(9) C c remplissage des variables internes c accumulated plastic strain kappa at the beginning of the step c in tension and in compression PARAHOT3(idimpara3-37,NPT) = STATEV(1) PARAHOT3(idimpara3-36,NPT) = STATEV(2) c damage variable in tension and in compression PARAHOT3(idimpara3-41,NPT) = STATEV(3) PARAHOT3(idimpara3-40,NPT) = STATEV(4) c plastic strains do i=1,6 PARAHOT3(idimpara3-25+i,NPT) = STATEV(4+i) end do c transient creep strain at the beginning of the step do i=1,6 PARAHOT3(17+i,NPT) = STATEV(10+i) end do c effective stress in compression (used for transient creep strain calculation) do i=1,6 PARAHOT3(23+i,NPT) = STATEV(16+i) end do do i=1,6 PARAHOT3(29+i,NPT) = PARAHOT3(23+i,NPT) end do c effective stress after convergence (used for sigf3d forward calculation) do i=1,6 PARAHOT3(idimpara3-7+i,NPT) = STATEV(23+i) end do c remplissage des contraintes c remplissage du vecteur contraintes et du vecteur c increment de deformation if (.not.lcp) then do i=1,6 PARAHOT3(idimpara3-13+i,NPT) = STRESS(i) sigma6v(i)= STRESS(i) deps6v(i)= DSTRAN(i) end do else do i=1,4 PARAHOT3(idimpara3-13+i,NPT) = STRESS(i) sigma6v(i)= STRESS(i) deps6v(i)= DSTRAN(i) end do PARAHOT3(idimpara3-8,NPT)=0.d0 PARAHOT3(idimpara3-7,NPT)=0.d0 sigma6v(5) = 0.d0 sigma6v(6) = 0.d0 deps6v(5) = 0.d0 deps6v(6) = 0.d0 endif c remplissage des deformations au debut du pas t c et en fin de pas if (.not.lcp) then do i=1,6 PARAHOT3(idimpara3-31+i,NPT) = STRAN(i) PARAHOT3(idimpara3-19+i,NPT) = STRAN(i) + DSTRAN(i) end do else do i=1,4 PARAHOT3(idimpara3-31+i,NPT) = STRAN(i) PARAHOT3(idimpara3-19+i,NPT) = STRAN(i) + DSTRAN(i) end do endif c TEMPERATURES c ******************************************************************************** c probleme pour vecteur STATEV - STATEV(30,31,32) ne passe pas d'un pas à l'autre c => il faut passer le vecteur STATEV à 32 composantes au lieu de 29 c ******************************************************************************** c actualisation de la temperature au pas courant **** TEMP = TEMP + DTEMP **** PARAHOT3(idimpara3-39,NPT) = TEMP PARAHOT3(idimpara3-39,NPT) = TETA2 c actualisation éventuelle de la temperature maximale atteinte **** if (TEMP.gt.STATEV(30)) STATEV(30) = TEMP if (TETA2.gt.STATEV(30)) STATEV(30) = TETA2 PARAHOT3(idimpara3-38,NPT) = STATEV(30) c phi function PARAHOT3(17,NPT) = STATEV(31) c Thermal strain PARAHOT3(idimpara3,NPT) = STATEV(32) c appel a mater2 NTOT = 1 CALL MATER2TG (cmat,NTOT, lerror,idimpara3, parahot3, npttot3, & ddsdde, sigma6v, deps6v,lcp,NOEL,NPT,Kinc) * print *, ' sortie de mater2tg' CALL MATER3TG (cmat,npt,idimpara3,parahot3,npttot3,lcp, . NOEL,NPT) * print *, ' sortie de mater3tg' C sauvetage des resultats c remplissage des contraintes if (.not.lcp) then do i=1,6 STRESS(i) = sigma6v(i) end do else do i=1,4 STRESS(i) = sigma6v(i) end do endif c if ((NOEL.eq.18).and.(NPT.eq.23)) then * print *,'GERNAY CONTRAINTES EN FIN DE PAS stress = ', stress c endif c actualisation des variables internes c accumulated plastic strain kappa at the beginning of the step c in tension and in compression STATEV(1) = PARAHOT3(idimpara3-37,NPT) STATEV(2) = PARAHOT3(idimpara3-36,NPT) c damage variable in tension and in compression STATEV(3) = PARAHOT3(idimpara3-41,NPT) STATEV(4) = PARAHOT3(idimpara3-40,NPT) c plastic strains do i=1,6 STATEV(4+i) = PARAHOT3(idimpara3-25+i,NPT) end do c transient creep strain at the beginning of the step do i=1,6 STATEV(10+i) = PARAHOT3(17+i,NPT) end do c effective stress in compression (used for transient creep strain calculation) do i=1,6 STATEV(16+i) = PARAHOT3(23+i,NPT) end do c equivalent strain STATEV(23) = STATEV(1) + STATEV(2) c effective stress after convergence (used for sigf3d forward calculation) do i=1,6 STATEV(23+i) = PARAHOT3(idimpara3-7+i,NPT) end do c maximum temperature STATEV(30) = PARAHOT3(idimpara3-38,NPT) c phi function STATEV(31) = PARAHOT3(17,NPT) c Thermal strain STATEV(32) = PARAHOT3(idimpara3,NPT) c do iloc=1,NTENS do iloc=1,4 test = ABS(stress(iloc)) if (test.gt.10.d20) then print *,'>>>>>>>>>>>>>>>>>>>> problem stress = ', . stress print *,'in NOEL = , NPT = ', IB , IGAU KINC = -2 endif enddo if (KINC.eq.-2) then print *,'in NOEL = , NPT = ', IB , IGAU endif * * TRANSFERT DE L'ETAT FINAL * DO I=1,NSTRS SIGF(I) = STRESS(I) ENDDO * DO I=1,NVARI VARF(I) = STATEV(I) ENDDO * RETURN END