1 : C MODCLB SOURCE CHAT 05/01/13 01:47:35 5004 2 : C 3 : C 4 : C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5 : C % % 6 : C % INTEGRATION DU MODCLB UNILATERAL AVEC FERMETURE LINEAIRE % 7 : C % DE FISSURES EN CONTRAINTES PLANES % 8 : C % (MODCLB NON ASSOCIE) % 9 : C % A L"AIDE DE LA METHODE DE NEWTON % 10 : C % % 11 : C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 12 : C 13 : SUBROUTINE MODCLB(WRK0,NMATT,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,Y2, 14 : 1 JSIGN,SIG1,SIG2,CONV) 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 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 DELSIG 42 : INTEGER NDIM 43 : PARAMETER (NDIM=6) 44 : REAL*8 JACOB(NDIM,NDIM),J15,J16,J25,J26 45 : REAL*8 DELAFO(NDIM) 46 : REAL*8 X(NDIM),DX(NDIM),X0(NDIM),XX,XIG1,XIG2 47 : REAL*8 ERREUR,TEST 48 : REAL*8 H1,H2,HT,HT2,trf 49 : REAL*8 ALFA1,ALFA2 50 : REAL*8 SIGY 51 : REAL*8 UN,DEUX 52 : REAL*8 ZMAX 53 : INTEGER ITER ,I,DIMPB 54 : LOGICAL LDOM1,LDOM2 55 : SEGMENT WRK0 56 : REAL*8 XMAT(NMATT) 57 : ENDSEGMENT 58 : UN=1.D0 59 : DEUX=2.D0 60 : ZMAX=1.e+04 61 : E = XMAT(1) 62 : NU = XMAT(2) 63 : Y01 = XMAT(5) 64 : Y02 = XMAT(6) 65 : A1 = XMAT(7) 66 : A2 = XMAT(8) 67 : B1 = XMAT(9) 68 : B2 = XMAT(10) 69 : BETA1 = XMAT(11) 70 : BETA2 = XMAT(12) 71 : SIGF = XMAT(13) 72 : SIGY =DSQRT (2.D0*E*Y01) 73 : * DELSIG=SIGY*1.D-4 74 : DELSIG=0.d0 75 : * print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2' 76 : * print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2 77 : TEST=1.D-12 78 : ERREUR=1.D0 79 : ITER=0 80 : X(1)=SIG1 81 : X(2)=SIG2 82 : X(3)=Y1 83 : X(4)=Y2 84 : X(5)=Z1 85 : X(6)=Z2 86 : DO I=1, NDIM 87 : X0(I)=X(I) 88 : END DO 89 : C 90 : C 91 : C 92 : JACOB(1,1) = -UN 93 : JACOB(1,2) = 0.D0 94 : JACOB(1,3) = 0.D0 95 : JACOB(1,4) = 0.D0 96 : JACOB(2,1) = 0.D0 97 : JACOB(2,2) = -UN 98 : JACOB(2,3) = 0.D0 99 : JACOB(2,4) = 0.D0 100 : JACOB(3,3) = -UN 101 : JACOB(3,4) = 0.D0 102 : JACOB(3,6) = 0.D0 103 : JACOB(4,3) = 0.D0 104 : JACOB(4,4) = -UN 105 : JACOB(4,5) = 0.D0 106 : JACOB(5,1) = 0.D0 107 : JACOB(5,2) = 0.D0 108 : JACOB(5,4) = 0.D0 109 : JACOB(5,5) = -UN 110 : JACOB(5,6) = 0.D0 111 : JACOB(6,1) = 0.D0 112 : JACOB(6,2) = 0.D0 113 : JACOB(6,3) = 0.D0 114 : JACOB(6,5) = 0.D0 115 : JACOB(6,6) = -UN 116 : alfa1=UN 117 : alfa2=UN 118 : 119 : DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 50) 120 : ITER=ITER+1 121 : IF (X(1) .GE. DELSIG) THEN 122 : H1=UN 123 : ELSE IF (X(1) .LE. -DELSIG) THEN 124 : H1=0.D0 125 : ELSE 126 : H1=(X(1)+DELSIG)/(DEUX*DELSIG) 127 : END IF 128 : IF (X(2) .GE. DELSIG) THEN 129 : H2=UN 130 : ELSE IF (X(2) .LE. -DELSIG) THEN 131 : H2=0.D0 132 : ELSE 133 : H2=(X(2)+DELSIG)/(DEUX*DELSIG) 134 : END IF 135 : c 136 : c smooth transition on ht ht2 137 : c 138 : trf = x(1) + x(2) 139 : IF(trf.GE.delsig)THEN 140 : HT=UN 141 : HT2=0.D0 142 : else if(trf.ge.(-delsig)) then 143 : ht = (trf+delsig)/DEUX/delsig 144 : ht2 = UN-ht 145 : else if(trf.ge.(-sigf+delsig)) then 146 : ht = 0.D+00 147 : ht2 = UN 148 : else if(trf.ge.(-sigf-delsig)) then 149 : ht = 0.D+00 150 : c ht2 =UN- (trf-(sigf-delsig))/(DEUX*sigf) 151 : ht2 = (trf+sigf+delsig)/DEUX/delsig 152 : else 153 : ht = 0.D+00 154 : ht2 = 0.D+00 155 : endif 156 : LDOM1= (X(3) .GT. YLIM1 .AND. X(3).GT. Y01 157 : 1 .AND. X(5) .LT. ZMAX) 158 : LDOM2= (X(4) .GT. YLIM2 .AND. X(4).GT. Y02 159 : 1 .AND. X(6) .LT. ZMAX) 160 : * print*,'y1',x(3),'ylim1',ylim1 161 : * print*,'ldom1',ldom1,'ldom2',ldom2 162 : C =====> PROBLEME ELASTIQUE 163 : C 164 : C 165 : C CALCUL DU JACOBIEN 166 : C 167 : C 168 : JACOB(3,1) = H1*X(1)*X(5)*X(5)/E 169 : C 170 : JACOB(3,2) = H2*X(2)*X(5)*X(5)/E 171 : C 172 : JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E 173 : C 174 : JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E 175 : C 176 : C CALCUL DU SECOND MEMBRE 177 : C 178 : XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN) 179 : 1 *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2) 180 : 2 *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+ 181 : 3 BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1* 182 : 4 (X(6)*(H2-UN)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+ 183 : 5 (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)- 184 : 6 X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+ 185 : 7 (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)* 186 : 8 (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5)) 187 : 9 *BETA1*H1)*HT2) 188 : C 189 : XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN) 190 : 1 *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)* 191 : 2 X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+ 192 : 3 (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+ 193 : 4 (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))* 194 : 5 BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1) 195 : 6 *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2* 196 : 7 (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))* 197 : 8 BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)* 198 : 9 (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2) 199 : 200 : DELAFO(1) =X(1)-XIG1 201 : C 202 : DELAFO(2) =X(2)-XIG2 203 : C 204 : XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E) 205 : DELAFO(3) = X(3)-XX 206 : C 207 : XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) ) 208 : 1 *X(6)*X(6)/(DEUX*E) 209 : DELAFO(4) = X(4)-XX 210 : C 211 : IF (LDOM1) THEN 212 : C 213 : C =====> ENDOMMAGEMENT EN TRACTION 214 : C 215 : J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2 216 : & ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H 217 : & 2)-E*EPSP1*H2)*SIGF+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+ 218 : & H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6 219 : & )*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1 220 : & *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN 221 : & )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU 222 : & X*X(5)*H2)-X(6)*H2)*SIGF-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H 223 : & 2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5 224 : & )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+ 225 : & BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1* 226 : & (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*( 227 : & X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X 228 : & (5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS 229 : & P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2* 230 : & (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BE 231 : & TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX) 232 : & +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 233 : C 234 : C 235 : J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6 236 : & )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X( 237 : & 6)-UN)*BETA2-E*EPSP2)*H1)*SIGF+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT 238 : & 2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2 239 : & ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT 240 : & 2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX) 241 : & )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX 242 : & *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF-DEUX*BETA1*HT2*NU+(BETA1*(- 243 : & DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H 244 : & T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN) 245 : & *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X( 246 : & 6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN) 247 : & *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)* 248 : & X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF+(B 249 : & ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2 250 : & +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X( 251 : & 6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA 252 : & 1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U 253 : & N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 254 : C 255 : JACOB(1,5) = J15 256 : JACOB(2,5) = J25 257 : JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E 258 : JACOB(4,5) = 0.D0 259 : JACOB(5,5) = -1.D0 260 : JACOB(5,1) = 0.D0 261 : JACOB(5,2) = 0.D0 262 : JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN) 263 : JACOB(5,4) = 0.D0 264 : DELAFO(5) = X(5)-((A1*(X(3)-Y01))**B1+UN) 265 : IF (LDOM2) THEN 266 : C 267 : C =====> ENDOMMAGEMENT EN COMPRESSION 268 : C 269 : J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2) 270 : 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF/ 271 : 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)* 272 : 3 X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X 273 : 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5) 274 : 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5 275 : 6 )*H2+2*X(6)*(H2-1))*SIGF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE 276 : 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X 277 : 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5) 278 : 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+ 279 : : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6 280 : ; )*(H2-1)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE 281 : < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6 282 : = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1 283 : > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+( 284 : ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 285 : C 286 : J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B 287 : 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1- 288 : 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X 289 : 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2- 290 : 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5) 291 : 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X 292 : 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF+(BETA1* 293 : 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H 294 : 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6) 295 : 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6 296 : : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2 297 : ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS 298 : < P2+(X(6)**2-X(6))*BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5) 299 : = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2* 300 : > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5)) 301 : ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2) 302 : @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 303 : C 304 : C 305 : JACOB(1,6) = J16 306 : JACOB(2,6) = J26 307 : JACOB(3,6) = 0.D0 308 : JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2) 309 : & *X(2)*X(2))*X(6)/E 310 : JACOB(5,6) = 0.D0 311 : JACOB(6,6) = -1.D0 312 : JACOB(6,1) = 0.D0 313 : JACOB(6,2) = 0.D0 314 : JACOB(6,3) = 0.D0 315 : * JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0)) 316 : JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0) 317 : JACOB(6,5) = 0.D0 318 : * DELAFO(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0) 319 : DELAFO(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0) 320 : DIMPB=6 321 : ELSE 322 : DIMPB=5 323 : END IF 324 : ELSE 325 : C 326 : C =====> PAS D'ENDOMMAGEMENT EN TRACTION 327 : C 328 : IF (LDOM2) THEN 329 : C 330 : C =====> ENDOMMAGEMENT EN COMPRESSION 331 : C 332 : J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2) 333 : 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF/ 334 : 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)* 335 : 3 X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X 336 : 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5) 337 : 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5 338 : 6 )*H2+2*X(6)*(H2-1))*SIGF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE 339 : 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X 340 : 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5) 341 : 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+ 342 : : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6 343 : ; )*(H2-1)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE 344 : < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6 345 : = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1 346 : > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+( 347 : ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 348 : C 349 : J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B 350 : 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1- 351 : 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X 352 : 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2- 353 : 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5) 354 : 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X 355 : 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF+(BETA1* 356 : 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H 357 : 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6) 358 : 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6 359 : : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2 360 : ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS 361 : < P2+(X(6)**2-X(6))*BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5) 362 : = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2* 363 : > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5)) 364 : ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2) 365 : @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2 366 : C 367 : C 368 : 369 : JACOB(1,5) = J16 370 : JACOB(2,5) = J26 371 : JACOB(3,5) = 0.D0 372 : JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2) 373 : & *X(2)*X(2))*X(6)/E 374 : JACOB(5,5) = -1.D0 375 : JACOB(5,1) = 0.D0 376 : JACOB(5,2) = 0.D0 377 : JACOB(5,3) = 0.D0 378 : * JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0)) 379 : JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0) 380 : JACOB(5,5) = -1.D0 381 : * DELAFO(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0) 382 : DELAFO(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0) 383 : DIMPB=5 384 : ELSE 385 : DIMPB=4 386 : END IF 387 : END IF 388 : C 389 : C RESOLUTION 390 : C 390 : CALL GAUSCL(JACOB,DX,DELAFO,NDIM,DIMPB,0) 391 : C 392 : C CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION 393 : C 394 : ERREUR=0.D0 395 : DO I=1,2 396 : ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY) 397 : X(I)=X(I)+DX(I) 398 : END DO 399 : X(3)=X(3)+DX(3) 400 : ERREUR=ERREUR+DX(3)*DX(3)/(Y01*Y01) 401 : X(4)=X(4)+DX(4) 402 : ERREUR=ERREUR+DX(4)*DX(4)/(Y02*Y02) 403 : IF (LDOM1) THEN 404 : IF ((X(5)+DX(5)).LT.Z1) THEN 405 : DX(5)=Z1-X(5) 406 : END IF 407 : IF ( X(5) .GT. ZMAX)THEN 408 : DX(5)=ZMAX + 0.1D0 - X(5) 409 : END IF 410 : X(5)=X(5)+DX(5) 411 : ERREUR=ERREUR+DX(5)*DX(5)/(X(5)*X(5)) 412 : IF (LDOM2) THEN 413 : IF ((X(6)+DX(6)).LT.Z2) THEN 414 : DX(6)=Z2-X(6) 415 : END IF 416 : IF ( X(6) .GT. ZMAX)THEN 417 : DX(6)=ZMAX + 0.1D0 - X(6) 418 : END IF 419 : X(6)=X(6)+DX(6) 420 : ERREUR=ERREUR+DX(6)*DX(6)/(X(6)*X(6)) 421 : END IF 422 : ELSE IF (LDOM2) THEN 423 : IF ((X(6)+DX(5)).LT.Z2) THEN 424 : DX(5)=Z2-X(6) 425 : END IF 426 : IF ( X(6) .GT. ZMAX)THEN 427 : DX(5)=ZMAX + 0.1D0 - X(6) 428 : END IF 429 : X(6)=X(6)+DX(5) 430 : ERREUR=ERREUR+DX(5)*DX(5)/(X(6)*X(6)) 431 : END IF 432 : IF (ERREUR .GT. 1.D+100)THEN 433 : CONV=.FALSE. 434 : print*,'erreur dans modele',erreur 435 : print*,'iteration num',iter 436 : print*,'delafo',delafo 437 : RETURN 438 : END IF 439 : C 440 : C QUELQUES IMPRESSIONS EN CAS D'ERREUR 441 : C 442 : if (.not. conv )then 443 : PRINT*,'ITERATION ',ITER 444 : print*,'LDOM1 :',LDOM1,' LDOM2 :',LDOM2 445 : PRINT*,'EPS1= ',EPSP1,' EPS2= ',EPSP2 446 : PRINT*,'SIG1:',X(1) ,' DSIG1:',DX(1),'ERR :',DELAFO(1) 447 : PRINT*,'SIG2:',X(2) ,' DSIG2:',DX(2),'ERR :',DELAFO(2) 448 : PRINT*,'Y1 :',X(3) ,' DY1 :',DX(3),'ERR :',DELAFO(3) 449 : PRINT*,'Y2 :',X(4) ,' DY2 :',DX(4),'ERR :',DELAFO(4) 450 : PRINT*,'Z1 :',X(5) ,' DZ1 :',DX(5),'ERR :',DELAFO(5) 451 : PRINT*,'Z2 :',X(6) ,' DZ2 :',DX(6),'ERR :',DELAFO(6) 452 : PRINT*,'H1',H1,'H2',H2,ht,ht2 453 : PRINT* 454 : write(*,200)jacob 455 : 200 format(6(6(2x,e11.4),/)) 456 : END IF 457 : 458 : END DO 459 : IF (ERREUR .GT. TEST ) THEN 460 : CONV=.FALSE. 461 : ELSE 462 : CONV=.TRUE. 463 : SIG1= X(1) 464 : SIG2= X(2) 465 : Y1 = X(3) 466 : Y2 = X(4) 467 : Z1 = X(5) 468 : Z2 = X(6) 469 : IF (SIG1 .GE. 0.D0) THEN 470 : JSIGN(1)=1 471 : ELSE 472 : JSIGN(1)=0 473 : END IF 474 : IF (SIG2.GE. 0.D0)THEN 475 : JSIGN(2)=1 476 : ELSE 477 : JSIGN(2)=0 478 : END IF 479 : JSIGN(3)=1 480 : END IF 481 : RETURN 482 : END 483 : 484 :
© Cast3M 2003 - All rights reserved.
Disclaimer