1 : C DPRAN3D SOURCE AM 12/02/27 23:39:38 7289 2 : SUBROUTINE dpran3d(parahot3,idimpara3,H66t,rt1,rc1,esigma, 3 : . deps,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,n,lnoconv, 4 : . istep,lcp,lbroyden,H66,lfulldamage) 5 : 6 : c parahot3 is the array that contains the material parameters and state variables 7 : c idimpara3 is the number of parameters for parahot3 8 : c H66t is the tangent matrix 9 : c rt1 and rc1 are the vectors used in DAMTG3D for the calculation of the tangent matrix 10 : c esigma is the effective stress vector 11 : c deps is the incremental strains vector 12 : c lg and lerror are logical 13 : c ntot is the integration point (global number for the whole structure): counter 14 : c n, istep, lnoconv, lbroyden and lfulldamage are used for the algorithmic strategy (step subdivider) 15 : c H66 is the elastic matrix 16 : 17 : c =========================================================================================== 18 : c ! Plastic Box of the 3D plastic-damage concrete model by T. Gernay ! 19 : c ! Yield functions of Drucker Prager in compression and Rankine in tension ! 20 : c ! Calculate the effective stresses and the effective elastoplastic tangent ! 21 : c ! modulus Dt (implicit process) ! 22 : c ! ! 23 : c ! input ! 24 : c ! DEPS : Incremental strain vector (6 components) ! 25 : c ! ! 26 : c ! output ! 27 : c ! H66t : tangent stiffness matrix (6x6) ! 28 : c ! ! 29 : c ! input and output ! 30 : c ! PARAHOT3 : the first columns contain the material properties at elevated temperature ! 31 : c ! the last columns contain strains, stresses, .... ! 32 : c ! ! 33 : c ! ESIGMA : Effective stress vector (6 components) ! 34 : c ! ! 35 : c =========================================================================================== 36 : 37 : IMPLICIT REAL*8 (A-B,D-H,O-Z) 38 : implicit integer (I-K,M,N) 39 : implicit logical (L) 40 : implicit character*10 (C) 41 : 42 : dimension deps(i6) 43 : c Incremental strain since the end of the previous time step (input) 44 : c i.e. since the last converged point 45 : dimension H66t(i6,i6) 46 : c Tangent matrix (output) 47 : dimension parahot3(idimpara3) 48 : c parahot3: Material prop. and various info. at elevated temp. (input and output) 49 : c parahot3(idimpara3-37,ntot): equivalent plastic strain in tension at the last converged step 50 : c parahot3(idimpara3-35,ntot): equivalent plastic strain in tension at the end of the current step 51 : c parahot3(idimpara3-36,ntot): equivalent plastic strain in compression at the last converged step 52 : c parahot3(idimpara3-34,ntot): equivalent plastic strain in compression at the end of the current step 53 : c parahot3(idimpara3-30:idimpara3-25,ntot) = strain at the end of the previous time step 54 : c parahot3(idimpara3-24:idimpara3-19,ntot) = plastic strains at the beginning of dpran3D 55 : c parahot3(idimpara3-18:idimpara3-13,ntot) = strain at the end of dpran3D 56 : c parahot3(idimpara3- 6:idimpara3- 1,ntot) = effective stress at the end of dpran3D 57 : c parahot3(idimpara3,ntot) = thermal strain 58 : dimension esigma(i6) 59 : c Effective stress (output) 60 : **** dimension EPSMEC(i6),EPSSIG(i6),esigmae(i6),esigmaf(i6) 61 : dimension EPSMEC(6),EPSSIG(6),esigmae(6),esigmaf(6) 62 : **** dimension H66(i6,i6),x(i2),d2gt(i6,i6),trav66b(i6,i6),trav6(i6) 63 : dimension H66(6,6),x(2),d2gt(6,6),trav66b(6,6),trav6(6) 64 : **** dimension trav66(i6,i6),Da(i6,i6),dft(i6) 65 : dimension trav66(6,6),Da(6,6),dft(6) 66 : **** dimension dfc(i6),dgc(i6),dgt(i6) 67 : dimension dfc(6),dgc(6),dgt(6) 68 : **** dimension trav1(1),DftD(i1,i6),DfcD(i1,i6),Ddgt(i6),Ddgc(i6) 69 : dimension trav1(1),DftD(1,6),DfcD(1,6),Ddgt(6),Ddgc(6) 70 : **** dimension trav16(i1,i6),rn1(i1,i6),rn2(i1,i6),rn3(i1,i6) 71 : dimension trav16(1,6),rn1(1,6),rn2(1,6),rn3(1,6) 72 : **** dimension rt1(i1,i6),rc1(i1,i6),vloc1(i1),vloc(i6) 73 : dimension rt1(1,6),rc1(1,6),vloc1(1),vloc(6) 74 : **** dimension Ddgtn1(i6,i6),Ddgtn2(i6,i6),Ddgcn3(i6,i6) 75 : dimension Ddgtn1(6,6),Ddgtn2(6,6),Ddgcn3(6,6) 76 : **** dimension Ddgcn4(i6,i6),rn4(i1,i6),dplas(6),d2gc(6,6),vloc2(6) 77 : dimension Ddgcn4(6,6),rn4(1,6),dplas(6),d2gc(6,6),vloc2(6) 78 : 79 : 80 : c Deviatoric projection operator P2 (stressdeviator = P2 * stress) 81 : c -------------------------------- 82 : dimension P2(6,6) 83 : data P2 / 2.0d0, -1.0d0, -1.0d0, 0.0d0, 0.0d0, 0.0d0, 84 : . -1.0d0, 2.0d0, -1.0d0, 0.0d0, 0.0d0, 0.0d0, 85 : . -1.0d0, -1.0d0, 2.0d0, 0.0d0, 0.0d0, 0.0d0, 86 : . 0.0d0, 0.0d0, 0.0d0, 6.0d0, 0.0d0, 0.0d0, 87 : . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 6.0d0, 0.0d0, 88 : . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 6.0d0 / 89 : 90 : 91 : c Matrix of unity I, or U 92 : c ----------------------- 93 : dimension U(6,6) 94 : data U / 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 95 : . 0.0d0, 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 96 : . 0.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0, 0.0d0, 97 : . 0.0d0, 0.0d0, 0.0d0, 1.0d0, 0.0d0, 0.0d0, 98 : . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0, 0.0d0, 99 : . 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1.0d0 / 100 : 101 : c Vector pi2 102 : c ----------------------- 103 : dimension pi2(6) 104 : data pi2 / 1.0d0, 105 : . 1.0d0, 106 : . 1.0d0, 107 : . 0.0d0, 108 : . 0.0d0, 109 : . 0.0d0 / 110 : 111 : 112 : 113 : ***** 114 : * MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS) 115 : IF(I1.GT.1) PRINT *, ' DPRAND3D - ERREUR I1 = ', I1, ' > 1 ' 116 : IF(I6.GT.6) PRINT *, ' DPRAND3D - ERREUR I6 = ', I6, ' > 6 ' 117 : ***** 118 : 119 : lerror = .false. 120 : lnoconv = .false. 121 : 122 : P2(1,1) = 2.0d0 123 : P2(2,2) = 2.0d0 124 : P2(3,3) = 2.0d0 125 : i7 = 7 126 : i8 = 8 127 : i9 = 9 128 : i10=10 129 : i11=11 130 : i12=12 131 : i13=13 132 : i14=14 133 : i15=15 134 : i16=16 135 : i17=17 136 : i18=18 137 : i19=19 138 : i20=20 139 : i21=21 140 : i22=22 141 : i23=23 142 : i24=24 143 : i25=25 144 : i26=26 145 : i27=27 146 : i28=28 147 : i29=29 148 : i30=30 149 : i31=31 150 : i32=32 151 : i33=33 152 : i34=34 153 : i35=35 154 : i36=36 155 : i37=37 156 : i38=38 157 : i39=39 158 : 159 : r0 = 0. 160 : r1 = 1. 161 : r2 = 2. 162 : r10 = 10. 163 : 164 : 165 : Eo = parahot3(i1) 166 : poison = parahot3(i2) 167 : fc = parahot3(i3) 168 : ft = parahot3(i4) 169 : c Strain at peak stress 170 : epsc1 = parahot3(i5) 171 : c Compressive limit of elasticity 172 : fc0 = parahot3(i6) 173 : c Biaxial compressive strength 174 : fb = parahot3(i7) 175 : c Dilatancy parameter 176 : ag = parahot3(i8) 177 : c Damage to peak stress 178 : dc = parahot3(i9) 179 : c Compressive dissipated energy parameter 180 : rxc = parahot3(i10) 181 : c Crack energy in tension 182 : gt = parahot3(i11) 183 : c Model parameters 184 : rkappa1 = parahot3(i13) 185 : ac = parahot3(i14) 186 : bc = parahot3(i15) 187 : 188 : c small strain for starting the return mapping algorithm 189 : tol = 0.0001*epsc1 190 : 191 : c Hardening parameter for tension at the beginning of the step 192 : rkappait= parahot3(idimpara3-37) 193 : c Hardening parameter for compression at the beginning of the step 194 : rkappaic= parahot3(idimpara3-36) 195 : 196 : c =================================================================== 197 : c ! ! 198 : c ! CALCULATION OF THE MATRIX OF ELASTICITY ! 199 : c ! ======================================= ! 200 : c =================================================================== 201 : 202 : c Matrix of elasticity, (6x6, full 3D) 203 : c -------------------- 204 : c [Eo*(1-nu)/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) 0 0 0 ] 205 : c [Eo*nu/((1+nu)*(1-2*nu)) Eo*(1-nu)/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) 0 0 0 ] 206 : c H66 = [Eo*nu/((1+nu)*(1-2*nu)) Eo*nu/((1+nu)*(1-2*nu)) Eo*(1-nu)/((1+nu)*(1-2*nu)) 0 0 0 ] 207 : c [0 0 0 Eo/(1+nu) 0 0 ] 208 : c [0 0 0 0 Eo/(1+nu) 0 ] 209 : c [0 0 0 0 0 Eo/(1+nu) ] 210 : 211 : c Matrix of elasticity, (3x3, plane stress) 212 : c -------------------- 213 : c [ 1 nu 0 0 0 0] 214 : c H66 = Eo/(1-nu*nu) [nu 1 0 0 0 0] 215 : c [ 0 0 0 0 0 0] 216 : c [ 0 0 0 (1-nu)/2 0 0] 217 : c [ 0 0 0 0 0 0] 218 : c [ 0 0 0 0 0 0] 219 : c++ 220 : 221 : do jloc=i1,i6 222 : do iloc=i1,i6 223 : H66(iloc,jloc) = r0 224 : end do 225 : end do 226 : 227 : if (.not.lcp) then 228 : c 3D model 229 : rloc = Eo/((r1+poison)*(r1-2.0d0*poison)) 230 : H66(i1,i1) = rloc*(r1-poison) 231 : H66(i2,i2) = rloc*(r1-poison) 232 : H66(i3,i3) = rloc*(r1-poison) 233 : H66(i1,i2) = rloc*poison 234 : H66(i1,i3) = rloc*poison 235 : H66(i2,i3) = rloc*poison 236 : H66(i2,i1) = rloc*poison 237 : H66(i3,i1) = rloc*poison 238 : H66(i3,i2) = rloc*poison 239 : rloc2 = Eo/(2.0d0*(r1+poison)) 240 : H66(i4,i4) =rloc2 241 : H66(i5,i5) =rloc2 242 : H66(i6,i6) =rloc2 243 : else 244 : c PLANE STRESS 245 : rloc = Eo/(r1-poison*poison) 246 : H66(i1,i1) = rloc 247 : H66(i2,i1) = poison * (rloc) 248 : H66(i1,i2) = poison * (rloc) 249 : H66(i2,i2) = rloc 250 : H66(4,4) = rloc * ((1.0d0-poison)/(2.0d0)) 251 : deps(3) = 0.d0 252 : deps(5) = 0.d0 253 : deps(6) = 0.d0 254 : endif 255 : 256 : c ============================================================= 257 : c ! ! 258 : c ! ELASTIC STRESS PREDICTOR ! 259 : c ! ======================== ! 260 : c ============================================================= 261 : 262 : c Total mechanical strain: EPSMEC[i,j] = epsmec[i-1,conv] + deps[i,j] - epspl[i-1,j] 263 : c epsplastic has to be subtracted because Sigma initial = 0. 264 : do iloc=i1,i6 265 : EPSMEC(iloc) = parahot3(idimpara3-31+iloc) 266 : . + ( (DEPS(iloc) * istep)/n ) 267 : . - parahot3(idimpara3-25+iloc) 268 : c Calculation of the INSTANTANEOUS STRESS-DEPENDENT STRAIN 269 : EPSSIG(iloc) = EPSMEC(iloc) - parahot3(17+iloc) 270 : end do 271 : c The effective stresses are calculated as functions of the instantaneous 272 : c stress-dependent strain eps,sig because the transient creep strain is explicitly 273 : c calculated apart from these constitutive relationships. Still in parahot, the mechanical 274 : c strains are saved, not the instantaneous stress-dependent strains. 275 : 276 : 277 : if ( (.not.lnoconv).and.(n.gt.1) ) then 278 : c We are in a subdivision of the step (for algorithmic reasons) 279 : c The plastic strain considered in epsmec has to be adjusted 280 : rkappait = parahot3(idimpara3-35) 281 : rkappaic = parahot3(idimpara3-34) 282 : x1 = rkappait - parahot3(idimpara3-37) 283 : x2 = rkappaic - parahot3(idimpara3-36) 283 : call dplas3D(x1,x2,P2,pi2,ag,fc,fb,parahot3,idimpara3, 284 : . dplas,lcp) 285 : do iloc=i1,i6 286 : EPSSIG(iloc) = EPSMEC(iloc)-dplas(iloc)-parahot3(17+iloc) 287 : end do 288 : endif 289 : 290 : c Elastic stress predictor: eff sigma[i,j] = Hel x EPSSIG[i,j] 290 : call mulab(H66,EPSSIG,ESIGMAE,i6,i6,i6,i1) 291 : 292 : * 293 : * print *,' prediction contraintes = ',esigmae 294 : * 295 : 296 : c ============================================================== 297 : c ! ! 298 : c ! CONCRETE COMPLETELY DAMAGED ! 299 : c ! =========================== ! 300 : c ============================================================== 301 : dc = parahot3(idimpara3-40) 302 : dt = parahot3(idimpara3-41) 303 : 304 : if ((dc.gt.0.98).or.(lfulldamage).or.(dt.gt.0.98)) then 305 : c Concrete completely damaged 306 : do iloc=1,6 307 : ESIGMA(iloc) = 0.d0 308 : end do 309 : do iloc=1,6 310 : parahot3(idimpara3-13+iloc) = esigma(iloc) 311 : end do 312 : do jloc=1,6 313 : do iloc=1,6 314 : H66t(iloc,jloc) = 0.d0 315 : end do 316 : end do 317 : c store the mechanical strain at the end of this iteration 318 : c mechanical strain = epsmec + epspl 319 : do iloc=i1,i6 320 : parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + 321 : . parahot3(idimpara3-25+iloc) 322 : end do 323 : return 324 : endif 325 : 326 : c ========================================================================== 327 : c ! ! 328 : c ! CALCULATE THE PLASTICITY FUNCTIONS WITH THE ELASTIC STRESS PREDICTOR ! 329 : c ! ==================================================================== ! 330 : c ========================================================================== 331 : 332 : c Calculation of the hardening functions 333 : tau_compr = fhardcc(rkappaic,fc,fc0,rkappa1, 334 : . ac,bc,lerror) 335 : tau_tension = fhardct(rkappait,ft) 336 : if (lerror) then 337 : write(2,*) '>Error returned from compressive hardening 338 : . to DPRAN3d' 339 : return 340 : endif 341 : 342 : c Calculation of the yield functions for the elastic stress predictor 343 : alpha = (fb-fc)/(2.*fb-fc) 344 : f_compr = DPfunc(esigmae,alpha) 345 : 346 : f_tension = Ranfunc(esigmae,i3,i6,parahot3, 347 : . idimpara3,lcp) 348 : 349 : c++ 350 : f_compr = f_compr - ( (1.0d0-alpha) * tau_compr ) 351 : f_tension = f_tension - tau_tension 352 : 353 : c ======================================================= 354 : c ! ! 355 : c ! TEST ELASTIC OR PLASTIC STEP ! 356 : c ! ============================ ! 357 : c ======================================================= 358 : 359 : c The plasticity functions must be satisfied with a precision of 10 Pa, i.e. 0.000010 N/mm² 360 : if ((f_compr.lt.-r10).and.(f_tension.lt.-r10)) then 361 : 362 : C *************** 363 : c ! ELASTICITY ! 364 : C *************** 365 : 366 : do iloc=i1,i6 367 : ESIGMA(iloc) = ESIGMAE(iloc) 368 : end do 369 : 370 : do jloc=i1,i6 371 : do iloc=i1,i6 372 : H66T(iloc,jloc) = H66(iloc,jloc) 373 : end do 374 : end do 375 : c store the stress at the end of this iteration 376 : do iloc=i1,i6 377 : parahot3(idimpara3-7+iloc) = ESIGMAE(iloc) 378 : end do 379 : 380 : c The equivalent plastic strain in tension has not changed 381 : parahot3(idimpara3-35) = parahot3(idimpara3-37) 382 : c The equivalent plastic strain in compression has not changed 383 : parahot3(idimpara3-34) = parahot3(idimpara3-36) 384 : 385 : c vectors for the calculation of H66t in subroutine damtg3d 386 : do iloc=i1,i6 387 : rt1(i1,iloc) = r0 388 : rc1(i1,iloc) = r0 389 : end do 390 : 391 : else 392 : 393 : C *************** 394 : c ! PLASTICITY ! 395 : C *************** 396 : 397 : * PRINT *, ' >>>>>>>> ON ENTRE EN PLASTICITE' 398 : 399 : c RETURN MAPPING 400 : c -------------- 401 : c Calculation of Dlambdac and Dlambdat, the increments of the eq. plastic strain which brings 402 : c back on the plasticity surfaces 403 : c We have to solve the following non linear system of equations 404 : c rkt * function_t(Dlambdac,Dlambdat) + (1-rkt) * Dlambdat = 0 405 : c rkc * function_c(Dlambdac,Dlambdat) + (1-rkc) * Dlambdac = 0 406 : c with rkc and rkt = 1 or 0, depending whether the related surface of plasticity is violated 407 : 408 : if (f_compr.lt.-r10) then 409 : c only the tension surface is violated by the trial stress 410 : Dlambdac = r0 411 : rkc = r0 412 : rkt = r1 413 : 414 : else if (f_tension.lt.-r10) then 415 : c only the compression surface is violated by the trial stress 416 : Dlambdat = r0 417 : rkc = r1 418 : rkt = r0 419 : 420 : else 421 : c Both surfaces are violated by the trial stress 422 : rkc = r1 423 : rkt = r1 424 : 425 : endif 426 : 427 : c X=0 at the beginning of every SAFIR iteration, X being the vector (Dlambdat,Dlambdac) 428 : X(i1) = r0 429 : X(i2) = r0 430 : 431 : c =================================================================================== 432 : c ! ! 432 : c ! CALL TO SUBROUTINE ALGOPLAS TO SOLVE THE NONLINEAR SYSTEM IN DLAMBDAT, DLAMBDAC ! 433 : c ! =============================================================================== ! 434 : c =================================================================================== 435 : 436 : if (.not.lcp) then 437 : c 3D MODEL 438 : if (.not.lbroyden) then 438 : call Algoplas1(rkt,rkc,esigmae,H66,P2,pi2,alpha, 439 : . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, 440 : . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, 441 : . parahot3,idimpara3,deps,lnoconv,lcp,U) 442 : else 442 : call Algoplas2(rkt,rkc,esigmae,H66,P2,pi2,alpha, 443 : . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, 444 : . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, 445 : . parahot3,idimpara3,deps,lnoconv,lcp,U) 446 : endif 447 : if (lnoconv) then 448 : return 449 : endif 450 : 451 : if (lerror) then 452 : write(2,*) 'lerror is returned from ALGOPLAS to 453 : . DPRAN3D. Ntot =',ntot 454 : return 455 : endif 456 : 457 : c The increment of plastic multipliers Dlambdat, Dlambdac have been found 458 : 459 : c ========================================================== 460 : c ! ! 461 : c ! UPDATE OF THE HARDENING PARAMETERS AND THE STRESSES ! 462 : c ! =================================================== ! 463 : c ========================================================== 464 : 465 : Dlambdat = x(i1) 466 : Dlambdac = x(i2) 467 : 468 : if (n.gt.1) then 469 : c we are currently in a step subdivision (for algorithmic reason) 470 : c we take the value saved at the end of the previous step-subdivision 471 : rkappait = parahot3(idimpara3-35) 472 : rkappaic = parahot3(idimpara3-34) 473 : else 474 : c we take the value at the previous (converged) time step 475 : rkappait= parahot3(idimpara3-37) 476 : rkappaic= parahot3(idimpara3-36) 477 : endif 478 : 479 : c update of the hardening parameters 480 : rkappat = rkappait+Dlambdat 481 : rkappac = rkappaic+Dlambdac 482 : 483 : c Update of the stresses: esigmaf = F(Dlambdat,Dlambdac (and esigmae)) 484 : ag = parahot3(8) 484 : call sigf3D_implicit(esigmae,Dlambdat,Dlambdac,esigmaf,H66, 485 : c *************** 486 : . P2,pi2,ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp, 487 : . esigmae,U) 488 : 489 : if (lerror) then 490 : return 491 : endif 492 : 493 : else 494 : 495 : c PLANE STRESS MODEL 496 : if (.not.lbroyden) then 496 : call Algoplas1cp(rkt,rkc,esigmae,H66,P2,pi2,alpha, 497 : . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, 498 : . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, 499 : . parahot3,idimpara3,deps,esigmaf,rkappat,rkappac, 500 : . lnoconv,lcp,U) 501 : else 501 : call Algoplas2cp(rkt,rkc,esigmae,H66,P2,pi2,alpha, 502 : . rkappait,rkappaic,x,fc,fb,ft,Eo,poison,lg,ntot, 503 : . fc0,rkappa1,ac,bc,ag,lerror,i1,i2,i3,i4,i5,i6, 504 : . parahot3,idimpara3,deps,esigmaf,rkappat,rkappac, 505 : . lnoconv,lcp,U) 506 : endif 507 : 508 : if (lnoconv) then 509 : return 510 : endif 511 : 512 : if (lerror) then 513 : write(2,*) 'lerror is returned from ALGOPLAS to 514 : . DPRAN3D. Ntot =',ntot 515 : return 516 : endif 517 : 518 : Dlambdat = x(i1) 519 : Dlambdac = x(i2) 520 : if (n.gt.1) then 521 : c we are currently in a step subdivision (for algorithmic reason) 522 : c we take the value saved at the end of the previous step-subdivision 523 : rkappait = parahot3(idimpara3-35) 524 : rkappaic = parahot3(idimpara3-34) 525 : else 526 : c we take the value at the previous (converged) time step 527 : rkappait= parahot3(idimpara3-37) 528 : rkappaic= parahot3(idimpara3-36) 529 : endif 530 : 531 : c update of the hardening parameters 532 : rkappat = rkappait+Dlambdat 533 : rkappac = rkappaic+Dlambdac 534 : 535 : endif 536 : 537 : c store the stress at the end of this iteration 538 : do iloc=i1,i6 539 : ESIGMA(iloc) = ESIGMAF(iloc) 540 : end do 541 : do iloc=i1,i6 542 : parahot3(idimpara3-i7+iloc) = ESIGMAF(iloc) 543 : end do 544 : 545 : c store the accumulated plastic strains at the end of this iteration 546 : parahot3(idimpara3-35) = rkappat 547 : parahot3(idimpara3-34) = rkappac 548 : 549 : c THE EFFECTIVE STRESSES HAVE BEEN FOUND 550 : c ************************************** 551 : 552 : c ================================================================ 553 : c ! TANGENT MATRIX. (effective!) ! 554 : c ! ---------------------------- ! 555 : c ! Here, we calculate the effective elastoplastic tangent ! 556 : c ! modulus Dt. As this is a damage-plastic model, we have ! 557 : c ! then to correct this modulus Dt with the damage terms to ! 558 : c ! obtain the (nominal) consistent algorithmic tangent modulus ! 559 : c ! Ct. This will be done in DAMTG3D. Here, we calculate Dt. ! 560 : c ================================================================ 561 : 562 : c APEX of Drucker-Prager surface 563 : alpha = (fb-fc)/(2.*fb-fc) 564 : f_compr = DPfunc(esigmae,alpha) 565 : if ((dabs(f_compr).le.100.).and.(rkc.eq.1)) then 566 : do jloc=i1,i6 567 : do iloc=i1,i6 568 : H66t(iloc,jloc) = 0.d0 569 : end do 570 : enddo 571 : do iloc=i1,i6 572 : parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + 573 : . parahot3(idimpara3-25+iloc) 574 : end do 575 : RETURN 576 : endif 577 : 578 : if (lcp) then 579 : c recalculate the full 3D matrix 580 : do jloc=i1,i6 581 : do iloc=i1,i6 582 : H66(iloc,jloc) = r0 583 : end do 584 : end do 585 : rloc = Eo/((r1+poison)*(r1-2.0d0*poison)) 586 : H66(i1,i1) = rloc*(r1-poison) 587 : H66(i2,i2) = rloc*(r1-poison) 588 : H66(i3,i3) = rloc*(r1-poison) 589 : H66(i1,i2) = rloc*poison 590 : H66(i1,i3) = rloc*poison 591 : H66(i2,i3) = rloc*poison 592 : H66(i2,i1) = rloc*poison 593 : H66(i3,i1) = rloc*poison 594 : H66(i3,i2) = rloc*poison 595 : rloc2 = Eo/(2.0d0*(r1+poison)) 596 : H66(i4,i4) =rloc2 597 : H66(i5,i5) =rloc2 598 : H66(i6,i6) =rloc2 599 : endif 600 : if (lcp) then 601 : tau_tension = fhardct(rkappat,ft) 602 : psi1 = tau_tension - (esigmaf(1)+esigmaf(2))/2 603 : if (abs(psi1).le.fc*1.e-5) then 604 : do jloc=1,6 605 : do iloc=1,6 606 : H66t(iloc,jloc) = 0 607 : enddo 608 : enddo 609 : do iloc=i1,i6 610 : parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + 611 : . parahot3(idimpara3-25+iloc) 612 : end do 613 : RETURN 614 : endif 615 : endif 616 : 617 : c First, calculation of Da, the effective algorithmic modulus 618 : c ************************************************** 619 : 619 : call ddsigrank(esigmaf,d2gt,i1,i2,i3,i4,i5,i6,lcp, 620 : . parahot3,idimpara3) 621 : 622 : do jloc=i1,i6 623 : do iloc=i1,i6 624 : TRAV66b(iloc,jloc) = Dlambdat*d2gt(iloc,jloc) 625 : c [6;6] 626 : end do 627 : end do 627 : call mulAB(H66,TRAV66b,TRAV66b,i6,i6,i6,i6) 628 : c [6;6] = [6;6] x [6;6] 629 : 630 : tau_c = fhardcc(rkappac,fc,fc0,rkappa1, 631 : . ac,bc,lerror) 632 : if (lerror) then 633 : write(2,*) '>Error returned from compressive hardening 634 : . to DPRAN3d' 635 : return 636 : endif 636 : call mulATB(pi2,esigmaf,vloc1,i6,i1,i6,i1) 637 : c [1] = [1;6] x [6;1] 638 : psi2 = (1d0-alpha)*tau_c - alpha*vloc1(i1) 639 : 640 : c First come the second order terms. 640 : call mulab(P2,ESIGMAF,TRAV6,i6,i6,i6,i1) 641 : c [6;1] = [6;6] x [6;1] 642 : 643 : c t t t 644 : c T66 = T6 T6 = P2 s s P2 644 : call mulABT(trav6,trav6,trav66,i6,i1,i6,i1) 645 : c [6;6] = [6;1] x [1;6] 646 : test = 1.E+6 647 : psi2lim=max(1.E-7*test,1.d-3) 648 : if (abs(psi2).le.psi2lim) then 649 : do jloc=i1,i6 650 : do iloc=i1,i6 651 : TRAV66(iloc,jloc) = r0 652 : end do 653 : end do 654 : else 655 : do jloc=i1,i6 656 : do iloc=i1,i6 657 : TRAV66(iloc,jloc) = (P2(iloc,jloc)-TRAV66(iloc,jloc) 658 : . /(r2*psi2*psi2))/(r2*psi2) 659 : c [6;6] 660 : end do 661 : end do 662 : endif 663 : 664 : do jloc=i1,i6 665 : do iloc=i1,i6 666 : TRAV66(iloc,jloc) = Dlambdac*TRAV66(iloc,jloc) 667 : end do 668 : end do 668 : call mulAB(H66,TRAV66,TRAV66,i6,i6,i6,i6) 669 : c **** ***** [6;6]=[6;6]x[6;6] 670 : do jloc=i1,i6 671 : do iloc=i1,i6 672 : TRAV66(iloc,jloc) = U(iloc,jloc)+TRAV66(iloc,jloc) 673 : c TRAV66(iloc,jloc)=U(iloc,jloc)+TRAV66(iloc,jloc) 674 : . +TRAV66b(iloc,jloc) 675 : c TRAV66b requires the calculation of d2gt - we prefer not to take it 676 : end do 677 : end do 678 : c t t -1 679 : c Da = ( I + Dlambdat H66 d2gt + Dlambdac H66 (P2/2psi2)-(P2 s s P2)/(4psi2³) ) H66 680 : 680 : call invert9(TRAV66,i6,i6) 681 : c **** ******* [6;6] 681 : call mulab(TRAV66,H66,Da,i6,i6,i6,i6) 682 : c [6;6]=[6;6]x[6;6] 683 : 684 : c The effective algorithmic modulus Da has been obtained 685 : c Now, calculation of the effective elastoplastic tangent modulus Dt 686 : c ****************************************************************** 687 : c Now come the first order terms 688 : 689 : dtau1 = dfhardct(rkappat) 690 : dtau2 = dfhardcc(rkappac,fc,fc0,rkappa1,ac,bc) 691 : 692 : c Calculation of n1, n2, n3, n4 693 : 693 : call dsigRank(esigmaf,dgt,i3,i6,parahot3,idimpara3, 694 : c **** ******** [6;1] 695 : . rkappait,lcp) 696 : rloc = 0.d0 697 : do iloc=1,6 698 : rloc = rloc + (dgt(iloc))**2 699 : enddo 700 : if (rloc.eq.0.d0) rkt = 0.d0 701 : do iloc=i1,i6 702 : dft(iloc) = rkt*dgt(iloc) 703 : c [6;1] = [1] x [6;1] 704 : end do 705 : 706 : e1 = r1-rkt*(r1+dtau1) 707 : c [1] 708 : e2 = r1-rkc*(r1+dtau2) 709 : c [1] 710 : 711 : fc = parahot3(i3) 712 : c Compressive strength 713 : fb = parahot3(i7) 714 : c Biaxial compressive strength 715 : 715 : call dsigDrucker(fc,fb,esigmaf,dfc,alpha,P2, 716 : c **** **************** [6;1] 717 : . pi2,i1,i6,lcp) 718 : do iloc=i1,i6 719 : dfc(iloc) = rkc*dfc(iloc) 720 : c [6;1] = [1]x[6;1] 721 : end do 722 : ag = parahot3(i8) 722 : call dsiggc(esigmaf,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp) 723 : c **** **************** [6;1] 724 : 724 : call mulATB(dft,Da,trav16,i6,i1,i6,i6) 725 : c [1;6] = [1;6]x[6;6] 726 : do iloc=i1,i6 727 : dftD(i1,iloc) = trav16(i1,iloc) 728 : end do 729 : 729 : call mulab(Da,dgt,Ddgt,i6,i6,i6,i1) 730 : c [6;1] = [6;6] x [6;1] 731 : 731 : call mulATB(dfc,Da,trav16,i6,i1,i6,i6) 732 : c [1;6] = [1;6]x[6;6] 733 : do iloc=i1,i6 734 : dfcD(i1,iloc) = trav16(i1,iloc) 735 : end do 736 : 736 : call mulab(Da,dgc,Ddgc,i6,i6,i6,i1) 737 : c [6;1] = [6;6]x[6;1] 738 : 739 : c n1 739 : call mulATB(dft,Ddgc,trav1,i6,i1,i6,i1) 740 : c **** ****** [1]=[1;6]x[6;1] 741 : dftDdgc = trav1(i1) 742 : c [1] 743 : do iloc=i1,i6 744 : rn1(i1,iloc) = -dftDdgc*dfcD(i1,iloc) 745 : end do 746 : c [1;6]=[1]x[1;6] 747 : 748 : c n2 [1]=[1;6]x[6;1] 748 : call mulATB(dfc,Ddgc,trav1,i6,i1,i6,i1) 749 : c **** ****** 750 : dfcDdgc = trav1(i1) 751 : do iloc=i1,i6 752 : rn2(i1,iloc) = (dfcDdgc-e2)*dftD(i1,iloc) 753 : end do 754 : 755 : c n3 [1]=[1;6]x[6;1] 755 : call mulATB(dft,Ddgt,trav1,i6,i1,i6,i1) 756 : c **** ****** 757 : dftDdgt = trav1(i1) 758 : do iloc=i1,i6 759 : rn3(i1,iloc) = (dftDdgt-e1)*dfcD(i1,iloc) 760 : end do 761 : 762 : c n4 [1]=[1;6]x[6;1] 762 : call mulATB(dfc,Ddgt,trav1,i6,i1,i6,i1) 763 : c **** ****** 764 : dfcDdgt = trav1(i1) 765 : do iloc=i1,i6 766 : rn4(i1,iloc) = -(dfcDdgt)*dftD(i1,iloc) 767 : end do 768 : 769 : c denominator [1] 770 : Denom = (dftDdgt-e1)*(dfcDdgc-e2) - (dfcDdgt*dftDdgc) 771 : 772 : c Calculation of rt1 and rc1 that are used in Mater 2 for calculation of Ct 773 : do iloc=i1,i6 774 : rt1(i1,iloc) = (rn1(i1,iloc)+rn2(i1,iloc))/denom 775 : rc1(i1,iloc) = (rn3(i1,iloc)+rn4(i1,iloc))/denom 776 : end do 777 : 778 : do iloc=i1,i6 779 : vloc(iloc) = rn1(i1,iloc) 780 : end do 780 : call mulABT(Ddgt,vloc,Ddgtn1,i6,i1,i6,i1) 781 : c **** ****** [6;6]=[6;1]x[1;6] 782 : do iloc=i1,i6 783 : vloc(iloc) = rn2(i1,iloc) 784 : end do 784 : call mulABT(Ddgt,vloc,Ddgtn2,i6,i1,i6,i1) 785 : c **** ****** [6;6]=[6;1]x[1;6] 786 : do iloc=i1,i6 787 : vloc(iloc) = rn3(i1,iloc) 788 : end do 788 : call mulABT(Ddgc,vloc,Ddgcn3,i6,i1,i6,i1) 789 : c **** ****** [6;6]=[6;1]x[1;6] 790 : do iloc=i1,i6 791 : vloc(iloc) = rn4(i1,iloc) 792 : end do 792 : call mulABT(Ddgc,vloc,Ddgcn4,i6,i1,i6,i1) 793 : c **** ****** [6;6]=[6;1]x[1;6] 794 : 795 : do jloc=i1,i6 796 : do iloc=i1,i6 797 : TRAV66(iloc,jloc) = Ddgtn1(iloc,jloc) + Ddgtn2(iloc,jloc) 798 : . + Ddgcn3(iloc,jloc) + Ddgcn4(iloc,jloc) 799 : c [6;6] 800 : TRAV66(iloc,jloc) = TRAV66(iloc,jloc)/Denom 801 : c [6;6] 802 : end do 803 : end do 804 : 805 : do jloc=i1,i6 806 : do iloc=i1,i6 807 : H66t(iloc,jloc) = Da(iloc,jloc) - TRAV66(iloc,jloc) 808 : end do 809 : end do 810 : c END OF CALCULATION OF EFFECTIVE TANGENT MATRIX 811 : c ************************************************** 812 : 813 : endif 814 : c END OF PLASTICITY 815 : 816 : c store the mechanical strain at the end of this iteration 817 : c mechanical strain = epsmec + epspl 818 : do iloc=i1,i6 819 : parahot3(idimpara3-19+iloc) = EPSMEC(iloc) + 820 : . parahot3(idimpara3-25+iloc) 821 : end do 822 : 823 : RETURN 824 : END 825 : 826 :
© Cast3M 2003 - All rights reserved.
Disclaimer