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 <= 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