1 : C BEHAV3 SOURCE CHAT 05/01/12 21:38:27 5004 2 : SUBROUTINE BEHAV3(SIGR,DSTRN,DEPS1,DEPS2,IPLA, 3 : & SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, 4 : & BETJEF,VISCO,NECH0,NECH1) 5 : C 6 : C ================================================================== 7 : C 8 : C MODELE DE PLASTICITE EN TRAITEMENT PRE ET POST PIC 9 : C Deux criteres de traction compression: 10 : C Rankine et Von Mises 11 : C 12 : C ================================================================== 13 : C CE SOUS-PROGRAMME EST APPELE DANS "BONE". 14 : C 15 : IMPLICIT INTEGER(I-N) 16 : IMPLICIT REAL*8(A-H,O-Z) 17 : DIMENSION DFSIG1(4),DFSIG2(4),SIGF(4),DSIGT(4),VEC1(4),VEC2(4) 18 : DIMENSION V1(4),AC1(4),AC2(4),SIGE(4),SIGP(4) 19 : DIMENSION DLAM1(2),DLAM0(2),SIGR(4),IS1(2) 20 : DIMENSION FCRI0(2),FCRI1(2),DEPSI(4),DSTRN(4),SIGRV(4) 21 : DIMENSION D2FSI1(4,4),D2FSI2(4,4),DEP(4,4),P1(4,4),P2(4,4) 22 : DIMENSION D(4,4),A(4,4),AI(4,4),DJAC0(4,4),DJI(4,4),AH(4,4) 23 : C 24 : * 25 : SEGMENT BETJEF 26 : REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF, 27 : & TCON,DPSTF1,DPSTF2,TETA,PDT,TP0 28 : INTEGER ICT,ICC,IMOD,IVIS,ITR, 29 : & ISIM,IBB,IGAU,IZON 30 : ENDSEGMENT 31 : SEGMENT VISCO 32 : REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV 33 : ENDSEGMENT 34 : SEGMENT NECH0 35 : REAL*8 DT,DC,ALFG,S0,ENDO 36 : ENDSEGMENT 37 : SEGMENT NECH1 38 : REAL*8 ENDL 39 : ENDSEGMENT 40 : * 41 : * 42 : * COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF, 43 : * & TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR, 44 : * & ISIM,IBB,IGAU,IZON 45 : C 46 : MNSTRS=NSTRS 47 : C 48 : IAPEX=0 49 : Ft=PALF*COLI 50 : Rb=COLI 51 : PRB=1.D-10 52 : PRB2=1.D-6 53 : ITER=1 54 : IEC1=0 55 : IEC2=0 56 : IBROY=0 57 : ITANG=0 58 : CRIMAX=0.D0 59 : SEQ = 0.D0 60 : SEQ1 = 0.D0 61 : SEQ2 = 0.D0 61 : CALL ZERO(SIGE,4,1) 61 : CALL ZERO(V1,4,1) 61 : CALL ZERO(P1,4,4) 61 : CALL ZERO(P2,4,4) 61 : CALL ZERO(D,4,4) 62 : C 63 : DO 10 I=1,NSTRS 64 : SIGE(I)=SIGF(I) 65 : 10 CONTINUE 66 : C 67 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 68 : AD=YOUN/(1.D0-XNU*XNU) 69 : D(1,1)=AD 70 : D(2,2)=D(1,1) 71 : D(3,3)=AD*(1.D0-XNU)/2.D0 72 : D(1,2)=AD*XNU 73 : D(2,1)=D(1,2) 74 : ENDIF 75 : C 76 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 77 : ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU)) 78 : D(1,1)=ADD*(1.D0-XNU) 79 : D(2,2)=D(1,1) 80 : D(3,3)=D(1,1) 81 : D(1,2)=ADD*XNU 82 : D(2,1)=D(1,2) 83 : D(1,3)=D(1,2) 84 : D(2,3)=D(1,2) 85 : D(3,1)=D(1,2) 86 : D(3,2)=D(1,2) 87 : D(4,4)=0.5*ADD*(1.D0-2.D0*XNU) 88 : ENDIF 89 : C 90 : C ************ Le point est fissure pour 1ere fois *********** 91 : C 92 : IF (IFIS.EQ.0) THEN 92 : CALL PRINC(SIGF,V1,NSTRS) 93 : DEPS1=0.D0 94 : SIG1=Ft 95 : IF (V1(1).GT.Ft) THEN 96 : IFIS=1 97 : ENDIF 98 : ENDIF 99 : IF (IPLA.EQ.0) THEN 100 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 100 : CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF) 101 : ENDIF 102 : IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN 102 : CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF) 103 : ENDIF 104 : DEPS2=0.D0 105 : SIG2=Rb*AA 106 : IF (SEQ.GT.SIG2) THEN 107 : IPLA=1 108 : ENDIF 109 : ENDIF 110 : C 111 : 12 CONTINUE 112 : C 112 : CALL PRINC(SIGE,V1,NSTRS) 113 : TETA=V1(4) 114 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 114 : CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF) 115 : ENDIF 116 : IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN 116 : CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF) 117 : ENDIF 118 : FCRI0(1)=V1(1)-SIG1 119 : FCRI0(2)=SEQ-SIG2 120 : IF (DABS(FCRI0(1)).GT.DABS(FCRI0(2))) THEN 121 : CRIMAX=DABS(100.D0*FCRI0(1)) 122 : ELSE 123 : CRIMAX=DABS(100.D0*FCRI0(2)) 124 : ENDIF 125 : IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN 126 : IEC1=1 127 : IEC2=1 128 : ENDIF 129 : IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN 130 : IEC1=1 131 : IEC2=0 132 : ENDIF 133 : IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN 134 : IEC1=0 135 : IEC2=1 136 : ENDIF 137 : C 138 : 15 CONTINUE 139 : C 140 : DO 16 I=1,NSTRS 141 : SIGF(I)=SIGE(I) 142 : 16 CONTINUE 143 : C 144 : IF (IEC1.EQ.1.AND.IEC2.EQ.1) THEN 145 : DO 17 I=1,NSTRS 146 : DO 17 J=1,NSTRS 147 : DEP(I,J)=D(I,J) 148 : 17 CONTINUE 149 : GOTO 100 150 : ENDIF 151 : C 152 : IF (IEC1.EQ.1.AND.IEC2.EQ.0) THEN 152 : CALL BEHAV2(SIGR,DSTRN,DEPS1,DEPS2,IPLA 153 : & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, 154 : & BETJEF,VISCO,NECH0,NECH1) 155 : GOTO 100 156 : ENDIF 157 : C 158 : IF (IEC1.EQ.0.AND.IEC2.EQ.1) THEN 158 : CALL BEHAV1(SIGR,DSTRN,DEPS1,DEPS2,IPLA 159 : & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP, 160 : & BETJEF,VISCO,NECH0,NECH1) 161 : GOTO 100 162 : ENDIF 163 : C 164 : DK1=DEPS1 165 : DK2=DEPS2 166 : DLAM0(1)=0.D0 167 : DLAM0(2)=0.D0 168 : C 169 : 18 CONTINUE 170 : C 171 : C ************ Preparation des criteres ******************** 172 : C 173 : C ------------------ Von Mises ----------------------------- 174 : IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN 175 : A1=0.D0 176 : A2=1.D0 177 : ENDIF 178 : C ----------------- Drucker Prager ------------------------- 179 : IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN 180 : A1=(BETA-1.D0)/(2*BETA-1.D0) 181 : A2=BETA/(2*BETA-1.D0) 182 : ENDIF 183 : C ---------------------------------------------------------- 184 : C 185 : C ---------------- direction de traction ------------------ 186 : C 187 : C--------------- Determination du parametre d'ecrouissage ---------- 188 : C 189 : IF(IVIS.LE.2) THEN 189 : CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF) 190 : ELSE 190 : CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1) 191 : ENDIF 192 : C 193 : C------------------------------------------------------------ 194 : PI=4.D0*DATAN(1.D0) 195 : PHIC=V1(4)*(PI/180.D0) 196 : COSA=DCOS(PHIC) 197 : SINA=DSIN(PHIC) 198 : C 199 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 200 : DFSIG1(1)=COSA*COSA 201 : DFSIG1(2)=SINA*SINA 202 : DFSIG1(3)=2.D0*SINA*COSA 203 : DO 19 I=1,NSTRS 204 : DO 19 J=1,NSTRS 205 : P1(I,J)=0.D0 206 : 19 CONTINUE 207 : P1(1,1)=1.D0/2.D0 208 : P1(1,2)=-1.D0/2.D0 209 : P1(2,1)=-1.D0/2.D0 210 : P1(2,2)=1.D0/2.D0 211 : P1(3,3)=2.D0 212 : ENDIF 213 : C 214 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 215 : DFSIG1(1)=COSA*COSA 216 : DFSIG1(2)=SINA*SINA 217 : DFSIG1(3)=0.D0 218 : DFSIG1(4)=2.D0*SINA*COSA 219 : DO 21 I=1,NSTRS 220 : DO 21 J=1,NSTRS 221 : P1(I,J)=0.D0 222 : 21 CONTINUE 223 : P1(1,1)=1.D0/2.D0 224 : P1(1,2)=-1.D0/2.D0 225 : P1(2,1)=-1.D0/2.D0 226 : P1(2,2)=1.D0/2.D0 227 : P1(4,4)=2.D0 228 : ENDIF 229 : C 230 : DO 20 I=1,NSTRS 231 : AC1(I)=0.D0 232 : DO 20 J=1,NSTRS 233 : AC1(I)=AC1(I)+D(I,J)*DFSIG1(J) 234 : 20 CONTINUE 235 : C 236 : C ---------------- direction de compression ---------------- 237 : C 238 : IF(IVIS.LE.2) THEN 238 : CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF) 239 : ELSE 239 : CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) 240 : ENDIF 241 : C 242 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 243 : SX=SIGF(1) 244 : SY=SIGF(2) 245 : SXY=SIGF(3) 246 : R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY 247 : R=SQRT(0.5*R1) 248 : IF (DABS(R).LE.1.D-10) THEN 249 : WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3' 250 : STOP 251 : ENDIF 252 : DFSIG2(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2 253 : DFSIG2(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2 254 : DFSIG2(3)=3.D0*SXY/R/A2 255 : P2(1,1)=2.D0/A2/A2 256 : P2(1,2)=-1.D0/A2/A2 257 : P2(1,3)=0.D0/A2/A2 258 : P2(2,1)=-1.D0/A2/A2 259 : P2(2,2)=2.D0/A2/A2 260 : P2(2,3)=0.D0/A2/A2 261 : P2(3,1)=0.D0/A2/A2 262 : P2(3,2)=0.D0/A2/A2 263 : P2(3,3)=6.D0/A2/A2 264 : ENDIF 265 : C 266 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 267 : SX=SIGF(1) 268 : SY=SIGF(2) 269 : SZ=SIGF(3) 270 : SXY=SIGF(4) 271 : TX=(SX-SY)*(SX-SY) 272 : TY=(SX-SZ)*(SX-SZ) 273 : TZ=(SY-SZ)*(SY-SZ) 274 : R1=TX+TY+TZ+6*SXY*SXY 275 : R=SQRT(0.5*R1) 276 : IF (DABS(R).LE.1.D-10) THEN 277 : WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3' 278 : STOP 279 : ENDIF 280 : DFSIG2(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2 281 : DFSIG2(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2 282 : DFSIG2(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2 283 : DFSIG2(4)=3.D0*SXY/R/A2 284 : P2(1,1)=2.D0/A2/A2 285 : P2(1,2)=-1.D0/A2/A2 286 : P2(1,3)=-1.D0/A2/A2 287 : P2(1,4)=0.D0/A2/A2 288 : P2(2,1)=-1.D0/A2/A2 289 : P2(2,2)=2.D0/A2/A2 290 : P2(2,3)=-1.D0/A2/A2 291 : P2(2,4)=0.D0/A2/A2 292 : P2(3,1)=-1.D0/A2/A2 293 : P2(3,2)=-1.D0/A2/A2 294 : P2(3,3)=2.D0/A2/A2 295 : P2(3,4)=0.D0/A2/A2 296 : P2(4,1)=0.D0/A2/A2 297 : P2(4,2)=0.D0/A2/A2 298 : P2(4,3)=0.D0/A2/A2 299 : P2(4,4)=6.D0/A2/A2 300 : ENDIF 301 : C 302 : DO 25 I=1,NSTRS 303 : AC2(I)=0.D0 304 : DO 25 J=1,NSTRS 305 : AC2(I)=AC2(I)+D(I,J)*DFSIG2(J) 306 : 25 CONTINUE 307 : C 308 : C ************ Preparation du jacobien ******************** 309 : C 310 : F1DF1=0.D0 311 : DO 30 J=1,NSTRS 312 : F1DF1=F1DF1+DFSIG1(J)*AC1(J) 313 : 30 CONTINUE 314 : F1DF2=0.D0 315 : DO 31 J=1,NSTRS 316 : F1DF2=F1DF2+DFSIG1(J)*AC2(J) 317 : 31 CONTINUE 318 : F2DF2=0.D0 319 : DO 32 J=1,NSTRS 320 : F2DF2=F2DF2+DFSIG2(J)*AC2(J) 321 : 32 CONTINUE 322 : F2DF1=0.D0 323 : DO 33 J=1,NSTRS 324 : F2DF1=F2DF1+DFSIG2(J)*AC1(J) 325 : 33 CONTINUE 326 : C 327 : DJAC0(1,1)=-(F1DF1+PAECT) 328 : DJAC0(1,2)=-(F1DF2) 329 : DJAC0(2,2)=-(F2DF2+PAEC) 330 : DJAC0(2,1)=-(F2DF1) 331 : C 332 : DO 43 I=1,2 333 : DO 43 J=1,2 334 : DJI(I,J)=DJAC0(I,J) 335 : 43 CONTINUE 336 : C 336 : CALL INVMA2(DJI,2,ISING) 337 : IF (ISING.EQ.1) THEN 338 : WRITE(*,*)'MATRICE DJI singuliere ds BEHAV3' 339 : ENDIF 340 : C 341 : C ************ Debut des iterations internes ******************* 342 : C 343 : 40 CONTINUE 344 : C 345 : C *************** Determination de DK et DLAM ****************** 346 : C 347 : C 348 : C 349 : DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2) 350 : DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2) 351 : C 352 : IF (DLAM1(1).LE.0.D0) IEC1=1 353 : IF (DLAM1(2).LE.0.D0) IEC2=1 354 : IF (IEC1.EQ.1.OR.IEC2.EQ.1) THEN 355 : C WRITE(*,*)'Dans BEHAV3, DLAMDA1 est negatif:',DLAM1(1) 356 : C WRITE(*,*)'Dans BEHAV3, DLAMDA2 est negatif:',DLAM1(2) 357 : C WRITE(*,*)'A l iteration :',ITER 358 : GOTO 15 359 : ENDIF 360 : C 361 : DK1=DEPS1+DLAM1(1) 362 : DK2=DEPS2+DLAM1(2) 363 : C 364 : IF(IVIS.LE.2) THEN 364 : CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF) 364 : CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF) 365 : ELSE 365 : CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1) 365 : CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) 366 : ENDIF 367 : C 368 : C ************** Determination de DPHI ********************* 369 : C 370 : IF (IMOD.EQ.1) THEN 371 : AD1=YOUN/(1.D0-XNU) 372 : coef=2.D0*SEQ2/(2.D0*SEQ2+DLAM1(2)*AD1) 373 : DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1)) 374 : TO=SIGF(3) 375 : DPHI2=SEQ2 376 : ENDIF 377 : IF (IMOD.EQ.2) THEN 378 : DG=YOUN/2.D0/(1.D0+XNU) 379 : AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU) 380 : coef1=(SEQ2+2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG) 381 : coef2=(2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG) 382 : CAS1=(SIGE(1)+SIGE(2)-0.5*DLAM1(1)*2.D0*AD1) 383 : CAS2=SIGE(3)-DLAM1(1)*XNU*AD1 384 : DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2) 385 : TO=SIGF(4) 386 : DPHI2=SEQ2 387 : ENDIF 388 : IF (IMOD.EQ.3) THEN 389 : DPHI2=SEQ2-A1/A2*(SIGF(1)+SIGF(2)) 390 : AD1=YOUN/(1.D0-XNU) 391 : coef=2.D0*DPHI2/(2.D0*DPHI2+DLAM1(2)*AD1) 392 : SOM=A1/A2*DLAM1(2)*2.D0*AD1 393 : DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1)-SOM) 394 : TO=SIGF(3) 395 : ENDIF 396 : IF (IMOD.EQ.4) THEN 397 : SOM=SIGE(1)+SIGE(2)+SIGE(3) 398 : DPHI2=A1*SOM/A2 399 : DPHI2=DPHI2-A1*DLAM1(1)*ADD*(1.D0+XNU)/A2 400 : DPHI2=DPHI2-A1*A1*DLAM1(2)*3.D0*ADD*(1.D0+XNU)/A2/A2 401 : DPHI2=SEQ2-DPHI2 402 : DG=YOUN/2.D0/(1.D0+XNU)/A2/A2 403 : AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU) 404 : AD2=YOUN/(1.D0-2.D0*XNU) 405 : coef1=(DPHI2+2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG) 406 : coef2=(2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG) 407 : CAS1=(SIGE(1)+SIGE(2)-DLAM1(1)*AD1-2.D0*A1/A2*DLAM1(2)*AD2) 408 : CAS2=SIGE(3)-DLAM1(1)*XNU*AD1-A1/A2*DLAM1(2)*AD2 409 : DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2) 410 : TO=SIGF(4) 411 : ENDIF 412 : C IF (DPHI1.LT.0.D0) THEN 413 : C WRITE(*,*)'ATTENTION DPHI1 NEGATIF' 414 : C WRITE(*,*)'DPHI1=',DPHI1 415 : C ENDIF 416 : C IF (DPHI2.LT.0.D0) THEN 417 : C WRITE(*,*)'ATTENTION DPHI2 ds BEHAV3 NEGATIF' 418 : C WRITE(*,*)'DPHI2=',DPHI2 419 : C ENDIF 420 : C 421 : C--------------- Cas de l'apex ----------------------------- 422 : C 423 : IF (DABS(DPHI1).LE.10E-10.AND. 424 : *DABS(DPHI2).GT.10E-10) THEN 425 : IAPEX=1 426 : C WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX 427 : C WRITE(*,*)'Dans l element',IBB 428 : C WRITE(*,*)'et au point d intégration',IGAU 429 : DO 50 I=1,NSTRS 430 : DO 50 J=1,NSTRS 431 : AI(I,J)=0.D0 432 : 50 CONTINUE 433 : AI(1,1)=0.5 434 : AI(1,2)=0.5 435 : AI(2,1)=AI(1,2) 436 : AI(2,2)=AI(1,1) 437 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0 438 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0 439 : GOTO 75 440 : ENDIF 441 : C 442 : C---------- Cas du critere reduit a un point ---------------- 443 : C 444 : IF (DABS(DPHI2).LE.10E-10) THEN 445 : IAPEX=2 446 : C WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX 447 : C WRITE(*,*)'Dans l element',IBB 448 : C WRITE(*,*)'et au point d intégration',IGAU 449 : DO 52 I=1,NSTRS 450 : DO 52 J=1,NSTRS 451 : AI(I,J)=0.D0 452 : 52 CONTINUE 453 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 454 : DO 53 I=1,3 455 : DO 53 J=1,3 456 : AI(I,J)=1.D0/3.D0 457 : 53 CONTINUE 458 : ENDIF 459 : GOTO 75 460 : ENDIF 461 : C 462 : C ************** Mise a jour des contraintes *************** 463 : C 464 : C ---------------- calcul de la matrice A ------------------ 465 : C 466 : DO 60 I=1,NSTRS 467 : DO 60 J=1,NSTRS 468 : A(I,J)=0.D0 469 : 60 CONTINUE 470 : C 471 : DG=YOUN/2.D0/(1.D0+XNU) 472 : C 473 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 474 : DW1=AD*(2.D0-XNU) 475 : DW2=AD*(2.D0*XNU-1.D0) 476 : A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1 477 : A(1,1)=A(1,1)+(DLAM1(2)*DW1)/2.D0/DPHI2/A2/A2 478 : A(2,2)=A(1,1) 479 : A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1 480 : A(1,2)=A(1,2)+(DLAM1(2)*DW2)/2.D0/DPHI2/A2/A2 481 : A(2,1)=A(1,2) 482 : A(3,3)=(DLAM1(1)*DG)/DPHI1 483 : A(3,3)=A(3,3)+1.D0+(DLAM1(2)*3.D0*DG)/DPHI2/A2/A2 484 : ENDIF 485 : C 486 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 487 : A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1 488 : A(1,1)=A(1,1)+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2 489 : A(2,2)=A(1,1) 490 : A(3,3)=1.D0+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2 491 : A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1 492 : A(1,2)=A(1,2)-(DLAM1(2)*DG)/DPHI2/A2/A2 493 : A(2,1)=A(1,2) 494 : A(1,3)=-(DLAM1(2)*DG)/DPHI2/A2/A2 495 : A(3,1)=A(1,3) 496 : A(3,2)=A(1,3) 497 : A(2,3)=A(1,3) 498 : A(4,4)=1.D0+(DLAM1(1)*DG)/DPHI1 499 : A(4,4)=A(4,4)+(DLAM1(2)*DG*3.D0)/DPHI2/A2/A2 500 : ENDIF 501 : C 502 : C -------------- invertion de la matrice A ----------------- 503 : C 504 : DO 70 I=1,NSTRS 505 : DO 70 J=1,NSTRS 506 : AI(I,J)=A(I,J) 507 : 70 CONTINUE 508 : C 508 : CALL INVMA2(AI,NSTRS,ISING) 509 : IF (ISING.EQ.1) THEN 510 : WRITE(*,*)'MATRICE AI singuliere ds BEHAV3' 511 : ENDIF 512 : C 513 : C -------------- mise a jour des contraintes ------------ 514 : C 515 : 75 CONTINUE 516 : VEC1(1)=AD1 517 : VEC1(2)=AD1 518 : VEC1(3)=0.D0 519 : VEC1(4)=0.D0 520 : IF(IMOD.EQ.3) THEN 521 : VEC2(1)=AD1*(1.D0+XNU) 522 : VEC2(2)=AD1*(1.D0+XNU) 523 : VEC2(3)=0.D0 524 : VEC2(4)=0.D0 525 : ENDIF 526 : IF(IMOD.EQ.4) THEN 527 : VEC2(1)=AD1*(1.D0+XNU) 528 : VEC2(2)=AD1*(1.D0+XNU) 529 : VEC2(3)=AD1*(1.D0+XNU) 530 : VEC2(4)=0.D0 531 : ENDIF 532 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 533 : VEC2(1)=0.D0 534 : VEC2(2)=0.D0 535 : VEC2(3)=0.D0 536 : VEC2(4)=0.D0 537 : ENDIF 538 : C 539 : DO 80 I=1,NSTRS 540 : DEPSI(I)=SIGE(I)-0.5*DLAM1(1)*VEC1(I) 541 : DEPSI(I)=DEPSI(I)-A1*DLAM1(2)*VEC2(I)/A2 542 : 80 CONTINUE 543 : C 544 : DO 90 I=1,NSTRS 545 : SIGF(I)=0.0D+00 546 : DO 90 J=1,NSTRS 547 : SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 548 : 90 CONTINUE 549 : C 550 : C ******** Verification des criteres **************** 551 : C 551 : CALL PRINC(SIGF,V1,NSTRS) 552 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 552 : CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF) 553 : ENDIF 554 : IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN 554 : CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF) 555 : ENDIF 556 : FCRI1(1)=V1(1)-SEQ1 557 : IF (IAPEX.EQ.2) FCRI1(1)=0.D0 558 : FCRI1(2)=SEQ-SEQ2 559 : C 560 : IF (IBROY.EQ.0.AND.(DABS(FCRI1(1)).GE.CRIMAX.OR 561 : * .DABS(FCRI1(2)).GE.CRIMAX)) THEN 562 : C WRITE(*,*)'****************************************' 563 : C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN' 564 : C WRITE(*,*)'on passe donc a la secante' 565 : C WRITE(*,*)'Dans l element',IBB 566 : C WRITE(*,*)'et au point d intégration',IGAU 567 : C WRITE(*,*)'CRIMAX=',CRIMAX 568 : C WRITE(*,*)'****************************************' 569 : ITER=ITR 570 : ENDIF 571 : C 572 : C ******* Compteur sur la methode de resolution **** 573 : C 574 : IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN 575 : IBROY=1 576 : ITANG=1 577 : ITER=1 578 : IAPEX=0 579 : IEC1=0 580 : IEC2=0 581 : GOTO 12 582 : ENDIF 583 : C 584 : C ******* non convergence ************************** 585 : C 586 : IF ((DABS(FCRI1(1)).GT.PRB.OR.DABS(FCRI1(2)).GT.PRB) 587 : *.AND.ITER.LT.ITR) THEN 588 : IF (IBROY.EQ.0) THEN 588 : CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1) 589 : DLAM0(1)=DLAM1(1) 590 : DLAM0(2)=DLAM1(2) 591 : FCRI0(1)=FCRI1(1) 592 : FCRI0(2)=FCRI1(2) 593 : IF (ITER.GE.(ITR-1)) THEN 594 : C WRITE(*,*)'***********************' 595 : C WRITE(*,*)'BROYDEN n a pas aboutit' 596 : C WRITE(*,*)'ITER=',ITER 597 : C WRITE(*,*)'FCRIT=',FCRI0(1) 598 : C WRITE(*,*)'FCRIC=',FCRI0(2) 599 : C WRITE(*,*)'Dans l element',IBB 600 : C WRITE(*,*)'et au point d intégration',IGAU 601 : C WRITE(*,*)'***********************' 602 : ENDIF 603 : ITER=ITER+1 604 : GOTO 40 605 : ENDIF 606 : IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN 607 : DLAM0(1)=DLAM1(1) 608 : DLAM0(2)=DLAM1(2) 609 : FCRI0(1)=FCRI1(1) 610 : FCRI0(2)=FCRI1(2) 610 : CALL PRINC(SIGF,V1,NSTRS) 611 : IF (ITER.GE.(ITR-5)) THEN 612 : C WRITE(*,*)'ITER=',ITER 613 : C WRITE(*,*)'FCRIT=',FCRI0(1) 614 : C WRITE(*,*)'FCRIC=',FCRI0(2) 615 : ENDIF 616 : ITER=ITER+1 617 : GOTO 18 618 : ENDIF 619 : ENDIF 620 : IF (ITER.GE.ITR.AND.(DABS(FCRI1(1)).GT.PRB2.OR. 621 : * DABS(FCRI1(2)).GT.PRB2)) THEN 622 : WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3' 623 : WRITE(*,*)'Dans l element',IBB 624 : WRITE(*,*)'et au point d intégration',IGAU 625 : WRITE(*,*)'FCRIT=',FCRI0(1) 626 : WRITE(*,*)'FCRIC=',FCRI0(2) 627 : WRITE(*,*)'IPLA=',IPLA 628 : WRITE(*,*)'IFIS=',IFIS 629 : C STOP 630 : ENDIF 631 : C 632 : C **************** Fin des iterations internes ******************* 633 : C 634 : C 635 : C ************* Calcul de la DEP **************** 636 : C 637 : IF (IAPEX.EQ.0) THEN 637 : CALL DERI2(D2FSI1,DPHI1,P1,SIGF,NSTRS,BETJEF) 637 : CALL DERI2(D2FSI2,DPHI2,P2,SIGF,NSTRS,BETJEF) 637 : CALL LADEP1(SIGF,DEP,PAECT,PAEC,NSTRS,IFOU, 638 : & DFSIG1,DFSIG2,D2FSI1,D2FSI2,DLAM1,BETJEF) 639 : ELSE IF (IAPEX.EQ.1) THEN 640 : DO 95 I=1,NSTRS 641 : DO 95 J=1,NSTRS 642 : AH(I,J)=0.D0 643 : DO 95 K=1,NSTRS 644 : AH(I,J)=AH(I,J)+AI(I,K)*D(K,J) 645 : 95 CONTINUE 646 : DFSIG1(1)=0.5*DSQRT(2.D0) 647 : DFSIG1(2)=0.5*DSQRT(2.D0) 648 : DFSIG1(3)=0.D0 649 : IF (IMOD.EQ.2) DFSIG1(4)=0.D0 649 : CALL CREDEP(AH,DFSIG1,PAEC,NSTRS,DEP,BETJEF) 650 : ELSE IF (IAPEX.EQ.2) THEN 651 : DO 96 I=1,NSTRS 652 : DO 96 J=1,NSTRS 653 : DEP(I,J)=0.D0 654 : DO 96 K=1,NSTRS 655 : DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J) 656 : 96 CONTINUE 657 : ENDIF 658 : IAPEX=0 659 : C DO 97 I=1,NSTRS 660 : C DO 97 J=1,NSTRS 661 : C DEP(I,J)=D(I,J) 662 : C 97 CONTINUE 663 : C 664 : C *********************************************** 665 : C 666 : DEPS1=DK1 667 : DEPS2=DK2 668 : SIG1=SEQ1 669 : SIG2=SEQ2 670 : IF(IVIS.LE.2) THEN 670 : CALL INDICA(DK1,DK2,IFIS,IPLA,3,BETJEF) 671 : ELSE 671 : CALL INDIC1(DK1,DK2,IFIS,IPLA,3,BETJEF,NECH0,NECH1) 672 : ENDIF 673 : C ********************** CALCUL VISCOPLASTIQUE *********************** 674 : IF (IVIS.EQ.1) THEN 674 : CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2, 675 : & SIGP,SIGRV,DSTRN,D,BETJEF,VISCO) 676 : ENDIF 677 : C ******************************************************************** 678 : 100 CONTINUE 679 : IEC1=0 680 : IEC2=0 681 : CRIMAX=0.D0 682 : C 683 : RETURN 684 : END 685 : 686 : 687 : 688 :
© Cast3M 2003 - All rights reserved.
Disclaimer