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