1 : C BEHAV2 SOURCE PV 07/11/23 21:15:30 5978 2 : SUBROUTINE BEHAV2(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 Un seul critere de compression: Von Mises ou Drucker-Prager 10 : C 11 : C ================================================================== 12 : C CE SOUS-PROGRAMME EST APPELE DANS "BONE". 13 : C 14 : IMPLICIT INTEGER(I-N) 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : DIMENSION DFSIG(4),SIGF(4),DSIGT(4),VEC1(4),DEPSI(4),DSTRN(4) 17 : DIMENSION SIGR(4),AC(4),SIGE(4),SIGRV(4),SIGP(4) 18 : DIMENSION D(4,4),DEP(4,4),D2FSIG(4,4),P2(4,4),DP(4,4) 19 : DIMENSION A(4,4),AI(4,4),AH(4,4) 20 : * 21 : SEGMENT BETJEF 22 : REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF, 23 : & TCON,DPSTF1,DPSTF2,TETA,PDT,TP0 24 : INTEGER ICT,ICC,IMOD,IVIS,ITR, 25 : & ISIM,IBB,IGAU,IZON 26 : ENDSEGMENT 27 : * 28 : SEGMENT VISCO 29 : REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV 30 : ENDSEGMENT 31 : SEGMENT NECH0 32 : REAL*8 DT,DC,ALFG,S0,ENDO 33 : ENDSEGMENT 34 : SEGMENT NECH1 35 : REAL*8 ENDL 36 : ENDSEGMENT 37 : * 38 : * 39 : C 40 : * COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF, 41 : * & TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR, 42 : * & ISIM,IBB,IGAU,IZON 43 : C 44 : IAPEX=0 45 : PRB=1.D-10 46 : PRB2=1.D-6 47 : ITER=1 48 : IBROY=0 49 : ITANG=0 50 : Rb=COLI 51 : Ft=PALF*COLI 52 : IEC2=0 53 : CRIMAX=0.D0 54 : SEQ = 0.D0 55 : SEQ2 = 0.D0 55 : CALL ZERO(SIGE,4,1) 55 : ***** CALL ZERO(V1,4,1) 55 : CALL ZERO(P2,4,4) 55 : CALL ZERO(D,4,4) 56 : C 57 : DO 10 I=1,NSTRS 58 : SIGE(I)=SIGF(I) 59 : 10 CONTINUE 60 : C 61 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 62 : AD=YOUN/(1.D0-XNU*XNU) 63 : D(1,1)=AD 64 : D(2,2)=D(1,1) 65 : D(3,3)=AD*(1.D0-XNU)/2.D0 66 : D(1,2)=AD*XNU 67 : D(2,1)=D(1,2) 68 : ENDIF 69 : C 70 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 71 : ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU)) 72 : D(1,1)=ADD*(1.D0-XNU) 73 : D(2,2)=D(1,1) 74 : D(3,3)=D(1,1) 75 : D(1,2)=ADD*XNU 76 : D(2,1)=D(1,2) 77 : D(1,3)=D(1,2) 78 : D(2,3)=D(1,2) 79 : D(3,1)=D(1,2) 80 : D(3,2)=D(1,2) 81 : D(4,4)=0.5*ADD*(1.D0-2.D0*XNU) 82 : ENDIF 83 : C 84 : C ************ Le point est comprime pour 1ere fois *********** 85 : C 86 : IF (IPLA.EQ.0) THEN 87 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 87 : CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF) 88 : ENDIF 89 : IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN 89 : CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF) 90 : ENDIF 91 : DEPS2=0.D0 92 : SIG2=Rb*AA 93 : IF (SEQ.GT.SIG2) THEN 94 : IPLA=1 95 : ENDIF 96 : ENDIF 97 : IF (IFIS.EQ.0) THEN 98 : DEPS1=0.D0 99 : SIG1=Ft 100 : ENDIF 101 : C 102 : 8 CONTINUE 103 : C 104 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 104 : CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF) 105 : ENDIF 106 : IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN 106 : CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF) 107 : ENDIF 108 : FCRIC0=SEQ-SIG2 109 : CRIMAX=DABS(100.D0*FCRIC0) 110 : IF (FCRIC0.LE.0.D0) THEN 111 : IEC2=1 112 : ENDIF 113 : DK=DEPS2 114 : DLAM0=0.D0 115 : C 116 : 12 CONTINUE 117 : C 118 : DO 15 I=1,NSTRS 119 : SIGF(I)=SIGE(I) 120 : 15 CONTINUE 121 : C 122 : IF (IEC2.EQ.1) THEN 123 : DO 16 I=1,NSTRS 124 : DO 16 J=1,NSTRS 125 : DEP(I,J)=D(I,J) 126 : 16 CONTINUE 127 : GOTO 100 128 : ENDIF 129 : C 130 : 9 CONTINUE 131 : C 132 : IF(IVIS.LE.2) THEN 132 : CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF) 133 : ELSE 133 : CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) 134 : ENDIF 135 : C 136 : C ************ Preparation des criteres ******************** 137 : C 138 : C ------------------ Von Mises ----------------------------- 139 : IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN 140 : A1=0.D0 141 : A2=1.D0 142 : ENDIF 143 : C ----------------- Drucker Prager ------------------------- 144 : IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN 145 : A1=(BETA-1.D0)/(2*BETA-1.D0) 146 : A2=BETA/(2*BETA-1.D0) 147 : ENDIF 148 : C ---------------------------------------------------------- 149 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 150 : SX=SIGF(1) 151 : SY=SIGF(2) 152 : SXY=SIGF(3) 153 : R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY 154 : R=SQRT(0.5*R1) 155 : IF (DABS(R).LE.1.D-10) THEN 156 : WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2' 157 : STOP 158 : ENDIF 159 : DFSIG(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2 160 : DFSIG(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2 161 : DFSIG(3)=3.D0*SXY/R/A2 162 : P2(1,1)=2.D0/A2/A2 163 : P2(1,2)=-1.D0/A2/A2 164 : P2(1,3)=0.D0/A2/A2 165 : P2(2,1)=-1.D0/A2/A2 166 : P2(2,2)=2.D0/A2/A2 167 : P2(2,3)=0.D0/A2/A2 168 : P2(3,1)=0.D0/A2/A2 169 : P2(3,2)=0.D0/A2/A2 170 : P2(3,3)=6.D0/A2/A2 171 : ENDIF 172 : C 173 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 174 : SX=SIGF(1) 175 : SY=SIGF(2) 176 : SZ=SIGF(3) 177 : SXY=SIGF(4) 178 : TX=(SX-SY)*(SX-SY) 179 : TY=(SX-SZ)*(SX-SZ) 180 : TZ=(SY-SZ)*(SY-SZ) 181 : R1=TX+TY+TZ+6*SXY*SXY 182 : R=SQRT(0.5*R1) 183 : IF (DABS(R).LE.1.D-10) THEN 184 : WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2' 185 : ENDIF 186 : DFSIG(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2 187 : DFSIG(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2 188 : DFSIG(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2 189 : DFSIG(4)=3.D0*SXY/R/A2 190 : P2(1,1)=2.D0/A2/A2 191 : P2(1,2)=-1.D0/A2/A2 192 : P2(1,3)=-1.D0/A2/A2 193 : P2(1,4)=0.D0/A2/A2 194 : P2(2,1)=-1.D0/A2/A2 195 : P2(2,2)=2.D0/A2/A2 196 : P2(2,3)=-1.D0/A2/A2 197 : P2(2,4)=0.D0/A2/A2 198 : P2(3,1)=-1.D0/A2/A2 199 : P2(3,2)=-1.D0/A2/A2 200 : P2(3,3)=2.D0/A2/A2 201 : P2(3,4)=0.D0/A2/A2 202 : P2(4,1)=0.D0/A2/A2 203 : P2(4,2)=0.D0/A2/A2 204 : P2(4,3)=0.D0/A2/A2 205 : P2(4,4)=6.D0/A2/A2 206 : ENDIF 207 : C 208 : DO 25 I=1,NSTRS 209 : AC(I)=0.D0 210 : DO 25 J=1,NSTRS 211 : AC(I)=AC(I)+D(I,J)*DFSIG(J) 212 : 25 CONTINUE 213 : C 214 : C ************ Preparation du jacobien *********************** 215 : C 216 : FDF=0.D0 217 : DO 30 J=1,NSTRS 218 : FDF=FDF+DFSIG(J)*AC(J) 219 : 30 CONTINUE 220 : C 221 : DJAC0=-(FDF+PAEC) 222 : C 223 : C ************ Debut des iterations internes ******************* 224 : C 225 : 40 CONTINUE 226 : C 227 : C *************** Determination de DK et DLAM ****************** 228 : C 229 : C 230 : DLAM1=DLAM0-FCRIC0/DJAC0 231 : C 232 : IF (DLAM1.LT.0.D0) IEC2=1 233 : IF (IEC2.EQ.1) THEN 234 : C WRITE(*,*)'Dans BEHAV2, DLAMDA est negatif:',DLAM1 235 : C WRITE(*,*)'A l iteration :',ITER 236 : GOTO 12 237 : ENDIF 238 : C 239 : DK=DEPS2+DLAM1 240 : C 241 : IF(IVIS.LE.2) THEN 241 : CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF) 242 : ELSE 242 : CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1) 243 : ENDIF 244 : C 245 : C ************** Determination de DPHI2 ********************* 246 : C 247 : IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN 248 : DPHI2=SEQ2 249 : ENDIF 250 : IF (IMOD.EQ.3) THEN 251 : SOM=SIGF(1)+SIGF(2)+SIGF(3) 252 : DPHI2=SEQ2-A1*SOM/A2 253 : ENDIF 254 : IF (IMOD.EQ.4) THEN 255 : SOM=SIGE(1)+SIGE(2)+SIGE(3) 256 : DPHI2=A1*SOM/A2-A1*A1*DLAM1*3.D0*ADD*(1.D0+XNU)/A2/A2 257 : DPHI2=SEQ2-DPHI2 258 : ENDIF 259 : C IF (DPHI2.LT.0.D0) THEN 260 : C WRITE(*,*)'ATTENTION DPHI2 ds BEHAV2 NEGATIF' 261 : C WRITE(*,*)'DPHI2=',DPHI2 262 : C ENDIF 263 : C 264 : C---------- Cas du critere reduit a un point ---------------- 265 : C 266 : IF (DABS(DPHI2).LE.10E-10) THEN 267 : IAPEX=1 268 : C WRITE(*,*)'IAPEX ds BEHAV2 =',IAPEX 269 : DO 50 I=1,NSTRS 270 : DO 50 J=1,NSTRS 271 : AI(I,J)=0.D0 272 : 50 CONTINUE 273 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 274 : DO 53 I=1,3 275 : DO 53 J=1,3 276 : AI(I,J)=1.D0/3.D0 277 : 53 CONTINUE 278 : ENDIF 279 : GOTO 75 280 : ENDIF 281 : C 282 : C ************** Mise a jour des contraintes *************** 283 : C 284 : C ---------------- calcul de la matrice A ------------------ 285 : C 286 : DO 60 I=1,NSTRS 287 : DO 60 J=1,NSTRS 288 : A(I,J)=0.D0 289 : 60 CONTINUE 290 : C 291 : DG=YOUN/2.D0/(1.D0+XNU) 292 : C 293 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 294 : DW1=AD*(2.D0-XNU) 295 : DW2=AD*(2.D0*XNU-1.D0) 296 : A(1,1)=1.D0+(DLAM1*DW1)/2.D0/DPHI2/A2/A2 297 : A(2,2)=A(1,1) 298 : A(1,2)=(DLAM1*DW2)/2.D0/DPHI2/A2/A2 299 : A(2,1)=A(1,2) 300 : A(3,3)=1.D0+(DLAM1*3.D0*DG)/DPHI2/A2/A2 301 : ENDIF 302 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 303 : A(1,1)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2 304 : A(2,2)=A(1,1) 305 : A(3,3)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2 306 : A(1,2)=-(DLAM1*DG)/DPHI2/A2/A2 307 : A(2,1)=A(1,2) 308 : A(1,3)=-(DLAM1*DG)/DPHI2/A2/A2 309 : A(3,1)=A(1,3) 310 : A(3,2)=A(1,3) 311 : A(2,3)=A(1,3) 312 : A(4,4)=1.D0+(DLAM1*DG*3.D0)/DPHI2/A2/A2 313 : ENDIF 314 : C 315 : C -------------- invertion de la matrice A ----------------- 316 : C 317 : DO 70 I=1,NSTRS 318 : DO 70 J=1,NSTRS 319 : AI(I,J)=A(I,J) 320 : 70 CONTINUE 320 : CALL INVMA2(AI,NSTRS,ISING) 321 : IF (ISING.EQ.1) THEN 322 : WRITE(*,*)'MATRICE AI singuliere ds BEHAV2' 323 : WRITE(*,*)'NSTRS=',NSTRS 324 : ENDIF 325 : C 326 : C -------------- mise a jour des contraintes -------------- 327 : C 328 : 75 CONTINUE 329 : IF(IMOD.EQ.3) THEN 330 : VEC1(1)=AD*(1.D0+XNU) 331 : VEC1(2)=AD*(1.D0+XNU) 332 : VEC1(3)=0.D0 333 : VEC1(4)=0.D0 334 : ENDIF 335 : IF(IMOD.EQ.4) THEN 336 : VEC1(1)=ADD*(1.D0+XNU) 337 : VEC1(2)=ADD*(1.D0+XNU) 338 : VEC1(3)=ADD*(1.D0+XNU)*2.D0*XNU 339 : VEC1(4)=0.D0 340 : ENDIF 341 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 342 : VEC1(1)=0.D0 343 : VEC1(2)=0.D0 344 : VEC1(3)=0.D0 345 : VEC1(4)=0.D0 346 : ENDIF 347 : C 348 : DO 80 I=1,NSTRS 349 : DEPSI(I)=SIGE(I)-A1*DLAM1*VEC1(I)/A2 350 : 80 CONTINUE 351 : C 352 : DO 90 I=1,NSTRS 353 : SIGF(I)=0.0D+00 354 : DO 90 J=1,NSTRS 355 : SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 356 : 90 CONTINUE 357 : C 358 : C ************ Verification du critere ******************** 359 : C 360 : IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN 360 : CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF) 361 : ENDIF 362 : IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN 362 : CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF) 363 : ENDIF 364 : FCRIC1=SEQ-SEQ2 365 : C 366 : IF (IBROY.EQ.0.AND.DABS(FCRIC1).GE.CRIMAX) THEN 367 : C WRITE(*,*)'****************************************' 368 : C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN' 369 : C WRITE(*,*)'on passe donc a la secante' 370 : C WRITE(*,*)'Dans l element',IBB 371 : C WRITE(*,*)'et au point d intégration',IGAU 372 : C WRITE(*,*)'CRIMAX=',CRIMAX 373 : C WRITE(*,*)'****************************************' 374 : ITER=ITR 375 : ENDIF 376 : C 377 : C ******* Compteur sur la methode de resolution **** 378 : C 379 : IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN 380 : IBROY=1 381 : ITANG=1 382 : IAPEX=0 383 : IEC2=0 384 : ITER=1 385 : GOTO 8 386 : ENDIF 387 : C 388 : C ******* non convergence ************************** 389 : C 390 : IF (DABS(FCRIC1).GT.PRB.AND.ITER.LT.ITR) THEN 391 : IF (IBROY.EQ.0) THEN 392 : DJAC0=(FCRIC0-FCRIC1)/(DLAM0-DLAM1) 393 : DLAM0=DLAM1 394 : FCRIC0=FCRIC1 395 : ITER=ITER+1 396 : IF (ITER.GE.(ITR-1)) THEN 397 : C WRITE(*,*)'***********************' 398 : C WRITE(*,*)'BROYDEN n a pas aboutit' 399 : C WRITE(*,*)'Dans l element',IBB 400 : C WRITE(*,*)'et au point d intégration',IGAU 401 : C WRITE(*,*)'ITER=',ITER 402 : C WRITE(*,*)'FCRI=',FCRIC1 403 : C WRITE(*,*)'***********************' 404 : ENDIF 405 : GOTO 40 406 : ENDIF 407 : IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN 408 : DLAM0=DLAM1 409 : FCRIC0=FCRIC1 410 : ITER=ITER+1 411 : IF (ITER.GE.(ITR-5)) THEN 412 : C WRITE(*,*)'ITER=',ITER 413 : C WRITE(*,*)'FCRI=',FCRIC1 414 : ENDIF 415 : GOTO 9 416 : ENDIF 417 : ENDIF 418 : IF (ITER.GE.ITR.AND.DABS(FCRIC1).GT.PRB2) THEN 419 : WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV2' 420 : WRITE(*,*)'FCRI=',FCRIC1 421 : WRITE(*,*)'Dans l element',IBB 422 : WRITE(*,*)'et au point d intégration',IGAU 423 : WRITE(*,*)'IPLA=',IPLA 424 : WRITE(*,*)'IFIS=',IFIS 425 : C STOP 426 : ENDIF 427 : C 428 : C ************* Calcul de la DEP **************** 429 : C 430 : IF (IAPEX.EQ.0) THEN 430 : CALL DERI2(D2FSIG,DPHI2,P2,SIGF,NSTRS,BETJEF) 430 : CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU, 431 : & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF) 432 : ELSE IF (IAPEX.EQ.1) THEN 433 : C WRITE(*,*)'ELAS2',IAPEX 434 : DO 96 I=1,NSTRS 435 : DO 96 J=1,NSTRS 436 : DEP(I,J)=0.D0 437 : DO 96 K=1,NSTRS 438 : DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J) 439 : 96 CONTINUE 440 : ENDIF 441 : IAPEX=0 442 : CRIMAX=0.D0 443 : * DO 97 I=1,NSTRS 444 : * DO 97 J=1,NSTRS 445 : * DEP(I,J)=D(I,J) 446 : * 97 CONTINUE 447 : C 448 : C *********************************************** 449 : C 450 : DEPS2=DK 451 : SIG2=SEQ2 452 : IF(IVIS.LE.2) THEN 452 : CALL INDICA(0.D0,DK,IFIS,IPLA,2,BETJEF) 453 : ELSE 453 : CALL INDIC1(0.D0,DK,IFIS,IPLA,2,BETJEF,NECH0,NECH1) 454 : ENDIF 455 : C ********************** CALCUL VISCOPLASTIQUE *********************** 456 : IF (IVIS.EQ.1) THEN 456 : CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2, 457 : & SIGP,SIGRV,DSTRN,D,BETJEF,VISCO) 458 : ENDIF 459 : C ******************************************************************** 460 : 100 CONTINUE 461 : IEC2=0 462 : C 463 : RETURN 464 : END 465 : 466 : 467 : 468 : 469 : 470 :
© Cast3M 2003 - All rights reserved.
Disclaimer