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