1 : C CODJFD SOURCE KICH 09/09/11 21:16:46 6505 2 : C MODJFD SOURCE AM 00/12/13 21:41:28 4045 3 : SUBROUTINE CODJFD(WRK52,WRK53,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,Y2, 4 : 1 JSIGN,SIG1,SIG2,CONV,DT1) 5 : C 6 : C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7 : C % % 8 : C % INTEGRATION DU MODJFD UNILATERAL AVEC FERMETURE LINEAIRE % 9 : C % DE FISSURES EN CONTRAINTES PLANES % 10 : C % (MODCLB NON ASSOCIE) % 11 : C % A L"AIDE DE LA METHODE DE NEWTON % 12 : C % AJOUT VICOENDOMMAGEMENT J.F. DUBE % 13 : C % % 14 : C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15 : C 16 : C 17 : C 18 : IMPLICIT INTEGER(I-N) 19 : -INC DECHE 20 : C VARIABLES D"ENTREES 21 : C 22 : C 23 : REAL*8 EPSP1,EPSP2 24 : C 25 : C VARIABLES DE SORTIE 26 : C 27 : C 28 : REAL*8 SIG1,SIG2,Z1,Z2,Y1,Y2,DT1 29 : INTEGER JSIGN(3) 30 : LOGICAL CONV 31 : C 32 : C 33 : C PARAMETRES D"EVOLUTION 34 : C 35 : C 36 : REAL*8 E,NU,A1,A2,B1,B2,Y01,Y02,SIGF1 37 : REAL*8 BETA1,BETA2,YLIM1,YLIM2 38 : C 39 : C 40 : C VARIABLES INTERNES 41 : C 42 : C 43 : REAL*8 Y001,Y002,DY01,DY02,hdp1,hdp2 44 : REAL*8 DELTASIG 45 : INTEGER NDIM 46 : PARAMETER (NDIM=6) 47 : REAL*8 JACOB(NDIM,NDIM),J15,J16,J25,J26 48 : REAL*8 DELAFON(NDIM) 49 : REAL*8 X(NDIM),DX(NDIM),X0(NDIM),XX,XIG1,XIG2 50 : REAL*8 ERREUR,TEST 51 : REAL*8 H1,H2,HT,HT2,trf 52 : REAL*8 ALFA1,ALFA2 53 : REAL*8 SIGY 54 : REAL*8 UN,DEUX 55 : REAL*8 ZMAX 56 : INTEGER ITER ,I,DIMPB 57 : LOGICAL LDOM1,LDOM2,evol1,evol2 58 : REAL*8 M1,M2,ALPHA1,ALPHA2 59 : REAL*8 x3,x4,x5,x6,dx5,dx6 60 : UN=1.D0 61 : DEUX=2.D0 62 : ZMAX=1.d+04 63 : E = XMAT(1) 64 : NU = XMAT(2) 65 : Y001 = XMAT(5) 66 : Y002 = XMAT(6) 67 : A1 = XMAT(7) 68 : A2 = XMAT(8) 69 : B1 = XMAT(9) 70 : B2 = XMAT(10) 71 : BETA1 = XMAT(11) 72 : BETA2 = XMAT(12) 73 : SIGF1 = XMAT(13) 74 : M1 = XMAT(14) 75 : M2 = XMAT(15) 76 : ALPHA1= XMAT(16) 77 : ALPHA2= XMAT(17) 78 : SIGY =DSQRT (2.D0*E*Y001) 79 : * DELTASIG=SIGY*1.D-4 80 : DELTASIG=0.d0 81 : * print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2' 82 : * print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2 83 : TEST=1.D-12 84 : ERREUR=1.D0 85 : ITER=0 86 : X(1)=SIG1 87 : X(2)=SIG2 88 : X(3)=Y1 89 : X(4)=Y2 90 : X(5)=Z1 91 : X(6)=Z2 92 : DO I=1, NDIM 93 : X0(I)=X(I) 94 : END DO 95 : x5=Z1 96 : x6=z2 97 : C 98 : C 99 : C 100 : JACOB(1,1) = -1.d0 101 : JACOB(1,2) = 0.D0 102 : JACOB(1,3) = 0.D0 103 : JACOB(1,4) = 0.D0 104 : JACOB(1,5) = 0.D0 105 : JACOB(1,6) = 0.D0 106 : JACOB(2,1) = 0.D0 107 : JACOB(2,2) = -1.d0 108 : JACOB(2,3) = 0.D0 109 : JACOB(2,4) = 0.D0 110 : JACOB(2,5) = 0.D0 111 : JACOB(2,6) = 0.D0 112 : JACOB(3,1) = 0.D0 113 : JACOB(3,2) = 0.D0 114 : JACOB(3,3) = -1.d0 115 : JACOB(3,4) = 0.D0 116 : JACOB(3,5) = 0.D0 117 : JACOB(3,6) = 0.D0 118 : JACOB(4,1) = 0.D0 119 : JACOB(4,2) = 0.D0 120 : JACOB(4,3) = 0.D0 121 : JACOB(4,4) = -1.d0 122 : JACOB(4,5) = 0.D0 123 : JACOB(4,6) = 0.D0 124 : JACOB(5,1) = 0.D0 125 : JACOB(5,2) = 0.D0 126 : JACOB(5,3) = 0.D0 127 : JACOB(5,4) = 0.D0 128 : JACOB(5,5) = -1.d0 129 : JACOB(5,6) = 0.D0 130 : JACOB(6,1) = 0.D0 131 : JACOB(6,2) = 0.D0 132 : JACOB(6,3) = 0.D0 133 : JACOB(6,4) = 0.D0 134 : JACOB(6,5) = 0.D0 135 : JACOB(6,6) = -1.d0 136 : 137 : DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 100) 138 : ITER=ITER+1 139 : IF (X(1) .GE. DELTASIG) THEN 140 : H1=UN 141 : ELSE IF (X(1) .LE. -DELTASIG) THEN 142 : H1=0.D0 143 : ELSE 144 : H1=(X(1)+DELTASIG)/(DEUX*DELTASIG) 145 : END IF 146 : IF (X(2) .GE. DELTASIG) THEN 147 : H2=UN 148 : ELSE IF (X(2) .LE. -DELTASIG) THEN 149 : H2=0.D0 150 : ELSE 151 : H2=(X(2)+DELTASIG)/(DEUX*DELTASIG) 152 : END IF 153 : c 154 : c smooth transition on ht ht2 155 : c 156 : trf = x(1) + x(2) 157 : IF(trf.GE.deltasig)THEN 158 : HT=UN 159 : HT2=0.D0 160 : else if(trf.ge.(-deltasig)) then 161 : ht = (trf+deltasig)/DEUX/deltasig 162 : ht2 = UN-ht 163 : else if(trf.ge.(-sigf+deltasig)) then 164 : ht = 0.D+00 165 : ht2 = UN 166 : else if(trf.ge.(-sigf-deltasig)) then 167 : ht = 0.D+00 168 : c ht2 =UN- (trf-(sigf-deltasig))/(DEUX*sigf) 169 : ht2 = (trf+sigf+deltasig)/DEUX/deltasig 170 : else 171 : ht = 0.D+00 172 : ht2 = 0.D+00 173 : endif 174 : c print*,'DTjfd',dt1 175 : dy01=(ddim(x(5),z1)/(x(5)*z1*dt1)) 176 : hdp1=1.0 177 : if (dy01 .gt. 0.0) hdp1=0.0 178 : y01=y001+m1*dy01**alpha1 179 : dy02=(ddim(x(6),z2)/(x(6)*z2*dt1)) 180 : hdp2=1.0 181 : if (dy02 .gt. 0.0) hdp2=0.0 182 : y02=y002+m2*dy02**alpha2 183 : c print*,'iter=',iter,' dy01=',dy01,'dy02=',dy02 184 : 185 : ldom1=.false. 186 : evol1=.false. 187 : if (X(3).GT. Y01 .AND. X(5) .LT. ZMAX) then 188 : if ( (1.+(A1*(-Y01 + X(3)))**B1-z1) .gt. 1.d-06) then 189 : evol1=.true. 190 : ldom1=.true. 191 : x5=z1 192 : endif 193 : endif 194 : ldom2=.false. 195 : evol2=.false. 196 : if (X(4).GT. Y02 .AND. X(6) .LT. ZMAX) then 197 : if ( (1.+(A2*(-Y02 + X(4)))**B2-z2) .gt. 1.d-06) then 198 : evol2=.true. 199 : ldom2=.true. 200 : x6=z2 201 : endif 202 : endif 203 : c print*,'y1',x(3),'ylim1',ylim1 204 : c print*,'z1',x(5),'x5',x5 205 : c print*,'ldom1',ldom1,'ldom2',ldom2 206 : C =====> PROBLEME ELASTIQUE 207 : C 208 : C 209 : C CALCUL DU JACOBIEN 210 : C 211 : C 212 : JACOB(3,1) = H1*X(1)*X(5)*X(5)/E 213 : C 214 : JACOB(3,2) = H2*X(2)*X(5)*X(5)/E 215 : C 216 : JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E 217 : C 218 : JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E 219 : C 220 : C CALCUL DU SECOND MEMBRE 221 : C 222 : XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN) 223 : 1 *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2) 224 : 2 *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+ 225 : 3 BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1* 226 : 4 (X(6)*(H2-UN)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+ 227 : 5 (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)- 228 : 6 X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+ 229 : 7 (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)* 230 : 8 (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5)) 231 : 9 *BETA1*H1)*HT2) 232 : C 233 : XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN) 234 : 1 *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)* 235 : 2 X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+ 236 : 3 (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+ 237 : 4 (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))* 238 : 5BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1) 239 : 6 *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2* 240 : 7 (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))* 241 : 8 BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)* 242 : 9 (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2) 243 : 244 : DELAFON(1) =X(1)-XIG1 245 : C 246 : DELAFON(2) =X(2)-XIG2 247 : C 248 : XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E) 249 : DELAFON(3) = X(3)-XX 250 : C 251 : XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) ) 252 : 1 *X(6)*X(6)/(DEUX*E) 253 : DELAFON(4) = X(4)-XX 254 : C 255 : IF (LDOM1) THEN 256 : C 257 : C =====> ENDOMMAGEMENT EN TRACTION 258 : C 259 : J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2 260 : & ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H 261 : & 2)-E*EPSP1*H2)*SIGF1+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+ 262 : & H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6 263 : & )*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1 264 : & *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN 265 : & )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU 266 : & X*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H 267 : & 2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5 268 : & )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+ 269 : & BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1* 270 : & (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*( 271 : & X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X 272 : & (5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS 273 : & P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2* 274 : & (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BE 275 : & TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX) 276 : & +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 277 : C 278 : C 279 : J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6 280 : & )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X( 281 : & 6)-UN)*BETA2-E*EPSP2)*H1)*SIGF1+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT 282 : & 2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2 283 : & ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT 284 : & 2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX) 285 : & )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX 286 : & *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(- 287 : & DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H 288 : & T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN) 289 : & *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X( 290 : & 6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN) 291 : & *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)* 292 : & X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF1+(B 293 : & ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2 294 : & +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X( 295 : & 6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA 296 : & 1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U 297 : & N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 298 : C 299 : JACOB(1,5) = J15 300 : JACOB(2,5) = J25 301 : JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E 302 : JACOB(4,5) = 0.D0 303 : JACOB(5,5) = -1.D0 304 : & -A1*B1*(A1*(X(3)-Y01))**(B1-UN) 305 : & *alpha1*m1*dy01**alpha1/(dy01+hdp1)/(x(5)**2*dt1) 306 : JACOB(5,1) = 0.D0 307 : JACOB(5,2) = 0.D0 308 : JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN) 309 : JACOB(5,4) = 0.D0 310 : DELAFON(5) = X(5)-((A1*(X(3)-Y01))**B1+UN) 311 : 312 : IF (LDOM2) THEN 313 : C 314 : C =====> ENDOMMAGEMENT EN COMPRESSION 315 : C 316 : J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2) 317 : 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/ 318 : 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)* 319 : 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X 320 : 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5) 321 : 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5 322 : 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE 323 : 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X 324 : 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5) 325 : 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+ 326 : : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6 327 : ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE 328 : < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6 329 : = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1 330 : > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+( 331 : ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 332 : C 333 : J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B 334 : 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1- 335 : 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X 336 : 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2- 337 : 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5) 338 : 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X 339 : 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1* 340 : 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H 341 : 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6) 342 : 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6 343 : : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2 344 : ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS 345 : < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5) 346 : = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2* 347 : > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5)) 348 : ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2) 349 : @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 350 : C 351 : C 352 : JACOB(1,6) = J16 353 : JACOB(2,6) = J26 354 : JACOB(3,6) = 0.D0 355 : JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2) 356 : & *X(2)*X(2))*X(6)/E 357 : JACOB(5,6) = 0.D0 358 : JACOB(6,6) = -1.D0 359 : & -A2*B2*(A2*(X(4)-Y02))**(B2-UN) 360 : & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1) 361 : JACOB(6,1) = 0.D0 362 : JACOB(6,2) = 0.D0 363 : JACOB(6,3) = 0.D0 364 : * JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0)) 365 : JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0) 366 : JACOB(6,5) = 0.D0 367 : * DELAFON(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0) 368 : DELAFON(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0) 369 : DIMPB=6 370 : ELSE 371 : DIMPB=5 372 : END IF 373 : ELSE 374 : C 375 : C =====> PAS D'ENDOMMAGEMENT EN TRACTION 376 : C 377 : IF (LDOM2) THEN 378 : C 379 : C =====> ENDOMMAGEMENT EN COMPRESSION 380 : C 381 : J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2) 382 : 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/ 383 : 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)* 384 : 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X 385 : 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5) 386 : 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5 387 : 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE 388 : 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X 389 : 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5) 390 : 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+ 391 : : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6 392 : ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE 393 : < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6 394 : = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1 395 : > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+( 396 : ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 397 : C 398 : J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B 399 : 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1- 400 : 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X 401 : 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2- 402 : 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5) 403 : 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X 404 : 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1* 405 : 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H 406 : 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6) 407 : 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6 408 : : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2 409 : ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS 410 : < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5) 411 : = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2* 412 : > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5)) 413 : ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2) 414 : @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 415 : C 416 : C 417 : 418 : JACOB(1,5) = J16 419 : JACOB(2,5) = J26 420 : JACOB(3,5) = 0.D0 421 : JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2) 422 : & *X(2)*X(2))*X(6)/E 423 : JACOB(5,5) = -1.D0 424 : & -A2*B2*(A2*(X(4)-Y02))**(B2-UN) 425 : & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1) 426 : JACOB(5,1) = 0.D0 427 : JACOB(5,2) = 0.D0 428 : JACOB(5,3) = 0.D0 429 : * JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0)) 430 : JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0) 431 : * DELAFON(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0) 432 : DELAFON(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0) 433 : DIMPB=5 434 : ELSE 435 : DIMPB=4 436 : END IF 437 : END IF 438 : C 439 : C RESOLUTION 440 : C 440 : CALL GAUSCL(JACOB,DX,DELAFON,NDIM,DIMPB,0) 441 : C 442 : C CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION 443 : C 444 : ERREUR=0.D0 445 : DO I=1,2 446 : ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY) 447 : X(I)=X(I)+DX(I) 448 : if (x(i) .ge. 0.) then 449 : x(i)=dmin1(1.d+15,x(i)) 450 : else 451 : x(i)=dmax1(-1.d+15,x(i)) 452 : endif 453 : END DO 454 : X(3)=X(3)+DX(3) 455 : X3=DMAX1(X(3),Y001) 456 : X3=Dmin1(dabs(X3),1.d+07) 457 : ERREUR=ERREUR+DX(3)*DX(3)/(x3*x3) 458 : X(4)=X(4)+DX(4) 459 : X4=DMAX1(X(4),Y002) 460 : X4=Dmin1(dabs(X4),1.d+07) 461 : ERREUR=ERREUR+DX(4)*DX(4)/(x4*x4) 462 : if (.not. ldom1) then 463 : dx5=(x(5)-x5)*1.0/4.0 464 : x(5)=dx5+x5 465 : x(5)=dmax1(z1,x(5)) 466 : ERREUR=ERREUR+DX5*DX5/(z1*z1) 467 : endif 468 : if (.not. ldom2) then 469 : dx6=(x(6)-x6)/4. 470 : x(6)=dx6+x6 471 : x(6)=dmax1(z2,x(6)) 472 : ERREUR=ERREUR+DX6*DX6/(z2*z2) 473 : endif 474 : IF (LDOM1) THEN 475 : IF ((X(5)+DX(5)).LT.Z1) THEN 476 : DX(5)=(Z1-X(5)) 477 : 478 : c if (evol1) dx(5)=dx(5)/2. 479 : 480 : END IF 481 : IF ( X(5) .GT. ZMAX)THEN 482 : DX(5)=ZMAX + 0.1D0 - X(5) 483 : END IF 484 : X(5)=X(5)+DX(5) 485 : ERREUR=ERREUR+DX(5)*DX(5)/(z1*z1) 486 : IF (LDOM2) THEN 487 : IF ((X(6)+DX(6)).LT.Z2) THEN 488 : DX(6)=(Z2-X(6)) 489 : 490 : c if (evol2) dx(6)=dx(6)/2. 491 : 492 : END IF 493 : IF ( X(6) .GT. ZMAX)THEN 494 : DX(6)=ZMAX + 0.1D0 - X(6) 495 : END IF 496 : X(6)=X(6)+DX(6) 497 : ERREUR=ERREUR+DX(6)*DX(6)/(z2*z2) 498 : END IF 499 : ELSE IF (LDOM2) THEN 500 : IF ((X(6)+DX(5)).LT.Z2) THEN 501 : DX(5)=(Z2-X(6)) 502 : 503 : c if (evol2) dx(5)=dx(5)/2. 504 : 505 : END IF 506 : IF ( X(6) .GT. ZMAX)THEN 507 : DX(5)=ZMAX + 0.1D0 - X(6) 508 : END IF 509 : X(6)=X(6)+DX(5) 510 : ERREUR=ERREUR+DX(5)*DX(5)/(z2*z2) 511 : END IF 512 : 513 : 514 : IF (ERREUR .GT. 1.D+100)THEN 515 : CONV=.FALSE. 516 : print*,'erreur dans modele',erreur 517 : print*,'iteration num',iter 518 : print*,'delafon',delafon 519 : print*,z1,z2,x3,x4,sigy 520 : print*,dx 521 : print*,'y02=',y02,'y002=',y002 522 : RETURN 523 : END IF 524 : C 525 : C QUELQUES IMPRESSIONS EN CAS D'ERREUR 526 : C 527 : if (.not. conv )then 528 : PRINT*,'ITERATION ',ITER,'ERREUR',erreur 529 : print*,'LDOM1 :',LDOM1,' LDOM2 :',LDOM2 530 : PRINT*,'EPS1= ',EPSP1,' EPS2= ',EPSP2 531 : PRINT*,'SIG1:',X(1) ,' DSIG1:',DX(1),'ERR :',DELAFON(1) 532 : PRINT*,'SIG2:',X(2) ,' DSIG2:',DX(2),'ERR :',DELAFON(2) 533 : PRINT*,'Y1 :',X(3) ,' DY1 :',DX(3),'ERR :',DELAFON(3) 534 : PRINT*,'Y2 :',X(4) ,' DY2 :',DX(4),'ERR :',DELAFON(4) 535 : PRINT*,'Z1 :',X(5) ,' DZ1 :',DX(5),'ERR :',DELAFON(5) 536 : PRINT*,'Z2 :',X(6) ,' DZ2 :',DX(6),'ERR :',DELAFON(6) 537 : PRINT*,'H1',H1,'H2',H2,ht,ht2 538 : PRINT* 539 : write(*,200)jacob 540 : END IF 541 : 542 : END DO 543 : IF (ERREUR .GT. TEST ) THEN 544 : CONV=.FALSE. 545 : ELSE 546 : CONV=.TRUE. 547 : SIG1= X(1) 548 : SIG2= X(2) 549 : Y1 = X(3) 550 : Y2 = X(4) 551 : Z1 = X(5) 552 : Z2 = X(6) 553 : IF (SIG1 .GE. 0.D0) THEN 554 : JSIGN(1)=1 555 : ELSE 556 : JSIGN(1)=0 557 : END IF 558 : IF (SIG2.GE. 0.D0)THEN 559 : JSIGN(2)=1 560 : ELSE 561 : JSIGN(2)=0 562 : END IF 563 : JSIGN(3)=1 564 : END IF 565 : RETURN 566 : 200 format(6(6(2x,e11.4),/)) 567 : END 568 : 569 : 570 : 571 : 572 : 573 : 574 :
© Cast3M 2003 - All rights reserved.
Disclaimer