1 : C CFLUE2 SOURCE KICH 09/09/11 21:15:53 6505 2 : SUBROUTINE CFLUE2(wrk52,wrk53,wrk54,wrk2,wrk3, 3 : & IB,IGAU,NBPGAU,NBGMAT,NELMAT,IFOURB) 4 : * 5 : * modele `lemaitre fluage` Lemaitre et Chaboche pp.289,411 6 : * voir aussi syntheses RUPTHER SEMT/LISN/RT/99-124/A 7 : * remarque ep et sigma' sont colineaires donc ep et (Sigma0 + K depst)' le sont 8 : * Runge_Kutta d ordre 2 : Sigma1 = 1/2 YK0 + 1/2 YK1 9 : * YK0 = seq0 * (Sigma0 + K depst) / seqtot 10 : * etablit increment de r et D 11 : * YK1 = seq1 * (YK0 + K depst) / seqtot 12 : * etablit increment de r et D 13 : * puis moyenne 14 : IMPLICIT INTEGER(I-N) 15 : IMPLICIT REAL*8 (A-H,O-Z) 16 : -INC CCOPTIO 17 : -INC DECHE 18 : * 19 : * 20 : * 21 : SEGMENT WRK2 22 : REAL*8 TRAC(LTRAC) 23 : ENDSEGMENT 24 : * 25 : SEGMENT WRK3 26 : REAL*8 WORK(LW),WORK2(LW2) 27 : ENDSEGMENT 28 : * 29 : dimension YK0(8),YK1(8),spri0(8),delep0(8),spri1(8),delep1(8), 30 : &DIV(8) 31 : 32 : 33 : d0 = var0(3) 34 : r0 = var0(2) 35 : 36 : * cas isotrope 37 : ip1 = 4 38 : * 39 : youn0 = xmat0(1) 40 : sigy0 = xmat0(ip1+1) 41 : xn0 = xmat0(ip1+2) 42 : xm0 = xmat0(ip1+3) 43 : gk0 = xmat0(ip1+4) 44 : alp0 = xmat0(ip1+5) 45 : blp0 = xmat0(ip1+6) 46 : rini0 = xmat0(ip1+7) 47 : a0 = xmat0(ip1+8) 48 : pk0 = xmat0(ip1+9) 49 : x2mu0= xmat0(1)/(1.+xmat0(2)) 50 : 51 : if (d0.le.0.) d0 = 0.D0 52 : if (1. - d0.le.0.) then 53 : do ic = 1,nstrs 54 : sigf(ic) = 0. 55 : enddo 56 : 57 : depse = 0.d0 58 : do is = 1, nstrs 59 : defp(is) = depst(is) 60 : depse = depse + defp(is)*defp(is) 61 : enddo 62 : 63 : varf(1) = var0(1) + dsqrt(depse*0.5) 64 : varf(2) = var0(2) 65 : varf(3) = 1. 66 : goto 1002 67 : endif 68 : 69 : if (r0.lt.0.) then 70 : c write(6,*) 'erreur valeur r ', r0 71 : moterr(1:16) = conm 72 : moterr(17:24) = 'CFLUE2-3' 72 : call erreur(943) 73 : return 74 : endif 75 : if (r0.eq.0.) r0 = 1.e-10 76 : 77 : if(ib.eq.1.and.igau.eq.1) then 78 : * write(6,*) 't0' , youn0, sigy0, xn0 ,xm0 ,gk0,pk0 79 : endif 80 : 81 : youn1 = xmat(1) 82 : sigy1 = xmat(ip1+1) 83 : xn1 = xmat(ip1+2) 84 : xm1 = xmat(ip1+3) 85 : gk1 = xmat(ip1+4) 86 : alp1 = xmat(ip1+5) 87 : blp1 = xmat(ip1+6) 88 : rini1 = xmat(ip1+7) 89 : a1 = xmat(ip1+8) 90 : pk1 = xmat(ip1+9) 91 : x2mu1= xmat(1)/(1.+xmat(2)) 92 : 93 : if(ib.eq.1.and.igau.eq.1) then 94 : * write(6,*) 't1' , youn1, sigy1, xn1 ,xm1 ,gk1,pk1 95 : endif 96 : * 97 : delt = tempf - temp0 98 : * gerer la variable cachee avec les variations de temperature 99 : 100 : C---------CARACTERISTIQUES GEOMETRIQUES--------------------------------- 101 : C 102 : C COQUES 103 : C 104 : ALFAH=1.D0 105 : IF(MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9) THEN 106 : EP1=work(1) 107 : IF(MFR.NE.5) ALFAH=work(2)**2 108 : ENDIF 109 : C---------COQUES AVEC CT------------------------------------------------ 110 : C ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES 111 : 112 : IF(MFR.EQ.9) THEN 113 : NCONT=6 114 : ELSE 115 : NCONT=NSTRS 116 : ENDIF 117 : 118 : * calcul YK0 119 : * determine direction et sens de eplas 120 : * calcul increments de contrainte / utilise xmatf 121 : * remarque on utilise les caracteristiques elastiques a la date finale 122 : 123 : do ic = 1,valmat(/1) 124 : xmatf(ic) = valmat(ic) 125 : * completer pour le cas anisotrope 126 : if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1. - d0) 127 : enddo 128 : 128 : CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR, 129 : & N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST, 130 : & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR, 131 : & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD) 132 : 133 : IF(IRTD.NE.1) THEN 134 : KERRE=69 135 : GOTO 1010 136 : ENDIF 137 : 138 : DO I=1,NSTRS 139 : DSIGT(I)=SIG0(I) + dsigt(i) 140 : ENDDO 141 : 142 : C---------CAS DES POUTRES----------------------------------------------- 143 : 144 : IF(MFR.EQ.7) THEN 145 : DIV(1)=1.D0/work(4) 146 : DIV(2)=1.D0 147 : DIV(3)=1.D0 148 : DIV(4)=work(10)/work(1) 149 : DIV(5)=work(11)/work(2) 150 : DIV(6)=work(12)/work(3) 151 : IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4)) 152 : IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4)) 153 : IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4)) 154 : DO I=1,NCONT 155 : DSIGT(I)= DSIGT(I)*DIV(I) 156 : ENDDO 157 : ENDIF 158 : 159 : C----------------------------------------------------------------------- 160 : C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ 161 : C----------------------------------------------------------------------- 162 : seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH) 163 : 164 : if (seqtot.gt.0.) then 165 : seq0 = seqtot 166 : else 167 : * contrainte spherique 168 : do ic = 1,nstrs 169 : YK0(ic) = dsigt(ic) 170 : enddo 171 : 172 : do is = 1, nstrs 173 : delep0(is) = 0. 174 : enddo 175 : 176 : delr0 = 0. 177 : deld0 = 0. 178 : goto 500 179 : endif 180 : 181 : * point fixe pour determiner le multiplicateur de sigtot'/seqtot 182 : icaz = 1 183 : do ipfx = 1,10 184 : if(ib.eq.1.and.igau.eq.1) then 185 : c write(6,*) 'pt fixe' , ipfx, seq0,seqtot,icaz 186 : endif 187 : goto (70,80) icaz 188 : 70 continue 189 : * fonction 190 : if (gk1.gt.0..and.xm1.gt.0..and.r0.gt.0..and.(1.-d0).gt.0.) then 191 : delr0 = seq0/(1.-d0)/gk1/(r0 ** (1/xm1)) 192 : else 193 : delr0 = 0. 194 : endif 195 : delr0 = max(delr0,0.d0) 196 : c cas delr0 nul ? 197 : if (xn1.gt.0..and.x2mu1.gt.0.) then 198 : delr0 = delr0 ** xn1 199 : else 200 : moterr(1:16) = conm 201 : moterr(17:24) = 'CFLUE2-1' 201 : call erreur(943) 202 : return 203 : endif 204 : xmult = 1.5 * delr0 * x2mu1 * delt / (1. - d0) / (1. - d0) 205 : if(ib.eq.1.and.igau.eq.1) then 206 : c write(6,*) 'delr0' , delr0 207 : endif 208 : seq01 = seqtot - xmult 209 : goto 90 210 : 211 : * fonction associee 212 : 80 continue 213 : xmult = seqtot - seq0 214 : c en principe 1. - d0 > 0 215 : if (1. - d0.lt.0..or.delt.lt.0..or.x2mu1.le.0.) then 216 : * write(6,*) 'erreur parametres 1' , delt, x2mu1, d0 217 : moterr(1:16) = conm 218 : moterr(17:24) = 'CFLUE2-2' 218 : call erreur(943) 219 : return 220 : endif 221 : delr0 = xmult / delt * (1.- d0) * (1. - d0) /x2mu1 / 1.5 222 : if (xn1.gt.0) then 223 : delr0 = delr0 ** (1/xn1) 224 : * write(6,*) ' delr0 ', delr0 225 : else 226 : c write(6,*) 'erreur parametres 2' , xn1 227 : moterr(1:16) = conm 228 : moterr(17:24) = 'CFLUE2-4' 228 : call erreur(943) 229 : return 230 : endif 231 : if (gk1.gt.0..and.xm1.gt.0.) then 232 : seq01 = delr0 * (r0 ** (1/xm1)) * gk1 * (1. - d0) 233 : else 234 : c write(6,*) 'erreur parametres 3' , gk1,xm1 235 : moterr(1:16) = conm 236 : moterr(17:24) = 'CFLUE2-5' 236 : call erreur(943) 237 : return 238 : endif 239 : 240 : 90 continue 241 : if (ipfx.eq.1) then 242 : if (seq01.lt.0) then 243 : icaz = 2 244 : seq01 = seqtot/2. 245 : else if (seq01.eq.0.) then 246 : icaz =1 247 : seq01 = seqtot/2. 248 : endif 249 : endif 250 : varseq = abs((seq01 - seq0) / seq0) 251 : if (ib.eq.1.and.igau.eq.1) then 252 : c write(6,*) 'variation relative', varseq 253 : endif 254 : seq0 = seq01 255 : if (seq0 .lt.0.) then 256 : c write(6,*) 'erreur point fixe', seqtot, seq0 257 : moterr(1:16) = conm 258 : moterr(17:24) = 'CFLUE2-6' 258 : call erreur(943) 259 : return 260 : endif 261 : 262 : enddo 263 : 264 : 100 if (seqtot - seq0 .lt.0.) then 265 : c write(6,*) 'erreur point fixe', seqtot, seq0 266 : moterr(1:16) = conm 267 : moterr(17:24) = 'CFLUE2-7' 267 : call erreur(943) 268 : return 269 : endif 270 : 271 : *a revoir pour alp0 et beta0 . 272 : xki0 = seq0 273 : if ((xki0*a1).gt.0.and.(1. - d0).gt.0.. 274 : &and.pk1.gt.0.and.a1.ne.0..and.rini1.gt.0.)then 275 : deld0 = ((xki0 /a1) ** rini1) * ((1. - d0) ** (pk1 * (-1.))) 276 : else 277 : deld0 = 0. 278 : endif 279 : 280 : if(ib.eq.1.and.igau.eq.1) then 281 : * write(6,*) 'deld0' , deld0 282 : endif 283 : 284 : trsig0 = (dsigt(1) + dsigt(2) + dsigt(3)) * 0.3333333333 285 : spri0(1) = dsigt(1) - trsig0 286 : spri0(2) = dsigt(2) - trsig0 287 : spri0(3) = dsigt(3) - trsig0 288 : do is = 4,nstrs 289 : spri0(is) = dsigt(is) 290 : enddo 291 : 292 : if (deld0*delt.lt.1.) then 293 : 294 : if (seq0.gt.0.and.(1. - d0).gt.0..and.seqtot - seq0.gt.0.) then 295 : coep0 = (seqtot - seq0) * (1. - d0) / x2mu1/ seqtot 296 : else 297 : coep0 = 0.d0 298 : endif 299 : do is = 1, nstrs 300 : delep0(is) = spri0(is) * coep0 301 : enddo 302 : 303 : do ic = 1,nstrs 304 : YK0(ic) = spri0(ic) * seq0 / seqtot 305 : if (ic.le.3) YK0(ic) = YK0(ic) + trsig0 * seq0/seqtot 306 : enddo 307 : 308 : else 309 : 310 : do is = 1, nstrs 311 : delep0(is) = depst(is) 312 : enddo 313 : 314 : do ic = 1,nstrs 315 : YK0(ic) = 0. 316 : enddo 317 : endif 318 : 319 : 500 continue 320 : * calcul YK1 321 : * gerer la variable cachee avec les variations de temperature 322 : r1 = r0 + delr0 * delt 323 : if (r1.lt.0.) then 324 : c write(6,*) 'erreur valeur r ', r1 325 : moterr(1:16) = conm 326 : moterr(17:24) = 'CFLUE2-8' 326 : call erreur(943) 327 : return 328 : endif 329 : if (r1.eq.0.) r1 = 1.e-10 330 : 331 : d1 = var0(3) + deld0*delt 332 : if (d1.le.0.) d1 = 0.D0 333 : if (d1.ge.1.) then 334 : do ic = 1,nstrs 335 : YK1(ic) = 0. 336 : enddo 337 : deld1 = 0. 338 : delr1 = 0. 339 : do is = 1, nstrs 340 : delep1(is)= depst(is) 341 : enddo 342 : goto 1000 343 : endif 344 : 345 : do ic = 1,valmat(/1) 346 : xmatf(ic) = valmat(ic) 347 : * completer pour le cas anisotrope 348 : if (ic.eq.1) xmatf(ic) = xmatf(ic) * (1. - d1) 349 : enddo 350 : 350 : CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,xmatf,VALCAR, 351 : & N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST, 352 : & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR, 353 : & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD) 354 : 355 : IF(IRTD.NE.1) THEN 356 : KERRE=69 357 : GOTO 1010 358 : ENDIF 359 : 360 : DO I=1,NSTRS 361 : DSIGT(I)=SIG0(I) + dsigt(i) 362 : ENDDO 363 : C---------CAS DES POUTRES----------------------------------------------- 364 : 365 : IF(MFR.EQ.7) THEN 366 : DIV(1)=1.D0/work(4) 367 : DIV(2)=1.D0 368 : DIV(3)=1.D0 369 : DIV(4)=work(10)/work(1) 370 : DIV(5)=work(11)/work(2) 371 : DIV(6)=work(12)/work(3) 372 : IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4)) 373 : IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4)) 374 : IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4)) 375 : DO I=1,NCONT 376 : DSIGT(I)= DSIGT(I)*DIV(I) 377 : ENDDO 378 : ENDIF 379 : 380 : C----------------------------------------------------------------------- 381 : C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ 382 : C----------------------------------------------------------------------- 383 : seqtot=VONMIS0(dsigt,NSTRS,MFR,IFOURB,EP1,ALFAH) 384 : 385 : 386 : if (seqtot.gt.0.) then 387 : seq1 = seqtot 388 : else 389 : * contrainte spherique 390 : do ic = 1,nstrs 391 : YK1(ic) = dsigt(ic) 392 : enddo 393 : 394 : do is = 1, nstrs 395 : delep1(is) = 0. 396 : enddo 397 : 398 : delr1 = 0. 399 : deld1 = 0. 400 : goto 1000 401 : endif 402 : 403 : * point fixe pour determiner le multiplicateur de sigtot'/seqtot 404 : icaz = 1 405 : do ipfx = 1,10 406 : if(ib.eq.1.and.igau.eq.1) then 407 : c write(6,*) 'pt fixe' , ipfx, seq1,seqtot,icaz 408 : endif 409 : goto (570,580) icaz 410 : 411 : 570 continue 412 : * fonction 413 : if (gk1.gt.0..and.xm1.gt.0..and.r1.gt.0..and.(1.-d1).gt.0.) then 414 : delr1 = seq1/(1.-d1)/gk1/(r1 ** (1/xm1)) 415 : else 416 : delr1 = 0. 417 : endif 418 : delr1 = max(delr1,0.d0) 419 : c cas delr1 nul ? 420 : if (xn1.gt.0.) then 421 : delr1 = delr1 ** xn1 422 : else 423 : c write(6,*) 'erreur parametres 2' , xn1 424 : moterr(1:16) = conm 425 : moterr(17:24) = 'CFLUE2-9' 425 : call erreur(943) 426 : return 427 : endif 428 : if(ib.eq.1.and.igau.eq.1) then 429 : * write(6,*) 'delr1' , delr1 430 : endif 431 : xmult = 1.5 * delr1 * x2mu1 / (1.- d1) / (1. - d1) * delt 432 : seq11 = seqtot - xmult 433 : goto 590 434 : 435 : * fonction associee 436 : 580 continue 437 : xmult = seqtot - seq1 438 : c en principe 1. - d1 > 0 439 : if (1. - d1.lt.0..or.delt.lt.0..or.x2mu1.le.0.) then 440 : c write(6,*) 'erreur parametres 1' , delt, x2mu1, d1 441 : moterr(1:16) = conm 442 : moterr(17:24) = 'CFLUE210' 442 : call erreur(943) 443 : return 444 : endif 445 : delr1 = xmult / delt * (1.- d1) * (1. - d1) /x2mu1 / 1.5 446 : if (xn1.gt.0) then 447 : delr1 = delr1 ** (1/xn1) 448 : * write(6,*) ' delr1 ', delr1 449 : else 450 : c write(6,*) 'erreur parametres 2' , xn1 451 : moterr(1:16) = conm 452 : moterr(17:24) = 'CFLUE211' 452 : call erreur(943) 453 : return 454 : endif 455 : if (gk1.gt.0..and.xm1.gt.0.) then 456 : seq11 = delr1 * (r1 ** (1/xm1)) * gk1 * (1. - d1) 457 : else 458 : c write(6,*) 'erreur parametres 3' , gk1,xm1 459 : moterr(1:16) = conm 460 : moterr(17:24) = 'CFLUE212' 460 : call erreur(943) 461 : return 462 : endif 463 : 464 : 590 continue 465 : if (ipfx.eq.1) then 466 : if (seq11.lt.0) then 467 : icaz = 2 468 : seq11 = seqtot/2. 469 : else if (seq11.eq.0.) then 470 : icaz =1 471 : seq11 = seqtot/2. 472 : endif 473 : endif 474 : varseq = abs((seq11 - seq1) / seq1) 475 : if(ib.eq.1.and.igau.eq.1) then 476 : c write(6,*) 'variation relative', varseq 477 : endif 478 : seq1 = seq11 479 : if (seq1 .lt.0.) then 480 : c write(6,*) 'erreur point fixe', seqtot, seq1 481 : moterr(1:16) = conm 482 : moterr(17:24) = 'CFLUE213' 482 : call erreur(943) 483 : return 484 : endif 485 : 486 : 487 : enddo 488 : 489 : 490 : if (seqtot - seq1 .lt.0.) then 491 : c write(6,*) 'erreur point fixe', seqtot, seq1 492 : moterr(1:16) = conm 493 : moterr(17:24) = 'CFLUE214' 493 : call erreur(943) 494 : return 495 : endif 496 : 497 : 498 : *a revoir pour alp1 et bta1 . 499 : xki1 = seq1 500 : if ((xki1*a1).gt.0.and.(1. - d1).gt.0.. 501 : &and.pk1.gt.0.) then 502 : deld1 = ((xki1 /a1) ** rini1) * ((1. - d1) ** (pk1 * (-1.))) 503 : else 504 : deld1 = 0. 505 : endif 506 : if(ib.eq.1.and.igau.eq.1) then 507 : * write(6,*) 'deld1' , deld1 508 : endif 509 : 510 : trsig1 = (dsigt(1) + dsigt(2) + dsigt(3)) * 0.3333333333 511 : spri1(1) = dsigt(1) - trsig1 512 : spri1(2) = dsigt(2) - trsig1 513 : spri1(3) = dsigt(3) - trsig1 514 : do is = 4,nstrs 515 : spri1(is) = dsigt(is) 516 : enddo 517 : 518 : if (deld1*delt.lt.1) then 519 : if (seq1.gt.0.and.(1. - d1).gt.0..and.seqtot - seq1.gt.0.) then 520 : coep1 = (seqtot - seq1) * (1. - d1) / x2mu1/ seqtot 521 : else 522 : coep1 = 0. 523 : endif 524 : do is = 1, nstrs 525 : delep1(is) = spri1(is) * coep1 526 : enddo 527 : 528 : do ic = 1,nstrs 529 : YK1(ic) = spri1(ic) * seq1 / seqtot 530 : if (ic.le.3) YK1(ic) = YK1(ic) + trsig1 * seq1 / seqtot 531 : enddo 532 : 533 : else 534 : 535 : do is = 1, nstrs 536 : delep1(is) = depst(is) 537 : enddo 538 : 539 : do ic = 1,nstrs 540 : YK1(ic) = 0. 541 : enddo 542 : 543 : endif 544 : 545 : c au final 546 : 1000 continue 547 : do ic = 1,nstrs 548 : sigf(ic) = 0.5*YK0(ic) + 0.5*YK1(ic) 549 : enddo 550 : 551 : depse = 0.d0 552 : do is = 1, nstrs 553 : defp(is) = 0.5*(delep0(is) + delep1(is)) 554 : depse = depse + defp(is)*defp(is) 555 : enddo 556 : 557 : varf(1) = var0(1) + dsqrt(depse*0.5) 558 : vtd = var0(3) + (deld0 + deld1)*0.5*delt 559 : if (vtd.le.1.) then 560 : varf(3) = vtd 561 : else 562 : varf(3) = 1. 563 : endif 564 : varf(2) = (r0 + r1)*0.5 565 : 566 : 567 : 1002 continue 568 : if(ib.eq.1.and.igau.eq.1) then 569 : * write(6,*) 't0' ,sigf(3),yk0(3),yk1(3),varf(3),depst(3) 570 : endif 571 : 572 : 1010 continue 573 : RETURN 574 : END 575 : 576 : 577 : 578 : 579 :
© Cast3M 2003 - All rights reserved.
Disclaimer