C HUJEUX SOURCE FANDEUR 22/05/02 21:15:24 11359 SUBROUTINE HUJEUX(SIG0,NSTRS,DEPST,VAR0,NVARI,XMAT,NCOMAT,XCAR, . SIGF,VARF,DEFP,PRECIS,MFR,KERRE) C---------------------------------------------------------------------- C PLASTICITE MODELE HUJEUX C C ENTREES C SIG0(NSTRS) = CONTRAINTES INITIALES C NSTRS = NOMBRE DE CONTRAINTES C DEPST(NSTRS) = INCREMENT DE DEFORMATIONS TOTALES C VAR0(NVARI) = VARIABLES INTERNES DEBUT C VAR0(1) = EPSE C VAR0(2) = TRAP C XMAT(NCOMAT) = COMPOSANTES DE MATERIAU C XMAT(1) = YOUN C XMAT(2) = NU C XMAT(3) = RHO C XMAT(4) = ALPH C XMAT(5) = PI0 C XMAT(6) = COHE C XMAT(7) = M C XMAT(8) = BETA C XMAT(9) = HUA C XMAT(10) = HUB C XMAT(11) = HUR C XMAT(12) = P1 C XMAT(13) = N C XMAT(14) = E1 C NCOMAT = NOMBRE DE COMPOSANTES DE MATERIAU C XCAR(ICARA) = CARACTERISTIQUES C MFR = NUMERO DE LA FORMULATION DE L'ELEMENT FINI C = 1 MASSIF C C SORTIES C SIGF(NSTRS) = CONTRAINTES FINALES C VARF(NVARI) = VARIABLES INTERNES FINALES C DEFP(NSTRS) = DEFORMATIONS PLASTIQUES C KERRE = 0 TOUT OK C 1 SI DLAMBDA NEGATIF C 2 NOMBRE MAX D ITERATIONS INTERNES DEPASSE C 21 ON NE TROUVE PAS L INTERSECTION AVEC LA SURFACE DE CHARGE C 22 SIG0 A L EXTERIEUR DE LA SURFACE DE CHARGE C C----------------------------------------------------------------------- C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO C DIMENSION SIG0(*),DEPST(*),VAR0(*),XMAT(*),XCAR(*) DIMENSION SIGF(*),VARF(*),DEFP(*) DIMENSION STOT(6),XINV(2),SIGEL(6) DIMENSION DFDS(6),DGDS(6),COVNMS(1),DEFPLA(6),DDEFPL(6) REAL*8 DDOT EXTERNAL DDOT DATA ITMAX/25/ ALFAH=0.D0 COVNMS(1) = 0.D0 * * TEST SUR IFOUR : ON ACCEPTE 2D DP, AXI, 3D * KERRE=0 IF(IFOUR.NE.0.AND.IFOUR.NE.2.AND.IFOUR.NE.-1) THEN KERRE=99 RETURN ENDIF * * QUELQUES INITIALISATIONS * JEBOUC=0 IIMPI0=IIMPI 2020 JEBOUC=JEBOUC+1 * IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77881) (XMAT(I),I=1,NCOMAT) 77881 FORMAT(1X,' XMAT=',1PE12.5/) WRITE(IOIMP,77882) (SIG0(I),I=1,NSTRS) 77882 FORMAT(1X,' SIG0 ' /(6(1X,1PE12.5))/) WRITE(IOIMP,77883) (VAR0(I),I=1,NVARI) 77883 FORMAT(1X,' VAR0 ' /(6(1X,1PE12.5))/) WRITE(IOIMP,77884) (DEPST(I),I=1,NSTRS) 77884 FORMAT(1X,' DEPST ' /(6(1X,1PE12.5))/) ENDIF * * ON RECUPERE LES CONSTANTES DU MATERIAU * YOUN = XMAT(1) YOUE1= XMAT(14) XNU = XMAT(2) PI0 = XMAT(5) COHE = XMAT(6) GM = XMAT(7) BETA = XMAT(8) HUA = XMAT(9) HUB = XMAT(10) HUR = XMAT(11) P1 = XMAT(12) ZN = XMAT(13) * * TESTS SUR LES DONNEES * IF(COHE.EQ.0.D0) COHE = 1.D-10 * PI0 TRASI0 = COHE - TRACE(SIG0)/3.D0 IF(TRASI0.LT.0.D0.OR.HUB.LE.0.D0) THEN KERRE=22 RETURN ENDIF IF(ZN.LT.0.D0.OR.ZN.GE.1.D0) THEN KERRE=22 RETURN ENDIF * E2=YOUE1/(P1**ZN) E1=E2*(TRASI0**ZN) G2=E1/(1.D0+XNU) G=G2*0.5D0 UNSURE=1.D0/E1 BULK2=E2/(3.D0*(1.D0-2.D0*XNU)) * * ON CALCULE LES CONTRAINTES FINALES EN ELASTIQUE * TRACEP= TRACE(DEPST) ZN1=1.D0-ZN UNSURN=1.D0/ZN1 FAC= TRASI0 - (TRASI0**ZN1 - ZN1*BULK2*TRACEP)**UNSURN STOT(1)= SIG0(1)+G2*(DEPST(1)-TRACEP/3.D0) +FAC STOT(2)= SIG0(2)+G2*(DEPST(2)-TRACEP/3.D0) +FAC STOT(3)= SIG0(3)+G2*(DEPST(3)-TRACEP/3.D0) +FAC STOT(4)= SIG0(4)+G*DEPST(4) IF(IFOUR.EQ.2) THEN STOT(5)= SIG0(5)+G*DEPST(5) STOT(6)= SIG0(6)+G*DEPST(6) ELSE STOT(5)= 0.D0 STOT(6)= 0.D0 ENDIF * IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77875) TRACEP 77875 FORMAT(1X,' TRACEP = ',1PE12.5/) WRITE(IOIMP,77886) (STOT(I),I=1,NSTRS) 77886 FORMAT(1X,' STOT ' /(6(1X,1PE12.5))/) ENDIF * * ON CALCULE LE CRITERE * PI= (PI0+COHE)*EXP(-BETA*VAR0(2)) - COHE IF(HUA.EQ.0.D0) THEN XKSI=1.D0+HUR ELSE XKSI=VAR0(1)/(HUA+VAR0(1)) + HUR ENDIF CALL HUJCRI(STOT,PI,XKSI,XMAT,XINV,Y,SELAS) * IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77887) XINV(1),XINV(2),Y,SELAS,PI,XKSI 77887 FORMAT(1X,' XINV(1)= ',1PE12.5,2X,'XINV(2)= ',1PE12.5, . 2X,'Y= ',1PE12.5,2X,'SELAS= ',1PE12.5,2X,'PI= ',1PE12.5/ . 2X,'KSI=',1PE12.5/) ENDIF * PETI=1.1D00*PRECIS*Y ZERO=0.D0 CALL EPSPRE(Y,ZERO,PETI,ITRY) IF(ITRY.EQ.0.AND.Y.GT.ZERO) GO TO 4 * * SORTIE ELASTIQUE * VARF(1)=VAR0(1) VARF(2)=VAR0(2) DO 3 IBA=1,NSTRS SIGF(IBA)=STOT(IBA) DEFP(IBA)=0.D0 3 CONTINUE IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77003) 77003 FORMAT(1X,' SORTIE ELASTIQUE ' /) ENDIF RETURN * * ECOULEMENT PLASTIQUE * 4 CONTINUE ITER =0 DEPS =0.D0 DEPSI=0.D0 CALL SHIFTD(STOT,SIGEL,6) CALL ZDANUL(DEFPLA,6) * * ------------------------------- * | LES ITERATIONS INTERNES | * ------------------------------- 555 CONTINUE ITER=ITER+1 IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77888) ITER 77888 FORMAT(1X,' >>>>>>>>>>> HUJEUX - ITER =',I4/) ENDIF * CALL HUJFLO(IFOUR,SIGEL,XINV,PI,XKSI,XMAT,DFDS,DFDPI, . DFDKSI,DGDS,HV,HQ) IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77889) (DFDS(IJ),IJ=1,NSTRS) 77889 FORMAT(1X,' DFDS '/(6(1X,1PE12.5))/) WRITE(IOIMP,77899) (DGDS(IJ),IJ=1,NSTRS) 77899 FORMAT(1X,' DGDS '/(6(1X,1PE12.5))/) WRITE(IOIMP,77890) DFDPI,DFDKSI,HV,HQ 77890 FORMAT(1X,' DFDPI= ',1PE12.5,2X,'DFDKSI= ',1PE12.5/ . 1X,' HV=',1PE12.5,2X,'HQ=',1PE12.5/) ENDIF * * CALCUL DE D*DGDS * ABSI=ABS(COHE-XINV(1)) BULKGT=(BULK2*(ABSI**ZN) - G2/3.D0)*TRACE(DGDS) STOT(1)=G2*DGDS(1)+BULKGT STOT(2)=G2*DGDS(2)+BULKGT STOT(3)=G2*DGDS(3)+BULKGT STOT(4)=G*DGDS(4) IF(IFOUR.EQ.2) THEN STOT(5)=G*DGDS(5) STOT(6)=G*DGDS(6) ELSE STOT(5)=0.D0 STOT(6)=0.D0 ENDIF * IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77891) (STOT(IJ),IJ=1,NSTRS) 77891 FORMAT(1X,' D*DFDS '/(6(1X,1PE12.5))/) ENDIF * DFTDDG=DDOT(NSTRS,DFDS,1,STOT,1) DEPSI=Y/(DFTDDG - DFDPI*HV - DFDKSI*HQ +1.D-20) DEPS = DEPS+DEPSI IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77892) DEPSI,DEPS 77892 FORMAT(1X,' DEPSI=',1PE12.5,2X,'DEPS=',1PE12.5/) ENDIF IF(DEPS.LT.0.D0) GO TO 6 * * ON INCREMENTE * DO 66 I=1,NSTRS DEFPLA(I)=DEFPLA(I)+DGDS(I)*DEPSI STOT(I)=DEPST(I)-DEFPLA(I) 66 CONTINUE IF(IFOUR.NE.2) THEN DEFPLA(5)=0.D0 DEFPLA(6)=0.D0 STOT(5)=0.D0 STOT(6)=0.D0 ENDIF * * DEVIATEUR DE DEFPLA * DO 67 I=1,6 IF(I.LE.3) THEN DDEFPL(I)=DEFPLA(I)-TRACE(DEFPLA)/3.D0 ELSE DDEFPL(I)=DEFPLA(I) ENDIF 67 CONTINUE IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77991) (DEFPLA(IJ),IJ=1,NSTRS) 77991 FORMAT(1X,' DEFPLA '/(6(1X,1PE12.5))/) WRITE(IOIMP,73991) (DDEFPL(IJ),IJ=1,NSTRS) 73991 FORMAT(1X,' DDEFPL '/(6(1X,1PE12.5))/) ENDIF * TRACEP= TRACE(STOT) FAC = TRASI0 - (TRASI0**ZN1 - ZN1*BULK2*TRACEP)**UNSURN SIGEL(1)= SIG0(1)+G2*(STOT(1)-TRACEP/3.D0) +FAC SIGEL(2)= SIG0(2)+G2*(STOT(2)-TRACEP/3.D0) +FAC SIGEL(3)= SIG0(3)+G2*(STOT(3)-TRACEP/3.D0) +FAC SIGEL(4)= SIG0(4)+G*STOT(4) IF(IFOUR.EQ.2) THEN SIGEL(5)= SIG0(5)+G*STOT(5) SIGEL(6)= SIG0(6)+G*STOT(6) ELSE SIGEL(5)= 0.D0 SIGEL(6)= 0.D0 ENDIF * PI= (PI0+COHE)*EXP(-BETA*(VAR0(2)+TRACE(DEFPLA))) - COHE EPSE=VONEPS(DDEFPL,1,ALFAH,COVNMS)+VAR0(1) IF(HUA.EQ.0.D0) THEN XKSI=1.D0+HUR ELSE XKSI=EPSE/(HUA+EPSE) + HUR ENDIF * * TEST DE CONVERGENCE * CALL HUJCRI(SIGEL,PI,XKSI,XMAT,XINV,Y,SELAS) STST=ABS(Y/SELAS) * IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77893) (SIGEL(IJ),IJ=1,NSTRS) 77893 FORMAT(1X,' NEW SIG '/(6(1X,1PE12.5))/) WRITE(IOIMP,77894) PI,Y,SELAS,STST 77894 FORMAT(1X,' PI= ',1PE12.5,2X,'Y= ',1PE12.5,2X, . 'SELAS= ',1PE12.5,2X,'STST= ',1PE12.5/) ENDIF * IF(STST.LE.PRECIS) GO TO 8 IF(ITER.GT.ITMAX) GO TO 7 GO TO 555 * 6 CONTINUE KERRE=1 GO TO 9 7 CONTINUE KERRE=2 9 CONTINUE IF (JEBOUC.EQ.1) THEN IIMPI=48 GO TO 2020 ELSE IIMPI=IIMPI0 ENDIF * 8 CONTINUE * * REMPLISSAGE * VARF(1)=EPSE VARF(2)=VAR0(2)+TRACE(DEFPLA) DO 12 IBA=1,NSTRS SIGF(IBA)=SIGEL(IBA) DEFP(IBA)=DEFPLA(IBA) 12 CONTINUE * IF(IIMPI.EQ.48) THEN WRITE(IOIMP,77895) (STOT(IJ),IJ=1,NSTRS) 77895 FORMAT(1X,' DEF_ELA '/(6(1X,1PE12.5))/) WRITE(IOIMP,77896) (DEFP(IJ),IJ=1,NSTRS) 77896 FORMAT(1X,' DEF_PLA '/(6(1X,1PE12.5))/) WRITE(IOIMP,77897) VARF(1),VARF(2) 77897 FORMAT(1X,' VARF(1)= ',1PE12.5,2X,'VARF(2)= ',1PE12.5/) ENDIF * 1000 RETURN END