1 : C TUFIEN SOURCE CHAT 05/01/13 03:55:36 5004 2 : SUBROUTINE TUFIEN(XM,XP,DELTAM,DELTAP,THETA,XJP,PRECIS,THETA0, 3 : & XM0,XP0,XJ1C,T,RAYOM,EPAI,YOUN) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : C======================================================================= 7 : C CE SOUS PROGRAMME SERT A CALCULER LE NOUVEL ANGLE DE LA 8 : C FISSURE DANS LE CAS D UN INCREMENT DE CHARGE POUR 9 : C LEQUEL IL Y A PROPAGATION EN ELASTICITE 10 : C SI ON PLASTIFIE ON RENVOIT LA PARTIE DE L INCREMENT 11 : C QU IL VA FALLOIR ECOULER 12 : C IL EST APPELE PAR TUFPLA 13 : C 14 : C ENTREES:XM VALEUR DU MOMENT AU DEBUT DU PAS 15 : C XP VALEUR DE L EFFORT NORMAL AU DEBUT DU PAS 16 : C DELTAM INCREMENT DE MOMENT 17 : C DELTAP INCREMENT D EFFORT NORMAL 18 : C THETA ANGLE DE LA FISSURE AU DEBUT DU PAS 19 : C XJP VALEUR DE LA PARTIE PLASTIQUE DE J 20 : C 21 : C SORTIES:XM VALEUR DU MOMENT A LA FIN DU PAS 22 : C XP VALEUR DE L EFFORT NORMAL A LA FIN DU PAS 23 : C DELTAM PARTIE DE L INCREMENT DE MOMENT RESTANT 24 : C A ECOULER EN PLASTICITE 25 : C DELTAP PARTIE DE L INCREMENT D EFFORT NORMAL 26 : C RESTANT A ECOULER EN PLASTICITE 27 : C THETA ANGLE DE LA FISSURE A LA FIN DU PAS 28 : C======================================================================= 29 : -INC CCREEL 30 : -INC CCOPTIO 31 : DIMENSION XX(2,2),YY(2,2) 32 : XM1=XM+DELTAM 33 : XP1=XP+DELTAP 34 : C 35 : C ON CALCULE LE NOUVEL ANGLE DE LA FISSURE PAR LA METHODE 36 : C DE NEWTON EN APPROXIMANT LA DERIVEE 37 : C 38 : X0=THETA 39 : X1=THETA+XPI/180.D0 40 : NBIT=0 41 : 1 CONTINUE 42 : F0=CRIT2(XM1,XP1,X0,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN) 43 : F1=CRIT2(XM1,XP1,X1,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN) 44 : FP=(F1-F0)/(X1-X0) 45 : X2=X1-F1/FP 46 : IF(X2.LT.THETA0) THEN 47 : IF(IIMPI.EQ.999)THEN 48 : WRITE(IOIMP,111) 49 : 111 FORMAT(1X,'L ANGLE EST INFERIEUR A L ANGLE INITIAL') 50 : ENDIF 51 : DELTAM=0.D0 52 : DELTAP=0.D0 53 : RETURN 54 : ENDIF 55 : IF(X2.GT.360.D0) THEN 56 : IF(IIMPI.EQ.999)THEN 57 : WRITE(IOIMP,222) 58 : 222 FORMAT(1X,'PAS DE SOLUTION REALISTE') 59 : ENDIF 60 : DELTAM=0.D0 61 : DELTAP=0.D0 62 : RETURN 63 : ENDIF 64 : NBIT=NBIT+1 65 : CRI=ABS((X2-X1)/X1) 66 : F2=CRIT2(XM1,XP1,X2,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN) 67 : CRU=ABS(F2) 68 : IF(CRI.GE.PRECIS.OR.CRU.GE.PRECIS) THEN 69 : X0=X1 70 : X1=X2 71 : GOTO 1 72 : ELSE 73 : THETA=X2 74 : IF(IIMPI.EQ.999)THEN 75 : WRITE(IOIMP,*)THETA 76 : WRITE(IOIMP,*)'NBRE ITERATIONS',NBIT 77 : VER=CRIT2(XM1,XP1,THETA,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN) 78 : WRITE(IOIMP,900)VER 79 : 900 FORMAT(1X,'VERIF RESOLUTION',1X,F12.5) 80 : ENDIF 81 : C 82 : C THETA EST LE NOUVEL ANGLE DE LA FISSURE 83 : C 84 : ENDIF 85 : C 86 : C ON TESTE POUR SAVOIR SI ON EST TOUJOURS ELASTIQUE 87 : C 87 : CALL TUFIC1(XM1,XP1,THETA,IRETOU,PRECIS,XM0,XP0) 88 : C 89 : C SI OUI ON A TOUT ECOULE ET ON MET DELTAM ET DELTAP 90 : C A ZERO 91 : C 92 : IF(IRETOU.EQ.0) THEN 93 : XM=XM+DELTAM 94 : XP=XP+DELTAP 95 : DELTAM=0.D0 96 : DELTAP=0.D0 97 : RETURN 98 : C 99 : C SINON ON CHERCHE L ANGLE QUI CORRESPOND A LA FRONTIERE 100 : C DU DOMAINE D ELASTICITE 101 : C ON RESOUD UN SYSTEME DU TYPE F(A,B)=G(A,B)=0 D UNE 102 : C MANIERE ANALOGUE A LA METHODE DE NEWTON 103 : C 104 : ELSE 105 : IF(IIMPI.EQ.999)THEN 106 : WRITE(IOIMP,901) 107 : 901 FORMAT(1X,'MAIS ON PLASTIFIE') 108 : ENDIF 109 : NBIT=0 110 : A1=0.D0 111 : B1=1.D0 112 : XM1=XM+A1*DELTAM 113 : XP1=XP+A1*DELTAP 114 : THETA1=B1*THETA 115 : A2=1.D0 116 : B2=0.9D0 117 : XM2=XM+A2*DELTAM 118 : XP2=XP+A2*DELTAP 119 : THETA2=B2*THETA 120 : 2 CONTINUE 121 : F11=CRIT1(XM1,XP1,THETA1,XM0,XP0) 122 : F21=CRIT1(XM2,XP2,THETA1,XM0,XP0) 123 : XX(1,1)=(F21-F11)/(A2-A1) 124 : F12=CRIT1(XM1,XP1,THETA2,XM0,XP0) 125 : XX(1,2)=(F12-F11)/(B2-B1) 126 : F22=CRIT1(XM2,XP2,THETA2,XM0,XP0) 127 : G11=CRIT2(XM1,XP1,THETA1,XJP,XJ1C,T,RAYOM,THETA0,EPAI, 128 : & YOUN) 129 : G21=CRIT2(XM2,XP2,THETA1,XJP,XJ1C,T,RAYOM,THETA0,EPAI, 130 : & YOUN) 131 : XX(2,1)=(G21-G11)/(A2-A1) 132 : G12=CRIT2(XM1,XP1,THETA2,XJP,XJ1C,T,RAYOM,THETA0,EPAI, 133 : & YOUN) 134 : XX(2,2)=(G12-G11)/(B2-B1) 135 : G22=CRIT2(XM2,XP2,THETA2,XJP,XJ1C,T,RAYOM,THETA0,EPAI, 136 : & YOUN) 137 : DET=XX(1,1)*XX(2,2)-XX(1,2)*XX(2,1) 138 : YY(1,1)=XX(2,2)/DET 139 : YY(1,2)=-XX(1,2)/DET 140 : YY(2,1)=-XX(2,1)/DET 141 : YY(2,2)=XX(1,1)/DET 142 : A3=A2-YY(1,1)*F22-YY(1,2)*G22 143 : B3=B2-YY(2,1)*F22-YY(2,2)*G22 144 : XM3=XM+A3*DELTAM 145 : XP3=XP+A3*DELTAP 146 : THETA3=B3*THETA 147 : NBIT=NBIT+1 148 : CRA=ABS((A3-A2)/A2) 149 : CRO=ABS((B3-B2)/B2) 150 : CRU=CRA+CRO 151 : F3=CRIT1(XM3,XP3,THETA3,XM0,XP0) 152 : G3=CRIT2(XM3,XP3,THETA3,XJP,XJ1C,T,RAYOM,THETA0,EPAI, 153 : & YOUN) 154 : CRI=(F3**2+G3**2)**0.5D0 155 : IF(CRI.GE.PRECIS.OR.CRU.GE.PRECIS) THEN 156 : A1=A2 157 : B1=B2 158 : A2=A3 159 : B2=B3 160 : XM1=XM2 161 : XP1=XP2 162 : THETA1=THETA2 163 : XM2=XM3 164 : XP2=XP3 165 : THETA2=THETA3 166 : GOTO 2 167 : ELSE 168 : C 169 : C ON A ENFIN TROUVE A ET B QUI CONVIENNENT 170 : C 171 : A=A3 172 : B=B3 173 : ENDIF 174 : XM=XM+A*DELTAM 175 : XP=XP+A*DELTAP 176 : DELTAM=(1.D0-A)*DELTAM 177 : DELTAP=(1.D0-A)*DELTAP 178 : THETA=B*THETA 179 : VER1=CRIT1(XM,XP,THETA,XM0,XP0) 180 : VER2=CRIT2(XM,XP,THETA,XJP,XJ1C,T,RAYOM,THETA0,EPAI, 181 : & YOUN) 182 : IF(IIMPI.EQ.999)THEN 183 : WRITE(IOIMP,903)VER1,VER2 184 : 903 FORMAT(1X,'REVERIF RESO',1X,2F12.5) 185 : WRITE(IOIMP,*)'NBRE ITERATIONS',NBIT 186 : ENDIF 187 : RETURN 188 : ENDIF 189 : END 190 : 191 : 192 : 193 : 194 :
© Cast3M 2003 - All rights reserved.
Disclaimer