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