1 : C SYCO12 SOURCE YURI 11/12/09 21:15:01 7199 2 : C SYCO12 SOURCE FANDEUR 11/04/12 21:15:09 6938 3 : 4 : SUBROUTINE SYCO12(wrk52,wrk53,wrk54,wrk2,wrk3, 5 : & IB,IGAU,NBPGAU,LTRAC,IFOURB,iecou,xecou) 6 : 7 : C----------------------------------------------------------------------- 8 : C ECOULEMENT PLASTIQUE POUR UN POINT 9 : C ALGORITHME ORTIZ ET SIMO 10 : C 11 : C EN ENTREE : 12 : C 13 : C SIG0 CONTRAINTES AU DEBUT DU PAS 14 : C DEPST INCREMENT DE DEFORMATIONS TOTALES 15 : C ( THERMIQUE ENLEVEE ) 16 : C VAR0 VARIABLES INTERNES DEDUT DU PAS 17 : C VAREX0 VARIABLES EXTERNES DEBUT DU PAS 18 : C VAREXF VARIABLES EXTERNES FIN DU PAS 19 : C XMAT COEFFICIENTS DU MATERIAU 20 : C PRECIS PRECISION POUR ITERATIONS INTERNES 21 : C WORK DES CARACTERISTIQUES 22 : C TRAC COURBE DE TRACTION 23 : C MFR1 INDICE DE FORMULATION 24 : C NSTRSS NOMBRE DE CONTRAINTES CA2000 25 : C INPLAS NUMERO DU MODELE DE PLASTICITE 26 : C DDAUX = MATRICE DE HOOKE ELASTIQUE 27 : C CMATE = NOM DU MATERIAU 28 : C VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU 29 : C VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES 30 : C N2EL = NBRE D ELEMENTS DANS SEGMENT DE HOOKE 31 : C N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE 32 : C IFOU = OPTION DE CALCUL 33 : C IB = NUMERO DE L ELEMENT COURANT 34 : C IGAU = NUMERO DU POINT COURANT 35 : C EPAIST= EPAISSEUR 36 : C NBPGAU= NBRE DE POINTS DE GAUSS 37 : C MELE = NUMERO DE L ELEMENT FINI 38 : C NPINT = NBRE DE POINTS D INTEGRATION 39 : C NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES 40 : C NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES 41 : C SECT = SECTION 42 : C LHOOK = TAILLE DE LA MATRICE DE HOOKE 43 : C TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES 44 : C UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE 45 : C 46 : C EN SORTIE : 47 : C 48 : C SIGF CONTRAINTES A LA FIN DU PAS 49 : C VARF VARIABLES INTERNES FIN DU PAS 50 : C DEFP INCR. DE DEFORMATIONS PLASTIQUES 51 : C KERRE CODE D'ERREUR 52 : C = 0 SI TOUT OK 53 : C = 99 SI FORMULATION NON DISPONIBLE 54 : C EN PLASTICITE 55 : C = 1 SI DEPS EST NEGATIF 56 : C = 2 SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE 57 : C 58 : C IFOUR : OPTION DE CALCUL 59 : C 60 : C IFOUR = -3 DEFORMATION PLANE GENERALISEE 61 : C IFOUR = -2 CONTRAINTES PLANES 62 : C IFOUR = -1 DEFORMATIONS PLANES 63 : C IFOUR = 0 AXISYMETRIQUE 64 : C IFOUR = 1 SERIE DE FOURIER 65 : C IFOUR = 2 TRIDIMENSIONNEL 66 : C 67 : C MFR1 : NUMERO DE LA FORMULATION ELEMENT FINI 68 : C 69 : C MFR1 = 1 MASSIF 70 : C MFR1 = 63 XFEM 71 : C 72 : C ISOTROPE 73 : C INPLAS : NUMERO DU MODELE DE VISCOPLASTICITE 74 : C 75 : C INPLAS = 153 SYCO1 76 : C INPLAS = 154 SYCO2 77 : C 78 : C----------------------------------------------------------------------- 79 : C CONVENTION DE REMPLISSAGE DES MEMOIRES 80 : C----------------------------------------------------------------------- 81 : C 82 : C XMAT(1) MODULE D'YOUNG 83 : C XMAT(2) COEFFICIENT DE POISSON 84 : C 85 : C TRAC(1 A 2*LTRAC) COURBE DE TRACTION 86 : C 87 : C MODELE ISOTROPE 88 : C VAR0(1) EPSE 89 : C VAR0(2) VP 90 : C 91 : C 92 : C----------------------------------------------------------------------- 93 : 94 : IMPLICIT INTEGER(I-N) 95 : IMPLICIT REAL*8(A-H,O-Z) 96 : 97 : -INC CCOPTIO 98 : -INC DECHE 99 : 100 : SEGMENT IECOU 101 : INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7, 102 : 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16, 103 : 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24, 104 : 3 icow25,icow26,icow27,icow28,icow29,icow30,icow31, 105 : 4 icow32,icow33,NSTRSS,MFR1, NBGMAT,NELMAT,icow38, 106 : 5 icow39,icow40,icow41,icow42,icow43,icow44 107 : INTEGER icow45,icow46,icow47,icow48,icow49,icow50, 108 : . icow51,icow52,icow53,icow54,icow55,icow56 109 : . icow57,icow58 110 : ENDSEGMENT 111 : * 112 : SEGMENT WRK2 113 : REAL*8 TRAC(LTRAC) 114 : ENDSEGMENT 115 : * 116 : SEGMENT WRK3 117 : REAL*8 WORK(LW),WORK2(LW2) 118 : ENDSEGMENT 119 : * 120 : * Commun XECOU: sert de fourre-tout pour les initialisations 121 : * de réels 122 : * 123 : SEGMENT XECOU 124 : * COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00 125 : REAL*8 DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00 126 : ENDSEGMENT 127 : C COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00 128 : C----------------------------------------------------------------------- 129 : 130 : DIMENSION SIG(130),EPS(130) 131 : DIMENSION S(8),SX(8),DS(8),DSIG(8),SPHER(8),SPHER1(8),SPHER2(8) 132 : DIMENSION DSPHER1(8),DSPHER2(8),DEPSE(8),DEPSP(8),DDEPSE(8) 133 : DIMENSION F(8),W1(8),W2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8) 134 : DIMENSION CRIGI(12) 135 : 136 : NCONT=NSTRSS 137 : C 138 : C--------PAS DE TEMPS ----------------------- 139 : TSYC9 = DT 140 : itracb=1 141 : 2222 continue 142 : 143 : C----------------------------------------------------------------------- 144 : DO I=1,NSTRSS 145 : S(I)=0.D0 146 : SPHER(I)=0.D0 147 : SPHER1(I)=0.D0 148 : SPHER2(I)=0.D0 149 : ENDDO 150 : KERRE=0 151 : YUNG=XMAT(1) 152 : XNU=XMAT(2) 153 : 154 : C---------CARACTERISTIQUES GEOMETRIQUES--------------------------------- 155 : C 156 : C 157 : EP1=1.D0 158 : ALFAH=1.D0 159 : 160 : IF(INPLAS.EQ.153) THEN 161 : C write(6,*) 'INPLAS', inplas 162 : C write (6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5) 163 : C write (6,*) 'XMAT', xmat(6), xmat(7), xmat(8) 164 : PSYC9=XMAT(6) 165 : DSYC9=XMAT(7) 166 : ENDIF 167 : IF(INPLAS.EQ.154) THEN 168 : C write(6,*) 'INPLAS', inplas 169 : C write (6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5) 170 : C write (6,*) 'XMAT', xmat(6), xmat(7), xmat(8),xmat(9),xmat(10) 171 : PSYC9=XMAT(6) 172 : ASYC9=XMAT(7) 173 : BSYC9=XMAT(8) 174 : CSYC9=XMAT(9) 175 : ENDIF 176 : 177 : DO I=1,LTRAC 178 : SIG(I)=TRAC(2*I-1) 179 : EPS(I)=TRAC(2*I) 180 : ENDDO 181 : 182 : EPST=VAR0(1) 183 : C 184 : C VP=VAR0(2) 185 : C if(VP.LE.1.d-4) VP = 1.d-4 186 : C VP0 = 1.D-4 187 : C IF(INPLAS.EQ.153) THEN 188 : C deno1 = 1. + ( (VP0/DSYC9)**(1./PSYC9) ) 189 : C deno2 = 1. + ( (VP/DSYC9)**(1./PSYC9) ) 190 : C else if(INPLAS.EQ.154) then 191 : C HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) 192 : C deno1 = 1. + (HSYC9 * ((VP0)**(1./PSYC9))) 193 : C deno2 = 1. + (HSYC9 * ((VP)**(1./PSYC9))) 194 : C endif 195 : C 196 : 197 : C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- 198 : C-xxxx pour le materiau syco1 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- 199 : C DSYC9 parametre D de la loi Symonds Cowper 200 : C PSYC9 parametre P de la loi Symonds Cowper 201 : C VP valeur minimum de la vitesse de deformation equivalente 202 : IF(INPLAS.EQ.153) THEN 203 : VP0 = 1.D-4 204 : VP = VP0 205 : else if(INPLAS.EQ.154) then 206 : C VP valeur minimum de la vitesse de deformation equivalente 207 : VP0 = 1.D-4 208 : VP = VP0 209 : endif 210 : C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- 210 : CALL CALSIG(DEPST,DDAUX,NSTRSS,CMATE,VALMAT,VALCAR, 211 : & N2EL,N2PTEL,MFR1,IFOURB,IB,IGAU,EPAIST, 212 : & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR, 213 : & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD) 214 : 215 : IF(IRTD.NE.1) THEN 216 : KERRE=69 217 : GOTO 1000 218 : ENDIF 219 : 220 : IF (ifourb.eq.-2) then 221 : * en cas de contraintes planes on repasse en 3D 222 : do iu=1,6 223 : do ju=1,6 224 : DDA(iu,ju)=0.d0 225 : enddo 226 : enddo 227 : cte_cp = xnu / (1.d0 - xnu) 228 : aux= YUNG / ((1.d0+xnu)*(1.d0-2.d0*xnu)) 229 : aux1= aux * xnu 230 : aux2= aux * (1.d0-xnu) 231 : gege = Yung / (2.d0 *(1.d0 +xnu)) 232 : DDA(1,1)=AUX2 233 : DDA(2,1)=AUX1 234 : DDA(1,2)=AUX1 235 : DDA(2,2)=AUX2 236 : DDA(3,3)=aux2 237 : DDA(1,3)=aux1 238 : DDA(2,3)=aux1 239 : DDA(3,1)=aux1 240 : DDA(3,2)=aux1 241 : 242 : DDA(4,4)=gege 243 : DDA(5,5)=gege 244 : DDA(6,6)=GEGE 245 : ENDIF 246 : * 247 : DO I=1,NSTRSS 248 : S(I)=SIG0(I)+DSIGT(I) 249 : SIGB(I)=S(I) 250 : SX(I)=S(I)-SPHER(I) 251 : ENDDO 252 : 253 : 254 : 255 : C----------------------------------------------------------------------- 256 : C CALCUL DE LA LIMITE ELASTIQUE SI 257 : C----------------------------------------------------------------------- 258 : 259 : C 259 : CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI) 260 : c 261 : IF(IBI.EQ.1) THEN 262 : KERRE=75 263 : GOTO 1000 264 : ENDIF 265 : 266 : C----------------------------------------------------------------------- 267 : C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ 268 : C----------------------------------------------------------------------- 269 : * attention en contraintes planes on se declare en 3D 270 : * rien besoin de faire dans vonmis0 car ifourb n'est pas utilisé 271 : * et les contraintes sont dimensionnées à 6 272 : * if(igau.eq.1) then 273 : * write(6,*) ' en entrée sx ',sx(1),sx(2),sx(3) 274 : * endif 275 : C 276 : C SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) 277 : C calcul de la contrainte de Von Mises 278 : SOM1=0.D0 279 : DO I=4,NSTRSS 280 : SOM1=SOM1+SX(I)*SX(I) 281 : ENDDO 282 : 283 : SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3) 284 : & -SX(2)*SX(3)+3.D0*SOM1 285 : SEQ=SQRT(ABS(SEQ)) 286 : C 287 : C----------------------------------------------------------------------- 288 : 289 : C LE CRITERE EST-IL VERIFIE 290 : C----------------------------------------------------------------------- 291 : C if(itracb.eq.1) then 292 : C write(6,*)' en entrée si seq', si, seq 293 : C endif 294 : C 295 : c test de plasticite sur SI statique 296 : C SI = SI * deno2/deno1 297 : C 298 : PHI=SEQ-SI 299 : NITER=0 300 : PETI=1.1D0*0.5D0*PRECIS*SEQ 300 : CALL EPSPRE(SEQ,SI,PETI,ITRY) 301 : C 302 : IF((ITRY.EQ.1).OR.(SEQ.LE.SI)) THEN 303 : * rien a faire on n'a pas plastifié 304 : DO I=1,NCONT 305 : SIGF(I)=S(I) 306 : DEFP(I)=0.D0 307 : ENDDO 308 : 309 : VARF(1)=VAR0(1) 310 : VARF(2)=VP0 311 : RETURN 312 : ENDIF 313 : 314 : C----------------------------------------------------------------------- 315 : C ON A PLASTIFIE 316 : C----------------------------------------------------------------------- 317 : C 318 : cri0= si * 1.d-8 319 : PHI0=PHI 320 : SI0 = SI 321 : RR=0.D0 322 : 323 : DO I=1,NCONT 324 : DSIG(I)=0.D0 325 : DEPSP(I)=0.D0 326 : DSPHER1(I)=0.D0 327 : DSPHER2(I)=0.D0 328 : ENDDO 329 : 330 : C----------------------------------------------------------------------- 331 : C DEBUT DE LA BOUCLE D'ITERATIONS INTERNES 332 : C----------------------------------------------------------------------- 333 : sx1in=sx(1) 334 : sx2in=sx(2) 335 : sx3in= sx(3) 336 : s1in=s(1) 337 : s2in=s(2) 338 : 339 : iderin=0 340 : 10 NITER=NITER+1 341 : 342 : phi= seq - si 343 : 344 : 345 : C---------CALCUL DE W1=DF/D(SIGMA)-------------------------------------- 346 : 347 : C---------ELEMENTS MASSIFS---------------------------------------------- 348 : 349 : 350 : F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0 351 : F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0 352 : F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0 353 : DO I=4,NSTRSS 354 : F(I)=SX(I) 355 : ENDDO 356 : DO I=1,3 357 : W1(I)=1.5D0*F(I)/SEQ 358 : Z1(I)=W1(I) 359 : ENDDO 360 : DO I=4,NCONT 361 : W1(I)=3.D0*F(I)/SEQ 362 : Z1(I)=1.5D0*F(I)/SEQ 363 : ENDDO 364 : 365 : 366 : if(ifourb.eq.-2) then 367 : DO I=1,NCONT 368 : r_z = 0.D0 369 : DO J=1,NCONT 370 : r_z=r_z+DDA(I,J)*W1(J) 371 : ENDDO 372 : W2(I)=r_z 373 : ENDDO 374 : else 375 : DO I=1,NCONT 376 : r_z = 0.D0 377 : DO J=1,NCONT 378 : r_z=r_z+DDAUX(I,J)*W1(J) 379 : ENDDO 380 : W2(I)=r_z 381 : ENDDO 382 : endif 383 : COEF=0.D0 384 : DO I=1,NCONT 385 : COEF=COEF+W1(I)*W2(I) 386 : ENDDO 387 : 388 : C----------------------------------------------------------------------- 389 : C ECROUISSAGE ISOTROPE 390 : C----------------------------------------------------------------------- 391 : 391 : CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC,1,IBI) 392 : * write(6,*) ' pente epst', pente,epst 393 : IF(IBI.EQ.1) THEN 394 : KERRE=75 395 : GOTO 1000 396 : ENDIF 397 : 398 : RP=PENTE 399 : C=0.D0 400 : 401 : DENOM=COEF+C 402 : C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- 403 : C---------CALCUL des termes en sup dus a la viscosité---------------- 404 : C il s'agit de -qsi1 = -d(phi)/d(°p) = +d(SY)/d(°p) 405 : c (où °p = dp/dt cad la vitesse inelastiq equivalente) 406 : IF(INPLAS.EQ.153) THEN 407 : deno1 = 1. + ( (VP0/DSYC9)**(1./PSYC9) ) 408 : deno2 = 1. + ( (VP/DSYC9)**(1./PSYC9) ) 409 : 410 : qsi1 = (SI/(PSYC9*DSYC9))*((VP/DSYC9)**( (1./PSYC9)-1.)) 411 : qsi1 = qsi1 / (deno2 * TSYC9) 412 : DENOM = DENOM + qsi1 413 : DENOM = DENOM + (RP * deno2/deno1) 414 : C 415 : ELSE IF (INPLAS.EQ.154) THEN 416 : HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) 417 : deno1 = 1. + (HSYC9 * ((VP0)**(1./PSYC9))) 418 : deno2 = 1. + (HSYC9 * ((VP)**(1./PSYC9))) 419 : 420 : qsi1 = (SI*HSYC9/PSYC9) * (VP**((1./PSYC9)-1.)) 421 : qsi1 = qsi1 / (deno2 * TSYC9) 422 : DENOM = DENOM + qsi1 423 : DENOM = DENOM + (RP * deno2/deno1) 424 : C 425 : qsi2 = (SI/deno2)*(VP**(1./PSYC9))*((ASYC9-HSYC9)/CSYC9) 426 : DENOM = DENOM + qsi2 427 : C 428 : ENDIF 429 : C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx- 430 : C23456789012345678901234567890123456789012345678901234567890123456789012 431 : DELTA=PHI/DENOM 432 : DMU=C*DELTA/SEQ 433 : C if(itracb.eq.1) then 434 : C write(6,*)'pente coef denom delta dmu niter si seq' 435 : C write(6,*)pente ,coef , denom ,delta ,dmu,niter,si,seq 436 : C endif 437 : 438 : DO I=1,NCONT 439 : DSIG(I)=-DELTA*W2(I) 440 : DSPHER1(I)=DMU*SX(I) 441 : ENDDO 442 : 443 : * Cas des contraintes planes en massif 444 : if(ifourb.eq.-2) then 445 : 446 : bb= abs(dsig(3)+ s(3) ) 447 : 448 : C if(itracb.eq.1) then 449 : C write(6,*) ' iderin bb ' , iderin , bb 450 : C endif 451 : 452 : r_z = dsig(3) * cte_cp 453 : sx(3)= sx3in - dsig(3) 454 : sx(1)= sx1in - r_z 455 : sx(2)= sx2in - r_z 456 : C 457 : C SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) 458 : C calcul de la contrainte de Von Mises 459 : SOM1=0.D0 460 : DO I=4,NSTRSS 461 : SOM1=SOM1+SX(I)*SX(I) 462 : ENDDO 463 : 464 : SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3) 465 : & -SX(2)*SX(3)+3.D0*SOM1 466 : SEQ=SQRT(ABS(SEQ)) 467 : C 468 : s(3)= - dsig(3) 469 : s(1)= s1in - r_z 470 : s(2)= s2in - r_z 471 : 472 : cri0= si * 1.d-8 473 : if( bb.gt.cri0) then 474 : 475 : if(iderin.eq.0) then 476 : niter=niter - 1 477 : endif 478 : 479 : C if(itracb.eq.1) then 480 : C write(6,* ) ' cri0 bb seq delta ',cri0,bb,seq , delta 481 : C endif 482 : 483 : iderin = iderin + 1 484 : 485 : ** write(6,*) ' iderin ' , iderin 486 : 487 : if(iderin.gt.50) then 488 : write(6,*) 'on iderin superieur a 50' 489 : write(ioimp,*) ' probleme dans iterations internes' 490 : KERRE=2 491 : GO TO 1000 492 : endif 493 : 494 : C23456789012345678901234567890123456789012345678901234567890123456789012 495 : C write(6,*)'cri0 bb niter iderin',cri0,bb,niter,iderin,dsig(3),s(3) 496 : go to 10 497 : endif 498 : 499 : DMU=C*DELTA/SEQ 500 : DO I=1,NCONT 501 : DSPHER1(I)=DMU*SX(I) 502 : ENDDO 503 : endif 504 : 505 : iderin=0 506 : C 507 : DR=RP*DELTA 508 : RR=RR+DR 509 : EPST=EPST+DELTA 510 : 511 : C mise à jour des contraintes 512 : DO I=1,NCONT 513 : S(I)=S(I)+DSIG(I) 514 : SPHER1(I)=SPHER1(I)+DSPHER1(I) 515 : SPHER2(I)=SPHER2(I)+DSPHER2(I) 516 : SPHER(I)=SPHER1(I)+SPHER2(I) 517 : SX(I)=S(I)-SPHER(I) 518 : ENDDO 519 : if(ifourb.eq.-2) then 520 : s(3)=0.d0 521 : endif 522 : C if(igau.eq.1) then 523 : C write(6,*) ' niter',niter 524 : C write(6,*) ' dsig' , dsig(1),dsig(2),dsig(3),dsig(4) 525 : C write(6,*) ' s' , s(1),s(2),s(3),s(4) 526 : C write(6,*) ' sx' , sx(1),sx(2),sx(3),sx(4) 527 : C write(6,*) ' spher' , spher(1),spher(2),spher(3),spher(4) 528 : C endif 529 : 530 : C SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH) 531 : C calcul de la contrainte de Von Mises 532 : SOM1=0.D0 533 : DO I=4,NSTRSS 534 : SOM1=SOM1+SX(I)*SX(I) 535 : ENDDO 536 : 537 : SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3) 538 : & -SX(2)*SX(3)+3.D0*SOM1 539 : SEQ=SQRT(ABS(SEQ)) 540 : C 541 : C---------CONTRAINTES PLANES-------------------------------------------- 542 : 543 : IF(IFOURB.EQ.-2) THEN 544 : 545 : F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0 546 : F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0 547 : F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0 548 : DO I=4,NSTRSS 549 : F(I)=SX(I) 550 : ENDDO 551 : DO I=1,3 552 : W1(I)=1.5D0*F(I)/SEQ 553 : ENDDO 554 : DO I=4,NSTRSS 555 : W1(I)=3.D0*F(I)/SEQ 556 : ENDDO 557 : 558 : 559 : DO I=1,NSTRSS 560 : DEPSP(I)=DEPSP(I)+DELTA*W1(I) 561 : ENDDO 562 : ENDIF 563 : 564 : C----------------------------------------------------------------------- 565 : C TEST 566 : C CALCUL DE LA NOUVELLE VALEUR DE PHI 567 : C----------------------------------------------------------------------- 567 : CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI) 568 : C****************************************************** 569 : c-xxxx calcul de la valeur (i+1) de H, de la vitesse de def VP, 570 : c et de SI en prenant en compte la viscosité 571 : IF(INPLAS.EQ.153) THEN 572 : C materiau SYCO1 573 : deno1 = 1. + ( (VP0/DSYC9)**(1./PSYC9) ) 574 : VP = MAX( (VP+(DELTA/TSYC9)) , VP0 ) 575 : C deno2 = 1. + ( (VP/DSYC9)**(1./PSYC9) ) 576 : SI = SI * (1. + ( (VP/DSYC9)**(1./PSYC9) )) / deno1 577 : C 578 : else if (INPLAS.EQ.154) then 579 : C materiau SYCO2 580 : HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9)) 581 : deno1 = 1. + (HSYC9 * ((VP0)**(1./PSYC9))) 582 : VP = MAX( (VP+(DELTA/TSYC9)) , VP0 ) 583 : C deno2 = 1. + (HSYC9 * ((VP)**(1./PSYC9))) 584 : SI = SI * (1. + (HSYC9 * (VP**(1./PSYC9)))) / deno1 585 : endif 586 : C****************************************************** 587 : 588 : 589 : PHI=SEQ-SI 590 : 591 : PETI=1.D-7 592 : APHI=ABS(PHI) 593 : APHI0=ABS(PHI0) 594 : TEST=PETI*APHI0 595 : 596 : 597 : IF(NITER.GT.200) THEN 598 : write(*,*) 'on a bien planté ici,car APHI=',APHI 599 : write(*,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI 600 : TEST2 = 1.E-6*APHI0 601 : if(APHI.gt.TEST2) then 602 : KERRE=2 603 : GO TO 1000 604 : else 605 : TEST=TEST2 606 : endif 607 : ELSEIF(NITER.eq.101) THEN 608 : write(*,*) 'on aurait du planté ici,car APHI=',APHI 609 : write(*,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI 610 : write(*,*) 'mais on va continuer 50pas de plus' 611 : ENDIF 612 : 613 : * if(igau.eq.1) then 614 : * write(6 ,*)'test conver niter phi seq si testphi0',niter, 615 : * $ phi,seq,si,test,aphi0 616 : * write(6,*) ' rr si0' , rr, si0 617 : * endif 618 : 619 : IF(APHI.LE.TEST) THEN 620 : 621 : 622 : C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES------------------- 623 : 624 : IF(IFOURB.NE.-2) THEN 625 : DO I=1,NCONT 626 : DS(I)=S(I)-SIGB(I) 627 : ENDDO 627 : C CALL EPSIG0(DS,DDEPSE,MFR1,IFOURB,YUNG,XNU,WORK,NSTRSS) 628 : C copie d'un bout de EPSIG0 629 : DO I=1,NSTRSS 630 : DDEPSE(I)=0.D0 631 : ENDDO 632 : AUX1=1.D0/YUNG 633 : AUX2=2.D0*(1.D0+XNU) 634 : AUX3=AUX1*AUX2 635 : DDEPSE(1)=(DS(1)-XNU*(DS(2)+DS(3)))*AUX1 636 : DDEPSE(2)=(DS(2)-XNU*(DS(1)+DS(3)))*AUX1 637 : DDEPSE(3)=(DS(3)-XNU*(DS(1)+DS(2)))*AUX1 638 : DO I=4,NSTRSS 639 : DDEPSE(I)=AUX3*DS(I) 640 : ENDDO 641 : C fin copie d'un bout de EPSIG0 642 : DO I=1,NCONT 643 : DEPSE(I)=DEPST(I)+DDEPSE(I) 644 : DEPSP(I)=DEPST(I)-DEPSE(I) 645 : ENDDO 646 : ENDIF 647 : 648 : DO I=1,NSTRSS 649 : SIGF(I)=S(I) 650 : DEFP(I)=DEPSP(I) 651 : ENDDO 652 : 653 : VARF(1)=EPST 654 : VARF(2)=VP 655 : 656 : KERRE=0 657 : RETURN 658 : 659 : ELSE 660 : sx1in=sx(1) 661 : sx2in=sx(2) 662 : s1in=s(1) 663 : s2in=s(2) 664 : 665 : GOTO 10 666 : ENDIF 667 : 668 : 1000 CONTINUE 669 : RETURN 670 : END 671 : 672 : 673 : 674 : 675 :
© Cast3M 2003 - All rights reserved.
Disclaimer