C OUGLOV SOURCE BR232186 16/09/12 12:43:41 9078 SUBROUTINE OUGLOV(XMAT,VAR0,VARF,SIG0,SIGF,DEPST,IFOUR) C C====&===1=========2=========3=========4=========5=========6=========7== C Commentaires : Subroutine permettant de mettre en oeuvre le C modele OUGLOVA pour representer C le comportement de l'acier corrodé C C Auteurs : R. Paredes C : CEA-DEN/DANS/DM2S/SEMT/EMSI C : Romili.Paredes@cea.fr C C Date : Aout 2016 C====&===1=========2=========3=========4=========5=========6=========7== C C-----DECLARATION GENERALE---------------------------------------------- C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C C-----DECLARATION DES VARIABLES----------------------------------------- C DIMENSION XMAT(*),DEPST(*),VAR0(*),SIGF(*),VARF(*) REAL*8 UNIT3(3,3), DFEDPHI(3) REAL*8 EPSR(3), XEPSP(3,3), XEPSF(3,3), XEPSE(3,3), SIG(3,3) REAL*8 DFDSG(3,3), DPHDSG(3,3), DEPSP(3,3), DVSIG(3,3), EDPH(3,3) C-----Parametre pour le nombre d iteration locales internes IMAX = 1000 C C-----DEFINITION DE LA MATRICE UNITE D ORDRE 3-------------------------- C DO I = 1,3 DO J = 1,3 IF (I.EQ.J) THEN UNIT3(I,J) = 1.0D0 ELSE UNIT3(I,J) = 0.0D0 ENDIF ENDDO ENDDO C C-----DONNEES MATERIAUX------------------------------------------------- C XYG = XMAT(1) XNU = XMAT(2) XSIGY = XMAT(5) XK = XMAT(6) XM = XMAT(7) XTC = XMAT(8) XDC = XMAT(9) XLA = (XNU*XYG)/((1.0D0 - 2.0D0*XNU)*(1.0D0 + XNU)) XMU = XYG/(2.0D0*(1.0D0 + XNU)) C C-----CORROSION--------------------------------------------------------- C C-----Deformation plastique a rupture IF (XTC.LE.15.0D0) THEN EPSR(1) = -0.0111D0*XTC + 0.2345D0 ELSE EPSR(1) = -0.0006D0*XTC + 0.051D0 ENDIF C-----Coefficient de contraction XNUSTAR = 0.5D0 - (XSIGY/(XYG*EPSR(1)))*(0.5D0 - XNU) EPSR(2) = -EPSR(1)*XNUSTAR EPSR(3) = -EPSR(1)*XNUSTAR C-----Seuil de rupture et d'endommagement TEPSR = EPSR(1)**2.0D0 + EPSR(2)**2.0D0 + EPSR(3)**2.0D0 XPR = ((2.0D0/3.0D0)*TEPSR)**(0.5D0) XPD = 0.8D0*XPR C C-----ETAT INITIAL------------------------------------------------------ C XD = VAR0(1) XR = VAR0(2) XP = VAR0(3) XZT = VAR0(4) XNRUP = VAR0(5) C----------------------------------------------------------------------- C-----CALCUL TRIDIMENSIONNEL-------------------------------------------- C----------------------------------------------------------------------- IF (IFOUR.EQ.2) THEN C--------Deformation plastique XEPSP(1,1) = VAR0(6) XEPSP(2,2) = VAR0(7) XEPSP(3,3) = VAR0(8) XEPSP(1,2) = VAR0(9)/2.0D0 XEPSP(1,3) = VAR0(10)/2.0D0 XEPSP(2,3) = VAR0(11)/2.0D0 XEPSP(2,1) = XEPSP(1,2) XEPSP(2,3) = XEPSP(2,3) XEPSP(3,1) = XEPSP(1,3) XEPSP(3,2) = XEPSP(2,3) C--------Deformation totale XEPSF(1,1) = DEPST(1) + VAR0(12) XEPSF(2,2) = DEPST(2) + VAR0(13) XEPSF(3,3) = DEPST(3) + VAR0(14) XEPSF(1,2) = (DEPST(4) + VAR0(15))/2.0D0 XEPSF(1,3) = (DEPST(5) + VAR0(16))/2.0D0 XEPSF(2,3) = (DEPST(6) + VAR0(17))/2.0D0 XEPSF(2,1) = XEPSF(1,2) XEPSF(2,3) = XEPSF(2,3) XEPSF(3,1) = XEPSF(1,3) XEPSF(3,2) = XEPSF(2,3) C--------Deformation elastique DO I = 1,3 DO J = 1,3 XEPSE(I,J) = XEPSF(I,J) - XEPSP(I,J) ENDDO ENDDO C--------Prediction elastique TREPSE = XEPSE(1,1) + XEPSE(2,2) + XEPSE(3,3) DO I = 1,3 DO J = 1,3 AUX1 = XLA*TREPSE*UNIT3(I,J) AUX2 = 2.0D0*XMU*XEPSE(I,J) SIG(I,J) = (1.0D0 - XD)*(AUX1 + AUX2) ENDDO ENDDO CALL J2(SIG,SIGEQ) CALL DEVIAT(SIG,DVSIG) C--------Boucle d'endommagement DO I = 1,IMAX C-----------Evaluation de la fonction seuil FP = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY) FP0 = FP C-----------Boucle de plasticite IF (FP.GT.1.0D0) THEN DO J = 1,IMAX C-----------------Calcul des derives DO K = 1,3 DO L = 1,3 AUX1 = ((3.0D0/2.0D0)/(SIGEQ*(1.0D0 - XD))) DFDSG(K,L) = AUX1*DVSIG(K,L) ENDDO ENDDO DO K = 1,3 DO L = 1,3 DPHDSG(K,L) = DFDSG(K,L) ENDDO ENDDO C-----------------Calcul du multiplicateur plastique TRDPHDSG = DPHDSG(1,1) + DPHDSG(2,2) + DPHDSG(3,3) DO K = 1,3 DO L = 1,3 AUX1 = XLA*TRDPHDSG*UNIT3(K,L) AUX2 = 2.0D0*XMU*DPHDSG(K,L) EDPH(K,L) = (AUX1 + AUX2) ENDDO ENDDO DO K = 1,3 DFEDPHI(K) = 0.0D0 DO L = 1,3 AUX1 = DFDSG(K,L)*EDPH(L,K) DFEDPHI(K) = DFEDPHI(K) + AUX1 ENDDO ENDDO AUX4 = DFEDPHI(1) + DFEDPHI(2) + DFEDPHI(3) AUX5 = XK/XM AUX6 = SIGEQ/(XK*(1.0D0 - XD)) AUX7 = XSIGY/XK AUX8 = (1.0D0 - XM) XLAMBDA = FP/(AUX4 + AUX5*((AUX6 + AUX7)**AUX8)) C-----------------Mise a jour des variables DO K = 1,3 DO L = 1,3 SIG(K,L) = SIG(K,L) - XLAMBDA*EDPH(K,L) ENDDO ENDDO CALL J2(SIG,SIGEQ) CALL DEVIAT(SIG,DVSIG) DO K = 1,3 DO L = 1,3 XEPSP(K,L) = XEPSP(K,L) + XLAMBDA*DPHDSG(K,L) ENDDO ENDDO XP = XP + XLAMBDA/(1.0D0 - XD); XR = XK*(XP**(1.0D0/XM)) C-----------------Reevaluation de la fonction seuil FP = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY) CRIT1 = DABS(FP/FP0) CRIT2 = DABS(XLAMBDA) IF ((CRIT1.LE.1.0D-8).OR.(CRIT2.LE.1.0D-10)) THEN GOTO 910 ENDIF ENDDO C-----Sortie de la boucle de plasticité 910 CONTINUE C--------------Evaluation du seuil d'endommagement FD = XP - (XPD + XZT) IF (I.EQ.1) THEN FD0 = FD ENDIF CRIT = DABS(FD/FD0) IF ((FD.GT.0.0D0).AND.(CRIT.GT.1.0D-5)) THEN C-----------------Calcul de l'endommagement XD = (XDC/(XPR - XPD))*(XP - XPD) XZT = XP - XPD C-----------------Endommagement critique IF (XD.GT.XDC) THEN XD = XDC ENDIF ELSE GOTO 920 ENDIF IF (CRIT.LE.1.0D-5) THEN GOTO 920 ENDIF ELSE GOTO 920 ENDIF ENDDO C-----Sortie de la boucle d'endommagement 920 CONTINUE C--------Rupture IF (XP.GE.XPR) THEN XNRUP = 1 XP = XPR DO I = 1,3 DO J = 1,3 SIG(I, J) = 1.0D0 ENDDO ENDDO ENDIF C C-----SOCKAGE EN SORTIE------------------------------------------------- C C--------Les variables internes VARF(1) = XD VARF(2) = XR VARF(3) = XP VARF(4) = XZT VARF(5) = XNRUP VARF(6) = XEPSP(1,1) VARF(7) = XEPSP(2,2) VARF(8) = XEPSP(3,3) VARF(9) = XEPSP(1,2)*2.0D0 VARF(10) = XEPSP(1,3)*2.0D0 VARF(11) = XEPSP(2,3)*2.0D0 VARF(12) = XEPSF(1,1) VARF(13) = XEPSF(2,2) VARF(14) = XEPSF(3,3) VARF(15) = XEPSF(1,2)*2.0D0 VARF(16) = XEPSF(1,3)*2.0D0 VARF(17) = XEPSF(2,3)*2.0D0 C--------Les contraintes SIGF(1) = SIG(1,1) SIGF(2) = SIG(2,2) SIGF(3) = SIG(3,3) SIGF(4) = SIG(1,2) SIGF(5) = SIG(1,3) SIGF(6) = SIG(2,3) C----------------------------------------------------------------------- C-----CALCUL DEFORMATION PLANE------------------------------------------ C ----------------------------------------------------------------------- ELSEIF (IFOUR.EQ.-1) THEN C--------Deformation plastique XEPSP(1,1) = VAR0(6) XEPSP(2,2) = VAR0(7) XEPSP(3,3) = VAR0(8) XEPSP(1,2) = VAR0(9)/2.0D0 XEPSP(1,3) = 0.0D0 XEPSP(2,3) = 0.0D0 XEPSP(2,1) = XEPSP(1,2) XEPSP(2,3) = XEPSP(2,3) XEPSP(3,1) = XEPSP(1,3) XEPSP(3,2) = XEPSP(2,3) C--------Deformation totale XEPSF(1,1) = DEPST(1) + VAR0(10) XEPSF(2,2) = DEPST(2) + VAR0(11) XEPSF(3,3) = DEPST(3) + VAR0(12) XEPSF(1,2) = (DEPST(4) + VAR0(13))/2.0D0 XEPSF(1,3) = 0.0D0 XEPSF(2,3) = 0.0D0 XEPSF(2,1) = XEPSF(1,2) XEPSF(2,3) = XEPSF(2,3) XEPSF(3,1) = XEPSF(1,3) XEPSF(3,2) = XEPSF(2,3) C--------Deformation elastique DO I = 1,3 DO J = 1,3 XEPSE(I,J) = XEPSF(I,J) - XEPSP(I,J) ENDDO ENDDO C--------Deformation plane XEPSE(1,3) = 0.0D0 XEPSE(2,3) = 0.0D0 XEPSE(3,1) = 0.0D0 XEPSE(3,2) = 0.0D0 XEPSE(3,3) = 0.0D0 C--------Prediction elastique TREPSE = XEPSE(1,1) + XEPSE(2,2) + XEPSE(3,3) DO I = 1,3 DO J = 1,3 AUX1 = XLA*TREPSE*UNIT3(I,J) AUX2 = 2.0D0*XMU*XEPSE(I,J) SIG(I,J) = (1.0D0 - XD)*(AUX1 + AUX2) ENDDO ENDDO CALL J2(SIG,SIGEQ) CALL DEVIAT(SIG,DVSIG) C--------Boucle d'endommagement DO I = 1,IMAX C-----------Evaluation de la fonction seuil FP = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY) FP0 = FP C-----------Boucle de plasticite IF (FP.GT.1.0D0) THEN DO J = 1,IMAX C-----------------Calcul des derives DO K = 1,3 DO L = 1,3 AUX1 = ((3.0D0/2.0D0)/(SIGEQ*(1.0D0 - XD))) DFDSG(K,L) = AUX1*DVSIG(K,L) ENDDO ENDDO DO K = 1,3 DO L = 1,3 DPHDSG(K,L) = DFDSG(K,L) ENDDO ENDDO C-----------------Calcul du multiplicateur plastique TRDPHDSG = DPHDSG(1,1) + DPHDSG(2,2) + DPHDSG(3,3) DO K = 1,3 DO L = 1,3 AUX1 = XLA*TRDPHDSG*UNIT3(K,L) AUX2 = 2.0D0*XMU*DPHDSG(K,L) EDPH(K,L) = (AUX1 + AUX2) ENDDO ENDDO DO K = 1,3 DFEDPHI(K) = 0.0D0 DO L = 1,3 AUX1 = DFDSG(K,L)*EDPH(L,K) DFEDPHI(K) = DFEDPHI(K) + AUX1 ENDDO ENDDO AUX4 = DFEDPHI(1) + DFEDPHI(2) + DFEDPHI(3) AUX5 = XK/XM AUX6 = SIGEQ/(XK*(1.0D0 - XD)) AUX7 = XSIGY/XK AUX8 = (1.0D0 - XM) XLAMBDA = FP/(AUX4 + AUX5*((AUX6 + AUX7)**AUX8)) C-----------------Mise a jour des variables DO K = 1,3 DO L = 1,3 SIG(K,L) = SIG(K,L) - XLAMBDA*EDPH(K,L) ENDDO ENDDO CALL J2(SIG,SIGEQ) CALL DEVIAT(SIG,DVSIG) DO K = 1,3 DO L = 1,3 XEPSP(K,L) = XEPSP(K,L) + XLAMBDA*DPHDSG(K,L) ENDDO ENDDO XP = XP + XLAMBDA/(1.0D0 - XD); XR = XK*(XP**(1.0D0/XM)) C-----------------Reevaluation de la fonction seuil FP = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY) CRIT1 = DABS(FP/FP0) CRIT2 = DABS(XLAMBDA) IF ((CRIT1.LE.1.0D-10).OR.(CRIT2.LE.1.0D-10)) THEN GOTO 930 ENDIF ENDDO C-----Sortie de la boucle de plasticité 930 CONTINUE C--------------Evaluation du seuil d'endommagement FD = XP - (XPD + XZT) IF (I.EQ.1) THEN FD0 = FD ENDIF CRIT = DABS(FD/FD0) IF ((FD.GT.0.0D0).AND.(CRIT.GT.1.0D-5)) THEN C-----------------Calcul de l'endommagement XD = (XDC/(XPR - XPD))*(XP - XPD) XZT = XP - XPD C-----------------Endommagement critique IF (XD.GT.XDC) THEN XD = XDC ENDIF ELSE GOTO 940 ENDIF IF (CRIT.LE.1.0D-5) THEN GOTO 940 ENDIF ELSE GOTO 940 ENDIF ENDDO C-----Sortie de la boucle d'endommagement 940 CONTINUE C-----Rupture IF (XP.GE.XPR) THEN XNRUP = 1 XP = XPR DO I = 1,3 DO J = 1,3 SIG(I, J) = 1.0D0 ENDDO ENDDO ENDIF C C-----SOCKAGE EN SORTIE------------------------------------------------- C C--------Les variables internes VARF(1) = XD VARF(2) = XR VARF(3) = XP VARF(4) = XZT VARF(5) = XNRUP VARF(6) = XEPSP(1,1) VARF(7) = XEPSP(2,2) VARF(8) = XEPSP(3,3) VARF(9) = XEPSP(1,2)*2.0D0 VARF(10) = XEPSF(1,1) VARF(11) = XEPSF(2,2) VARF(12) = XEPSF(3,3) VARF(13) = XEPSF(1,2)*2.0D0 C--------Les contraintes SIGF(1) = SIG(1,1) SIGF(2) = SIG(2,2) SIGF(3) = SIG(3,3) SIGF(4) = SIG(1,2) C----------------------------------------------------------------------- C-----CALCUL CONTRAINTES PLANES----------------------------------------- C----------------------------------------------------------------------- ELSEIF (IFOUR.EQ.-2) THEN C--------Deformation plastique XEPSP(1,1) = VAR0(6) XEPSP(2,2) = VAR0(7) XEPSP(3,3) = VAR0(8) XEPSP(1,2) = VAR0(9)/2.0D0 XEPSP(1,3) = 0.0D0 XEPSP(2,3) = 0.0D0 XEPSP(2,1) = XEPSP(1,2) XEPSP(2,3) = XEPSP(2,3) XEPSP(3,1) = XEPSP(1,3) XEPSP(3,2) = XEPSP(2,3) C--------Deformation totale XEPSF(1,1) = DEPST(1) + VAR0(10) XEPSF(2,2) = DEPST(2) + VAR0(11) XEPSF(3,3) = DEPST(3) + VAR0(12) XEPSF(1,2) = (DEPST(4) + VAR0(13))/2.0D0 XEPSF(1,3) = 0.0D0 XEPSF(2,3) = 0.0D0 XEPSF(2,1) = XEPSF(1,2) XEPSF(2,3) = XEPSF(2,3) XEPSF(3,1) = XEPSF(1,3) XEPSF(3,2) = XEPSF(2,3) C-------Modification de la déformation 33 XEPSF(3,3) = VAR0(14) C--------Boucle contrainte plane DO M = 1,IMAX C-----------Deformation elastique DO I = 1,3 DO J = 1,3 XEPSE(I,J) = XEPSF(I,J) - XEPSP(I,J) ENDDO ENDDO C-----------Prediction elastique TREPSE = XEPSE(1,1) + XEPSE(2,2) + XEPSE(3,3) DO I = 1,3 DO J = 1,3 AUX1 = XLA*TREPSE*UNIT3(I,J) AUX2 = 2.0D0*XMU*XEPSE(I,J) SIG(I,J) = (1.0D0 - XD)*(AUX1 + AUX2) ENDDO ENDDO CALL J2(SIG,SIGEQ) CALL DEVIAT(SIG,DVSIG) C-----------Boucle d'endommagement DO I = 1,IMAX C--------------Evaluation de la fonction seuil FP = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY) FP0 = FP C--------------Boucle de plasticite IF (FP.GT.1.0D0) THEN DO J = 1,IMAX C--------------------Calcul des derives DO K = 1,3 DO L = 1,3 AUX1 = ((3.0D0/2.0D0)/(SIGEQ*(1.0D0 - XD))) DFDSG(K,L) = AUX1*DVSIG(K,L) ENDDO ENDDO DO K = 1,3 DO L = 1,3 DPHDSG(K,L) = DFDSG(K,L) ENDDO ENDDO C--------------------Calcul du multiplicateur plastique TRDPHDSG = DPHDSG(1,1) + DPHDSG(2,2) + DPHDSG(3,3) DO K = 1,3 DO L = 1,3 AUX1 = XLA*TRDPHDSG*UNIT3(K,L) AUX2 = 2.0D0*XMU*DPHDSG(K,L) EDPH(K,L) = (AUX1 + AUX2) ENDDO ENDDO DO K = 1,3 DFEDPHI(K) = 0.0D0 DO L = 1,3 AUX1 = DFDSG(K,L)*EDPH(L,K) DFEDPHI(K) = DFEDPHI(K) + AUX1 ENDDO ENDDO AUX4 = DFEDPHI(1) + DFEDPHI(2) + DFEDPHI(3) AUX5 = XK/XM AUX6 = SIGEQ/(XK*(1.0D0 - XD)) AUX7 = XSIGY/XK AUX8 = (1.0D0 - XM) XLAMBDA = FP/(AUX4 + AUX5*((AUX6 + AUX7)**AUX8)) C--------------------Mise a jour des variables DO K = 1,3 DO L = 1,3 SIG(K,L) = SIG(K,L) - XLAMBDA*EDPH(K,L) ENDDO ENDDO CALL J2(SIG,SIGEQ) CALL DEVIAT(SIG,DVSIG) DO K = 1,3 DO L = 1,3 XEPSP(K,L) = XEPSP(K,L) + XLAMBDA*DPHDSG(K,L) ENDDO ENDDO XP = XP + XLAMBDA/(1.0D0 - XD); XR = XK*(XP**(1.0D0/XM)) C--------------------Reevaluation de la fonction seuil FP = (SIGEQ/(1.0D0 - XD)) - (XR + XSIGY) CRIT1 = DABS(FP/FP0) CRIT2 = DABS(XLAMBDA) IF ((CRIT1.LE.1.0D-10).OR.(CRIT2.LE.1.0D-10)) THEN GOTO 950 ENDIF ENDDO C-----Sortie de la boucle de plasticité 950 CONTINUE C-----------------Evaluation du seuil d'endommagement FD = XP - (XPD + XZT) IF (I.EQ.1) THEN FD0 = FD ENDIF CRIT = DABS(FD/FD0) IF ((FD.GT.0.0D0).AND.(CRIT.GT.1.0D-5)) THEN C--------------------Calcul de l'endommagement XD = (XDC/(XPR - XPD))*(XP - XPD) XZT = XP - XPD C--------------------Endommagement critique IF (XD.GT.XDC) THEN XD = XDC ENDIF ELSE GOTO 960 ENDIF IF (CRIT.LE.1.0D-5) THEN GOTO 960 ENDIF ELSE GOTO 960 ENDIF ENDDO C-----Sortie de la boucle d'endommagement 960 CONTINUE C-----------Contraintes Planes CRIT = SIG(3,3)/SIG(2,2) IF (CRIT.GT.1.0D-8) THEN XEPSF(3,3) = XEPSF(3,3) - SIG(3,3)/XYG ELSE GOTO 970 ENDIF ENDDO C-----Sortie de la boucle de contraintes planes 970 CONTINUE C-----Rupture IF (XP.GE.XPR) THEN XNRUP = 1 XP = XPR DO I = 1,3 DO J = 1,3 SIG(I, J) = 1.0D0 ENDDO ENDDO ENDIF C C-----SOCKAGE EN SORTIE------------------------------------------------- C C--------Les variables internes VARF(1) = XD VARF(2) = XR VARF(3) = XP VARF(4) = XZT VARF(5) = XNRUP VARF(6) = XEPSP(1,1) VARF(7) = XEPSP(2,2) VARF(8) = XEPSP(3,3) VARF(9) = XEPSP(1,2)*2.0D0 VARF(10) = XEPSF(1,1) VARF(11) = XEPSF(2,2) VARF(12) = XEPSF(3,3) VARF(13) = XEPSF(1,2)*2.0D0 VARF(14) = XEPSF(3,3) C--------Les contraintes SIGF(1) = SIG(1,1) SIGF(2) = SIG(2,2) SIGF(3) = SIG(3,3) SIGF(4) = SIG(1,2) ENDIF C-----Fin de l integration RETURN END