1 : C CCOIN0 SOURCE FANDEUR 11/04/12 21:15:09 6938 2 : 3 : SUBROUTINE CCOIN0(wrk52,wrk53,wrk54,wrk2,wrk3, 4 : & IB,IGAU,NBPGAU,LTRAC,IFOURB,iecou) 5 : 6 : C----------------------------------------------------------------------- 7 : C ECOULEMENT PLASTIQUE POUR UN POINT 8 : C ALGORITHME ORTIZ ET SIMO 9 : C 10 : C EN ENTREE : 11 : C 12 : C SIG0 CONTRAINTES AU DEBUT DU PAS 13 : C DEPST INCREMENT DE DEFORMATIONS TOTALES 14 : C ( THERMIQUE ENLEVEE ) 15 : C VAR0 VARIABLES INTERNES DEDUT DU PAS 16 : C VAREX0 VARIABLES EXTERNES DEBUT DU PAS 17 : C VAREXF VARIABLES EXTERNES FIN DU PAS 18 : C XMAT COEFFICIENTS DU MATERIAU 19 : C PRECIS PRECISION POUR ITERATIONS INTERNES 20 : C WORK DES CARACTERISTIQUES 21 : C TRAC COURBE DE TRACTION 22 : C MFR1 INDICE DE FORMULATION 23 : C NSTRSS NOMBRE DE CONTRAINTES CA2000 24 : C INPLAS NUMERO DU MODELE DE PLASTICITE 25 : C DDAUX = MATRICE DE HOOKE ELASTIQUE 26 : C CMATE = NOM DU MATERIAU 27 : C VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU 28 : C VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES 29 : C N2EL = NBRE D ELEMENTS DANS SEGMENT DE HOOKE 30 : C N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE 31 : C IFOU = OPTION DE CALCUL 32 : C IB = NUMERO DE L ELEMENT COURANT 33 : C IGAU = NUMERO DU POINT COURANT 34 : C EPAIST= EPAISSEUR 35 : C NBPGAU= NBRE DE POINTS DE GAUSS 36 : C MELE = NUMERO DE L ELEMENT FINI 37 : C NPINT = NBRE DE POINTS D INTEGRATION 38 : C NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES 39 : C NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES 40 : C SECT = SECTION 41 : C LHOOK = TAILLE DE LA MATRICE DE HOOKE 42 : C TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES 43 : C UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE 44 : C 45 : C EN SORTIE : 46 : C 47 : C SIGF CONTRAINTES A LA FIN DU PAS 48 : C VARF VARIABLES INTERNES FIN DU PAS 49 : C DEFP INCR. DE DEFORMATIONS PLASTIQUES 50 : C KERRE CODE D'ERREUR 51 : C = 0 SI TOUT OK 52 : C = 99 SI FORMULATION NON DISPONIBLE 53 : C EN PLASTICITE 54 : C = 1 SI DEPS EST NEGATIF 55 : C = 2 SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE 56 : C 57 : C IFOUR : OPTION DE CALCUL 58 : C 59 : C IFOUR = -3 DEFORMATION PLANE GENERALISEE 60 : C IFOUR = -2 CONTRAINTES PLANES 61 : C IFOUR = -1 DEFORMATIONS PLANES 62 : C IFOUR = 0 AXISYMETRIQUE 63 : C IFOUR = 1 SERIE DE FOURIER 64 : C IFOUR = 2 TRIDIMENSIONNEL 65 : C 66 : C MFR1 : NUMERO DE LA FORMULATION ELEMENT FINI 67 : C 68 : C MFR1 = 1 MASSIF 69 : C MFR1 = 3 COQUE 70 : C MFR1 = 5 COQUE EPAISSE 71 : C MFR1 = 7 POUTRE 72 : C MFR1 = 9 COQUE AVEC CISAILLEMENT TRANSVERSE 73 : C 74 : C INPLAS : NUMERO DU MODELE DE PLASTICITE 75 : C 76 : C INPLAS = 1 PARFAIT 77 : C INPLAS = 4 CINEMATIQUE 78 : C INPLAS = 5 ISOTROPE 79 : C INPLAS = 7 CHABOCHE1 80 : C INPLAS = 12 CHABOCHE2 81 : C 82 : C----------------------------------------------------------------------- 83 : C CONVENTION DE REMPLISSAGE DES MEMOIRES 84 : C----------------------------------------------------------------------- 85 : C 86 : C XMAT(1) MODULE D'YOUNG 87 : C XMAT(2) COEFFICIENT DE POISSON 88 : C 89 : C TRAC(1 A 2*LTRAC) COURBE DE TRACTION 90 : C WORK( " +1) ALFAH POUR COQUES PLASTICITE GLOBALE 91 : C WORK( " +..) DONNEES POUR CRITERE POUTRES GLOBALES 92 : C 93 : C MODELE ISOTROPE 94 : C VAR0(1) EPS* 95 : C 96 : C MODELE CINEMATIQUE LINEAIRE 97 : C VAR0(1) EPS* 98 : C VAR0(2) A VAR0(1+NSTRSS) DEPLACEMENT DE LA SPHERE 99 : C 100 : C MODELE CHABOCHE 101 : C XMAT(5) .... COEFFICIENTS A,C,... 102 : C VAR0(1) EPS* 103 : C VAR0(2) A VAR0(1+NSTRSS) DEPLACEMENT DE LA SPHERE 1 104 : C VAR0(2+NSTRSS) A VAR0(1+2*NSTRSS) " " " " 2 105 : C 106 : C----------------------------------------------------------------------- 107 : 108 : IMPLICIT INTEGER(I-N) 109 : IMPLICIT REAL*8(A-H,O-Z) 110 : 111 : -INC CCOPTIO 112 : -INC DECHE 113 : 114 : SEGMENT IECOU 115 : INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7, 116 : 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16, 117 : 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24, 118 : 3 icow25,icow26,icow27,icow28,icow29,icow30,icow31, 119 : 4 icow32,icow33,NSTRSS,MFR1, NBGMAT,NELMAT,icow38, 120 : 5 icow39,icow40,icow41,icow42,icow43,icow44 121 : INTEGER icow45,icow46,icow47,icow48,icow49,icow50, 122 : . icow51,icow52,icow53,icow54,icow55,icow56 123 : . icow57,icow58 124 : ENDSEGMENT 125 : * 126 : SEGMENT WRK2 127 : REAL*8 TRAC(LTRAC) 128 : ENDSEGMENT 129 : * 130 : SEGMENT WRK3 131 : REAL*8 WORK(LW),WORK2(LW2) 132 : ENDSEGMENT 133 : 134 : DIMENSION SIG(130),EPS(130) 135 : DIMENSION S(8),SX(8),DS(8),DSIG(8),SPHER(8),SPHER1(8),SPHER2(8) 136 : DIMENSION DSPHER1(8),DSPHER2(8),DEPSE(8),DEPSP(8),DDEPSE(8) 137 : DIMENSION F(8),W1(8),W2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8) 138 : DIMENSION CRIGI(12) 139 : 140 : C---------COQUES AVEC CT------------------------------------------------ 141 : C ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES 142 : IF(MFR1.EQ.9) THEN 143 : NCONT=6 144 : ELSE 145 : NCONT=NSTRSS 146 : ENDIF 147 : itracb=0 148 : 2222 continue 149 : 150 : C----------------------------------------------------------------------- 151 : DO I=1,NSTRSS 152 : S(I)=0.D0 153 : SPHER(I)=0.D0 154 : SPHER1(I)=0.D0 155 : SPHER2(I)=0.D0 156 : ENDDO 157 : KERRE=0 158 : YUNG=XMAT(1) 159 : XNU=XMAT(2) 160 : 161 : C---------CARACTERISTIQUES GEOMETRIQUES--------------------------------- 162 : C 163 : C COQUES 164 : C 165 : ALFAH=1.D0 166 : IF(MFR1.EQ.3.OR.MFR1.EQ.5.OR.MFR1.EQ.9) THEN 167 : EP1=WORK(1) 168 : IF(MFR1.NE.5) ALFAH=WORK(2)**2 169 : ENDIF 170 : 171 : C---------COEFFICIENTS CHABOCHE----------------------------------------- 172 : 173 : IF(INPLAS.EQ.7) THEN 174 : A1=XMAT(5) 175 : C1=XMAT(6) 176 : R0=XMAT(7) 177 : PSI=XMAT(8) 178 : OME=XMAT(9) 179 : RM=XMAT(10) 180 : B=XMAT(11) 181 : A2=0.D0 182 : C2=0.D0 183 : ELSE IF(INPLAS.EQ.12) THEN 184 : A1=XMAT(5) 185 : C1=XMAT(6) 186 : A2=XMAT(7) 187 : C2=XMAT(8) 188 : R0=XMAT(9) 189 : PSI=XMAT(10) 190 : OME=XMAT(11) 191 : RM=XMAT(12) 192 : B=XMAT(13) 193 : ELSE 194 : DO I=1,LTRAC 195 : SIG(I)=TRAC(2*I-1) 196 : EPS(I)=TRAC(2*I) 197 : ENDDO 198 : ENDIF 199 : 200 : EPST=VAR0(1) 201 : C---------ECROUISSAGE CINEMATIQUE--------------------------------------- 202 : 203 : IF(INPLAS.EQ.4.OR.INPLAS.EQ.7.OR.INPLAS.EQ.12) THEN 204 : DO I=1,NSTRSS 205 : SPHER1(I)=VAR0(I+1) 206 : ENDDO 207 : IF(INPLAS.EQ.12) THEN 208 : DO I=1,NSTRSS 209 : SPHER2(I)=VAR0(NSTRSS+1+I) 210 : ENDDO 211 : ENDIF 212 : DO I=1,NSTRSS 213 : SPHER(I)=SPHER1(I)+SPHER2(I) 214 : ENDDO 215 : C----------------------------------------------------------------------- 216 : C PREDICTEUR ELASTIQUE 217 : C S : CONTRAINTE 218 : C SPHER : VARIABLE D'ECROUISSAGE CINEMATIQUE 219 : C SX = S - X 220 : C----------------------------------------------------------------------- 221 : 222 : * en elastique non lineaire on annule les contraintes initiales 223 : * mais on cumule les epsilons elastiques 224 : ELSE IF(INPLAS.EQ.87) THEN 225 : EPST=0.D0 226 : DO I=1,NSTRSS 227 : SIG0(I)=0.D0 228 : DEPST(I)=DEPST(I) + VAR0(I+1) 229 : ENDDO 230 : ENDIF 231 : 231 : CALL CALSIG(DEPST,DDAUX,NSTRSS,CMATE,VALMAT,VALCAR, 232 : & N2EL,N2PTEL,MFR1,IFOURB,IB,IGAU,EPAIST, 233 : & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR, 234 : & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD) 235 : 236 : IF(IRTD.NE.1) THEN 237 : KERRE=69 238 : GOTO 1000 239 : ENDIF 240 : 241 : IF (mfr.eq.1.and. ifourb.eq.-2) then 242 : * en cas de contraintes planes on repasse en 3D 243 : do iu=1,6 244 : do ju=1,6 245 : DDA(iu,ju)=0.d0 246 : enddo 247 : enddo 248 : cte_cp = xnu / (1.d0 - xnu) 249 : aux= YUNG / ((1.d0+xnu)*(1.d0-2.d0*xnu)) 250 : aux1= aux * xnu 251 : aux2= aux * (1.d0-xnu) 252 : gege = Yung / (2.d0 *(1.d0 +xnu)) 253 : DDA(1,1)=AUX2 254 : DDA(2,1)=AUX1 255 : DDA(1,2)=AUX1 256 : DDA(2,2)=AUX2 257 : DDA(3,3)=aux2 258 : DDA(1,3)=aux1 259 : DDA(2,3)=aux1 260 : DDA(3,1)=aux1 261 : DDA(3,2)=aux1 262 : 263 : DDA(4,4)=gege 264 : DDA(5,5)=gege 265 : DDA(6,6)=GEGE 266 : ELSE IF ((MFR.EQ.3.AND.NPINT.EQ.0).OR.MFR.EQ.9) THEN 267 : AUX=YUNG/(1.D0-XNU*XNU) 268 : AUX1=AUX*XNU 269 : DO I=1,NCONT 270 : DO J=1,NCONT 271 : DDAUX(I,J)=0.D0 272 : ENDDO 273 : ENDDO 274 : C 275 : C CAS TRIDIMENSIONNEL ET FOURIER 276 : C 277 : IF(IFOURB.EQ.2.OR.IFOURB.EQ.1) THEN 278 : 279 : GEGE=0.5D0*YUNG/(1.D0+XNU) 280 : DDAUX(1,1)=AUX 281 : DDAUX(2,1)=AUX1 282 : DDAUX(1,2)=AUX1 283 : DDAUX(2,2)=AUX 284 : DDAUX(3,3)=GEGE 285 : DDAUX(4,4)=AUX 286 : DDAUX(5,4)=AUX1 287 : DDAUX(4,5)=AUX1 288 : DDAUX(5,5)=AUX 289 : DDAUX(6,6)=GEGE 290 : C 291 : C CAS AXISYMETRIQUE ET DEFORMATIONS PLANES 292 : C 293 : ELSE IF(IFOURB.EQ.0.OR.IFOURB.EQ.-1.OR.IFOURB.EQ.-3) THEN 294 : 295 : DDAUX(1,1)=AUX 296 : DDAUX(2,1)=AUX1 297 : DDAUX(1,2)=AUX1 298 : DDAUX(2,2)=AUX 299 : DDAUX(3,3)=AUX 300 : DDAUX(4,3)=AUX1 301 : DDAUX(3,4)=AUX1 302 : DDAUX(4,4)=AUX 303 : C 304 : C CAS CONTRAINTES PLANES 305 : C 306 : ELSE IF(IFOURB.EQ.-2) THEN 307 : DDAUX(1,1)=YUNG 308 : DDAUX(3,3)=YUNG 309 : ENDIF 310 : ENDIF 311 : * 312 : DO I=1,NSTRSS 313 : S(I)=SIG0(I)+DSIGT(I) 314 : SIGB(I)=S(I) 315 : SX(I)=S(I)-SPHER(I) 316 : ENDDO 317 : * if( igau.eq.1.and.itracb.eq.1) then 318 : * write(6,*) ' entre dans cconi0' 319 : * write(6,*) ' sig0' , sig0(1),sig0(2),sig0(3),sig0(4) 320 : * write(6,*) ' dsigt' , dsigt(1),dsigt(2),dsigt(3),dsigt(4) 321 : * write(6,*) ' s' , s(1),s(2),s(3),s(4) 322 : * write(6,*) ' sx' , sx(1),sx(2),sx(3),sx(4) 323 : * write(6,*) ' spher' , spher(1),spher(2),spher(3),spher(4) 324 : * endif 325 : 326 : C---------CAS DES POUTRES----------------------------------------------- 327 : 328 : IF(MFR1.EQ.7) THEN 329 : DIV(1)=1.D0/WORK(4) 330 : DIV(2)=1.D0 331 : DIV(3)=1.D0 332 : DIV(4)=WORK(10)/WORK(1) 333 : DIV(5)=WORK(11)/WORK(2) 334 : DIV(6)=WORK(12)/WORK(3) 335 : IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(WORK(1)*WORK(4)) 336 : IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(WORK(2)*WORK(4)) 337 : IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(WORK(3)*WORK(4)) 338 : DO I=1,NCONT 339 : S(I)=S(I)*DIV(I) 340 : SX(I)=SX(I)*DIV(I) 341 : ENDDO 342 : ENDIF 343 : 344 : C----------------------------------------------------------------------- 345 : C CALCUL DE LA LIMITE ELASTIQUE SI 346 : C----------------------------------------------------------------------- 347 : 348 : IF(INPLAS.EQ.1.OR.INPLAS.EQ.4) THEN 349 : SI=TRAC(1) 350 : ELSE IF(INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN 350 : CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI) 351 : IF(IBI.EQ.1) THEN 352 : KERRE=75 353 : GOTO 1000 354 : ENDIF 355 : C* Modele de CHABOCHE (prise en compte ecrouissage) 356 : ELSE IF (INPLAS.EQ.7 .OR. INPLAS.EQ.12) THEN 357 : RMmRR = (RM - R0) * EXP(-B*EPST) 358 : SI = RM - RMmRR 359 : ENDIF 360 : 361 : C----------------------------------------------------------------------- 362 : C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ 363 : C----------------------------------------------------------------------- 364 : * attention en contraintes planes on se declare en 3D 365 : * rien besoin de faire dans vonmis0 car ifourb n'est pas utilisé 366 : * et les contraintes sont dimensionnées à 6 367 : * if(igau.eq.1) then 368 : * write(6,*) ' en entrée sx ',sx(1),sx(2),sx(3) 369 : * endif 370 : SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) 371 : 372 : C----------------------------------------------------------------------- 373 : C LE CRITERE EST-IL VERIFIE 374 : C----------------------------------------------------------------------- 375 : * if(itracb.eq.1) then 376 : * write(6,*)' en entrée si seq', si, seq 377 : * endif 378 : PHI=SEQ-SI 379 : NITER=0 380 : PETI=1.1D0*0.5D0*PRECIS*SEQ 380 : CALL EPSPRE(SEQ,SI,PETI,ITRY) 381 : IF((ITRY.EQ.1).OR.(SEQ.LE.SI)) THEN 382 : * rien a faire on n'a pas plastifié 383 : IF(MFR1.EQ.7) THEN 384 : DO I=1,NCONT 385 : S(I)=S(I)/DIV(I) 386 : ENDDO 387 : ENDIF 388 : DO I=1,NCONT 389 : SIGF(I)=S(I) 390 : DEFP(I)=0.D0 391 : ENDDO 392 : IF(MFR1.EQ.9) THEN 393 : DEFP(7)=0.D0 394 : DEFP(8)=0.D0 395 : SIGF(7)=S(7) 396 : SIGF(8)=S(8) 397 : ENDIF 398 : 399 : VARF(1)=VAR0(1) 400 : IF(INPLAS.EQ.4.OR.INPLAS.EQ.7) THEN 401 : DO I=1,NSTRSS 402 : VARF(I+1)=VAR0(I+1) 403 : ENDDO 404 : ELSE IF(INPLAS.EQ.12) THEN 405 : DO I=1,2*NSTRSS 406 : VARF(I+1)=VAR0(I+1) 407 : ENDDO 408 : ELSE IF(INPLAS.EQ.87) THEN 409 : DO I=1,NSTRSS 410 : VARF(I+1)=DEPST(I) 411 : ENDDO 412 : ENDIF 413 : RETURN 414 : ENDIF 415 : 416 : C----------------------------------------------------------------------- 417 : C ON A PLASTIFIE 418 : C----------------------------------------------------------------------- 419 : cri0= si * 1.e-8 420 : PHI0=PHI 421 : SI0=SI 422 : RR=0.D0 423 : 424 : DO I=1,NCONT 425 : DSIG(I)=0.D0 426 : DEPSP(I)=0.D0 427 : DSPHER1(I)=0.D0 428 : DSPHER2(I)=0.D0 429 : ENDDO 430 : 431 : C----------------------------------------------------------------------- 432 : C DEBUT DE LA BOUCLE D'ITERATIONS INTERNES 433 : C----------------------------------------------------------------------- 434 : sx1in=sx(1) 435 : sx2in=sx(2) 436 : sx3in= sx(3) 437 : s1in=s(1) 438 : s2in=s(2) 439 : 440 : iderin=0 441 : 10 NITER=NITER+1 442 : 443 : phi= seq - si 444 : c* if(itracb.eq.1) then 445 : c* write(6,*) ' niter seq si phi ' , niter,seq,si ,phi 446 : c* endif 447 : 448 : C---------CALCUL DE W1=DF/D(SIGMA)-------------------------------------- 449 : 450 : C---------ELEMENTS MASSIFS---------------------------------------------- 451 : 452 : IF(MFR1.EQ.1) THEN 453 : 454 : F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0 455 : F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0 456 : F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0 457 : DO I=4,NSTRSS 458 : F(I)=SX(I) 459 : ENDDO 460 : DO I=1,3 461 : W1(I)=1.5D0*F(I)/SEQ 462 : Z1(I)=W1(I) 463 : ENDDO 464 : DO I=4,NCONT 465 : W1(I)=3.D0*F(I)/SEQ 466 : Z1(I)=1.5D0*F(I)/SEQ 467 : ENDDO 468 : c* if(itracb.eq.1) then 469 : c* write(6,*)' W1' , w1(1), w1(2), w1(3), w1(4) 470 : c* endif 471 : 472 : C---------COQUES MINCES------------------------------------------------- 473 : 474 : ELSE IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN 475 : 476 : IF(IFOURB.GE.1) THEN 477 : W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1) 478 : W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1) 479 : W1(3)=3.D0*SX(3)/(SEQ*EP1) 480 : W1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(5))/(SEQ*EP1*EP1) 481 : W1(5)=3.D0*WORK(2)*(2.D0*SX(5)-SX(4))/(SEQ*EP1*EP1) 482 : W1(6)=18.D0*WORK(2)*SX(6)/(SEQ*EP1*EP1) 483 : ELSE 484 : W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1) 485 : W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1) 486 : W1(3)=3.D0*WORK(2)*(2.D0*SX(3)-SX(4))/(SEQ*EP1*EP1) 487 : W1(4)=3.D0*WORK(2)*(2.D0*SX(4)-SX(3))/(SEQ*EP1*EP1) 488 : ENDIF 489 : 490 : 491 : C---------COQUES EPAISSES----------------------------------------------- 492 : ELSE IF(MFR1.EQ.5) THEN 493 : F(1)=(2.D0*SX(1)-SX(2))/3.D0 494 : F(2)=(2.D0*SX(2)-SX(1))/3.D0 495 : DO I=3,5 496 : F(I)=SX(I) 497 : ENDDO 498 : DO I=1,2 499 : W1(I)=1.5D0*F(I)/SEQ 500 : Z1(I)=W1(I) 501 : ENDDO 502 : DO I=3,5 503 : W1(I)=3.D0*F(I)/SEQ 504 : Z1(I)=1.5D0*F(I)/SEQ 505 : ENDDO 506 : 507 : C---------POUTRES------------------------------------------------------- 508 : 509 : ELSE IF(MFR1.EQ.7) THEN 510 : 511 : DO I=1,NCONT 512 : DO J=1,NCONT 513 : DDA(I,J)=0.D0 514 : ENDDO 515 : ENDDO 516 : DDA(1,1)=YUNG 517 : DDA(4,4)=0.5*YUNG/(1.+XNU) 518 : DDA(5,5)=YUNG 519 : DDA(6,6)=YUNG 520 : W1(1)=SX(1)/SEQ 521 : W1(2)=0.D0 522 : W1(3)=0.D0 523 : W1(4)=SX(4)/SEQ 524 : W1(5)=SX(5)/SEQ 525 : W1(6)=SX(6)/SEQ 526 : ENDIF 527 : 528 : IF(MFR1.EQ.7) THEN 529 : DO I=1,NCONT 530 : r_z = 0.D0 531 : DO J=1,NCONT 532 : r_z=r_z+DDA(I,J)*W1(J) 533 : ENDDO 534 : W2(I)=r_z 535 : ENDDO 536 : ELSE 537 : IF(mfr.eq.1.and. ifourb.eq.-2) then 538 : DO I=1,NCONT 539 : r_z = 0.D0 540 : DO J=1,NCONT 541 : r_z=r_z+DDA(I,J)*W1(J) 542 : ENDDO 543 : W2(I)=r_z 544 : ENDDO 545 : else 546 : DO I=1,NCONT 547 : r_z = 0.D0 548 : DO J=1,NCONT 549 : r_z=r_z+DDAUX(I,J)*W1(J) 550 : ENDDO 551 : W2(I)=r_z 552 : ENDDO 553 : endif 554 : ENDIF 555 : COEF=0.D0 556 : DO I=1,NCONT 557 : COEF=COEF+W1(I)*W2(I) 558 : ENDDO 559 : 560 : C----------------------------------------------------------------------- 561 : C PLASTICITE PARFAITE, ECROUISSAGE ISOTROPE ET CINEMATIQUE ZIEGLER 562 : C----------------------------------------------------------------------- 563 : 564 : IF(INPLAS.EQ.1.OR.INPLAS.EQ.4.OR.INPLAS.EQ.5 565 : $ .OR.INPLAS.EQ.87) THEN 565 : CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC,1,IBI) 566 : * write(6,*) ' pente epst', pente,epst 567 : IF(IBI.EQ.1) THEN 568 : KERRE=75 569 : GOTO 1000 570 : ENDIF 571 : 572 : IF(INPLAS.EQ.1) THEN 573 : RP=0.D0 574 : C=0.D0 575 : ELSE IF(INPLAS.EQ.4) THEN 576 : RP=0.D0 577 : C=PENTE 578 : ELSE IF(INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN 579 : RP=PENTE 580 : C=0.D0 581 : ENDIF 582 : 583 : DENOM=COEF+C+RP 584 : DELTA=PHI/DENOM 585 : DMU=C*DELTA/SEQ 586 : * if(itracb.eq.1) then 587 : * write(6,*)'pente coef denom delta dmu ' 588 : * write(6,*)pente ,coef , denom ,delta ,dmu 589 : * endif 590 : 591 : DO I=1,NCONT 592 : DSIG(I)=-DELTA*W2(I) 593 : DSPHER1(I)=DMU*SX(I) 594 : ENDDO 595 : 596 : * Cas des contraintes planes en massif 597 : if(mfr.eq.1.and.ifourb.eq.-2) then 598 : 599 : * if(itracb.eq.1) then 600 : * write(6,*)'niter dsig delta',niter,dsig(1),dsig(2),dsig(3),delta 601 : * endif 602 : 603 : bb= abs(dsig(3)+ s(3) ) 604 : * if(itracb.eq.1) then 605 : * write(6,*) ' iderin bb ' , iderin , bb 606 : * endif 607 : r_z = dsig(3) * cte_cp 608 : sx(3)= sx3in - dsig(3) 609 : sx(1)= sx1in - r_z 610 : sx(2)= sx2in - r_z 611 : SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) 612 : s(3)= - dsig(3) 613 : s(1)= s1in - r_z 614 : s(2)= s2in - r_z 615 : if( bb.gt.cri0) then 616 : if(iderin.eq.0) then 617 : niter=niter - 1 618 : endif 619 : * if(itracb.eq.1) then 620 : * write(6,* ) ' cri0 bb seq delta ',cri0,bb,seq , delta 621 : * endif 622 : iderin=iderin+1 623 : ** write(6,*) ' iderin ' , iderin 624 : if(iderin.gt.50) then 625 : write(ioimp,*) ' probleme dans iterations internes' 626 : KERRE=2 627 : GO TO 1000 628 : endif 629 : go to 10 630 : endif 631 : DMU=C*DELTA/SEQ 632 : DO I=1,NCONT 633 : DSPHER1(I)=DMU*SX(I) 634 : ENDDO 635 : endif 636 : iderin=0 637 : DP=DELTA 638 : DR=RP*DP 639 : 640 : ELSE 641 : 642 : C----------------------------------------------------------------------- 643 : C MODELE DE CHABOCHE 644 : C----------------------------------------------------------------------- 645 : 646 : C---------UNIQUEMENT POUR LES ELEMENTS MASSIFS-------------------------- 647 : 648 : XPRO1=0.D0 649 : XPRO2=0.D0 650 : DO I=1,NCONT 651 : XPRO1=XPRO1+W1(I)*SPHER1(I) 652 : XPRO2=XPRO2+W1(I)*SPHER2(I) 653 : ENDDO 654 : 655 : FIP=1.D0+(PSI-1.D0)*EXP(-OME*EPST) 656 : 657 : DENOM=COEF+(A1*C1+A2*C2)*FIP-C1*XPRO1-C2*XPRO2+B*RMmRR 658 : DELTA=PHI/DENOM 659 : 660 : DO I=1,NCONT 661 : DSIG(I)=-DELTA*W2(I) 662 : DSPHER1(I)=(2.D0*A1*FIP*Z1(I)/3.D0-SPHER1(I))*C1*DELTA 663 : DSPHER2(I)=(2.D0*A2*FIP*Z1(I)/3.D0-SPHER2(I))*C2*DELTA 664 : ENDDO 665 : 666 : DR=B* RMmRR *DELTA 667 : DP=DELTA 668 : ENDIF 669 : 670 : RR=RR+DR 671 : EPST=EPST+DP 672 : 673 : IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN 674 : IF(IFOURB.GE.1) THEN 675 : DO I=1,3 676 : DSIG(I)=DSIG(I)*EP1 677 : ENDDO 678 : DO I=4,6 679 : DSIG(I)=EP1*EP1*DSIG(I)/(6.D0*WORK(2)) 680 : ENDDO 681 : ELSE 682 : DSIG(1)=DSIG(1)*EP1 683 : DSIG(2)=DSIG(2)*EP1 684 : DSIG(3)=EP1*EP1*DSIG(3)/(6.D0*WORK(2)) 685 : DSIG(4)=EP1*EP1*DSIG(4)/(6.D0*WORK(2)) 686 : ENDIF 687 : ENDIF 688 : C mise à jour des contraintes 689 : DO I=1,NCONT 690 : S(I)=S(I)+DSIG(I) 691 : SPHER1(I)=SPHER1(I)+DSPHER1(I) 692 : SPHER2(I)=SPHER2(I)+DSPHER2(I) 693 : SPHER(I)=SPHER1(I)+SPHER2(I) 694 : SX(I)=S(I)-SPHER(I) 695 : ENDDO 696 : if(ifourb.eq.-2.and.mfr.eq.1) then 697 : s(3)=0.d0 698 : endif 699 : * if(igau.eq.1) then 700 : * write(6,*) ' niter',niter 701 : * write(6,*) ' dsig' , dsig(1),dsig(2),dsig(3) 702 : * write(6,*) ' s' , s(1),s(2),s(3) 703 : * write(6,*) ' sx' , sx(1),sx(2),sx(3) 704 : * write(6,*) ' spher' , spher(1),spher(2),spher(3) 705 : * endif 706 : 707 : SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) 708 : 709 : C---------CONTRAINTES PLANES-------------------------------------------- 710 : 711 : IF(IFOURB.EQ.-2) THEN 712 : 713 : IF(MFR1.EQ.1) THEN 714 : F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0 715 : F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0 716 : F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0 717 : DO I=4,NSTRSS 718 : F(I)=SX(I) 719 : ENDDO 720 : DO I=1,3 721 : W1(I)=1.5D0*F(I)/SEQ 722 : ENDDO 723 : DO I=4,NSTRSS 724 : W1(I)=3.D0*F(I)/SEQ 725 : ENDDO 726 : 727 : ELSE IF(MFR1.EQ.3.OR.MFR1.EQ.9) THEN 728 : AUX=EP1*EP1*EP1*EP1 729 : W1(1)=(2.D0*SX(1)-SX(2))/(2.D0*SEQ*EP1*EP1) 730 : W1(2)=(2.D0*SX(2)-SX(1))/(2.D0*SEQ*EP1*EP1) 731 : W1(3)=18.D0*ALFAH*(2.D0*SX(3)-SX(4))/(SEQ*AUX) 732 : W1(4)=18.D0*ALFAH*(2.D0*SX(4)-SX(3))/(SEQ*AUX) 733 : ENDIF 734 : 735 : DO I=1,NSTRSS 736 : DEPSP(I)=DEPSP(I)+DELTA*W1(I) 737 : ENDDO 738 : ENDIF 739 : 740 : C----------------------------------------------------------------------- 741 : C TEST 742 : C CALCUL DE LA NOUVELLE VALEUR DE PHI 743 : C----------------------------------------------------------------------- 744 : IF(INPLAS.EQ.5.OR.INPLAS.EQ.87) THEN 744 : CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI) 745 : C* Modele de CHABOCHE (prise en compte ecrouissage) 746 : ELSE IF (INPLAS.EQ.7 .OR. INPLAS.EQ.12) THEN 747 : RMmRR = (RM - R0) * EXP(-B*EPST) 748 : SI = RM - RMmRR 749 : ELSE 750 : SI=RR+SI0 751 : ENDIF 752 : PHI=SEQ-SI 753 : 754 : PETI=1.E-7 755 : APHI=ABS(PHI) 756 : APHI0=ABS(PHI0) 757 : TEST=PETI*APHI0 758 : 759 : IF(NITER.GT.50) THEN 760 : if(itracb.eq.0) then 761 : itracb=1 762 : go to 2222 763 : endif 764 : KERRE=2 765 : GO TO 1000 766 : ENDIF 767 : * if(igau.eq.1) then 768 : * write(6 ,*)'test conver niter phi seq si testphi0',niter, 769 : * $ phi,seq,si,test,aphi0 770 : * write(6,*) ' rr si0' , rr, si0 771 : * endif 772 : IF(APHI.LE.TEST) THEN 773 : 774 : IF(MFR1.EQ.7) THEN 775 : DO I=1,NCONT 776 : S(I)=S(I)/DIV(I) 777 : ENDDO 778 : ENDIF 779 : 780 : C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES------------------- 781 : 782 : IF(IFOURB.NE.-2) THEN 783 : DO I=1,NCONT 784 : DS(I)=S(I)-SIGB(I) 785 : ENDDO 785 : CALL EPSIG0(DS,DDEPSE,MFR1,IFOURB,YUNG,XNU,WORK,NSTRSS) 786 : DO I=1,NCONT 787 : DEPSE(I)=DEPST(I)+DDEPSE(I) 788 : DEPSP(I)=DEPST(I)-DEPSE(I) 789 : ENDDO 790 : ENDIF 791 : 792 : DO I=1,NSTRSS 793 : SIGF(I)=S(I) 794 : DEFP(I)=DEPSP(I) 795 : ENDDO 796 : 797 : C---------COQUES AVEC CISAILLEMENT TRANSVERSE--------------------------- 798 : 799 : IF(MFR1.EQ.9) THEN 800 : DEFP(7)=0.D0 801 : DEFP(8)=0.D0 802 : SIGF(7)=SIGB(7) 803 : SIGF(8)=SIGB(8) 804 : ENDIF 805 : 806 : VARF(1)=EPST 807 : IF(INPLAS.EQ.4.OR.INPLAS.EQ.7.OR.INPLAS.EQ.12) THEN 808 : DO I=1,NSTRSS 809 : VARF(I+1)=SPHER1(I) 810 : ENDDO 811 : IF(INPLAS.EQ.12) THEN 812 : DO I=1,NSTRSS 813 : VARF(NSTRSS+1+I)=SPHER2(I) 814 : ENDDO 815 : ENDIF 816 : ENDIF 817 : IF(INPLAS.EQ.87) THEN 818 : DO I=1,NSTRSS 819 : VARF(1+I)=DEPST(I) 820 : ENDDO 821 : * write(6,*) ' sortie ecoin0' 822 : * write(6,*) 'depst',(depst(I),i=1,nstrs) 823 : * write(6,*) 'depsp',(depsp(i),i=1,nstrs) 824 : ENDIF 825 : KERRE=0 826 : RETURN 827 : 828 : ELSE 829 : sx1in=sx(1) 830 : sx2in=sx(2) 831 : s1in=s(1) 832 : s2in=s(2) 833 : * if(igau.eq.1) then 834 : * write(6,*) ' fin iter', niter,' nouvel sx',sx(1),sx(2),sx(3) 835 : * endif 836 : GOTO 10 837 : ENDIF 838 : C ENDIF 839 : 840 : 1000 CONTINUE 841 : RETURN 842 : END 843 : 844 :
© Cast3M 2003 - All rights reserved.
Disclaimer