1 : C BEHAV1 SOURCE PV 07/11/23 21:15:27 5978 2 : SUBROUTINE BEHAV1(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 POST PIC 9 : C Un seul critere de traction: Rankine 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),SIGRV(4) 17 : DIMENSION V1(4),AC(4),SIGE(4),SIGP(4),SIGR(4),DSTRN(4) 18 : DIMENSION D2FSIG(4,4),DEP(4,4),P1(4,4),D(4,4),DP(4,4),A(4,4) 19 : DIMENSION 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 : C 38 : * COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF, 39 : * & TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR, 40 : * & ISIM,IBB,IGAU,IZON 41 : C 42 : C 43 : IAPEX=0 44 : PRB=1.D-10 45 : PRB2=1.D-6 46 : ITER=1 47 : ITANG=0 48 : IBROY=0 49 : Ft=PALF*COLI 50 : CRIMAX=0.D0 51 : SEQ = 0.D0 52 : SEQ1 = 0.D0 52 : CALL ZERO(SIGE,4,1) 52 : CALL ZERO(V1,4,1) 52 : CALL ZERO(P1,4,4) 52 : CALL ZERO(D,4,4) 53 : C 54 : DO 10 I=1,NSTRS 55 : SIGE(I)=SIGF(I) 56 : 10 CONTINUE 57 : C 58 : C DO 11 I=1,NSTRS 59 : C DO 11 J=1,NSTRS 60 : C D(I,J)=0.D0 61 : C 11 CONTINUE 62 : C 63 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 64 : AD=YOUN/(1.D0-XNU*XNU) 65 : D(1,1)=AD 66 : D(2,2)=D(1,1) 67 : D(3,3)=AD*(1.D0-XNU)/2.D0 68 : D(1,2)=AD*XNU 69 : D(2,1)=D(1,2) 70 : ENDIF 71 : C 72 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 73 : ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU)) 74 : D(1,1)=ADD*(1.D0-XNU) 75 : D(2,2)=D(1,1) 76 : D(3,3)=D(1,1) 77 : D(1,2)=ADD*XNU 78 : D(2,1)=D(1,2) 79 : D(1,3)=D(1,2) 80 : D(2,3)=D(1,2) 81 : D(3,1)=D(1,2) 82 : D(3,2)=D(1,2) 83 : D(4,4)=0.5*ADD*(1.D0-2.D0*XNU) 84 : ENDIF 85 : C 86 : C ************ Le point est fissure pour 1ere fois *********** 87 : C 88 : IF (IFIS.EQ.0) THEN 88 : CALL PRINC(SIGF,V1,NSTRS) 89 : DEPS1=0.D0 90 : SIG1=Ft 91 : IF (V1(1).GT.Ft) THEN 92 : IFIS=1 93 : ENDIF 94 : ENDIF 95 : IF (IPLA.EQ.0) THEN 96 : DEPS2=0.D0 97 : SIG2=COLI*AA 98 : ENDIF 99 : C 100 : 8 CONTINUE 101 : C 102 : DO 4 I=1,NSTRS 103 : SIGF(I)=SIGE(I) 104 : 4 CONTINUE 105 : C 105 : CALL PRINC(SIGF,V1,NSTRS) 106 : TETA=V1(4) 107 : DK=DEPS1 108 : DLAM0=0.D0 109 : FCRI0=V1(1)-SIG1 110 : CRIMAX=DABS(100.D0*FCRI0) 111 : C 112 : IF (FCRI0.LT.0.D0) THEN 113 : DO 45 I=1,NSTRS 114 : DO 45 J=1,NSTRS 115 : DEP(I,J)=D(I,J) 116 : 45 CONTINUE 117 : GOTO 100 118 : ENDIF 119 : C 120 : C ************ Traitement du point deja fissure ******************** 121 : C 122 : 9 CONTINUE 123 : PI=4.D0*DATAN(1.D0) 124 : PHIC=V1(4)*(PI/180.D0) 125 : COSA=DCOS(PHIC) 126 : SINA=DSIN(PHIC) 127 : C------------------------------------------------------------------- 128 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 129 : DFSIG(1)=COSA*COSA 130 : DFSIG(2)=SINA*SINA 131 : DFSIG(3)=2.D0*SINA*COSA 132 : DO 12 I=1,NSTRS 133 : DO 12 J=1,NSTRS 134 : P1(I,J)=0.D0 135 : 12 CONTINUE 136 : P1(1,1)=1.D0/2.D0 137 : P1(1,2)=-1.D0/2.D0 138 : P1(2,1)=-1.D0/2.D0 139 : P1(2,2)=1.D0/2.D0 140 : P1(3,3)=2.D0 141 : ENDIF 142 : C------------------------------------------------------------------- 143 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 144 : DFSIG(1)=COSA*COSA 145 : DFSIG(2)=SINA*SINA 146 : DFSIG(3)=0.D0 147 : DFSIG(4)=2.D0*SINA*COSA 148 : DO 13 I=1,NSTRS 149 : DO 13 J=1,NSTRS 150 : P1(I,J)=0.D0 151 : 13 CONTINUE 152 : P1(1,1)=1.D0/2.D0 153 : P1(1,2)=-1.D0/2.D0 154 : P1(2,1)=-1.D0/2.D0 155 : P1(2,2)=1.D0/2.D0 156 : P1(4,4)=2.D0 157 : ENDIF 158 : C------------------------------------------------------------------- 159 : DO 20 I=1,NSTRS 160 : AC(I)=0.D0 161 : DO 20 J=1,NSTRS 162 : AC(I)=AC(I)+D(I,J)*DFSIG(I) 163 : 20 CONTINUE 164 : C 165 : FDF=0.D0 166 : DO 30 J=1,NSTRS 167 : FDF=FDF+DFSIG(J)*AC(J) 168 : 30 CONTINUE 169 : C--------------- Determination du parametre d'ecrouissage ---------- 170 : C 171 : IF(IVIS.LE.2) THEN 171 : CALL UNICOU(DK,PAEC,1,SEQ,BETJEF) 172 : ELSE 172 : CALL UNICO1(DK,PAEC,1,SEQ,BETJEF,NECH0,NECH1) 173 : ENDIF 174 : C 175 : C------------------------------------------------------------------- 176 : DJAC0=-(PAEC+FDF) 177 : C 178 : C ************ Debut des iterations internes *************** 179 : C 180 : 40 CONTINUE 181 : C 182 : C *************** Determination de DK ********************* 183 : C 184 : DLAM1=-FCRI0/DJAC0+DLAM0 185 : DK=DEPS1+DLAM1 186 : C IF (DLAM1.LT.0.D0) THEN 187 : C WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1 188 : C WRITE(*,*)'A l iteration :',ITER 189 : C ENDIF 190 : C 191 : C--------------- Estimation contrainte quivalente ---------- 192 : C 193 : IF(IVIS.LE.2) THEN 193 : CALL UNICOU(DK,PAEC,1,SEQ1,BETJEF) 194 : ELSE 194 : CALL UNICO1(DK,PAEC,1,SEQ1,BETJEF,NECH0,NECH1) 195 : ENDIF 196 : C 197 : C ************** Determination de DPHI ********************* 198 : C 199 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 200 : AD1=YOUN/(1.D0-XNU) 201 : TO=SIGF(3) 202 : ENDIF 203 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 204 : AD1=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU)) 205 : TO=SIGF(4) 206 : ENDIF 207 : C 208 : DPHI=SEQ1-0.5*(SIGE(1)+SIGE(2))+0.5*AD1*DLAM1 209 : C IF (ITER.EQ.1) THEN 210 : C DPHI=0.25*(SIGF(1)-SIGF(2))*(SIGF(1)-SIGF(2)) 211 : C DPHI=DPHI+TO*TO 212 : C DPHI=DSQRT(DPHI) 213 : C ENDIF 214 : IF (DPHI.LT.0.D0) THEN 215 : C WRITE(*,*)'ATTENTION DPHI NEGATIF' 216 : C WRITE(*,*)'DPHI=',DPHI 217 : ENDIF 218 : C 219 : C--------------- Cas de l'apex ----------------------------- 220 : C 221 : IF (DABS(DPHI).LE.10E-10) THEN 222 : IAPEX=1 223 : C WRITE(*,*)'IAPEX ds BEHAV1 =',IAPEX 224 : C WRITE(*,*)'Dans l element',IBB 225 : C WRITE(*,*)'et au point d intégration',IGAU 226 : DO 50 I=1,NSTRS 227 : DO 50 J=1,NSTRS 228 : AI(I,J)=0.D0 229 : 50 CONTINUE 230 : AI(1,1)=0.5 231 : AI(1,2)=0.5 232 : AI(2,1)=AI(1,2) 233 : AI(2,2)=AI(1,1) 234 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0 235 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0 236 : GOTO 75 237 : ENDIF 238 : C 239 : C ************** Mise a jour des contraintes *************** 240 : C 241 : C ---------------- calcul de la matrice A ------------------ 242 : C 243 : DO 60 I=1,NSTRS 244 : DO 60 J=1,NSTRS 245 : A(I,J)=0.D0 246 : 60 CONTINUE 247 : C 248 : DG=YOUN/2.D0/(1.D0+XNU) 249 : C 250 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 251 : A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI 252 : A(2,2)=A(1,1) 253 : A(1,2)=-(A(1,1)-1.D0) 254 : A(2,1)=A(1,2) 255 : A(3,3)=1.D0+(DLAM1*DG)/DPHI 256 : ELSE 257 : A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI 258 : A(2,2)=A(1,1) 259 : A(1,2)=-(A(1,1)-1.D0) 260 : A(2,1)=A(1,2) 261 : A(3,3)=1.D0 262 : A(4,4)=1.D0+(DLAM1*DG)/DPHI 263 : ENDIF 264 : C 265 : C -------------- invertion de la matrice A ----------------- 266 : C 267 : DO 70 I=1,NSTRS 268 : DO 70 J=1,NSTRS 269 : AI(I,J)=A(I,J) 270 : 70 CONTINUE 270 : CALL INVMA2(AI,NSTRS,ISING) 271 : IF (ISING.EQ.1) THEN 272 : WRITE(*,*)'MATRICE AI singuliere ds BEHAV1' 273 : ENDIF 274 : C 275 : C -------------- mise a jour des contraintes ------------ 276 : C 277 : 75 CONTINUE 278 : IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN 279 : VEC1(1)=AD1 280 : VEC1(2)=AD1 281 : VEC1(3)=0.D0 282 : VEC1(4)=0.D0 283 : ENDIF 284 : IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN 285 : VEC1(1)=AD1 286 : VEC1(2)=AD1 287 : VEC1(3)=AD1*2.D0*XNU 288 : VEC1(4)=0.D0 289 : ENDIF 290 : C 291 : DO 80 I=1,NSTRS 292 : DEPSI(I)=SIGE(I)-0.5*DLAM1*VEC1(I) 293 : 80 CONTINUE 294 : C 295 : DO 90 I=1,NSTRS 296 : SIGF(I)=0.0D+00 297 : DO 90 J=1,NSTRS 298 : SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J) 299 : 90 CONTINUE 300 : C 301 : C ******** Verification du critere **************** 302 : C 302 : CALL PRINC(SIGF,V1,NSTRS) 303 : FCRI1=V1(1)-SEQ1 304 : C 305 : IF (IBROY.EQ.0.AND.DABS(FCRI1).GE.CRIMAX) THEN 306 : C WRITE(*,*)'****************************************' 307 : C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN' 308 : C WRITE(*,*)'on passe donc a la secante' 309 : C WRITE(*,*)'Dans l element',IBB 310 : C WRITE(*,*)'et au point d intégration',IGAU 311 : C WRITE(*,*)'CRIMAX=',CRIMAX 312 : C WRITE(*,*)'****************************************' 313 : ITER=ITR 314 : ENDIF 315 : C 316 : C ******* Compteur sur la methode de resolution **** 317 : C 318 : IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN 319 : IBROY=1 320 : ITANG=1 321 : ITER=1 322 : IAPEX=0 323 : GOTO 8 324 : ENDIF 325 : C 326 : C ******* non convergence ************************** 327 : C 328 : IF (DABS(FCRI1).GT.PRB.AND.ITER.LT.ITR) THEN 329 : IF (IBROY.EQ.0) THEN 330 : DJAC1=(FCRI0-FCRI1)/(DLAM0-DLAM1) 331 : DLAM0=DLAM1 332 : DJAC0=DJAC1 333 : FCRI0=FCRI1 334 : ITER=ITER+1 335 : IF (ITER.GE.(ITR-1)) THEN 336 : C WRITE(*,*)'***********************' 337 : C WRITE(*,*)'BROYDEN n a pas aboutit' 338 : C WRITE(*,*)'Dans l element',IBB 339 : C WRITE(*,*)'et au point d intégration',IGAU 340 : C WRITE(*,*)'ITER=',ITER 341 : C WRITE(*,*)'FCRI=',FCRI1 342 : C WRITE(*,*)'***********************' 343 : ENDIF 344 : GOTO 40 345 : ENDIF 346 : IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN 347 : DLAM0=DLAM1 348 : FCRI0=FCRI1 348 : CALL PRINC(SIGF,V1,NSTRS) 349 : ITER=ITER+1 350 : IF (ITER.GE.(ITR-5)) THEN 351 : WRITE(*,*)'ITER=',ITER 352 : WRITE(*,*)'FCRI=',FCRI1 353 : ENDIF 354 : GOTO 9 355 : ENDIF 356 : ENDIF 357 : IF (ITER.GE.ITR.AND.DABS(FCRI1).GT.PRB2) THEN 358 : WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV1' 359 : WRITE(*,*)'FCRI=',FCRI1 360 : WRITE(*,*)'Dans l element',IBB 361 : WRITE(*,*)'et au point d intégration',IGAU 362 : WRITE(*,*)'IPLA=',IPLA 363 : WRITE(*,*)'IFIS=',IFIS 364 : C STOP 365 : ENDIF 366 : C 367 : C ************* Calcul de la DEP **************** 368 : C 369 : IF (IAPEX.EQ.0) THEN 369 : CALL DERI2(D2FSIG,DPHI,P1,SIGF,NSTRS,BETJEF) 369 : CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU, 370 : & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF) 371 : ELSE 372 : C WRITE(*,*)'ELAS',IAPEX 373 : DO 95 I=1,NSTRS 374 : DO 95 J=1,NSTRS 375 : AH(I,J)=0.D0 376 : DO 95 K=1,NSTRS 377 : AH(I,J)=AH(I,J)+AI(I,K)*D(K,J) 378 : 95 CONTINUE 379 : DFSIG(1)=0.5*DSQRT(2.D0) 380 : DFSIG(2)=0.5*DSQRT(2.D0) 381 : DFSIG(3)=0.D0 382 : DFSIG(4)=0.D0 382 : CALL CREDEP(AH,DFSIG,PAEC,NSTRS,DEP,BETJEF) 383 : ENDIF 384 : IAPEX=0 385 : CRIMAX=0.D0 386 : C 387 : C *********************************************** 388 : C 389 : IF (DK.LT.0.D0) THEN 390 : WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1 391 : WRITE(*,*)'A l iteration :',ITER 392 : C STOP 393 : ENDIF 394 : DEPS1=DK 395 : SIG1=V1(1) 396 : IF (V1(1).LT.0.D0) V1(1)=0.D0 396 : CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF) 397 : IF(IVIS.LE.2) THEN 397 : CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF) 398 : ELSE 398 : CALL INDIC1(DK,0.0D0,IFIS,IPLA,1,BETJEF,NECH0,NECH1) 399 : ENDIF 400 : C ********************** CALCUL VISCOPLASTIQUE *********************** 401 : IF (IVIS.EQ.1) THEN 401 : CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2, 402 : & SIGP,SIGRV,DSTRN,D,BETJEF,VISCO) 403 : ENDIF 404 : C ******************************************************************** 405 : 100 CONTINUE 406 : C 407 : RETURN 408 : END 409 : 410 : 411 : 412 : 413 : 414 :
© Cast3M 2003 - All rights reserved.
Disclaimer