1 : C CODJF2 SOURCE KICH 09/09/11 21:16:45 6505 2 : C MODJF2 SOURCE AM 00/12/13 21:41:10 4045 3 : SUBROUTINE CODJF2(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 MODJFD2 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 DICHOTOMIE % 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 : INTEGER IMAX 61 : REAL*8 pas1,erreurd1,xd1,xg1,pas2,erreurd2,xd2,xg2 62 : REAL*8 erreur_max,fctd1,fctd2 63 : logical passage1,passage2 64 : UN=1.D0 65 : DEUX=2.D0 66 : ZMAX=1.d+04 67 : E = XMAT(1) 68 : NU = XMAT(2) 69 : Y001 = XMAT(5) 70 : Y002 = XMAT(6) 71 : A1 = XMAT(7) 72 : A2 = XMAT(8) 73 : B1 = XMAT(9) 74 : B2 = XMAT(10) 75 : BETA1 = XMAT(11) 76 : BETA2 = XMAT(12) 77 : SIGF1 = XMAT(13) 78 : M1 = XMAT(14) 79 : M2 = XMAT(15) 80 : ALPHA1= XMAT(16) 81 : ALPHA2= XMAT(17) 82 : SIGY =DSQRT (2.D0*E*Y001) 83 : * DELTASIG=SIGY*1.D-4 84 : DELTASIG=0.d0 85 : * print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2' 86 : * print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2 87 : TEST=1.D-12 88 : ERREUR=1.D0 89 : ITER=0 90 : X(1)=SIG1 91 : X(2)=SIG2 92 : X(3)=Y1 93 : X(4)=Y2 94 : X(5)=Z1 95 : X(6)=Z2 96 : DO I=1, NDIM 97 : X0(I)=X(I) 98 : END DO 99 : x5=Z1 100 : x6=z2 101 : C 102 : C 103 : C 104 : JACOB(1,1) = -1.d0 105 : JACOB(1,2) = 0.D0 106 : JACOB(1,3) = 0.D0 107 : JACOB(1,4) = 0.D0 108 : JACOB(1,5) = 0.D0 109 : JACOB(1,6) = 0.D0 110 : JACOB(2,1) = 0.D0 111 : JACOB(2,2) = -1.d0 112 : JACOB(2,3) = 0.D0 113 : JACOB(2,4) = 0.D0 114 : JACOB(2,5) = 0.D0 115 : JACOB(2,6) = 0.D0 116 : JACOB(3,1) = 0.D0 117 : JACOB(3,2) = 0.D0 118 : JACOB(3,3) = -1.d0 119 : JACOB(3,4) = 0.D0 120 : JACOB(3,5) = 0.D0 121 : JACOB(3,6) = 0.D0 122 : JACOB(4,1) = 0.D0 123 : JACOB(4,2) = 0.D0 124 : JACOB(4,3) = 0.D0 125 : JACOB(4,4) = -1.d0 126 : JACOB(4,5) = 0.D0 127 : JACOB(4,6) = 0.D0 128 : JACOB(5,1) = 0.D0 129 : JACOB(5,2) = 0.D0 130 : JACOB(5,3) = 0.D0 131 : JACOB(5,4) = 0.D0 132 : JACOB(5,5) = -1.d0 133 : JACOB(5,6) = 0.D0 134 : JACOB(6,1) = 0.D0 135 : JACOB(6,2) = 0.D0 136 : JACOB(6,3) = 0.D0 137 : JACOB(6,4) = 0.D0 138 : JACOB(6,5) = 0.D0 139 : JACOB(6,6) = -1.d0 140 : Erreur_max=1.0d-10 141 : pas1=1.0 142 : Xd1=z1 143 : Xg1=z1 144 : imax=0 145 : Erreurd1=1.0 146 : passage1=.true. 147 : Do while ( (Erreurd1 .GT. Erreur_max) .and. (pas1 .ge. 1.d-6) ) 148 : x(5)=xd1+pas1 149 : 150 : c print*,'*** Z1 ***' 151 : c print*,x(5) 152 : 153 : ** Calcul de Y01 modifie par l'effet de vitesse 154 : c dp1=(ddim(x(5),z1)/(x(5)*z1*dt1)) 155 : c y01=y001+m1*dp1**alpha1 156 : 157 : pas2=1.0 158 : Xd2=z2 159 : Xg2=x(6) 160 : Erreurd2=1.0 161 : passage2=.true. 162 : X(4)=Y2 163 : Do while ((Erreurd2 .GT. Erreur_max) .and. (pas2 .ge. 1.d-6)) 164 : imax=imax+1 165 : if (imax .ge. 10000) then 166 : print*,'pb dans jfd2' 167 : conv=.false. 168 : return 169 : end if 170 : x(6)=xd2+pas2 171 : 172 : c print*,'*** Z2 ***' 173 : c print*,x(6) 174 : 175 : ** Calcul de Y02 modifie par l'effet de vitesse 176 : c dp2=(ddim(x(6),z2)/(x(6)*z2*dt1)) 177 : c y02=y002+m2*dp2**alpha2 178 : ** Calcul de sigma et Y1 et Y2 179 : ERREUR=1.D0 180 : ITER=0 181 : X(1)=Sig1 182 : X(2)=Sig2 183 : DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 100) 184 : ITER=ITER+1 185 : IF (X(1) .GE. DELTASIG) THEN 186 : H1=UN 187 : ELSE IF (X(1) .LE. -DELTASIG) THEN 188 : H1=0.D0 189 : ELSE 190 : H1=(X(1)+DELTASIG)/(DEUX*DELTASIG) 191 : END IF 192 : IF (X(2) .GE. DELTASIG) THEN 193 : H2=UN 194 : ELSE IF (X(2) .LE. -DELTASIG) THEN 195 : H2=0.D0 196 : ELSE 197 : H2=(X(2)+DELTASIG)/(DEUX*DELTASIG) 198 : END IF 199 : c 200 : c smooth transition on ht ht2 201 : c 202 : trf = x(1) + x(2) 203 : IF(trf.GE.deltasig)THEN 204 : HT=UN 205 : HT2=0.D0 206 : else if(trf.ge.(-deltasig)) then 207 : ht = (trf+deltasig)/DEUX/deltasig 208 : ht2 = UN-ht 209 : else if(trf.ge.(-sigf+deltasig)) then 210 : ht = 0.D+00 211 : ht2 = UN 212 : else if(trf.ge.(-sigf-deltasig)) then 213 : ht = 0.D+00 214 : c ht2 =UN- (trf-(sigf-deltasig))/(DEUX*sigf) 215 : ht2 = (trf+sigf+deltasig)/DEUX/deltasig 216 : else 217 : ht = 0.D+00 218 : ht2 = 0.D+00 219 : endif 220 : c print*,'DTjfd',dt1 221 : ** Calcul de Y01 modifie par l'effet de vitesse 222 : dy01=(ddim(x(5),z1)/(x(5)*z1*dt1)) 223 : hdp1=1.0 224 : if (dy01 .gt. 0.0) hdp1=0.0 225 : y01=y001+m1*dy01**alpha1 226 : ** Calcul de Y02 modifie par l'effet de vitesse 227 : dy02=(ddim(x(6),z2)/(x(6)*z2*dt1)) 228 : hdp2=1.0 229 : if (dy02 .gt. 0.0) hdp2=0.0 230 : y02=y002+m2*dy02**alpha2 231 : c print*,'dans clbeco iter=',iter,' y01=',y01,'y02=',y02 232 : 233 : ldom1=.false. 234 : evol1=.false. 235 : ldom2=.false. 236 : evol2=.false. 237 : c print*,'y1',x(3),'ylim1',ylim1 238 : c print*,'z1',x(5),'x5',x5 239 : c print*,'ldom1',ldom1,'ldom2',ldom2 240 : C =====> PROBLEME ELASTIQUE 241 : C 242 : C 243 : C CALCUL DU JACOBIEN 244 : C 245 : C 246 : JACOB(3,1) = H1*X(1)*X(5)*X(5)/E 247 : C 248 : JACOB(3,2) = H2*X(2)*X(5)*X(5)/E 249 : C 250 : JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E 251 : C 252 : JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E 253 : C 254 : C CALCUL DU SECOND MEMBRE 255 : C 256 : XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN) 257 : 1 *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2) 258 : 2 *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+ 259 : 3 BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1* 260 : 4 (X(6)*(H2-UN)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+ 261 : 5 (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)- 262 : 6 X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+ 263 : 7 (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)* 264 : 8 (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5)) 265 : 9 *BETA1*H1)*HT2) 266 : C 267 : XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN) 268 : 1 *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)* 269 : 2 X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+ 270 : 3 (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+ 271 : 4 (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))* 272 : 5 BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1) 273 : 6 *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2* 274 : 7 (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))* 275 : 8 BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)* 276 : 9 (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2) 277 : 278 : DELAFON(1) =X(1)-XIG1 279 : C 280 : DELAFON(2) =X(2)-XIG2 281 : C 282 : XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E) 283 : DELAFON(3) = X(3)-XX 284 : C 285 : XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) ) 286 : 1 *X(6)*X(6)/(DEUX*E) 287 : DELAFON(4) = X(4)-XX 288 : C 289 : IF (LDOM1) THEN 290 : C 291 : C =====> ENDOMMAGEMENT EN TRACTION 292 : C 293 : J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2 294 : & ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H 295 : & 2)-E*EPSP1*H2)*SIGF1+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+ 296 : & H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6 297 : & )*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1 298 : & *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN 299 : & )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU 300 : & X*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H 301 : & 2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5 302 : & )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+ 303 : & BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1* 304 : & (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*( 305 : & X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X 306 : & (5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS 307 : & P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2* 308 : & (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BE 309 : & TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX) 310 : & +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 311 : C 312 : C 313 : J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6 314 : & )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X( 315 : & 6)-UN)*BETA2-E*EPSP2)*H1)*SIGF1+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT 316 : & 2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2 317 : & ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT 318 : & 2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX) 319 : & )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX 320 : & *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(- 321 : & DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H 322 : & T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN) 323 : & *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X( 324 : & 6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN) 325 : & *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)* 326 : & X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF1+(B 327 : & ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2 328 : & +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X( 329 : & 6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA 330 : & 1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U 331 : & N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 332 : C 333 : JACOB(1,5) = J15 334 : JACOB(2,5) = J25 335 : JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E 336 : JACOB(4,5) = 0.D0 337 : JACOB(5,5) = -1.D0 338 : & -A1*B1*(A1*(X(3)-Y01))**(B1-UN) 339 : & *alpha1*m1*dy01**alpha1/(dy01+hdp1)/(x(5)**2*dt1) 340 : JACOB(5,1) = 0.D0 341 : JACOB(5,2) = 0.D0 342 : JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN) 343 : JACOB(5,4) = 0.D0 344 : DELAFON(5) = X(5)-((A1*(X(3)-Y01))**B1+UN) 345 : 346 : IF (LDOM2) THEN 347 : C 348 : C =====> ENDOMMAGEMENT EN COMPRESSION 349 : C 350 : J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2) 351 : 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/ 352 : 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)* 353 : 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X 354 : 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5) 355 : 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5 356 : 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE 357 : 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X 358 : 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5) 359 : 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+ 360 : : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6 361 : ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE 362 : < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6 363 : = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1 364 : > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+( 365 : ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 366 : C 367 : J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B 368 : 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1- 369 : 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X 370 : 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2- 371 : 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5) 372 : 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X 373 : 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1* 374 : 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H 375 : 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6) 376 : 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6 377 : : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2 378 : ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS 379 : < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5) 380 : = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2* 381 : > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5)) 382 : ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2) 383 : @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 384 : C 385 : C 386 : JACOB(1,6) = J16 387 : JACOB(2,6) = J26 388 : JACOB(3,6) = 0.D0 389 : JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2) 390 : & *X(2)*X(2))*X(6)/E 391 : JACOB(5,6) = 0.D0 392 : JACOB(6,6) = -1.D0 393 : & -A2*B2*(A2*(X(4)-Y02))**(B2-UN) 394 : & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1) 395 : JACOB(6,1) = 0.D0 396 : JACOB(6,2) = 0.D0 397 : JACOB(6,3) = 0.D0 398 : * JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0)) 399 : JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0) 400 : JACOB(6,5) = 0.D0 401 : * DELAFON(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0) 402 : DELAFON(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0) 403 : DIMPB=6 404 : ELSE 405 : DIMPB=5 406 : END IF 407 : ELSE 408 : C 409 : C =====> PAS D'ENDOMMAGEMENT EN TRACTION 410 : C 411 : IF (LDOM2) THEN 412 : C 413 : C =====> ENDOMMAGEMENT EN COMPRESSION 414 : C 415 : J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2) 416 : 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/ 417 : 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)* 418 : 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X 419 : 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5) 420 : 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5 421 : 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE 422 : 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X 423 : 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5) 424 : 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+ 425 : : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6 426 : ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE 427 : < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6 428 : = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1 429 : > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+( 430 : ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 431 : C 432 : J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B 433 : 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1- 434 : 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X 435 : 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2- 436 : 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5) 437 : 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X 438 : 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1* 439 : 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H 440 : 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6) 441 : 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6 442 : : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2 443 : ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS 444 : < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5) 445 : = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2* 446 : > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5)) 447 : ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2) 448 : @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 449 : C 450 : C 451 : 452 : JACOB(1,5) = J16 453 : JACOB(2,5) = J26 454 : JACOB(3,5) = 0.D0 455 : JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2) 456 : & *X(2)*X(2))*X(6)/E 457 : JACOB(5,5) = -1.D0 458 : & -A2*B2*(A2*(X(4)-Y02))**(B2-UN) 459 : & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1) 460 : JACOB(5,1) = 0.D0 461 : JACOB(5,2) = 0.D0 462 : JACOB(5,3) = 0.D0 463 : * JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0)) 464 : JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0) 465 : * DELAFON(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0) 466 : DELAFON(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0) 467 : DIMPB=5 468 : ELSE 469 : DIMPB=4 470 : END IF 471 : END IF 472 : C 473 : C RESOLUTION 474 : C 474 : CALL GAUSCL(JACOB,DX,DELAFON,NDIM,DIMPB,0) 475 : C 476 : C CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION 477 : C 478 : ERREUR=0.D0 479 : do i=1,ndim 480 : erreur=erreur+(delafon(i)*delafon(i)) 481 : enddo 482 : ERREUR=0.D0 483 : DO I=1,2 484 : ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY) 485 : X(I)=X(I)+DX(I) 486 : c if (x(i) .ge. 0.) then 487 : c x(i)=dmin1(1.d+15,x(i)) 488 : c else 489 : c x(i)=dmax1(-1.d+15,x(i)) 490 : c endif 491 : END DO 492 : X(3)=X(3)+DX(3) 493 : X3=DMAX1(X(3),Y01) 494 : c X3=Dmin1(dabs(X3),1.d+07) 495 : ERREUR=ERREUR+DX(3)*DX(3)/(x3*x3) 496 : X(4)=X(4)+DX(4) 497 : X4=DMAX1(X(4),Y02) 498 : c X4=Dmin1(dabs(X4),1.d+07) 499 : ERREUR=ERREUR+DX(4)*DX(4)/(x4*x4) 500 : c if (.not. ldom1) then 501 : c dx5=(x(5)-x5)*1.0/4.0 502 : c x(5)=dx5+x5 503 : c x(5)=dmax1(z1,x(5)) 504 : c ERREUR=ERREUR+DX5*DX5/(z1*z1) 505 : c endif 506 : c if (.not. ldom2) then 507 : c dx6=(x(6)-x6)/4. 508 : c x(6)=dx6+x6 509 : c x(6)=dmax1(z2,x(6)) 510 : c ERREUR=ERREUR+DX6*DX6/(z2*z2) 511 : c endif 512 : IF (LDOM1) THEN 513 : IF ((X(5)+DX(5)).LT.Z1) THEN 514 : DX(5)=(Z1-X(5)) 515 : END IF 516 : IF ( X(5) .GT. ZMAX)THEN 517 : DX(5)=ZMAX + 0.1D0 - X(5) 518 : END IF 519 : X(5)=X(5)+DX(5) 520 : ERREUR=ERREUR+DX(5)*DX(5)/(z1*z1) 521 : IF (LDOM2) THEN 522 : IF ((X(6)+DX(6)).LT.Z2) THEN 523 : DX(6)=(Z2-X(6)) 524 : END IF 525 : IF ( X(6) .GT. ZMAX)THEN 526 : DX(6)=ZMAX + 0.1D0 - X(6) 527 : END IF 528 : X(6)=X(6)+DX(6) 529 : ERREUR=ERREUR+DX(6)*DX(6)/(z2*z2) 530 : END IF 531 : ELSE IF (LDOM2) THEN 532 : IF ((X(6)+DX(5)).LT.Z2) THEN 533 : DX(5)=(Z2-X(6)) 534 : END IF 535 : IF ( X(6) .GT. ZMAX)THEN 536 : DX(5)=ZMAX + 0.1D0 - X(6) 537 : END IF 538 : X(6)=X(6)+DX(5) 539 : ERREUR=ERREUR+DX(5)*DX(5)/(z2*z2) 540 : END IF 541 : 542 : 543 : IF (ERREUR .GT. 1.D+100)THEN 544 : CONV=.FALSE. 545 : print*,'erreur dans modele',erreur 546 : print*,'iteration num',iter 547 : print*,'delafon',delafon 548 : RETURN 549 : END IF 550 : C 551 : C QUELQUES IMPRESSIONS EN CAS D'ERREUR 552 : C 553 : if (.not. conv )then 554 : PRINT*,'ITERATION ',ITER,'ERREUR',erreur 555 : print*,'LDOM1 :',LDOM1,' LDOM2 :',LDOM2 556 : PRINT*,'EPS1= ',EPSP1,' EPS2= ',EPSP2 557 : PRINT*,'SIG1:',X(1) ,' DSIG1:',DX(1),'ERR :',DELAFON(1) 558 : PRINT*,'SIG2:',X(2) ,' DSIG2:',DX(2),'ERR :',DELAFON(2) 559 : PRINT*,'Y1 :',X(3) ,' DY1 :',DX(3),'ERR :',DELAFON(3) 560 : PRINT*,'Y2 :',X(4) ,' DY2 :',DX(4),'ERR :',DELAFON(4) 561 : PRINT*,'Z1 :',X(5) ,' DZ1 :',DX(5),'ERR :',DELAFON(5) 562 : PRINT*,'Z2 :',X(6) ,' DZ2 :',DX(6),'ERR :',DELAFON(6) 563 : PRINT*,'H1',H1,'H2',H2,ht,ht2 564 : PRINT* 565 : write(*,200)jacob 566 : END IF 567 : 568 : END DO 569 : 570 : *** pour z2 **** 571 : ******************** 572 : fctd2=0. 573 : If ( x(4) .Ge. Y02 ) Then 574 : ** si z converge est > a x(6) alors fctd est < 0 575 : Fctd2 =-( 1.0+(A2*(-Y02 + X(4)))**B2-X(6) ) 576 : if (fctd2.ge.0.0 .and. x(6).eq.z2) fctd2 =0.0 577 : Else 578 : if (x(4) .le. y002 ) then 579 : fctd2=x(6)-z2 580 : c x(6)=z2 581 : else 582 : Fctd2 =x(6)-z2 583 : endif 584 : Endif 585 : If ( Fctd2 .lt. 0.0) then 586 : Xg2=Xd2 587 : if (passage2) pas2=pas2*2.0 588 : Xd2=Xd2+pas2 589 : Else 590 : Xd2=Xg2 591 : Pas2=Pas2/2.0 592 : Xd2=Xd2+pas2 593 : passage2=.false. 594 : Endif 595 : Erreurd2=(fctd2*fctd2)/(z2*z2) 596 : if (x(6) .gt. zmax) then 597 : x(6)=zmax 598 : Erreurd2=0.0 599 : endif 600 : 601 : c print*,'Erreurd2' 602 : c print*,Erreurd2,x(4),y02,y002 603 : 604 : ENDDO 605 : 606 : *** pour z1 **** 607 : ******************** 608 : fctd1=0.0 609 : If ( x(3) .Ge. Y01 ) Then 610 : ** si z converge est > a x(5) alors fctd1 est < 0 611 : Fctd1 =-( 1.0+(A1*(-Y01 + X(3)))**B1-X(5) ) 612 : if (fctd1.ge.0.0 .and. x(5).eq.z1) fctd1=0.0 613 : Else 614 : if (x(3) .le. y001) then 615 : fctd1=x(5)-z2 616 : x(5)=z1 617 : else 618 : Fctd1 =x(5)-z2 619 : endif 620 : Endif 621 : If ( Fctd1 .lt. 0.0) then 622 : Xg1=Xd1 623 : if (passage1) pas1=pas1*2.0 624 : Xd1=Xd1+pas1 625 : Else 626 : Xd1=Xg1 627 : Pas1=Pas1/2.0 628 : Xd1=Xd1+pas1 629 : passage1=.false. 630 : Endif 631 : Erreurd1=(Fctd1*Fctd1)/(z1*z1) 632 : if (x(5) .gt. zmax) then 633 : x(5)=zmax 634 : Erreurd1=0.0 635 : endif 636 : 637 : c print*,'Erreurd1' 638 : c print*,Erreurd1,x(3),y02,Y002 639 : 640 : ENDDO 641 : 642 : IF (ERREUR .GT. TEST ) THEN 643 : CONV=.FALSE. 644 : ELSE 645 : CONV=.TRUE. 646 : SIG1= X(1) 647 : SIG2= X(2) 648 : Y1 = X(3) 649 : Y2 = X(4) 650 : Z1 = X(5) 651 : Z2 = X(6) 652 : IF (SIG1 .GE. 0.D0) THEN 653 : JSIGN(1)=1 654 : ELSE 655 : JSIGN(1)=0 656 : END IF 657 : IF (SIG2.GE. 0.D0)THEN 658 : JSIGN(2)=1 659 : ELSE 660 : JSIGN(2)=0 661 : END IF 662 : JSIGN(3)=1 663 : END IF 664 : RETURN 665 : 200 format(6(6(2x,e11.4),/)) 666 : END 667 : 668 : 669 : 670 : 671 : 672 : 673 : 674 : 675 :
© Cast3M 2003 - All rights reserved.
Disclaimer