C DP_SOL SOURCE CB215821 17/11/30 21:16:00 9639 C DP_SOIL SOURCE AF221230 12/10/04 21:15:18 7520 SUBROUTINE DP_SOL (XMAT,VAR0,VARF,SIG0,SIGF,DEPST,XCAR) C C====&===1=========2=========3=========4=========5=========6=========7== C Commentaires : Loi Drucker Prager avec ecrouissage non lineaire et C loi non associée C Traits : - C C Auteur : A. Frau (Dr - Ing.) - CEA/DEN/DANS/DM2S/SEMT/EMSI C====&===1=========2=========3=========4=========5=========6=========7== C C----DECLARATION GENERALES---------------------------------------------- C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO C C----DECLARATION SEGMENTS----------------------------------------------- C REAL*8 XMAT(*),SIG0(*),SIGF(*),VAR0(*),VARF(*),DEPST(*),XCAR(*) REAL*8 XNU,XYOUN,XSY REAL*8 XMU,XLAM,XALPHA,XSIGI,XBET1,XGAM1,XDEL1 REAL*8 XDEEPP(6),XEPP0(6) REAL*8 XSIGT(6),XCHIT,XQT,XFHIT,XDSIGT(6),XNDSIGT REAL*8 XNUT(6),XPIT(6) REAL*8 XTOL1,XVPL1 REAL*8 XGAMN,XSIGN(6),XEPPN(6),XCHIN,XQN,XFHIN,XDSIGN(6),XNDSIGN REAL*8 XGAMN1,XSIGN1(6),XEPPN1(6),XCHIN1,XQN1,XFHIN1 REAL*8 XDFDG,XDEGAM,XDSIGN1(6),XNDSIGN1,XTSIGN1 REAL*8 XGAMF,XCHIF,XQF,XEPPF(6),XVALT1,XVALT2 C C----CONTROLE OPTION CALCUL -------------------------------------------- IF ((IFOUR.NE.2).AND.(IFOUR.NE.0).AND.(IFOUR.NE.-1)) THEN WRITE(IOIMP,10000) 10000 FORMAT('Option calcul mauvais') WRITE(IOIMP,10001) 10001 FORMAT(' seulement:') WRITE(IOIMP,10002) 10002 FORMAT(' 3D MODE PLAN') WRITE(IOIMP,10003) 10003 FORMAT(' 2D MODE AXIS') WRITE(IOIMP,10004) 10004 FORMAT(' 2D MODE PLAN') STOP ENDIF C C----MISE EN DONNEES---------------------------------------------------- C C Module de young XYOUN = XMAT(1) C Coef Poisson XNU = XMAT(2) C alpha XALPHA = XMAT(5) C Sigma_y XSY = XMAT(6) C Sig_inf XSIGI = XMAT(7) C Beta XBET1 = XMAT(8) C Gamma XGAM1 = XMAT(9) C Delta XDEL1 = XMAT(10) C C----VARIABLES INTERNES------------------------------------------------- c Coef mu elastique XMU = (0.5D0)*((XYOUN)/(1.0D0 + XNU)) c Coef lambda elastique XLAM = ((XNU)/(1.0D0-((2.0D0)*(XNU)))) C C----ACTUALISATION DEFORMATION------------------------------------------ C C C----PREDICTION ELASTIQUE----------------------------------------------- C c gamma_test = 0 donc le delta_plastique = 0 c Initialissation Eps_pla_initiales DO 11 I1=1,6 XDEEPP(I1) = 0.0D0 XEPP0(I1) = VAR0(I1) 11 CONTINUE c c Calcul de la variable q_test = q_ini -> gamma_test=0 XCHIT = VAR0(7) c c Calcul de la variable q_test = q_ini -> gamma_test=0 XQT = ((-1.0D0)*(XSIGI - XSY))*(1.0D0 & - EXP(((-1.0D0)*(XBET1))*(XCHIT))) c c Calcul Sigma_test a partir du step initiale DO 12 I1=1,6 IF (I1.LE.3) THEN XSIGT(I1) = SIG0(I1) + & ((2.0D0)*(XMU))*(DEPST(I1) + ((XLAM)*(TRACE(DEPST)))) ELSE IF (IFOUR.EQ.2) THEN XSIGT(I1) = SIG0(I1) + ((1.0D0)*(XMU))*(DEPST(I1)) ELSE IF (I1.EQ.4) THEN XSIGT(I1) = SIG0(I1) + ((1.0D0)*(XMU))*(DEPST(I1)) ELSE XSIGT(I1) = 0.0D0 ENDIF ENDIF ENDIF 12 CONTINUE C C Deviatoire de Sigma_test DO 13 I1=1,6 IF (I1.LE.3) THEN XDSIGT(I1) = XSIGT(I1) - ((TRACE(XSIGT)/(3.0D0))) ELSE IF (IFOUR.EQ.2) THEN XDSIGT(I1) = XSIGT(I1) ELSE IF (I1.EQ.4) THEN XDSIGT(I1) = XSIGT(I1) ELSE XDSIGT(I1) = 0.0D0 ENDIF ENDIF ENDIF 13 CONTINUE c c Calcul de la norme partie deviatoire XNDSIGT = 0.0D0 DO 14 I1=1,6 IF (I1.LE.3) THEN XNDSIGT = XNDSIGT + (XDSIGT(I1)*XDSIGT(I1)) ELSE XNDSIGT = XNDSIGT + ((2.0D0)*(XDSIGT(I1)*XDSIGT(I1))) ENDIF 14 CONTINUE XNDSIGT = ((XNDSIGT)**(0.5D0)) C c Calcul du tensuer pi XPIT(1) = 1.0D0 XPIT(2) = 1.0D0 XPIT(3) = 1.0D0 XPIT(4) = 0.0D0 XPIT(5) = 0.0D0 XPIT(6) = 0.0D0 XPIT = (((3.0D0)**(0.5D0))/(3.0D0))*(XPIT) C c Calcul ((3**0.5)*(trace(sigma_test))) XTSIGT = 0.0D0 DO 23 I1=1,6 XTSIGT = XTSIGT + (XPIT(I1)*(XSIGT(I1))) 23 CONTINUE c c Calcul du vecteur nu_test IF ((ABS(XNDSIGT/XTSIGT)).LT.(1.0E-8)) THEN XNUT(1) = (1.0D0)/(3.0D0) XNUT(2) = (1.0D0)/(3.0D0) XNUT(3) = (-2.0D0)/(3.0D0) XNUT(4) = 0.0D0 XNUT(5) = 0.0D0 XNUT(6) = 0.0D0 XVPL1 = 0.0D0 ELSE XNUT(1) = XDSIGT(1)/XNDSIGT XNUT(2) = XDSIGT(2)/XNDSIGT XNUT(3) = XDSIGT(3)/XNDSIGT XNUT(4) = XDSIGT(4)/XNDSIGT XNUT(5) = XDSIGT(5)/XNDSIGT XNUT(6) = XDSIGT(6)/XNDSIGT XVPL1 = 1.0D0 ENDIF c Calcul du critere fhi_test XFHIT = XNDSIGT + (((2.0D0)**(0.5))*(XTSIGT))*(XALPHA) & - (((2.0D0/3.0D0)**(0.5D0))*(XSY - XQT)) C C-----TEST--------------------------------------------- C IF (XFHIT.LT.0) THEN GOTO 101 c fhi_test < 0.0D0 pas de plasticité ELSE IF (XFHIT.GE.0) THEN GOTO 102 c fhi_test >= 0.0D0 on est en plasticité END IF C C-----CALCUL ELASTIQUE--------------------------------------------- C 101 CONTINUE c c si fhi_test < 0 gamma=0 DO 15 I1=1,6 SIGF(I1) = XSIGT(I1) 15 CONTINUE C c pas des variations des variables internes DO 16 I1=1,8 IF (I1.EQ.8) THEN VARF(I1) = 0.0D0 ELSE VARF(I1) = VAR0(I1) ENDIF 16 CONTINUE RETURN C C-----CALCUL PLASTIQUE--------------------------------------------- C 102 CONTINUE c c si fhi_test >= 0 determination de gamma c Nombre max interaction NINT1 = 10000 C Tollerance XTOL1 = 1.E-8 XFHIN = XFHIT XQN = XQT XGAMN = 0.0D0 XCHIN = XCHIT c DO 17 I1=1,6 XEPPN(I1) = XEPP0(I1) XSIGN(I1) = XSIGT(I1) 17 CONTINUE c DO 110 II1=1,NINT1 IF (II1.EQ.NINT1) THEN c PRINT*,'!!!!! ATTENTION !!!!' c PRINT*,'Nombre des iteractions' c PRINT*,' maximales depasse' c STOP ENDIF XDFDG = (((2.0D0)*(XMU*XVPL1))*(XDEL1)) & + ((2.0D0)/(3.0D0))*(((XSIGI - XSY)*( & EXP(((-1.0D0)*(XBET1))*(XCHIN))))*(XBET1)) + & ((4.0D0)*((XALPHA)*(XGAM1)))* & (1.0D0 + ((3.0D0)*(XLAM)))*(XMU) XDEGAM = XFHIN/XDFDG XGAMN1 = XGAMN + XDEGAM XCHIN1 = XCHIN + (((2.0D0/3.0D0)**(0.5D0))*(XDEGAM)) XQN1 = ((-1.0D0)*(XSIGI - XSY))*(1.0D0 & - EXP(((-1.0D0)*(XBET1))*(XCHIN1))) IF ((XGAMN1).LT.(0.0D0)) THEN WRITE(IOIMP,10005) 10005 FORMAT('Error....multiplicateur negative') STOP ENDIF DO 18 I1=1,6 IF (I1.LE.3) THEN XDEEPP(I1) = ((((XDEL1*XVPL1))*(XNUT(I1)))+ & (((XGAM1))*((2.0D0)**(0.5D0)))*(XPIT(I1)))*XDEGAM ELSE IF (IFOUR.EQ.2) THEN XDEEPP(I1) = ((((XDEL1*XVPL1))*(XNUT(I1)))+ & (((XGAM1))*((2.0D0)**(0.5D0)))*(XPIT(I1)))*XDEGAM ELSE IF (I1.EQ.4) THEN XDEEPP(I1) = ((((XDEL1*XVPL1))*(XNUT(I1)))+ & (((XGAM1))*((2.0D0)**(0.5D0)))*(XPIT(I1)))*XDEGAM ELSE XDEEPP(I1) = 0.0D0 ENDIF ENDIF ENDIF IF (I1.LE.3) THEN XEPPN1(I1) = XEPPN(I1) + XDEEPP(I1) ELSE IF (IFOUR.EQ.2) THEN XEPPN1(I1) = XEPPN(I1) + XDEEPP(I1) ELSE IF (I1.EQ.4) THEN XEPPN1(I1) = XEPPN(I1) + XDEEPP(I1) ELSE XEPPN1(I1) = 0.0D0 ENDIF ENDIF ENDIF 18 CONTINUE DO 21 I1=1,6 IF (I1.LE.3) THEN XSIGN1(I1) = XSIGN(I1) - & ((2.0D0)*(XMU))*(XDEEPP(I1) & + ((XLAM)*(TRACE(XDEEPP)))) ELSE IF (IFOUR.EQ.2) THEN XSIGN1(I1) = XSIGN(I1) - & ((1.0D0)*(XMU))*(XDEEPP(I1)) ELSE IF (I1.EQ.4) THEN XSIGN1(I1) = XSIGN(I1) - & ((1.0D0)*(XMU))*(XDEEPP(I1)) ELSE XSIGN1(I1) = 0.0D0 ENDIF ENDIF ENDIF 21 CONTINUE DO 19 I1=1,6 IF (I1.LE.3) THEN XDSIGN1(I1) = XSIGN1(I1) & - ((TRACE(XSIGN1)/(3.0D0))) ELSE IF (IFOUR.EQ.2) THEN XDSIGN1(I1) = XSIGN1(I1) ELSE IF (I1.EQ.4) THEN XDSIGN1(I1) = XSIGN1(I1) ELSE XDSIGN1(I1) = 0.0D0 ENDIF ENDIF ENDIF 19 CONTINUE XNDSIGN1 = 0.0D0 DO 20 I1=1,6 IF (I1.LE.3) THEN XNDSIGN1 = XNDSIGN1 & + (XDSIGN1(I1)*XDSIGN1(I1)) ELSE XNDSIGN1 = XNDSIGN1 & + ((2.0D0)*(XDSIGN1(I1)*XDSIGN1(I1))) ENDIF 20 CONTINUE XNDSIGN1 = ((XNDSIGN1)**(0.5D0)) XTSIGN1 = XSIGN1(1) + XSIGN1(2) + XSIGN1(3) XTSIGN1 = ((XTSIGN1)/((3.0D0)**(0.5D0))) XFHIN1 = XNDSIGN1 + & (((2.0D0)**(0.5))*(XTSIGN1))*(XALPHA) & - (((2.0D0/3.0D0)**(0.5D0))*(XSY - XQN1)) c IF (ABS((XFHIN1)/(XFHIT)).LT.(XTOL1)) THEN GOTO 103 ELSE XFHIN = XFHIN1 XGAMN = XGAMN1 XCHIN = XCHIN1 XQN = XQN1 XEPPN(1) = XEPPN1(1) XEPPN(2) = XEPPN1(2) XEPPN(3) = XEPPN1(3) XEPPN(4) = XEPPN1(4) XEPPN(5) = XEPPN1(5) XEPPN(6) = XEPPN1(6) XSIGN(1) = XSIGN1(1) XSIGN(2) = XSIGN1(2) XSIGN(3) = XSIGN1(3) XSIGN(4) = XSIGN1(4) XSIGN(5) = XSIGN1(5) XSIGN(6) = XSIGN1(6) ENDIF 110 CONTINUE 103 CONTINUE XCHIF = XCHIN1 XGAMF = XGAMN1 SIGF(1) = XSIGN1(1) SIGF(2) = XSIGN1(2) SIGF(3) = XSIGN1(3) SIGF(4) = XSIGN1(4) SIGF(5) = XSIGN1(5) SIGF(6) = XSIGN1(6) VARF(1) = XEPPN1(1) VARF(2) = XEPPN1(2) VARF(3) = XEPPN1(3) VARF(4) = XEPPN1(4) VARF(5) = XEPPN1(5) VARF(6) = XEPPN1(6) VARF(7) = XCHIF VARF(8) = XGAMF RETURN END