1 : C SJONL3 SOURCE CHAT 06/06/01 21:20:08 5450 2 : SUBROUTINE SJONL3(sigi, depst, vari, xmat, tras, jtras , 3 : > trat, jtrat, trac, jtrac, sigf, varf, 4 : > depsp, ireturn) 5 : c --- SJONL3 --- 6 : c computing plastic flow in joint element 7 : c under 3D plane (strain) conditions 8 : c 9 : c beat release 24-7-00 ylp. 10 : c sigi IN initial stress 11 : c depst IN total strain increment 12 : c vari IN initial internal variable 13 : c xmat IN material properties 14 : c tras IN shear plastic softening behaviour 15 : c jtras INnumber of points 16 : c trat IN tensile plastic softening behaviour 17 : c jtrat IN number of points 18 : c trac IN compression plastic behaviour 19 : c jtrac IN number of points 20 : c sigf OUT final stress 21 : c varf OUT final internal variable 22 : c depsp OUT plastic strain increment 23 : c ireturn OUT routine control integer 24 : 25 : IMPLICIT REAL*8(A-H,K-Z) 26 : 27 : IMPLICIT INTEGER(I-J) 28 : dimension sigi(*), depst(*), vari(*), 29 : > sigf(*), depsp(*), varf(*) 30 : dimension ddepst(3), ddepsp(3) 31 : dimension xmat(*) 32 : dimension tras(2,jtras), trat(2,jtrat), trac(2,jtrac) 33 : c 1 --> strain 34 : c 2 --> stress 35 : c_______________________________________________________________________ 36 : c 37 : * save iloop 38 : * data iloop/0/ 39 : c print *,'debug -- sjonl3d -- check everything' 40 : c 41 : c do i=1,3 42 : c print *,'depst ',depst(i) 43 : c print *,'vari ', vari(i) 44 : c print *,'sigi ', sigi(i) 45 : c print *,'sigf ', sigf(i) 46 : c print *,'varf ',varf(i) 47 : c print *,'depst ',depst(i) 48 : c enddo 49 : c do i=1,6 50 : c print *,'xmat ', xmat(i) 51 : c enddo 52 : c print *,'tras ', tras ,' jtras ', jtras 53 : c print *,'trat ', trat ,' jtrat ', jtrat 54 : c print *,'trac ', trac ,' jtrac ', jtrac 55 : c 56 : c print *,'ireturn ', ireturn 57 : * iloop = iloop + 1 58 : c look at xmat(6) and deduce elastic/ plastic... 59 : c ...coupled/uncoupled variables. 60 : c all variables Jxxx are integer used as logical with the 61 : c following rule (0 is true, 1 is false) 62 : c JCP : use plastic law in compression 63 : c JSP : use plastic law in shear 64 : c JTP : use plastic law in tension 65 : c JST : couple shear and tension 66 : c JSC : couple shear and compression 66 : call sjoncpl(xmat(6),JCP,JSP,JTP,JST,JSC) 67 : * print *,'DEBUG --> entree dans sjonl2' 68 : * print *,'DepsN DepsS ' ,depst(2), depst(1) 69 : * set some usefull variables 70 : ZERO = 0.0d0 71 : ONE = 1.0d0 72 : c set the initial location of the edge of the cone 73 : ETi = xmat(5) 74 : c index definition 75 : Ipress = 2 76 : Ishear = 1 77 : Itotal = 3 78 : 79 : Ishea1 = 1 80 : Ishea2 = 2 81 : Inorm = 3 82 : Icompr = 4 83 : 84 : Istrain = 1 85 : Istress = 2 86 : c number of increments 87 : INCR = 5 88 : c______________________________________________________________________ 89 : c 90 : c elastic initialisation of the joint stiffness --> K 91 : c K[z] 92 : c with z = t tensile ; s shear ; c compression 93 : Kt0 = xmat(Ipress) 94 : Ks0 = xmat(Ishear) 95 : Kc0 = Kt0 96 : c print *,'Deps - tension ',depst(Ipress) 97 : c print *,'Deps - shear ',depst(Ishear) 98 : c split the displacement increment in few smaller increments 99 : ddepst(Inorm) = depst(Inorm)/INCR 100 : ddepst(Ishea1) = depst(Ishea1)/INCR 101 : ddepst(Ishea2) = depst(Ishea2)/INCR 102 : c STARTING MAIN LOOP 103 : c [0 104 : do iflag =1 , INCR 105 : c______________________________________________________________________ 106 : c 107 : c initial residual strength 108 : c RTS --> tensile 109 : c RSS --> shear 110 : c 111 : c internal variables initialisation --> Alpha 112 : c A[xyz] 113 : c with x = i initial ; f final ; k intermediate computational value 114 : c y = p plastic ; t total 115 : c z = t tensile ; s shear 116 : Aipt = vari(Ipress) 117 : c \--> initial plastic tensile strain 118 : Aips = vari(Ishear) 119 : c \--> initial plastic shear strain 120 : Aipc = vari(Icompr) 121 : c \--> initial plastic compression strain 122 : c______________________________________________________________________ 123 : c 124 : c initial residual strength 125 : c RTS --> tensile 126 : c RSS --> shear 127 : c RCS --> compression 128 : c compute the tensile stress RTS corresponding to the initial 129 : c plastic tensile strain with evolution 129 : call yofxcu(Aipt ,trat, jtrat, RTS, dRTS, ireturn) 130 : c \--> Aipt initial plastic tensile strain 131 : c RTS OUT 132 : c dRTS = dTkt/dAkps useless value 133 : c [1 134 : if (ireturn .ne. 0) return 135 : c \--> ERROR 136 : c 1] 137 : c compute the initial cohesion residual 137 : call yofxcu(Aips, tras, jtras, RSS, dRSS, ireturn) 138 : c [1 139 : if (ireturn .ne. 0) return 140 : c \--> ERROR 141 : c 1] 142 : c comput the initial compression residual 142 : call yofxcu(Aipc, trac, jtrac, RCS, dRCS, ireturn) 143 : c [1 144 : if (ireturn .ne. 0) return 145 : c \--> ERROR 146 : c 1] 147 : c______________________________________________________________________ 148 : c 149 : 150 : Aitt = vari(Itotal) 151 : c \--> initial total tensile strain 152 : c compute final total tensile strain... 153 : Aftt = Aitt + ddepst(Inorm) 154 : c ... and the initial tensile strength... 155 : TT0 = trat(Istress,1) 156 : c ... and the elastic compressive strength 157 : TC0 = trac(Istress,1) 158 : c print *,'TC0 ',TC0 159 : c check one or two things if entering sjonl2 for first time 160 : c [1 161 : * mis en commentaire par TC 162 : * if ( iloop .eq. 1) then 163 : c [2 164 : * if ( TT0 .eq. ZERO ) then 165 : * print *,'cut-off tensile strength is zero' 166 : * print *,'denied option' 167 : * return 168 : * endif 169 : c 2] 170 : c [2 171 : * if ( TT0 .gt. ETi ) then 172 : * print *,'WARNING cut-off ouside the friction cone...' 173 : * print *,'... tensile strength at the edge of the cone' 174 : * print *,'set to cut-off tensile strength' 175 : * return 176 : * endif 177 : c 2] 178 : * endif 179 : * fin de mise en commentaire par TC 180 : c 1] 181 : c initialization of the damage variables 182 : Dtensi = 1.0d0 183 : Dshear = 1.0d0 184 : Dcompr = 1.0d0 185 : c WARNING the damage variables are understand here as : 186 : c D = 1 means no damage 187 : c D = 0 means full damage... 188 : c which is exactly the opposite of the usual definition !!! 189 : c compute the predicted elastic stress state 190 : c using the secant stiffness 191 : c 192 : c first, look at the tension modulus... 193 : c maximal (initial) tensile stress initialisation 194 : c N.B.: it is thus assumed that the behaviour in tension is elastic 195 : c plastic, and that the peak value features the end of the elastic 196 : c regime. 197 : c check whether tension or compression region is reached 198 : c [1 199 : if ( Aftt .gt. ZERO ) then 200 : c \--> tensile region reached... 201 : c [2 202 : if ( JTP .ne. 0 ) then 203 : c \--> use elastic law 204 : Kt = Kt0 205 : else 206 : c \--> use damaged elastic predictor 207 : c [3 208 : if ( RTS .eq. ZERO ) then 209 : c \--> zero residual tensile strength 210 : Kt = ZERO 211 : c \--> tensile modulus set to zero 212 : else 213 : c \--> still got residual tensile strength 214 : Kt = RTS / ( Aipt + TT0 / Kt0 ) 215 : endif 216 : c 3] 217 : c update tensile damage variable to enable coupling effect 218 : c [3 219 : if ( JST .eq. 0) Dtensi = RTS / TT0 220 : c 3] 221 : endif 222 : c 2] 223 : c print *,'Dtensi ' ,Dtensi 224 : JCT = 0 225 : c JCT is an integer value used as a logical to control the location 226 : c of the stress state with respect to tension or compression: 227 : c tension --> JCT = 0 228 : c compression --> JCT = 1 229 : else 230 : c \--> compression region reached... 231 : c [2 232 : if ( JCP .ne. 0 ) then 233 : c \--> use elastic law 234 : Kc = Kc0 235 : else 236 : c \--> use damaged elastic predictor 237 : c [3 238 : if ( RCS .eq. ZERO ) then 239 : Kc = ZERO 240 : else 241 : Kc = RCS / (Aipc + TC0 / Kc0 ) 242 : endif 243 : c 3] 244 : c update compression damage variable 245 : c Dcompr = Kc / Kc0 --> too fast 246 : c ... try something else... 247 : endif 248 : c 2] 249 : JCT = 1 250 : endif 251 : c ]1 252 : c ... then, look at the shear modulus 253 : TS0 = tras(Istress,1) 254 : c [1 255 : if (JSP .ne. 0 ) then 256 : Ks = Ks0 257 : else 258 : c [2 259 : if ( RSS .eq. ZERO ) then 260 : c \--> zero residual shear strength 261 : c or opened bonds 262 : Ks = ZERO 263 : c \--> shear modulus set to zero 264 : c Such case shall only appear if the bonds are opened. 265 : c otherwise, residual cohesion is generally observed 266 : else 267 : c \--> still got residual shear strength 268 : Ks = RSS / (Aips + TS0 / Ks0) 269 : endif 270 : c update the shear damage variable 271 : if (JST .eq. 0) Dshear = RSS / TS0 272 : c(Ks / Ks0) 273 : c print *,'Dshear ' , Dshear 274 : c 2] 275 : endif 276 : c 1] 277 : 278 : c [1 279 : if (JCT .eq. 0) then 280 : c check if previous stress state was in compression 281 : c [2 282 : if ( (Aitt * Aftt) .lt. ZERO ) then 283 : Tft = Kt * Aftt * Dtensi 284 : else 285 : Tft = sigi(Inorm) + ( ddepst(Inorm) * Kt * Dshear) 286 : endif 287 : c 2] 288 : Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dtensi) 289 : Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dtensi) 290 : Tfs0 = sqrt((Tfs1**2)+(Tfs2**2)) 291 : else 292 : c check if previous stress state was in tension 293 : c [2 294 : if ( (Aitt * Aftt) .lt. ZERO ) then 295 : Tfc = Kc * Aftt * Dcompr 296 : else 297 : c Tfc = sigi(Ipress) + ( ddepst(Ipress) * Kc * Dshear) 298 : Tfc = sigi(Inorm) + ( ddepst(Inorm) * Kc ) 299 : endif 300 : c 2] 301 : Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dcompr) 302 : Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dcompr) 303 : Tfs0 = sqrt((Tfs1**2)+(Tfs2**2)) 304 : endif 305 : c 1] 306 : 307 : c print *,'stress state elastic prediction' 308 : c print *,'tension compression shear' 309 : c print *,Tft, Tfc, Tfs1, Tfs2 310 : c______________________________________________________________________ 311 : c 312 : c before checking if any kind of yielding occured, update the 313 : c final internal variables in the case of simple elastic regime. 314 : Afpt = Aipt 315 : Afps = Aips 316 : Afpc = Aipc 317 : c______________________________________________________________________ 318 : c 319 : c now, check yielding 320 : c N.B.: note that each mode might be activated independently 321 : c without any restriction!!!! 322 : c first, look at the tension/compression mode... 323 : c [1 324 : if ( JCT .eq. 0 ) then 325 : c [2 326 : if ( JTP .eq. 0 ) then 327 : c [3 328 : if ( abs(Tft) .gt. RTS*Dshear ) then 329 : c \--> YIELDING IN TENSION 330 : c print *,'yielding in tension...' 331 : c Afpt = Aftt - TT0 / Kt0 332 : Afpt = Aipt + (Tft-RTS*Dshear)/Kt 333 : c print *,'Afpt ',Afpt 334 : c \--> updated plastic tensile strain 335 : * print *,'updated plastic tensile strain ', Afpt 336 : c compute the updated residual tensile strength... 336 : call yofxcu(Afpt, trat, jtrat, RTS ,dRTS, ireturn) 337 : c \--> Afpt updated plastic tensile strain 338 : c RTS updated residual tensile strength 339 : c dTkt dummy 340 : Tft = RTS * Dshear 341 : c \--> compute the final tensile stress 342 : c print *,'residual tensile strength ',Tft 343 : endif 344 : c 3] 345 : endif 346 : c 2] 347 : else 348 : c [2 349 : if ( JCP .eq. 0 ) then 350 : c if ( abs(Tfc) .gt. abs(RCS)*Dshear) then 351 : c [3 352 : if ( abs(Tfc) .gt. abs(RCS) ) then 353 : c \--> YIELDING IN COMPRESSION 354 : c print *,'yielding in compression' 355 : c Afpc = abs(Aftt) - TC0 / Kc0 356 : c Afpc = Aipc + ((Tfc*-1.0d0)-RCS*Dshear)/Kc 357 : Afpc = Aipc + ((Tfc*(-1.0d0))-RCS)/Kc 358 : c \--> updated plastic compression strain 359 : c print *,'updated plastic compression strain ', Afpc 360 : c compute the updated residual compression strength 360 : call yofxcu(Afpc, trac, jtrac, RCS , dRCS, ireturn) 361 : c then, go back to the compression domain (negative stress) 362 : Tfc = RCS * (-1.0d0) 363 : c print *,'residual compressive strength ',RCS 364 : endif 365 : c 3] 366 : endif 367 : c 2] 368 : endif 369 : c 1] 370 : c ... finally, look at the shear mode... 371 : c compute the tensile value at the edge using an homothetic transform. 372 : c the cut-off value plays the role of a drag-stress 373 : c ETf = RTS + ETi - TT0 374 : c ... definitely not relevant!!!! 375 : c ... try something else... like a simple homothetic rule related to 376 : c the location of the residual cut-off with regard to the initial 377 : c location of the edge of the friction cone 378 : ETf = ETi * RTS / TT0 379 : c the behaviour in shear is assumed to follow a 380 : c homothetic rule with respect to the behaviour in pure shear 381 : c the coefficient of this homothetic transform is given by... 382 : c [1 383 : if ( JSP .eq. 0 ) then 384 : c [2 385 : if ( ETf .ne. ZERO ) then 386 : c [3 387 : if ( JCT .eq. 0 ) then 388 : c \--> shear yielding with tension 389 : Chom = (ETf - Tft)/(ETf)*Dtensi 390 : c print *,'Coef. hom. ', Chom 391 : c WARNING the original version used 392 : c Chom = (RTS - Tft)/(RTS +Pnor) 393 : c xmat(5) or PNOR is a simple numerical value used is order to 394 : c avoid division by zero if the residual cohesion vanishes. 395 : c print *,'???', Dtensi, ONE 396 : c [4 397 : if ( Dtensi .eq. ONE) then 398 : c \--> no tensile damage 399 : c [5 400 : if ( abs(Tfs0) .gt. (RSS*Chom) ) then 401 : c \--> YIELDING IN SHEAR 402 : c print *,'YIELDING IN SHEAR' 403 : c compute the total plastic shear strain... 404 : Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks 405 : c print *,'plasti shear strain ', Afps 406 : c ...and the final cohesion residual 406 : call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn) 407 : c [6 408 : if (ireturn .ne. 0) return 409 : c \--> ERROR 410 : c 6] 411 : Tfs = sign(Chom*RSS , Tfs) 412 : if (Tfs0 .eq. ZERO) then 413 : Tfs1 = ZERO 414 : Tfs2 = ZERO 415 : else 416 : Tfs1= Tfs*Tfs1/Tfs0 417 : Tfs2= Tfs*Tfs2/Tfs0 418 : endif 419 : c print *,'updated shear stress ',Tfs 420 : c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 421 : endif 422 : c 5] 423 : else 424 : c \--> already got some damage due to tension 425 : Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dtensi) 426 : c print *, Afpstmp 427 : c [5 428 : if ( Afpstmp .gt. ZERO ) then 428 : call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn) 429 : c [6 430 : if (ireturn .ne. 0) return 431 : c \--> ERROR 432 : c 6] 433 : c print *,'RSSTMP', RSSTMP 434 : c [6 435 : if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then 436 : c \--> YIELDING IN SHEAR 437 : Afps = Afpstmp 438 : print *,'YIELDING IN SHEAR 2' 439 : RSS = RSSTMP 440 : Tfs = sign(Chom*RSS , Tfs) 441 : if (Tfs0 .eq. ZERO) then 442 : Tfs1 = ZERO 443 : Tfs2 = ZERO 444 : else 445 : Tfs1= Tfs*Tfs1/Tfs0 446 : Tfs2= Tfs*Tfs2/Tfs0 447 : endif 448 : print *,'updated shear stress ',Tfs 449 : print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 450 : endif 451 : c 6] 452 : endif 453 : c 5] 454 : endif 455 : c 4] 456 : else 457 : c \--> check shear yielding in compression region 458 : Chom = (ETf - Tfc)/(ETf)*Dcompr 459 : c print *,'Coef. hom. ', Chom 460 : c print *,'???', Dcompr, ONE 461 : c [4 462 : if ( Dcompr .eq. ONE) then 463 : c \--> no crushing damage 464 : c [5 465 : if ( abs(Tfs0) .gt. (RSS*Chom) ) then 466 : c \--> YIELDING IN SHEAR 467 : c print *,'YIELDING IN SHEAR' 468 : c compute the total plastic shear strain... 469 : Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks 470 : c print *,'plasti shear strain ', Afps 471 : c ...and the final cohesion residual 471 : call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn) 472 : c [6 473 : if (ireturn .ne. 0) return 474 : c \--> ERROR 475 : c 6] 476 : Tfs = sign(Chom*RSS , Tfs) 477 : if (Tfs0 .eq. ZERO) then 478 : Tfs1 = ZERO 479 : Tfs2 = ZERO 480 : else 481 : Tfs1= Tfs*Tfs1/Tfs0 482 : Tfs2= Tfs*Tfs2/Tfs0 483 : endif 484 : c print *,'updated shear stress ',Tfs 485 : c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 486 : endif 487 : c 5] 488 : else 489 : c \--> already got some damage due to compression 490 : Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dcompr) 491 : c print *, Afpstmp 492 : c [5 493 : if ( Afpstmp .gt. ZERO ) then 493 : call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn) 494 : c [6 495 : if (ireturn .ne. 0) return 496 : c \--> ERROR 497 : c 6] 498 : c print *,'RSSTMP', RSSTMP 499 : c [6 500 : if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then 501 : c \--> YIELDING IN SHEAR 502 : Afps = Afpstmp 503 : c print *,'YIELDING IN SHEAR 2' 504 : RSS = RSSTMP 505 : Tfs = sign(Chom*RSS , Tfs) 506 : if (Tfs0 .eq. ZERO) then 507 : Tfs1 = ZERO 508 : Tfs2 = ZERO 509 : else 510 : Tfs1= Tfs*Tfs1/Tfs0 511 : Tfs2= Tfs*Tfs2/Tfs0 512 : endif 513 : c print *,'updated shear stress ',Tfs 514 : c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2 515 : endif 516 : c 6] 517 : endif 518 : c 5] 519 : endif 520 : c 4] 521 : endif 522 : c 3] 523 : else 524 : c \--> the edge has reach the origin... no homothetic rule 525 : c can be computed ! 526 : print *,'EDGE ZERO -- ERROR' 527 : endif 528 : c 2] 529 : endif 530 : c 1] 531 : c______________________________________________________________________ 532 : c 533 : c final stage --> updating the output variables... 534 : c ... first the final stress... 535 : c [1 536 : if ( JCT .eq. 0 ) then 537 : sigf(Inorm) = Tft 538 : sigi(Inorm) = Tft 539 : else 540 : sigf(Inorm) = Tfc 541 : sigi(Inorm) = Tfc 542 : endif 543 : c 1] 544 : sigf(Ishea1) = Tfs1 545 : sigi(Ishea1) = Tfs1 546 : sigf(Ishea2) = Tfs2 547 : sigi(Ishea2) = Tfs2 548 : c ... then the final internal variables... 549 : varf(Ipress) = Afpt 550 : varf(Ishear) = Afps 551 : varf(Icompr) = Afpc 552 : varf(Itotal) = Aftt 553 : c ... and finally, the plastic strain increments. 554 : c [1 555 : if ( JCT .eq. 0 ) then 556 : ddepsp(Inorm) = varf(Ipress) - vari(Ipress) 557 : else 558 : ddepsp(Inorm) = varf(Icompr) - vari(Icompr) 559 : endif 560 : c 1] 561 : ddepsp(Ishea1) = varf(Ishear) - vari(Ishear) 562 : ddepsp(Ishea2) = varf(Ishear) - vari(Ishear) 563 : vari(Ipress) = Afpt 564 : vari(Ishear) = Afps 565 : vari(Icompr) = Afpc 566 : vari(Itotal) = Aftt 567 : depsp(Inorm) = depsp(Inorm) + ddepsp(Inorm) 568 : depsp(Ishea1) = depsp(Ishea1) + ddepsp(Ishea1) 569 : depsp(Ishea2) = depsp(Ishea2) + ddepsp(Ishea2) 570 : c END OF MAIN LOOP 571 : c 0] 572 : enddo 573 : c print *,'output variables :' 574 : c print *,' --> stress state - normal shear' 575 : c print *, sigf(Inorm), sigf(Ishea1), sigf(Ishea2) 576 : c print *,'internal variables - normal - comp - shear - total' 577 : c print *, varf(Ipress), varf(Icompr),varf(Ishear), varf(Itotal) 578 : c print *,'plastic strain - normal - shear' 579 : c print *, depsp(Inorm), depsp(Ishea1),depsp(Ishea2) 580 : c And, that's all folks !!!! 581 : c print *,'sjonl2 ended normaly...' 582 : end 583 : 584 : 585 :
© Cast3M 2003 - All rights reserved.
Disclaimer