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